text:
type GLFramebuffer
    id         ::GLuint
    color      ::Texture{RGBA{UFixed8}, 2}
    objectid   ::Texture{Vec{2, GLushort}, 2}
    depth      ::Texture{Float32, 2}
    postprocess::RenderObject
end
Base.size(fb::GLFramebuffer) = size(fb.color) # it's guaranteed, that they all have the same size

loadshader(name) = load(joinpath(dirname(@__FILE__), name))

function draw_fullscreen(vao_id)
    glBindVertexArray(vao_id)
    glDrawArrays(GL_TRIANGLES, 0, 3)
    glBindVertexArray(0)
end
immutable PostprocessPrerender
end
@compat function (sp::PostprocessPrerender)()
    glDisable(GL_DEPTH_TEST)
    glDisable(GL_BLEND)
    glDisable(GL_STENCIL_TEST)
    glStencilMask(0xff)
    glDisable(GL_CULL_FACE)
end
"""
Creates a postprocessing render object.
This will transfer the pixels from the color texture of the Framebuffer
to the screen and while at it, it can do some postprocessing (not doing it right now):
E.g fxaa anti aliasing, color correction etc.
"""
function postprocess(color::Texture, framebuffer_size)
    shader = LazyShader(
        loadshader("fxaa.vert"),
        loadshader("fxaa_combine.frag")
    )
    data = Dict{Symbol, Any}(
        :color_texture => color
    )
    robj = RenderObject(data, shader, PostprocessPrerender(), nothing)
    robj.postrenderfunction = () -> draw_fullscreen(robj.vertexarray.id)
    robj
end

function attach_framebuffer{T}(t::Texture{T, 2}, attachment)
    glFramebufferTexture2D(GL_FRAMEBUFFER, attachment, GL_TEXTURE_2D, t.id, 0)
end

function GLFramebuffer(fb_size)
    render_framebuffer = glGenFramebuffers()
    glBindFramebuffer(GL_FRAMEBUFFER, render_framebuffer)

    buffersize      = tuple(value(fb_size)...)
    color_buffer    = Texture(RGBA{UFixed8},    buffersize, minfilter=:nearest, x_repeat=:clamp_to_edge)
    objectid_buffer = Texture(Vec{2, GLushort}, buffersize, minfilter=:nearest, x_repeat=:clamp_to_edge)
    depth_buffer    = Texture(Float32, buffersize,
        internalformat = GL_DEPTH_COMPONENT32F,
        format         = GL_DEPTH_COMPONENT,
        minfilter=:nearest, x_repeat=:clamp_to_edge
    )
    attach_framebuffer(color_buffer, GL_COLOR_ATTACHMENT0)
    attach_framebuffer(objectid_buffer, GL_COLOR_ATTACHMENT1)
    attach_framebuffer(depth_buffer, GL_DEPTH_ATTACHMENT)

    p  = postprocess(color_buffer, fb_size)
    fb = GLFramebuffer(render_framebuffer, color_buffer, objectid_buffer, depth_buffer, p)
    glBindFramebuffer(GL_FRAMEBUFFER, 0)
    fb
end

function Base.resize!(fb::GLFramebuffer, window_size)
    ws = tuple(window_size...)
    if ws!=size(fb) && all(x->x>0, window_size)
        buffersize = tuple(window_size...)
        resize_nocopy!(fb.color, buffersize)
        resize_nocopy!(fb.objectid, buffersize)
        resize_nocopy!(fb.depth, buffersize)
    end
    nothing
end

immutable MonitorProperties
    name::Compat.UTF8String
    isprimary::Bool
    position::Vec{2, Int}
    physicalsize::Vec{2, Int}
    videomode::GLFW.VidMode
    videomode_supported::Vector{GLFW.VidMode}
    dpi::Vec{2, Float64}
    monitor::Monitor
end

function MonitorProperties(monitor::Monitor)
    name 		 = GLFW.GetMonitorName(monitor)
    isprimary 	 = GLFW.GetPrimaryMonitor() == monitor
    position	 = Vec{2, Int}(GLFW.GetMonitorPos(monitor)...)
    physicalsize = Vec{2, Int}(GLFW.GetMonitorPhysicalSize(monitor)...)
    videomode 	 = GLFW.GetVideoMode(monitor)

    dpi			 = Vec(videomode.width * 25.4, videomode.height * 25.4) ./ Vec{2, Float64}(physicalsize)
    videomode_supported = GLFW.GetVideoModes(monitor)

    MonitorProperties(name, isprimary, position, physicalsize, videomode, videomode_supported, dpi, monitor)
end

immutable GLContext
    window::GLFW.Window
    framebuffer::GLFramebuffer
end

type Screen
    name 		::Symbol
    area 		::Signal{SimpleRectangle{Int}}
    parent 		::Screen
    children 	::Vector{Screen}
    inputs 		::Dict{Symbol, Any}
    renderlist 	::Tuple # a tuple of specialized renderlists

    hidden 		::Bool
    color       ::RGBA{Float32}

    cameras 	::Dict{Symbol, Any}

    glcontext   ::GLContext

    function Screen(
            name        ::Symbol,
            area        ::Signal{SimpleRectangle{Int}},
            parent 		::Screen,
            children 	::Vector{Screen},
            inputs 		::Dict{Symbol, Any},
            renderlist 	::Tuple,
            hidden 		::Bool,
            color       ::Colorant,
            cameras 	::Dict{Symbol, Any},
            context     ::GLContext
        )
        new(
            name, area, parent,
            children, inputs, renderlist,
            hidden, RGBA{Float32}(color), cameras,
            context
        )
    end

    function Screen(
            name        ::Symbol,
            area        ::Signal{SimpleRectangle{Int}},
            children    ::Vector{Screen},
            inputs      ::Dict{Symbol, Any},
            renderlist  ::Tuple,
            hidden      ::Bool,
            color       ::Colorant,
            cameras     ::Dict{Symbol, Any},
            context     ::GLContext
        )
        parent = new()
        new(
            name, area, parent,
            children, inputs, renderlist,
            hidden, RGBA{Float32}(color), cameras,
            context
        )
    end
end


text:
module EchoviewEvr

using Filetimes
using Dates

export regions, polygon, polygons

"""
An Echoview region.

See [documentation](http://bit.ly/2uH0O4a).
"""
struct Region
    name::String
    version::String
    classification::String
    boundingrectangle
    notes::String
    detectionsettings::String
    points
    regiontype::String
end

function readheader(f)
    line = readline(f)
    # Note that the file actually starts with a byte order mark U+FEFF
    @assert startswith(line, "\ufeffEVRG") "Invalid EVRG header."
    line
end

function readmultiline(f)
    n = parse(Int, readline(f))
    lines = String[]
    for i = 1:n
        line = readline(f)
        push!(lines, line)
    end
    join(lines)
end

function readpoints(f)
    line = readline(f)
    fields = split(line)
    # BUG points = fields[1:2:end-1]
    points = fields[1:end-1]
    regiontype = fields[end]
    return points, regiontype
end

function readblankline(f)
    line = readline(f)
    @assert length(line) < 1
end

function readregion(f)
    readblankline(f)

    line = readline(f)
    fields = split(line)
    version = fields[1]
    @assert version=="13"
    pointcount = fields[2]
    regionid = fields[3]
    selected = fields[4]
    @assert selected == "0"
    regioncreationtype = fields[5]
    dummy = fields[6]
    @assert dummy == "-1"
    boundingrectanglecalculated = fields[7]
    boundingrectangle = []
    if boundingrectanglecalculated == "1"
        # Not enitirely consistent with the documentation here?
        datep1 = fields[8]
        timep1 = fields[9]
        depthp1 = fields[10]
        datep2 = fields[11]
        timep2 = fields[12]
        depthp2 = fields[13]
        boundingrectangle = [datep1, timep1, depthp1,
                              datep2, timep2, depthp2]
    end

    notes = readmultiline(f)
    detectionsettings = readmultiline(f)
    classification = readline(f)
    points, regiontype = readpoints(f)
    name = readline(f)

    Region(name, version, classification, boundingrectangle,
           notes, detectionsettings, points, regiontype)

end

function regions(f::IOStream)

    function _it(chn1)
        readheader(f)
        n = parse(Int, readline(f))
        for i in 1:n
            region = readregion(f)
            put!(chn1, region)
        end
    end

    return Channel(_it, ctype=Region)
end

"""
    regions(filename::AbstractString)

`regions` returns an iterator over regions in the given Echoview EVR
file.

"""
function regions(filename::AbstractString)
    regions([filename])
end

"""
    load(filename::AbstractString)

Loads an Echoview region file.
"""
function load(filename::AbstractString)
    collect(regions([filename]))
end

"""
    regions(filenames::Vector{String})

`regions` returns an iterator over regions in the given Echoview EVR
files.

"""
function regions(filenames::Vector{String})

    function _it(chn1)
        for filename in filenames
            open(filename) do f
                for region in regions(f)
                    put!(chn1, region)
                end
            end
        end
    end

    return Channel(_it, ctype=Region)
end

function makefiletime(date, time)

    CCYY = parse(Int,date[1:4])
    MM = parse(Int, date[5:6])
    DD = parse(Int, date[7:8])

    HH = parse(Int, time[1:2])
    mm = parse(Int, time[3:4])
    SS = parse(Int, time[5:6])

    ssss = div(parse(Int, time[7:10]), 10)

    d = DateTime(CCYY, MM, DD, HH, mm, SS)

    ft = Filetimes.filetime(d)

    ft = ft + ssss *10000
end

function polygon(region::Region)
    p = region.points
    d = p[1:3:end]
    t = p[2:3:end]
    y = parse.(Float64,p[3:3:end])

    x = makefiletime.(d,t)

    return x,y
end

function polygons(ps::Vector{Region})
    polygon.(ps)
end

function polygons(x)
    polygons(collect(regions(x)))
end

end # module

ext: jl


text:
using Shapley
using Test

@testset "Shapley.jl" begin
    # Write your tests here.
end


text:
"""
Accumulator-like data structure for storing frequencies of CJK words (although other tokens
can be stored as well).

You generally don't need to explicitly call this struct's constructor yourself; rather, use
the `charfreq` function.
"""
struct CJKFrequency{S <: AbstractString, C <: Number}
    freq::Accumulator{S, C}
    size::Ref{C}
    function CJKFrequency(frequencies=Dict{String, Int}())
        key_type = eltype(keys(frequencies))
        val_type = eltype(values(frequencies))
        new{key_type, val_type}(counter(frequencies), Ref(sum(values(frequencies))))
    end
end

# mutation
function DataStructures.inc!(cf::CJKFrequency, key, count=1)
    inc!(cf.freq, key, count)
    cf.size[] += count
end
function DataStructures.dec!(cf::CJKFrequency, key, count=1)
    if key in keys(cf.freq)
        cf.size[] -= cf.freq[key] - max(0, cf.freq[key] - count)
        dec!(cf.freq, key, min(count, cf.freq[key]))
    end
end
function DataStructures.reset!(cf::CJKFrequency, key)
    key in keys(cf.freq) && (size[] -= cf.freq[key])
    reset!(cf.freq, key)
end

# collections operations
Base.isempty(cf::CJKFrequency) = isempty(cf.freq)
function Base.empty!(cf::CJKFrequency)
     empty!(cf.freq)
     size = 0
     return cf
end
Base.in(item, cf::CJKFrequency) = in(item, cf.freq)

# dictionaries
Base.haskey(cf::CJKFrequency, key) = haskey(cf.freq, key)
Base.keys(cf::CJKFrequency) = keys(cf.freq)
Base.values(cf::CJKFrequency) = values(cf.freq)

# indexing
Base.getindex(cf::CJKFrequency, i) = getindex(cf.freq, i)
Base.firstindex(cf::CJKFrequency) = firstindex(cf.freq)
Base.lastindex(cf::CJKFrequency) = lastindex(cf.freq)

# iteration
Base.iterate(cf::CJKFrequency) = iterate(cf.freq)
Base.iterate(cf::CJKFrequency, state) = iterate(cf.freq, state)
Base.length(cf::CJKFrequency) = length(cf.freq)
Base.size(cf::CJKFrequency) = cf.size[]

"""
    charfreq(text)
    charfreq(charfreq_type)

Create a character frequency mapping from either text or load it from a default location for
pre-specified character frequency datasets (e.g. `SimplifiedLCMC`, `SimplifiedJunDa`, etc.).

## Examples
When creating a character frequency from text, this method behaves almost exactly like
`DataStructures.counter` except that the return value always has type `CharacterFrequency`
(`Accumulator{String, Int}`).
```julia-repl
julia> text = split("王老师性格内向，沉默寡言，我除在课外活动小组“文学研究会”听过他一次报告，并听-邓知识渊博，是“老师的老师”外，对他一无所知。所以，研读他的作", "");

julia> charfreq(text)
CJKFrequency{SubString{String}, Int64}(Accumulator(除 => 1, 报 => 1, 是 => 1, 知 => 2, 并 => 1, 性 => 1, ， => 6, 言 => 1, 邓 => 1, 外 => 2, 所 => 2, 对 => 1, 动 => 1, 寡 => 1, 。 => 1, 渊 => 1, 学 => 1, - => 1, 听 => 2, 我 => 1, 次 => 1, 一 => 2, 读 => 1, 作 => 1, 格 => 1, “ => 2, 博 => 1, 课 => 1, 老 => 3, 会 => 1, 告 => 1, 无 => 1, 活 => 1, 组 => 1, 内 => 1, 师 => 3, 的 => 2, 小 => 1, 文 => 1, 默 => 1, 究 => 1, 过 => 1, 在 => 1, 以 => 1, ” => 2, 研 => 2, 他 => 3, 向 => 1, 沉 => 1, 王 => 1), Base.RefValue{Int64}(71))
```

See the documentation for individual character frequency dataset structs for examples of the
second case.
"""
function charfreq end
charfreq(tokens) = CJKFrequency(counter(tokens))
charfreq(text::AbstractString) = CJKFrequency(counter(split(text, "")))


text:
using Pkg
using PackageCompiler

function build_julia_image(project_path::String)

    Pkg.activate(joinpath(project_path, "julia_envs", "Xpress"))
    Pkg.instantiate(verbose=true)

    # must point pycall to the python path we want it to use
    ENV["PYTHON"] = joinpath(project_path, "env", "bin", "python")
    ENV["PROJECT_PATH"] = project_path

    # TODO create REopt module to attach to Base for system image
    if Sys.islinux()
        ext = ".so"
    elseif Sys.isapple()
        ext = ".dylib"
    elseif Sys.iswindows()
        ext = ".dll"
    else
        error("Unsupported operating system")
    end

    PackageCompiler.create_sysimage(
        [:AxisArrays, :JuMP, :MathOptInterface, :PyCall, :Xpress, :MutableArithmetics],
        sysimage_path=joinpath(project_path, "julia_envs", "Xpress", "JuliaXpressSysimage" * ext),
        precompile_execution_file=joinpath(project_path, "julia_envs", "Xpress", "precompile.jl")
    )
end


text:
"""
`DeepONet(architecture_branch::Tuple, architecture_trunk::Tuple,
                        act_branch = identity, act_trunk = identity;
                        init_branch = Flux.glorot_uniform,
                        init_trunk = Flux.glorot_uniform,
                        bias_branch=true, bias_trunk=true)`
`DeepONet(branch_net::Flux.Chain, trunk_net::Flux.Chain)`

Create an (unstacked) DeepONet architecture as proposed by Lu et al.
arXiv:1910.03193

The model works as follows:

x --- branch --
               |
                -⊠--u-
               |
y --- trunk ---

Where `x` represents the input function, discretely evaluated at its respective sensors.
So the ipnut is of shape [m] for one instance or [m x b] for a training set.
`y` are the probing locations for the operator to be trained. It has shape [N x n] for
N different variables in the PDE (i.e. spatial and temporal coordinates) with each n distinct evaluation points.
`u` is the solution of the queried instance of the PDE, given by the specific choice of parameters.

Both inputs `x` and `y` are multiplied together via dot product Σᵢ bᵢⱼ tᵢₖ.

You can set up this architecture in two ways:

1. By Specifying the architecture and all its parameters as given above. This always creates
 `Dense` layers for the branch and trunk net and corresponds to the DeepONet proposed by Lu et al.

2. By passing two architectures in the form of two Chain structs directly. Do this if you want more
flexibility and e.g. use an RNN or CNN instead of simple `Dense` layers.

Strictly speaking, DeepONet does not imply either of the branch or trunk net to be a simple
 DNN. Usually though, this is the case which is why it's treated as the default case here.

# Example

Consider a transient 1D advection problem ∂ₜu + u ⋅ ∇u = 0, with an IC u(x,0) = g(x).
We are given several (b = 200) instances of the IC, discretized at 50 points each and want
 to query the solution for 100 different locations and times [0;1].

That makes the branch input of shape [50 x 200] and the trunk input of shape [2 x 100]. So the
 input for the branch net is 50 and 100 for the trunk net.

# Usage

```julia
julia> model = DeepONet((32,64,72), (24,64,72))
DeepONet with
branch net: (Chain(Dense(32, 64), Dense(64, 72)))
Trunk net: (Chain(Dense(24, 64), Dense(64, 72)))

julia> model = DeepONet((32,64,72), (24,64,72), σ, tanh; init_branch=Flux.glorot_normal, bias_trunk=false)
DeepONet with
branch net: (Chain(Dense(32, 64, σ), Dense(64, 72, σ)))
Trunk net: (Chain(Dense(24, 64, tanh; bias=false), Dense(64, 72, tanh; bias=false)))

julia> branch = Chain(Dense(2,128),Dense(128,64),Dense(64,72))
Chain(
  Dense(2, 128),                        # 384 parameters
  Dense(128, 64),                       # 8_256 parameters
  Dense(64, 72),                        # 4_680 parameters
)                   # Total: 6 arrays, 13_320 parameters, 52.406 KiB.

julia> trunk = Chain(Dense(1,24),Dense(24,72))
Chain(
  Dense(1, 24),                         # 48 parameters
  Dense(24, 72),                        # 1_800 parameters
)                   # Total: 4 arrays, 1_848 parameters, 7.469 KiB.

julia> model = DeepONet(branch,trunk)
DeepONet with
branch net: (Chain(Dense(2, 128), Dense(128, 64), Dense(64, 72)))
Trunk net: (Chain(Dense(1, 24), Dense(24, 72)))
```
"""
struct DeepONet{T1, T2}
    branch_net::T1
    trunk_net::T2
end

# Declare the function that assigns Weights and biases to the layer
function DeepONet(architecture_branch::Tuple, architecture_trunk::Tuple,
                        act_branch = identity, act_trunk = identity;
                        init_branch = Flux.glorot_uniform,
                        init_trunk = Flux.glorot_uniform,
                        bias_branch=true, bias_trunk=true)

    @assert architecture_branch[end] == architecture_trunk[end] "Branch and Trunk net must share the same amount of nodes in the last layer. Otherwise Σᵢ bᵢⱼ tᵢₖ won't work."

    # To construct the subnets we use the helper function in subnets.jl
    # Initialize the branch net
    branch_net = construct_subnet(architecture_branch, act_branch;
                                    init=init_branch, bias=bias_branch)
    # Initialize the trunk net
    trunk_net = construct_subnet(architecture_trunk, act_trunk;
                                    init=init_trunk, bias=bias_trunk)

    return DeepONet{typeof(branch_net),typeof(trunk_net)}(branch_net, trunk_net)
end

Flux.@functor DeepONet

#= The actual layer that does stuff
x is the input function, evaluated at m locations (or m x b in case of batches)
y is the array of sensors, i.e. the variables of the output function
with shape (N x n) - N different variables with each n evaluation points =#
function (a::DeepONet)(x::AbstractArray, y::AbstractVecOrMat)
    # Assign the parameters
    branch, trunk = a.branch_net, a.trunk_net

    #= Dot product needs a dim to contract
    However, we perform the transformations by the NNs always in the first dim
    so we need to adjust (i.e. transpose) one of the inputs,
    which we do on the branch input here =#
    return branch(x)' * trunk(y)
end

# Sensors stay the same and shouldn't be batched
(a::DeepONet)(x::AbstractArray, y::AbstractArray) =
  throw(ArgumentError("Sensor locations fed to trunk net can't be batched."))

# Print nicely
function Base.show(io::IO, l::DeepONet)
    print(io, "DeepONet with\nbranch net: (",l.branch_net)
    print(io, ")\n")
    print(io, "Trunk net: (", l.trunk_net)
    print(io, ")\n")
end


text:
PCA = @load PCA pkg="MultivariateStats"

"""
    shuffle_split(X, Y, Z, p_split)

    Creates a stratified partition with p_train of data in the test set.
    Stratification occurs along feature Y.
"""
function shuffle_split(X::DataFrame, Y::Vector, Z::DataFrame, p_split::Float64=0.7)
    train, test = partition(eachindex(Y), p_split, stratify=Y, shuffle=true)
    Xtr, Ytr, Ztr = X[train,:], coerce(Y[train], Finite), Z[train,:]
    Xte, Yte, Zte = X[test,:], coerce(Y[test], Finite), Z[test,:]
    return Xtr, Ytr, Ztr, Xte, Yte, Zte
end

"""
    split_prune_bootstrap(X, Y, Z, train_p, prune_p, bootstrap_samples)

    Runs the initial shuffle-splitting (with train_p proportion in the training set), followed by
    pruning of variables with more than prune_p proportion of data missing, and finally followed by
    backdoor causal bootstrap resampling on the discrete confounder Z.
"""
function split_prune_bootstrap(
    X::DataFrame,
    Y::Vector,
    Z::Union{DataFrame, Nothing},
    train_p::Float64=0.8,
    prune_p::Float64=0.6,
    bootstrap_samples::Int=1000)

    Xₜᵣ, Yₜᵣ, Zₜᵣ, Xₜₑ, Yₜₑ, Zₜₑ = shuffle_split(X, Y, Z, train_p)
    Xₜᵣ = select(Xₜᵣ, prune_missingness(Xₜᵣ, prune_p))
    if bootstrap_samples > 0
        Xₜᵣ, Yₜᵣ, Zₜᵣ = discrete_causal_bootstrapping(Xₜᵣ, Yₜᵣ, Zₜᵣ, bootstrap_samples)
    end
    return Xₜᵣ, Yₜᵣ, Zₜᵣ, Xₜₑ, Yₜₑ, Zₜₑ
end

"""
    classification_analysis(X, Y, Z, label, n_splits, split_p, prune_p, bootstrap_samples)

    Runs a classification analysis
"""
function classification_analysis(
    model::Any,
    X::DataFrame,
    Y::Vector{Float64},
    Z::Union{DataFrame, Nothing},
    label::String="Classification Analysis",
    n_splits::Int=100,
    split_p::Float64=0.5,
    prune_p::Float64=0.6,
    bootstrap_samples::Int=1000)

    res = DataFrame()
    roc_curves = []
    predictions = DataFrame()
    for i ∈ 1:n_splits
        Xₜᵣ,Yₜᵣ,Zₜᵣ,Xₜₑ,Yₜₑ,Zₜₑ = split_prune_bootstrap(X, Y, Z, split_p, prune_p, bootstrap_samples)
        ℳ = model(Xₜᵣ, Yₜᵣ) |> fit!
        Y′ = predict(ℳ, Xₜₑ)

        # Store results
        res = store_classification_results(res, i, Y′, Yₜₑ)
        predictions = store_classification_predictions(predictions, i, label, Zₜₑ, Yₜₑ, Y′)
        fprs, tprs, ts = roc(Y′, coerce(Yₜₑ, OrderedFactor))
        push!(roc_curves, [fprs, tprs])
    end

    return res, predictions, roc_curves
end

"""
    pca_plot(X,U,pc_i,pc_j; palette, alpha, title, dpi)

Performs principal components analysis and plots PCs.
"""
function pca_plot(
    X::DataFrame,
    U::Vector,
    pc_i::Int64=1,
    pc_j::Int64=2;
    palette::Symbol=:seaborn_bright,
    alpha::Float64=0.5,
    title::String="",
    dpi::Int64=350)

    pipe = @pipeline FeatureSelector FillImputer Standardizer OneHotEncoder PCA(maxoutdim=10)
    ℳ = machine(pipe, X) |> fit!
    Z′ = ℳ(X)

    fig = plot(palette=:seaborn_bright,
        title=title,
        xlabel=L"PC_1",
        ylabel=L"PC_2",
        size=(400,400))
    for u ∈ unique(U)
        fig = scatter!(
            Z′[U .== u, pc_i],
            Z′[U .== u, pc_j],
            alpha=0.5, label=u,
            dpi=dpi)
    end
    return fig
end


text:
#  Achdou, Han, Lasry, Lions, Moll (2015)
#  Heterogeneous Agent Models in Continuous Time

##############################################################################
##
## Type
##
##############################################################################

type AiyagariProblem
    π::Float64               # productivity
    γ::Float64               # CRRA utility with parameter gamma
    α::Float64               # Production function F = K^α * L^(1-α)
    δ::Float64               # Capital depreciation
    ρ::Float64               # discount rate
    σ2::Float64
    θ::Float64
    zmean::Float64
    z::Vector{Float64}       # productivity vector
    a::Vector{Float64}       # wealth vector
    invΔ::Float64            # 1/δ in HznB algorithm
end

function AiyagariProblem(;
                         π::Float64 = 0.0,
                         γ::Float64 = 2.0,
                         α::Float64 = 0.35,
                         δ::Float64 = 0.1,
                         ρ::Float64 = 0.05,
                         σ2::Float64 = (0.10)^2,
                         θ::Float64 = 0.3,
                         zmean::Float64 = 1.0,
                         zn::Int = 20,
                         zmin::Float64 = 0.5,
                         zmax::Float64 = 1.5,
                         amin::Float64 = -1.0,
                         amax::Float64 = 50.0,
                         an::Int = 50,
                         invΔ::Float64 = 1e-3
                         )

    a = collect(linspace(amin,amax,an))
    z = collect(linspace(zmin,zmax,zn))
    AiyagariProblem(π, γ, α, δ, ρ, σ2, θ, zmean, z, a, invΔ)
end

##############################################################################
##
## Aiyagari Solution
##
##############################################################################
type AiyagariSolution{T}
    V::Vector{T}       # initial value function
    g::Vector{T}      # distribution
    K::T               # initial aggregate capital.
    r::T
    w::T
end

function AiyagariSolution{T}(ap::AiyagariProblem; K::T = 3.8)
    # initial value function
    π = ap.π ; α = ap.α ; δ = ap.δ ; z = ap.z ; a = ap.a ; γ = ap.γ ; ρ = ap.ρ
    an = length(a)
    zn = length(z)
    r = exp(π) * α * K^(α-1) - δ
    w = exp(π) * (1-α) * K^α
    V = Array(T, an*zn)
    ij = zero(Int)
    for zi in 1:zn, ai in 1:an
        ij += 1
        V[ij] = (w * z[zi] + r * a[ai])^(1-γ)/((1-γ)*ρ)
    end

    # storage arrays
    g = Array(T, an*zn)
    AiyagariSolution(V, g, K, r, w)
end

##############################################################################
##
## Aiyagari Storage
##
##############################################################################

type AiyagariArrays{T}
    C::Base.SparseMatrix.SparseMatrixCSC{T, Int}
    B::Base.SparseMatrix.SparseMatrixCSC{T, Int}
    A::Base.SparseMatrix.SparseMatrixCSC{T, Int}
    b::Vector{T}
    u::Vector{T}
end

function AiyagariArrays{T}(ap::AiyagariProblem, as::AiyagariSolution{T})
    π = ap.π ; α = ap.α ; δ = ap.δ ; z = ap.z ; a = ap.a ; γ = ap.γ ; ρ = ap.ρ; σ2 = ap.σ2 ; θ = ap.θ ; zmean = ap.zmean
    invdz = 1/(z[2] - z[1])
    invda = 1/(a[2] - a[1])
    an = length(a)
    zn = length(z)

    # create 3 sparse matrices of the right form filled with zeros
    A = spdiagm(
        (ones(T, an*zn), ones(T, an*zn-1), ones(T, an*zn-1), ones(T, (zn-1)*an), ones(T, (zn-1)*an)),
        (0, 1, -1, an, -an)
    )
    fill!(nonzeros(A), zero(T))
    B = deepcopy(A)
    C = deepcopy(A)

    # fill up matrix C
    Cvals = nonzeros(C)
    Crows = rowvals(C)
    ij = zero(Int)
    @inbounds for zi in 1:zn, ai in 1:an
        ij += 1
        krange = nzrange(C, ij)
        rows = Crows[krange]

        if zi > 1
            current =  0.5 * σ2 * invdz^2
            index = searchsortedfirst(rows, ij - an)
            Cvals[krange[index]] += current
            index = searchsortedfirst(rows, ij)
            Cvals[krange[index]] -= current
        end

        if zi < zn
            current = θ * (zmean - z[zi]) * invdz + 0.5 * σ2 * invdz^2
            index = searchsortedfirst(rows, ij + an)
            Cvals[krange[index]] += current
            index = searchsortedfirst(rows, ij)
            Cvals[krange[index]] -= current
        end
    end

    # b such that Ag = b in plank
    b = fill(zero(T), an*zn)
    b[1] = .1

    u = similar(b)

    AiyagariArrays(C, B, A, b, u)
end

##############################################################################
##
## Solve hamilton-jacobi-bellman
##
##############################################################################

    # update A and u
function  update_Au!(ap::AiyagariProblem, aa::AiyagariArrays, as::AiyagariSolution)
    γ = ap.γ  ;  a = ap.a ; z = ap.z ; invΔ = ap.invΔ ;  ρ = ap.ρ ;     z = ap.z ; a = ap.a ;
    C = aa.C  ; A = aa.A ; u = aa.u ; B = aa.B ;
    an = length(a)
    zn = length(z)
    invda = (an-1)/(a[end]-a[1])
    invγ = 1/γ
    inv1γ = 1.0/(1.0-γ)
    Cvals = nonzeros(C)
    Avals = nonzeros(A)
    copy!(Avals, Cvals)
    ij = zero(Int)
    @inbounds for zi in 1:zn, ai in 1:an
        ij += 1
        krange = nzrange(A, ij)
        rows = rowvals(A)[krange]
        if ai < an
            ∂V = (as.V[ij+1] - as.V[ij]) * invda
        else
            # state constraint boundary condition
            ∂V = (as.w * z[zi] + as.r * a[end])^(-γ)
        end
        saving = as.w * z[zi] + as.r * a[ai] - ∂V^(-invγ)
        if saving > 0 # case of positive drift
            if ai < an
                current = saving * invda
                index = searchsortedfirst(rows, ij)
                Avals[krange[index]] -= current
                index = searchsortedfirst(rows, ij + 1)
                Avals[krange[index]] += current
            end
        else
            if ai > 1
                ∂V = (as.V[ij] - as.V[ij-1]) * invda
            else
                # state constraint boundary condition
                ∂V =  (as.w * z[zi] + as.r * a[1])^(-γ)
            end
            saving = as.w * z[zi] + as.r * a[ai] - ∂V^(-invγ)
            if saving < 0 # case of negative drift
                if ai > 1
                    current = saving * invda
                    index = searchsortedfirst(rows, ij - 1)
                    Avals[krange[index]] -= current
                    index = searchsortedfirst(rows, ij)
                    Avals[krange[index]] += current
                end
            else
                ∂V = (as.w * z[zi] + as.r * a[ai])^(-γ)
            end
        end

        # update u
        u[ij] = inv1γ * ∂V^(1-invγ) + invΔ * as.V[ij]
    end
end

function update_B!(ap::AiyagariProblem, aa::AiyagariArrays, as::AiyagariSolution)
    γ = ap.γ  ;  a = ap.a ; z = ap.z ; invΔ = ap.invΔ ;  ρ = ap.ρ ;     z = ap.z ; a = ap.a ;
    C = aa.C  ; A = aa.A ; u = aa.u ; B = aa.B ;
    an = length(a)
    zn = length(z)
    Bvals = nonzeros(B)
    Brows = rowvals(B)
    Avals = nonzeros(A)
    ij = zero(Int)
    @inbounds for zi in 1:zn, ai in 1:an
        ij += 1
        for k in nzrange(B, ij)
            # loop over elements in the column ij
            row = Brows[k]
            Bvals[k] = - Avals[k]
            if row == ij
                Bvals[k] += invΔ + ρ
            end
        end
    end
end

function solve_hjb!(ap::AiyagariProblem,
                    aa::AiyagariArrays,
                    as::AiyagariSolution;
                    maxit::Int = 100,
                    crit::Float64 = 1e-6,
                    verbose::Bool = true)
    # update price vectors

    for iter in 1:maxit
        # update newV
        update_Au!(ap, aa, as)
        update_B!(ap, aa, as)
        newV = aa.B' \ aa.u
        # check convergence
        distance = chebyshev(newV, as.V)
        if distance < crit
            if verbose
                println("hjb solved : $(iter) iterations")
            end
            break
        else
            # update V using newV
            (newV, as.V) = (as.V, newV)
        end
    end
end

##############################################################################
##
## Solve fokker-planck equation
##
##############################################################################

# solve fokker-planck
function solve_fp!(ap::AiyagariProblem, aa::AiyagariArrays, as::AiyagariSolution)
    A = aa.A
    # change first row of matrix A
    Arows = rowvals(A)
    Avals = nonzeros(A)
    @inbounds for ij in 1:size(A, 2)
        for k in nzrange(A, ij)
            if Arows[k] == 1
                Avals[k] = zero(Float64)
            end
        end
    end
    Avals[1] = one(Float64)

    # solve system Ag = b
    as.g = A \ aa.b
    scale!(as.g, 1/sum(as.g))
end

##############################################################################
##
## Solve equilibrium by iterating on K
##
##############################################################################
function update_rw!(ap::AiyagariProblem, as::AiyagariSolution)
    as.r = exp(ap.π) * ap.α * as.K^(ap.α-1) - ap.δ
    as.w = exp(ap.π) * (1-ap.α) * as.K^ap.α
end

# compute agregate capital
function sum_capital{T}(a::Vector{Float64}, g::Vector{T})
    numerator = zero(Float64)
    denominator = zero(Float64)
    an = length(a)
    zn = div(length(g), length(a))
    ij = zero(Int)
    @inbounds for zi in 1:zn
        @simd for ai in 1:an
            ij += 1
            current = g[ij]
            numerator += current * a[ai]
            denominator += current
        end
    end
    return numerator / denominator
end

# solve
function solve!(ap::AiyagariProblem, aa::AiyagariArrays, as::AiyagariSolution;
                maxit::Int  = 100,      # maximum number of iterations in the HznB loop
                maxitK::Int = 100,      # maximum number of iterations in the K loop
                crit::Float64 = 1e-6,   # criterion HznB loop
                critK::Float64 = 1e-5,  # criterion K loop
                relax::Float64 = 0.99,  # relaxation parameter
                verbose = true)
    for iter in 1:maxitK
        if verbose
            println("Main loop iteration ",iter)
        end
        update_rw!(ap, as)

        # solve hamilton-jacobi-bellman equation
        solve_hjb!(ap, aa, as, maxit = maxit, crit = crit, verbose = verbose)

        # solve fokker-planck equation
        solve_fp!(ap, aa, as)

        # Update agregate capital
        newK = sum_capital(ap.a, as.g)
        if verbose
            @show newK
        end
        if abs(as.K - newK) < critK
            break
        else
            # relaxation algorithm (to ensure convergence)
            as.K = relax * as.K + (1 - relax) * newK
        end
    end
    return as
end

function solve(ap::AiyagariProblem;
    maxit::Int  = 100,      # maximum number of iterations in the HznB loop
    maxitK::Int = 100,      # maximum number of iterations in the K loop
    crit::Float64 = 1e-6,   # criterion HznB loop
    critK::Float64 = 1e-5,  # criterion K loop
    relax::Float64 = 0.99,  # relaxation parameter
    verbose = true)
    as = AiyagariSolution(ap)
    aa = AiyagariArrays(ap, as)
    solve!(ap, aa, as; maxit = maxit, maxitK = maxitK, crit = crit, critK = critK, relax = relax, verbose = verbose)
end

# as = AiyagariSolution(ap)
# HJBFiniteDifference.solve_hjb!(ap, as)


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

    Copyright © 2019: Stockholm Environment Institute U.S.

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

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

if @isdefined GLPK
    @info "Testing scenario solution with GLPK."

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

        # Test with default outputs
        NemoMod.calculatescenario(dbfile; jumpmodel=Model(optimizer_with_attributes(GLPK.Optimizer, "presolve" => true)), numprocs=1, restrictvars=false, quiet = false)

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

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

        @test isapprox(testqry[1,:val], 3845.15711985937; atol=TOL)
        @test isapprox(testqry[2,:val], 146.552326874185; atol=TOL)
        @test isapprox(testqry[3,:val], 139.573639845721; atol=TOL)
        @test isapprox(testqry[4,:val], 132.927276043543; atol=TOL)
        @test isapprox(testqry[5,:val], 126.597405755756; atol=TOL)
        @test isapprox(testqry[6,:val], 120.568957862625; atol=TOL)
        @test isapprox(testqry[7,:val], 114.827578916785; atol=TOL)
        @test isapprox(testqry[8,:val], 109.359598968367; atol=TOL)
        @test isapprox(testqry[9,:val], 104.151999017492; atol=TOL)
        @test isapprox(testqry[10,:val], 99.1923800166593; atol=TOL)

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

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

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

        @test isapprox(testqry[1,:val], 3845.15711985937; atol=TOL)
        @test isapprox(testqry[2,:val], 146.552326874185; atol=TOL)
        @test isapprox(testqry[3,:val], 139.573639845721; atol=TOL)
        @test isapprox(testqry[4,:val], 132.927276043543; atol=TOL)
        @test isapprox(testqry[5,:val], 126.597405755756; atol=TOL)
        @test isapprox(testqry[6,:val], 120.568957862625; atol=TOL)
        @test isapprox(testqry[7,:val], 114.827578916785; atol=TOL)
        @test isapprox(testqry[8,:val], 109.359598968367; atol=TOL)
        @test isapprox(testqry[9,:val], 104.151999017492; atol=TOL)
        @test isapprox(testqry[10,:val], 99.1923800166593; atol=TOL)

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

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

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

        @test isapprox(testqry[1,:val], 3845.15711985937; atol=TOL)
        @test isapprox(testqry[2,:val], 146.552326874185; atol=TOL)
        @test isapprox(testqry[3,:val], 139.573639845721; atol=TOL)
        @test isapprox(testqry[4,:val], 132.927276043543; atol=TOL)
        @test isapprox(testqry[5,:val], 126.597405755756; atol=TOL)
        @test isapprox(testqry[6,:val], 120.568957862625; atol=TOL)
        @test isapprox(testqry[7,:val], 114.827578916785; atol=TOL)
        @test isapprox(testqry[8,:val], 109.359598968367; atol=TOL)
        @test isapprox(testqry[9,:val], 104.151999017492; atol=TOL)
        @test isapprox(testqry[10,:val], 99.1923800166593; atol=TOL)

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

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

        @test isapprox(testqry[1,:val], 3840.94032304097; atol=TOL)
        @test isapprox(testqry[2,:val], 459.294966842956; atol=TOL)
        @test isapprox(testqry[3,:val], 437.423777945669; atol=TOL)
        @test isapprox(testqry[4,:val], 416.594074233972; atol=TOL)
        @test isapprox(testqry[5,:val], 396.756261175212; atol=TOL)
        @test isapprox(testqry[6,:val], 377.863105881154; atol=TOL)
        @test isapprox(testqry[7,:val], 359.869624648718; atol=TOL)
        @test isapprox(testqry[8,:val], 342.732975855922; atol=TOL)
        @test isapprox(testqry[9,:val], 326.412357958021; atol=TOL)
        @test isapprox(testqry[10,:val], 310.868912340972; atol=TOL)

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

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

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

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

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

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

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

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

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

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

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

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

        @test isapprox(testqry[1,:val], 9786.56610437941; atol=TOL)
        @test isapprox(testqry[2,:val], 239.495268121409; atol=TOL)
        @test isapprox(testqry[3,:val], 228.090669080493; atol=TOL)
        @test isapprox(testqry[4,:val], 217.229187314711; atol=TOL)
        @test isapprox(testqry[5,:val], 206.884940299725; atol=TOL)
        @test isapprox(testqry[6,:val], 197.033276475929; atol=TOL)
        @test isapprox(testqry[7,:val], 187.650739500884; atol=TOL)
        @test isapprox(testqry[8,:val], 178.714990000842; atol=TOL)
        @test isapprox(testqry[9,:val], 170.204752381754; atol=TOL)
        @test isapprox(testqry[10,:val], 162.099764173099; atol=TOL)

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

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

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

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

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

        @test isapprox(testqry[1,:val], 12672.8117352623; atol=TOL)
        @test isapprox(testqry[2,:val], 2510.44571676115; atol=TOL)
        @test isapprox(testqry[3,:val], 1611.02249720726; atol=TOL)

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

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


text:
module UMFPACK

export UmfpackLU,
       decrement,
       decrement!,
       increment,
       increment!

import Base: (\), Ac_ldiv_B, At_ldiv_B, findnz, getindex, show, size

import ..LinAlg: A_ldiv_B!, Ac_ldiv_B!, At_ldiv_B!, Factorization, det, lufact, lufact!

include("umfpack_h.jl")
type MatrixIllConditionedException <: Exception
    message :: String
end

function umferror(status::Int)
     if status==UMFPACK_OK
         return
     elseif status==UMFPACK_WARNING_singular_matrix
         throw(MatrixIllConditionedException("Singular matrix"))
     elseif status==UMFPACK_WARNING_determinant_underflow
         throw(MatrixIllConditionedException("The determinant is nonzero but underflowed"))
     elseif status==UMFPACK_WARNING_determinant_overflow
         throw(MatrixIllConditionedException("The determinant overflowed"))
     elseif status==UMFPACK_ERROR_out_of_memory
         throw(MemoryError())
     elseif status==UMFPACK_ERROR_invalid_Numeric_object
         throw(ArgumentError("Invalid UMFPack numeric object"))
     elseif status==UMFPACK_ERROR_invalid_Symbolic_object
         throw(ArgumentError("Invalid UMFPack symbolic object"))
     elseif status==UMFPACK_ERROR_argument_missing
         throw(ArgumentError("A required argument to UMFPack is missing"))
     elseif status==UMFPACK_ERROR_n_nonpositive
         throw(BoundsError("The number of rows or columns of the matrix must be greater than zero"))
     elseif status==UMFPACK_ERROR_invalid_matrix
         throw(ArgumentError("Invalid matrix"))
     elseif status==UMFPACK_ERROR_different_pattern
         throw(ArgumentError("Pattern of the matrix changed"))
     elseif status==UMFPACK_ERROR_invalid_system
         throw(ArgumentError("Invalid sys argument provided to UMFPack solver"))
     elseif status==UMFPACK_ERROR_invalid_permutation
         throw(ArgumentError("Invalid permutation"))
     elseif status==UMFPACK_ERROR_file_IO
         throw(IOError())
     elseif status==UMFPACK_ERROR_ordering_failed
         error("The ordering method failed")
     elseif status==UMFPACK_ERROR_internal_error
         error("An internal error has occurred, of unknown cause")
     else
         error("Unknown error code: $status")
     end
end

macro isok(A)
    :(umferror($A))
end

function decrement!{T<:Integer}(A::AbstractArray{T})
    for i in 1:length(A) A[i] -= one(T) end
    A
end
decrement{T<:Integer}(A::AbstractArray{T}) = decrement!(copy(A))
function increment!{T<:Integer}(A::AbstractArray{T})
    for i in 1:length(A) A[i] += one(T) end
    A
end
increment{T<:Integer}(A::AbstractArray{T}) = increment!(copy(A))

typealias UMFVTypes Union(Float64,Complex128)
typealias UMFITypes Union(Int32,Int64)

## UMFPACK

# the control and info arrays
const umf_ctrl = Array(Float64, UMFPACK_CONTROL)
ccall((:umfpack_dl_defaults,:libumfpack), Void, (Ptr{Float64},), umf_ctrl)
const umf_info = Array(Float64, UMFPACK_INFO)

function show_umf_ctrl(level::Real)
    old_prt::Float64 = umf_ctrl[1]
    umf_ctrl[1] = float64(level)
    ccall((:umfpack_dl_report_control, :libumfpack), Void, (Ptr{Float64},), umf_ctrl)
    umf_ctrl[1] = old_prt
end
show_umf_ctrl() = show_umf_ctrl(2.)

function show_umf_info(level::Real)
    old_prt::Float64 = umf_ctrl[1]
    umf_ctrl[1] = float64(level)
    ccall((:umfpack_dl_report_info, :libumfpack), Void,
          (Ptr{Float64}, Ptr{Float64}), umf_ctrl, umf_info)
    umf_ctrl[1] = old_prt
end
show_umf_info() = show_umf_info(2.)

## Should this type be immutable?
type UmfpackLU{Tv<:UMFVTypes,Ti<:UMFITypes} <: Factorization{Tv}
    symbolic::Ptr{Void}
    numeric::Ptr{Void}
    m::Int
    n::Int
    colptr::Vector{Ti}                  # 0-based column pointers
    rowval::Vector{Ti}                  # 0-based row indices
    nzval::Vector{Tv}
end

function lufact{Tv<:UMFVTypes,Ti<:UMFITypes}(S::SparseMatrixCSC{Tv,Ti})
    S.m == S.n || error("argument matrix must be square")

    zerobased = S.colptr[1] == 0
    res = UmfpackLU(C_NULL, C_NULL, S.m, S.n,
                    zerobased ? copy(S.colptr) : decrement(S.colptr),
                    zerobased ? copy(S.rowval) : decrement(S.rowval),
                    copy(S.nzval))
    finalizer(res, umfpack_free_symbolic)
    umfpack_numeric!(res)
end

function lufact!{Tv<:UMFVTypes,Ti<:UMFITypes}(S::SparseMatrixCSC{Tv,Ti})
    S.m == S.n || error("argument matrix must be square")

    zerobased = S.colptr[1] == 0
    res = UmfpackLU(C_NULL, C_NULL, S.m, S.n,
                    zerobased ? S.colptr : decrement!(S.colptr),
                    zerobased ? S.rowval : decrement!(S.rowval),
                    S.nzval)
    finalizer(res, umfpack_free_symbolic)
    umfpack_numeric!(res)
end

function show(io::IO, f::UmfpackLU)
    println(io, "UMFPACK LU Factorization of a $(f.m)-by-$(f.n) sparse matrix")
    f.numeric != C_NULL && println(f.numeric)
end

## Wrappers for UMFPACK functions

for (sym_r,sym_c,num_r,num_c,sol_r,sol_c,det_r,det_z,lunz_r,lunz_z,get_num_r,get_num_z,itype) in
    (("umfpack_di_symbolic","umfpack_zi_symbolic",
      "umfpack_di_numeric","umfpack_zi_numeric",
      "umfpack_di_solve","umfpack_zi_solve",
      "umfpack_di_get_determinant","umfpack_zi_get_determinant",
      "umfpack_di_get_lunz","umfpack_zi_get_lunz","umfpack_di_get_numeric","umfpack_zi_get_numeric",:Int32),
     ("umfpack_dl_symbolic","umfpack_zl_symbolic",
      "umfpack_dl_numeric","umfpack_zl_numeric",
      "umfpack_dl_solve","umfpack_zl_solve",
      "umfpack_dl_get_determinant","umfpack_zl_get_determinant",
      "umfpack_dl_get_lunz","umfpack_zl_get_lunz","umfpack_dl_get_numeric","umfpack_zl_get_numeric",:Int64))
    @eval begin
        function umfpack_symbolic!{Tv<:Float64,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
            if U.symbolic != C_NULL return U end
            tmp = Array(Ptr{Void},1)
            @isok ccall(($sym_r, :libumfpack), Ti,
                           (Ti, Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Tv}, Ptr{Void},
                            Ptr{Float64}, Ptr{Float64}),
                           U.m, U.n, U.colptr, U.rowval, U.nzval, tmp,
                           umf_ctrl, umf_info)
            U.symbolic = tmp[1]
            return U
        end
        function umfpack_symbolic!{Tv<:Complex128,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
            if U.symbolic != C_NULL return U end
            tmp = Array(Ptr{Void},1)
            @isok ccall(($sym_c, :libumfpack), Ti,
                           (Ti, Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Void},
                            Ptr{Float64}, Ptr{Float64}),
                           U.m, U.n, U.colptr, U.rowval, real(U.nzval), imag(U.nzval), tmp,
                           umf_ctrl, umf_info)
            U.symbolic = tmp[1]
            return U
        end
        function umfpack_numeric!{Tv<:Float64,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
            if U.numeric != C_NULL return U end
            if U.symbolic == C_NULL umfpack_symbolic!(U) end
            tmp = Array(Ptr{Void}, 1)
            status = ccall(($num_r, :libumfpack), Ti,
                           (Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Void}, Ptr{Void},
                            Ptr{Float64}, Ptr{Float64}),
                           U.colptr, U.rowval, U.nzval, U.symbolic, tmp,
                           umf_ctrl, umf_info)
            status > 0 && throw(MatrixIllConditionedException(""))
            umferror(status)
            U.numeric = tmp[1]
            return U
        end
        function umfpack_numeric!{Tv<:Complex128,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
            if U.numeric != C_NULL return U end
            if U.symbolic == C_NULL umfpack_symbolic!(U) end
            tmp = Array(Ptr{Void}, 1)
            status = ccall(($num_c, :libumfpack), Ti,
                           (Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Void},
                            Ptr{Float64}, Ptr{Float64}),
                           U.colptr, U.rowval, real(U.nzval), imag(U.nzval), U.symbolic, tmp,
                           umf_ctrl, umf_info)
            status > 0 && throw(MatrixIllConditionedException(""))
            umferror(status)
            U.numeric = tmp[1]
            return U
        end
        function solve{Tv<:Float64,Ti<:$itype}(lu::UmfpackLU{Tv,Ti}, b::VecOrMat{Tv}, typ::Integer)
            umfpack_numeric!(lu)
            size(b,1)==lu.m || throw(DimensionMismatch(""))
            x = similar(b)
            joff = 1
            for k = 1:size(b,2)
                @isok ccall(($sol_r, :libumfpack), Ti,
                            (Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64},
                             Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
                            typ, lu.colptr, lu.rowval, lu.nzval, pointer(x,joff), pointer(b,joff), lu.numeric, umf_ctrl, umf_info)
                joff += size(b,1)
            end
            x
        end
        function solve{Tv<:Complex128,Ti<:$itype}(lu::UmfpackLU{Tv,Ti}, b::VecOrMat{Tv}, typ::Integer)
            umfpack_numeric!(lu)
            size(b,1)==lu.m || throw(DimensionMismatch(""))
            x = similar(b)
            n = size(b,1)
            br = Array(Float64, n)
            bi = Array(Float64, n)
            xr = Array(Float64, n)
            xi = Array(Float64, n)
            joff = 0
            for k = 1:size(b,2)
                for j = 1:n
                    bj = b[joff+j]
                    br[j] = real(bj)
                    bi[j] = imag(bj)
                end
                @isok ccall(($sol_c, :libumfpack),
                            Ti,
                            (Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64},
                             Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
                             Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
                            typ, lu.colptr, lu.rowval, real(lu.nzval), imag(lu.nzval),
                            xr, xi, br, bi,
                            lu.numeric, umf_ctrl, umf_info)
                for j = 1:n
                    x[joff+j] = complex(xr[j],xi[j])
                end
                joff += n
            end
            x
        end
        function det{Tv<:Float64,Ti<:$itype}(lu::UmfpackLU{Tv,Ti})
            mx = Array(Tv,1)
            @isok ccall(($det_r,:libumfpack), Ti,
                           (Ptr{Tv},Ptr{Tv},Ptr{Void},Ptr{Float64}),
                           mx, C_NULL, lu.numeric, umf_info)
            mx[1]
        end
        function det{Tv<:Complex128,Ti<:$itype}(lu::UmfpackLU{Tv,Ti})
            mx = Array(Float64,1)
            mz = Array(Float64,1)
            @isok ccall(($det_z,:libumfpack), Ti,
                           (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Void},Ptr{Float64}),
                           mx, mz, C_NULL, lu.numeric, umf_info)
            complex(mx[1], mz[1])
        end
        function umf_lunz{Tv<:Float64,Ti<:$itype}(lu::UmfpackLU{Tv,Ti})
            lnz = Array(Ti, 1)
            unz = Array(Ti, 1)
            n_row = Array(Ti, 1)
            n_col = Array(Ti, 1)
            nz_diag = Array(Ti, 1)
            @isok ccall(($lunz_r,:libumfpack), Ti,
                           (Ptr{Ti},Ptr{Ti},Ptr{Ti},Ptr{Ti},Ptr{Ti},Ptr{Void}),
                           lnz, unz, n_row, n_col, nz_diag, lu.numeric)
            (lnz[1], unz[1], n_row[1], n_col[1], nz_diag[1])
        end
        function umf_lunz{Tv<:Complex128,Ti<:$itype}(lu::UmfpackLU{Tv,Ti})
            lnz = Array(Ti, 1)
            unz = Array(Ti, 1)
            n_row = Array(Ti, 1)
            n_col = Array(Ti, 1)
            nz_diag = Array(Ti, 1)
            @isok ccall(($lunz_z,:libumfpack), Ti,
                           (Ptr{Ti},Ptr{Ti},Ptr{Ti},Ptr{Ti},Ptr{Ti},Ptr{Void}),
                           lnz, unz, n_row, n_col, nz_diag, lu.numeric)
            (lnz[1], unz[1], n_row[1], n_col[1], nz_diag[1])
        end
        function umf_extract{Tv<:Float64,Ti<:$itype}(lu::UmfpackLU{Tv,Ti})
            umfpack_numeric!(lu)        # ensure the numeric decomposition exists
            (lnz,unz,n_row,n_col,nz_diag) = umf_lunz(lu)
            Lp = Array(Ti, n_col + 1)
            Lj = Array(Ti, lnz) # L is returned in CSR (compressed sparse row) format
            Lx = Array(Tv, lnz)
            Up = Array(Ti, n_col + 1)
            Ui = Array(Ti, unz)
            Ux = Array(Tv, unz)
            P  = Array(Ti, n_row)
            Q  = Array(Ti, n_col)
            Rs = Array(Tv, n_row)
            @isok ccall(($get_num_r,:libumfpack), Ti,
                           (Ptr{Ti},Ptr{Ti},Ptr{Tv},
                            Ptr{Ti},Ptr{Ti},Ptr{Tv},
                            Ptr{Ti},Ptr{Ti},Ptr{Void},
                            Ptr{Ti},Ptr{Tv},Ptr{Void}),
                           Lp,Lj,Lx,
                           Up,Ui,Ux,
                           P, Q, C_NULL,
                           &0, Rs, lu.numeric)
            (transpose(SparseMatrixCSC(n_row,n_row,increment!(Lp),increment!(Lj),Lx)),
             SparseMatrixCSC(n_row,n_col,increment!(Up),increment!(Ui),Ux),
             increment!(P), increment!(Q), Rs)
        end
        function umf_extract{Tv<:Complex128,Ti<:$itype}(lu::UmfpackLU{Tv,Ti})
            umfpack_numeric!(lu)        # ensure the numeric decomposition exists
            (lnz,unz,n_row,n_col,nz_diag) = umf_lunz(lu)
            Lp = Array(Ti, n_col + 1)
            Lj = Array(Ti, lnz) # L is returned in CSR (compressed sparse row) format
            Lx = Array(Float64, lnz)
            Lz = Array(Float64, lnz)
            Up = Array(Ti, n_col + 1)
            Ui = Array(Ti, unz)
            Ux = Array(Float64, unz)
            Uz = Array(Float64, unz)
            P  = Array(Ti, n_row)
            Q  = Array(Ti, n_col)
            Rs = Array(Float64, n_row)
            @isok ccall(($get_num_z,:libumfpack), Ti,
                           (Ptr{Ti},Ptr{Ti},Ptr{Float64},Ptr{Float64},
                            Ptr{Ti},Ptr{Ti},Ptr{Float64},Ptr{Float64},
                            Ptr{Ti},Ptr{Ti},Ptr{Void}, Ptr{Void},
                            Ptr{Ti},Ptr{Float64},Ptr{Void}),
                           Lp,Lj,Lx,Lz,
                           Up,Ui,Ux,Uz,
                           P, Q, C_NULL, C_NULL,
                           &0, Rs, lu.numeric)
            (transpose(SparseMatrixCSC(n_row,n_row,increment!(Lp),increment!(Lj),complex(Lx,Lz))),
             SparseMatrixCSC(n_row,n_col,increment!(Up),increment!(Ui),complex(Ux,Uz)),
             increment!(P), increment!(Q), Rs)
        end
    end
end

### Solve with Factorization
A_ldiv_B!{T<:UMFVTypes}(lu::UmfpackLU{T}, b::VecOrMat{T}) = solve(lu, b, UMFPACK_A)
function A_ldiv_B!{Tlu<:Real,Tb<:Complex}(lu::UmfpackLU{Tlu}, b::Vector{Tb})
    r = solve(lu, [convert(Tlu,real(be)) for be in b], UMFPACK_A)
    i = solve(lu, [convert(Tlu,imag(be)) for be in b], UMFPACK_A)
    Tb[r[k]+im*i[k] for k = 1:length(r)]
end
A_ldiv_B!{Tlu<:UMFVTypes,Tb<:Number}(lu::UmfpackLU{Tlu}, b::StridedVecOrMat{Tb}) = A_ldiv_B!(lu, convert(Array{Tlu}, b))
A_ldiv_B!{Tlu<:UMFVTypes,Tb<:Number}(lu::UmfpackLU{Tlu}, b::AbstractVecOrMat{Tb}) = A_ldiv_B!(lu, convert(Array{Tlu}, b))

Ac_ldiv_B!{T<:UMFVTypes}(lu::UmfpackLU{T}, b::VecOrMat{T}) = solve(lu, b, UMFPACK_At)
function Ac_ldiv_B!{Tlu<:Real,Tb<:Complex}(lu::UmfpackLU{Tlu}, b::Vector{Tb})
    r = solve(lu, [convert(Float64,real(be)) for be in b], UMFPACK_At)
    i = solve(lu, [convert(Float64,imag(be)) for be in b], UMFPACK_At)
    Tb[r[k]+im*i[k] for k = 1:length(r)]
end
Ac_ldiv_B!{Tlu<:UMFVTypes,Tb<:Number}(lu::UmfpackLU{Tlu}, b::StridedVecOrMat{Tb}) = Ac_ldiv_B!(lu, convert(Array{Tlu}, b))
Ac_ldiv_B!{Tlu<:UMFVTypes,Tb<:Number}(lu::UmfpackLU{Tlu}, b::AbstractVecOrMat{Tb}) = Ac_ldiv_B!(lu, convert(Array{Tlu}, b))

At_ldiv_B!{T<:UMFVTypes}(lu::UmfpackLU{T}, b::VecOrMat{T}) = solve(lu, b, UMFPACK_Aat)
function At_ldiv_B!{Tlu<:Real,Tb<:Complex}(lu::UmfpackLU{Tlu}, b::Vector{Tb})
    r = solve(lu, [convert(Float64,real(be)) for be in b], UMFPACK_Aat)
    i = solve(lu, [convert(Float64,imag(be)) for be in b], UMFPACK_Aat)
    Tb[r[k]+im*i[k] for k = 1:length(r)]
end
At_ldiv_B!{Tlu<:UMFVTypes,Tb<:Number}(lu::UmfpackLU{Tlu}, b::StridedVecOrMat{Tb}) = At_ldiv_B!(lu, convert(Array{Tlu}, b))
At_ldiv_B!{Tlu<:UMFVTypes,Tb<:Number}(lu::UmfpackLU{Tlu}, b::AbstractVecOrMat{Tb}) = At_ldiv_B!(lu, convert(Array{Tlu}, b))

function getindex(lu::UmfpackLU, d::Symbol)
    L,U,p,q,Rs = umf_extract(lu)
    d == :L ? L :
    (d == :U ? U :
     (d == :p ? p :
      (d == :q ? q :
       (d == :Rs ? Rs :
        (d == :(:) ? (L,U,p,q,Rs) :
         throw(KeyError(d)))))))
end

## The C functions called by these Julia functions do not depend on
## the numeric and index types, even though the umfpack names indicate
## they do.  The umfpack_free_* functions can be called on C_NULL without harm.
function umfpack_free_symbolic(symb::Ptr{Void})
    tmp = [symb]
    ccall((:umfpack_dl_free_symbolic, :libumfpack), Void, (Ptr{Void},), tmp)
end
show_umf_info() = show_umf_info(2.)

function umfpack_free_symbolic(lu::UmfpackLU)
    if lu.symbolic == C_NULL return lu end
    umfpack_free_numeric(lu)
    umfpack_free_symbolic(lu.symbolic)
    lu.symbolic = C_NULL
    return lu
end

function umfpack_free_numeric(num::Ptr{Void})
    tmp = [num]
    ccall((:umfpack_dl_free_numeric, :libumfpack), Void, (Ptr{Void},), tmp)
end

function umfpack_free_numeric(lu::UmfpackLU)
    if lu.numeric == C_NULL return lu end
    umfpack_free_numeric(lu.numeric)
    lu.numeric = C_NULL
    return lu
end

function umfpack_report_symbolic(symb::Ptr{Void}, level::Real)
    old_prl::Float64 = umf_ctrl[UMFPACK_PRL]
    umf_ctrl[UMFPACK_PRL] = float64(level)
    @isok ccall((:umfpack_dl_report_symbolic, :libumfpack), Int,
                   (Ptr{Void}, Ptr{Float64}), symb, umf_ctrl)
    umf_ctrl[UMFPACK_PRL] = old_prl
end

umfpack_report_symbolic(symb::Ptr{Void}) = umfpack_report_symbolic(symb, 4.)

function umfpack_report_symbolic(lu::UmfpackLU, level::Real)
    umfpack_report_symbolic(umfpack_symbolic!(lu).symbolic, level)
end

umfpack_report_symbolic(lu::UmfpackLU) = umfpack_report_symbolic(lu.symbolic,4.)
function umfpack_report_numeric(num::Ptr{Void}, level::Real)
    old_prl::Float64 = umf_ctrl[UMFPACK_PRL]
    umf_ctrl[UMFPACK_PRL] = float64(level)
    @isok ccall((:umfpack_dl_report_numeric, :libumfpack), Int,
                   (Ptr{Void}, Ptr{Float64}), num, umf_ctrl)
    umf_ctrl[UMFPACK_PRL] = old_prl
end

umfpack_report_numeric(num::Ptr{Void}) = umfpack_report_numeric(num, 4.)
function umfpack_report_numeric(lu::UmfpackLU, level::Real)
    umfpack_report_numeric(umfpack_numeric!(lu).numeric, level)
end

umfpack_report_numeric(lu::UmfpackLU) = umfpack_report_numeric(lu,4.)

end # UMFPACK module

max_stars_repo_name: vchuravy/julia

max_issues_repo_name: vchuravy/julia

max_forks_repo_name: vchuravy/julia


text:
import RowEchelon

#=
Homework 2
Section 1.4
Problem 37
=#

# Create the Matrix to be evaluated
A = [ 7 2 -5 8 ;
	 -5 -3 4 -9 ;
	 6 10 -2 7 ;
	 -7 9 2 15 ]

# Store the columns for future reference
columns = 4

# Row reduce the Matrix
A_rref = RowEchelon.rref(A)

println("Original Matrix, A")
println(A)
println("Row Reduced Echelon Matrix:")
println(A_rref)

# Determine whether or not matrix has a pivot position in each row
hasAllZeroRow = false
index = 1
currentZeros = 0
while index <= columns
	global index
	global hasAllZeroRow
	global currentZeros

	# Parse through each entry of the rows, searching for an all 0 row.
	for entry in A_rref[index, :]
		global hasAllZeroRow
		global currentZeros

		if entry == 0
			# If here then:
			# Current entry is a zero, increment zero counter
			currentZeros += 1
		end
	end
	if currentZeros == columns
		# If here then:
		# A row with all zeroes was found
		hasAllZeroRow = true
	end
	currentZeros = 0
	index += 1
end

if hasAllZeroRow == true
	# If here then:
	# The matrix does not span R4
	println("The Matrix does not span R4")
else
	# If here then:
	# The matrix does span all elements in R4
	println("The Matrix does span R4")
end

println()
println()
exit(0)


text:
using TSne
using DataFrames

doc = """Use t-SNE to generate a PDF called myplot.pdf from an input CSV file. Default assumption is to have no header and no labels. If these are available in the CSV these must be given as arguments.

Usage:
  demo-csv.jl <filename>
  demo-csv.jl [--labelcol=<col>] <filename>
  demo-csv.jl haveheader [--labelcol=<col>] <filename>
  demo-csv.jl [--labelcolname=<colname>] <filename>
  demo-csv.jl haveheader [--labelcolname=<colname>] <filename>

Options:
  -h --help     Show this screen.
  --version     Show version.
  --filename=   Path to CSV file
  --noheader    The CSV file does not have a header row
  --nolabel    The CSV file does not have a label column

"""

"""
Normalize `A` columns, so that the mean and standard deviation
of each column are 0 and 1, resp.
"""
function rescale(A, dim::Integer=1)
    res = A .- mean(A, dim)
    res ./= map!(x -> x > 0.0 ? x : 1.0, std(A, dim))
    res
end

using DocOpt

arguments = docopt(doc, version=v"2.0.0")
dump(arguments)

if nothing==arguments["--labelcol"]
    lblcol = -1
else
    lblcol = parse(Int64, arguments["--labelcol"])
end

df = readtable(arguments["<filename>"], header = nothing!=arguments["haveheader"])

if nothing!=arguments["--labelcolname"]
    lblcol = findfirst(x -> x==symbol(arguments["--labelcolname"]), names(df))
end

println("Data is $df")
labels = lblcol > 0 ? df[:, lblcol] : nothing

dataset = df[filter(x -> x!=lblcol, 1:ncol(df)), :]
data = convert(Matrix{Float64}, dataset)
# Normalize the data, this should be done if there are large scale differences in the dataset
X = rescale(data)

# Run t-SNE
Y = tsne(X, 2, 50, 1000, 20.0)

using Gadfly
theplot = plot(x=Y[:,1], y=Y[:,2], color=labels)
draw(PDF("myplot.pdf", 8inch, 6inch), theplot)


text:
__precompile__(true)

module MarkovChainIntExperiment

using Reproduce
using JLD2
using Resampling
import Resampling.ExpUtils: algorithm_settings!, build_algorithm_dict
using Resampling.ExpUtils.MarkovChainUtil
using Random
using Statistics
import Flux: Descent
import ProgressMeter
import StatsBase

function exp_settings(as::ArgParseSettings = ArgParseSettings(exc_handler=Reproduce.ArgParse.debug_handler))
    @add_arg_table as begin
        # Basic Experiment Arguments
        "--exp_loc"
        arg_type=String
        required=true
        "--run"
        arg_type=Int64
        required=true
        "--seed"
        help = "Starting seed (initialize seed + run)"
        default = 0
        arg_type = Int64
        "--numiter"
        help = "Number of training iterations"
        default = 1000
        arg_type = Int64
        "--buffersize"
        help = "Size of buffer"
        default = 15000
        arg_type = Int64
        "--batchsize"
        help = "Size of batch"
        default = 16
        arg_type = Int64
        "--onpolicy"
        action=:store_true
        "--working"
        action=:store_true
    end

    @add_arg_table as begin
        "--alphas"
        arg_type=Float64
        nargs='+'
        required=true
    end

    @add_arg_table as begin
        "--beta"
        help="Interpolation between policies"
        arg_type=Float64
        required=true
    end

    # algorithm_settings!(as)
    env_settings!(as)

    return as
end

function mave(model::Resampling.TabularLayer{Array{Float64, 1}}, truth, args...)
    return mean(abs.(model.W .- truth))
end

function mave(model::Resampling.TabularLayer{Array{Float64, 2}}, truth, target_policy)
    return mean(abs.(model.W'*target_policy  .- truth))
end

function main_experiment(args::Vector{String})

    as = exp_settings()
    parsed = parse_args(args, as)
    save_loc=""
    if !parsed["working"]
        create_info!(parsed, parsed["exp_loc"]; filter_keys=["verbose", "working", "exp_loc"])
        save_loc = Reproduce.get_save_dir(parsed)
        if isfile(joinpath(save_loc, "results.jld2"))
            return
        end
    end

    run = parsed["run"]
    chain_size = parsed["chainsize"]
    num_iterations = parsed["numiter"]
    buffer_size = parsed["buffersize"]
    batch_size = parsed["batchsize"]
    α_arr = parsed["alphas"]
    β = parsed["beta"]

    μ, π = get_policies(parsed)
    if parsed["onpolicy"]
        μ = π
    end

    truth = Resampling.DynamicProgramming(MarkovChain(chain_size), π, 0.9)

    lu = BatchTD()
    value_arr = [Tabular(chain_size) for α in α_arr]
    local_error = zeros(length(α_arr), num_iterations)

    rng = MersenneTwister(parsed["seed"] + run)
    WER = get_experience_interpolation(buffer_size, μ, π, β*μ + (1.0-β)*π; rng=rng, chain_size=chain_size)

    for iter in 1:num_iterations

        samp_wer = sample(WER, batch_size; rng=rng)
        arg_wer = (samp_wer[:ρ], samp_wer[:s_t], samp_wer[:s_tp1],
                   samp_wer[:r], samp_wer[:γ_tp1], samp_wer[:terminal])

        for (α_idx, α) in enumerate(α_arr)
            opt = Descent(α)
            update!(value_arr[α_idx], opt, lu, arg_wer...;)
            local_error[α_idx, iter] = mave(value_arr[α_idx], truth, π)
        end
    end

    # println(value_function)

    if parsed["working"]
        return local_error
    else
        results = local_error
        @save joinpath(save_loc, "results.jld2") results
    end

end

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    main_experiment(ARGS)
    return 0
end

end


text:
export
    ZmqServer

mutable struct ZmqServer{T <: DFG.AbstractDFG}
    fg::T
    config::Dict{String, Any}
    isServerActive::Bool
    binding::String
end


text:
mutable struct Graph
    adjMatrix:: Array{Int64,2}
    pheromoneMatrix:: Array{Float64,2}
    foodIndex:: Int64
    nodesCount :: Int64
end

function generateGraph(nodeCount:: Int64, adjMatrix:: Array{Int64,2}, foodIndex::Int64)
    pheromoneMatrix = ones(Float64,nodeCount,nodeCount)
    graph = Graph(adjMatrix,pheromoneMatrix,foodIndex,nodeCount)
    return graph
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
#Residential Consumer Balancing Code
function annual_energy_balance_nmw(consumer::Residential, pvsys::PVSystem; print_output=false)

    results = Dict{Int64,Any}()

    d_d = Truncated(Normal(0.0, 0.5), -0.10, 0.10)
    d_s = Truncated(Normal(-0.2, 0.5), -0.10, 0.15)
    days_p_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    m_h = 0
    count = 0

    # Initial conditions
    carry_over = 0.0
    global_allowance = 0.0
    global_withdrawal = 0.0
    global_generation = 0.0
    energy_loss = 0.0

    #Loop for every Month of the year
    for (ix,m) in enumerate(days_p_month)

        #Get solar Data with houtly resolution for the month
        solar_month_base = pvsys.time_series[24*(m_h) + 1:24*(m_h+m)]
        solar_month = [(s+rand(d_s, 1)[1])*(s > 0.1) for s in solar_month_base]*pvsys.capacity
        m_h = m

        #Simulate month on an hourly basis
        #initial conditions.
        grid_energy = 0.0
        withdrawn_energy = 0.0
        consumer_energy = 0.0
        PV_energy = 0.0
        injection_grid = 0.0

        #Loop for every day of the month
        for i in 1:m
            var = rand(d_d, 24).+1.0
            daily_p = (consumer.load_curve.*var)*consumer.peak_power[ix]
            solar_day = solar_month[(((m-1)*24)+1):(((m-1)*24)+24)]

            #Loop for every hour of the day
            for t in 1:24

               consumer_energy += daily_p[t]

               PV_energy +=  solar_day[t]

               balance = daily_p[t] - solar_day[t]

               injection_grid += max(0.0, -1*balance)

               grid_energy += max(0.0,balance)

            end #end of hourly loop

        end #end of month loop

           global_generation += PV_energy #GenAcum
           global_allowance = 0.49*PV_energy
           global_withdrawal += min(global_allowance,injection_grid)
           energy_loss = max(0.0, injection_grid - global_allowance)
           carry_over = 0.0

           results[ix] = Dict("consumer_energy" => consumer_energy,
                               "PV_energy" => PV_energy,
                                "injection_grid" => injection_grid,
                                "withdrawn_energy" => min(global_allowance,injection_grid),
                            "grid_energy" => grid_energy,
                            "carry_over" => carry_over,
                            "global_generation" => global_generation,
                            "global_withdrawal" => global_withdrawal,
                            "global_allowance" => global_allowance,
                            "energy_loss" => energy_loss
                    )
    end

    if print_output
        print_residential(results)
    end

    return results

end

#Commercial and Industrial Consumer Balancing Code
function annual_energy_balance_nmw(consumer::CommIndus, pvsys::PVSystem; print_output=false)

    results = Dict{Int64,Any}()

    d_d = Truncated(Normal(0.0, 0.5), -0.10, 0.10)
    d_s = Truncated(Normal(-0.2, 0.5), -0.10, 0.15)
    days_p_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    m_h = 0
    count = 0

    # Initial conditions
    peak_power = 0.0
    peak_demand = 0.0
    carry_over = 0.0
    global_allowance = 0.0
    global_withdrawal = 0.0
    global_generation = 0.0
    energy_loss = 0.0

    #Loop for every Month of the year
    for (ix,m) in enumerate(days_p_month)

        #Get solar Data with houtly resolution for the month
        solar_month_base = pvsys.time_series[24*(m_h) + 1:24*(m_h+m)]
        solar_month = [(s+rand(d_s, 1)[1])*(s > 0.1) for s in solar_month_base]*pvsys.capacity
        m_h = m

        #Simulate month on an hourly basis
        #initial conditions.
        grid_energy = 0.0
        withdrawn_energy = 0.0
        consumer_energy = 0.0
        PV_energy = 0.0
        injection_grid = 0.0

        #Loop for every day of the month
        for i in 1:m
            var = rand(d_d, 24).+1.0
            daily_p = (consumer.load_curve.*var)*consumer.peak_power[ix]
            solar_day = solar_month[(((m-1)*24)+1):(((m-1)*24)+24)]

                #Loop for every hour of the day
                for t in 1:24

                   consumer_energy += daily_p[t]

                   PV_energy +=  solar_day[t]

                   balance = daily_p[t] - solar_day[t]

                   peak_power = max(peak_power, daily_p[t])

                   peak_demand = max(peak_demand, balance)

                   injection_grid += max(0.0, -1*balance)

                   grid_energy += max(0.0,balance)

                end #end of hourly loop

            end #end of month loop

           global_generation += PV_energy #GenAcum
           global_allowance = 0.49*PV_energy
           global_withdrawal += min(global_allowance,injection_grid)
           energy_loss = max(0.0, injection_grid - global_allowance)
           carry_over = 0.0

           results[ix] = Dict("consumer_energy" => consumer_energy,
                               "PV_energy" => PV_energy,
                                "injection_grid" => injection_grid,
                                "withdrawn_energy" => min(global_allowance,injection_grid),
                            "grid_energy" => grid_energy,
                            "peak_demand" => peak_demand,
                            "peak_power" => peak_power,
                            "carry_over" => carry_over,
                            "global_generation" => global_generation,
                            "global_withdrawal" => global_withdrawal,
                            "global_allowance" => global_allowance,
                            "energy_loss" => energy_loss
                    )

        peak_power = 0.0
        peak_demand = 0.0

    end

    if print_output
        print_commercial(results)
    end

    return results

end

#TMT Consumer Balancing Code
function annual_energy_balance_nmw(consumer::TMT, pvsys::PVSystem; print_output=false)

    results = Dict{Int64,Any}()

    d_d = Truncated(Normal(0.0, 0.5), -0.10, 0.10)
    d_s = Truncated(Normal(-0.2, 0.5), -0.10, 0.15)
    days_p_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    m_h = 0
    count = 0

    carry_over = 0.0
    global_allowance = 0.0
    global_withdrawal = 0.0
    global_generation = 0.0
    energy_loss = 0.0

    #Loop for every Month of the year
    for (ix,m) in enumerate(days_p_month)

        # Initial monthly conditions
        peak = [10, 11, 12, 17, 18, 19]
        peak_power_peak = 0.0
        peak_demand_peak = 0.0

        valley = [6, 7, 8, 9, 13, 14, 15, 16]
        peak_power_valley = 0.0
        peak_demand_valley = 0.0

        night = [20, 21, 22, 23, 24, 1, 2, 3, 4, 5]
        peak_power_night = 0.0
        peak_demand_night = 0.0

        #Get solar Data with houtly resolution for the month
        solar_month_base = pvsys.time_series[24*(m_h) + 1:24*(m_h+m)]
        solar_month = [(s+rand(d_s, 1)[1])*(s > 0.1) for s in solar_month_base]*pvsys.capacity
        m_h = m

        #Simulate month on an hourly basis
        #initial conditions.
        withdrawn_energy = 0.0
        consumer_energy_peak = 0.0
        consumer_energy_valley = 0.0
        consumer_energy_night = 0.0
        grid_energy_peak = 0.0
        grid_energy_valley = 0.0
        grid_energy_night = 0.0
        PV_energy = 0.0
        injection_grid = 0.0

        #Loop for every day of the month
        for i in 1:m
            var = rand(d_d, 24).+1.0
            daily_p = (consumer.load_curve.*var)*consumer.peak_power[ix]
            solar_day = solar_month[(((m-1)*24)+1):(((m-1)*24)+24)]

            #Loop for every hour of the day
            for t in 1:24

               (t in peak) ? consumer_energy_peak += daily_p[t] : true
               (t in valley) ? consumer_energy_valley += daily_p[t] : true
               (t in night) ? consumer_energy_night += daily_p[t] : true

               PV_energy +=  solar_day[t]

               #println(t,"--------")
               #println("PV"," ", solar_day[t])
               #println("---------")
               balance = daily_p[t] - solar_day[t]

               injection_grid += max(0.0, -1*balance)

               #Peak demand corresponds to the peak as measured by the distribution company

               #(t in peak) ? println("peak d"," ",  peak_demand_peak, " ", balance, " ", max(0.0, peak_demand_peak, max(0.0,balance)) ) : true
               (t in peak) ? peak_demand_peak = max(0.0, peak_demand_peak, max(0.0,balance)) : true

               #(t in valley) ? println("valley d", " ", peak_demand_valley," ",  balance, " ", max(0.0, peak_demand_valley, max(0.0,balance)) ) : true
               (t in valley) ? peak_demand_valley = max(0.0, peak_demand_valley, max(0.0,balance)) : true

               #(t in night) ? println("night d", " ", peak_demand_night, " ", balance, " ", max(0.0, peak_demand_night, max(0.0,balance)) ) : true
               (t in night) ? peak_demand_night = max(0.0, peak_demand_night, max(0.0,balance)) : true

               #Peak power corresponds to the peak of the consumer load

               #(t in peak) ? println("peak p", " ", peak_power_peak, " ", daily_p[t], " ", max(0.0, peak_power_peak, daily_p[t]) ) : true
               (t in peak) ? peak_power_peak = max(0.0, peak_power_peak, daily_p[t]) : true

               #(t in valley) ? println("valley p", " ", peak_power_valley, " ", daily_p[t], " ", max(0.0, peak_power_valley, daily_p[t])) : true
               (t in valley) ? peak_power_valley = max(0.0, peak_power_valley, daily_p[t]) : true

               #(t in night) ? println("night p", " ", peak_power_night, " ", daily_p[t], " ", max(0.0, peak_power_night, daily_p[t]) ) : true
               (t in night) ? peak_power_night = max(0.0, peak_power_night, daily_p[t]) : true

               available_energy = carry_over + max(0.0, -1*balance)

               withdrawal = 0.0

               if balance >= 0.0

                    withdrawal = available_energy - balance

                    if withdrawal > 0.0 #There is enough stored to meet energy demand in t

                        #Book Keeping

                        withdrawn_energy += balance

                        (t in peak) ? grid_energy_peak += 0.0 : true
                        (t in valley) ? grid_energy_valley += 0.0 : true
                        (t in night) ? grid_energy_night += 0.0 : true

                        carry_over = withdrawal

                    elseif withdrawal <= 0.0 #There is not enough stored to meet energy demand in t

                        withdrawn_energy += available_energy

                        carry_over = 0.0

                        (t in peak) ? grid_energy_peak -= withdrawal : true
                        (t in valley) ? grid_energy_valley -= withdrawal : true
                        (t in night) ? grid_energy_night -= withdrawal : true

                    end

               elseif balance < 0.0

                    carry_over  = available_energy

               end

            end #end of hourly loop

        end #end of month loop

           global_generation += PV_energy #GenAcum
           global_allowance = 0.49*PV_energy
           global_withdrawal += min(global_allowance,injection_grid)
           energy_loss = max(0.0, injection_grid - global_allowance)
           carry_over = 0.0

           results[ix] = Dict("consumer_energy_peak" => consumer_energy_peak,
                              "consumer_energy_valley" => consumer_energy_valley,
                              "consumer_energy_night" => consumer_energy_night,
                              "PV_energy" => PV_energy,
                              "injection_grid" => injection_grid,
                              "withdrawn_energy" => withdrawn_energy,
                              "grid_energy_peak" => grid_energy_peak,
                              "grid_energy_valley" => grid_energy_valley,
                              "grid_energy_night" => grid_energy_night,
                              "peak_demand_peak" => peak_demand_peak,
                              "peak_demand_valley" => peak_demand_valley,
                              "peak_demand_night" => peak_demand_night,
                              "peak_power_peak" => peak_power_peak,
                              "peak_power_valley" => peak_power_valley,
                              "peak_power_night" => peak_power_night,
                              "carry_over" => carry_over,
                              "global_generation" => global_generation,
                              "global_withdrawal" => global_withdrawal,
                              "global_allowance" => global_allowance,
                        )

    end

    if print_output
        print_TMT(results)
    end

    return results

end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
# # ======================
# # Example usage
# using CSV
# using Plots
#
# file_name = "sts_for_VP_test.csv"
# csv_matrix = CSV.read(file_name)[:,2:end]
# almost_not_csv_matrix = Matrix(csv_matrix)
#
# file_name2 = "spikes.csv"
# spikes = load_csv_file_to_array(file_name2)
#
# file_name3 = "th_chunks.csv"
# th_chunks = load_csv_file_to_array(file_name3)
#
# sts = generate_spike_matrix(spikes; th_chunks=th_chunks)
# VPd = [get_selfdist(s; n_chan=32, cost=60., dt=0.01) for s in sts]
#
#
# plot_set = Any[]
#     for matrix in VPd
#         push!(plot_set, heatmap(matrix, color=:Reds_9, colorbar=false, yflip = true))
#     end
#
# plot(plot_set[1], plot_set[2], plot_set[3],
#     plot_set[4], plot_set[5], plot_set[6],
#     plot_set[7], plot_set[8], plot_set[9],
#     layout=(1,9), size=(9*1200,1100), legend=false, colorbar=false)

"""
    get_selfdist(st_inp; n_chan=32, cost=60., dt=0.01)

Method for computing pair-wise spike distances from a range of spike trains.
Function copied from Mikolaj SETCOmodel

Inputs:
    st_inp: [2 x N] array with spike times and indices of neurons.
    N - number of spikes generated, 1st row - index of neuron generating given spikes, 2nd row - spike time.
    n_chan - number of neurons (default: 32)
    cost - cost parameter for VP spike distance, in ms (default: 60 ms)
    dt - simulation timestep, in ms (default: 0.01 ms -> 100 kHz)
Output:
    pc - [n_chan x n_chan] matrix containing pairwise VP spikes distances for each pair of neurons.
"""
function get_selfdist(st_inp; n_chan=32, cost=60., dt=0.01)

    unique_ids = unique(st_inp[1,:])
    sts_new = Any[]
    for id = unique_ids
        push!(sts_new, st_inp[2,findall(x->x==id, st_inp[1,:])])
    end

    # sts = [st_inp[0,st_inp[1,:]==i] for i in 1:n_chan]
    pc = zeros(n_chan, n_chan)

    for i = 1:n_chan, j = 1:n_chan
        pc[i,j] = spkd(sts_new[i], sts_new[j], dt/(cost))
    end
    return pc
end

# TODO Add test with MATLAB code run for some spike train and compare with
# results from this
"""
    spkd(s1, s2, cost)

Fast implementation of victor-purpura spike distance (faster than neo & elephant python packages)
Direct Python port of http://www-users.med.cornell.edu/~jdvicto/pubalgor.htmlself.
The below code was tested against the original implementation and yielded exact results.
All credits go to the authors of the original code.

Code was translated from Frotran to Matlab, from Matlab to Python, from
Python to Julia. It was veryfied with MATLAB code.

Input:
    s1,2: pair of vectors of spike times
    cost: cost parameter for computing Victor-Purpura spike distance.
    (Note: the above need to have the same units!)
Output:
    d: VP spike distance.
"""
function spkd(s1, s2, cost)
    nspi=length(s1)
    nspj=length(s2)

    # Why not to use this?
    if cost==0
       return d=abs(nspi-nspj)
    elseif cost==Inf
       return d=nspi+nspj
   end

    scr=zeros(nspi+1, nspj+1)

    # initialize margins with cost of adding a spike
    scr[:,1]=0:nspi
    scr[1,:]=0:nspj

    for i = 2:nspi+1, j = 2:nspj+1
        component1 = scr[i-1,j]+1
        component2 = scr[i,j-1]+1
        component3 = scr[i-1,j-1]+cost*abs(s1[i-1]-s2[j-1])
        scr[i,j] = min(component1, component2, component3)
	end

    d=scr[end,end]

    return d
end

"""
	generate_spike_matrix(spikes; th_chunks=[[0,0]])

Generates matrix of the form [2xN] array with spike times and indices of
neurons. N - number of spikes generated, 1st row - index of neuron generating
given spikes, 2nd row - spike time.

Resulting matrix is time sorted. Resulting matrix may be splint into fragments
by setting 'th_chunks' to a list of fragments in a way such that
'th_chunks[k,1]' is a starting index and 'th_chunks[k,2]' is an ending index of
k'th fragment.
"""
function generate_spike_matrix(spikes; th_chunks=[[0,0]], val_range=20:52)

    if th_chunks == [[0,0]]
        th_chunks[1,1] = 1
        th_chunks[1,2] = size(spikes,1)
    end

    spike_train_simplified = Any[]
    for i = 1:size(th_chunks,1)
        #1 Get a chunk of spike trains of interest
        syllable_spikes = spikes[th_chunks[i, 1]:th_chunks[i, 2], val_range]
        # find all spikes
        all_spikes = findall(x->x==1,syllable_spikes)
        # convert to [2xN] matrix
        total_spikes = length(all_spikes)
        sorted_spikes = zeros(Int, 2,total_spikes)
        for k = 1:total_spikes
            sorted_spikes[1,k] = all_spikes[k][2]
            sorted_spikes[2,k] = all_spikes[k][1]
        end

        push!(spike_train_simplified, sorted_spikes[:,sortperm(sorted_spikes[2,:])])
    end

    return spike_train_simplified
end


text:
using Documenter, IteratorAsArrayDisguises

makedocs(
    modules = [IteratorAsArrayDisguises],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Gustavo Goretkin",
    sitename = "IteratorAsArrayDisguises.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/goretkin/IteratorAsArrayDisguises.jl.git",
)


text:
function _autoapi(ex::Expr, options::AbstractArray{Any}, @nospecialize(__source__))
    except = undef
    var = nothing
    rc = nothing
    cast = nothing
    castexc = true
    for opt_ex in options
        @switch opt_ex begin
        @case :(except($v))
            except = v
        @case :(mangle($v))
            var = v
        @case :borrow2new
            rc = :(Py_XINCREF(ret))
        @case :(cast($f))
            cast = f
        @case :nocastexc
            castexc = false
        @case _
            Base.@warn "@autoapi: unknown option $(opt_ex)."
        end
    end

    sym, anns, returntype, argtypes, args =
        @switch ex begin

        @case :($sym($(argtypes...))::$returntype) && let except′ = undef end ||
              :($sym($(argtypes...))::$returntype != $except′)
            except === undef && (except = except′;)
            narg = length(argtypes)
            args = [Symbol("_arg$(i)") for i = 1:narg]
            anns = [:($(args[i]) :: $(argtypes[i])) for i = 1:narg]
            sym, anns, returntype, Expr(:tuple, argtypes...), args
        end
    var === nothing && (var = sym;)
    r = @q begin
        $__source__
        function $var(apis, $(anns...))::$returntype
            $__source__
            ret = ccall(apis.$sym, $returntype, $argtypes, $(args...))
            $rc
            return ret
        end
    end
    val = gensym(:val)
    if except !== undef
        push!(r.args, :($DIO.DIO_ExceptCode(::typeof($var)) = $except))
    end
    if nothing !== cast
        push!(r.args, :($DIO.DIO_HasCast(::typeof($var)) = true))
        push!(r.args, :($DIO.DIO_Cast(::typeof($var), $val::$returntype)::PyPtr = $cast($val)))
        push!(r.args, :($DIO.DIO_CastExc(::typeof($var)) = $castexc))
    end
    return r
end

"""
create an API from a C symbol
    ```
        @autoapi PyObject_Call(PyPtr, PyPtr)::PyPtr except(Py_NULL)
    ```
    expands to
    ```
    function PyObject_Call(apis, x::PyPtr, y::PyPtr)::PyPtr
        ccall(apis.PyObject_Call, PyPtr, (PyPtr, PyPtr), x, y)
    end

    DIO.DIO_ExceptCode(::typeof(PyObject_Call)) = Py_NULL
    ```
"""
macro autoapi(ex::Expr, @nospecialize(options...))
    esc(_autoapi(ex, Any[options...], __source__))
end

"""
Handling Python exception, the result may not be boxed.
"""
macro callexc(f, args...)
    tmp = gensym("callexc")
    r = @q begin
        $__source__
        $tmp = $f($(args...))
        if DIO_ExceptCode($f) === $tmp
            return Py_NULL
        else
            $tmp
        end
    end
    esc(r)
end


text:
function receive_select(f::Function, ch1::Channel, ch2::Channel)

    while true

        result = @select begin
            ch1 |> value => begin
                f(ch1, value)
            end
            ch2 |> value => begin
                f(ch2, value)
            end
        end

        if result isa Break
            break
        end
    end

end


text:
# See: ttps://julialang.github.io/PackageCompiler.jl/dev/examples/plots/

import Pkg; Pkg.activate(".")
using PackageCompiler

create_sysimage(:Plots,
                sysimage_path="sys_plots.so",
                precompile_execution_file="test.jl")

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
module hitbtc

using HTTP
using JSON
using Nettle
using Base64
using Dates
using Printf
using DataStructures

global const BASE_URL = "https://api.hitbtc.com"

publickey() = ENV["HITBTC-PUBLIC-KEY"]
secretkey() = ENV["HITBTC-SECRET-KEY"]

function public_request(
    method::String,
    request_path::String,
    query::OrderedDict
    )
    url = join([BASE_URL, request_path], "/")
    pairs = join(string.(keys(query),"=", values(query)), "&")
    curl = url * "?" * pairs
    response = HTTP.request(method, curl, query)
    json = JSON.parse(String(response.body))
    return json
end

function details()
    method = "GET"
    path = "api/2/public/symbol"

    result = public_request(method, path, OrderedDict())
    form =
        x -> (
            symbol = x["id"],
            base_currency = x["baseCurrency"],
            quote_currency = x["quoteCurrency"],
   )

    return map(form, result)
end

function stats_24hr(symbol::String = "ALL")
    method = "GET"
    path = "api/2/public/ticker"

    format = "yyyy-mm-ddTHH:MM:SS.sss"
    form =
        x -> (
            symbol = x["symbol"],
            first_price = x["open"],
            last_price = x["last"],
            high_price = x["high"],
            low_price = x["low"],
            base_volume = x["volume"],
            quote_volume = x["volumeQuote"],
            time = DateTime(x["timestamp"][1:length(x["timestamp"])-1], format),
        )
    if !isequal(symbol, "ALL")
        path = "api/2/public/ticker/$symbol"

        x = public_request(method, path, OrderedDict())
        return (
            symbol = x["symbol"],
            first_price = x["open"],
            last_price = x["last"],
            high_price = x["high"],
            low_price = x["low"],
            base_volume = x["volume"],
            quote_volume = x["volumeQuote"],
            time = DateTime(x["timestamp"][1:length(x["timestamp"])-1], format),
        )
    end
    result = public_request(method, path, OrderedDict())
    return map(form, result)
end
stats_24hr()

function order_book(symbol::String)

end  # module hitbtc