text:
# random initialization

function randinit(X, k::Integer; normalize::Bool=false, zeroh::Bool=false)
    p, n = size(X)
    T = eltype(X)

    W = rand(T, p, k)
    if normalize
        normalize1_cols!(W)
    end

    H = zeroh ? zeros(T, k, n) : rand(T, k, n)
    return W, H
end

# NNDSVD: Non-Negative Double Singular Value Decomposition
#
# Reference
# ----------
#   C. Boutsidis, and E. Gallopoulos. SVD based initialization: A head
#   start for nonnegative matrix factorization. Pattern Recognition, 2007.
#
function _nndsvd!(X, W, Ht, inith::Bool, variant::Int)

    p, n = size(X)
    k = size(W, 2)
    T = eltype(W)

    # compute SVD
    (U, s, V) = svd(X, thin=true)

    # main loop
    v0 = variant == 0 ? zero(T) :
         variant == 1 ? convert(T, mean(X)) : convert(T, mean(X) * 0.01)

    for j = 1:k
        x = view(U,:,j)
        y = view(V,:,j)
        xpnrm, xnnrm = posnegnorm(x)
        ypnrm, ynnrm = posnegnorm(y)
        mp = xpnrm * ypnrm
        mn = xnnrm * ynnrm

        vj = v0
        if variant == 2
            vj *= rand(T)
        end

        if inith
            if mp >= mn
                scalepos!(view(W,:,j), x, 1 / xpnrm, vj)
                scalepos!(view(Ht,:,j), y, s[j] * mp / ypnrm, vj)
            else
                ss = sqrt(s[j] * mn)
                scaleneg!(view(W,:,j), x, 1 / xnnrm, vj)
                scaleneg!(view(Ht,:,j), y, s[j] * mn / ynnrm, vj)
            end
        else
            if mp >= mn
                scalepos!(view(W,:,j), x, 1 / xpnrm, vj)
            else
                scaleneg!(view(W,:,j), x, 1 / xnnrm, vj)
            end
        end
    end
end

function nndsvd(X, k::Integer; zeroh::Bool=false, variant::Symbol=:std)

    p, n = size(X)
    T = eltype(X)
    ivar = variant == :std ? 0 :
           variant == :a   ? 1 :
           variant == :ar  ? 2 :
           throw(ArgumentError("Invalid value for variant"))

    W = Array{T,2}(p, k)
    H = Array{T,2}(k, n)
    if zeroh
        Ht = reshape(view(H,:,:), (n, k))
        _nndsvd!(X, W, Ht, false, ivar)
        fill!(H, 0)
    else
        Ht = Array{T,2}(n, k)
        _nndsvd!(X, W, Ht, true, ivar)
        for j = 1:k
            for i = 1:n
                H[j,i] = Ht[i,j]
            end
        end
    end
    return (W, H)
end

function posnegnorm(x::AbstractArray{T}) where T
    pn = zero(T)
    nn = zero(T)
    for i = 1:length(x)
        @inbounds xi = x[i]
        if xi > zero(T)
            pn += abs2(xi)
        else
            nn += abs2(xi)
        end
    end
    return (sqrt(pn), sqrt(nn))
end

function scalepos!(y, x, c::T, v0::T) where T<:Number
    @inbounds for i = 1:length(y)
        xi = x[i]
        if xi > zero(T)
            y[i] = xi * c
        else
            y[i] = v0
        end
    end
end

function scaleneg!(y, x, c::T, v0::T) where T<:Number
    @inbounds for i = 1:length(y)
        xi = x[i]
        if xi < zero(T)
            y[i] = - (xi * c)
        else
            y[i] = v0
        end
    end
end


text:
# depth limited
function dls(depth, goal, start)
    frontier = [createtree(start)];

    while !isempty(frontier)
        @timeit to "pop" node = pop!(frontier);

        if node.state == goal
            return node
        end
        if node.pathcost > depth
            continue;

        elseif !iscycle(node, 3)
            @timeit to "Moves" acts = moves(node.state);
            @timeit to "Expand" res = expansion(node.state);
            for (action, parent) in zip(acts, res)
                @timeit to "Node Creation" newnode = addnode(action, node, parent);
                @timeit to "Push" push!(frontier, newnode);
            end
        end
    end

    return nothing
end

# iterative
function idfs(goal, start; initdepth=1, limit=100, step=1)
    reset_timer!(to::TimerOutput)
    searchrange = initdepth:step:limit;
    for d in searchrange
            @timeit to "DLS" solnode = dls(d, goal, start)
        if !isnothing(solnode)
            displaysolution(solnode)
            return solnode
        end
    end

    println("We have gone too deep, no solutions found!  (limit=$limit)")
    return nothing
end

ext: jl

max_stars_repo_name: brice-allen/npuz

max_issues_repo_name: brice-allen/npuz

max_forks_repo_name: brice-allen/npuz


text:
import FEM.Nvec
import FEM.dNmatrix
import FEM.Jmatrix
import FEM.dNdxmatrix
import FEM.FENode2
import FEM.Point2

facts("FEM.Interpolator") do

context("FEM.Interpolator.LinTrigInterp") do
      interp = LinTrigInterp()

      node_1 = FENode2(1, [1.3, 0.5])
      node_2 = FENode2(2, [1.6, 1.0])
      node_3 = FENode2(3, [0.2, 0.9])

      loc_coords = Point2(0.5, 0.3)

      N = Nvec(interp, loc_coords)

      @fact norm(N - [0.5, 0.3, 1.0 - 0.5 - 0.3]) => roughly(0.0)

end # context

end # facts


text:
export
    AbstractCustomNGraph,
    AbstractCustomNDiGraph,
    AbstractCustomNTree,
    AbstractCustomNETree,
    AbstractCustomTree,

    get_graph,
    get_vtx_ids,
    get_vtx_map,
    get_nodes,

    get_vtx,
    get_vtx_id,
    get_node,
    get_parent,

    get_edge,
    replace_edge!,
    # delete_edge!,
    # add_custom_edge!,

    node_val,
    edge_val,
    node_id, # may cause an issue
    edge_source,
    edge_target,

    set_vtx_map!,
    insert_to_vtx_map!,

    replace_node!,
    add_node!,
    make_node,
    make_edge,
    add_child!,
    add_parent!,
    rem_node!,
    rem_nodes!

abstract type AbstractCustomGraph <: AbstractGraph{Int} end

"""
    abstract type AbstractCustomNGraph{G,N,ID} <: AbstractCustomGraph

An abstract Custom Graph type, with an underlying graph of type `G`, nodes of
type `N` with ids of type `ID`. All concrete subtypes `CG<:AbstractCustomNGraph`
must implement the following methods: `get_graph(g::CG)`,`get_vtx_ids(g::CG)`,
`get_vtx_map(g::CG)`, and `get_nodes(g::CG)`. These methods are implemented by
default if `g::CG` has the following fields:
- `graph     ::CG`
- `nodes     ::Vector{N}`
- `vtx_map   ::Dict{ID,Int}`
- `vtx_ids   ::Vector{ID}`

Abstract subtypes of `AbstractCustomNGraph{G,N,ID}` include:
- `AbstractCustomNEGraph{G,N,E,ID}` - a graph with custom edges of type `E`
"""
abstract type AbstractCustomNGraph{G,N,ID} <: AbstractCustomGraph end

"""
    abstract type AbstractCustomNEGraph{G,N,E,ID} <: AbstractCustomNGraph{G,N,ID}

An abstract Custom Graph type, with an underlying graph of type `G`, nodes of
type `N` with ids of type `ID`, and edges of type `E`. All concrete subtypes
`CG<:AbstractCustomNEGraph` must implement the required `AbstractCustomNGraph`
interface in addition to the following methods:
- `out_edges(g::CG)` : returns an integer-indexed forward adjacency list `fadj`
    such that `fadj[u::Int][v::Int]` contains the custom edge associated
    with `u → v`.
- `in_edges(g::CG)` : returns an integer-indexed backward adjacency list `badj`
    such that `badj[v::Int][u::Int]` contains the custom edge associated
    with `u → v`.
The above methods are implemented by default if `g::CG` has the following
fields:
- `outedges  ::Vector{Dict{Int,E}}`
- `inedges   ::Vector{Dict{Int,E}}`
"""
abstract type AbstractCustomNEGraph{G,N,E,ID} <: AbstractCustomNGraph{G,N,ID} end

const AbstractCustomNDiGraph{N,ID} = AbstractCustomNGraph{DiGraph,N,ID}
const AbstractCustomNEDiGraph{N,E,ID} = AbstractCustomNEGraph{DiGraph,N,E,ID}

"""
    abstract type AbstractCustomNTree{N,ID} <: AbstractCustomNDiGraph{N,ID}

Abstract custom graph type with tree edge structure.
"""
abstract type AbstractCustomNTree{N,ID} <: AbstractCustomNDiGraph{N,ID} end
abstract type AbstractCustomNETree{N,E,ID} <: AbstractCustomNEDiGraph{N,E,ID} end
const AbstractCustomTree = Union{AbstractCustomNTree,AbstractCustomNETree}

# Common interface
_graph_type(::AbstractCustomNGraph{G,N,ID}) where {G,N,ID}  = G
_node_type(::AbstractCustomNGraph{G,N,ID}) where {G,N,ID}   = N
_id_type(::AbstractCustomNGraph{G,N,ID}) where {G,N,ID}     = ID
"""
    get_graph(g::AbstractCustomNGraph{G,N,ID})

return the underlying graph of type `G`.
"""
get_graph(g::AbstractCustomGraph)   = g.graph
get_graph(g::AbstractGraph)   = g
"""
    get_vtx_ids(g::AbstractCustomNGraph{G,N,ID})

return the vector `Vector{ID}` of `g`'s unique vertex ids.
"""
get_vtx_ids(g::AbstractCustomNGraph) = g.vtx_ids
"""
    get_vtx_map(g::AbstractCustomNGraph{G,N,ID})

return a data structure (e.g, a 'Dict{ID,Int}') mapping node id to node index.
"""
get_vtx_map(g::AbstractCustomNGraph) = g.vtx_map
"""
    get_nodes(g::AbstractCustomNGraph{G,N,ID})

Return the vector `Vector{N}` of `g`'s nodes.
"""
get_nodes(g::AbstractCustomNGraph)   = g.nodes

Base.zero(g::G) where {G<:AbstractCustomNGraph} = G(graph=_graph_type(g)())

get_vtx(g::AbstractCustomNGraph,v::Int) = v
get_vtx(g::AbstractCustomNGraph{G,N,ID},id::ID) where {G,N,ID} = get(get_vtx_map(g), id, -1)
get_vtx(g::AbstractCustomNGraph{G,N,ID},node::N) where {G,N,ID} = get_vtx(g,node.id)
get_vtx_id(g::AbstractCustomNGraph,v::Int)             = get_vtx_ids(g)[v]
get_node(g::AbstractCustomNGraph,v) = get_nodes(g)[get_vtx(g,v)]

for op in [:edgetype,:ne,:nv,:vertices,:edges,:is_cyclic,:topological_sort_by_dfs,:is_directed,:is_connected]
    @eval LightGraphs.$op(g::AbstractCustomNGraph) = $op(get_graph(g))
end
for op in [:outneighbors,:inneighbors,:indegree,:outdegree,:has_vertex]
    @eval LightGraphs.$op(g::AbstractCustomNGraph,v::Int) = $op(get_graph(g),v)
    @eval LightGraphs.$op(g::AbstractCustomNGraph,id) = $op(g,get_vtx(g,id))
end
for op in [:bfs_tree]
    @eval LightGraphs.$op(g::AbstractCustomNGraph,v::Int;kwargs...) = $op(get_graph(g),v;kwargs...)
    @eval LightGraphs.$op(g::AbstractCustomNGraph,id;kwargs...) = $op(get_graph(g),get_vtx(g,id);kwargs...)
end
for op in [:has_edge] #,:add_edge!,:rem_edge!]
    @eval LightGraphs.$op(s::AbstractCustomNGraph,u,v) = $op(get_graph(s),get_vtx(s,u),get_vtx(s,v))
end

"""
    node_id(node)

Return the id of a node. Part of the required node interface for nodes in an
 `AbstractCustomNGraph`.
"""
node_id(node)       = node.id

"""
    node_val(node)

Return the value associated with a node. Part of the optional node interface for
nodes in an `AbstractCustomNGraph`.
"""
node_val(node)      = node.val

"""
    edge_source(edge)

Returns the ID of the source node of an edge. Part of the required interface for
edges in an `AbstractCustomNEGraph`.
"""
edge_source(edge)   = edge.src
"""
    edge_target(edge)

Returns the ID of the target node of an edge. Part of the required interface for
edges in an `AbstractCustomNEGraph`.
"""
edge_target(edge)   = edge.dst
"""
    edge_val(edge)

Returns the value associated with an edge. Part of the optional interface for
edges in an `AbstractCustomNEGraph`.
"""
edge_val(edge)      = edge.val

function set_vtx_map!(g::AbstractCustomNGraph,node,id,v::Int)
    @assert nv(g) >= v
    get_vtx_map(g)[id] = v
    get_nodes(g)[v] = node
end
function insert_to_vtx_map!(g::AbstractCustomNGraph,node,id,idx::Int=nv(g))
    push!(get_vtx_ids(g), id)
    push!(get_nodes(g), node)
    set_vtx_map!(g,node,id,idx)
end

# TODO code a convenient macro to implement this
#
# for op in node_accessor_interface
#     @eval $op(g::AbstractCustomNGraph,v) = $op(get_node(g,v))
#     @eval $op(g::AbstractCustomNGraph) = map(v->$op(get_node(g,v)), vertices(g))
# end
# for op in node_mutator_interface
#     @eval $op(g::AbstractCustomNGraph,v,val) = $op(get_node(g,v),val)
#     @eval $op(g::AbstractCustomNGraph,val) = begin
#         for v in vertices(g)
#             $op(get_node(g,v),val)
#         end
#     end
# end

"""
    replace_node!(g::AbstractCustomNGraph{G,N,ID},node::N,id::ID) where {G,N,ID}

Replace the current node associated with `id` with the new node `node`.
"""
function replace_node!(g::AbstractCustomNGraph{G,N,ID},node::N,id::ID) where {G,N,ID}
    v = get_vtx(g, id)
    @assert v != -1 "node id $(string(id)) is not in graph and therefore cannot be replaced"
    set_vtx_map!(g,node,id,v)
    node
end

"""
    add_node!(g::AbstractCustomNGraph{G,N,ID},node::N,id::ID) where {G,N,ID}

Add `node` to `g` with associated `id`.
"""
function add_node!(g::AbstractCustomNGraph{G,N,ID},node::N,id::ID) where {G,N,ID}
    @assert !has_vertex(g, id) "Trying to add node $(string(node)) with id $(string(id)) to g, but a node with id $(string(id)) already exists in g"
    add_vertex!(get_graph(g))
    insert_to_vtx_map!(g,node,id,nv(g))
    add_edge_lists!(g)
    node
end
"""
    make_node(g::AbstractCustomNGraph{G,N,ID},val,id)

Construct a node of type `N` from val and id. This method must be implemented
for whatever custom node type is used.
"""
function make_node end
add_node!(g::AbstractCustomNGraph{G,N,ID},val,id::ID) where {G,N,ID} = add_node!(g,make_node(g,val,id),id)
# add_node!(g::AbstractCustomNEGraph{G,N,E,ID},val,id::ID) where {G,N,E,ID} = add_node!(g,make_node(g,val,id),id)
add_node!(g::AbstractCustomNGraph,val,id) = add_node!(g,make_node(g,val,id))
replace_node!(g::AbstractCustomNGraph{G,N,ID},val,id::ID) where {G,N,ID} = replace_node!(g,make_node(g,val,id),id)

"""
    function add_child!(graph,parent,node,id)

add node `child` to `graph` with id `id`, then add edge `parent` → `child`
"""
function add_child!(g::AbstractCustomNGraph{G,N,ID},parent,child,id) where {G,N,ID}
    n = add_node!(g,child,id)
    if add_edge!(g,parent,id)
        return n
    else
        rem_node!(g,id)
        return nothing
    end
end

"""
    function add_parent!(graph,child,parent,id)

add node `parent` to `graph` with id `id`, then add edge `parent` → `child`
"""
function add_parent!(g::AbstractCustomNGraph{G,N,ID},child,parent::N,id::ID) where {G,N,ID}
    n = add_node!(g,parent,id)
    if add_edge!(g,id,child)
        return n
    else
        rem_node!(g,id)
        return nothing
    end
end
for op in [:add_child!,:add_parent!]
    @eval $op(g::AbstractCustomNGraph{G,N,ID},u,v::N) where {G,N,ID} = $op(g,u,v,v.id)
end

"""
    swap_with_end_and_delete!(vec,v)

Replaces `vec[v]` with `last(vec)`, then removes the last element from `vec`
"""
function swap_with_end_and_delete!(vec::Vector,v)
    vec[v] = last(vec)
    pop!(vec)
    return vec
end

"""
    rem_node!

removes a node (by id) from g.
Note about LightGraphs.rem_vertex!:
"internally the removal is performed swapping the vertices `v` and `nv(G)``, and
removing the last vertex `nv(G)` from the graph"
"""
function rem_node!(g::AbstractCustomNGraph{G,N,ID}, id::ID) where {G,N,ID}
    v = get_vtx(g, id)
    rem_vertex!(get_graph(g), v)
    swap_with_end_and_delete!(get_nodes(g),v)
    swap_with_end_and_delete!(get_vtx_ids(g),v)
    # deleteat!(get_nodes(g), v)
    delete!(get_vtx_map(g), id)
    # if v was not at end, need to update id_map too.
    if v <= nv(g)
        get_vtx_map(g)[get_vtx_ids(g)[v]] = v
    end
    # deleteat!(get_vtx_ids(g), v)
    # @assert length(get_vtx_ids(g)) == length(get_vtx_map(g)) == nv(get_graph(g))
    # for vtx in v:nv(get_graph(g))
    #     n_id = get_vtx_ids(g)[vtx]
    #     get_vtx_map(g)[n_id] = vtx
    # end
    delete_from_edge_lists!(g,v) # no effect except for AbstractCustomNEGraph
    g
end
rem_node!(g::AbstractCustomNGraph, v) = rem_node!(g,get_vtx_id(g,v))
function rem_nodes!(g::AbstractCustomNGraph, vtxs::Vector)
    node_ids = collect(map(v->get_vtx_id(g,v), vtxs))
    for id in node_ids
        rem_node!(g,id)
    end
    g
end

# Edge graph interface
_edge_type(::AbstractCustomNEGraph{G,N,E,ID}) where {G,N,E,ID} = E

"""
    in_edges(g::AbstractCustomNEGraph{G,N,E,ID})

Returns an integer-indexed backward adjacency list `badj` (e.g.,
`badj::Vector{Dict{Int,E}}`) such that `badj[v::Int][u::Int]` contains the
custom edge associated with `u → v`.
"""
in_edges(g::AbstractCustomNEGraph) = g.inedges
in_edges(g::AbstractCustomNEGraph, v) = in_edges(g)[get_vtx(g,v)]
in_edge(g::AbstractCustomNEGraph, v, u) = in_edges(g,v)[get_vtx(g,u)]
function set_edge!(g::AbstractCustomNEGraph{G,N,E,ID}, u,v, edge::E) where {G,N,E,ID}
    out_edges(g,u)[get_vtx(g,v)] = edge
    in_edges(g,v)[get_vtx(g,u)] = edge
    # TODO for undirected graph, insert reversed edge too?
    edge
end

"""
    out_edges(g::AbstractCustomNEGraph{G,N,E,ID})

Returns an integer-indexed forward adjacency list `fadj` (e.g.,
`fadj::Vector{Dict{Int,E}}`) such that `fadj[u::Int][v::Int]` contains the
custom edge associated with `u → v`.
"""
out_edges(g::AbstractCustomNEGraph) = g.outedges
out_edges(g::AbstractCustomNEGraph, v) = out_edges(g)[get_vtx(g,v)]
out_edge(g::AbstractCustomNEGraph, u, v) = out_edges(g,u)[get_vtx(g,v)]
function get_edge(g::AbstractCustomNEGraph,u,v)
    @assert has_edge(g,u,v) "Edge $u → $v does not exist."
    return out_edge(g,u,v)
end
get_edge(g::AbstractCustomNEGraph,edge) = get_edge(g,edge_source(edge),edge_target(edge))
LightGraphs.has_edge(g::AbstractCustomNEGraph,e) = has_edge(g,edge_source(e),edge_target(e))

add_edge_lists!(g::AbstractCustomNGraph) = g
function add_edge_lists!(g::AbstractCustomNEGraph{G,N,E,ID}) where {G,N,E,ID}
    push!(g.outedges, Dict{Int,E}())
    push!(g.inedges, Dict{Int,E}())
    return g
end
delete_from_edge_lists!(g::AbstractCustomNGraph,v::Int) = g
function delete_from_edge_lists!(g::AbstractCustomNEGraph,v::Int)
    swap_with_end_and_delete!(in_edges(g),v)
    swap_with_end_and_delete!(out_edges(g),v)
    # deleteat!(in_edges(g),v)
    # deleteat!(out_edges(g),v)
    return g
end

# TODO Make a GraphType with default constructible edges
function add_custom_edge!(g::AbstractCustomNEGraph{G,N,E,ID},u,v,edge::E) where {G,N,E,ID}
    if !is_legal_edge(g,u,v,edge)
        @warn "An edge $u → $v is illegal in $g"
        return false
    end
    if add_edge!(get_graph(g),get_vtx(g,u),get_vtx(g,v))
        set_edge!(g,u,v,edge)
        return true
    end
    @warn "Cannot add edge $u → $v. Does an edge already exist?"
    # println("Cannot add edge $u → $v. Does an edge already exist?")
    return false
end
add_custom_edge!(g::AbstractCustomNEGraph,edge) = add_custom_edge!(g,edge_source(edge),edge_target(edge),edge)
# LightGraphs.add_edge!(g::AbstractCustomNEGraph,args...) = add_custom_edge!(g,args...)

add_custom_edge!(g::AbstractCustomNGraph,edge) = add_custom_edge!(g,edge_source(edge),edge_target(edge))
add_custom_edge!(g::AbstractCustomNGraph,u,v) = add_edge!(get_graph(g),get_vtx(g,u),get_vtx(g,v))
add_custom_edge!(g::AbstractCustomNGraph,u,v,args...) = add_custom_edge!(g,u,v)
# LightGraphs.add_edge!(g::AbstractCustomNGraph,u,v,edge) = add_custom_edge!(g,u,v) # no custom edge type here
# LightGraphs.add_edge!(g::AbstractCustomNGraph,u,v) = add_custom_edge!(g,u,v) # no custom edge type here
LightGraphs.add_edge!(g::AbstractCustomGraph,args...) = add_custom_edge!(g,args...) # no custom edge type here

"""
    make_edge(g::G,u,v,val) where {G}

Construct an edge `u → v` of type `_edge_type(g)` based on val. Default behavior
is to throw an error.
"""
function make_edge(g::G,u,v,val) where {G}
    throw(ErrorException(
    """
    MethodError: `make_edge(g,u,v,val)` not implemented.
    To add an edge to a graph g::$G, either:
    - add the edge explicitly using `add_edge!(g,edge,[u,v])
    - make a default constructor for the edge type
    - implement `GraphUtils.make_edge(g,u,v,val)`
    """
    ))
end
make_edge(g,u,v) = make_edge(g,u,v,nothing)
function add_custom_edge!(g::AbstractCustomNEGraph,u,v,val)
    add_custom_edge!(g,u,v,make_edge(g,u,v,val))
end
function add_custom_edge!(g::AbstractCustomNEGraph,u,v)
    add_custom_edge!(g,u,v,make_edge(g,u,v))
end
# function LightGraphs.add_edge!(g::AbstractCustomNEGraph,u,v,val)
#     add_edge!(g,u,v,make_edge(g,u,v,val))
# end
# function LightGraphs.add_edge!(g::AbstractCustomNEGraph,u,v)
#     add_edge!(g,u,v,make_edge(g,u,v))
# end
function replace_edge!(g::AbstractCustomNEGraph{G,N,E,ID},u,v,edge::E) where {G,N,E,ID}
    if has_edge(g,u,v)
        set_edge!(g,u,v,edge)
        return true
    end
    @warn "graph does not have edge $u → $v"
    return false
end
replace_edge!(g::AbstractCustomNEGraph,u,v,val) = replace_edge!(g,u,v,make_edge(g,u,v,val))
function replace_edge!(g::AbstractCustomNEGraph{G,N,E,ID},old_edge,edge::E) where {G,N,E,ID}
    replace_edge!(g,edge_source(old_edge),edge_target(old_edge),edge)
end

function delete_edge!(g::AbstractCustomNEGraph, u, v)
    if rem_edge!(get_graph(g),get_vtx(g,u),get_vtx(g,v))
        delete!(out_edges(g,u),get_vtx(g,v))
        delete!(in_edges(g,v),get_vtx(g,u))
        return true
    end
    @warn "Cannot remove edge $u → $v. Does it exist?"
    return false
end
delete_edge!(g::AbstractCustomNGraph, u, v) = rem_edge!(get_graph(g),get_vtx(g,u),get_vtx(g,v))
delete_edge!(g,edge) = delete_edge!(g,edge_source(edge),edge_target(edge))
LightGraphs.rem_edge!(g::AbstractCustomGraph, args...) = delete_edge!(g,args...)
# LightGraphs.rem_edge!(g::AbstractCustomNEGraph, u, v) = delete_edge!(g,u,v)
# LightGraphs.rem_edge!(g::AbstractCustomNEGraph, edge) = rem_edge!(g,edge_source(edge),edge_target(edge))

# Tree interface
get_parent(g::AbstractCustomTree,v) = get(inneighbors(g,v),1,-1)
is_legal_edge(g,u,v) = true
is_legal_edge(g,u,v,e) = is_legal_edge(g,u,v)
is_legal_edge(g::AbstractCustomTree,u,v) = !(has_vertex(g,get_parent(g,v)) || get_vtx(g,u) == get_vtx(g,v))
# function LightGraphs.add_edge!(g::AbstractCustomNTree,u,v)
function add_custom_edge!(g::AbstractCustomTree,u,v)
    is_legal_edge(g,u,v)
    if !is_legal_edge(g,u,v)
        return false
    end
    add_edge!(get_graph(g),get_vtx(g,u),get_vtx(g,v))
end

################################################################################
################################ Concrete Types ################################
################################################################################

export
    CustomNode,
    CustomEdge,
    CustomNGraph,
    CustomNDiGraph,
    CustomNEGraph,
    CustomNEDiGraph,
    CustomNTree,
    CustomNETree,
    NGraph,
    NEGraph,
    NTree,
    NETree

"""
    CustomGraph

An example concrete subtype of `AbstractCustomNGraph`.
"""
@with_kw struct CustomNGraph{G<:AbstractGraph,N,ID} <: AbstractCustomNGraph{G,N,ID}
    graph               ::G                     = G()
    nodes               ::Vector{N}             = Vector{N}()
    vtx_map             ::Dict{ID,Int}          = Dict{ID,Int}()
    vtx_ids             ::Vector{ID}            = Vector{ID}() # maps vertex uid to actual graph node
end
const CustomNDiGraph{N,ID} = CustomNGraph{DiGraph,N,ID}

"""
    CustomNEGraph{G,N,E,ID}

Custom graph type with custom edge and node types.
"""
@with_kw struct CustomNEGraph{G,N,E,ID} <: AbstractCustomNEGraph{G,N,E,ID}
    graph               ::G                     = G()
    nodes               ::Vector{N}             = Vector{N}()
    vtx_map             ::Dict{ID,Int}          = Dict{ID,Int}()
    vtx_ids             ::Vector{ID}            = Vector{ID}()
    inedges             ::Vector{Dict{Int,E}}   = Vector{Dict{Int,E}}()
    outedges            ::Vector{Dict{Int,E}}   = Vector{Dict{Int,E}}()
end
const CustomNEDiGraph{N,E,ID} = CustomNEGraph{N,E,ID}

"""
    CustomNTree

An example concrete subtype of `AbstractCustomNTree`.
"""
@with_kw struct CustomNTree{N,ID} <: AbstractCustomNTree{N,ID}
    graph               ::DiGraph               = DiGraph()
    nodes               ::Vector{N}             = Vector{N}()
    vtx_map             ::Dict{ID,Int}           = Dict{ID,Int}()
    vtx_ids             ::Vector{ID}             = Vector{ID}()
end

"""
    CustomNETree{G,N,E,ID}

Custom tree type with custom edge and node types.
"""
@with_kw struct CustomNETree{N,E,ID} <: AbstractCustomNETree{N,E,ID}
    graph               ::DiGraph               = DiGraph()
    nodes               ::Vector{N}             = Vector{N}()
    vtx_map             ::Dict{ID,Int}          = Dict{ID,Int}()
    vtx_ids             ::Vector{ID}            = Vector{ID}()
    inedges             ::Vector{Dict{Int,E}}   = Vector{Dict{Int,E}}()
    outedges            ::Vector{Dict{Int,E}}   = Vector{Dict{Int,E}}()
end

"""
    CustomNode{N,ID}

A custom node type. Fields:
- `id::ID`
- `val::N`
"""
struct CustomNode{N,ID}
    id::ID
    val::N
end
CustomNode{N,ID}(id::ID) where {N,ID} = CustomNode{N,ID}(id,N())
function make_node(g::AbstractCustomNGraph{G,CustomNode{N,ID},ID},val::N,id) where {G,N,ID}
    _node_type(g)(id,val)
end
# function make_node(g::AbstractCustomNEGraph{G,CustomNode{N,ID},E,ID},val::N,id) where {G,N,E,ID}
#     _node_type(g)(id,val)
# end
"""
    CustomEdge{E,ID}

A custom node type. Fields:
- `id::ID`
- `val::E`
"""
struct CustomEdge{E,ID}
    src::ID
    dst::ID
    val::E
end
CustomEdge{E,ID}(id1::ID,id2::ID) where {E,ID} = CustomEdge{E,ID}(id1,id2,E())
function make_edge(g::AbstractCustomNEGraph{G,N,CustomEdge{E,ID},ID},u,v,val::E) where {G,N,E,ID}
    _edge_type(g)(
        get_vtx_id(g,get_vtx(g,u)),
        get_vtx_id(g,get_vtx(g,v)),
        val)
end

# """
#     abstract type CustomNGraph{G,N,ID} <: AbstractCustomNGraph{G,CustomNode{N,ID},ID}
#
# Abstract custom graph with nodes of type `CustomNode{N,ID}`.
# """
# abstract type CustomNGraph{G,N,ID} <: AbstractCustomNGraph{G,CustomNode{N,ID},ID} end
# """
#     abstract type CustomNEGraph{G,N,ID} <: AbstractCustomNEGraph{G,CustomEdge{E,ID},CustomNode{N,ID},ID}
#
# Abstract custom graph with nodes of type `CustomNode{N,ID}` and edges of type
# `CustomEdge{E,ID}`.
# """
# abstract type CustomNEGraph{G,N,E,ID} <: AbstractCustomNEGraph{G,CustomEdge{E,ID},CustomNode{N,ID},ID} end
const NGraph{G,N,ID} = CustomNGraph{G,CustomNode{N,ID},ID}
const NEGraph{G,N,E,ID} = CustomNEGraph{G,CustomNode{N,ID},CustomEdge{E,ID},ID}
const NTree{N,ID} = CustomNTree{CustomNode{N,ID},ID}
const NETree{N,E,ID} = CustomNETree{CustomNode{N,ID},CustomEdge{E,ID},ID}

for T in (:CustomNGraph,:CustomNTree)
    @eval begin
        function Base.convert(::Type{X},g::AbstractCustomNGraph{G,N,ID}) where {G,N,ID,X<:$T}
            X(
                deepcopy(get_graph(g)),
                deepcopy(get_nodes(g)),
                deepcopy(get_vtx_map(g)),
                deepcopy(get_vtx_ids(g))
            )
        end
    end
end
for T in (:CustomNEGraph,:CustomNETree)
    @eval begin
        function Base.convert(::Type{X},g::AbstractCustomNEGraph{G,N,E,ID}) where {G,N,E,ID,X<:$T}
            X(
                deepcopy(get_graph(g)),
                deepcopy(get_nodes(g)),
                deepcopy(get_vtx_map(g)),
                deepcopy(get_vtx_ids(g)),
                deepcopy(in_edges(g)),
                deepcopy(out_edges(g)),
            )
        end
    end
end


text:
module Pcap

include("offline.jl")
include("online.jl")
include("headers.jl")

end # module Pcap


text:
using  LinearAlgebra
const  IdentityMatrix = UniformScaling{Bool}

function clip(s::Real, min_thresh::Real, max_thresh::Real, min_new::Real = min_thresh, max_new::Real = max_thresh)
	s = ifelse(s < min_thresh, min_new, ifelse(s > max_thresh, max_new, s))
end

function scaled_norm(E::IdentityMatrix, v::Array, p::Real = 2)
	E.λ ? norm(v, p) : zero(eltype(v))
end

function scaled_norm(E::Diagonal, v::Array{T}, p::Real = 2) where{T}
	if p == 2
		return scaled_norm2(E, v)
	elseif p == Inf
		return scaled_norm_Inf(E, v)
	elseif p == 1
		return scaled_norm1(E, v)
	else
		throw(ArgumentError("bad norm specified"))
	end
end

function scaled_norm2(E::Diagonal, v::Array)
	sum_sq  = zero(eltype(v))
	for i = 1:length(v)
		sum_sq += (E.diag[i] * v[i])^2
	end
	return sqrt(sum_sq)
end

function scaled_norm_Inf(E::Diagonal, v::Array)
	norm  = zero(eltype(v))
	for i = 1:length(v)
		norm = max(norm, abs(E.diag[i] * v[i]))
	end
	return norm
end

function scaled_norm1(E::Diagonal, v::Array)
	norm  = zero(eltype(v))
	for i = 1:length(v)
		norm += abs(E.diag[i] * v[i])
	end
	return norm
end

function col_norms!(v::Array{Tf, 1},
	A::Matrix{Tf};
	reset::Bool = true) where{Tf <: AbstractFloat}

	if(reset) v.= 0 end

	for i = 1:size(A, 2)
		v[i] = max(v[i], norm(view(A, :, i), Inf))
	end
	return v
end

function col_norms!(v::Array{Tf, 1},
	A::SparseMatrixCSC{Tf,Ti}; reset::Bool = true) where{Tf <: AbstractFloat, Ti <: Integer}

	if(reset) v.= 0 end

	for i = 1:A.n
		@inbounds for j = A.colptr[i]:(A.colptr[i + 1] - 1)
		@inbounds v[i] = max(v[i], abs(A.nzval[j]))
	end
	end
	return v
end

function row_norms!(v::Array{Tf, 1},
	A::Matrix{Tf};
	reset::Bool = true) where{Tf <: AbstractFloat}

	if(reset) v.= 0 end

	for i = 1:size(A, 1)
		v[i] = max(v[i], norm(view(A, i, :), Inf))
	end
	return v
end

function row_norms!(v::Array{Tf, 1},
	A::SparseMatrixCSC{Tf, Ti};
	reset::Bool = true) where{Tf <: AbstractFloat, Ti <: Integer}

	if(reset) v.= 0 end

	@inbounds for i = 1:(A.colptr[end] - 1)
		@inbounds v[A.rowval[i]] = max(v[A.rowval[i]], abs(A.nzval[i]))
	end
	return v
end

function lmul!(L::Diagonal, M::SparseMatrixCSC)

	#NB : Same as:  @views M.nzval .*= D.diag[M.rowval]
	#but this way allocates no memory at all and
	#is marginally faster
	m, n = size(M)
	(m == length(L.diag)) || throw(DimensionMismatch())

	@inbounds for i = 1:(M.colptr[end] - 1)
	@inbounds M.nzval[i] *= L.diag[M.rowval[i]]
	end
	return M
end

lmul!(L::IdentityMatrix, M::AbstractMatrix) = L.λ ? M : M .= zero(eltype(M))

function rmul!(M::SparseMatrixCSC, R::Diagonal)

	m, n = size(M)
	(n == length(R.diag)) || throw(DimensionMismatch())

	@inbounds for i = 1:n, j = M.colptr[i]:(M.colptr[i + 1] - 1)
	@inbounds M.nzval[j] *= R.diag[i]
	end
	return M
end

rmul!(M::AbstractMatrix, R::IdentityMatrix) = R.λ ? R : R .= zero(eltype(R))

function lrmul!(L::Diagonal, M::SparseMatrixCSC, R::Diagonal)

	Mnzval  = M.nzval
	Mrowval = M.rowval
	Mcolptr = M.colptr
	Rd      = R.diag
	Ld      = L.diag

	m, n = size(M)
	(m == length(Ld) && n == length(Rd)) || throw(DimensionMismatch())

	@inbounds for i = 1:n, j = Mcolptr[i]:(Mcolptr[i + 1] - 1)
	@inbounds Mnzval[j] *= Ld[Mrowval[j]] * Rd[i]
	end
	return M
end

lrmul!(L::IdentityMatrix,
	M::AbstractMatrix,
	R::IdentityMatrix) = (L.λ && R.λ) ? M : M .= zero(eltype(M))

lrmul!(L::Diagonal,
	M::SparseMatrixCSC,
	R::IdentityMatrix) = R.λ ? lmul!(L, M) : M .= zero(eltype(M))

lrmul!(L::Diagonal,
	M::AbstractMatrix,
	R::Diagonal) = LinearAlgebra.lmul!(L, LinearAlgebra.rmul!(M, R))

lrmul!(L::Diagonal,
	M::AbstractMatrix,
	R::IdentityMatrix) = R.λ ? LinearAlgebra.lmul!(L, M) : M .= zero(eltype(M))

lrmul!(L::IdentityMatrix,
	M::AbstractMatrix,
	R::Diagonal) = L.λ ? LinearAlgebra.rmul!(M, R) : M .= zero(eltype(M))

"""
    symmetrize_upper!(A)

Symmetrizes the matrix A by calculating A = 0.5 * (A + A') but only performs the operation on the upper triangular part.
"""
function symmetrize_upper!(A::AbstractMatrix)
	n = size(A, 1)
	@assert(size(A, 1) == size(A, 2), "Matrix is not square.")
	@inbounds for j in 1:n, i in 1:j
 		A[i, j] = (A[i, j] + A[j, i]) / 2
	end
	nothing
end

"""
    symmetrize_full!(A)

Symmetrizes the matrix A by calculating A = 0.5 * (A + A') and storing the result in-place.
"""
function symmetrize_full!(A::AbstractMatrix)
	n = size(A, 1)
	@assert(size(A, 1) == size(A, 2), "Matrix is not square.")
	@inbounds for j in 1:n, i in 1:j
 		A[i, j] = (A[i, j] + A[j, i]) / 2
 		A[j, i] = A[i, j]
	end
	nothing
end

# this function assumes real symmetric X and only considers the upper triangular part
function is_pos_sem_def(X, tol)
    # set option 'N' to only compute eigenvalues, s is ordered from min to max
   s, U = LAPACK.syevr!('N', 'A', 'U', X, 0.0, 0.0, 0, 0, -1.0);
   return s[1] >= -tol
end

function is_neg_sem_def(X, tol)
    # set option 'N' to only compute eigenvalues, s is ordered from min to max
   s, U = LAPACK.syevr!('N', 'A', 'U', X, 0.0, 0.0, 0, 0, -1.0);
   return s[end] <= tol
end


text:
module Game

using CodeManager

function run()
  println("Create module TEST...")
  #eval(:(module TEST; x=zeros(Float32,128^3*300); end))

  CodeManager.include_module(@__MODULE__, joinpath(@__DIR__,"script.jl"))
  println("Done.")
end

function reload()
  println("Reload module Script in 5 Seconds...")
  sleep(5)
  CodeManager.safe_clean!(@__MODULE__, :Script)
  CodeManager.include_module(@__MODULE__, joinpath(@__DIR__,"script.jl"))
  println("Done.")
end

function cleanUp()
  println("Clean module Script in 5 Seconds...")
  sleep(5)
  CodeManager.safe_clean!(@__MODULE__, :Script)
  println("Done.")
end

end

max_stars_repo_name: Gilga/JLGLAssets

max_issues_repo_name: Gilga/JLGLAssets

max_forks_repo_name: Gilga/JLGLAssets


text:
using GAIO
using ForwardDiff
using StaticArrays
using BenchmarkTools

grid = LinRange(-1, 1, 300)
points = collect(Iterators.product(grid, grid))

f(x) = SA[1/2 - 2*1.4*x[1]^2 + x[2], 0.3*x[1]]

domain = Box((0.0, 0.0), (1.0, 1.0))
partition = RegularPartition(domain, 18)
point = ntuple(i->rand(2), 10000)
B = partition[point]

g1 = PointDiscretizedMap(f, domain, points)
g2 = AdaptiveBoxMap(f, domain)

@time g1B = g1(B)
g2B = g2(B)

# and for the Lorenz system

@inline function rk4(f, x, τ)
    τ½ = τ/2

    k = f(x)
    dx = @. k/6

    k = f(@. x+τ½*k)
    dx = @. dx + k/3

    k = f(@. x+τ½*k)
    dx = @. dx + k/3

    k = f(@. x+τ*k)
    dx = @. dx + k/6

    return @. x+τ*dx
end

function lorenz_v(x)
    s = 10.0
    rh = 28.0
    b = 0.4
    return s * (x[2] - x[1]), rh * x[1] - x[2] - x[1] * x[3], x[1] * x[2] - b * x[3]
end

function lorenz_f(x)
    h = 0.01
    n = 20

    for i = 1:n
        x = rk4(lorenz_v, x, h)
    end

    return x
end

grid = LinRange(-1, 1, 6)
points = collect(Iterators.product(grid, grid, grid))

domain = Box((0.0, 0.0, 27.0), (30.0, 30.0, 40.0))
partition = BoxPartition(domain, 27)

rh = 28.0
b = 0.4
x0 = (sqrt(b*(rh-1)), sqrt(b*(rh-1)), rh-1)

boxset = partition[x0]

g_adaptive = AdaptiveBoxMap(lorenz_f, domain)
g_points = PointDiscretizedMap(lorenz_f, domain, points)

point = ntuple(i->60.0*rand(3).-30.0, 1000000);
B = partition[point]; @show length(B)

@time g1B = g_points(B); @show length(g1B)
gaB = g_adaptive(B); @show length(gaB)
@show g_adaptive.image_points(point, point)


text:
#  SPI Read function.  This version returns at single byte at the address.
#  This will need to form an Avalon packet, just like the write function.

#module spi_read
#import SPIManager

#=
What is the data flow of the read and write functions?
This process is pretty strange, since Julia is calling C library functions.
Also there are functions to initialize and call the FTDI SPI device mixed in.
I will ignore the FTDI in this commentary.
The C code, which is derived from the Altera example code, is built into a
"shared library".  The library has functions which translate read and write
requests into SPI data.  It is a bit confusing, as there are "Avalon Packets",
physical Avalon transfers, and then the SPI data itself, which is what is
sent to the FTDI device.

The C library creates its own buffers (arrays) which contain the SPI data.
Julia can "wrap" these buffers with the "unsafe_wrap" function, which is the
mechanism for Julia to access the C library generated arrays.

Since the C buffer/array is "malloced", Julia needs to get the pointer.
This is done with the get_send_byte_buffer_pointer() function, which was added
to the Altera example code for this purpose.

Also needed is the size of the buffer.  This size is variable, as there are
substitutions made by the Altera code to make the physical Avalon transfer, due
to having to work around characters used for special purposes, for example
end of buffer marking.  This is some quite funky C code.

With the pointer and buffer size available, the "unsafe_wrap" function opens
the door for Julia to use the C buffer to write to the FTDI SPI device.

Writing is fairly straightforward, as the physical Avalon transfer is created,
and this buffer is transmitted by an FTDI SPI write.  So this requires and
address and data.

The read, however, is a little trickier.  The Avalon documentation is not clear
on this.  This should require an address, and if more than one byte is being
read, the length.

=#

export readbyte, printit, readavalon
#  Path to the shared library:  /mnt/ds216/embedded_design/FPGA/avalon_bus_stuff/data_to_avalon_packets_library/build/default/libavalonpacket.so.  Non-sequential read is 0x00, Sequential read is 0x04.
readavalon(addr, size, buff) = ccall((:transaction_channel_read, "/mnt/ds216/embedded_design/FPGA/avalon_bus_stuff/data_to_avalon_packets_library/build/default/libavalonpacket"), UInt32, (UInt32, UInt32, Ptr{UInt8}, UInt32), addr, size, buff, 0x04)
#  Note, the buff parameter is meaningless here.  It can be a dummy buffer, and it won't be useful.
#  In the original project, the SPI command would be used to populate this buffer.  Can't do that here, it is more of a multi-step process.
#packet_length() = ccall((:get_xmit_packet_length, "/mnt/ds216/embedded_design/FPGA/avalon_bus_stuff/data_to_avalon_packets_library/build/default/libavalonpacket"), UInt32,(),)

#  The function should return the retrieved data.
function readbyte(
            address
       )

    # Open the channel to the SPI device.
    SPIManager.SPI_GetNumChannels()
    SPIManager.SPI_InitChannel()

    avalonbuffer = zeros(UInt8, 100)  #  This should be more than adequate for 1 byte!

    #sizeToTransfer = packet_length()
    # println("transfer_data length is $transfer_data_length")
    sizeTransfered = Ref{Culong}()

    #  Need an equivalent to writeavalon, call it readavalon.
    bogusbuffer = zeros(UInt8, 1)
    # Transform to Avalon packet:
    #println("Prior to readavalon")
    #  This function is attempting to free invalid pointer.
    #  What is the buffer in the read case???
    readavalon(address, 1, bogusbuffer)
    #println("After readavalon")
    #  Now, get the length of the resultant Avalon packet:
    sizeToTransfer = SPIManager.packet_length()
    println("transfer_data length is $sizeToTransfer")
    buffersize = sizeToTransfer  + 12  #  This is the write packet plus maximum response length.
    #  This is the buffer to be transmitted on the SPI.
    avalonbuffer = unsafe_wrap(Array, SPIManager.get_send_byte_buffer_pointer(),buffersize, own = false)
    #  Create another buffer for the received data:
    readbuffer = zeros(UInt8, 100)
    #  Now, print out the buffer temporarily for debugging.
    for i in 1:12
      #println("avalonbuffer[$i] = $avalonbuffer[i]")
      Printf.@printf("avalonbuffer[%d] = %x\n", i, avalonbuffer[i])
    end

# FT_STATUS SPI_Read(FT_HANDLE handle, uint8 *buffer, uint32 sizeToTransfer, uint32 *sizeTransferred, uint32 transferOptions)

    #  This function needs to write to some buffer.                                                                             Culonglong               Ptr{UInt8}  Ptr(UInt8)    UInt32      Ref(Culong)     UInt32
       spiread = ccall((:SPI_ReadWrite, "libMPSSE"), UInt32, (Culonglong, Ptr{UInt8}, Ptr{UInt8}, UInt32, Ref{Culong}, UInt32), SPIManager.actualhandle, readbuffer, avalonbuffer, buffersize, sizeTransfered, 0)
        #  Return the first byte of the buffer?
       println("spiread = $spiread")

   #    for i in 1:24
   #    println("readbuffer is $readbuffer[i]")
   #    end

       println("The sizeTransfered is $sizeTransfered[]")
       for i in 1:100
         Printf.@printf("readbuffer[%d] = %x\n", i, readbuffer[i])
       end
       SPIManager.SPI_CloseChannel()
    end

max_stars_repo_name: Greg-R/spiavalonfpga

max_issues_repo_name: Greg-R/spiavalonfpga

max_forks_repo_name: Greg-R/spiavalonfpga


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

export PolyPairBasis

# TODO: allow PairBasis with different Pr for each z, z' combination

struct PolyPairBasis{TJ, NZ} <: IPBasis
   J::TJ
   zlist::SZList{NZ}
   bidx0::SMatrix{NZ,NZ,Int}
end

fltype(pB::PolyPairBasis) = fltype(pB.J)

Base.length(pB::PolyPairBasis) = length(pB.J) * (numz(pB) * (numz(pB) + 1)) ÷ 2
Base.length(pB::PolyPairBasis, z0::AtomicNumber) = length(pB, z2i(pB, z0))
Base.length(pB::PolyPairBasis, iz0::Integer) = length(pB.J)

zlist(pB::PolyPairBasis) = pB.zlist

function scaling(pB::PolyPairBasis, p)
   ww = zeros(Float64, length(pB))
   for iz0 = 1:numz(pB), iz = 1:numz(pB)
      idx0 = _Bidx0(pB, iz0, iz)
      for n = 1:length(pB.J)
         # TODO: very crude, can we do better?
         #       -> need a proper H2-orthogonbality?
         ww[idx0+n] = n^p
      end
   end
   return ww
end

PolyPairBasis(J::ScalarBasis, species) =
   PolyPairBasis( J, ZList(species; static=true) )

PolyPairBasis(J::ScalarBasis, zlist::SZList) =
   PolyPairBasis(J, zlist, get_bidx0(J, zlist))

function get_bidx0(J, zlist::SZList{NZ}) where {NZ}
   NJ = length(J)
   bidx0 = fill(zero(Int), (NZ, NZ))
   i0 = 0
   for i = 1:NZ, j = i:NZ
      bidx0[i,j] = i0
      bidx0[j,i] = i0
      i0 += NJ
   end
   return SMatrix{NZ, NZ, Int}(bidx0...)
end

==(B1::PolyPairBasis, B2::PolyPairBasis) =
      (B1.J == B2.J) && (B1.zlist == B2.zlist)

JuLIP.cutoff(pB::PolyPairBasis) = cutoff(pB.J)

write_dict(pB::PolyPairBasis) = Dict(
      "__id__" => "ACE_PolyPairBasis",
          "Pr" => write_dict(pB.J),
       "zlist" => write_dict(pB.zlist) )

read_dict(::Val{:SHIPs_PolyPairBasis}, D::Dict) =
   read_dict(Val{:ACE_PolyPairBasis}(), D)

read_dict(::Val{:ACE_PolyPairBasis}, D::Dict) =
      PolyPairBasis( read_dict(D["Pr"]), read_dict(D["zlist"]) )

alloc_temp(pB::PolyPairBasis, args...) = (
              J = alloc_B(pB.J),
          tmp_J = alloc_temp(pB.J)  )

alloc_temp_d(pB::PolyPairBasis, args...) =  (
             J = alloc_B( pB.J),
         tmp_J = alloc_temp(pB.J),
            dJ = alloc_dB(pB.J),
        tmpd_J = alloc_temp_d(pB.J)  )

"""
compute the zeroth index of the basis corresponding to the potential between
two species zi, zj; as precomputed in `PolyPairBasis.bidx0`
"""
_Bidx0(pB, zi, zj) = pB.bidx0[ z2i(pB, zi), z2i(pB, zj) ]
_Bidx0(pB, i::Integer, j::Integer) = pB.bidx0[ i, j ]

function energy(pB::PolyPairBasis, at::Atoms{T}) where {T}
   E = zeros(T, length(pB))
   tmp = alloc_temp(pB)
   for (i, j, R) in pairs(at, cutoff(pB))
      r = norm(R)
      evaluate!(tmp.J, tmp.tmp_J, pB.J, r)
      idx0 = _Bidx0(pB, at.Z[i], at.Z[j])
      for n = 1:length(pB.J)
         E[idx0 + n] += 0.5 * tmp.J[n]
      end
   end
   return E
end

function forces(pB::PolyPairBasis, at::Atoms{T}) where {T}
   F = zeros(JVec{T}, length(at), length(pB))
   tmp = alloc_temp_d(pB)
   for (i, j, R) in pairs(at, cutoff(pB))
      r = norm(R)
      evaluate_d!(tmp.J, tmp.dJ, tmp.tmpd_J, pB.J, r)
      idx0 = _Bidx0(pB, at.Z[i], at.Z[j])
      for n = 1:length(pB.J)
         F[i, idx0 + n] += 0.5 * tmp.dJ[n] * (R/r)
         F[j, idx0 + n] -= 0.5 * tmp.dJ[n] * (R/r)
      end
   end
   return [ F[:, iB] for iB = 1:length(pB) ]
end

function virial(pB::PolyPairBasis, at::Atoms{T}) where {T}
   V = zeros(JMat{T}, length(pB))
   tmp = alloc_temp_d(pB)
   for (i, j, R) in pairs(at, cutoff(pB))
      r = norm(R)
      evaluate_d!(tmp.J, tmp.dJ, tmp.tmpd_J, pB.J, r)
      idx0 = _Bidx0(pB, at.Z[i], at.Z[j])
      for n = 1:length(pB.J)
         V[idx0 + n] -= 0.5 * (tmp.dJ[n]/r) * R * R'
      end
   end
   return V
end


text:
println("First line of the script")
time_start = time()
using Test
using Logging
using SparseArrays
using ExSup
using Pio3d
using JLD2
@info "Before Distributed"
using Distributed
if nprocs() == 1
    addprocs(3)
end
# using Revise
# using ExSup
using LarSurf
# @everywhere using LarSurf

fn = "exp_surface_extraction_ircad_times.csv"

@info "before everywhere using"
@info "time from start: $(time()-time_start) [s]"
# using Plasm
# @everywhere using LarSurf
# @everywhere using Distributed

@info "after everywhere using, time from start: $(time()-time_start) [s]"

block_size = [64, 64, 64]

data_id = 1
show = false
taubin = true
taubin_n = 5
taubin_p1 = 0.4
taubin_p2 = -0.2

stepz = 1
stepxy = 1
do_crop = false
cropx = 200
cropy = 200
cropz = 400
# stepxy = 4
# block_size = [128, 128, 128]
# block_size = [128, 128, 128]
# block_size = [32, 32, 32]
# data_size1 = 128
# data_size1 = 256
# data_size1 = 512

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

data = LarSurf.Experiments.report_init_row(@__FILE__)
LarSurf.set_time_data(data)

# data["nprocs"] = nprocs()
# data["fcn"] = String(Symbol(fcni))
data["block size"] = block_size[1]

data["using done"] = time()-time_start
# segmentation = LarSurf.data234()
@info "Generate data..."
@info "time from start: $(time() - time_start) [s]"
mask_labels=["liver", "portalvein"]
# pth = Pio3d.datasets_join_path("medical/orig/3Dircadb1.$data_id/MASKS_DICOM/liver")

mask_label = "nrn"
threshold = 7000
# for mask_label in mask_labels

	pth = Pio3d.datasets_join_path("medical/orig/jatra_mikro_data/Nejlepsi_rozliseni_nevycistene")
	pth = Pio3d.datasets_join_path("medical/processed/corrosion_cast/nrn10.pklz")
	datap = Pio3d.read3d(pth)
	data3d_full = datap["data3d"]
	@info "raw data size=$(size(data3d_full))"
	data3d_full = data3d_full[1:stepz:end, 1:stepxy:end, 1:stepxy:end]
	if do_crop
		data3d_full = data3d_full[1:cropx, 1:cropy, 1:cropz]
	end
	    # round(size(data3d_full, 1) / target_size1)
	    # return data3d_full
	# segmentation = convert(Array{Int8, 3}, data3d_full .> threshold)
	segmentation = Int8.(data3d_full .> threshold)

	sz = [
		size(data3d_full, 1),
	 	size(data3d_full, 2),
	 	size(data3d_full, 3)
		]
	# @info "nvoxels=$(nvoxels)"
	voxelsize_mm = datap["voxelsize_mm"]
	voxelsize_mm[1] = voxelsize_mm[1] * stepz
	voxelsize_mm[2] = voxelsize_mm[2] * stepxy
	voxelsize_mm[3] = voxelsize_mm[3] * stepxy
	@info "voxelsize mm = $(voxelsize_mm), size = $(sz)"
	data = LarSurf.Experiments.report_add_data_info(data, segmentation, voxelsize_mm)
	# segmentation = LarSurf.generate_cube(data_size1; remove_one_pixel=true)
	@info "==== using done, data generated time from start: $(time() - time_start) [s]"
	data["data generated"] = time() - time_start

	@info "Setup..."
	setup_time = @elapsed LarSurf.lsp_setup(block_size;reference_time=time_start)
	println("setup time: $setup_time")
	@info "==== setup done, time from start: $(time() - time_start) [s]"
	data["setup done"] = time() - time_start
	# for wid in workers()
	#     # println("testing on $wid")
	#     ftr = @spawnat wid LarSurf._single_boundary3
	#     @test fetch(ftr) != nothing
	# end

	# @debug "Setup done"
	tmd = @timed larmodel = LarSurf.lsp_get_surface(segmentation; voxelsize=voxelsize_mm)
	val, tm, mem, gc = tmd
	println("Total time per $(mask_label): $tm")
	@info "==== time from start: $(time() - time_start) [s]"
	data["finished"] = time() - time_start
	ExSup.datetime_to_dict!(data)
	ExSup.add_to_csv(data, fn)

	V, FV = larmodel
	FVtri = LarSurf.triangulate_quads(FV)

	@JLD2.save "ircad_$(mask_label).jld2" V FV
	using ViewerGL
	println("=== ViewerGL readed")
	ViewerGL.VIEW([
	    ViewerGL.GLGrid(V, FVtri, ViewerGL.Point4d(1,1,1,0.1))
		ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
	])
	# FV = FVtri
	@JLD2.save "ircad_$(mask_label)_tri.jld2" V FVtri
	@info "smoothing"
	objlines = LarSurf.Lar.lar2obj(V, FVtri, "nrn_tri.obj")

	if taubin
		t = @elapsed Vs = LarSurf.Smoothing.smoothing_FV_taubin(V, FV, taubin_p1, taubin_p2, taubin_n)
		@info "smoothing taubin time", t
		# t = @elapsed FVtri = LarSurf.triangulate_quads(FV)

		if show
			ViewerGL.VIEW([
			    ViewerGL.GLGrid(Vs, FVtri, ViewerGL.Point4d(1, 0, 1, 0.1))
				ViewerGL.GLAxis(ViewerGL.Point3d(-1, -1, -1),ViewerGL.Point3d(1, 1, 1))
			])
		end
	else
		t = @elapsed Vs = LarSurf.Smoothing.smoothing_FV(V, FVtri, 0.6, 3)
		@info "smoothing time", t
	end
	@info "Smoothing numer of Vs: $(size(Vs))"
	@JLD2.save "ircad_$(mask_label)_sm.jld2" Vs FVtri
	# @JLD2.save "liver01tri.jld2" V FVtri
	objlines = LarSurf.Lar.lar2obj(Vs, FVtri, "nrn_tri_taubin.obj")
# end
# Plasm.view(val)


text:
@testset "Constructor" begin

using Base.Test
using DataValueArrays

# test Inner Constructor
@test_throws ArgumentError DataValueArray([1, 2, 3, 4], [true, false, true])

# test (A::AbstractArray, m::Array{Bool}) constructor
v = [1, 2, 3, 4]
dv = DataValueArray(v, fill(false, size(v)))

m = [1 2; 3 4]
dm = DataValueArray(m, fill(false, size(m)))

t = Array{Int}(2, 2, 2)
t[1:2, 1:2, 1:2] = 1

dt = DataValueArray(t, fill(false, size(t)))
dv = DataValueArray(v, [false, false, false, false])

y2 = DataValueArray([1, 2, 3, 4, 5, 6],
                    [true, false, false, false, false ,false])
@test isa(y2, DataValueVector{Int})
@test y2.isnull[1]

# test (::AbstractArray) constructor
dv = DataValueArray(v)
@test isa(dv, DataValueVector{Int})

y = DataValueArray([1, 2, 3, 4, 5, 6])
@test isa(y, DataValueVector{Int})

z = DataValueArray(1.:6.)
@test isa(z, DataValueVector{Float64})

# test (::Type{T}, dims::Dims) constructor
u1 = DataValueArray(Int, (5, ))
u2 = DataValueArray(Int, (2, 2))
u3 = DataValueArray(Int, (2, 2, 2))
@test isa(u1, DataValueVector{Int})
@test isa(u2, DataValueMatrix{Int})
@test isa(u3, DataValueArray{Int, 3})

# test (::Type{T}, dims::Int...) constructor
x1 = DataValueArray(Int, 2)
x2 = DataValueArray(Int, 2, 2)
x3 = DataValueArray(Int, 2, 2, 2)
@test isa(x1, DataValueVector{Int})
@test isa(x2, DataValueMatrix{Int})
@test isa(x3, DataValueArray{Int, 3})

# test DataValueArray{T}(dims::Dims)
d1, d2 = rand(1:100), rand(1:100)
X1 = DataValueArray{Int}((d1,))
X2 = DataValueArray{Int}((d1, d2))
@test isequal(X1, DataValueArray(Array{Int}((d1,)), fill(true, (d1,))))
@test isequal(X2, DataValueArray(Array{Int}((d1, d2)), fill(true, (d1, d2))))
for i in 1:5
    m = rand(3:5)
    dims = tuple([ rand(1:5) for i in 1:m ]...)
    X3 = DataValueArray{Int}(dims)
    @test isequal(X3, DataValueArray(Array{Int}(dims), fill(true, dims)))
end

# test DataValueArray{T}(dims::Int...)
d1, d2 = rand(1:100), rand(1:100)
X1 = DataValueArray{Int}(d1)
X2 = DataValueArray{Int}(d1, d2)
@test isequal(X1, DataValueArray(Array{Int}(d1), fill(true, d1)))
@test isequal(X2, DataValueArray(Array{Int}(d1, d2), fill(true, d1, d2)))
for i in 1:5
    m = rand(3:5)
    dims = [ rand(1:5) for i in 1:m ]
    X3 = DataValueArray{Int}(dims...)
    @test isequal(X3, DataValueArray(Array{Int}(dims...), fill(true, dims...)))
end

# test DataValueArray{T}(dims::Int...)
d1, d2 = rand(1:100), rand(1:100)
X1 = DataValueArray{Int,1}(d1)
X2 = DataValueArray{Int,2}(d1, d2)
@test isequal(X1, DataValueArray(Array{Int,1}(d1), fill(true, d1)))
@test isequal(X2, DataValueArray(Array{Int,2}(d1, d2), fill(true, d1, d2)))
for i in 1:5
    m = rand(3:5)
    dims = [ rand(1:5) for i in 1:m ]
    X3 = DataValueArray{Int,length(dims)}(dims...)
    @test isequal(X3, DataValueArray(Array{Int}(dims...), fill(true, dims...)))
end

# test (A::AbstractArray, ::Type{T}, ::Type{U}) constructor
z = DataValueArray([1, nothing, 2, nothing, 3], Int, Void)
@test isa(z, DataValueVector{Int})
@test z.isnull[2]
@test z.isnull[4]

# test (A::AbstractArray, ::Type{T}, na::Any) constructor
Z = DataValueArray([1, "na", 2, 3, 4, 5, "na"], Int, "na")
@test isa(Z, DataValueVector{Int})
@test Z.isnull == [false, true, false, false, false, false, true]

Y = DataValueArray([1, nothing, 2, 3, 4, 5, nothing], Int, Void)
@test isequal(Y, Z)

# test DataValueArray{T}()
X = DataValueArray{Int}()
@test isequal(size(X), ())

# test conversion from arrays, arrays of DataValues and DataValueArrays
miss1 = [false,false,false,false]
miss2 = [false,false,false,true]
for (a, miss) in zip(([1, 2, 3, 4],
                        # 1:4, # Currently does not work on 0.4, cf. JuliaLang/julia#16265
                        DataValue{Int}[DataValue(1), DataValue(2), DataValue(3), DataValue()],
                        DataValueArray(DataValue{Int}[DataValue(1), DataValue(2), DataValue(3), DataValue()])),
                        (miss1, miss2, miss2))
    @test isa(DataValueArray(a), DataValueArray{Int,1})
    @test isequal(DataValueArray(a), DataValueArray{Int,1}([1,2,3,4],miss))
    @test isa(DataValueArray{Int}(a), DataValueArray{Int,1})
    @test isequal(DataValueArray{Int}(a), DataValueArray{Int,1}([1,2,3,4],miss))
    @test isa(DataValueArray{Float64}(a), DataValueArray{Float64,1})
    @test isequal(DataValueArray{Float64}(a), DataValueArray{Float64,1}([1.0,2.0,3.0,4.0],miss))
    @test isa(DataValueArray{Int,1}(a), DataValueArray{Int,1})
    @test isequal(DataValueArray{Int,1}(a), DataValueArray{Int,1}([1,2,3,4],miss))
    @test isa(DataValueArray{Float64,1}(a), DataValueArray{Float64,1})
    @test isequal(DataValueArray{Float64,1}(a), DataValueArray{Float64,1}([1.0,2.0,3.0,4.0],miss))

    @test isa(convert(DataValueArray, a), DataValueArray{Int,1})
    @test isequal(convert(DataValueArray, a), DataValueArray{Int,1}([1,2,3,4],miss))
    @test isa(convert(DataValueArray{Int}, a), DataValueArray{Int,1})
    @test isequal(convert(DataValueArray{Int}, a), DataValueArray{Int,1}([1,2,3,4],miss))
    @test isa(convert(DataValueArray{Float64}, a), DataValueArray{Float64,1})
    @test isequal(convert(DataValueArray{Float64}, a), DataValueArray{Float64,1}([1.0,2.0,3.0,4.0],miss))
    @test isa(convert(DataValueArray{Int,1}, a), DataValueArray{Int,1})
    @test isequal(convert(DataValueArray{Int,1}, a), DataValueArray{Int,1}([1,2,3,4],miss))
    @test isa(convert(DataValueArray{Float64,1}, a), DataValueArray{Float64,1})
    @test isequal(convert(DataValueArray{Float64,1}, a), DataValueArray{Float64,1}([1.0,2.0,3.0,4.0],miss))
end

# test conversion from Array{DataValue} without element type (issue #145)
@test isa(DataValueArray(DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Any})
@test isa(DataValueArray{Int}(DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Int})
@test isa(DataValueArray{Float64}(DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Float64})

@test isa(convert(DataValueArray, DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Any})
@test isa(convert(DataValueArray{Int}, DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Int})
@test isa(convert(DataValueArray{Float64}, DataValue[DataValue(1), DataValue(2), DataValue(3), DataValue()]), DataValueArray{Float64})

# converting a DataValueArray to unqualified type DataValueArray should be no-op
m = rand(10:100)
A = rand(m)
M = rand(Bool, m)
X = DataValueArray(A, M)
@test X === convert(DataValueArray, X)

end


text:
using DataFrames, UnicodePlots, Test, Random, StatsBase

input = deserialize(joinpath(@__DIR__, "names_df.jls"))
df = input[:df]
col_dict = input[:col_dict]

# naive O(n^2) algorithm
function kendall_tau_distance(σ1, σ2)
    n = length(σ1)
    n == length(σ2) || throw(DimensionMismatch())
    D = 0
    for j = 1:n
        for i = 1:(j-1)
            D += ( (σ1[i] < σ1[j]) & (σ2[i] > σ2[j]) ) |
                ( (σ1[i] > σ1[j]) & (σ2[i] < σ2[j]) )
        end
    end
    return D
end

# Example borrowed from http://arxiv.org/abs/1905.02752
@test kendall_tau_distance([2, 4, 1, 3], [4, 1, 3, 2]) == 5
π = randperm(4)
@test kendall_tau_distance([2, 4, 1, 3][π], [4, 1, 3, 2][π]) == 5

function normalized_kendall_tau_distance(σ1, σ2)
    n = length(σ1)
    2*kendall_tau_distance(σ1, σ2) / (n * (n-1))
end
kendall_tau_coeff(x,y) = 1 - normalized_kendall_tau_distance(x,y)

function corrs(coeff, suffix)
    cols =  [v for (k, v) in pairs(col_dict) if k[2] == suffix ]
    [ coeff(df[!, c1] , df[!, c2])  for c1 in cols, c2 in cols]
end

all_cols = collect(values(col_dict))
sort!(all_cols)
corrs(coeff) = [ coeff(df[!, s1] , df[!, s2])   for s1 in all_cols, s2 in all_cols  ]

suffix = :sqrt_normalized
for coeff in (corkendall, corspearman, kendall_tau_coeff)
    cols =  [v for (k, v) in pairs(col_dict) if k.normalization == suffix ]
    @info cols
    @info "$coeff with $suffix" corrs(coeff, suffix) heatmap(corrs(coeff, suffix))
end

for coeff in (corkendall, corspearman, kendall_tau_coeff)
    @info all_cols
    @info "$coeff" corrs(coeff) heatmap(corrs(coeff))
end

for (k, v) in pairs(col_dict)
    k.normalization == :sqrt_normalized || continue
    if !isempty(k.params)
        k.params.ϵ > .5 && continue
        k.params.ρ > 5 && continue
    end
    n = size(df, 1)

    # Shuffle dataframe
    p = randperm(n)
    for c in propertynames(df)
        permute!(df[!, c], p)
    end

    top5 = sort!(df, v)[1:5, ["name1", "name2", v]]
    num = count(x -> x <= df[5, v], df[!, v])
    @info "$k; showing random top 5 with $(num) candidates to show" top5
end


text:
using Documenter, SwitchOnSafety

makedocs(
    sitename = "SwitchOnSafety",
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    # See https://github.com/JuliaOpt/JuMP.jl/issues/1576
    strict = true,
    pages = [
        "Index" => "index.md",
        "Invariant Sets" => "invariant.md",
        "Joint Spectral Radius" => "jsr.md"
    ]
)

deploydocs(
    repo   = "github.com/blegat/SwitchOnSafety.jl.git"
)


text:
struct Dual{V<:Real} <: Real
    value::V
    partials::AbstractVector{V}
end
partials(d::Dual) = d.partials
value(d::Dual) = d.value

Dual(v::V,ps...) where V<:Real = Dual{V}(v,collect(ps))

function Base.show(io::IO, d::Dual{V}) where V
    print(io, "Dual(", d.value)
    for p in d.partials
        print(io,",",p)
    end
    print(io, ")")
end

function jacobian(f, x::AbstractVector{T}) where T
    duals = map(1:length(x)) do i
        ps = map(1:length(x)) do j
            i==j ? T(1) : T(0)
        end
        Dual(x[i], ps)
    end
    mapreduce(r->reshape(partials(r),1,:), vcat, f(duals))
end

##########  RULES  #############################################################

Base.:+(a::Dual, b::Dual) = Dual(a.value + b.value, a.partials + b.partials)
Base.:*(a::Dual, b::Dual) = Dual(a.value * b.value, b.value * a.partials + a.value * b.partials)

function cprod(x::AbstractVector)
    y = similar(x)
    y[1] = x[1]
    for i in 2:length(y)
        y[i] = y[i-1]*x[i]
    end
    y
end

x = collect(1:4)
#csum(x) |> display

dx = [Dual(x[1], 1, 0, 0),
      Dual(x[2], 0, 1, 0),
      Dual(x[3], 0, 0, 1)]
cprod(dx) |> display
csum(dx)

jacobian(cprod, x)


text:
using ABM
using Random
using Distributions
mutable struct BoidData
    name::String
    pos::Real2D
    orientation::Real2D
end
function fstep(state::SimState, agent::Agent)
end
pos = Real2D(rand(Uniform(0, 100)), rand(Uniform(0, 100)))
d = BoidData("Boid", pos, Real2D(rand(Uniform(0, 100)),rand(Uniform(0, 100))))
a = Agent(fstep,d)

dict1 = Dict{Int64,Int64}()
array1 = Array{Int64,1}()

dict2 = Dict{Agent,Int64}()
array2 = Array{Agent,1}()

set1 = Set{Int64}()
set2 = Set{Agent}()

println("------------------------------")

@time for i in 1:1000000
    pos = Real2D(rand(Uniform(0, 100)), rand(Uniform(0, 100)))
    d = BoidData("Boid", pos, Real2D(rand(Uniform(0, 100)),rand(Uniform(0, 100))))
    a = Agent(fstep,d)
    push!(array2, a)
end

@time for i in 1:1000000
    pos = Real2D(rand(Uniform(0, 100)), rand(Uniform(0, 100)))
    d = BoidData("Boid", pos, Real2D(rand(Uniform(0, 100)),rand(Uniform(0, 100))))
    a = Agent(fstep,d)
    push!(dict2, a=>i)
end

@time for i in 1:1000000
    pos = Real2D(rand(Uniform(0, 100)), rand(Uniform(0, 100)))
    d = BoidData("Boid", pos, Real2D(rand(Uniform(0, 100)),rand(Uniform(0, 100))))
    a = Agent(fstep,d)
    push!(set2, a)
end

@time for x in values(dict2)
end

@time for x in values(set2)
end

@time for x in array2
end


text:
using PyPlot, PyCall
using Distributions
import StatsFuns.logsumexp
PyDict(matplotlib["rcParams"])["mathtext.fontset"] = "cm"
PyDict(matplotlib["rcParams"])["mathtext.rm"] = "serif"
PyDict(matplotlib["rcParams"])["lines.linewidth"] = 1.5
PyDict(matplotlib["rcParams"])["font.family"] = "TakaoPGothic"

function expt(a, b, sigma, Y, X, N_s)
    S = rand(Gamma(a, 1.0/b), N_s)
    C = mean([exp(sum(logpdf.(Normal(s, sigma), Y))) for s in S])
    curve = [exp(sum(logpdf.(Normal(mu, sigma), Y))) * pdf(Gamma(a, 1.0/b), mu) for mu in X]
    m = mean([s*exp(sum(logpdf.(Normal(s, sigma), Y)))/C for s in S])
    v = mean([(s-m)^2 * exp(sum(logpdf.(Normal(s, sigma), Y)))/C for s in S])
    return curve/C, m, v
end

X = linspace(-5, 10, 1000)

a = 2.0
b = 2.0
mu = 1.0
sigma=1.0

# data
N = 10
Y = rand(Normal(mu, sigma), N)

# calc posterior
N_s = 100000
posterior, m, v = expt(a, b, sigma, Y, X, N_s)

a_h = m^2 / v
b_h = m / v

figure()
plot(X, pdf(Normal(mu,sigma), X))
plot(X, pdf(Gamma(a,1.0/b), X))
plot(X, posterior)
plot(X, pdf(Gamma(a_h,1.0/b_h), X))
plot(Y, 0.02*ones(N), "o")
legend(["generator", "prior", "posterior", "approx", "samples"])
#legend(["データ生成分布", "事前分布", "事後分布", "近似分布", "データ"], fontsize=12)
xlim([-3, 6])
ylim([0, 1.8])


text:
struct UniformContinuous <: Distribution{Float64} end

const uniform_continuous = UniformContinuous()

"""
    uniform(low::Real, high::Real)

Sample a `Float64` from the uniform distribution on the interval [low, high].
"""
const uniform = uniform_continuous

function logpdf(::UniformContinuous, x::Real, low::Real, high::Real)
    (x >= low && x <= high) ? -log(high-low) : -Inf
end

function logpdf_grad(::UniformContinuous, x::Real, low::Real, high::Real)
    inv_diff = 1. / (high-low)
    (0., inv_diff, -inv_diff)
end

function random(::UniformContinuous, low::Real, high::Real)
    rand() * (high - low) + low
end

(::UniformContinuous)(low, high) = random(UniformContinuous(), low, high)
is_discrete(::UniformContinuous) = false

has_output_grad(::UniformContinuous) = true
has_argument_grads(::UniformContinuous) = (true, true)

export uniform_continuous, uniform


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

mutable struct VirtualNetworkPeering <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function VirtualNetworkPeering(;id=nothing)
        o = new()
        validate_property(VirtualNetworkPeering, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type VirtualNetworkPeering

const _property_map_VirtualNetworkPeering = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_VirtualNetworkPeering = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ VirtualNetworkPeering }) = collect(keys(_property_map_VirtualNetworkPeering))
Swagger.property_type(::Type{ VirtualNetworkPeering }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkPeering[name]))}
Swagger.field_name(::Type{ VirtualNetworkPeering }, property_name::Symbol) =  _property_map_VirtualNetworkPeering[property_name]

function check_required(o::VirtualNetworkPeering)
    true
end

function validate_property(::Type{ VirtualNetworkPeering }, name::Symbol, val)
end


text:
using BatteryCost
using Test

@testset "cost calculation" begin
    cell1 = cell()
    cost1 = BatteryCost.cost_default()
    cell3,cost3 = BatteryCost.convert_all(cell1,cost1,BatteryCost.mult)
    @test cost_calc(cell3,cost3)==(209.68524992517177, 1708.1761931921092)
end

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
############################################################
# joSAparallelToggles - parallel V/MV conversion operators #
############################################################

# includes
include("joAbstractSAparallelToggleOperator/constructors.jl")
include("joAbstractSAparallelToggleOperator/base_functions.jl")
include("joAbstractSAparallelToggleOperator/extra_functions.jl")


text:
using SparseInterp, Base.Test

G = NGrid([8,8,8,8,8])
A = rand(length(G))

@test all(isapprox.(G(A, G.grid), A))


text:
# Autogenerated wrapper script for Gumbo_jll for i686-w64-mingw32
export libgumbo

JLLWrappers.@generate_wrapper_header("Gumbo")
JLLWrappers.@declare_library_product(libgumbo, "libgumbo-1.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgumbo,
        "bin\\libgumbo-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Parameters
using ModelParameters

abstract type IHACRESNode <: NetworkNode end

Base.@kwdef mutable struct BilinearNode{P, A<:Real} <: IHACRESNode
    @network_node

    # https://wiki.ewater.org.au/display/SD41/IHACRES-CMD+-+SRG
    d::P = Param(200.0, bounds=(10.0, 550.0))  # flow threshold
    d2::P = Param(2.0, bounds=(0.0001, 10.0))   # flow threshold2
    e::P = Param(1.0, bounds=(0.1, 1.5))  # temperature to PET conversion factor
    f::P = Param(0.8, bounds=(0.01, 3.0))  # plant stress threshold factor (multiplicative factor of d)
    a::P = Param(0.9, bounds=(0.1, 10.0))  # quickflow storage coefficient == (1/tau_q)
    b::P = Param(0.1, bounds=(1e-3, 0.1))  # slowflow storage coefficent == (1/tau_s)

    # These changes require ihacres_nim to be changed as well...
    # a::A = Param(0.9, bounds=(0.1, 10.0))  # quickflow storage coefficient == exp(-1/tau_q)
    # b::A = Param(0.1, bounds=(10.0, 1000.0))  # slowflow storage coefficent == exp(-1/tau_s)
    # storage_threshold::A = Param(10.0, ...)  # optional threshold controlling bounds between `a` and `b`
    storage_coef::P = Param(2.9, bounds=(1e-10, 10.0))
    alpha::P = Param(0.1, bounds=(1e-5, 1 - 1/10^9))

    level_params::Array{P, 1} = [
        Param(-0.01, bounds=(-10.0, -0.01)),  # p1
        Param(0.8, bounds=(0.0, 1.5)),  # p2
        Param(4.5, bounds=(0.0, 20.0)), # p3
        Param(5.0, bounds=(1.0, 10.0)), # p4
        Param(0.35, bounds=(0.0, 1.0)), # p5
        Param(1.41, bounds=(-2.0, 2.0)), # p6
        Param(-1.45, bounds=(-2.5, 0.0)), # p7
        Param(6.75, bounds=(0.0, 10.0)), # p8
        Param(150.0, bounds=(50.0, 200.0)) # ctf
    ]

    storage::Array{A} = [100.0]
    quick_store::Array{A} = [0.0]
    slow_store::Array{A} = [0.0]
    outflow::Array{A} = []
    effective_rainfall::Array{A} = []
    et::Array{A} = []
    inflow::Array{A} = []
    level::Array{A} = []
    gw_store::Array{A} = [0.0]
end

# function prep_state!(node::IHACRESNode, sim_length::Int64)
#     ini_storage = node.storage[1]
#     node.storage = fill(0.0, sim_length+1)
#     node.storage[1] = ini_storage

#     ini_qs = node.quick_store[1]
#     node.quick_store = fill(0.0, sim_length+1)
#     node.quick_store[1] = ini_qs

#     ini_ss = node.slow_store[1]
#     node.slow_store = fill(0.0, sim_length+1)
#     node.slow_store[1] = ini_ss

#     node.outflow = fill(0.0, sim_length)
#     node.effective_rainfall = fill(0.0, sim_length)
#     node.et = fill(0.0, sim_length)
#     node.inflow = fill(0.0, sim_length)
#     node.level = fill(0.0, sim_length)

#     ini_gw = node.gw_store[1]
#     node.gw_store = fill(0.0, sim_length+1)
#     node.gw_store[1] = ini_gw
# end

function BilinearNode(name::String, spec::Dict)
    n = create_node(BilinearNode, name, spec["area"])
    node_params = copy(spec["parameters"])

    if haskey(spec, "level_params")
        n_lparams = n.level_params
        s_lparams = spec["level_params"]
        node_params["level_params"] = Param[
            Param(s_lparams[1], bounds=n_lparams[1].bounds)
            Param(s_lparams[2], bounds=n_lparams[2].bounds)
            Param(s_lparams[3], bounds=n_lparams[3].bounds)
            Param(s_lparams[4], bounds=n_lparams[4].bounds)
            Param(s_lparams[5], bounds=n_lparams[5].bounds)
            Param(s_lparams[6], bounds=n_lparams[6].bounds)
            Param(s_lparams[7], bounds=n_lparams[7].bounds)
            Param(s_lparams[8], bounds=n_lparams[8].bounds)
            Param(s_lparams[9], bounds=n_lparams[9].bounds)
        ]
    end

    node_params["storage"] = [node_params["initial_storage"]]
    delete!(node_params, "initial_storage")

    for (k, p) in node_params
        s = Symbol(k)
        if p isa String
            p = eval(Meta.parse(p))
        end

        try
            f = getfield(n, s)
            setfield!(n, s, Param(p, bounds=f.bounds))
        catch err
            msg = sprint(showerror, err, catch_backtrace())
            if occursin("no field bounds", string(msg))
                setfield!(n, s, p)
            else
                throw(err)
            end
        end
    end

    return n
end

"""
    BilinearNode(name::String, area::Float64, d::Float64, d2::Float64, e::Float64, f::Float64,
                a::Float64, b::Float64, s_coef::Float64, alpha::Float64,
                store::Float64, quick::Float64, slow::Float64, gw_store::Float64)

Create a IHACRES node that adopts the bilinear CMD module.
"""
function BilinearNode(name::String, area::Float64, d::Float64, d2::Float64, e::Float64, f::Float64,
                      a::Float64, b::Float64, s_coef::Float64, alpha::Float64,
                      store::Float64, quick::Float64, slow::Float64, gw_store::Float64)
    n = create_node(BilinearNode, name, area)

    update_params!(n, d, d2, e, f, a, b, s_coef, alpha)

    n.storage = [store]
    n.quick_store = [quick]
    n.slow_store = [slow]
    n.gw_store = [gw_store]

    return n
end

function update_state(s_node::IHACRESNode, storage::Float64, e_rainfall::Float64, et::Float64,
                      qflow_store::Float64, sflow_store::Float64, outflow::Float64,
                      level::Float64, gw_store::Float64)
    push!(s_node.storage, storage)
    push!(s_node.effective_rainfall, e_rainfall)
    push!(s_node.et, et)
    push!(s_node.outflow, outflow)

    push!(s_node.quick_store, qflow_store)
    push!(s_node.slow_store, sflow_store)
    push!(s_node.level, level)
    push!(s_node.gw_store, gw_store)
end

"""
    run_node!(node::IHACRESNode, climate::Climate, timestep::Int;
                  inflow=nothing, extraction=nothing, exchange=nothing)

Run a specific node for a specified time step.

# Arguments
- `node::IHACRESNode` :
- `climate::Climate` :
- `timestep::Int` : current time step
- `inflow::DataFrame` : Time series of inflows from any upstream node.
- `extraction::DataFrame` : Time series of water orders (expects column of `_releases`)
- `exchange::DataFrame` : Time series of groundwater flux
"""
function run_node!(node::IHACRESNode, climate::Climate, timestep::Int;
                   inflow=nothing, extraction=nothing, exchange=nothing)
    ts = timestep
    if checkbounds(Bool, node.outflow, ts)
        if node.outflow[ts] != undef
            # already ran for this time step so no need to run
            return node.outflow[ts], node.level[ts]
        end
    end

    node_name = node.name
    rain, et = climate_values(node, climate, ts)
    wo = timestep_value(ts, node_name, "releases", extraction)
    ex = timestep_value(ts, node_name, "exchange", exchange)
    in_flow = timestep_value(ts, node_name, "inflow", inflow)

    return run_timestep!(node, rain, et, ts; inflow=in_flow, extraction=wo, exchange=ex)
end

"""
    run_step!(s_node::BilinearNode,
              rain::Float64,
              evap::Float64,
              inflow::Float64,
              ext::Float64,
              gw_exchange::Float64;
              current_store::Union{Nothing, Float64}=nothing,
              quick_store::Union{Nothing, Float64}=nothing,
              slow_store::Union{Nothing, Float64}=nothing,
              gw_store::Union{Nothing, Float64}=nothing)::Tuple{Float64, Float64}

Run node with ET data to calculate outflow and update state.

# Arguments
- `s_node::BilinearNode` : IHACRESNode
- rain : rainfall for time step
- evap : evapotranspiration for time step
- inflow : inflow from previous node
- ext : irrigation and other water extractions
- gw_exchange : flux in ML where positive is contribution to stream, negative is infiltration
- current_store : replacement cmd state value (uses last known state if not provided)
- quick_store : replacement quick store state value
- slow_store : replacement slow store state value
- gw_store : replacement groundwater store state value

# Returns
- float, outflow from node [ML/day], stream level
"""
function run_step!(s_node::BilinearNode,
                    rain::Float64,
                    evap::Float64,
                    inflow::Float64,
                    ext::Float64,
                    gw_exchange::Float64;
                    current_store::Union{Nothing, Float64}=nothing,
                    quick_store::Union{Nothing, Float64}=nothing,
                    slow_store::Union{Nothing, Float64}=nothing,
                    gw_store::Union{Nothing, Float64}=nothing)::Tuple{Float64, Float64}

    if !isnothing(current_store)
        s_node.storage[end] = current_store
    end

    if !isnothing(quick_store)
        s_node.quick_store[end] = quick_store
    end

    if !isnothing(slow_store)
        s_node.slow_store[end] = slow_store
    end

    if !isnothing(gw_store)
        s_node.gw_store[end] = gw_store
    end

    current_store = s_node.storage[end]
    quick_store = s_node.quick_store[end]
    slow_store = s_node.slow_store[end]
    gw_store = s_node.gw_store[end]

    interim_results = [0.0, 0.0, 0.0]
    @ccall IHACRES.calc_ft_interim_cmd(interim_results::Ptr{Cdouble},
                                       current_store::Cdouble,
                                       rain::Cdouble,
                                       s_node.d::Cdouble,
                                       s_node.d2::Cdouble,
                                       s_node.alpha::Cdouble)::Cvoid

    (mf, e_rainfall, recharge) = interim_results

    et::Float64 = @ccall IHACRES.calc_ET(
        s_node.e::Cdouble,
        evap::Cdouble,
        mf::Cdouble,
        s_node.f::Cdouble,
        s_node.d::Cdouble
    )::Cdouble

    cmd::Float64 = @ccall IHACRES.calc_cmd(
        current_store::Cdouble,
        rain::Cdouble,
        et::Cdouble,
        e_rainfall::Cdouble,
        recharge::Cdouble
    )::Cdouble

    flow_results = [0.0, 0.0, 0.0]
    @ccall IHACRES.calc_ft_flows(
        flow_results::Ptr{Cdouble},
        quick_store::Cdouble,
        slow_store::Cdouble,
        e_rainfall::Cdouble,
        recharge::Cdouble,
        s_node.area::Cdouble,
        s_node.a::Cdouble,
        s_node.b::Cdouble
    )::Cvoid

    (nq_store, ns_store, outflow) = flow_results

    routing_res = [0.0, 0.0]
    @ccall IHACRES.routing(
        routing_res::Ptr{Cdouble},
        gw_store::Cdouble,
        s_node.storage_coef::Cdouble,
        inflow::Cdouble,
        outflow::Cdouble,
        ext::Cdouble,
        gw_exchange::Cdouble)::Cvoid

    (gw_store, outflow) = routing_res

    level_params = Array{Float64}(s_node.level_params)
    level::Float64 = @ccall IHACRES.calc_ft_level(
        outflow::Cdouble,
        level_params::Ptr{Cdouble}
    )::Cdouble

    push!(s_node.inflow, inflow)
    update_state(s_node, cmd, e_rainfall, et, nq_store, ns_store, outflow, level, gw_store)

    return outflow, level
end

"""
    run_timestep!(s_node::BilinearNode,
                  rain::Float64,
                  evap::Float64,
                  timestep::Union{Int64, Nothing};
                  inflow::Float64,
                  extraction::Float64,
                  exchange::Float64)::Tuple{Float64, Float64}

Run node for a given time step.
"""
function run_timestep!(s_node::BilinearNode,
                       rain::Float64,
                       evap::Float64,
                       timestep::Union{Int64, Nothing};
                       inflow::Float64=0.0,
                       extraction::Float64=0.0,
                       exchange::Float64=0.0)::Tuple{Float64, Float64}
    ts = timestep
    if !isnothing(ts)
        current_store = s_node.storage[ts]
        quick_store = s_node.quick_store[ts]
        slow_store = s_node.slow_store[ts]
        gw_store = s_node.gw_store[ts]
    else
        current_store = nothing
        quick_store = nothing
        slow_store = nothing
        gw_store = nothing
    end

    return run_step!(s_node, rain, evap, inflow, extraction, exchange;
                     current_store, quick_store, slow_store, gw_store)
end

"""
    param_info(node::IHACRESNode; with_level::Bool = true)::Tuple

Extract node parameter names, values, and bounds for IHACRESNode types.
"""
function param_info(node::IHACRESNode; with_level::Bool = false)::Tuple
    tmp = Model(node)
    values = collect(tmp[:val])
    bounds = collect(tmp[:bounds])
    param_names = collect(tmp[:fieldname])

    if with_level
        level_param_vals = map(x->x.val, node.level_params)
        append!(values, level_param_vals)

        level_param_bounds = map(x->x.bounds, node.level_params)
        append!(bounds, level_param_bounds)
    end

    return param_names, values, bounds
end

function run_node_with_temp!(sn::StreamfallNetwork, nid::Int64, climate::Climate;
                             inflow=nothing, extraction=nothing, exchange=nothing)
    node = sn[nid]
    timesteps = sim_length(climate)
    for ts in 1:timesteps
        run_node_with_temp!(node, climate, ts;
                            inflow=inflow, extraction=extraction, exchange=exchange)
    end

    return node.outflow, node.level
end

function run_node_with_temp!(node::BilinearNode, climate::Climate;
                             inflow=nothing, extraction=nothing, exchange=nothing)
    timesteps = sim_length(climate)
    for ts in 1:timesteps
        run_node_with_temp!(node, climate, ts;
                            inflow=inflow, extraction=extraction, exchange=exchange)
    end

    return node.outflow, node.level
end

function run_node_with_temp!(node::BilinearNode, climate::Climate, timestep::Int64;
                             inflow=nothing, extraction=nothing, exchange=nothing)
    ts = timestep
    P, T = climate_values(node, climate, ts)
    i = timestep_value(ts, node.name, "_inflow", inflow)
    ext = timestep_value(ts, node.name, "_extraction", extraction)
    flux = timestep_value(ts, node.name, "_exchange", exchange)
    run_node_with_temp!(node, P, T, i, ext, flux)
end

"""
    run_node_with_temp!(s_node::BilinearNode,
                        rain::Float64,
                        temp::Float64,
                        inflow::Float64,
                        ext::Float64,
                        gw_exchange::Float64=0.0;
                        current_store=nothing,
                        quick_store=nothing,
                        slow_store=nothing
                        gw_store=nothing)::Tuple{Float64, Float64}

Run node with temperature data to calculate outflow and update state.
"""
function run_node_with_temp!(s_node::BilinearNode,
                             rain::Float64,
                             temp::Float64,
                             inflow::Float64,
                             ext::Float64,
                             gw_exchange::Float64=0.0;
                             current_store=nothing,
                             quick_store=nothing,
                             slow_store=nothing,
                             gw_store=nothing)::Tuple{Float64, Float64}
    if !isnothing(current_store)
        s_node.storage[end] = current_store
    end

    if !isnothing(quick_store)
        s_node.quick_store[end] = quick_store
    end

    if !isnothing(slow_store)
        s_node.slow_store[end] = slow_store
    end

    if !isnothing(gw_store)
        s_node.gw_store[end] = gw_store
    end

    current_store = s_node.storage[end]
    quick_store = s_node.quick_store[end]
    slow_store = s_node.slow_store[end]
    gw_store = s_node.gw_store[end]

    interim_results = [0.0, 0.0, 0.0]
    @ccall IHACRES.calc_ft_interim_cmd(interim_results::Ptr{Cdouble},
                                       current_store::Cdouble,
                                       rain::Cdouble,
                                       s_node.d::Cdouble,
                                       s_node.d2::Cdouble,
                                       s_node.alpha::Cdouble)::Cvoid

    @assert any(isnan.(interim_results)) == false
    (mf, e_rainfall, recharge) = interim_results

    et::Float64 = @ccall IHACRES.calc_ET_from_T(
        s_node.e::Cdouble,
        temp::Cdouble,
        mf::Cdouble,
        s_node.f::Cdouble,
        s_node.d::Cdouble
    )::Cdouble

    cmd::Float64 = @ccall IHACRES.calc_cmd(
        mf::Cdouble,
        rain::Cdouble,
        et::Cdouble,
        e_rainfall::Cdouble,
        recharge::Cdouble
    )::Cdouble

    flow_results = [0.0, 0.0, 0.0]
    @ccall IHACRES.calc_ft_flows(
        flow_results::Ptr{Cdouble},
        quick_store::Cdouble,
        slow_store::Cdouble,
        e_rainfall::Cdouble,
        recharge::Cdouble,
        s_node.area::Cdouble,
        s_node.a::Cdouble,
        s_node.b::Cdouble
    )::Cvoid

    @assert any(isnan.(flow_results)) == false
    (nq_store, ns_store, outflow) = flow_results

    routing_res = [0.0, 0.0]
    @ccall IHACRES.routing(
        routing_res::Ptr{Cdouble},
        gw_store::Cdouble,
        s_node.storage_coef::Cdouble,
        inflow::Cdouble,
        outflow::Cdouble,
        ext::Cdouble,
        gw_exchange::Cdouble)::Cvoid

    @assert any(isnan.(routing_res)) == false
    (gw_store, outflow) = routing_res

    level_params = Array{Float64}(s_node.level_params)
    level::Float64 = @ccall IHACRES.calc_ft_level(
        outflow::Cdouble,
        level_params::Ptr{Cdouble}
    )::Cdouble

    push!(s_node.inflow, inflow)
    update_state(s_node, cmd, e_rainfall, et, nq_store, ns_store, outflow, level, gw_store)

    return outflow, level
end

"""
    update_params!(node::BilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64,
                   a::Float64, b::Float64, s_coef::Float64, alpha::Float64)::Nothing

Update model parameters.
"""
function update_params!(node::BilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64,
                        a::Float64, b::Float64, s_coef::Float64, alpha::Float64)::Nothing
    node.d = Param(d, bounds=node.d.bounds::Tuple)
    node.d2 = Param(d2, bounds=node.d2.bounds::Tuple)
    node.e = Param(e, bounds=node.e.bounds::Tuple)
    node.f = Param(f, bounds=node.f.bounds::Tuple)
    node.a = Param(a, bounds=node.a.bounds::Tuple)
    node.b = Param(b, bounds=node.b.bounds::Tuple)
    node.storage_coef = Param(s_coef, bounds=node.storage_coef.bounds::Tuple)
    node.alpha = Param(alpha, bounds=node.alpha.bounds::Tuple)

    return nothing
end

"""
    update_params!(node::BilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64,
                   a::Float64, b::Float64, s_coef::Float64, alpha::Float64,
                   p1::Float64, p2::Float64, p3::Float64, p4::Float64, p5::Float64, p6::Float64, p7::Float64, p8::Float64, CTF::Float64)::Nothing

Update all parameters.
"""
function update_params!(node::BilinearNode, d::Float64, d2::Float64, e::Float64, f::Float64,
                        a::Float64, b::Float64, s_coef::Float64, alpha::Float64,
                        p1::Float64, p2::Float64, p3::Float64, p4::Float64, p5::Float64,
                        p6::Float64, p7::Float64, p8::Float64, CTF::Float64)::Nothing
    node.d = Param(d, bounds=node.d.bounds::Tuple)
    node.d2 = Param(d2, bounds=node.d2.bounds::Tuple)
    node.e = Param(e, bounds=node.e.bounds::Tuple)
    node.f = Param(f, bounds=node.f.bounds::Tuple)
    node.a = Param(a, bounds=node.a.bounds::Tuple)
    node.b = Param(b, bounds=node.b.bounds::Tuple)
    node.storage_coef = Param(s_coef, bounds=node.storage_coef.bounds::Tuple)
    node.alpha = Param(alpha, bounds=node.alpha.bounds::Tuple)

    n_lparams = node.level_params
    node.level_params = [
        Param(p1, bounds=n_lparams[1].bounds::Tuple)
        Param(p2, bounds=n_lparams[2].bounds::Tuple)
        Param(p3, bounds=n_lparams[3].bounds::Tuple)
        Param(p4, bounds=n_lparams[4].bounds::Tuple)
        Param(p5, bounds=n_lparams[5].bounds::Tuple)
        Param(p6, bounds=n_lparams[6].bounds::Tuple)
        Param(p7, bounds=n_lparams[7].bounds::Tuple)
        Param(p8, bounds=n_lparams[8].bounds::Tuple)
        Param(CTF, bounds=n_lparams[9].bounds::Tuple)
    ]

    return nothing
end

# TODO
# update_bounds!()

"""
    reset!(s_node::IHACRESNode)::Nothing

Reset node. Clears all states back to their initial values.
"""
function reset!(s_node::IHACRESNode)::Nothing
    s_node.storage = [s_node.storage[1]]
    s_node.quick_store = [s_node.quick_store[1]]
    s_node.slow_store = [s_node.slow_store[1]]
    s_node.gw_store = [s_node.gw_store[1]]
    s_node.outflow = []
    s_node.level = []
    s_node.effective_rainfall = []
    s_node.et = []
    s_node.inflow = []

    return nothing
end


text:
module HKQM

using Reexport

using FastGaussQuadrature
using LinearAlgebra: dot, ⋅
using LinearAlgebra: cross, ×
using Polynomials
using ProgressMeter
using SpecialFunctions
using StaticArrays
using TensorOperations
using Tullio
@reexport using PeriodicTable
@reexport using Unitful
@reexport using UnitfulAtomic

# Submodule
include("submodules/ReferenceStates.jl")
using .ReferenceStates

import LinearAlgebra.dot
import LinearAlgebra.cross
import LinearAlgebra.normalize!

export AbstractElementGrid,
       AbstractOperator,
       CubicElementGrid,
       CubicElement,
       coulomb_correction,
       coulomb_tensor,
       CoulombTransformation,
       CoulombTransformationCombination,
       CoulombTransformationLocal,
       CoulombTransformationLog,
       CoulombTransformationLogLocal,
       coulombtransformation,
       DerivativeTensor,
       Element1D,
       elementsize,
       gaussian_coulomb_integral,
       gaussiandensity_self_energy,
       gausspoints,
       grid1d,
       normalize!,
       self_energy,
       test_accuracy,
       test_accuracy_new,
       test_accuracy_ad,
       test_kinetic_energy,
       transformation_tensor,
       transformation_tensor_alt,
       transformation_harrison,
       transformation_tensor_alt,
       xgrid,
       ygrid,
       zgrid

# Methods
export bracket,
       coulomb_operator,
       cross, ×,
       density_operator,
       dot, ⋅,
       get_center,
       helmholtz_equation,
       helmholtz_equation!,
       ketbra, ⋆,
       magnetic_current,
       momentum_operator,
       #nuclear_potential,
       optimal_coulomb_tranformation,
       para_magnetic_current,
       poisson_equation,
       poisson_equation!,
       position_operator,
       test_nuclear_potential,
       vector_potential,
       ω_tensor

# Abstract Types
export AbstractQuantumState

# Concrete Types
export CubicElementArray,
       DerivativeOperator,
       ElementGridVector,
       ElementVector,
       GradientOperator,
       HamiltonOperator,
       HamiltonOperatorFreeParticle,
       HamiltonOperatorMagneticField,
       LaplaceOperator,
       NuclearPotentialTensor,
       NuclearPotentialTensorCombination,
       NuclearPotentialTensorGaussian,
       NuclearPotentialTensorLog,
       NuclearPotentialTensorLogLocal,
       OperatorSum,
       ProjectionOperator,
       QuantumState,
       ScalarOperator,
       SlaterDeterminant,
       VectorOperator

include("abstract_types.jl")
include("elements.jl")
include("states.jl")
#include("bubbles.jl")
include("tensors.jl")
include("operators.jl")
include("integrations.jl")
include("accuracytests.jl")
include("scf.jl")

end


text:
using NNlib:softmax
using StatsBase:sample, Weights

"""
    struct LinearPolicy <: AbstractPolicy
        features::Array{Float64, 3}
        weights::Vector{Float64}
    end

The probability of each action is calculate by `features` and `weights` and then normalized by softmax.
"""
struct LinearPolicy <: AbstractPolicy
    features::Array{Float64, 3}
    weights::Vector{Float64}
end

(π::LinearPolicy)(s, ::Val{:dist}) = @view(π.features[s, :, :]) * π.weights |> softmax
(π::LinearPolicy)(s) = sample(Weights(π(s, Val(:dist)), 1.0))
(π::LinearPolicy)(s, a) = π(s, Val(:dist))[a]

function update!(π::LinearPolicy, e)
    π.weights .+= e
end


text:
import Dates

"""
    InvalidDayLDOM()

When the day calculated is invalid, move to the last day of the month. Will
use the "LDOM" short hand.
"""
struct InvalidDayLDOM <: InvalidDayConvention end
adjust(::InvalidDayLDOM, day, month, year) = Dates.Date(year, month, Dates.daysinmonth(year, month))
Base.show(io::IO, ::InvalidDayLDOM) = print(io, "LDOM")

"""
    InvalidDayFDONM()

When the day calculated is invalid, move to the first day of the next month.
Uses the "FDONM" short hand.
"""
struct InvalidDayFDONM <: InvalidDayConvention end
adjust(::InvalidDayFDONM, day, month, year) = month == 12 ? Dates.Date(year+1, 1, 1) : Dates.Date(year, month+1, 1)
Base.show(io::IO, ::InvalidDayFDONM) = print(io, "FDONM")

"""
    InvalidDayNDONM()

When the day calculated is invalid, move to the nth day of the next month where
n is the number of days past the last day of the month. Uses the "NDONM" short
hand.
"""
struct InvalidDayNDONM <: InvalidDayConvention end
function adjust(::InvalidDayNDONM, day, month, year)
    dayδ = day - Dates.daysinmonth(year, month)
    return month == 12 ? Dates.Date(year+1, 1, dayδ) : Dates.Date(year, month+1, dayδ)
end
Base.show(io::IO, ::InvalidDayNDONM) = print(io, "NDONM")

const INVALID_DAY_MAPPINGS = Dict(
    "LDOM" => InvalidDayLDOM(),
    "FDONM" => InvalidDayFDONM(),
    "NDONM" => InvalidDayNDONM()
)


text:
module RandomizedLasso

using Krylov: CgSolver, cg!, issolved
using Random
using LinearAlgebra
using StaticArrays
using Printf

include("utils.jl")
include("nystrom_pcg.jl")
include("admm.jl")

end


text:
export array_qudiff, prepare_init_state, LDEMSAlgHHL, bval, aval
export QuEuler, QuLeapfrog, QuAB2, QuAB3, QuAB4
export QuLDEMSProblem

using DiffEqBase
"""
    Based on : arxiv.org/abs/1010.2745v2

    * array_qudiff(N_t,N,h,A) - generates matrix for k-step solver
    * prepare_init_state(b,x,h,N_t) - generates inital states
    * solve_qudiff - solver

    x' = Ax + b

    * A - input matrix.
    * b - input vector.
    * x - inital vector
    * N - dimension of b (as a power of 2).
    * h - step size.
    * tspan - time span.
"""

"""
    LDEMSAlgHHL
    * step - step for multistep method
    * α - coefficients for xₙ
    * β - coefficent for xₙ'
"""
abstract type QuODEAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
#abstract type QuODEProblem{uType,tType,isinplace} <: DiffEqBase.AbstractODEProblem{uType,tType,isinplace} end
abstract type LDEMSAlgHHL <: QuODEAlgorithm end

struct QuLDEMSProblem{F,C,U,T} #<: QuODEProblem{uType,tType,isinplace}
    A::F
    b::C
    u0::U
    tspan::NTuple{2,T}

    #function QuLDEMSProblem(A,b,u0,tspan)
    #  new{typeof(u0),typeof(tspan),false,typeof(A),typeof(b)}(A,b,u0,tspan)
    #end
end

"""
    Explicit Linear Multistep Methods
"""
struct QuEuler{T}<:LDEMSAlgHHL
    step::Int
    α::Vector{T}
    β::Vector{T}

    QuEuler(::Type{T} = Float64) where {T} = new{T}(1,[1.0,],[1.0,])
end

struct QuLeapfrog{T}<:LDEMSAlgHHL
    step::Int
    α::Vector{T}
    β::Vector{T}

    QuLeapfrog(::Type{T} = Float64) where {T} = new{T}(2,[0, 1.0],[2.0, 0])
end
struct QuAB2{T}<:LDEMSAlgHHL
    step::Int
    α::Vector{T}
    β::Vector{T}

    QuAB2(::Type{T} = Float64) where {T} = new{T}(2,[1.0, 0], [1.5, -0.5])
end
struct QuAB3{T}<:LDEMSAlgHHL
    step::Int
    α::Vector{T}
    β::Vector{T}
    QuAB3(::Type{T} = Float64) where {T} = new{T}(3,[1.0, 0, 0], [23/12, -16/12, 5/12])
end
struct QuAB4{T}<:LDEMSAlgHHL
    step::Int
    α::Vector{T}
    β::Vector{T}

    QuAB4(::Type{T} = Float64) where {T} = new{T}(4,[1.0, 0, 0, 0], [55/24, -59/24, 37/24, -9/24])
end

function bval(alg::LDEMSAlgHHL,t,h,g::Function)
    b = zero(g(1))
    for i in 1:(alg.step)
        b += alg.β[i]*g(t-(i-1)*h)
    end
    return b
end

function aval(alg::LDEMSAlgHHL,t,h,g::Function)
    sz, = size(g(1))
    A = Array{ComplexF64}(undef,sz,(alg.step + 1)*sz)
    i_mat = Matrix{Float64}(I, size(g(1)))
    A[1:sz,sz*(alg.step) + 1:sz*(alg.step + 1)] = i_mat
    for i in 1:alg.step
        A[1:sz,sz*(i - 1) + 1: sz*i] = -1*(alg.α[alg.step - i + 1]*i_mat + h*alg.β[alg.step - i + 1]*g(t - (alg.step - i)*h))
    end
    return A
end

function prepare_init_state(tspan::NTuple{2, Float64},x::Vector,h::Float64,g::Function,alg::LDEMSAlgHHL)
    N_t = round(Int, (tspan[2] - tspan[1])/h + 1) #number of time steps
    N = nextpow(2,2*N_t + 1) # To ensure we have a power of 2 dimension for matrix
    sz, = size(g(1))
    init_state = zeros(ComplexF64,2*(N)*sz)
    #inital value
    init_state[1:sz] = x
    for i in 2:N_t
        b = bval(alg,h*(i - 1) + tspan[1],h,g)
        init_state[Int(sz*(i - 1) + 1):Int(sz*(i))] = h*b
    end
    return init_state
end

function array_qudiff(tspan::NTuple{2, Float64},h::Float64,g::Function,alg::LDEMSAlgHHL)
    sz, = size(g(1))
    i_mat = Matrix{Float64}(I, size(g(1)))
    N_t = round(Int, (tspan[2] - tspan[1])/h + 1) #number of time steps
    N = nextpow(2,2*N_t + 1) # To ensure we have a power of 2 dimension for matrix
    A_ = zeros(ComplexF64, N*sz, N*sz)
    # Generates First two rows
    @inbounds A_[1:sz, 1:sz] = i_mat
    @inbounds A_[sz + 1:2*sz, 1:sz] = -1*(i_mat + h*g(tspan[1]))
    @inbounds A_[sz + 1:2*sz,sz+1:sz*2] = i_mat
    #Generates additional rows based on k - step
    for i in 3:alg.step
        @inbounds A_[sz*(i - 1) + 1:sz*i, sz*(i - 3) + 1:sz*i] = aval(QuAB2(),(i-2)*h + tspan[1],h,g)
    end
    for i in alg.step + 1:N_t
        @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - alg.step - 1) + 1:sz*i] = aval(alg,(i - 2)*h + tspan[1],h,g)
    end
    #Generates half mirroring matrix
    for i in N_t + 1:N
        @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - 2) + 1:sz*(i - 1)] = -1*i_mat
        @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - 1) + 1:sz*i] = i_mat
    end
    A_ = [zero(A_) A_;A_' zero(A_)]
    return A_
end

function DiffEqBase.solve(prob::QuLDEMSProblem{F,C,U,T}, alg::LDEMSAlgHHL, dt = (prob.tspan[2]-prob.tspan[1])/100, n_reg::Int = 12) where {F,C,U,T}
    A = prob.A
    b = prob.b
    tspan = prob.tspan
    x = prob.u0

    mat = array_qudiff(tspan, dt, A, alg)
    state = prepare_init_state(tspan, x, dt, b, alg)
    λ = maximum(eigvals(mat))
    C_value = minimum(eigvals(mat) .|> abs)*0.01;
    mat = 1/(λ*2)*mat
    state = state*1/(2*λ) |> normalize!
    res = hhlsolve(mat,state, n_reg, C_value)
    res = res/λ
    return res
end;


text:
#=
Continuation of data, to e.g. fill missing values
=#
using SignalDecomposition
export sinusoidal_continuation

_numbertype(T) = _numbertype(eltype(T))
_numbertype(::Type{Union{T, Missing}}) where T = T
_numbertype(::Type{T}) where {T <: Real} = T

# TODO: This function needs to become more general to allow for time sampling
# other than monthly (currently explicitly assumed)

# TODO: perhaps it is worth considering a version that makes true time in seconds
# instead of months, which in the future would allow support of diurnal cycles
# TODO: YES! Make true time be seconds!

"""
    sinusoidal_continuation(T::ClimArray, fs = [1, 2]; Tmin = -Inf, Tmax = Inf)
Fill-in the missing values of spatiotemporal field `T`, by fitting sinusoidals
to the non-missing values, and then using the fitted sinusoidals for the missing values.

Frequencies are given per year (frequency 2 means 1/2 of a year).

`Tmin, Tmax` limits are used to clamp the result into this range (no clamping in the
default case).
"""
function sinusoidal_continuation(T, frequencies = [1.0, 2.0]; Tmin = -Inf, Tmax = Inf)
    E = _numbertype(T)
    lpv = Sinusoidal(E.(frequencies ./ DAYS_IN_YEAR))
    fullT = copy(T)
    # TODO: this must be extended to a general "true time" function
    truetime = Float32.(cumsum(daysinmonth.(dims(T, Time))))
    for i in otheridxs(T, Time)
        x = T[i...]
        any(ismissing, x) || continue # this space needs no correction
        mi = findall(!ismissing, x)
        x0 = Array{Float32}(x[mi])
        t0 = truetime[mi]
        sea, res = decompose(t0, x0, lpv)
        c = cosines(truetime, lpv)
        imi = findall(ismissing, x)
        x[imi] .= c[imi]
        fullT[i...] .= x
    end
    if !isinf(Tmin) && !isinf(Tmin)
        return clamp.(fullT, Tmin, Tmax)
    else
        return fullT
    end
end

function cosines(t, lpv::Sinusoidal)
    s = fill(lpv.A[1], length(t))
    for i in 2:length(lpv.A)
        @. s += lpv.A[i] * cos(2π * lpv.fs[i-1] * t + lpv.φ[i])
    end
    return s
end