text:
# 3-D plot commands

# surface plots
function surf(x::Coord,y::Coord,Z::Coord;
              handle::Handle = gnuplot_state.current,
              args...)

    # process arguments
    PA = PlotArgs()
    for k in keys(args)
        # substitute synonyms
        key = k
        for s in syn
            key ∈ s && (key = s[1]; break)
        end
        # store arguments
        for f in fieldnames(PlotArgs)
            f == key && (setfield!(PA, f, string(args[k])); break)
        end
    end

    # validation and defaults
    valid_3Dplotstyle(PA.plotstyle)
    valid_numeric(PA.linewidth)
    valid_linestyle(PA.linestyle)
    valid_pointtype(PA.pointtype)
    valid_numeric(PA.pointsize)
    valid_axis(PA.axis)
    valid_range(PA.xrange)
    valid_range(PA.yrange)
    valid_range(PA.zrange)

    if isa(Z, Matrix)
        if length(x) != Base.size(Z)[1] || length(y) != Base.size(Z)[2]
            throw(DimensionMismatch("invalid coordinates in surface plot."))
        end
    elseif isa(Z, Vector)
        if length(x) != length(y) || length(x) != length(Z)
            throw(DimensionMismatch("all coordinates must have the same dimension"))
        end
    else
        throw(DimensionMismatch("Z must be a matrix or a vector."))
    end

    term = config[:term][:terminal]
    PA.font == "" && (PA.font = TerminalDefaults[term][:font])
    PA.size == "" && (PA.size = TerminalDefaults[term][:size])
    PA.linewidth == "" && (PA.linewidth = "1")
    PA.background == "" && (PA.background = TerminalDefaults[term][:background])

    handle = figure(handle, redraw = false)
    clearfigure(handle)

    # create curve
    tc = TermConf(font       = PA.font,
                  size       = PA.size,
                  background = PA.background)

    ac = AxesConf(title      = PA.title,
                  xlabel     = PA.xlabel,
                  ylabel     = PA.ylabel,
                  zlabel     = PA.zlabel,
                  keyoptions = PA.keyoptions,
                  xrange     = PA.xrange,
                  yrange     = PA.yrange,
                  zrange     = PA.zrange,
                  xzeroaxis  = PA.xzeroaxis,
                  yzeroaxis  = PA.yzeroaxis,
                  zzeroaxis  = PA.zzeroaxis,
                  palette    = PA.palette)

    cc = CurveConf(plotstyle = PA.plotstyle,
                   legend    = PA.legend,
                   linecolor = PA.linecolor,
                   pointtype = PA.pointtype,
                   linewidth = PA.linewidth,
                   pointsize = PA.pointsize)

    c = Curve(x    = x,
              y    = y,
              Z    = Z,
              conf = cc)

    # push curve we just created to current figure
    push_figure!(handle,tc,ac,c,PA.gpcom)

    return gnuplot_state.figs[findfigure(handle)]
end
surf(x::Coord,y::Coord,f::Function;args...) = surf(x,y,meshgrid(x,y,f);args...)
surf(Z::Matrix;args...) = surf(1:size(Z)[2],1:size(Z)[1],Z;args...)

function surf!(x::Coord,y::Coord,Z::Coord;
               handle::Handle = gnuplot_state.current,
               args...)

    # process arguments
    PA = PlotArgs()
    for k in keys(args)
        for f in fieldnames(PlotArgs)
            k == f && setfield!(PA, f, string(args[k]))
        end
    end
    # validation and defaults
    valid_3Dplotstyle(PA.plotstyle)
    valid_linestyle(PA.linestyle)
    valid_pointtype(PA.pointtype)
    valid_numeric(PA.pointsize)
    valid_pointtype(PA.pointtype)
    valid_numeric(PA.linewidth)

    if isa(Z, Matrix)
        if length(x) != Base.size(Z)[1] || length(y) != Base.size(Z)[2]
            throw(DimensionMismatch("invalid coordinates in surface plot."))
        end
    elseif isa(Z, Vector)
        if length(x) != length(y) || length(x) != length(Z)
            throw(DimensionMismatch("all coordinates must have the same dimension"))
        end
    else
        throw(DimensionMismatch("Z must be a matrix or a vector."))
    end

    handles = gethandles()
    handle ∈ handles || error("Cannot append curve to non-existing handle")
    handle = figure(handle, redraw = false)

    # configure curve
    cc = CurveConf(plotstyle = PA.plotstyle,
                   legend    = PA.legend,
                   linecolor = PA.linecolor,
                   pointtype = PA.pointtype,
                   linewidth = PA.linewidth,
                   pointsize = PA.pointsize)

    c = Curve(x    = x,
              y    = y,
              Z    = Z,
              conf = cc)

    push_figure!(handle,c)

    return gnuplot_state.figs[findfigure(handle)]
end
surf!(x::Coord,y::Coord,f::Function;args...) = surf!(x,y,meshgrid(x,y,f);args...)
surf!(Z::Matrix;args...) = surf!(1:size(Z)[2],1:size(Z)[1],Z;args...)

function contour(x::Coord,y::Coord,Z::Coord;labels=true,args...)
    gp = """unset key
            set view map
            set contour base
            unset surface
            set cntrlabel font ",7"
            set cntrparam levels 10"""
    p = surf(x,y,Z;gpcom=gp,args...)
    labels && (p = surf!(x,y,Z;plotstyle="labels"))
    return p
end
contour(x::Coord,y::Coord,f::Function;args...) = contour(x,y,meshgrid(x,y,f);args...)
contour(Z::Matrix;args...) = contour(1:size(Z)[2],1:size(Z)[1],Z;args...)

scatter3(x::Coord,y::Coord,Z::Coord;args...) = surf(x,y,Z,plotstyle="points";args...)
scatter3!(x::Coord,y::Coord,Z::Coord;args...) = surf!(x,y,Z,plotstyle="points";args...)


text:
@testset "Events/Comments" begin
    evaluate(s) = collect(E.events(L.State(IOBuffer(s))))

    E = Events
    L = Events.Lexical

    CommentStart = E.CommentStart
    CommentEnd   = E.CommentEnd

    DC = E.DataContent
    ME = E.MarkupError

    @testset "Events/Comments (Positive)" begin
        @test evaluate("<!---->") == [ CommentStart(L.Location("a buffer", -1)),
                                       DC("", L.Location("a buffer", -1)),
                                       CommentEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<!--x-->") == [ CommentStart(L.Location("a buffer", -1)),
                                        DC("x", L.Location("a buffer", -1)),
                                        CommentEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<!--xxx-->") == [ CommentStart(L.Location("a buffer", -1)),
                                          DC("xxx", L.Location("a buffer", -1)),
                                          CommentEnd(L.Location("a buffer", -1)) ]

        @test (evaluate("<!--x--><!--y-->") == [ CommentStart(L.Location("a buffer", -1)),
                                                 DC("x", L.Location("a buffer", -1)),
                                                 CommentEnd(L.Location("a buffer", -1)),
                                                 CommentStart(L.Location("a buffer", -1)),
                                                 DC("y", L.Location("a buffer", -1)),
                                                 CommentEnd(L.Location("a buffer", -1)) ])

        @test (evaluate("<!-- Hello, World! -->") == [ CommentStart(L.Location("a buffer", -1)),
                                                       DC(" Hello, World! ", L.Location("a buffer", -1)),
                                                       CommentEnd(L.Location("a buffer", -1)) ])

        @test (evaluate("<!--é-->") == [ CommentStart(L.Location("a buffer", -1)),
                                         DC("é", L.Location("a buffer", -1)),
                                         CommentEnd(L.Location("a buffer", -1)) ])

        # This is from the XML specification, § 2.5.
        #
        @test (evaluate("<!-- declarations for <head> & <body> -->")
               == [ CommentStart(L.Location("a buffer", -1)),
                    DC(" declarations for <head> & <body> ", L.Location("a buffer", -1)),
                    CommentEnd(L.Location("a buffer", -1)) ])

        @test (evaluate("<!-- &e; -->") == [ CommentStart(L.Location("a buffer", -1)),
                                             DC(" &e; ", L.Location("a buffer", -1)),
                                             CommentEnd(L.Location("a buffer", -1)) ])
        @test (evaluate("<!-- <?target value?> -->") == [ CommentStart(L.Location("a buffer", -1)),
                                                          DC(" <?target value?> ", L.Location("a buffer", -1)),
                                                          CommentEnd(L.Location("a buffer", -1)) ])
    end

    @testset "Events/Comments (Negative ... no terminator)" begin
        @test (evaluate("<!--") == [ CommentStart(L.Location("a buffer", -1)),
                                     DC("", L.Location("a buffer", -1)),
                                     ME("ERROR: Expecting '-->' to end a comment.", L.Location("a buffer", -1)),
                                     CommentEnd(true, L.Location("a buffer", -1)) ])
        @test (evaluate("<!--x") == [ CommentStart(L.Location("a buffer", -1)),
                                      DC("x", L.Location("a buffer", -1)),
                                      ME("ERROR: Expecting '-->' to end a comment.", L.Location("a buffer", -1)),
                                      CommentEnd(true, L.Location("a buffer", -1)) ])
        @test (evaluate("<!--é") == [ CommentStart(L.Location("a buffer", -1)),
                                      DC("é", L.Location("a buffer", -1)),
                                      ME("ERROR: Expecting '-->' to end a comment.", L.Location("a buffer", -1)),
                                      CommentEnd(true, L.Location("a buffer", -1)) ])
        @test (evaluate("<!-- declarations for <head> & <body> ")
               == [ CommentStart(L.Location("a buffer", -1)),
                    DC(" declarations for <head> & <body> ", L.Location("a buffer", -1)),
                    ME("ERROR: Expecting '-->' to end a comment.", L.Location("a buffer", -1)),
                    CommentEnd(true, L.Location("a buffer", -1)) ])

        @test (evaluate("<!--x-") == [ CommentStart(L.Location("a buffer", -1)),
                                       DC("x-", L.Location("a buffer", -1)),
                                       ME("ERROR: Expecting '-->' to end a comment.", L.Location("a buffer", -1)),
                                       CommentEnd(true, L.Location("a buffer", -1)) ])
        @test (evaluate("<!--x--") == [ CommentStart(L.Location("a buffer", -1)),
                                        DC("x--", L.Location("a buffer", -1)),
                                        ME("ERROR: Expecting '-->' to end a comment.", L.Location("a buffer", -1)),
                                        CommentEnd(true, L.Location("a buffer", -1)) ])
    end

    @testset "Events/Comments (Negative ... nested --)" begin
        @test (evaluate("<!-- -- -->")
               == [ CommentStart(L.Location("a buffer", -1)),
                    DC(" ", L.Location("a buffer", -1)),
                    ME("ERROR: '--' is not allowed inside a comment.", L.Location("a buffer", -1)),
                    CommentEnd(true, L.Location("a buffer", -1)),
                    DC(" ", true, L.Location("a buffer", -1)),
                    DC("--", false, L.Location("a buffer", -1)),
                    DC(">", false, L.Location("a buffer", -1)) ])
        @test (evaluate("<!------>")
               == [ CommentStart(L.Location("a buffer", -1)),
                    DC("", L.Location("a buffer", -1)),
                    ME("ERROR: '--' is not allowed inside a comment.", L.Location("a buffer", -1)),
                    CommentEnd(true, L.Location("a buffer", -1)),
                    DC("--", false, L.Location("a buffer", -1)),
                    DC(">", false, L.Location("a buffer", -1)) ])

        # This is from the XML specification, § 2.5.
        #
        @test (evaluate("<!-- B+, B, or B--->")
               == [ CommentStart(L.Location("a buffer", -1)),
                    DC(" B+, B, or B", L.Location("a buffer", -1)),
                    ME("ERROR: '--' is not allowed inside a comment.", L.Location("a buffer", -1)),
                    CommentEnd(true, L.Location("a buffer", -1)),
                    DC("-", false, L.Location("a buffer", -1)),
                    DC(">", false, L.Location("a buffer", -1)) ])
    end
end

ext: jl

max_stars_repo_name: jlegare/NativeXML

max_issues_repo_name: jlegare/NativeXML

max_forks_repo_name: jlegare/NativeXML


text:
function numlattices(rcutmax, pbc, a, b, c)

    d = (length(c) < 3) ? 2 : 3

    m = 0
    n = 0
    p = 0

    if pbc[1] == 1
        m = Int64(round(rcutmax/a[1]))
    end
    if pbc[2] == 1
        n = Int64(round(rcutmax/b[2]))
    end
    if d==3
        if pbc[3] == 1
            p = Int64(round(rcutmax/c[3]))
        end
    end

    return m, n, p
end

function referenceorigins(m, n, p, d)

xref = zeros(d, (2*m+1)*(2*n+1)*(2*p+1))
for i = 1:(2*p+1) # z
    for j = 1:(2*n+1) # y
        for k = 1:(2*m+1) # x
            t = k + (2*m+1)*(j-1) + (2*m+1)*(2*n+1)*(i-1)
            xref[1, t] = k-m-1;
            xref[2, t] = j-n-1;
            if d==3
                xref[3, t] = i-p-1;
            end
        end
    end
end

return xref

end

function latticeorigins(a, b, c, m, n, p)

    d = (length(c) < 3) ? 2 : 3

    xo = referenceorigins(m, n, p, d)

    if d==2
        B2R, R2B = domainmaps(a, b);
    else
        B2R, R2B = domainmaps(a, b, c);
    end

    latorigins = R2B*xo;

    return latorigins
end

function latticecoords(x, rcutmax, pbc, a, b, c)

    m, n, p = numlattices(rcutmax, pbc, a, b, c)

    latorigins = latticeorigins(a, b, c, m, n, p)
    ind = (m+1) + (2*m+1)*(n) + (2*m+1)*(2*n+1)*(p)

    d = size(x,1)
    nx = size(x,2)
    nl = size(latorigins,2)  # number of lattices

    y = zeros(d, nx*nl)
    for j = 1:nx
        for k = 1:d
            y[k,j] = x[k,j]
        end
    end
    q = nx;
    for i = 1:nl
        if i != ind
            for j = 1:nx
                q = q + 1
                for k = 1:d
                    y[k,q] = latorigins[k,i] + x[k,j]
                end
            end
        end
    end

    alist = zeros(Int64,nx*nl)
    for i = 1:nl
        for j = 1:nx
            alist[j + nx*(i-1)] = j
        end
    end

    # get
    neighi, neighj, neighnum = neighborlist(y, rcutmax);

    neighi = neighi[neighi .<= nx];
    ne = length(neighi);

    # list of neighbors for each atom i in x
    neighlist = neighj[1:ne];

    # number of neighbors for each atom i in x
    neighnum = neighnum[1:nx];
    neighnum = [0; cumsum(neighnum)];

    return y, alist, neighlist, neighnum
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# DBSCAN Clustering
#
#   References:
#
#       Martin Ester, Hans-peter Kriegel, Jörg S, and Xiaowei Xu
#       A density-based algorithm for discovering clusters
#       in large spatial databases with noise. 1996.
#

"""
    DbscanResult <: ClusteringResult

The output of [`dbscan`](@ref) function (distance matrix-based implementation).

# Fields
 - `seeds::Vector{Int}`: indices of cluster starting points
 - `assignments::Vector{Int}`: vector of clusters indices, where each point was assigned to
 - `counts::Vector{Int}`: cluster sizes (number of assigned points)
"""
mutable struct DbscanResult <: ClusteringResult
    seeds::Vector{Int}          # starting points of clusters, size (k,)
    assignments::Vector{Int}    # assignments, size (n,)
    counts::Vector{Int}         # number of points in each cluster, size (k,)
end

"""
    DbscanCluster

DBSCAN cluster returned by [`dbscan`](@ref) function (point coordinates-based
implementation)

# Fields
 * `size::Int`: number of points in a cluster (core + boundary)
 * `core_indices::Vector{Int}`: indices of points in the cluster *core*
 * `boundary_indices::Vector{Int}`: indices of points on the cluster *boundary*
"""
struct DbscanCluster
    size::Int                      # number of points in cluster
    core_indices::Vector{Int}      # core points indices
    boundary_indices::Vector{Int}  # boundary points indices
end

## main algorithm

"""
    dbscan(D::AbstractMatrix, eps::Real, minpts::Int) -> DbscanResult

Perform DBSCAN algorithm using the distance matrix `D`.

# Arguments
The following options control which points would be considered
*density reachable*:
  - `eps::Real`: the radius of a point neighborhood
  - `minpts::Int`: the minimum number of neighboring points (including itself)
     to qualify a point as a density point.
"""
function dbscan(D::AbstractMatrix{T}, eps::Real, minpts::Int) where T<:Real
    # check arguments
    n = size(D, 1)
    size(D, 2) == n || throw(ArgumentError("D must be a square matrix ($(size(D)) given)."))
    n >= 2 || throw(ArgumentError("At least two data points are required ($n given)."))
    eps > 0 || throw(ArgumentError("eps must be a positive value ($eps given)."))
    minpts >= 1 || throw(ArgumentError("minpts must be positive integer ($minpts given)."))

    # invoke core algorithm
    _dbscan(D, convert(T, eps), minpts, 1:n)
end

function _dbscan(D::AbstractMatrix{T}, eps::T, minpts::Int, visitseq::AbstractVector{Int}) where T<:Real
    n = size(D, 1)

    # prepare
    seeds = Int[]
    counts = Int[]
    assignments = zeros(Int, n)
    visited = zeros(Bool, n)
    k = 0

    # main loop
    for p in visitseq
        if assignments[p] == 0 && !visited[p]
            visited[p] = true
            nbs = _dbs_region_query(D, p, eps)
            if length(nbs) >= minpts
                k += 1
                cnt = _dbs_expand_cluster!(D, k, p, nbs, eps, minpts, assignments, visited)
                push!(seeds, p)
                push!(counts, cnt)
            end
        end
    end

    # make output
    return DbscanResult(seeds, assignments, counts)
end

## key steps

function _dbs_region_query(D::AbstractMatrix{T}, p::Int, eps::T) where T<:Real
    n = size(D,1)
    nbs = Int[]
    dists = view(D,:,p)
    for i = 1:n
        @inbounds if dists[i] < eps
            push!(nbs, i)
        end
    end
    return nbs::Vector{Int}
end

function _dbs_expand_cluster!(D::AbstractMatrix{T},        # distance matrix
                              k::Int,                      # the index of current cluster
                              p::Int,                      # the index of seeding point
                              nbs::Vector{Int},            # eps-neighborhood of p
                              eps::T,                      # radius of neighborhood
                              minpts::Int,                 # minimum number of neighbors of a density point
                              assignments::Vector{Int},    # assignment vector
                              visited::Vector{Bool}) where T<:Real       # visited indicators
    assignments[p] = k
    cnt = 1
    while !isempty(nbs)
        q = popfirst!(nbs)
        if !visited[q]
            visited[q] = true
            qnbs = _dbs_region_query(D, q, eps)
            if length(qnbs) >= minpts
                for x in qnbs
                    if assignments[x] == 0
                        push!(nbs, x)
                    end
                end
            end
        end
        if assignments[q] == 0
            assignments[q] = k
            cnt += 1
        end
    end
    return cnt
end

"""
    dbscan(points::AbstractMatrix, radius::Real,
           [leafsize], [min_neighbors], [min_cluster_size]) -> Vector{DbscanCluster}

Cluster `points` using the DBSCAN (density-based spatial clustering of
applications with noise) algorithm.

# Arguments
 - `points`: the ``d×n`` matrix of points. `points[:, j]` is a
   ``d``-dimensional coordinates of ``j``-th point
 - `radius::Real`: query radius

Optional keyword arguments to control the algorithm:
 - `leafsize::Int` (defaults to 20): the number of points binned in each
   leaf node in the `KDTree`
 - `min_neighbors::Int` (defaults to 1): the minimum number of a *core* point
   neighbors
 - `min_cluster_size::Int` (defaults to 1): the minimum number of points in
   a valid cluster

# Example
``` julia
points = randn(3, 10000)
# DBSCAN clustering, clusters with less than 20 points will be discarded:
clusters = dbscan(points, 0.05, min_neighbors = 3, min_cluster_size = 20)
```
"""
function dbscan(points::AbstractMatrix, radius::Real; leafsize::Int = 20, kwargs ...)
    kdtree = KDTree(points; leafsize=leafsize)
    return _dbscan(kdtree, points, radius; kwargs ...)
end

# An implementation of DBSCAN algorithm that keeps track of both the core and boundary points
function _dbscan(kdtree::KDTree, points::AbstractMatrix, radius::Real;
                 min_neighbors::Int = 1, min_cluster_size::Int = 1)
    dim, num_points = size(points)
    num_points <= dim && throw(ArgumentError("points has $dim rows and $num_points columns. Must be a D x N matric with D < N"))
    0 <= radius || throw(ArgumentError("radius $radius must be ≥ 0"))
    1 <= min_neighbors || throw(ArgumentError("min_neighbors $min_neighbors must be ≥ 1"))
    1 <= min_cluster_size || throw(ArgumentError("min_cluster_size $min_cluster_size must be ≥ 1"))

    clusters = Vector{DbscanCluster}()
    visited = falses(num_points)
    cluster_selection = falses(num_points)
    core_selection = falses(num_points)
    to_explore = Vector{Int}()
    adj_list = Vector{Int}()
    for i = 1:num_points
        visited[i] && continue
        push!(to_explore, i) # start a new cluster
        fill!(core_selection, false)
        fill!(cluster_selection, false)
        while !isempty(to_explore)
            current_index = popfirst!(to_explore)
            visited[current_index] && continue
            visited[current_index] = true
            append!(adj_list, inrange(kdtree, points[:, current_index], radius))
            cluster_selection[adj_list] .= true
            # if a point doesn't have enough neighbors it is not a 'core' point and its neighbors are not added to the to_explore list
            if (length(adj_list) - 1) < min_neighbors
                empty!(adj_list)
                continue # query returns the query point as well as the neighbors
            end
            core_selection[current_index] = true
            update_exploration_list!(adj_list, to_explore, visited)
        end
        cluster_size = sum(cluster_selection)
        min_cluster_size <= cluster_size && accept_cluster!(clusters, core_selection, cluster_selection, cluster_size)
    end
    return clusters
end

"""
    update_exploration_list!(adj_list, exploration_list, visited) -> adj_list

Update the queue for expanding the cluster.

# Arguments
- `adj_list::Vector{Int}`: indices of the neighboring points to move to queue
- `exploration_list::Vector{Int}`: the indices that will be explored in the future
- `visited::BitVector`: a flag indicating whether a point has been explored already
"""
function update_exploration_list!(adj_list::Vector{T}, exploration_list::Vector{T},
                                  visited::BitVector) where T <: Integer
    for j in adj_list
        visited[j] && continue
        push!(exploration_list, j)
    end
    empty!(adj_list)
end

"""
    accept_cluster!(clusters, core_selection, cluster_selection) -> clusters

Accept cluster and update the clusters list.

# Arguments
- `clusters::Vector{DbscanCluster}`: a list of the accepted clusters
- `core_selection::Vector{Bool}`: selection of the core points of the cluster
- `cluster_selection::Vector{Bool}`: selection of all the cluster points
"""
function accept_cluster!(clusters::Vector{DbscanCluster}, core_selection::BitVector,
                         cluster_selection::BitVector, cluster_size::Int)
    core_idx = findall(core_selection) # index list of the core members
    boundary_selection = cluster_selection .& (~).(core_selection) #TODO change to .~ core_selection
                                                                            # when dropping 0.5
    boundary_idx = findall(boundary_selection) # index list of the boundary members
    push!(clusters, DbscanCluster(cluster_size, core_idx, boundary_idx))
end


text:
#=
ndgrid.jl

I know that using ndgrid is not the "julian way" but when porting code
from matlab sometimes it is easier just to use it instead of refactoring
to do some other way.

2019-06-18, Jeff Fessler, University of Michigan
=#

using Test: @test

"""
`(xx,yy) = ndgrid(x::AbstractVector{<:Any}, y::AbstractVector{<:Any})`
todo - improve?
"""
function ndgrid(
		x::AbstractVector{<:Any},
		y::AbstractVector{<:Any})

	tmp = Iterators.product(x, y)
	return [p[1] for p in tmp], [p[2] for p in tmp]
end

"""
`(xx,yy,zz) = ndgrid(x::AbstractVector{<:Any}, y::..., z::...)`
todo - improve?
"""
function ndgrid(
		x::AbstractVector{<:Any},
		y::AbstractVector{<:Any},
		z::AbstractVector{<:Any})

	tmp = Iterators.product(x, y, z)
	return [p[1] for p in tmp], [p[2] for p in tmp], [p[3] for p in tmp]
end

"""
`ndgrid(:test)`
self test
"""
function ndgrid(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))

	@test ndgrid(1:3, 2:4) == ([1 1 1; 2 2 2; 3 3 3], [2 3 4; 2 3 4; 2 3 4])
	@test size(ndgrid(1:3, 4:5, 6:9)[3]) == (3,2,4)
	true
end


text:
using ReinforcementLearning, ReinforcementLearningEnvironments, Flux, BenchmarkTools, CuArrays

b = CircularArrayBuffer{Float64}(84, 84, 10)

function show_msg(desc, x)
    println(repeat("=", 50))
    println(desc)
    println(repeat("=", 50))
    display(x)
    println(repeat("=", 50))
    println()
end

show_msg("push a frame into buffer", @benchmark(push!($b, $(rand(84, 84)))))
show_msg("view a frame", @benchmark(select_frame($b, $(rand(1:10)))))
show_msg("view frames", @benchmark(select_frame($b, $(rand(1:10, 32)))))
show_msg("send frames to gpu", @benchmark(gpu(select_frame($b, $(rand(1:10, 32))))))
show_msg("consecutive view a batch", @benchmark(consecutive_view($b, $(rand(1:5, 32)), 4)))
show_msg("get an element", @benchmark(getindex($b, $(rand(1:84)), $(rand(1:84)), $(rand(1:10)))))

m = Dense(84*84, 512)

show_msg("view a frame and do matrix mul", @benchmark($m(vec(select_frame($b, $(rand(1:10)))))))
show_msg("view a frame and do matrix mul (use CartesianIndex)", @benchmark($m(vec(view($b, :, :, $(rand(1:10)))))))
show_msg("view frames and do matrix mul", @benchmark($m(reshape(select_frame($b, $(rand(1:10, 32))), :, 32))))
show_msg("get a frame and do matrix mul", @benchmark($m(vec(getindex($b, :, :, $(rand(1:10)))))))
show_msg("do matrix mul", @benchmark($m(reshape($b, :, 10))))
show_msg("do matrix mul(for comparison)", @benchmark($m(reshape($(rand(84, 84, 10)), :, 10))))

#=
==================================================
push a frame into buffer
==================================================
BenchmarkTools.Trial:
  memory estimate:  128 bytes
  allocs estimate:  4
  --------------
  minimum time:     1.429 μs (0.00% GC)
  median time:      1.449 μs (0.00% GC)
  mean time:        1.465 μs (0.00% GC)
  maximum time:     7.480 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     10
==================================================

==================================================
view a frame
==================================================
BenchmarkTools.Trial:
  memory estimate:  64 bytes
  allocs estimate:  1
  --------------
  minimum time:     8.779 ns (0.00% GC)
  median time:      11.001 ns (0.00% GC)
  mean time:        19.588 ns (21.90% GC)
  maximum time:     4.445 μs (99.28% GC)
  --------------
  samples:          10000
  evals/sample:     999
==================================================

==================================================
view frames
==================================================
BenchmarkTools.Trial:
  memory estimate:  496 bytes
  allocs estimate:  6
  --------------
  minimum time:     101.667 ns (0.00% GC)
  median time:      104.926 ns (0.00% GC)
  mean time:        138.549 ns (22.18% GC)
  maximum time:     4.617 μs (95.96% GC)
  --------------
  samples:          10000
  evals/sample:     940
==================================================

==================================================
send frames to gpu
==================================================
BenchmarkTools.Trial:
  memory estimate:  553.31 KiB
  allocs estimate:  39
  --------------
  minimum time:     119.322 μs (0.00% GC)
  median time:      123.279 μs (0.00% GC)
  mean time:        140.553 μs (10.86% GC)
  maximum time:     2.255 ms (87.05% GC)
  --------------
  samples:          10000
  evals/sample:     1
==================================================

==================================================
consecutive view a batch
==================================================
BenchmarkTools.Trial:
  memory estimate:  4.59 KiB
  allocs estimate:  80
  --------------
  minimum time:     1.377 μs (0.00% GC)
  median time:      1.474 μs (0.00% GC)
  mean time:        2.001 μs (19.49% GC)
  maximum time:     575.531 μs (99.30% GC)
  --------------
  samples:          10000
  evals/sample:     10
==================================================

==================================================
get an element
==================================================
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     3.617 ns (0.00% GC)
  median time:      3.635 ns (0.00% GC)
  mean time:        3.657 ns (0.00% GC)
  maximum time:     19.967 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000
==================================================

==================================================
view a frame and do matrix mul
==================================================
BenchmarkTools.Trial:
  memory estimate:  31.98 KiB
  allocs estimate:  6
  --------------
  minimum time:     316.350 μs (0.00% GC)
  median time:      338.291 μs (0.00% GC)
  mean time:        347.118 μs (0.54% GC)
  maximum time:     11.323 ms (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1
==================================================

==================================================
view a frame and do matrix mul (use CartesianIndex)
==================================================
BenchmarkTools.Trial:
  memory estimate:  33.30 KiB
  allocs estimate:  51
  --------------
  minimum time:     378.907 μs (0.00% GC)
  median time:      403.400 μs (0.00% GC)
  mean time:        413.936 μs (0.51% GC)
  maximum time:     11.349 ms (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1
==================================================

==================================================
view frames and do matrix mul
==================================================
BenchmarkTools.Trial:
  memory estimate:  1010.80 KiB
  allocs estimate:  13
  --------------
  minimum time:     3.211 ms (0.00% GC)
  median time:      3.261 ms (0.00% GC)
  mean time:        3.345 ms (1.07% GC)
  maximum time:     14.094 ms (0.00% GC)
  --------------
  samples:          1489
  evals/sample:     1
==================================================

==================================================
get a frame and do matrix mul
==================================================
BenchmarkTools.Trial:
  memory estimate:  87.95 KiB
  allocs estimate:  39
  --------------
  minimum time:     379.192 μs (0.00% GC)
  median time:      401.962 μs (0.00% GC)
  mean time:        414.151 μs (1.11% GC)
  maximum time:     11.369 ms (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1
==================================================

==================================================
do matrix mul
==================================================
BenchmarkTools.Trial:
  memory estimate:  316.27 KiB
  allocs estimate:  17
  --------------
  minimum time:     1.406 ms (0.00% GC)
  median time:      1.452 ms (0.00% GC)
  mean time:        1.492 ms (0.85% GC)
  maximum time:     12.403 ms (0.00% GC)
  --------------
  samples:          3323
  evals/sample:     1
==================================================

==================================================
do matrix mul(for comparison)
==================================================
BenchmarkTools.Trial:
  memory estimate:  315.95 KiB
  allocs estimate:  8
  --------------
  minimum time:     1.198 ms (0.00% GC)
  median time:      1.262 ms (0.00% GC)
  mean time:        1.291 ms (0.97% GC)
  maximum time:     12.275 ms (0.00% GC)
  --------------
  samples:          3838
  evals/sample:     1
==================================================
=#


text:
"""The minimum mass considered."""
const mlow = 3.0

"""The maximum mass considered."""
const mhigh = 20.0

"""The reference redshift, at which rates are evaluated."""
const zref = 0.4

function _fm(m, a1, a2, mb)
    x = m/mb
    if m < mlow || m > mhigh
        -Inf
    elseif m < mb
        a1*log(x)
    else
        a2*log(x)
    end
end

"""
    make_log_dN(a1, a2, mb, mu_q, sigma_q)

Produces a function that computes `(m1, m2) -> log_dNdm1dm2` for the population
model.

The population model is

```math
m_1 m_2 \\frac{\\mathrm{d}N}{\\mathrm{d} m_1 \\mathrm{d} m_2 \\mathrm{d} V \\mathrm{d} t} = R f_m\\left( m_1 \\right) f_m\\left( m_2 \\right) \\exp\\left( -\\frac{\\left(q - \\mu_q\\right)^2}{2 \\sigma_q^2} + \\frac{\\left(1-\\mu_q\\right)^2}{2 \\sigma_q^2} \\right)^2
```
where
```math
f_m\\left( m \\right) = \\begin{cases}
m^{\\alpha_1} & m_\\mathrm{min} < m < m_b \\
m^{\\alpha_2} & m_b \\leq m < m_\\mathrm{max}
\\end{cases}
```
is the "common" black hole mass function, and a Gaussian-shaped pairing function
in ``q`` is used with parameters ``mu_q`` and ``sigma_q``.  The function
returned is normalized so that ``R`` measures the volumetric merger rate per
natural log mass squared at `m1 = m2 = m_b`.
"""

square(x) = x*x
function make_log_dN(a1, a2, mb, mu_q, sigma_q)
    function log_dN(m1, m2)
        q = m2/m1
        _fm(m1, a1, a2, mb) + _fm(m2, a1, a2, mb) - log(m1) - log(m2) - 0.5*square((q-mu_q)/sigma_q) + 0.5*square((1-mu_q)/sigma_q)
    end
    log_dN
end

"""
    make_dNdm1(a1, a2, mb, mu_q, sigma_q)

Returns a function that computes the marginal density (unnormalized) in `m1`.
"""
function make_dNdm1(a1, a2, mb, mu_q, sigma_q)
    log_dN = make_log_dN(a1, a2, mb, mu_q, sigma_q)
    m1 -> begin
        if m1 <= mlow
            zero(m1)
        elseif m1 > mhigh
            zero(m1)
        else
            n = Int(round(100*(log(m1)-log(mlow)))) + 2
            ms = exp.(range(log(mlow), stop=log(m1), length=n))
            dN = exp.(map(m2->log_dN(m1,m2), ms))
            trapz(ms, dN)
        end
    end
end

"""
    make_dNdm2(a1, a2, mb, beta)

Returns a function that computes the marginal density (unnormalized) in `m2`.
"""
function make_dNdm2(a1, a2, mb, mu_q, sigma_q)
    log_dN = make_log_dN(a1, a2, mb, mu_q, sigma_q)
    m2 -> begin
        if m2 < mlow
            zero(m2)
        elseif m2 > mhigh
            zero(m2)
        else
            n = Int(round(100*(log(mhigh) - log(m2)))) + 2
            ms = exp.(range(log(m2), stop=log(mhigh), length=n))
            dN = exp.(map(m1->log_dN(m1, m2), ms))
            trapz(ms, dN)
        end
    end
end

"""
Returns a function that computes the marginal (unnormalized) density in `q`.
"""
function make_dNdq(a1, a2, mb, mu_q, sigma_q)
    log_dN = make_log_dN(a1, a2, mb, mu_q, sigma_q)
    q -> begin
        if q < mlow/mhigh
            zero(q)
        elseif q > 1
            zero(q)
        else
            ml = mlow/q
            mh = mhigh
            n = Int(round(100*(log(mh)-log(ml)))) + 2
            ms = exp.(range(log(ml), stop=log(mh), length=n))
            dN = exp.(map(m1->log_dN(m1, q*m1) + log(m1), ms))
            trapz(ms, dN)
        end
    end
end

"""
    broken_pl_model(m1s, m2s, log_wts, m1s_sel, m2s_sel, log_pdraw, Ndraw)

Returns a Turing model for our broken power law mass function.

# Arguments

- `m1s`: An iterable containing iterables of the posterior samples for the
  source-frame primary mass for our observations.
- `m2s`: Same, but for secondary mass.
- `log_wts`: Iterable with iterables of the log of the prior weight assigned to
  each sample.
- `m1s_sel`: Samples of primary mass from the detected injections used to
  estimate detector sensitivity.
- `m2s_sel`: Secondary mass.
- `log_pdraw`: The log of the (normalized) density from which the injections
  were drawn.
- `Ndraw`: The number of injections drawn.

# Parameters

- `a1`: The low-mass power law slope.
- `a2`: The high-mass power law slope.
- `mb`: The break mass.
- `beta`: The pairing power law exponent.
- `R`: (generated) The volumetric merger rate per natural log mass squared at
  `m1 = m2 = mb`.
- `Neff_sel`: The effective number of samples (see [Farr
  (2019)](https://iopscience.iop.org/article/10.3847/2515-5172/ab1d5f)) in the
  Monte-Carlo selection integral.
"""
@model function broken_pl_model(m1s, m2s, log_wts, m1s_sel, m2s_sel, log_pdraw, Ndraw)
    Nobs = length(m1s)
    a1 ~ Uniform(0, 10)
    a2 ~ Uniform(-10, 0)
    mb ~ Uniform(sqrt(30), sqrt(200)) # Peaks somewhere around 10.
    mu_q ~ Uniform(0.2, 1.0)
    sigma_q ~ Uniform(0.1, 2)

    log_dN = make_log_dN(a1, a2, mb, mu_q, sigma_q)

    log_wts = map(m1s, m2s, log_wts) do mm1, mm2, log_wwt
        map(mm1, mm2, log_wwt) do m1, m2, log_wt
            log_dN(m1, m2) - log_wt
        end
    end

    Neff_samps = map(log_wts) do lws
        exp(2*logsumexp(lws) - logsumexp(2 .* lws))
    end

    log_likes = map(log_wts) do lw
        logsumexp(lw) - log(length(lw))
    end

    Turing.@addlogprob! sum(log_likes)

    log_wts_sel = map(m1s_sel, m2s_sel, log_pdraw) do m1, m2, logp
        log_dN(m1, m2) - logp
    end
    log_mu = logsumexp(log_wts_sel) - log(Ndraw)
    log_s2 = logsubexp(logsumexp(2 .* log_wts_sel) - 2.0*log(Ndraw), 2*log_mu - log(Ndraw))
    Neff_sel = exp(2*log_mu - log_s2)

    Turing.@addlogprob! -Nobs*log_mu

    mu_R = Nobs/exp(log_mu)
    sigma_R = sqrt(Nobs) / exp(log_mu)

    R = rand(Normal(mu_R, sigma_R))

    return (R = R, Neff_sel = Neff_sel, Neff_samps = Neff_samps)
end


text:
export ConvDiffFEMParam, getConvDiffFEMParam

"""
type ConvDiffFEMParam <: ForwardProbType

description of stationary convection diffusion forward problem

(- sig*Laplacian + v dot GRAD) u = f
with boundary conditions du/dn = 0 on Omega1, u=gd on Omega2

Here: Use a FEM discretization.

Construct an instance using getConvDiffFEMParam(M,kwargs...)

Fields:
	M         - describes computational mesh
	A         - stiffness matrix
	MassConst - mass matrix for piecewise constant rhs
	P         - receiver matrix
	sig       - scalar, >0
	v         - velocity (assumed to be spatially homogeneous)
	ubc       - field from Dirichlet boundary conditions
	rbc       - right hand side contribution from Dirichlet bc
	iddir     - indices of nodes at which we have Dirichlet BCs
	idint     - unconstrained nodes
	Fields    - stores fields, i.e., u
	Ainv      - factorization of PDE
"""
mutable struct ConvDiffFEMParam <: ForwardProbType
	M :: RegularMesh
	A :: SparseMatrixCSC{Float64} # stiffness matrix
	MassConst::SparseMatrixCSC{Float64} # mass matrix for piecewise constant source
	P :: AbstractArray{Float64} # receiver matrix
	sig ::Float64        # viscosity
	v::Array{Float64}    # velocity
	ubc::Array{Float64}  # field from Dirichlet boundary condition
	rbc::Array{Float64}  # right hand side from Dirichlet boundary condition
	iddir::Array         # indices of nodes on Dirichlet boundary
	idint::Array         # indices of nodes not on Dirichlet boudnary
	Fields::Array{Float64,1}
	Ainv::AbstractSolver
end

"""
function getConvDiffFEMParam

constructor for ConvDiffFEMParam

Required Input:
	M         - describes computational mesh

Keyword Arguments:
	sig       - scalar, >0
	v         - velocity (assumed to be spatially homogeneous)
	gd        - function for Dirichlet boundary conditions
	bc        - Tuple, defining the boundary conditions
	P         - receiver matrix
	Fields    - stores fields, i.e., u
	Ainv      - factorization of PDE
"""
function getConvDiffFEMParam(M::RegularMesh;   v=Matrix{Int}(I,M.dim,1), gd::Function=X->zeros(size(X,1)),
	                                           sig::Number=0.01,bc=(:dir,:neu,:neu,:neu,:neu,:neu), P=Diagonal(ones(prod(M.n.+1))),
											   Fields=zeros(0),Ainv=getJuliaSolver())

	if M.dim==2
		Mass, Mass_const, A = getFEMMatrices2D(M,sig=sig,v=v)
	else
		Mass, Mass_const, A = getFEMMatrices3D(M,sig=sig,v=v)
	end

	iddir = getBoundaryIndicesFEM(M, bc)  # Dirichlet nodes
	idint = setdiff(1:size(A,1),iddir) # interior nodes
	Abc = A[idint,iddir]
	A = A[idint,idint]
	Mass_const = Mass_const[idint,:]

	xd         = getNodalGrid(M)[iddir,:]
	ud         = gd(xd)
	ubc        = zeros(prod(M.n.+1));
	ubc[iddir] = ud
	rbc        = Abc*ud;

	return ConvDiffFEMParam(M,A,Mass_const,P,sig,v,ubc,rbc,iddir,idint,Fields,Ainv)
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
function Plots.plot(d::CircuitSimulation.TimeDomainData)
    scale = (1,"");
    mint = d.ts[2]
    if mint < 1e-15
        scale = (1e15, "f")
    elseif mint < 1e-12
        scale = (1e12, "p")
    elseif mint < 1e-9
        scale = (1e9, "n")
    elseif mint < 1e-6
        scale = (1e6, "m")
    elseif mint < 1e-3
        scale = (1e3, "m")
    end
    b = Plots.plot(
                   xaxis=(@sprintf("Time [%ss]",scale[2])),
                   size=(1200,800),
                   yaxis=("Time-domain Function"),
                   leg=false,
                  )
    Plots.plot!(b,d.ts*scale[1], d.vals,)
end

function Plots.plot(d::FrequencyDomainData; style=:magphase)
    scale = (1,"");
    maxf = maximum(d.fs)
    if maxf > 1e15
        scale = (1e-15, "P")
    elseif maxf > 1e12
        scale = (1e-12, "T")
    elseif maxf > 1e9
        scale = (1e-9, "G")
    elseif maxf > 1e6
        scale = (1e-6, "M")
    elseif maxf > 1e3
        scale = (1e-3, "k")
    end
    b = Plots.plot(
                   xaxis=(@sprintf("Frequency [%sHz]",scale[2])),
                   size=(1200,800),
                   layout=Plots.@layout([a; b]),
                   leg=[false false],
                  )
    if style==:magphase
        Plots.plot!(b,fftshift(d.fs)*scale[1],
                    hcat(fftshift(abs.(d.vals)),fftshift(rad2deg.(angle.(d.vals)))),
                    title=["Magnitude" "Phase"],
                    ylabel=["Magnitude (L1 norm)" "Phase [deg]"],
                   )
    elseif style==:realimag
        Plots.plot!(b,fftshift(d.fs)*scale[1],
                    hcat(fftshift(real.(d.vals)),fftshift(imag.(d.vals))),
                    title=["Real" "Imag"],
                    ylabel=["Real" "Imag"],
                   )
    else
        error("'style' must be either :magphase or :realimag.")
    end
end


text:
@testset "VanillaSVDD" begin

    dummy_data = [1.0 2.0 3.0; 4.0 5.0 6.0]
    vanilla_svdd = SVDD.VanillaSVDD(dummy_data)

    @testset "create" begin
        @test vanilla_svdd.state == SVDD.model_created
        @test_throws SVDD.ModelStateException SVDD.predict(vanilla_svdd, dummy_data)
    end

    @testset "initialize" begin
        init_strategy = SVDD.FixedParameterInitialization(MLKernels.GaussianKernel(0.5), 0.5)

        SVDD.initialize!(vanilla_svdd, init_strategy)
        @test vanilla_svdd.state == SVDD.model_initialized
        @test vanilla_svdd.data == dummy_data
        @test vanilla_svdd.C ≈ 0.5
        @test size(vanilla_svdd.K) == (3,3)
        @test all(diag(vanilla_svdd.K) .≈ 1.0)
    end

    @testset "fit" begin
        SVDD.fit!(vanilla_svdd, TEST_SOLVER)
        @test vanilla_svdd.state == SVDD.model_fitted
        @test length(vanilla_svdd.alpha_values) == size(dummy_data,2)
        @test vanilla_svdd.R > 0.0
        @test vanilla_svdd.alpha_values == SVDD.get_alpha_prime(vanilla_svdd)
    end

    @testset "predict" begin
        predictions = SVDD.predict(vanilla_svdd, dummy_data)
        @test length(predictions) == size(dummy_data, 2)
        @test predictions[1] == SVDD.predict(vanilla_svdd, hcat([1.0, 4.0]))[1]
    end

    @testset "params" begin
         @test all(map(x->SVDD.is_valid_param_value(vanilla_svdd, Val{:C}, x), [0.0, 0.1, 1.0]))
         @test !any(map(x->SVDD.is_valid_param_value(vanilla_svdd, Val{:C}, x), [-0.1, 1.1]))
     end
end


text:
# type for mortality, lapse, etc
abstract type Decrement end


text:
using MDBM
using LinearAlgebra, Arpack
function fooeig(x)
       mymx=Matrix{Float64}(I,5,5);
       mymx[2,1]=x;
       mymx[1,2]=x+5;
       abs(eigs(mymx, nev = 2)[1][1])-3
end
eig_problem=MDBM_Problem(fooeig,[Axis(-10:10)])

println(" X solution  //   Function val")
for k=1:12
solve!(eig_problem,1)
Xsol=getinterpolatedsolution(eig_problem)
print(Xsol[1])
print("  //   ")
println(maximum(abs.(fooeig.(Xsol[1]))))
end


text:
using Test

using RealNeuralNetworks.NeuralNets
using RealNeuralNetworks.Neurons

using CSV

@testset "test NeuralNets module..." begin
    neuronId2neuron = Dict{Int, Neuron}()
    neuronId2neuron[77625] = Neurons.load(joinpath(__DIR__, "../asset/77625.swc.bin"))
    neuronId2neuron[77641] = Neurons.load(joinpath(__DIR__, "../asset/77641.swc.bin"))

    df = CSV.load(joinpath(__DIR__, "../asset/77265.pre.synapses.csv")
    #Neurons.attach_pre_synapses()
    net = NeuralNet(neuronId2neuron)
end


text:
################################################################################
## Chi squared
################################################################################

#=---------------
Evaluate
---------------=#
function evaluate(dist::ChiSquared, a::T, b::T) where T <: AbstractFloat
    return (a-b)^2/(2*b)
end

function evaluate(dist::ChiSquared, a::AbstractVector{T}) where T <: AbstractFloat
    r = zero(T)
    @simd for i in eachindex(a)
        @inbounds r += (a[i] - 1.0)^2/2.0
    end
    return r
end

function evaluate(dist::ChiSquared, a::AbstractVector{T}, b::AbstractVector{T}) where T <: AbstractFloat
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    r = zero(T)
    @simd for i in eachindex(a, b)
        @inbounds r += (a[i]-b[i])^2/(2*b[i])
    end
    return r
end

#=---------------
gradient
---------------=#
function gradient(dist::ChiSquared, a::T, b::T) where T <: AbstractFloat
    return (a/b) - 1.0
end

function gradient(dist::ChiSquared, a::T) where T <: AbstractFloat
    return a - 1.0
end

function gradient!(u::Vector{T}, dist::ChiSquared, a::AbstractVector{T}, b::AbstractVector{T}) where T <: AbstractFloat
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    @simd for i = eachindex(a, b)
        ai = a[i]
        bi = b[i]
        @inbounds u[i] = (ai/bi)-1.0
    end
    u
end

function gradient!(u::Vector{T}, dist::ChiSquared, a::AbstractVector{T}) where T <: AbstractFloat
    @simd for i = eachindex(a)
        @inbounds u[i] = a[i]-1.0
    end
    u
end

#=---------------
hessian
---------------=#
function hessian(dist::ChiSquared, a::T, b::T) where T <: AbstractFloat
    return 1/b
end

function hessian(dist::ChiSquared, a::T) where T <: AbstractFloat
    return 1.0
end

function hessian!(u::Vector{T}, dist::ChiSquared, a::AbstractVector{T}, b::AbstractVector{T}) where T <: AbstractFloat
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    @simd for i = eachindex(a)
         @inbounds u[i] = 1/b[i]
    end
    u
end

function hessian!(u::Vector{T}, dist::ChiSquared, a::AbstractVector{T}) where T <: AbstractFloat
    @simd for i = eachindex(a)
        @inbounds u[i] = 1.0
    end
    u
end


text:
"Discrete filtering toolbox"
module DiscreteFiltering

using ApproxFun: Fun, chebyshevt, integrate, Interval, (..)
using ForwardDiff: derivative
using IntervalSets: ±
using LinearAlgebra: Diagonal, I, mul!, factorize, ldiv!, lu
using MLJLinearModels: LassoRegression, RidgeRegression, fit
using OrdinaryDiffEq: OrdinaryDiffEq
using OrdinaryDiffEq: ODEFunction, ODEProblem, QNDF, RK4, DiffEqArrayOperator, LinearExponential
using Parameters: @unpack
using SparseArrays: dropzeros!, sparse, spdiagm, spzeros
using NonNegLeastSquares: nonneg_lsq

# Domain
include("domain/domain.jl")
include("domain/get_npoint.jl")
include("domain/discretize.jl")

# Filter
include("filter/filter.jl")
include("filter/apply_filter.jl")
include("filter/apply_filter_int.jl")
include("filter/apply_filter_extend.jl")

# Matrix assembly
include("matrices/advection_matrix.jl")
include("matrices/diffusion_matrix.jl")
include("matrices/interpolation_matrix.jl")
include("matrices/filter_matrix.jl")
include("matrices/filter_matrix_meshwidth.jl")
include("matrices/reconstruction_matrix.jl")
include("matrices/reconstruction_matrix_meshwidth.jl")
include("matrices/get_W_R.jl")
include("matrices/fit_Cbar.jl")

# Equations
include("equations/equations.jl")
include("equations/solve.jl")
include("equations/solve_adbc.jl")

# Utils
include("utils/ridge.jl")
include("utils/sum_of_sines.jl")

# Domain
export ClosedIntervalDomain, PeriodicIntervalDomain, discretize, get_npoint

# Filter
export IdentityFilter
export TopHatFilter
export ConvolutionalFilter
export GaussianFilter
export apply_filter
export apply_filter_int
export apply_filter_extend

# Matrix assembly
export advection_matrix
export diffusion_matrix
export filter_matrix
export filter_matrix_meshwidth
export reconstruction_matrix
export reconstruction_matrix_meshwidth
export get_W_R
export fit_Cbar

# Equations
export AdvectionEquation, DiffusionEquation, BurgersEquation, solve, solve_adbc

# Utils
export sum_of_sines

end


text:
import Base: maximum, minimum, Fix2
export CircularListArray

# Implement interface functions directly on
include("circularlist_abstractarray.jl")

struct CircularListArray{T} <: AbstractArray{T,1}
    list::List{T}
    dict::Dict{T,Node{T}}
end
function CircularListArray(x; capacity=length(x))
    list = circularlist(x; capacity=capacity)
    dict = Dict(zip(x, list.nodes))
    CircularListArray(list, dict)
end
size(cla::CircularListArray) = size(cla.list)
getindex(cla::CircularListArray, x) = getindex(cla.list, x)
setindex!(cla::CircularListArray, v, i::Int) = setindex!(cla.list, v, i)
setindex!(cla::CircularListArray, X, I::UnitRange) = setindex!(cla.list, X, I)
IndexStyle(::CircularListArray) = IndexLinear()
function deleteat!(cla::CircularListArray, x)
    pop!.((cla.dict,), cla.list[x])
    deleteat!(cla.list, x)
end
function insert!(cla::CircularListArray, i::Integer, x)
    # TODO move dict logic here
    insert!(cla.list, i, x)
end
function insert!(cla::CircularListArray{T}, node::Node{T}, x) where T
    insert!(cla.list, node, x)
    for i = 1:length(x)
        cla.dict[node.next.data] = node.next
        node = node.next
    end
end
copy(cla::CircularListArray) = CircularListArray(copy(cla.list))
function findfirst(predicate::Function, cla::CircularListArray)
    findfirst(predicate, cla.list)
end
function findfirstnode( predicate::Base.Fix2{F,X}, cla::CircularListArray) where {F, X}
    if predicate.f ∈ ( isequal, (==) ,(===) )
        node = cla.dict[ predicate.x ]
    else
        findfirst(predicate, cla.list)
    end
end
maximum(cla::CircularListArray) = maximum(cla.list)
minimum(cla::CircularListArray) = minimum(cla.list)


text:
@testset "onda.annotation deprecations" begin
    root = mktempdir()
    possible_recordings = (uuid4(), uuid4(), uuid4())
    annotations = Annotation[Annotation(recording=rand(possible_recordings),
                                        id=uuid4(),
                                        span=TimeSpan(Second(rand(0:30)), Second(rand(31:60))),
                                        a=randstring('a':'z', 10),
                                        b=rand(Int, 1),
                                        c=rand(3)) for i in 1:50]
    annotations_file_path_1 = joinpath(root, "test-1.onda.annotations.arrow")
    annotations_file_path_2 = joinpath(root, "test-2.onda.annotations.arrow")
    cols = Tables.columns(annotations)
    io = IOBuffer()
    write_annotations(annotations_file_path_1, cols)
    Arrow.write(annotations_file_path_2, cols)
    write_annotations(io, cols)
    seekstart(io)
    for roundtripped in (read_annotations(annotations_file_path_1; validate_schema=false),
                         read_annotations(annotations_file_path_1; validate_schema=true),
                         read_annotations(annotations_file_path_2; validate_schema=false),
                         read_annotations(annotations_file_path_2; validate_schema=true),
                         Onda.materialize(read_annotations(io)),
                         read_annotations(seekstart(io); validate_schema=true))
        roundtripped = collect(Tables.rows(roundtripped))
        @test length(roundtripped) == length(annotations)
        for (r, a) in zip(roundtripped, annotations)
            @test NamedTuple(a) == NamedTuple(r)
            @test NamedTuple(a) == NamedTuple(Annotation(r))
        end
    end
    x = first(annotations)
    y = @compat Annotation(x.recording, x.id, x.span; x.a, x.b, x.c)
    @test x == y
    @test (@test_deprecated setproperties(y; a='+')) == Annotation(Tables.rowmerge(y; a='+'))
    df = DataFrame(annotations)
    @test Onda.gather(:recording, df) == Legolas.gather(:recording, df)

    names = (:recording, :id, :span)
    types = Tuple{Union{UInt128,UUID},Union{UInt128,UUID},Union{Onda.NamedTupleTimeSpan,TimeSpan}}
    @test (@test_deprecated Onda.validate_annotation_schema(nothing)) === nothing
    @test (@test_deprecated Onda.validate_annotation_schema(Tables.Schema(names, types))) === nothing
    @test_throws ArgumentError Onda.validate_annotation_schema(Tables.Schema((:x, :y), (Any, Any)))
end

@testset "onda.signal deprecations" begin
    root = mktempdir()
    possible_recordings = (uuid4(), uuid4(), uuid4())
    possible_sample_types = (UInt8, UInt16, UInt32, UInt64, Int8, Int16, Int32, Int64, Float32, Float64)
    signals = Signal[Signal(recording=rand(possible_recordings),
                            file_path=joinpath(root, "x_$(i)_file"),
                            file_format=rand(("lpcm", "lpcm.zst")),
                            span=TimeSpan(Second(rand(0:30)), Second(rand(31:60))),
                            kind="x_$i",
                            channels=["a_$i", "b_$i", "c_$i"],
                            sample_unit="unit_$i",
                            sample_resolution_in_unit=rand((0.25, 1)),
                            sample_offset_in_unit=rand((-0.25, 0.25)),
                            sample_type=rand(possible_sample_types),
                            sample_rate=rand((128, 50.5)),
                            a=join(rand('a':'z', 10)),
                            b=rand(Int, 1),
                            c=rand(3)) for i in 1:50]
    signals_file_path_1 = joinpath(root, "test-1.onda.signals.arrow")
    signals_file_path_2 = joinpath(root, "test-2.onda.signals.arrow")
    io = IOBuffer()
    write_signals(signals_file_path_1, signals)
    Arrow.write(signals_file_path_2, signals)
    write_signals(io, signals)
    seekstart(io)
    io2 = IOBuffer()
    write_signals(io2, signals; file=false)
    seekstart(io2)
    for roundtripped in (read_signals(signals_file_path_1; validate_schema=false),
                         read_signals(signals_file_path_1; validate_schema=true),
                         read_signals(signals_file_path_2; validate_schema=false),
                         read_signals(signals_file_path_2; validate_schema=true),
                         (tmp=read_signals(io); @test_deprecated Onda.materialize(tmp)),
                         (tmp=read_signals(io2); @test_deprecated Onda.materialize(tmp)),
                         read_signals(seekstart(io); validate_schema=true))
        roundtripped = collect(Tables.rows(roundtripped))
        @test length(roundtripped) == length(signals)
        for (r, s) in zip(roundtripped, signals)
            @test NamedTuple(s) == NamedTuple(r)
            @test NamedTuple(s) == NamedTuple(Signal(r))
        end
    end
    x = first(signals)
    y = @compat Signal(x.recording, x.file_path, x.file_format, x.span, x.kind, x.channels, x.sample_unit,
                       x.sample_resolution_in_unit, x.sample_offset_in_unit, x.sample_type, x.sample_rate;
                       x.a, x.b, x.c)
    @test x == y
    x = Onda.extract_samples_info(x)
    z = @test_deprecated (@compat Signal(x; y.recording, y.file_path, y.file_format, y.span, y.a, y.b, y.c))
    @test y == z
    @test (@test_deprecated setproperties(z; sample_rate=1.0)) == Signal(Tables.rowmerge(z; sample_rate=1.0))
    y = SamplesInfo(x.kind, x.channels, x.sample_unit, x.sample_resolution_in_unit, x.sample_offset_in_unit, x.sample_type, x.sample_rate)
    @test x == y
    @test (@test_deprecated setproperties(y; sample_rate=1.0)) == SamplesInfo(Tables.rowmerge(y; sample_rate=1.0))
    @test isnothing(@test_deprecated Onda.validate(y))
    df = DataFrame(signals)
    @test (@test_deprecated Onda.gather(:recording, df)) == Legolas.gather(:recording, df)

    names = (:recording, :file_path, :file_format, :span, :kind, :channels, :sample_unit, :sample_resolution_in_unit, :sample_offset_in_unit, :sample_type, :sample_rate)
    types = Tuple{Union{UInt128,UUID},Any,AbstractString,Union{Onda.NamedTupleTimeSpan,TimeSpan},AbstractString,AbstractVector{<:AbstractString},AbstractString,
                  Onda.LPCM_SAMPLE_TYPE_UNION,Onda.LPCM_SAMPLE_TYPE_UNION,AbstractString,Onda.LPCM_SAMPLE_TYPE_UNION}
    @test (@test_deprecated Onda.validate_signal_schema(nothing)) === nothing
    @test (@test_deprecated Onda.validate_signal_schema(Tables.Schema(names, types))) === nothing
    @test_throws ArgumentError Onda.validate_signal_schema(Tables.Schema((:x, :y), (Any, Any)))
end


text:
###### customizable show ########

const DECIMALS = 4

##### showing missing values ######

const NAN       = "NaN"
const NA        = "NA"
const BLACKHOLE = "\u2B24"
const DOTCIRCLE = "\u25CC"
const QUESTION  = "\u003F"
# const UNICORN   = "\u1F984" # doesn't render well

const MISSING = NAN


text:
module CircularArrays

export CircularArray, getindex, setindex!, advance_head!, reset_head!

# load files
include("circulararray.jl")

end # module


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

JLLWrappers.@generate_main_file_header("URCU")
JLLWrappers.@generate_main_file("URCU", UUID("aa747835-a391-587f-982f-064ff03f7d29"))
end  # module URCU_jll


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBasis

#    "Poi1" => (0, 1),
"1 node discrete point element",
type Poi1 <: AbstractBasis
end

function get_basis(::Element{Poi1}, ::Any, ::Any)
    return [1]
end

function get_dbasis(::Element{Poi1}, ::Any, ::Any)
    return [0]
end

function (::Element{Poi1})(::Any, ::Float64, ::Type{Val{:detJ}})
    return 1.0
end

function get_integration_order(::Poi1)
    return 1
end

function get_integration_points(::Poi1, ::Int64)
    return [ (1.0, (0.0, )) ]
end

function size(::Type{Poi1})
    return (0, 1)
end

function length(::Type{Poi1})
    return 1
end

function FEMBasis.get_reference_element_coordinates(::Type{Poi1})
    Vector{Float64}[[0.0]]
end

function inside(::Union{Type{Seg2}, Type{Seg3}, Type{Quad4}, Type{Quad8},
                        Type{Quad9}, Type{Pyr5}, Type{Hex8}, Type{Hex20},
                        Type{Hex27}}, xi)
    return all(-1.0 .<= xi .<= 1.0)
end

function inside(::Union{Type{Tri3}, Type{Tri6}, Type{Tri7}, Type{Tet4}, Type{Tet10}}, xi)
    return all(xi .>= 0.0) && (sum(xi) <= 1.0)
end

function get_reference_coordinates{B}(::Element{B})
    return get_reference_element_coordinates(B)
end


text:
function Sim_Model(Cell,Input,Def,Tk,SList,SOC,A0,B0,C0,D0,t)
    """
    Function to simulate generated reduced-order models.

    Sim_Model(Cell,Input,Def,Tk,SOC,A,B,C,D)

    """
    #Determine time span and allocate arrays
    tlength = size(Input,1)

    A = Array{Float64}(undef,size(A0[1]))
    B = Array{Float64}(undef,size(B0[1]))
    C = Array{Float64}(undef,size(C0[1]))
    D = Array{Float64}(undef,size(D0[1]))

    #Selecting SS Models
    for i in 1:length(SList)
        if SOC == SList[i] # Fix this, round?
            A = A0[i]
            B = B0[i]
            C = C0[i]
            D = D0[i]
        end
    end

    #Capturing Indices
    tfstr = Array{String}(undef,0,1)
    for i in 1:length(Cell.Transfer.tfs)
        t1 = t2 = Array{String}(undef,0,1)
        for j in 1:length(Cell.Transfer.Locs[i])
            t1 ="$(Cell.Transfer.tfs[i])_$(Cell.Transfer.Elec[i])"
            t2 = [t2; t1]
        end
        tfstr = [tfstr; t2]
    end

    CeInd = findall(isequal("C_e_Na"), tfstr)
    ϕ_ẽInd = findall(isequal("Phi_e_Na"), tfstr)
    CsePosInd = findall(isequal("C_se_Pos"), tfstr)
    CseNegInd = findall(isequal("C_se_Neg"), tfstr)
    ϕ_sNegInd = findall(isequal("Phi_s_Neg"), tfstr)
    ϕ_sPosInd = findall(isequal("Phi_s_Pos"), tfstr)
    ϕ_sePosInd = findall(isequal("Phi_se_Pos"), tfstr)
    ϕ_seNegInd = findall(isequal("Phi_se_Neg"), tfstr)
    FluxNegInd = findall(isequal("Flux_Neg"), tfstr)
    FluxPosInd = findall(isequal("Flux_Pos"), tfstr)

    #@infiltrate cond = true

    CeNeg = Cell.Transfer.Locs[1].<=Cell.Neg.L
    CeNegOffset = Cell.Transfer.Locs[1].<Cell.Neg.L
    CeSep = Cell.Transfer.Locs[1].<=Cell.Neg.L+Cell.Sep.L
    CeSepOffset = Cell.Transfer.Locs[1].<Cell.Neg.L+Cell.Sep.L
    CePos = Cell.Transfer.Locs[1].<=Cell.Const.Ltot
    CeNegInd = findall(CeNeg .== 1)
    CeSepInd = findall(CeSep.-CeNegOffset .== 1)
    CePosInd = findall(CePos.-CeSepOffset .== 1)

    csegain_neg = C[CseNegInd[1][1],1] #First Column in C Array (zeros column)
    csegain_pos = C[CsePosInd[1][1],1] #First Column in C Array (zeros column)

    #Memory Allocation
    θ_neg = Array{Float64}(undef,tlength,1) .= 0.
    θ_pos = Array{Float64}(undef,tlength,1) .= 0.
    jeq_neg = Array{Float64}(undef,tlength,1) .= 0.
    jeq_pos = Array{Float64}(undef,tlength,1) .= 0.
    x = Array{Float64}(undef,tlength+1,size(A,1)) .= 0.
    y = Array{Float64}(undef,tlength,size(C,1)) .= 0.
    Cse_Neg = Array{Float64}(undef,tlength,size(CseNegInd,1)) .= 0.
    Cse_Pos = Array{Float64}(undef,tlength,size(CsePosInd,1)) .= 0.
    Ce = Array{Float64}(undef,tlength,size(CeInd,1)) .= Cell.Const.ce0
    η0 = Array{Float64}(undef,tlength,1) .= 0.
    η_neg = Array{Float64}(undef,tlength,size(FluxNegInd,1)) .= 0.
    ηL = Array{Float64}(undef,tlength,1) .= 0.
    η_pos = Array{Float64}(undef,tlength,size(FluxPosInd,1)) .= 0.
    ϕ_ẽ1 = Array{Float64}(undef,tlength,size(ϕ_ẽInd,1)) .= 0.
    ϕ_ẽ2 = Array{Float64}(undef,tlength,size(CeInd,1)) .= 0.
    ϕ_se_neg_0 = Array{Float64}(undef,tlength,1) .= 0. #Replace with length of ϕ_seNegInd @ zero
    jNeg = Array{Float64}(undef,tlength,size(FluxNegInd,1)) .= 0.
    j0 = Array{Float64}(undef,tlength,1) .= 0.
    jPos = Array{Float64}(undef,tlength,size(FluxPosInd,1)) .= 0.
    jL = Array{Float64}(undef,tlength,1) .= 0.
    j0_CC_neg = Array{Float64}(undef,tlength,1) .= 0.
    Rtot_neg = Array{Float64}(undef,tlength,1) .= 0.
    Rtot_pos = Array{Float64}(undef,tlength,1) .= 0.
    j0_CC_pos = Array{Float64}(undef,tlength,1) .= 0.
    Uocp_Neg = Array{Float64}(undef,tlength,1) .= 0.
    Uocp_Pos = Array{Float64}(undef,tlength,1) .= 0.
    Cell_V = Array{Float64}(undef,tlength,1) .= 0.
    ϕ_e = Array{Float64}(undef,tlength,size(CeInd,1)) .= 0.
    Cell_SOC = Array{Float64}(undef,tlength,1) .= 0
    Iapp = Array{Float64}(undef,tlength+1,1) .= 0

    #Defining SOC
    SOC_Neg = SOC * (Cell.Neg.θ_100-Cell.Neg.θ_0) + Cell.Neg.θ_0
    SOC_Pos = SOC * (Cell.Pos.θ_100-Cell.Pos.θ_0) + Cell.Pos.θ_0
    θ_neg[1] = SOC_Neg
    θ_pos[1] = SOC_Pos
    Cell_SOC[1] = (SOC_Neg-Cell.Neg.θ_0)/(Cell.Neg.θ_100-Cell.Neg.θ_0)

    #Loop through time
    #Compute dependent variables (voltage, flux, etc.)
    for i in 0:(tlength-1)
        cs_neg_avg = x[i+1,1] * csegain_neg + SOC_Neg * Cell.Neg.cs_max < 0. ? 0. : x[i+1,1] * csegain_neg + SOC_Neg * Cell.Neg.cs_max #Zero if < 0
        cs_pos_avg = x[i+1,1] * csegain_pos + SOC_Pos * Cell.Pos.cs_max < 0. ? 0. : x[i+1,1] * csegain_pos + SOC_Pos * Cell.Pos.cs_max #Zero if < 0

        if cs_neg_avg > Cell.Neg.cs_max
            cs_neg_avg = Cell.Neg.cs_max
        end

        if cs_pos_avg > Cell.Pos.cs_max
            cs_pos_avg = Cell.Pos.cs_max
        end

        #Reaction Rates
        if Cell.Const.CellTyp == "Doyle_94"
            k_neg = Cell.Neg.k_norm#/Cell.Neg.cs_max/Cell.Const.ce0^(1-Cell.Neg.α)
            k_pos = Cell.Pos.k_norm#/Cell.Pos.cs_max/Cell.Const.ce0^(1-Cell.Pos.α)
            jeq_neg[i+1] = k_neg*sqrt(cs_neg_avg*Cell.Const.ce0*(Cell.Neg.cs_max-cs_neg_avg))
            jeq_pos[i+1] = k_pos*sqrt(cs_pos_avg*Cell.Const.ce0*(Cell.Pos.cs_max-cs_pos_avg))
        else
            #jeq_neg = Cell.Neg.k_norm*Cell.Neg.cs_max*(Cell.Const.ce0*(cs_neg_avg/Cell.Neg.cs_max*(1-cs_neg_avg/Cell.Neg.cs_max)))^(1-Cell.Neg.α)
            #jeq_pos = Cell.Pos.k_norm*Cell.Pos.cs_max*(Cell.Const.ce0*(cs_pos_avg/Cell.Pos.cs_max*(1-cs_pos_avg/Cell.Pos.cs_max)))^(1-Cell.Pos.α)
            k_neg = Cell.Neg.k_norm
            k_pos = Cell.Pos.k_norm
            jeq_neg[i+1] = Cell.Neg.k_norm*(Cell.Const.ce0*cs_neg_avg*(Cell.Neg.cs_max-cs_neg_avg))^(1-Cell.Neg.α)
            jeq_pos[i+1] = Cell.Pos.k_norm*(Cell.Const.ce0*cs_pos_avg*(Cell.Pos.cs_max-cs_pos_avg))^(1-Cell.Pos.α)
        end

        θ_neg[i+1] = cs_neg_avg/Cell.Neg.cs_max
        θ_pos[i+1] = cs_pos_avg/Cell.Pos.cs_max
        Cell_SOC[i+1] = (θ_neg[i+1]-Cell.Neg.θ_0)/(Cell.Neg.θ_100-Cell.Neg.θ_0)

        javg_neg = Iapp[i+1]/(Cell.Neg.as*F*Cell.Neg.L*Cell.Const.CC_A)
        javg_pos = Iapp[i+1]/(Cell.Pos.as*F*Cell.Pos.L*Cell.Const.CC_A)

        Arr_Factor = ((1/Cell.Const.T_ref)-(1/Tk[i+1]))/R
        κneg = Cell.Const.κf(mean(Ce[i+1,CeNegInd]))*exp(Cell.Const.Ea_κ*Arr_Factor)
        κpos = Cell.Const.κf(mean(Ce[i+1,CePosInd]))*exp(Cell.Const.Ea_κ*Arr_Factor)
        κsep = Cell.Const.κf(mean(Ce[i+1,CeSepInd]))*exp(Cell.Const.Ea_κ*Arr_Factor)
        σ_neg = Cell.Neg.σ*exp(Cell.Const.Ea_κ*Arr_Factor)
        σ_pos = Cell.Pos.σ*exp(Cell.Const.Ea_κ*Arr_Factor)
        κ_eff_Neg = κneg*(Cell.Neg.ϵ_e^(Cell.Neg.κ_brug))
        κ_eff_Sep = κsep*(Cell.Sep.ϵ_e^(Cell.Sep.κ_brug))
        κ_eff_Pos = κpos*(Cell.Pos.ϵ_e^(Cell.Pos.κ_brug))
        σ_eff_Neg = σ_neg*Cell.Neg.ϵ_s^Cell.Neg.σ_brug #Effective Conductivity Neg
        σ_eff_Pos = σ_pos*Cell.Pos.ϵ_s^Cell.Pos.σ_brug #Effective Conductivity Pos

        #Resistances
        #Rtot_neg[i] = (Tk[i]*R)/(F^2*mean([abs(jeq_neg[i]) abs(javg_neg)]))+Cell.Neg.RFilm
        Rtot_neg[i+1] = (Tk[i+1]*R)/(F^2*sqrt(jeq_neg[i+1]^2+javg_neg^2/4))+Cell.Neg.RFilm
        #Rtot_pos[i] = (Tk[i]*R)/(F^2*mean([abs(jeq_pos[i]) abs(javg_pos)]))+Cell.Pos.RFilm
        Rtot_pos[i+1] = (Tk[i+1]*R)/(F^2*sqrt(jeq_pos[i+1]^2+javg_pos^2/4))+Cell.Pos.RFilm

        #Condensing Variable
        ν_neg = @. Cell.Neg.L*sqrt((Cell.Neg.as*(1/κ_eff_Neg+1/σ_eff_Neg))/Rtot_neg[i+1])
        ν_pos = @. Cell.Pos.L*sqrt((Cell.Pos.as*(1/κ_eff_Pos+1/σ_eff_Pos))/Rtot_pos[i+1])

        #Relinearise dependent on ν, σ, κ
        D = D_Linear(Cell, ν_neg, ν_pos, σ_eff_Neg, κ_eff_Neg, σ_eff_Pos, κ_eff_Pos, κ_eff_Sep)

        #Interpolate C & D Matrices
        C = interp(C0,SList,Cell_SOC[i+1])
        #D = interp(D0,SList,Cell_SOC[i+1])
        #SS Output
        y[i+1,:] = C*x[i+1,:] + D*Iapp[i+1]

        #Concentrations & Force electrode concentration maximum
        Cse_Neg[i+1,:] = (SOC_Neg.*Cell.Neg.cs_max .+ y[i+1,CseNegInd]) > ones(size(Cse_Neg,2))*Cell.Neg.cs_max ? ones(size(Cse_Neg,2))*Cell.Neg.cs_max : (SOC_Neg.*Cell.Neg.cs_max .+ y[i+1,CseNegInd])
        Cse_Pos[i+1,:] = (SOC_Pos.*Cell.Pos.cs_max .+ y[i+1,CsePosInd]) > ones(size(Cse_Pos,2))*Cell.Pos.cs_max ? ones(size(Cse_Pos,2))*Cell.Pos.cs_max : (SOC_Pos.*Cell.Pos.cs_max .+ y[i+1,CsePosInd])
        Ce[i+1,:] = @. Cell.Const.ce0 + y[i+1,CeInd]

        #Potentials
        Uocp_Neg[i+1] = Cell.Const.Uocp("Neg",Cse_Neg[i+1,1]/Cell.Neg.cs_max)
        Uocp_Pos[i+1] = Cell.Const.Uocp("Pos",Cse_Pos[i+1,1]/Cell.Pos.cs_max)
        ϕ_se_neg_0[i+1] = y[i+1, ϕ_seNegInd[1]] + Uocp_Neg[i+1] #Location 0
        ϕ_ẽ1[i+1,:] = y[i+1,ϕ_ẽInd]
        ϕ_ẽ2[i+1,:] = @. ((Tk[i+1]*2*R*(1-Cell.Const.tpf(Ce[i+1,:])))/F)*(log(Ce[i+1,:]/Ce[i+1,1]))
        ϕ_e[i+1,:] = @. [0; ϕ_ẽ1[i+1,:]]+ϕ_ẽ2[i+1,:]-ϕ_se_neg_0[i+1]

        #Flux
        jNeg[i+1,:] = y[i+1,FluxNegInd]
        j0[i+1] = y[i+1,FluxNegInd[1]]
        jPos[i+1,:] = y[i+1,FluxPosInd]
        jL[i+1] = y[i+1,FluxPosInd[1]]

        #j0_CC_neg[i+1] = findmax([eps(); (k_neg*(Cell.Neg.cs_max-Cse_Neg[i+1,1])^(1-Cell.Neg.α))*((Cse_Neg[i+1,1]^Cell.Neg.α)*(Ce[i+1,1]^(1-Cell.Neg.α)))])[1]
        j0_neg = findmax([ones(size(Cse_Neg,2))*eps() (k_neg.*(Cse_Neg[i+1,:].^Cell.Neg.α).*(Ce[i+1,1].^(1-Cell.Neg.α))).*(Cell.Neg.cs_max.-Cse_Neg[i+1,:]).^(1-Cell.Neg.α)], dims=2)[1]

        η0[i+1] = Tk[i+1]*2*R/F*asinh(j0[i+1]/(2*j0_neg[1]))
        η_neg[i+1,:] = @. (Tk[i+1]*2*R)/F*asinh((jNeg[i+1,:])/(2*j0_neg))

        #j0_CC_pos[i+1] = findmax([eps(); (k_pos*(Cell.Pos.cs_max-Cse_Pos[i+1,1])^(1-Cell.Pos.α))*((Cse_Pos[i+1,1]^Cell.Pos.α)*(Ce[i+1,end]^(1-Cell.Pos.α)))])[1]
        j0_pos = findmax([ones(size(Cse_Pos,2))*eps() (k_pos.*(Cse_Pos[i+1,:].^Cell.Pos.α).*(Ce[i+1,1].^(1-Cell.Pos.α))).*(Cell.Pos.cs_max.-Cse_Pos[i+1,:]).^(1-Cell.Pos.α)], dims=2)[1]
        ηL[i+1] = (Tk[i+1]*2*R)/F*asinh(jL[i+1]/(2*j0_pos[1]))
        η_pos[i+1,:] = @. (Tk[i+1]*2*R)/F*asinh(jPos[i+1,:]/(2*j0_pos))

        #Cell Voltage
        Cell_V[i+1] = @. (Uocp_Pos[i+1]-Uocp_Neg[i+1]) + (ηL[i+1]-η0[i+1]) + (ϕ_ẽ1[i+1,end]+ϕ_ẽ2[i+1,end]) + (Cell.Pos.RFilm*jL[i+1]-Cell.Neg.RFilm*j0[i+1])*F

        #ϕ_s
        ϕ_s_neg = y[i+1,ϕ_sNegInd]
        ϕ_s_pos = @. y[i+1,ϕ_sPosInd] + Cell_V[i+1]

        #Interpolate A Matrix
        A = interp(A0,SList,Cell_SOC[i+1])
        B = interp(B0,SList,Cell_SOC[i+1])

        #Update States
        x[i+2,:] = A*x[i+1,:] + B*Iapp[i+1]

        if Def =="Power"
            Iapp[i+2] = Input[i+1,1]/Cell_V[i+1]
        else
            Iapp[i+2] = Input[i+1,1]
        end

    end

return Cell_V, Ce, jNeg, jPos, Rtot_neg, Rtot_pos, η0, ηL, η_neg, η_pos, ϕ_ẽ1, ϕ_ẽ2, Uocp_Neg, Uocp_Pos, ϕ_e, Cse_Neg, Cse_Pos, Cell_SOC, jeq_neg, jeq_pos, j0, jL, t
end

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
module FEDW
using Printf: @printf,@sprintf
using SparseArrays: spzeros, sparse
using QuantumLattices: AbstractFockOperator,Operators,FOperator,MAJORANA,CREATION,rank
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Mathematics.VectorSpaces: VectorSpace
using QuantumLattices.Mathematics.Combinatorics: Combinations

import QuantumLattices: kind,id,dimension
import QuantumLattices.Mathematics.AlgebraOverFields: idtype
import ExactDiagonalization.FED:FockBasis

export NCBasis,bkron,matrix,boplus

"""
      NCBasis{I}(vnstate::Vector{Int},vnpart::Vector{Vector{T}}) where {I<:Unsigned,T<:Tuple{Vararg{Int}}}

Generic fermionic/hard-core-bosonic basis, which uses conserved number of particles.
 e.g.
    vnstate=[m1,m2]; vnpart=[[(n1,n2),(n3,n4)],[(l1,l2),(l3,l4)]] denotes \$[(C_{m_1}^{n_1}⊕C_{m_1}^{n_2})⊗(C_{m_2}^{n_3}⊕C_{m_2}^{n_4})]⊕[(C_{m_1}^{l_1}⊕C_{m_1}^{l_2})⊗(C_{m_2}^{l_3}⊕C_{m_2}^{l_4})]\$ Hilbert space
============================================================================================================
vnstate=[4,6] denotes four and six single-particle states at the 1st and 2nd space,respectively.
vnpart=[[(2,),(3,)],[(1,),(1,)]];vnpart[1] denotes two and three particles at the 1st and 2nd space,respectively.
Thus, there are binomial(4,2)*binomial(6,3)+ binomial(4,1)*binomial(6,1) states.
"1110001010" denotes that "1010" ("111000") belongs to 1st (2nd) space.

#Examples:
```
julia> vnstate=[4,6];
julia> vnpart=[[(2,),(3,)]];
julia> NCBasis{UInt}(vnstate,vnpart).table.|>bitstring
120-element Vector{String}:
 "0000000000000000000000000000000000000000000000000000000001110011"
 "0000000000000000000000000000000000000000000000000000000010110011"
 "0000000000000000000000000000000000000000000000000000000100110011"
 ⋮
 "0000000000000000000000000000000000000000000000000000001101001100"
 "0000000000000000000000000000000000000000000000000000001110001100"
```
"""
struct NCBasis{I<:Unsigned} <: FockBasis{I}
    vnstate::Vector{Int}
    vnpart::Vector{Vector{Tuple{Vararg{Int}}}}
    table::Vector{I}
    sortflag::Bool
end
Base.kron(v1::Vector{NTuple{N,Int}},v2::Vector{NTuple{M,Int}}) where {N,M} = collect((v1[i]...,v2[j]...) for i=1:length(v1) for j=1:length(v2))
Base.kron(v1::Vector{T1},v2::Vector{T}) where {T<:Tuple{Vararg{Int}},T1<:Tuple{Vararg{Int}}} = collect((v1[i]...,v2[j]...) for i=1:length(v1) for j=1:length(v2))
Base.issorted(bs::NCBasis) = bs.sortflag

function NCBasis{I}(vnstate::Vector{Int},vnpart::Vector{Vector{T}}) where {I<:Unsigned,T<:Tuple{Vararg{Int}}}
    #@assert size(vnstate)==size(vnpart) "NCBasis error: size of vnstate ($(size(vnstate))) == size of vnpart ($(size(vnpart))) "
    table0=Vector{I}(undef,0)
    for k = 1:size(vnpart,1)
        @assert length(vnstate)==length(vnpart[k]) "NCBasis error: length of vnstate ($(length(vnstate))) == length of vnpart[$(k)] ($(length(vnpart[k]))) "
        @assert ~in(Tuple{}(),vnpart[k]) "NCBasis error: Tuple{}() in vnpart[$(k)]"
        vst=[0,cumsum(vnstate)...]
        temp = Vector{Vector{Tuple{Vararg{Int}}}}(undef,length(vnstate))
        for i=length(vnstate):-1:1
            temp[i]=vcat(collect.( Combinations{vnpart[k][i][j]}( vst[i]+1:vst[i+1]) for j=1:length(vnpart[k][i] ) )... )
        end
        tabl=reduce(kron,temp)
        for j=1:length(tabl)
            push!(table0,bitgen(tabl[j],I))
        end
    end
    return NCBasis{I}(vnstate,vnpart,sort(table0),true)
end

NCBasis{I}() where I<:Unsigned = NCBasis{I}([0],[[(0,)]])

function Base.show(io::IO,bs::NCBasis)
    @printf io "%s:\n" repr(bs)
    for i=1:dimension(bs)
        @printf io "  %s\n" string(bs[i],base=2,pad=sum(bs.vnstate))
    end
end
Base.repr(bs::NCBasis)=@sprintf "%s(%s)" nameof(typeof(bs)) join(id(bs),",")

"""
    bitgen(st::NTuple{N,Int},::Type{I}) where {N,I<:Integer}
Get a integer by NTuple `st`.

#Examples:
```
julia> bitgen((0,3),UInt)|>bitstring
"0000000000000000000000000000000000000000000000000000000000000100"
julia> bitgen(Tuple{}(),UInt)|>bitstring
"0000000000000000000000000000000000000000000000000000000000000000"
```
"""
function bitgen(st::NTuple{N,Int},::Type{I}) where {N,I<:Integer}
     eye,temp=one(I),zero(I)
    for j=1:length(st)
        temp=temp|eye<<(st[j]-1)
    end
    return temp
end

function dimension(bs::NCBasis)
        res0=0
        for k =1:size(bs.vnpart,1)
            res=1
            for i=1:length(bs.vnstate)
                res*=sum( binomial(bs.vnstate[i],bs.vnpart[k][i][j]) for j=1:length(bs.vnpart[k][i]) )
            end
            res0+=res
        end
        return res0
end

Base.getindex(bs::NCBasis,i::Int)=bs.table[i]

Base.searchsortedfirst(bs::NCBasis{I},i::I) where I<:Unsigned=findfirst(isequal(i),bs.table)
kind(::Type{<:NCBasis})=:nc
idtype(::Type{<:NCBasis}) = Tuple{  Vector{Int}, Vector{Vector{Tuple{Vararg{Int}}}} }
id(bs::NCBasis)=(bs.vnstate,bs.vnpart)

"""
   bkron(b1::I,n1::Int,b2::I,n2::Int) where I<:Integer -> I

Construct a state |b1>⊗|b2> by state |b1> and state |b2>. n1(n2) is the number of single-particle state |b1>(|b2>).

#Examples
```
julia> string(bkron(UInt(3),2,UInt(6),3),base=2,pad=5)
"11011"
```
"""
function bkron(b1::I,n1::Int,b2::I,n2::Int) where I<:Unsigned
            @assert n1+n2 <= length(bitstring(b1)) "bkron error: n1($(n1))+n2($(n2))<=length(bitstring(b1))($(length(bitstring(b1))))"
            b2<<n1|b1
        end

"""
   bkron(bs1::NCBasis{I},bs2::NCBasis{I}) where I<:Unsigned -> NCBasis{I}

Construct a basis {|b1>}⊗{|b2>} by basis {|b1>} and basis {|b2>}. If the `table`s of bs1 and bs2 are sorted, the `table` of result is sorted.
"""
function bkron(bs1::NCBasis{I},bs2::NCBasis{I}) where I<:Unsigned
            @assert size(bs2.vnstate,1)==size(bs1.vnstate,1) "bkron error: length of bs2.vnstate [$(length(bs2.vnstate))]== length of bs1.vnstate [$(length(bs1.vnstate))]"
            n2=length(bs2.vnpart)
            n1=length(bs1.vnpart)
            nvp=Vector{Vector{Tuple{Vararg{Int}}}}(undef,n1*n2)
            nvst=vcat(bs1.vnstate,bs2.vnstate)
            for  i= 1:n2
                for j=1:n1
                    nvp[j+(i-1)*n1]=vcat(bs1.vnpart[i],bs2.vnpart[j])
                end
            end
            table=collect(bkron(bs1.table[i],sum(bs1.vnstate),bs2.table[j],sum(bs2.vnstate)) for i=1:length(bs1.table) for j=1:length(bs2.table))
            #sortflag = bs1.sortflag && bs2.sortflag
            sortflag = false
            NCBasis{I}(nvst,nvp,table,sortflag)
        end

"""
   boplux(bs1::NCBasis{I},bs2::NCBasis{I}) where I<:Unsigned -> NCBasis{I}

Construct a basis {|b1>}⊕{|b2>} by basis {|b1>} and basis {|b2>}. The `table` is not sorted.
"""
function boplus(bs1::NCBasis{I},bs2::NCBasis{I}) where I<:Unsigned
        @assert bs2.vnstate==bs1.vnstate "boplus error: bs2.vnstate [$(bs2.vnstate)]== bs1.vnstate [$(bs1.vnstate)]"
        nvst=bs1.vnstate
        nvp=vcat(bs1.vnpart,bs2.vnpart)
        table=vcat(bs1.table,bs2.table)
        NCBasis{I}(nvst,nvp,table,false)
end

"""
    matrix(operator::AbstractFockOperator,basisbra::FockBasis,basisket::FockBasis,table) -> SparseMatrixCSC{valtype(operator),Int}
    matrix(operator::AbstractFockOperator,basis::FockBasis,table) -> SparseMatrixCSC{valtype(operator),Int}
    matrix(opts::Operators,basisbra::FockBasis,basisket::FockBasis,table,dtype::Type{<:Number}=Complex{Float}) -> SparseMatrixCSC{dtype,Int}
    matrix(opts::Operators,basis::FockBasis,table,dtype::Type{<:Number}=Complex{Float}) -> SparseMatrixCSC{dtype,Int}

Get the CSC-formed sparse matrix representation of an operator / a set of operators.
"""
function matrix(operator::AbstractFockOperator,basisbra::FockBasis,basisket::FockBasis, table)
    #@assert rank(operator)%2==0 "matrix error: two bases are needed for odd ranked Fock operators."
    ndata,data,indices,indptr=0,zeros(valtype(operator),dimension(basisket)),zeros(Int,dimension(basisbra)),zeros(Int,dimension(basisket))
    eye,temp,flag=one(basisket|>eltype),zeros(basisket|>eltype,rank(operator)+1),true
    seqs = NTuple{rank(operator), Int}(table[id(operator)[i]]-1 for i = 1:rank(operator))
    nambus = NTuple{rank(operator), Bool}(oid.index.nambu==CREATION for oid in id(operator))

    for i=1:dimension(basisket)
        flag=true
        temp[1]=basisket[i]
        for j=1:rank(operator)
            (temp[j]&eye<<seqs[end+1-j]>0)==nambus[end+1-j] && (flag=false; break)
            temp[j+1]=nambus[end+1-j] ? temp[j]|eye<<seqs[end+1-j] : temp[j]&~(eye<<seqs[end+1-j])
        end
        if flag
            (ndata+=1)
            nsign=0
            isa(operator,FOperator) && for j=1:rank(operator) for k=0:seqs[end+1-j]-1 (temp[j]&eye<<k>0 && (nsign+=1)) end end
            findind = searchsortedfirst(basisbra,temp[end])
            if findind == nothing
                ndata -= 1
            else
                indices[ndata]=findind
                data[ndata]=(-1)^nsign*operator.value
                indptr[ndata]=i
            end
        end

    end
    if ndata==0
        return spzeros(valtype(operator),dimension(basisbra),dimension(basisket))
    else
        return sparse(indices[1:ndata],indptr[1:ndata],data[1:ndata],dimension(basisbra),dimension(basisket))
    end
end

matrix(operator::AbstractFockOperator,basis::FockBasis,table) = matrix(operator,basis,basis,table)

function matrix(opts::Operators,basisbra::FockBasis,basisket::FockBasis,table,dtype::Type{<:Number}=Complex{Float})
    result=spzeros(dtype,dimension(basisbra),dimension(basisket))
    for opt in values(opts) result=result+matrix(opt,basisbra,basisket,table) end
    return result
end
matrix(opts::Operators,basis::FockBasis,table,dtype::Type{<:Number}=Complex{Float}) = matrix( opts,basis,basis,table,dtype)

end #module


text:
import Printf.@sprintf

# 9 tuning games.
const TUNING_SUITE = [
    "BeamRider",
    "DemonAttack",
    "DoubleDunk",
    "IceHockey",
    "MsPacman",
    "Pooyan",
    "RoadRunner",
    "Robotank",
    "Zaxxon",
]

# 36 testing games.
const TESTING_SUITE = [
    "Alien",
    "Amidar",
    "Assault",
    "Asterix",
    "Atlantis",
    "BankHeist",
    "BattleZone",
    "Boxing",
    "Breakout",
    "Carnival",
    "Centipede",
    "ChopperCommand",
    "CrazyClimber",
    "Enduro",
    "FishingDerby",
    "Freeway",
    "Frostbite",
    "Gopher",
    "Gravitar",
    "Hero",
    "Jamesbond",
    "Kangaroo",
    "Krull",
    "KungFuMaster",
    "NameThisGame",
    "Phoenix",
    "Pong",
    "Qbert",
    "Riverraid",
    "Seaquest",
    "SpaceInvaders",
    "StarGunner",
    "TimePilot",
    "UpNDown",
    "VideoPinball",
    "WizardOfWor",
    "YarsRevenge",
]

# Total of 45 games.
const ALL = cat(TUNING_SUITE, TESTING_SUITE, dims=1)

function fetch_rl_unplugged_atari(src, dest)
    try run(`which gsutil`) catch x throw("gsutil not found, install gsutil to proceed further") end

    run(`gsutil -m cp $src $dest`)
    return dest
end

num_shards = 100

function rl_unplugged_atari_init()
    for game in ALL
        for run in 1:5
            for index in 0:99
                register(
                    DataDep(
                        "rl-unplugged-atari-$game-$run-$index",
                        """
                        Dataset: RL Unplugged atari
                        Credits: https://arxiv.org/abs/2006.13888
                        Url: https://github.com/deepmind/deepmind-research/tree/master/rl_unplugged
                        Authors: Caglar Gulcehre, Ziyu Wang, Alexander Novikov, Tom Le Paine,
                        Sergio Gómez Colmenarejo, Konrad Zolna, Rishabh Agarwal,
                        Josh Merel, Daniel Mankowitz, Cosmin Paduraru, Gabriel
                        Dulac-Arnold, Jerry Li, Mohammad Norouzi, Matt Hoffman,
                        Ofir Nachum, George Tucker, Nicolas Heess, Nando deFreitas
                        Year: 2020

                        Data accompanying [RL Unplugged: Benchmarks for Offline Reinforcement Learning].
                        The dataset is generated by running an online DQN agent and recording transitions
                        from its replay during training with sticky actions Machado et al., 2018. As stated
                        in Agarwal et al., 2020, for each game we use data from five runs with 50 million
                        transitions each. States in each transition include stacks of four frames to be able
                        to do frame-stacking with our baselines. We release datasets for 46 Atari games.
                        For details on how the dataset was generated, please refer to the paper.
                        Atari is a standard RL benchmark. We recommend you to try offline RL methods
                        on Atari if you are interested in comparing your approach to other state of the
                        art offline RL methods with discrete actions.
                        """,
                        "gs://rl_unplugged/atari/$game/"*@sprintf("run_%i-%05i-of-%05i", run, index, num_shards);
                        fetch_method = fetch_rl_unplugged_atari
                    )
                )
            end
        end
    end
end


text:
using Compat
using Compat.Test

using Yao
using Yao.Zoo
using Yao.Blocks

witness_vec = zeros(1<<4)
witness_vec[[0, 3, 5, 10, 12, 15].+1] = 1
score_func(out_probs::Vector) = witness_vec'*out_probs

n = 4
depth = 4
pairs = [1=>2, 3=>4, 2=>3, 4=>1]
generator = diff_circuit(n, depth, pairs)

loss_func = () -> -score_func(apply!(zero_state(nqubits(witness_vec)), generator) |> probs)

import Base: gradient
function gradient(rots::Vector{<:RotationGate})
    ptb = perturb(loss_func, rots, π/2)
    (ptb[:,1] - ptb[:,2])/2
end

function train(generator::AbstractBlock, g_learning_rate::Real, niter::Int)
    rots = collect_rotblocks(generator)
    for i in 1:niter
        ggrad = gradient(rots)
        dispatch!(-, generator, ggrad.*g_learning_rate)
        if i%5==1 println("Step $i, loss = $(loss_func())") end
    end
end

train(generator, 0.1, 100)

using PyPlot
pl = apply!(zero_state(n), generator)|>probs
bar(0:1<<n-1, pl)
ylabel("p", fontsize=14)
savefig("/home/leo/jcode/Yao.jl/docs/src/assets/figures/trained22bs.png", dpi=300)

# biggest 3 indices
using PyPlot
using Yao.Intrinsics
configs = (pl |> sortperm)[end-2:end] .- 1 |> bitarray(n)

function show_configs(configs)
    M = size(configs, 2)
    for c in 1:M
        subplot(100+M*10+c)
        imshow(reshape(configs[:,c], 2, 2), vmin=0, vmax=1)
        axis("off")
    end
end

#show_configs(configs)
figure(figsize=(6,2))
show_configs([0, 3, 5, 10, 12, 15] |> bitarray(n))
savefig("/home/leo/jcode/Yao.jl/docs/src/assets/figures/22bs.png", dpi=300)

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
## Optimization Model for dtu42380a1
## Mar 11, 2020

using JuMP, Gurobi, CSV

include("./get_data.jl")

model = Model(with_optimizer(Gurobi.Optimizer, Presolve=0, OutputFlag=0));

## PRODUCTION VARIABLES
@variable(model, x[1:L,1:P,1:T] >= 0)
## SOURCE VARIABLES (PREBLENDING AND PACKAGING)
@variable(model, y[1:L,1:S,1:P,1:T] >= 0);
## MACHINE FLOW VARIABLES (ALLOCATION OF REGULAR TIME AND OVERTIME FOR BLENDING AND PACKAGING MACHINES)
@variable(model, mBR[1:L, 1:P, 1:T] >= 0);
@variable(model, mBO[1:L, 1:P, 1:T] >= 0);
@variable(model, mPR[1:L, 1:P, 1:T] >= 0);
@variable(model, mPO[1:L, 1:P, 1:T] >= 0);
## INVENTORY VARIABLES
@variable(model, I[1:L,1:P,1:T] >= 0);
## DEMAND NOT MET VARIABLES
@variable(model, delta[1:L,1:P,1:T] >= 0);

## OBJECTIVE FUNCTION
@objective(model, Min, sum( #=
=#           sum(y[i,s,p,t]*c[s] for s = 1:S)    #=  material cost
=#         + mBR[i,p,t]*cBR                      #= regular time blending machine cost
=#         + mBO[i,p,t]*cBO                      #= overtime blending machine cost
=#         + mPR[i,p,t]*cPR                      #= regular time packaging machine cost
=#         + mPO[i,p,t]*cPO                      #= overtime time packaging machine cost
=#         + I[i,p,t]*h[p]                       #= inventory cost
=#         + PI*delta[i,p,t]                     #= unsatisfied demand penalty
=#         for i = 1:L, p = 1:P, t = 1:T)
);

## PROCESS CONSTRAINTS - material and machine inflow must match (according to
##   the requirements) product outflow
@constraint(model, [i = 1:L, p = 1:P, t = 1:T,s = 1:S],
  y[i,s,p,t] == pi[p,s]*x[i,p,t]);
@constraint(model, [i = 1:L, p = 1:P, t = 1:T],
  mBR[i,p,t] + mBO[i,p,t] ==  muB[p]*x[i,p,t]);
@constraint(model, [i = 1:L, p = 1:P, t = 1:T],
  mPR[i,p,t] + mPO[i,p,t] ==  muP[p]*x[i,p,t]);
## CAPACITY CONSTRAINTS
@constraint(model, [i = 1:L, t = 1:T],
  sum(mBR[i,p,t] for p = 1:P) <= 24*5*60*60);
@constraint(model, [i = 1:L, t = 1:T],
  sum(mPR[i,p,t] for p = 1:P) <= 24*5*60*60);
@constraint(model, [i = 1:L, t = 1:T],
  sum(mBO[i,p,t] for p = 1:P) <= 24*1*60*60);
@constraint(model, [i = 1:L, t = 1:T],
  sum(mPO[i,p,t] for p = 1:P) <= 24*1*60*60);
## STORAGE CONSTRAINTS - inflow (production and start-of-period inventory) must
##   equal outflow (demand and end-of-period inventory) at each DC
@constraint(model, [i = 1:L, p = 1:P, t = 1  ],
  x[i,p,t] + I0[i,p] == I[i,p,t] + D[i,p,t] - delta[i,p,t]);
@constraint(model, [i = 1:L, p = 1:P, t = 2:T],
  x[i,p,t] + I[i,p,t-1] == I[i,p,t] + D[i,p,t] - delta[i,p,t]);

@constraint(model, [i = 1:L, p = 1:P, t = 1:T],
  I[i,p,t] >= ss[i,p,t]);

optimize!(model)

include("./export_result.jl")


text:
function main(n)
  n = round(Int, n / 2 * 2)
  a = rand(0:100_000_000, n, n)
  b = a
  c = a * b
  v = round(Int, n/2) + 1
  println(c[v, v])
end

function when_isMainModule()
  n = 100
  if length(ARGS) >= 1
    n = parse(Int, ARGS[1])
  end
  main(n)
end

when_isMainModule()

#########
# Results on MacOS + i5-5257U (Broadwell dual-core mobile 2.7GHz, turbo 3.1)
# Input 1500x1500 random large int64 matrix

# Nim 0.17.3 (devel) NO openMP. Compilation option: "-d:release --passc:-march=native"
# Julia v6.0
# Python 3.5.2 + Numpy 1.12

# Nim: 1.66s, 90 MB
# Julia: 4.49s, 185.2 MB
# Python Numpy: 9.49s, 55.8 MB

#########
# Results on Linux + E3-1230v5 (Skylake quad-core 3.4 GHz, turbo 3.8)
# Input 1500x1500 random large int64 matrix

# Nim 0.17.3 (devel) WITH OpenMP. Compilation option: "-d:release --passc:-march=native -d:openmp"
# Julia v6.0
# Python 3.6.2 + Numpy 1.12 compiled from source with march=native and openBLAS

# Nim + OpenMP: 0.36s, 55.5 MB
# Julia: 3.11s, 207.6 MB
# Python Numpy: 8.03s, 58.9 MB

max_stars_repo_name: dylanagreen/Arraymancer

max_issues_repo_name: dylanagreen/Arraymancer

max_forks_repo_name: dylanagreen/Arraymancer


text:
export solve_vf_terminal!, solve_vf_infill!, zero!

zero!(x::AbstractArray{T}) where {T<:Real} = fill!(x, zero(T))

function solve_vf_terminal!(EV::AbstractArray3)
    @views zero!(EV[:,:,end])
    @views zero!(EV[:,:,exploratory_terminal(wp)])
end

function solve_vf_terminal!(EV::AbstractArray3, dEV::AbstractArray4, dEVσ::AbstractArray3, wp::AbstractUnitProblem)
    @views zero!(EV[:,:,end])
    @views zero!(dEV[:,:,:,end])
    @views zero!(dEVσ[:,:,end])

    exp_trm = exploratory_terminal(wp)
    @views zero!(EV[  :,:,  exp_trm])
    @views zero!(dEV[ :,:,:,exp_trm])
    @views zero!(dEVσ[:,:,  exp_trm])
end

solve_vf_terminal!(evs::dcdp_Emax, wp::AbstractUnitProblem) = solve_vf_terminal!(evs.EV, evs.dEV, evs.dEVσ, wp)
solve_vf_terminal!(evs::dcdp_Emax, prim::dcdp_primitives) = solve_vf_terminal!(evs.EV, evs.dEV, evs.dEVσ, prim.wp)

# ---------------------------------------------

function solve_vf_infill!(evs::dcdp_Emax, t::dcdp_tmpvars, p::dcdp_primitives, θt::AbstractVector, σ::Real, dograd::Bool, itype::Tuple; maxit0::Integer=40, maxit1::Integer=20, vftol::Real=1e-9)

    EV         = evs.EV
    dEV        = evs.dEV
    wp         = p.wp
    Πz         = p.Πz
    β          = p.β

    nz,nψ,nS = size(EV)
    dmaxp1 = _dmax(wp)+1

    # ------------------------ size checks ----------------------------------

    (nz,nψ,dmaxp1) == size(t.ubVfull)         || throw(DimensionMismatch())
    (nz,nz) == size(t.IminusTEVp) == size(Πz) || throw(DimensionMismatch())
    (nz,nψ) == size(t.tmp) == size(t.lse)     || throw(DimensionMismatch())

    if dograd
        nθ = size(dEV,3)
        (nz,nψ,nθ,nS)     == size(dEV)        || throw(DimensionMismatch())
        (nz,nψ,nθ,dmaxp1) == size(t.dubVfull) || throw(DimensionMismatch())
    end

    # ------------------------ compute things ----------------------------------

    for i in ind_inf(wp)
        idxd, idxs, horzn = dp1space(wp,i), collect(sprimes(wp,i)), _horizon(wp,i)

        tvw = dcdp_tmpvars_view(t,idxd)
        ubV  = tvw.ubVfull
        dubV = tvw.dubVfull
        @views EV0 = EV[:,:,i]
        @views dEV0 = dEV[:,:,:,i]

        fillflows!(ubV, flow, p, θt, σ, i, itype...)
        @views ubV .+= β .* EV[:,:,idxs]

        if dograd
            fillflows_grad!(dubV, flowdθ, p, θt, σ, i, itype...)
            fill!(dEV0, 0.0)
            @views dubV .+= β .* dEV[:,:,:,idxs]
        end

        if horzn == :Finite
            if dograd
                vfit!(EV0, dEV0, tvw, p)
            else
                vfit!(EV0, tvw, p)
            end

        elseif horzn == :Infinite
            converged, iter, bnds =  solve_inf_vfit_pfit!(EV0, tvw, p; vftol=vftol, maxit0=maxit0, maxit1=maxit1)
            converged || @warn "Did not converge at state $i after $iter pfit. McQueen-Porteus bnds: $bnds. θt = $θt, σ = $σ"

            if dograd
                # TODO: only allows 0-payoff if no action
                ubV[:,:,1] .= β .* EV0
                gradinf!(dEV0, tvw, p)   # note: destroys ubV & dubV
            end
        else
            throw(error("i = $i, horzn = $horzn but must be :Finite or :Infinite"))
        end
    end
end


text:
const total_statements = Ref(0)
const succesful = Ref(0)

macro compile(statement)
    global succesful,total_statements
    total_statements[] += 1
    printstyled("compiling statement $(total_statements[]) \r",color = :reverse)
    try
        succesful[] +=  Base.eval(__module__,statement)
    catch error
        printstyled("\nWarning: ",bold = true , color = :yellow);
        printstyled("parsing error\n");
        @show error
        printstyled("$(__source__.file):$(__source__.line)\n\n",color=:cyan)
    end
end

reset_count() = begin
    succesful[] = 0
    total_statements[] = 0;
end

compile_summary() = begin
    printstyled("\n\nDONE: ",bold = true , color = :yellow);
    print("$(succesful[]) out of $(total_statements[]) succesfully compiled\n");
    printstyled("\nCreating sysimg...\n",bold = true , color = :red);
end


text:
include("roman_numerals.jl")

struct AtomicStage
    element::String
    name::String
    stage::Int
    stage_rn::RomanNumeral
    n_levels::Int
    g::Vector{Int}
    χ::typeof([1.0, 2.0]u"J")
    χ_ion::Unitful.Energy
    function AtomicStage(element, stage, g, χ, χ_ion)
        n_levels = length(χ)
        stage_rn = RomanNumeral(stage)
        stage = convert(Int, stage_rn)  # Ensures it can be called with string or Int
        name = string(element, "_", repr(stage_rn))
        new(element, name, stage, stage_rn, n_levels, g, χ, χ_ion)
    end
end