text:
## Simple graph manipulations (nodes/edges/attributes) ##

#=
##      Operations with nodes
=#

function _addnode!(g::Graph{N,E}) where {N,E}
    @_inline_meta
    @assert(nodecount(g) < typemax(N), "node ids exhausted type $N")
    push!(g.nodes, Node{N,E}())
    return nodecount(g)
end

"""
    addnode!(g::Graph[; attrs...])

Add one node to the graph. If `g` has node attributes, you can set
their values for the new node through keywords.

# Example

```julia-repl
julia> addnode!(g);
julia> g[:, :color] = "red";
julia> addnode!(g, color="blue");
julia> g[:, :color]
 ...
 "red"
 "blue"
```
"""
function addnode!(g::Graph; attrs...)
    n = _addnode!(g)
    for (k, v) = attrs
        @assert(haskey(g.nodeattrs, k), "node attribute $k not found")
        g.nodeattrs[k, n] = v
    end
    return n
end

"""
    addnodes!(g::Graph, count::Integer[; attrs...])

Add several nodes to the graph. If `g` has node attributes, you can set
their values for the new nodes through keywords.

# Example

```julia-repl
julia> addnodes!(g, 3, color=["red", "green", "blue"]);
julia> g[:, :color]
 ...
 "red"
 "green"
 "blue"
```
"""
function addnodes!(g::Graph, count::Integer)
    @assert(count >= 0, "can't add a negative number of nodes")
    for i = 1:count
        _addnode!(g)
    end
    return nodecount(g) - count + 1:nodecount(g)
end

"""
    remnode!(g::Graph, n::Integer)

Remove node `u` from the graph.
"""
function remnode!(g::Graph{N}, n::Integer) where {N}
    for p = fwd_ptrs(g.nodes[n])    # remove pointers to this node; works correctly with self-loops!
        delete_ptr!(rev_ptrs(g.nodes[p.node]), N(n), p.id)
    end
    for p = rev_ptrs(g.nodes[n])
        delete_ptr!(fwd_ptrs(g.nodes[p.node]), N(n), p.id)
    end
    g.edgecount -= ptr_length(g.nodes[n], Both)
    if n < nodecount(g)     # move the last node to n if not already there
        m = nodecount(g)
        for p = fwd_ptrs(g.nodes[m])    # reset ptrs from m to n
            swap_ptr!(rev_ptrs(g.nodes[p.node]), N(m), p.id, N(n))
        end
        for p = rev_ptrs(g.nodes[m])
            if p.node != n
                swap_ptr!(fwd_ptrs(g.nodes[p.node]), N(m), p.id, N(n))
            else
                swap_ptr!(fwd_ptrs(g.nodes[m]), N(m), p.id, N(n))
            end
        end
        g.nodes[n] = pop!(g.nodes);
    else
        pop!(g.nodes);
    end
end

remnodes!(g::Graph, it) = for i = sort(it, rev=true); remnode!(g, i) end

# FIXME: teach these operations work with attributes nicely

#=
##      Operations with edges
=#

"""
    hasedge(g::Graph, e::Edge)
    hasedge(g::Graph, n::Integer, m::Integer)

Check if graph `g` has an edge between nodes `n` and `m`.
"""
hasedge(g::Graph{N,E,D}, e::Edge) where {N,E,D<:Directed} =
    has_ptr(fwd_ptrs(g.nodes[e.source], D), e.target)

hasedge(g::Graph{N,E,D}, e::Edge) where {N,E,D<:Undirected} =
    has_ptr(fwd_ptrs(g.nodes[e.source], D), e.target) ||
    has_ptr(rev_ptrs(g.nodes[e.source], D), e.target)

hasedge(g::Graph{N,E,D}, n::Integer, m::Integer) where {N,E,D<:Directed} =
    has_ptr(fwd_ptrs(g.nodes[n], D), N(m))

hasedge(g::Graph{N,E,D}, n::Integer, m::Integer) where {N,E,D<:Undirected} =
    has_ptr(fwd_ptrs(g.nodes[n], D), N(m)) ||
    has_ptr(rev_ptrs(g.nodes[n], D), N(m))

"""
    new_ptr_pair!(g::Graph, n, m)

Create a pair of `NodePtr`s to be used for adding a new edge.
"""
function new_ptr_pair!(g::Graph{N,E}, n::N, m::N) where {N,E<:Integer}
    @_inline_meta
    @assert(g.edgemaxid < typemax(E), "edge ids exhausted type $E")
    g.edgemaxid += one(E)
    return NodePtr(m, g.edgemaxid), NodePtr(n, g.edgemaxid)
end

function new_ptr_pair!(g::Graph{N,E}, n::N, m::N) where {N,E<:Nothing}
    @_inline_meta
    return NodePtr(m, nothing), NodePtr(n, nothing)
end

function _addedge!(g::MultiGraph{N,E,D}, n::Integer, m::Integer) where {N,E,D}
    fwd, rev = new_ptr_pair!(g, N(n), N(m))
    add_ptr!(fwd_ptrs(g.nodes[n], D), fwd)
    add_ptr!(rev_ptrs(g.nodes[m], D), rev)
    g.edgecount += 1
    return g.edgecount
end

# TODO: ensure that simple graphs can't have selfloops + make undirected graph checks

function _addedge!(g::SimpleGraph{N,E,D}, n::Integer, m::Integer) where {N,E,D}
    fwd, rev = new_ptr_pair!(g, N(n), N(m))
    add_one_ptr!(fwd_ptrs(g.nodes[n], D), fwd) &&
    add_one_ptr!(rev_ptrs(g.nodes[m], D), rev) ||
    error("graph already contains edge $n → $m")
    g.edgecount += 1
    return g.edgecount
end

"""
    addedge!{N,E}(g::Graph{N,E}, n::Integer, m::Integer[; attrs...])

Add an edge between nodes `n` and `m` in graph `g`.
If `g` has edge attributes, you can set
their values for the new edge through keywords `attrs...`.

# Example

```julia-repl
julia> addedge!(g, 1, 2, color="red");
julia> g[:, :, :color]
 ...
 "red"
```
"""
function addedge!(g::Graph, n::Integer, m::Integer; attrs...)
    e = _addedge!(g, n, m)
    for (k, v) = attrs
        @assert(haskey(g.edgeattrs, k), "edge attribute $k not found")
        g.edgeattrs[k, e] = v
    end
    return e
end

addedge!(g::Graph, pair::Pair{T}) where {T<:Integer} = addedge!(g, pair...)

const DyadTuple = Tuple{T,T,Vararg{T}} where {T<:Integer}
const DyadTupleVector = AbstractVector{T} where {T<:DyadTuple}

"""
    addedges!(g::Graph, es)

Add multiple edges to the graph. There are 3 ways to define them:

* Integer tuple, e.g., `(1, 2)` or `(1, 2, 3, 1)`. In second case,
3 edges will be created: `1 → 2`, `2 → 3`, and `3 → 1`.

* Vector of tuples, e.g., `[(1, 2), (3, 4), (4, 3, 2)]`.

* Generator returning integer 2-tuples, e.g., `((i, i + 1) for i = 1:10)`.

# Example

```julia-repl
julia> g = Graph(4);
julia> addedges!(g, (1, 2, 3, 4, 1))
```
"""
function addedges!(g::Graph, es::DyadTuple)
    for i = 1:length(es) - 1
        addedge!(g, es[i], es[i + 1])
    end
end

function addedges!(g::Graph, es::DyadTupleVector)
    for e = es
        addedges!(g, e)
    end
end

addedges!(g::Graph, it) = for i = collect(it); addedge!(g, i...) end

"""
    remedge!(g::Graph, e::Edge)
    remedge!(g::Graph, u::Integer, v::Integer)

Remove edge `e` from graph `g`.
"""
function remedge!(g::Graph{N,E,D}, e::Edge) where {N,E,D<:Directed}
    delete_ptr!(fwd_ptrs(g.nodes[e.source], D), e.target, e.id) &&
    delete_ptr!(rev_ptrs(g.nodes[e.target], D), e.source, e.id) ||
    error("edge $e not found")
    g.edgecount -= 1;
end

function remedge!(g::Graph{N,E,D}, e::Edge) where {N,E,D<:Undirected}
    delete_ptr!(fwd_ptrs(g.nodes[e.source], D), e.target, e.id) &&
    delete_ptr!(rev_ptrs(g.nodes[e.target], D), e.source, e.id) ||
    delete_ptr!(fwd_ptrs(g.nodes[e.target], D), e.source, e.id) &&
    delete_ptr!(rev_ptrs(g.nodes[e.source], D), e.target, e.id) ||
    error("edge $e not found")
    g.edgecount -= 1;
end

function remedge!(g::Graph{N,E,D}, n::Integer, m::Integer) where {N,E,D<:Directed}
    delete_ptr!(fwd_ptrs(g.nodes[n], D), N(m)) &&
    delete_ptr!(rev_ptrs(g.nodes[m], D), N(n)) ||
    error("edge $(Edge(n, m, 0, true)) not found")
    g.edgecount -= 1;
end

function remedge!(g::Graph{N,E,D}, n::Integer, m::Integer) where {N,E,D<:Undirected}
    delete_ptr!(fwd_ptrs(g.nodes[n], D), N(m)) &&
    delete_ptr!(rev_ptrs(g.nodes[m], D), N(n)) ||
    delete_ptr!(fwd_ptrs(g.nodes[m], D), N(n)) &&
    delete_ptr!(rev_ptrs(g.nodes[n], D), N(m)) ||
    error("edge $(Edge(n, m, 0, false)) not found")
    g.edgecount -= 1;
end

remedges!(g::Graph, it) = for x = collect(it); remedge!(g, x); end

# FIXME: allow to use list of self-loops for bulk edge deletion

#=
##      Graph traversal
=#

"""
    neighbors(g::Graph, n::Integer)

Ids of all neighbors of node `n`.
"""
neighbors(g::Graph, n::Integer) = PtrView(g, n, Both, get_node_id)

"""
    outneighbors(g::Graph, n::Integer)

Ids of nodes towards which node `n` has an edge.
"""
outneighbors(g::Graph, n::Integer) = PtrView(g, n, Forward, get_node_id)

"""
    inneighbors(g::Graph, n::Integer)

Ids of nodes having an edge to node `n`.
"""
inneighbors(g::Graph, n::Integer) = PtrView(g, n, Reverse, get_node_id)

"""
    edgeids(g::Graph, [n::Integer, [m::Integer]])

Iterate over ids of (all, out-, or in-) edges starting at node `n` and
optionally ending in `m`.
"""
edgeids(g::Graph) = Base.flatten(PtrView(g, n, Forward, get_edge_id) for n = nodes(g))

edgeids(g::Graph, n::Integer) = PtrView(g, n, Both, get_edge_id)
outedgeids(g::Graph, n::Integer) = PtrView(g, n, Forward, get_edge_id)
inedgeids(g::Graph, n::Integer) = PtrView(g, n, Reverse, get_edge_id)

# FIXME: edgeids
# edgeids(g::Graph{N,E,D}, n::Integer, m::Integer) where {N,E,D} = EdgeRangeIDView(g.nodes[n], m, D)

"""
    edges(g::Graph, [n::Integer, [m::Integer]])

Iterate over all edges in the graph, or, if `n` is provided, edges that are
adjacent to node `n`. If both `n` and `m` are provided, list edges between those
two nodes.
"""
edges(g::Graph) = EdgeIter(g)
edges(g::Graph, n::Integer) = PtrView(g, n, Both, get_edge)

# FIXME: implement next 4 methods
# edges(g::Graph{N,E,D}, n::Integer, m::Integer) where {N,E,D} = EdgeRangeView(g.nodes[n], n, m, D, Forward)
# edges(g::Graph{N,E,D}, n::Integer, m::UnitRange) where {N,E,D} = EdgeRangeView(g.nodes[n], n, m, D, Forward)
# edges(g::Graph{N,E,D}, n::UnitRange, m::Integer) where {N,E,D} = EdgeRangeView(g.nodes[m], m, n, D, Reverse)
# edges(g::Graph{N,E,D}, n::UnitRange, m::UnitRange) where {N,E,D} = Base.Flatten(EdgeRangeView(g.nodes[i], i, m, Forward, D) for i = n)

"""
    outedges(g::Graph, n::Integer)

Iterate over all edges that start at node `n`.
"""
outedges(g::Graph, n::Integer) = PtrView(g, n, Forward, get_edge)
outedges(g::Graph, n::AbstractVector) = Base.flatten(PtrView(g, i, Forward, get_edge) for i = n)

"""
    inedges(g::Graph, n::Integer)

Iterate over all edges that finish at node `n`.
"""
inedges(g::Graph, n::Integer) = PtrView(g, n, Reverse, get_edge)
inedges(g::Graph, n::AbstractVector) = Base.Flatten(PtrView(g, i, Reverse, get_edge) for i = n)

#=
##      Indexing `Graph` objects
=#

# Nodes
getindex(g::Graph, ::Colon) = nodes(g)
getindex(g::Graph, n) = g.nodes[n]

# Edges
getindex(g::Graph, ::Colon, ::Colon) = edges(g)
getindex(g::Graph, n, ::Colon) = outedges(g, n)
getindex(g::Graph, ::Colon, n) = inedges(g, n)
# getindex(g::Graph, n, m) = edges(g, n, m)

# Attributes
getindex(d::AttributeDict, s::Symbol, i) = d[s][i]

function setindex!(d::AttributeDict, v, s::Symbol, i)
    if typeof(d[s]) <: ConstantAttribute
        d[s] = SparseAttribute(d[s])
    elseif typeof(d[s]) <: SparseAttribute && density(d[s]) > .3
        d[s] = DenseAttribute(d[s])
    end
    d[s][i] = v
end

addnodeattr!(g::Graph, s::Symbol, v) = (g.nodeattrs[s] = nodeattr(g, v); nothing)
addnodeattrs!(g::Graph; attrs...) = for (k, v) = attrs; addnodeattr!(g, k, v) end

addedgeattr!(g::Graph, s::Symbol, v) = (g.edgeattrs[s] = edgeattr(g, v); nothing)
addedgeattrs!(g::Graph; attrs...) = for (k, v) = attrs; addedgeattr!(g, k, v) end

hasnodeattr(g::Graph, s::Symbol) = haskey(g.nodeattrs, s)
hasedgeattr(g::Graph, s::Symbol) = haskey(g.edgeattrs, s)

# Node attributes
setindex!(g::Graph, v, ::Colon, s::Symbol) = addnodeattr!(g, s, v)
setindex!(g::Graph, v, i, s::Symbol) = g.nodeattrs[s, i] = v

getindex(g::Graph, ::Colon, s::Symbol) = g.nodeattrs[s]
getindex(g::Graph, i, s::Symbol) = g.nodeattrs[s, i]

# Edge atrributes
setindex!(g::Graph, v, ::Colon, ::Colon, s::Symbol) = addedgeattr!(g, s, v)
setindex!(g::Graph, v, i::Integer, ::Colon, s::Symbol) = g.edgeattrs[s, collect(outedgeids(g, i))] = v
setindex!(g::Graph, v, ::Colon, i::Integer, s::Symbol) = g.edgeattrs[s, collect(inedgeids(g, i))] = v
# setindex!(g::Graph, v, i::Integer, j::Integer, s::Symbol) = g.edgeattrs[s, edgeids(g, i, j)] = v
setindex!(g::Graph, v, e::Edge, s::Symbol) = g.edgeattrs[s, e.id] = v

getindex(g::Graph, ::Colon, ::Colon, s::Symbol) = g.edgeattrs[s]
getindex(g::Graph, i::Integer, ::Colon, s::Symbol) = g.edgeattrs[s, collect(outedgeids(g, i))]
getindex(g::Graph, ::Colon, i::Integer, s::Symbol) = g.edgeattrs[s, collect(inedgeids(g, i))]
# getindex(g::Graph, i::Integer, j::Integer, s::Symbol) = g.edgeattrs[s][edgeids(g, i, j)]
getindex(g::Graph, e::Edge, s::Symbol) = g.edgeattrs[s, e.id]

#=
##      Convenience syntax
=#

maxid(edges::Nothing) = 0
maxid(edges::DyadTuple) = maximum(edges)
maxid(edges::DyadTupleVector) = maximum(maximum(i) for i = edges)
maxid(edges) = error("invalid specification of edges")

"""
    Graph()
    Graph(kwargs...)
    Graph(n[; kwargs...])
    Graph(edges[, n; kwargs...])

Construct a `Graph` object with `n` nodes and specified `edges`.
The following keyword arguments control the data structure being used:

* `directed=true` — whether constructed graph is directed,
* `sumple=false` — whether graph is simple (multi-edges and self-loops are prohibited),
* `TNode=UInt32` — type used for node IDs, should be `Integer`,
* `TEdge=UInt32` — type used for edge IDs, should be `Integer` or `Nothing`.
"""
function Graph(; directed=true, simple=false, TNode=UInt32, TEdge=UInt32)
    D = directed ? Forward : Both
    M = simple ? Simple : Multi
    return Graph{TNode,TEdge,D,M}()
end

function Graph(n::Integer; kwargs...)
    g = Graph(; kwargs...)
    addnodes!(g, n)
    return g
end

function Graph(edges::Union{Tuple,Vector}, n::Integer=maxid(edges); kwargs...)
    g = Graph(; kwargs...)
    addnodes!(g, n)
    addedges!(g, edges)
    return g
end


text:
using Ejemplo
using Test

@testset "Ejemplo.jl" begin
    @test f(2, 1) == 7
    @test f(2, 3) == 13
    @test f(1, 3) == 11
end

@testset "Derivada de f" begin
    @test ∂ₓf(124.24, 245.245) == 2
end

ext: jl


text:
#this example is modified from https://numericalenvironmental.wordpress.com/2018/05/26/a-steady-state-variably-saturated-flow-model-in-vertical-cross-section-a-finite-difference-approach-using-julia/
import Ipopt
using JuMP
import MathOptInterface

include("utilities.jl")
include("inputdeck.jl")
include("ex.jl")

function req_jump(Ks, neighbors::Array, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
	neighborsdict = Dict{Int, Set{Int}}()
	neighbor2neighborindex = Dict{Pair{Int, Int}, Int}()
	for i = 1:size(coords, 2)
		neighborsdict[i] = Set{Int}()
	end
	for (i, (node1, node2)) in enumerate(neighbors)
		push!(neighborsdict[node1], node2)
		push!(neighborsdict[node2], node1)
		neighbor2neighborindex[node1=>node2] = i
		neighbor2neighborindex[node2=>node1] = i
	end
	req_jump(Ks, neighborsdict, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs, neighbor2neighborindex)
end

function req_jump(Ks, neighbors_dict::Dict, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs, n2ni)
	model = Model(with_optimizer(Ipopt.Optimizer, max_cpu_time=5 * 60.0))
	psi0 = -ones(size(coords, 2))
	psi0[dirichletnodes] = zeros(length(dirichletnodes))
	@variable(model, psi[j=1:size(coords, 2)], start=psi0[j])
	@constraint(model, psi[dirichletnodes] .== dirichletpsis[dirichletnodes])
	dirichletnodesset = Set(dirichletnodes)
	for i = filter(x->!(x in dirichletnodesset), 1:size(coords, 2))
		fixedneighbors = Int[]
		freeneighbors = Int[]
		for n in neighbors_dict[i]
			if n in dirichletnodesset
				push!(fixedneighbors, n)
			else
				push!(freeneighbors, n)
			end
		end
		j = freeneighbors[1]
		@NLconstraint(model, 0 == Qs[i]
					  + sum(2 / (1 / ifelse(psi[i] < 0, (1 - abs(alphas[n2ni[i => nbor]] * psi[i]) ^ (Ns[n2ni[i => nbor]] - 1) * (1 + abs(alphas[n2ni[i => nbor]] * psi[i]) ^ Ns[n2ni[i => nbor]]) ^ -((Ns[n2ni[i => nbor]] - 1) / Ns[n2ni[i => nbor]])) ^ 2 / (1 + abs(alphas[n2ni[i => nbor]] * psi[i]) ^ Ns[n2ni[i => nbor]]) ^ ((Ns[n2ni[i => nbor]] - 1) / (2 * Ns[n2ni[i => nbor]])), 1.0) + 1 / (ifelse(psi[nbor] < 0, (1 - abs(alphas[n2ni[i => nbor]] * psi[nbor]) ^ (Ns[n2ni[i => nbor]] - 1) * (1 + abs(alphas[n2ni[i => nbor]] * psi[nbor]) ^ Ns[n2ni[i => nbor]]) ^ -((Ns[n2ni[i => nbor]] - 1) / Ns[n2ni[i => nbor]])) ^ 2 / (1 + abs(alphas[n2ni[i => nbor]] * psi[nbor]) ^ Ns[n2ni[i => nbor]]) ^ ((Ns[n2ni[i => nbor]] - 1) / (2 * Ns[n2ni[i => nbor]])), 1.0))) * Ks[n2ni[i=>nbor]] * (psi[nbor] + coords[2, nbor] - psi[i] - coords[2, i]) * areasoverlengths[n2ni[i=>nbor]] for nbor in freeneighbors)
					  + sum(2 / (1 / ifelse(psi[i] < 0, (1 - abs(alphas[n2ni[i => nbor]] * psi[i]) ^ (Ns[n2ni[i => nbor]] - 1) * (1 + abs(alphas[n2ni[i => nbor]] * psi[i]) ^ Ns[n2ni[i => nbor]]) ^ -((Ns[n2ni[i => nbor]] - 1) / Ns[n2ni[i => nbor]])) ^ 2 / (1 + abs(alphas[n2ni[i => nbor]] * psi[i]) ^ Ns[n2ni[i => nbor]]) ^ ((Ns[n2ni[i => nbor]] - 1) / (2 * Ns[n2ni[i => nbor]])), 1.0) + 1 / (ifelse(dirichletpsis[nbor] < 0, (1 - abs(alphas[n2ni[i => nbor]] * dirichletpsis[nbor]) ^ (Ns[n2ni[i => nbor]] - 1) * (1 + abs(alphas[n2ni[i => nbor]] * dirichletpsis[nbor]) ^ Ns[n2ni[i => nbor]]) ^ -((Ns[n2ni[i => nbor]] - 1) / Ns[n2ni[i => nbor]])) ^ 2 / (1 + abs(alphas[n2ni[i => nbor]] * dirichletpsis[nbor]) ^ Ns[n2ni[i => nbor]]) ^ ((Ns[n2ni[i => nbor]] - 1) / (2 * Ns[n2ni[i => nbor]])), 1.0))) * Ks[n2ni[i=>nbor]] * (dirichletpsis[nbor] + coords[2, nbor] - psi[i] - coords[2, i]) * areasoverlengths[n2ni[i=>nbor]] for nbor in fixedneighbors))
	end
	@objective(model, Min, 0)
	JuMP.optimize!(model)
	return JuMP.value.(psi), model
end

psi_jump, model = req_jump(Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
res(psi) = req_residuals(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)
jac(psi) = req_psi(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs)

callback(psi_jump, res(psi_jump), jac(psi_jump), 0)


text:
"""
Created on Fri 15 May 2020

@author: Yoann Pradat

Performs NMF rank selection by the method of factor stability across sligth perturbations of the V matrix.
"""

using Distances
using Distributed
using Distributions
using Random
using PyCall

@everywhere include("./src/VNMF.jl")
@everywhere using .VNMF

include(pwd() * "/demo/src/simulation/simulate.jl")

#### folder for saving plots and other files
folder   = "./demo/plot/simulation/rank_stability_a"
mkpath(folder)

#### # 1. SIMULATE DATA
#### # ################################################################################################################

#### fix rng
rng = MersenneTwister(123)

#### simulate clusters of columns (individuals) and clusters of rows (variables)
W, H, lims_K, lims_N, active_clusters = simulate_WH_clustered(
    F              = 96,
    N              = 200,
    K              = 5,
    n_clusters_N   = 5,
    dist_W         = truncated(Normal(0, 3), 0, Inf),
    dist_H         = Poisson(1000),
    H_cluster_mode = :a,
    rng            = rng
)

#### get col cluster indices
col_indices = zeros(Int, size(H,2))
for i in 1:(length(lims_N)-1)
    range_indices = (lims_N[i]+1):lims_N[i+1]
    col_indices[range_indices] .= i
end

#### simulate noise
E = rand(rng, Poisson(2), size(W,1), size(H, 2))

#### V input matrix (observations)
V = round.(W * H) + E

py"""
import sys
if "." not in sys.path:
    sys.path.append(".")

from demo.src.util.plot_VWH import *

def make_plot_V(V, col_indices, filepath):

    row_groups_colors = {
        "C>A": "cornflowerblue",
        "C>G": "lightseagreen",
        "C>T": "coral",
        "T>A": "sandybrown",
        "T>C": "forestgreen",
        "T>G": "deeppink"
    }
    fig = plot_V(V, col_indices)
    fig.savefig(filepath, bbox_inches="tight")
"""

filepath = joinpath(folder, "plot_V_true.pdf")
py"make_plot_V"(V, col_indices, filepath)

#### remove cluster order by shuffling
rand_order_V = randperm(rng, size(V,2))
V_obs = V[:, rand_order_V]
col_indices_obs = col_indices[rand_order_V]

filepath = joinpath(folder, "plot_V_obs.pdf")
py"make_plot_V"(V_obs, col_indices_obs, filepath)

#### save V to test on SigProfiler
folder_results   = "./demo/results/simulation/rank_stability_a"
mkpath(folder_results)

using DelimitedFiles
writedlm(joinpath(folder_results, "V_a_simulated.txt"), V, '\t')

#### # 2. SELECT ALGORITHM
#### # ################################################################################################################

#### params for rank selection
rs_params = RSParams(
    K_min  = 2,
    K_max  = 8,
    n_iter = 10,
    pert   = :Multinomial,
    seed   = 123
)

#### nmf algorithm
nmf_global_params = NMFParams(
    init          = :random,
    dist          = Uniform(0, 1),
    max_iter      = 20_000,
    stopping_crit = :cost,
    stopping_tol  = 1e-6,
    verbose       = false,
)

nmf_local_params = MUParams(
    β       = 2,
    div     = :β,
    scale_W = false,
    alg     = :mu
)

nmf = NMF(
    solver        = nmf_MU,
    global_params = nmf_global_params,
    local_params  = nmf_local_params
)

#### run rank selection by stability procedure
@time rs_results = rank_by_stability(V_obs, rs_params, nmf)

##### REPRODUCE FOR RANK = 5

nmf_global_params = NMFParams(
    init          = :random,
    dist          = Uniform(0, 1),
    rank          = 5,
    max_iter      = 20_000,
    stopping_crit = :cost,
    stopping_tol  = 1e-6,
    verbose       = false,
)

nmf_local_params = MUParams(
    β       = 2,
    div     = :β,
    scale_W = true,
    alg     = :mu
)

nmf = NMF(
    solver        = nmf_MU,
    global_params = nmf_global_params,
    local_params  = nmf_local_params
)

#### build grid
list_V_pert = VNMF._get_list_perturbed_matrix(V, 50, :Multinomial, rng=rng)
list_nmf    = Vector{NMF}()
for (rank, seed) in zip(5:5, 5:5)
    #### using a different seed for each run of NMF
    nmf_rank = copy(nmf)
    nmf_rank.global_params.rank = rank
    nmf_rank.global_params.seed = seed

    push!(list_nmf, nmf_rank)
end

grid = vec(collect(Iterators.product(list_V_pert, list_nmf)))

#### run in parallel
grid_nmf_results = @showprogress pmap(grid) do g
    VNMF._rank_by_stability_one(g)
end

W_all = hcat([nmf_results.W for nmf_results in grid_nmf_results]...)
H_all = vcat([nmf_results.H for nmf_results in grid_nmf_results]...)

#### recover parameters of the algorithm
n_fac = nmf_global_params.rank
n_nmf = size(grid_nmf_results, 1)
n_tot = n_fac * n_nmf


text:
#=

  Interpolation tests

   Unit test to test interpolation software.  The origArray is
   expanded by interpFactor to interpolateArray.  The origArray is
   filled with 1.0(s), interpolation should create a larger
   array, also filled with 1.0(s) if successful.

=#

include("../src/Interpolate.jl")

using .Interpolate, Test

#  lower/upper bounds for 2D and 3D  arrays
idims = Array{Int64}(undef, 6)
odims = Array{Int64}(undef, 6)
ϵ  = 1e-5

# Helper function to declare the size of a particular dimension as opposed to subscripts
indexdim(a, b) = b - a + 1

#------------------------------------------------------------------
# Test the 2D interpolation.
#------------------------------------------------------------------
@testset "2D Interpolation" begin

for ii = 3: 7 # create these 2D arrays
  for interpfactor = 1: 3 # create these interpFactor(s)

    odims[1] = 1
    odims[2] = ii
    odims[3] = 1
    odims[4] = ii + 1 # non - sqare matrix

    # Allocate and populate the originalArray with 1.0
    origarray2d = ones(Float64, ii, ii +1)

    # Calculate dimensions for interpolated array
    interpolateCalculateSpace2D!(odims, interpfactor, idims)

    # Create the interpolateArray
    interparray2d = Array{Float64, 2}(undef, indexdim(idims[1],idims[2]), indexdim(idims[3],idims[4]))
    # populate the interparray with 2.0 (2.0 # = 1.0)
    interparray2d .= 2.0

    # Interpolate the interpolateArray.
    interpolateArray2D!(origarray2d, odims, interparray2d, idims, interpfactor)

    # The inperpolateArray should be all 1.0, + / - ϵ.
    # Print failures.
    for j = idims[3]: idims[4]
      for k = idims[1]: idims[2]

        @test abs(interparray2d[k, j] - 1.0) < ϵ
        if abs(interparray2d[k, j] - 1.0) > ϵ
          println("Failure 2D")
          println( k, " , ", j , " , ", interparray2d[k, j], " , ",  interparray2d[k, j] - 1.0)
          break
        end

      end
    end

  end
end

end # testset

#------------------------------------------------------------------
# Test the 3D interpolation.
#------------------------------------------------------------------
@testset "3D Interpolation " begin

for kk = 2: 3
  for ii = 3: 7 # create these 2D arrays
    for interpfactor = 1: 5 # create these interpFactor(s)

        # Allocate and populate the origArray3D with ones
        origarray3d = ones(Float64, ii, ii +1, kk)

        odims[1] = 1
        odims[2] = ii
        odims[3] = 1
        odims[4] = ii + 1 # non - sqare matrix
        odims[5] = 1
        odims[6] = kk

    # Calculate dimensions for interpolated array
    interpolateCalculateSpace3D!(odims, interpfactor, idims)
    # Create the interpolateArray
    interparray3d = Array{Float64, 3}(undef, indexdim(idims[1],idims[2]), indexdim(idims[3],idims[4]), indexdim(idims[5],idims[6]))
    # populate the interparray with 2.0 (2.0 # = 1.0)
    interparray3d .= 2.0

    # Interpolate the interpolateArray.
    interpolateArray3D!(origarray3d, odims, interparray3d, idims, interpfactor)

    # The inperpolateArray should be all 1.0, + / - ϵ.
    # Print failures.
    for i = idims[5]: idims[6]
      for j = idims[3]: idims[4]
        for k = idims[1]: idims[2]

          @test abs(interparray3d[k, j, i] - 1.0) < ϵ
            if abs(interparray3d[k, j, i] - 1.0) > ϵ
            println("failure 3d")
            println( "k, j, i = ", k, " , ",j, " , ", i, " , ", interparray3d[k, j, i], " , ", interparray3d[k, j, i] - 1.0)
            break
          end

        end
      end
    end

    end
  end
end

end #testset

return true


text:
module GraphBuilding

export build_graph
export weights_voronoi, weights_cutoff
export inverse_square, exp_decay

using PyCall
using ChemistryFeaturization
using Serialization
using Xtals
using NearestNeighbors
#rc[:paths][:crystals] = @__DIR__ # so that Xtals.jl knows where things are

# options for decay of bond weights with distance...
# user can of course write their own as well
inverse_square(x) = x^-2.0
exp_decay(x) = exp(-x)

"""
Build graph from a file storing a crystal structure (currently supports anything Xtals.jl can read in). Returns the weight matrix and elements used for constructing an `AtomGraph`.

# Arguments
## Required Arguments
- `file_path::String`: Path to ASE-readable file containing a molecule/crystal structure

## Keyword Arguments
- `use_voronoi::bool`: if true, use Voronoi method for neighbor lists, if false use cutoff method

    (The rest of these parameters are only used if `use_voronoi==false`)

- `cutoff_radius::Real=8.0`: cutoff radius for atoms to be considered neighbors (in angstroms)
- `max_num_nbr::Integer=12`: maximum number of neighbors to include (even if more fall within cutoff radius)
- `dist_decay_func::Function=inverse_square`: function to determine falloff of graph edge weights with neighbor distance
"""
function build_graph(
    file_path::String;
    use_voronoi::Bool = false,
    cutoff_radius::Real = 8.0,
    max_num_nbr::Integer = 12,
    dist_decay_func::Function = inverse_square,
)
    c = Crystal(abspath(file_path))
    atom_ids = String.(c.atoms.species)

    if use_voronoi
        @info "Note that building neighbor lists and edge weights via the Voronoi method requires the assumption of periodic boundaries. If you are building a graph for a molecule, you probably do not want this..."
        s = pyimport_conda("pymatgen.core.structure", "pymatgen", "conda-forge")
        struc = s.Structure.from_file(file_path)
        weight_mat = weights_voronoi(struc)
        return weight_mat, atom_ids, struc
    else
        build_graph(
            c;
            cutoff_radius = cutoff_radius,
            max_num_nbr = max_num_nbr,
            dist_decay_func = dist_decay_func,
        )
    end

end

"""
Build graph from a Crystal object. Currently only supports the "cutoff" method of neighbor list/weight calculation (not Voronoi).
This dispatch exists to support autodiff of graph-building.

# Arguments
## Required Arguments
- `crys::Crystal`: Crystal object representing the atomic geometry from which to build a graph

## Keyword Arguments
- `cutoff_radius::Real=8.0`: cutoff radius for atoms to be considered neighbors (in angstroms)
- `max_num_nbr::Integer=12`: maximum number of neighbors to include (even if more fall within cutoff radius)
- `dist_decay_func::Function=inverse_square`: function to determine falloff of graph edge weights with neighbor distance
"""
function build_graph(
    crys::Crystal;
    cutoff_radius::Real = 8.0,
    max_num_nbr::Integer = 12,
    dist_decay_func::Function = inverse_square,
)

    is, js, dists = neighbor_list(crys; cutoff_radius = cutoff_radius)
    weight_mat = weights_cutoff(
        is,
        js,
        dists;
        max_num_nbr = max_num_nbr,
        dist_decay_func = dist_decay_func,
    )
    return weight_mat, String.(crys.atoms.species), crys
end

"""
Build graph using neighbor number cutoff method adapted from original CGCNN.

!!! note
    `max_num_nbr` is a "soft" max, in that if there are more of the same distance as the last, all of those will be added.
"""
function weights_cutoff(is, js, dists; max_num_nbr = 12, dist_decay_func = inverse_square)
    # sort by distance
    ijd = sort([t for t in zip(is, js, dists)], by = t -> t[3])

    # initialize neighbor counts
    num_atoms = maximum(is)
    local nb_counts = Dict(i => 0 for i = 1:num_atoms)
    local longest_dists = Dict(i => 0.0 for i = 1:num_atoms)

    # iterate over list of tuples to build edge weights...
    # note that neighbor list double counts so we only have to increment one counter per pair
    weight_mat = zeros(Float32, num_atoms, num_atoms)
    for (i, j, d) in ijd
        # if we're under the max OR if it's at the same distance as the previous one
        if nb_counts[i] < max_num_nbr || isapprox(longest_dists[i], d)
            weight_mat[i, j] += dist_decay_func(d)
            longest_dists[i] = d
            nb_counts[i] += 1
        end
    end

    # average across diagonal, just in case
    weight_mat = 0.5 .* (weight_mat .+ weight_mat')

    # normalize weights
    weight_mat = weight_mat ./ maximum(weight_mat)
end

"""
Build graph using neighbors from faces of Voronoi polyedra and weights from areas. Based on the approach from https://github.com/ulissigroup/uncertainty_benchmarking/blob/aabb407807e35b5fd6ad06b14b440609ae09e6ef/BNN/data_pyro.py#L268
"""
function weights_voronoi(struc)
    num_atoms = size(struc)[1]
    sa = pyimport_conda("pymatgen.analysis.structure_analyzer", "pymatgen", "conda-forge")
    vc = sa.VoronoiConnectivity(struc)
    conn = vc.connectivity_array
    weight_mat = zeros(Float32, num_atoms, num_atoms)
    # loop over central atoms
    for atom_ind = 1:size(conn)[1]
        # loop over neighbor atoms
        for nb_ind = 1:size(conn)[2]
            # loop over each possible PBC image for chosen image
            for image_ind = 1:size(conn)[3]
                # only add as neighbor if atom is not current center one AND there is connectivity to image
                if (atom_ind != image_ind) && (conn[atom_ind, nb_ind, image_ind] != 0)
                    weight_mat[atom_ind, nb_ind] +=
                        conn[atom_ind, nb_ind, image_ind] / maximum(conn[atom_ind, :, :])
                end
            end
        end
    end

    # average across diagonal (because neighborness isn't strictly symmetric in the way we're defining it here)
    weight_mat = 0.5 .* (weight_mat .+ weight_mat')

    # normalize weights
    weight_mat = weight_mat ./ maximum(weight_mat)
end

"""
Find all lists of pairs of atoms in `crys` that are within a distance of `cutoff_radius` of each other, respecting periodic boundary conditions.

Returns as is, js, dists to be compatible with ASE's output format for the analogous function.
"""
function neighbor_list(crys::Crystal; cutoff_radius::Real = 8.0)
    n_atoms = crys.atoms.n

    # make 3 x 3 x 3 supercell and find indices of "middle" atoms
    # as well as index mapping from outer -> inner
    supercell = replicate(crys, (3, 3, 3))

    # check for size of cutoff radius relative to size of cell
    min_celldim = min(crys.box.a, crys.box.b, crys.box.c)
    if cutoff_radius >= min_celldim
        @warn "Your cutoff radius is quite large relative to the size of your unit cell. This may cause issues with neighbor list generation, and will definitely cause a very dense graph. To avoid issues, I'm setting it to be approximately equal to the smallest unit cell dimension."
        cutoff_radius = 0.99 * min_celldim
    end

    # todo: try BallTree, also perhaps other leafsize values
    #tree = BruteTree(sc.atoms.coords.xf, PeriodicEuclidean([1.0, 1.0, 1.0]))
    tree = BruteTree(Cart(supercell.atoms.coords, supercell.box).x)

    is_raw = 13*n_atoms+1:14*n_atoms
    js_raw =
        inrange(tree, Cart(supercell.atoms.coords[is_raw], supercell.box).x, cutoff_radius)

    index_map(i) = (i - 1) % n_atoms + 1 # I suddenly understand why some people dislike 1-based indexing

    # this looks horrifying but it does do the right thing...
    #ijraw_pairs = [p for p in Iterators.flatten([Iterators.product([p for p in zip(is_raw, js_raw)][n]...) for n in 1:4]) if p[1]!=p[2]]
    split1 = map(zip(is_raw, js_raw)) do x
        return [
            p for p in [(x[1], [j for j in js if j != x[1]]...) for js in x[2]] if
            length(p) == 2
        ]
    end
    ijraw_pairs = [(split1...)...]
    get_pairdist((i, j)) = distance(supercell.atoms, supercell.box, i, j, false)
    dists = get_pairdist.(ijraw_pairs)
    is = index_map.([t[1] for t in ijraw_pairs])
    js = index_map.([t[2] for t in ijraw_pairs])

    return is, js, dists
end

# TODO: graphs from SMILES via OpenSMILES.jl

end


text:
function create_cache_analysis(analyzer, mesh::TreeMesh{3},
                               equations, dg::DG, cache,
                               RealT, uEltype)

  # pre-allocate buffers
  # We use `StrideArray`s here since these buffers are used in performance-critical
  # places and the additional information passed to the compiler makes them faster
  # than native `Array`s.
  u_local = StrideArray(undef, uEltype,
                        StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
  u_tmp1  = StrideArray(undef, uEltype,
                        StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
  u_tmp2  = StrideArray(undef, uEltype,
                        StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
  x_local = StrideArray(undef, RealT,
                        StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
  x_tmp1  = StrideArray(undef, RealT,
                        StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
  x_tmp2  = StrideArray(undef, RealT,
                        StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))

  return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2)
end

function create_cache_analysis(analyzer, mesh::CurvedMesh{3},
                               equations, dg::DG, cache,
                               RealT, uEltype)

  # pre-allocate buffers
  # We use `StrideArray`s here since these buffers are used in performance-critical
  # places and the additional information passed to the compiler makes them faster
  # than native `Array`s.
  u_local = StrideArray(undef, uEltype,
                        StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
  u_tmp1  = StrideArray(undef, uEltype,
                        StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
  u_tmp2  = StrideArray(undef, uEltype,
                        StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
  x_local = StrideArray(undef, RealT,
                        StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
  x_tmp1  = StrideArray(undef, RealT,
                        StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
  x_tmp2  = StrideArray(undef, RealT,
                        StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
  jacobian_local = StrideArray(undef, RealT,
                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
  jacobian_tmp1  = StrideArray(undef, RealT,
                               StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
  jacobian_tmp2  = StrideArray(undef, RealT,
                               StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))

  return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2, jacobian_local, jacobian_tmp1, jacobian_tmp2)
end

function calc_error_norms(func, u, t, analyzer,
                          mesh::TreeMesh{3}, equations, initial_condition,
                          dg::DGSEM, cache, cache_analysis)
  @unpack vandermonde, weights = analyzer
  @unpack node_coordinates = cache.elements
  @unpack u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2 = cache_analysis

  # Set up data structures
  l2_error   = zero(func(get_node_vars(u, equations, dg, 1, 1, 1, 1), equations))
  linf_error = copy(l2_error)

  # Iterate over all elements for error calculations
  for element in eachelement(dg, cache)
    # Interpolate solution and node locations to analysis nodes
    multiply_dimensionwise!(u_local, vandermonde, view(u,                :, :, :, :, element), u_tmp1, u_tmp2)
    multiply_dimensionwise!(x_local, vandermonde, view(node_coordinates, :, :, :, :, element), x_tmp1, x_tmp2)

    # Calculate errors at each analysis node
    volume_jacobian_ = volume_jacobian(element, mesh, cache)

    for k in eachnode(analyzer), j in eachnode(analyzer), i in eachnode(analyzer)
      u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j, k), t, equations)
      diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j, k), equations)
      l2_error += diff.^2 * (weights[i] * weights[j] * weights[k] * volume_jacobian_)
      linf_error = @. max(linf_error, abs(diff))
    end
  end

  # For L2 error, divide by total volume
  total_volume_ = total_volume(mesh)
  l2_error = @. sqrt(l2_error / total_volume_)

  return l2_error, linf_error
end

function calc_error_norms(func, u, t, analyzer,
                          mesh::CurvedMesh{3}, equations, initial_condition,
                          dg::DGSEM, cache, cache_analysis)
  @unpack vandermonde, weights = analyzer
  @unpack node_coordinates, inverse_jacobian = cache.elements
  @unpack u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2, jacobian_local, jacobian_tmp1, jacobian_tmp2 = cache_analysis

  # Set up data structures
  l2_error   = zero(func(get_node_vars(u, equations, dg, 1, 1, 1, 1), equations))
  linf_error = copy(l2_error)
  total_volume = zero(real(mesh))

  # Iterate over all elements for error calculations
  for element in eachelement(dg, cache)
    # Interpolate solution and node locations to analysis nodes
    multiply_dimensionwise!(u_local, vandermonde, view(u,                :, :, :, :, element), u_tmp1, u_tmp2)
    multiply_dimensionwise!(x_local, vandermonde, view(node_coordinates, :, :, :, :, element), x_tmp1, x_tmp2)
    multiply_scalar_dimensionwise!(jacobian_local, vandermonde, inv.(view(inverse_jacobian, :, :, :, element)), jacobian_tmp1, jacobian_tmp2)

    # Calculate errors at each analysis node
    @. jacobian_local = abs(jacobian_local)

    for k in eachnode(analyzer), j in eachnode(analyzer), i in eachnode(analyzer)
      u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j, k), t, equations)
      diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j, k), equations)
      l2_error += diff.^2 * (weights[i] * weights[j] * weights[k] * jacobian_local[i, j, k])
      linf_error = @. max(linf_error, abs(diff))
      total_volume += weights[i] * weights[j] * weights[k] * jacobian_local[i, j, k]
    end
  end

  # For L2 error, divide by total volume
  l2_error = @. sqrt(l2_error / total_volume)

  return l2_error, linf_error
end

function integrate_via_indices(func::Func, u,
                               mesh::TreeMesh{3}, equations, dg::DGSEM, cache,
                               args...; normalize=true) where {Func}
  @unpack weights = dg.basis

  # Initialize integral with zeros of the right shape
  integral = zero(func(u, 1, 1, 1, 1, equations, dg, args...))

  # Use quadrature to numerically integrate over entire domain
  for element in eachelement(dg, cache)
    volume_jacobian_ = volume_jacobian(element, mesh, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      integral += volume_jacobian_ * weights[i] * weights[j] * weights[k] * func(u, i, j, k, element, equations, dg, args...)
    end
  end

  # Normalize with total volume
  if normalize
    integral = integral / total_volume(mesh)
  end

  return integral
end

function integrate_via_indices(func::Func, u,
                               mesh::CurvedMesh{3}, equations, dg::DGSEM, cache,
                               args...; normalize=true) where {Func}
  @unpack weights = dg.basis

  # Initialize integral with zeros of the right shape
  integral = zero(func(u, 1, 1, 1, 1, equations, dg, args...))
  total_volume = zero(real(mesh))

  # Use quadrature to numerically integrate over entire domain
  for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, k, element]))
      integral += volume_jacobian * weights[i] * weights[j] * weights[k] * func(u, i, j, k, element, equations, dg, args...)
      total_volume += volume_jacobian * weights[i] * weights[j] * weights[k]
    end
  end

  # Normalize with total volume
  if normalize
    integral = integral / total_volume
  end

  return integral
end

function integrate(func::Func, u,
                   mesh::Union{TreeMesh{3},CurvedMesh{3}},
                   equations, dg::DGSEM, cache; normalize=true) where {Func}
  integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, j, k, element, equations, dg
    u_local = get_node_vars(u, equations, dg, i, j, k, element)
    return func(u_local, equations)
  end
end

function analyze(::typeof(entropy_timederivative), du, u, t,
                 mesh::Union{TreeMesh{3},CurvedMesh{3}}, equations, dg::DG, cache)
  # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ
  integrate_via_indices(u, mesh, equations, dg, cache, du) do u, i, j, k, element, equations, dg, du
    u_node  = get_node_vars(u,  equations, dg, i, j, k, element)
    du_node = get_node_vars(du, equations, dg, i, j, k, element)
    dot(cons2entropy(u_node, equations), du_node)
  end
end

function analyze(::Val{:l2_divb}, du, u, t,
                 mesh::TreeMesh{3}, equations::IdealGlmMhdEquations3D,
                 dg::DG, cache)
  integrate_via_indices(u, mesh, equations, dg, cache, cache, dg.basis.derivative_matrix) do u, i, j, k, element, equations, dg, cache, derivative_matrix
    divb = zero(eltype(u))
    for l in eachnode(dg)
      divb += ( derivative_matrix[i, l] * u[6, l, j, k, element] +
                derivative_matrix[j, l] * u[7, i, l, k, element] +
                derivative_matrix[k, l] * u[7, i, j, l, element] )
    end
    divb *= cache.elements.inverse_jacobian[element]
    divb^2
  end |> sqrt
end

function analyze(::Val{:linf_divb}, du, u, t,
                 mesh::TreeMesh{3}, equations::IdealGlmMhdEquations3D,
                 dg::DG, cache)
  @unpack derivative_matrix, weights = dg.basis

  # integrate over all elements to get the divergence-free condition errors
  linf_divb = zero(eltype(u))
  for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      divb = zero(eltype(u))
      for l in eachnode(dg)
        divb += ( derivative_matrix[i, l] * u[6, l, j, k, element] +
                  derivative_matrix[j, l] * u[7, i, l, k, element] +
                  derivative_matrix[k, l] * u[7, i, j, l, element] )
      end
      divb *= cache.elements.inverse_jacobian[element]
      linf_divb = max(linf_divb, abs(divb))
    end
  end

  return linf_divb
end


text:
using PackageCompiler
create_app(".", "deps/VarianceComponentsHDFELinux", force=true, precompile_statements_file="deps/precompilation_statements.jl")


text:
@testset "getCurrentBlockRange" begin
    doc = settestdoc("ab")

    res = (LanguageServer.Position(0, 0), LanguageServer.Position(0, 2), LanguageServer.Position(0, 2))

    @test LanguageServer.julia_getCurrentBlockRange_request(LanguageServer.VersionedTextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), 0, LanguageServer.Position(0, 0)), server, server.jr_endpoint) == res

    @test LanguageServer.julia_getCurrentBlockRange_request(LanguageServer.VersionedTextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), 0, LanguageServer.Position(0, 1)), server, server.jr_endpoint) == res

    @test LanguageServer.julia_getCurrentBlockRange_request(LanguageServer.VersionedTextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("testdoc"), 0, LanguageServer.Position(0, 2)), server, server.jr_endpoint) == res
end


text:
"""
    TriggerJitter(trigger, offset)

or

    TriggerOffset(offset)

A trigger operation that apply [jitter](https://en.wikipedia.org/wiki/Jitter) to instant when a job should be triggered.

Addition `+` and substraction `-` are implemented so it's possible to define a new trigger using

    Trigger("H") + TriggerJitter(Date.Minute(3))

to be able to run a job every hour with a random jitter of 3 minutes.
This is same as:

    TriggerJitter(Trigger("H"), Date.Minute(3))

Randomize `next_dt_fire` by adding or subtracting a random value (the jitter). If the
resulting DateTime is in the past, returns the initial `next_dt_fire` without jitter.

next_dt_fire - jitter <= result <= next_dt_fire + jitter
"""
struct TriggerJitter
    trigger::AbstractTrigger
    jitter

    TriggerJitter() = error("ToDo")
end


text:
function check_num_nodes(g::GNNGraph, x::AbstractArray)
    @assert g.num_nodes == size(x, ndims(x))
end

sort_edge_index(eindex::Tuple) = sort_edge_index(eindex...)

function sort_edge_index(u, v)
    uv = collect(zip(u, v))
    p = sortperm(uv) # isless lexicographically defined for tuples
    return u[p], v[p]
end

cat_features(x1::Nothing, x2::Nothing) = nothing
cat_features(x1::AbstractArray, x2::AbstractArray) = cat(x1, x2, dims=ndims(x1))
cat_features(x1::Union{Number, AbstractVector}, x2::Union{Number, AbstractVector}) =
    cat(x1, x2, dims=1)

function cat_features(x1::NamedTuple, x2::NamedTuple)
    sort(collect(keys(x1))) == sort(collect(keys(x2))) ||
        @error "cannot concatenate feature data with different keys"

    NamedTuple(k => cat_features(getfield(x1,k), getfield(x2,k)) for k in keys(x1))
end

# Turns generic type into named tuple
normalize_graphdata(data::Nothing; kws...) = NamedTuple()

normalize_graphdata(data; default_name::Symbol, kws...) =
    normalize_graphdata(NamedTuple{(default_name,)}((data,)); default_name, kws...)

function normalize_graphdata(data::NamedTuple; default_name, n, duplicate_if_needed=false)
    # This had to workaround two Zygote bugs with NamedTuples
    # https://github.com/FluxML/Zygote.jl/issues/1071
    # https://github.com/FluxML/Zygote.jl/issues/1072

    if n == 1
        # If last array dimension is not 1, add a new dimension.
        # This is mostly usefule to reshape globale feature vectors
        # of size D to Dx1 matrices.
        function unsqz(v)
            if v isa AbstractArray && size(v)[end] != 1
                v = reshape(v, size(v)..., 1)
            end
            v
        end

        data = NamedTuple{keys(data)}(unsqz.(values(data)))
    end

    sz = map(x -> x isa AbstractArray ? size(x)[end] : 0, data)

    if duplicate_if_needed
        # Used to copy edge features on reverse edges
        @assert all(s -> s == 0 ||  s == n || s == n÷2, sz)

        function duplicate(v)
            if v isa AbstractArray && size(v)[end] == n÷2
                v = cat(v, v, dims=ndims(v))
            end
            v
        end

        data = NamedTuple{keys(data)}(duplicate.(values(data)))
    else
        @assert all(s -> s == 0 ||  s == n, sz)
    end
    return data
end

ofeltype(x, y) = convert(float(eltype(x)), y)


text:
function rel_energy_err_sde(sol)

    NQ = ndims(sol.q)

    if NQ==2

        en_ref  = 0.5 * ( sol.q[1,0]^2 + sol.q[2,0]^2 )
        en_last = 0.5 * ( sol.q[1,end]^2 + sol.q[2,end]^2 )

        return abs((en_last-en_ref)/en_ref)

    elseif NQ==3

        en_ref  = 0.5 * ( sol.q[1,0,:].^2 + sol.q[2,0,:].^2 )
        en_last = 0.5 * ( sol.q[1,end,:].^2 + sol.q[2,end,:].^2 )

        return maximum(abs.( (en_last .- en_ref) ./ en_ref ))

    end
end

function rel_energy_err_psde(sol)

    NQ = ndims(sol.q)

    if NQ==2

        en_ref  = 0.5 * ( sol.q[1,0]^2 + sol.p[1,0]^2 )
        en_last = 0.5 * ( sol.q[1,end]^2 + sol.p[1,end]^2 )

        return abs((en_last-en_ref)/en_ref)

    elseif NQ==3

        en_ref  = 0.5 * ( sol.q[1,0,:].^2 + sol.p[1,0,:].^2 )
        en_last = 0.5 * ( sol.q[1,end,:].^2 + sol.p[1,end,:].^2 )

        return maximum(abs.( (en_last .- en_ref) ./ en_ref ))

    end
end


text:
"""
**Paired Differences Index**

    pdi{T<:Number}(x::Array{T, 1})

This function will range the values of each row, so that the strongest link has
a value of one. This works for deterministic and quantitative networks.

> Poisot, T., Canard, E., Mouquet, N., Hochberg, M.E., 2012. A comparative study
> of ecological specialization estimators. Methods in Ecology and Evolution 3,
> 537–544. doi:10.1111/j.2041-210X.2011.00174.x

"""
function pdi(x::Vector{T}) where {T<:Number}
    s = reverse(sort(x))
    s = s ./ maximum(s)
    p = s[1] .- s
    return sum(p[2:end])/(length(s)-1)
end

"""
**Resource range**

    specificity(N::DeterministicNetwork)

Measure of specificity in a deterministic network. This returns a value between
0 and 1, where 1 indicates maximal specificity.

> Poisot, T., Canard, E., Mouquet, N., Hochberg, M.E., 2012. A comparative study
> of ecological specialization estimators. Methods in Ecology and Evolution 3,
> 537–544. doi:10.1111/j.2041-210X.2011.00174.x

"""
function specificity(N::BinaryNetwork)
    A = map(Int64, N.A)
    return Dict(zip(species(N; dims=1), vec(mapslices(pdi, A; dims=2))))
end

"""
**Paired Differences Index**

    specificity(N::QuantitativeNetwork)

Measure of specificity in a quantitative network. This returns a value between 0
and 1, where 1 indicates maximal specificity. Note that the PDI is measured
species-wise, and the maximal interaction strength of every species is set to 1.

> Poisot, T., Canard, E., Mouquet, N., Hochberg, M.E., 2012. A comparative study
> of ecological specialization estimators. Methods in Ecology and Evolution 3,
> 537–544. doi:10.1111/j.2041-210X.2011.00174.x

"""
function specificity(N::QuantitativeNetwork)
    return Dict(zip(species(N; dims=1), vec(mapslices(pdi, N.A; dims=2))))
end


text:
import 'lib_src/print.jl'
import 'lib_src/array.jl'
def main {
    a = []

    a:10 = 'doge'
    a:10 <- 'kitty'

    println(a)
}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
include("post_process_defs.jl")

### Parse Args ###

# Get
if length(ARGS) >= 1
  AWS_BUCKET = ARGS[1]
else
  println("usage: julia pullAndGenResults.jl <AWS_BUCKET>")
  exit(1)
end

# Create dir to do post processing
TMP_DIR = "TMP_DIR/"
mkpath(TMP_DIR)

RESULTS_DIR = split(read(`aws s3 ls $(AWS_BUCKET)/`, String))
RESULTS_DIR = filter(d -> d != "PRE", RESULTS_DIR)
RESULTS_DIR = filter(d -> !occursin("EXCLUDE_ME", d), RESULTS_DIR)
# RESULTS_DIR = filter(d -> occursin("K21", d) || occursin("missmech", d), RESULTS_DIR)

successes = []
failures = []

@time for dir in sort(RESULTS_DIR, rev=true)
  println("Processing: $dir")
  aws_dir = "$(AWS_BUCKET)/$(dir)"
  tmp_dir = "$(TMP_DIR)/$(dir)"

  try
    # Pull the dataset to TMP_DIR
    println("Pull data set from AWS to $(tmp_dir)...")
    cmd = `aws s3 sync $aws_dir $tmp_dir`
    println(cmd)
    run(cmd)

    # Post process
    println("Post processing...")
    path_to_mcmc = "$(tmp_dir)/output.jld2"
    @time post_process(path_to_mcmc)

    # Send back to AWS
    println("Sending results to S3...")
    run(`aws s3 sync $tmp_dir $aws_dir`)

    # Remove all files
    println("Remove results locally...")
    run(`rm -rf $(tmp_dir)`)

    append!(successes, [tmp_dir])
  catch
    append!(failures, [tmp_dir])
  end
end

println("Successes")
[println(s) for s in successes];
println()

println("Failures")
[println(f) for f in failures];
println()


text:
using Parameters
using LinearAlgebra
using StaticArrays
using Dierckx
using Calculus

abstract type AbstractDomain end
export AbstractDomain

abstract type AbstractConstantKDomain <: AbstractDomain  end
export AbstractConstantKDomain

abstract type AbstractVariableKDomain <: AbstractDomain end
export AbstractVariableKDomain

@with_kw struct ConstantTPDomain{N<:AbstractPhase,S<:Integer,W<:Real, W2<:Real, I<:Integer, Q<:AbstractArray} <: AbstractConstantKDomain
    phase::N
    interfaces::Array{AbstractInterface,1} = Array{AbstractInterface,1}()
    indexes::Q #assumed to be in ascending order
    constantspeciesinds::Array{S,1}
    T::W
    P::W
    kfs::Array{W,1}
    krevs::Array{W,1}
    efficiencyinds::Array{I,1}
    Gs::Array{W,1}
    rxnarray::Array{UInt16,2}
    mu::W = 0.0
    diffusivity::Array{W,1} = Array{Float64,1}()
    jacobian::Array{W,2} = Array{Float64,2}(undef,(0,0))
    sensitivity::Bool = false
    jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
    t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
end
function ConstantTPDomain(;phase::E2,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,X2},constantspecies::Array{X3,1}=Array{String,1}(),
    sparse::Bool=false,sensitivity::Bool=false) where {E<:Real,E2<:AbstractPhase,Q<:AbstractInterface,W<:Real,X,X2,X3}
    #set conditions and initialconditions
    T = 0.0
    P = 0.0
    if sensitivity
        y0 = zeros(length(phase.species)*(1+length(phase.species)+length(phase.reactions)))
    else
        y0 = zeros(length(phase.species))
    end
    spnames = [x.name for x in phase.species]
    for (key,val) in initialconds
        if key == "T"
            T = val
        elseif key == "P"
            P = val
        else
            ind = findfirst(isequal(key),spnames)
            @assert typeof(ind)<: Integer  "$key not found in species list: $spnames"
            y0[ind] = val
        end
    end

    @assert T != 0.0
    @assert P != 0.0
    ns = y0
    N = sum(ns)

    if length(constantspecies) > 0
        spcnames = getfield.(phase.species,:name)
        constspcinds = [findfirst(isequal(k),spcnames) for k in constantspecies]
    else
        constspcinds = Array{Int64,1}()
    end
    efficiencyinds = [rxn.index for rxn in phase.reactions if typeof(rxn.kinetics)<:AbstractFalloffRate && length(rxn.kinetics.efficiencies) > 0]
    Gs = calcgibbs(phase,T)
    if :solvent in fieldnames(typeof(phase)) && typeof(phase.solvent) != EmptySolvent
        mu = phase.solvent.mu(T)
    else
        mu = 0.0
    end
    if phase.diffusionlimited
        diffs = getfield.(phase.species,:diffusion)(T=T,mu=mu,P=P)
    else
        diffs = Array{typeof(T),1}()
    end
    C = P/(R*T)
    V = N*R*T/P
    kfs,krevs = getkfkrevs(phase=phase,T=T,P=P,C=C,N=N,ns=ns,Gs=Gs,diffs=diffs,V=V)
    if sparse
        jacobian=spzeros(typeof(T),length(phase.species),length(phase.species))
    else
        jacobian=zeros(typeof(T),length(phase.species),length(phase.species))
    end
    rxnarray = getreactionindices(phase)
    return ConstantTPDomain(phase,interfaces,SVector(phase.species[1].index,phase.species[end].index),constspcinds,
        T,P,kfs,krevs,efficiencyinds,Gs,rxnarray,mu,diffs,jacobian,sensitivity,MVector(false),MVector(0.0)), y0
end
export ConstantTPDomain

@with_kw struct ConstantVDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,Q<:AbstractArray} <: AbstractVariableKDomain
    phase::N
    interfaces::Array{AbstractInterface,1} = Array{AbstractInterface,1}()
    indexes::Q #assumed to be in ascending order
    constantspeciesinds::Array{S,1}
    V::W
    rxnarray::Array{UInt16,2}
    jacobian::Array{W,2}
    sensitivity::Bool = false
    jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
    t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
end
function ConstantVDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,E},constantspecies::Array{X2,1}=Array{String,1}(),
    sparse::Bool=false,sensitivity::Bool=false) where {E,X,X2,Z<:IdealGas,Q<:AbstractInterface}

    #set conditions and initialconditions
    T = 0.0
    P = 0.0
    V = 0.0
    ns = zeros(length(phase.species))
    spnames = [x.name for x in phase.species]
    for (key,val) in initialconds
        if key == "T"
            T = val
        elseif key == "P"
            P = val
        elseif key == "V"
            V = val
        else
            ind = findfirst(isequal(key),spnames)
            @assert typeof(ind)<: Integer  "$key not found in species list: $spnames"
            ns[ind] = val
        end
    end
    @assert V != 0.0 || (T != 0.0 && P != 0.0)
    N = sum(ns)
    if V == 0.0
        V = N*R*T/P
    elseif T == 0.0
        T = P*V/(R*N)
    elseif P == 0.0
        P = N*R*T/V
    else
        throw(error("ConstantVDomain overspecified with T,P and V"))
    end
    if sensitivity
        y0 = vcat(ns,T,zeros((length(ns)+1)*(length(ns)+length(phase.reactions))))
    else
        y0 = vcat(ns,T)
    end
    if length(constantspecies) > 0
        spcnames = getfield.(phase.species,:name)
        constspcinds = [findfirst(isequal(k),spcnames) for k in constantspecies]
    else
        constspcinds = Array{Int64,1}()
    end
    if sparse
        jacobian=zeros(typeof(T),length(phase.species)+1,length(phase.species)+1)
    else
        jacobian=zeros(typeof(T),length(phase.species)+1,length(phase.species)+1)
    end
    rxnarray = getreactionindices(phase)
    return ConstantVDomain(phase,interfaces,SVector(phase.species[1].index,phase.species[end].index,phase.species[end].index+1),constspcinds,
    V,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0)), y0
end
export ConstantVDomain

@with_kw struct ParametrizedTPDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,Q<:AbstractArray} <: AbstractVariableKDomain
    phase::N
    interfaces::Array{AbstractInterface,1} = Array{AbstractInterface,1}()
    indexes::Q #assumed to be in ascending order
    constantspeciesinds::Array{S,1}
    T::Function
    P::Function
    rxnarray::Array{UInt16,2}
    jacobian::Array{W,2}
    sensitivity::Bool = false
    jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
    t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
end
function ParametrizedTPDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,Any},constantspecies::Array{X2,1}=Array{String,1}(),
    sparse::Bool=false,sensitivity::Bool=false) where {X,X2,Z<:IdealGas,Q<:AbstractInterface}

    #set conditions and initialconditions
    T = 0.0
    P = 0.0
    V = 0.0
    ts = 0.0
    ns = zeros(length(phase.species))
    spnames = [x.name for x in phase.species]
    for (key,val) in initialconds
        if key == "T"
            T = val
        elseif key == "P"
            P = val
        elseif key == "V"
            V = val
        elseif key == "ts"
            ts = val
        else
            ind = findfirst(isequal(key),spnames)
            @assert typeof(ind) <: Integer  "$key not found in species list: $spnames"
            ns[ind] = val
        end
    end
    @assert V != 0.0 || (T != 0.0 && P != 0.0)
    if isa(T,AbstractArray)
        q = Spline1D(ts,T;k=3,s=1e-11)
        Tfcn(x::Float64) = q(x)
    elseif isa(T,Function)
        Tfcn = T
    else
        throw(error("ParametrizedTPDomain must take \"T\" as a function or if an array of times for \"ts\" is supplied as an array of volumes"))
    end
    if isa(P,AbstractArray)
        v = Spline1D(ts,P;k=3,s=1e-11)
        Pfcn(x::Float64) = v(x)
    elseif isa(P,Function)
        Pfcn = P
    else
        throw(error("ParametrizedTPDomain must take \"P\" as a function or if an array of times for \"ts\" is supplied as an array of volumes"))
    end

    N = sum(ns)
    V = N*R*Tfcn(0.0)/Pfcn(0.0)

    if sensitivity
        y0 = zeros(length(phase.species)*(length(phase.species)+1+length(phase.reactions)))
    else
        y0 = zeros(length(phase.species))
    end

    y0[phase.species[1].index:phase.species[end].index] = ns

    if length(constantspecies) > 0
        spcnames = getfield.(phase.species,:name)
        constspcinds = [findfirst(isequal(k),spcnames) for k in constantspecies]
    else
        constspcinds = Array{Int64,1}()
    end
    if sparse
        jacobian=zeros(typeof(V),length(phase.species)+1,length(phase.species)+1)
    else
        jacobian=zeros(typeof(V),length(phase.species)+1,length(phase.species)+1)
    end
    rxnarray = getreactionindices(phase)
    return ParametrizedTPDomain(phase,interfaces,SVector(phase.species[1].index,phase.species[end].index),constspcinds,
    Tfcn,Pfcn,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0)), y0
end
export ParametrizedTPDomain

@with_kw struct ParametrizedVDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,Q<:AbstractArray} <: AbstractVariableKDomain
    phase::N
    interfaces::Array{AbstractInterface,1} = Array{AbstractInterface,1}()
    indexes::Q #assumed to be in ascending order
    constantspeciesinds::Array{S,1}
    V::Function
    rxnarray::Array{UInt16,2}
    jacobian::Array{W,2}
    sensitivity::Bool = false
    jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
    t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
end
function ParametrizedVDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,Any},constantspecies::Array{X2,1}=Array{String,1}(),
    sparse::Bool=false,sensitivity::Bool=false) where {X,X2,E<:Real,Z<:IdealGas,Q<:AbstractInterface}

    #set conditions and initialconditions
    T = 0.0
    P = 0.0
    V = 0.0
    ts = Array{Float64,1}()
    ns = zeros(length(phase.species))
    spnames = [x.name for x in phase.species]
    @assert "V" in keys(initialconds)
    for (key,val) in initialconds
        if key == "T"
            T = val
        elseif key == "P"
            P = val
        elseif key == "V"
            V = val
        elseif key == "ts"
            ts = val
        else
            ind = findfirst(isequal(key),spnames)
            @assert typeof(ind)<: Integer  "$key not found in species list: $spnames"
            ns[ind] = val
        end
    end
    @assert isa(V,Function) || isa(V,AbstractArray)
    if isa(V,AbstractArray)
        q = Spline1D(ts,V;k=3,s=1e-11)
        Vfcn = f(x::Float64) = q(x)
    elseif isa(V,Function)
        Vfcn = V
    else
        throw(error("ParametrizedVDomain must take \"V\" as a function or if an array of times for \"ts\" is supplied as an array of volumes"))
    end
    N = sum(ns)
    if T == 0.0
        T = P*Vfcn(0.0)/(R*N)
    elseif P == 0.0
        P = N*R*T/Vfcn(0.0)
    else
        ns *= (P*Vfcn(0.0)/(R*T))/sum(ns) #automatically scale down moles if pressure specified
    end
    if sensitivity
        y0 = vcat(ns,T,zeros((length(ns)+1)*(length(ns)+length(phase.reactions))))
    else
        y0 = vcat(ns,T)
    end
    if length(constantspecies) > 0
        spcnames = getfield.(phase.species,:name)
        constspcinds = [findfirst(isequal(k),spcnames) for k in constantspecies]
    else
        constspcinds = Array{Int64,1}()
    end
    if sparse
        jacobian=zeros(typeof(T),length(phase.species)+1,length(phase.species)+1)
    else
        jacobian=zeros(typeof(T),length(phase.species)+1,length(phase.species)+1)
    end
    rxnarray = getreactionindices(phase)
    return ParametrizedVDomain(phase,interfaces,SVector(phase.species[1].index,phase.species[end].index,phase.species[end].index+1),constspcinds,
    Vfcn,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0)), y0
end
export ParametrizedVDomain
@with_kw struct ConstantTVDomain{N<:AbstractPhase,S<:Integer,W<:Real, W2<:Real, I<:Integer, Q<:AbstractArray} <: AbstractConstantKDomain
    phase::N
    interfaces::Array{AbstractInterface,1} = Array{AbstractInterface,1}()
    indexes::Q #assumed to be in ascending order
    constantspeciesinds::Array{S,1}
    T::W
    V::W
    kfs::Array{W,1}
    krevs::Array{W,1}
    efficiencyinds::Array{I,1}
    Gs::Array{W,1}
    rxnarray::Array{UInt16,2}
    mu::W = 0.0
    diffusivity::Array{W,1} = Array{Float64,1}()
    jacobian::Array{W,2} = Array{Float64,2}(undef,(0,0))
    sensitivity::Bool = false
    jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
    t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
end
function ConstantTVDomain(;phase::Z,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,E},constantspecies::Array{X2,1}=Array{String,1}(),
    sparse=false,sensitivity=false) where {E,X,X2, Z<:AbstractPhase,Q<:AbstractInterface,W<:Real}
    #set conditions and initialconditions
    T = 0.0
    V = 0.0
    P = 1.0e9
    if sensitivity
        y0 = zeros(length(phase.species)*(length(phase.species)+1+length(phase.reactions)))
    else
        y0 = zeros(length(phase.species))
    end
    spnames = [x.name for x in phase.species]
    for (key,val) in initialconds
        if key == "T"
            T = val
        elseif key == "P"
            P = val
        elseif key == "V"
            V = val
        else
            ind = findfirst(isequal(key),spnames)
            @assert typeof(ind)<: Integer  "$key not found in species list: $spnames"
            y0[ind] = val
        end
    end
    @assert T != 0.0
    @assert V != 0.0
    ns = y0
    N = sum(ns)

    if length(constantspecies) > 0
        spcnames = getfield.(phase.species,:name)
        constspcinds = [findfirst(isequal(k),spcnames) for k in constantspecies]
    else
        constspcinds = Array{Int64,1}()
    end
    efficiencyinds = [rxn.index for rxn in phase.reactions if typeof(rxn.kinetics)<:AbstractFalloffRate && length(rxn.kinetics.efficiencies) > 0]
    Gs = calcgibbs(phase,T)
    if :solvent in fieldnames(typeof(phase)) && typeof(phase.solvent) != EmptySolvent
        mu = phase.solvent.mu(T)
    else
        mu = 0.0
    end
    if phase.diffusionlimited
        diffs = [x(T=T,mu=mu,P=P) for x in getfield.(phase.species,:diffusion)]
    else
        diffs = Array{Float64,1}()
    end
    P = 1.0e9  #essentiallly assuming this is a liquid
    C = N/V
    kfs,krevs = getkfkrevs(phase=phase,T=T,P=P,C=C,N=N,ns=ns,Gs=Gs,diffs=diffs,V=V)
    if sparse
        jacobian=zeros(typeof(T),length(phase.species),length(phase.species))
    else
        jacobian=zeros(typeof(T),length(phase.species),length(phase.species))
    end
    rxnarray = getreactionindices(phase)
    return ConstantTVDomain(phase,interfaces,SVector(phase.species[1].index,phase.species[end].index),constspcinds,
        T,V,kfs,krevs,efficiencyinds,Gs,rxnarray,mu,diffs,jacobian,sensitivity,MVector(false),MVector(0.0)), y0
end
export ConstantTVDomain

@with_kw struct ParametrizedTConstantVDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,Q<:AbstractArray} <: AbstractVariableKDomain
    phase::N
    interfaces::Array{AbstractInterface,1} = Array{AbstractInterface,1}()
    indexes::Q #assumed to be in ascending order
    constantspeciesinds::Array{S,1}
    T::Function
    V::W
    rxnarray::Array{UInt16,2}
    jacobian::Array{W,2}
    sensitivity::Bool = false
    jacuptodate::MArray{Tuple{1},Bool,1,1}=MVector(false)
    t::MArray{Tuple{1},W2,1,1}=MVector(0.0)
end
function ParametrizedTConstantVDomain(;phase::IdealDiluteSolution,interfaces::Array{Q,1}=Array{EmptyInterface,1}(),initialconds::Dict{X,X3},constantspecies::Array{X2,1}=Array{String,1}(),
    sparse::Bool=false,sensitivity::Bool=false) where {X,X2,X3,Q<:AbstractInterface}
    #set conditions and initialconditions
    T = 0.0
    P = 0.0
    V = 0.0
    ts = Array{Float64,1}()
    ns = zeros(length(phase.species))
    spnames = [x.name for x in phase.species]
    @assert "V" in keys(initialconds)
    for (key,val) in initialconds
        if key == "T"
            T = val
        elseif key == "P"
            P = val
        elseif key == "V"
            V = val
        elseif key == "ts"
            ts = val
        else
            ind = findfirst(isequal(key),spnames)
            @assert typeof(ind)<: Integer  "$key not found in species list: $spnames"
            ns[ind] = val
        end
    end
    if isa(T,AbstractArray)
        q = Spline1D(ts,T;k=3,s=1e-11)
        Tfcn = f(x::Float64) = q(x)
    elseif isa(T,Function)
        Tfcn = T
    else
        throw(error("ParametrizedTConstantVDomain must take \"T\" as a function or if an array of times for \"ts\" is supplied as an array of volumes"))
    end
    N = sum(ns)
    if P == 0.0
        P = 1e8
    else
        throw(error("ParametrizedTConstantVDomain cannot specify P"))
    end
    if sensitivity
        y0 = zeros(length(phase.species)*(length(phase.species)+1+length(phase.reactions)))
    else
        y0 = zeros(length(phase.species))
    end
    y0[phase.species[1].index:phase.species[end].index] = ns
    if length(constantspecies) > 0
        spcnames = getfield.(phase.species,:name)
        constspcinds = [findfirst(isequal(k),spcnames) for k in constantspecies]
    else
        constspcinds = Array{Int64,1}()
    end
    if sparse
        jacobian=zeros(typeof(V),length(phase.species)+1,length(phase.species)+1)
    else
        jacobian=zeros(typeof(V),length(phase.species)+1,length(phase.species)+1)
    end
    rxnarray = getreactionindices(phase)
    return ParametrizedTConstantVDomain(phase,interfaces,SVector(phase.species[1].index,phase.species[end].index),constspcinds,
    Tfcn,V,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0)), y0
end
export ParametrizedTConstantVDomain

@inline function calcthermo(d::ConstantTPDomain{W,Y},y::J,t::Q) where {W<:IdealGas,Y<:Integer,J<:AbstractArray{Float64,1},Q<:Float64}
    if t != d.t[1]
        d.t[1] = t
        d.jacuptodate[1] = false
    end
    ns = y[d.indexes[1]:d.indexes[2]]
    N = sum(ns)
    V = N*d.T*R/d.P
    cs = ns./V
    C = N/V
    for ind in d.efficiencyinds #efficiency related rates may have changed
        d.kfs[ind],d.krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,d.Gs,d.diffusivity,V)
    end
    return ns,cs,d.T,d.P,V,C,N,d.mu,d.kfs,d.krevs,[],[],[],[],0.0
end

@inline function calcthermo(d::ConstantTPDomain{W,Y},y::J,t::Q) where {W<:IdealGas,Y<:Integer,J<:AbstractArray,Q<:Real}
    if t != d.t[1]
        if isa(t,Float64)
            d.t[1] = t
        end
        d.jacuptodate[1] = false
    end
    ns = y[d.indexes[1]:d.indexes[2]]
    N = sum(ns)
    V = N*d.T*R/d.P
    cs = ns./V
    C = N/V
    kfs = convert(typeof(y),copy(d.kfs))
    krevs = convert(typeof(y),copy(d.krevs))
    for ind in d.efficiencyinds #efficiency related rates may have changed
        kfs[ind],krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,d.Gs,d.diffusivity,V)
    end
    return ns,cs,d.T,d.P,V,C,N,d.mu,kfs,krevs,[],[],[],[],0.0
end

@inline function calcthermo(d::ConstantVDomain{W,Y},y::J,t::Q) where {W<:IdealGas,Y<:Integer,J<:AbstractArray,Q<:Real}
    if t != d.t[1]
        d.t[1] = t
        d.jacuptodate[1] = false
    end
    ns = y[d.indexes[1]:d.indexes[2]]
    T = y[d.indexes[3]]
    N = sum(ns)
    cs = ns./d.V
    C = N/d.V
    P = C*R*T
    Gs = zeros(length(d.phase.species))
    Us = zeros(length(d.phase.species))
    Cvave = 0.0
    @simd for i = 1:length(d.phase.species)
        @inbounds cpdivR,hdivRT,sdivR = calcHSCpdless(d.phase.species[i].thermo,T)
        @fastmath @inbounds Gs[i] = (hdivRT-sdivR)*R*T
        @fastmath @inbounds Us[i] = (hdivRT-1.0)*R*T
        @fastmath @inbounds Cvave += cpdivR*ns[i]
    end
    @fastmath Cvave *= R/N
    @fastmath Cvave -= R
    if d.phase.diffusionlimited
        diffs = getfield.(d.phase.species,:diffusion)(T=T,mu=mu,P=P)
    else
        diffs = Array{Float64,1}()
    end
    kfs,krevs = getkfkrevs(phase=d.phase,T=T,P=P,C=C,N=N,ns=ns,Gs=Gs,diffs=diffs,V=d.V)
    return ns,cs,T,P,d.V,C,N,0.0,kfs,krevs,[],Us,Gs,diffs,Cvave
end

@inline function calcthermo(d::ParametrizedVDomain{W,Y},y::J,t::Q) where {W<:IdealGas,Y<:Integer,J<:AbstractArray,Q<:Real}
    if t != d.t[1]
        d.t[1] = t
        d.jacuptodate[1] = false
    end
    V = d.V(t)
    ns = y[d.indexes[1]:d.indexes[2]]
    T = y[d.indexes[3]]
    N = sum(ns)
    cs = ns./V
    C = N/V
    P = C*R*T
    Gs = zeros(length(d.phase.species))
    Us = zeros(length(d.phase.species))
    Cvave = 0.0
    @simd for i = 1:length(d.phase.species)
        @inbounds cpdivR,hdivRT,sdivR = calcHSCpdless(d.phase.species[i].thermo,T)
        @fastmath @inbounds Gs[i] = (hdivRT-sdivR)*R*T
        @fastmath @inbounds Us[i] = (hdivRT-1.0)*R*T
        @fastmath @inbounds Cvave += cpdivR*ns[i]
    end
    @fastmath Cvave *= R/N
    @fastmath Cvave -= R
    if d.phase.diffusionlimited
        diffs = getfield.(d.phase.species,:diffusion)(T=T,mu=mu,P=P)
    else
        diffs = Array{Float64,1}()
    end
    kfs,krevs = getkfkrevs(phase=d.phase,T=T,P=P,C=C,N=N,ns=ns,Gs=Gs,diffs=diffs,V=V)
    return ns,cs,T,P,V,C,N,0.0,kfs,krevs,[],Us,Gs,diffs,Cvave
end

@inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y},y::J,t::Q) where {W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q<:Real}
    if t != d.t[1]
        d.t[1] = t
        d.jacuptodate[1] = false
    end
    V = d.V
    T = d.T(t)
    ns = y[d.indexes[1]:d.indexes[2]]
    N = sum(ns)
    cs = ns./V
    C = N/V
    P = C*R*T
    Gs = zeros(length(d.phase.species))
    mu = d.phase.solvent.mu(T)
    for i = 1:length(d.phase.species)
        @inbounds cpdivR,hdivRT,sdivR = calcHSCpdless(d.phase.species[i].thermo,T)
        @fastmath @inbounds Gs[i] = (hdivRT-sdivR)*R*T
    end
    if d.phase.diffusionlimited
        diffs = [x(T=T,mu=mu,P=P) for x in getfield.(d.phase.species,:diffusion)]
    else
        diffs = Array{Float64,1}()
    end
    kfs,krevs = getkfkrevs(phase=d.phase,T=T,P=P,C=C,N=N,ns=ns,Gs=Gs,diffs=diffs,V=V)
    return ns,cs,T,P,V,C,N,0.0,kfs,krevs,[],[],Gs,diffs,0.0
end

@inline function calcthermo(d::ParametrizedTPDomain{W,Y},y::J,t::Q) where {W<:IdealGas,Y<:Integer,J<:AbstractArray,Q<:Real}
    if t != d.t[1]
        d.t[1] = t
        d.jacuptodate[1] = false
    end
    T = d.T(t)
    @assert T < 10000.0
    P = d.P(t)
    ns = y[d.indexes[1]:d.indexes[2]]
    N = sum(ns)
    V = N*T*R/P
    cs = ns./V
    C = N/V
    P = C*R*T
    Gs = zeros(length(d.phase.species))
    Us = zeros(length(d.phase.species))
    Cvave = 0.0
    @simd for i = 1:length(d.phase.species)
        @inbounds cpdivR,hdivRT,sdivR = calcHSCpdless(d.phase.species[i].thermo,T)
        @fastmath @inbounds Gs[i] = (hdivRT-sdivR)*R*T
    end
    @fastmath Cvave *= R/N
    @fastmath Cvave -= R
    if d.phase.diffusionlimited
        diffs = getfield.(d.phase.species,:diffusion)(T=T,mu=mu,P=P)
    else
        diffs = Array{Float64,1}()
    end
    kfs,krevs = getkfkrevs(phase=d.phase,T=T,P=P,C=C,N=N,ns=ns,Gs=Gs,diffs=diffs,V=V)
    return ns,cs,T,P,V,C,N,0.0,kfs,krevs,[],[],Gs,diffs,0.0
end

@inline function calcthermo(d::ConstantTVDomain{W,Y},y::J,t::Q) where {W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q<:Real}
    if t != d.t[1]
        d.t[1] = t
        d.jacuptodate[1] = false
    end
    ns = y[d.indexes[1]:d.indexes[2]]
    N = sum(ns)
    cs = ns./d.V
    C = N/d.V
    P = 1.0e9
    return ns,cs,d.T,P,d.V,C,N,d.mu,d.kfs,d.krevs,[],[],[],[],0.0
end
export calcthermo

@inline function calcdomainderivatives!(d::Q,dydt::Array{Z7,1};t::Z10,T::Z4,P::Z9,Us::Array{Z,1},V::Z2,C::Z3,ns::Array{Z5,1},N::Z6,Cvave::Z8) where {Q<:AbstractDomain,Z10,Z9,Z8<:Real,Z7<:Real,W<:IdealGas,Y<:Integer,Z6,Z,Z2,Z3,Z4,Z5<:Real}
    for ind in d.constantspeciesinds #make dydt zero for constant species
        @inbounds dydt[ind] = 0.0
    end
end

@inline function calcdomainderivatives!(d::ConstantVDomain{W,Y},dydt::Array{K,1};t::Z10,T::Z4,P::Z9,Us::Array{Z,1},V::Z2,C::Z3,ns::Array{Z5,1},N::Z6,Cvave::Z7) where {Z10,Z9,W<:IdealGas,Z7<:Real,K<:Real,Y<:Integer,Z6,Z,Z2,Z3,Z4,Z5<:Real}
    @views @fastmath @inbounds dydt[d.indexes[3]] = -dot(Us,dydt[d.indexes[1]:d.indexes[2]])/(N*Cvave) #divide by V to cancel ωV to ω
    for ind in d.constantspeciesinds #make dydt zero for constant species
        @inbounds dydt[ind] = 0.0
    end
end

@inline function calcdomainderivatives!(d::ParametrizedVDomain{W,Y},dydt::Array{K,1};t::Z10,T::Z4,P::Z9,Us::Array{Z,1},V::Z2,C::Z3,ns::Array{Z5,1},N::Z6,Cvave::Z7) where {Z10,Z9,W<:IdealGas,Z7<:Real,K<:Real,Y<:Integer,Z6,Z,Z2,Z3,Z4,Z5<:Real}
    @views @fastmath @inbounds dydt[d.indexes[3]] = (-dot(Us,dydt[d.indexes[1]:d.indexes[2]])-P*Calculus.derivative(d.V,t))/(N*Cvave) #divide by V to cancel ωV to ω

    for ind in d.constantspeciesinds #make dydt zero for constant species
        @inbounds dydt[ind] = 0.0
    end
end
export calcdomainderivatives!

function getreactionindices(ig::Q) where {Q<:AbstractPhase}
    arr = zeros(UInt16,(6,length(ig.reactions)))
    for (i,rxn) in enumerate(ig.reactions)
        arr[1:length(rxn.reactantinds),i] = rxn.reactantinds
        arr[4:length(rxn.productinds)+3,i] = rxn.productinds
    end
    return arr
end
export getreactionindices


text:
for TYP in (:ReverseOrientation,:Reverse)
    WRAP = Meta.parse(string(TYP)*"Wrapper")
    @eval begin
        abstract type $TYP{T} <: Operator{T} end

        struct $WRAP{OS,T} <: Operator{T}
            op::OS
        end

        $WRAP(op::Operator) = $WRAP{typeof(op),eltype(op)}(op)
        convert(::Type{Operator{T}},op::$TYP) where {T} = $TYP{T}()
        convert(::Type{Operator{T}},op::$WRAP) where {T} = $WRAP(Operator{T}(op.op))::Operator{T}

        @wrapper $WRAP
    end
end


text:
module p3

# Return boundary conditions values and indices, ϕb and β, for problem 2,
# with N subdivisions
function boundaryconditions3(N::Integer)
    ϕb = zeros(N+1,N+1)                 # boundary conditions
    β  = zeros(Int,N+1,N+1)             # boundary condition locations

    # Conditions on the perimeter
    ϕb[:,1] = ϕb[:,N+1] = ϕb[1,:] = ϕb[N+1,:] = 0.0
    β[:,1]  = β[:,N+1]  = β[1,:]  = β[N+1,:]  = 1

    # Conditions on the interior
    # xmin < x < xmax  →  N*xmin + 1 < i < N*xmax + 1
    xmin = 0.5
    xmax = 0.75
    ymin = 0.625
    ymax = 0.875

    xs = int([ceil(N*xmin+1):floor(N*xmax+1)])
    ys = int([ceil(N*ymin+1):floor(N*ymax+1)])

    ϕb[xs,ys] = 0.0
    β[xs,ys]  = 1

    return ϕb, β
end

# Calculate the potential
function V(x::Real, y::Real, N::Integer, q::Real)

    a = b = 1:(N+1)
    a .-= (N*x + 0.5)           # avoid singularity
    b .-= (N*y + 0.5)
    a .*= a
    b .*= b
    return -N*q ./ √(a .+ b')

end

V(N::Integer, q::Real) = V(0.5, 0.25, N, q)

function minor(A::Matrix, β::Vector)

    size(A)[1] == size(A)[2] || error("Matrix must be square")

    N = sum(1 .- β)             # dimensions of new matrix

    B = Array(typeof(A[1]),N,N)

    ii = 1
    jj = 1

    for j in 1:size(A)[1]
        β[j] == 1 && continue
        for i in 1:size(A)[2]
            β[i] == 1 && continue
            B[ii,jj] = A[i,j]
            ii += 1
        end
        jj +=1
        ii = 1
    end

    return B
end

# repopulate a vector v with zeros for fixed parameters
function addwater(v::Vector, β::Vector)

    vf = Array(typeof(v[1]),length(β))
    ii = 1

    length(v) == sum(1-β) || error("v and β have different num free params")

    for i in 1:length(β)
        if β[i] == 1
            vf[i] = 0               # put zero for a fixed value
        elseif β[i] == 0
            vf[i] = v[ii]           # fill in the free value
            ii += 1                 # prep the next free value
        else
            error("β can only have 0 or 1 values")
        end
    end
    return vf

end

end


text:
@doc raw"""
    SIR(du,u,p,t)

Define classic SIR(suspected-infected-recovered) model.

Parameters: (birth rate, natural death rate, disaese induced death rate, population, infection rate, recovery rate)

```math
\begin{aligned}
& \frac{\rm{d}S}{\rm{dt}} = \Lambda -\beta S I/N - d S,\\
&\frac{\rm{d}I}{\rm{dt}} = \beta  S I/N - \gamma I - d I - \alpha I,\\
&\frac{\rm{d}R}{\rm{dt}} = \gamma I - d R,\\
\end{aligned}
```
"""
function SIRbasic(du, u, p, t)
    Λ, d, α, N, β, γ = p
    pop = (Λ > 0 & d > 0) ? Λ / d : N
    S, I, R = u
    du[1] = Λ - β * S * I / pop - d * S
    du[2] = β * S * I / pop - γ * I - d * I - α * I
    du[3] = γ * I - d * R
end

@doc raw"""
    SEIR(du,u,p,t)

Define classic SEIR(suspected-exposed-infected-recovered) model.

Parameters: (birth rate, natural death rate, disaese induced death rate, population, infection rate, recovery rate, incubation period,
exposed decreasing infection ratio)

```math
\begin{aligned}
& \frac{\rm{d}S}{\rm{dt}} = \Lambda -\beta S (I+k_{E}E)/N - d S,\\
&\frac{\rm{d}E}{\rm{dt}} = \beta  S (I+k_{E}E)/N - \sigma E -d E,\\
&\frac{\rm{d}I}{\rm{dt}} =\sigma E - \gamma I - d I - \alpha I,\\
&\frac{\rm{d}R}{\rm{dt}} = \gamma I - d R,\\
\end{aligned}
```
"""
function SEIRbasic(du, u, p, t)
    Λ, d, α, N, β, γ, σ, ke = p
    pop = (Λ > 0 & d > 0) ? Λ / d : N
    S, E, I, R = u
    du[1] = Λ - β * S * (I + ke * E) / pop - d * S
    du[2] = β * S * (I + ke * E) / pop - σ * E - d * E
    du[3] = σ * E - γ * I - d * I - α * I
    du[4] = γ * I - d * R
end

@doc raw"""
    SEIAR(du,u,p,t)

Define classic SEIR(suspected-exposed-infected-asymptomatic-recovered) model.

Parameters: (birth rate, natural death rate, disaese induced death rate, population, infection rate, recovery rate, incubation period,
exposed decreasing infection ratio, asymptomatic infection decreasing infection ratio, asymptomatic rate,recovery rate of asymptomatic)

```math
\begin{aligned}
& \frac{\rm{d}S}{\rm{dt}} = \Lambda -\beta S (I+k_{E}E+k_{A}A)/N - d S,\\
&\frac{\rm{d}E}{\rm{dt}} = \beta  S (I+k_{E}E++k_{A}A)/N - \sigma E -d E,\\
&\frac{\rm{d}I}{\rm{dt}} =(1-\rho)\sigma E - \gamma I - d I - \alpha I,\\
&\frac{\rm{d}A}{\rm{dt}} =\rho\sigma E - \gamma_{A} A - d A,\\
&\frac{\rm{d}R}{\rm{dt}} = \gamma I +\gamma_{A} A - d R,\\
\end{aligned}
```
"""
function SEIARbasic(du, u, p, t)
    Λ, d, α, N, β, γ, σ, ke, ka, ρ, γA= p
    pop = (Λ > 0 & d > 0) ? Λ / d : N
    S, E, I, A, R = u
    du[1] = Λ - β * S * (I + ke * E + ka * A) / pop - d * S
    du[2] = β * S * (I + ke * E + ka * A) / pop - σ * E - d * E
    du[3] = (1-ρ)*σ * E - γ * I - d * I - α * I
    du[3] = ρ * σ * E - γA * A - d * A
    du[4] = γ * I +  γA * A - d * R
end

@doc raw"""
    SIS(du,u,p,t)

Define classic SIR(suspected-infected-suspected) model.

Parameters: (population, infection rate, recovery rate)

```math
\begin{aligned}
& \frac{\rm{d}S}{\rm{dt}} = -\beta S I/N + \gamma I,\\
&\frac{\rm{d}I}{\rm{dt}} = \beta  S I/N - \gamma I,\\
\end{aligned}
```
"""
function SISbasic(du, u, p, t)
    pop, β, γ = p
    S, I = u
    du[1] = - β * S * I / pop + γ * I
    du[2] = β * S * I / pop - γ * I
end


text:
import GenDirectionalStats: vmf_rot3
import Rotations

struct GaussianVMFMixture <: Gen.Distribution{Pose} end
const gaussianVMF_mixture = GaussianVMFMixture()
function Gen.random(
    ::GaussianVMFMixture,
    centers::Array{Pose},
    positionStdev::Real,
    orientationKappa::Real,
)::Pose
    N = length(centers)
    i = categorical(fill(1.0/N, N))
    rot = Gen.random(vmf_rot3, Rotations.UnitQuaternion(1.0, 0.0, 0.0, 0.0), orientationKappa)
    offset = Pose(
        [
            Gen.normal(0, positionStdev),
            Gen.normal(0, positionStdev),
            Gen.normal(0, positionStdev),
        ],
        rot,
    )

    Pose(centers[i].pos + offset.pos, centers[i].orientation * offset.orientation)
end

function Gen.logpdf(
    ::GaussianVMFMixture,
    val::Pose,
    centers::Array{Pose},
    positionStdev::Real,
    orientationKappa::Real,
)
    N = length(centers)
    logpdfs = [
        logpdf(gaussianVMF, val, c, positionStdev, orientationKappa) for c in centers
    ]
    logsumexp(logpdfs) - log(N)
end
Gen.has_output_grad(::GaussianVMFMixture) = false
Gen.has_argument_grads(::GaussianVMFMixture) = (false, false, false)
function Gen.logpdf_grad(
    ::GaussianVMFMixture,
    center::Pose,
    positionStdev::Real,
    orientationKappa::Real,
)
    (nothing, nothing, nothing, nothing)
end

export gaussianVMF_mixture

struct GaussianVMFWeightedMixture <: Gen.Distribution{Pose} end
const gaussianVMF_weighted_mixture = GaussianVMFWeightedMixture()
function Gen.random(
    ::GaussianVMFWeightedMixture,
    centers::Array{Pose},
    weights::Array{<:Real},
    positionStdev::Real,
    orientationKappa::Real,
)::Pose
    N = length(centers)
    i = categorical(weights)
    rot = Gen.random(vmf_rot3, Rotations.UnitQuaternion(1.0, 0.0, 0.0, 0.0), orientationKappa)
    offset = Pose(
        [
            Gen.normal(0, positionStdev),
            Gen.normal(0, positionStdev),
            Gen.normal(0, positionStdev),
        ],
        rot,
    )

    Pose(centers[i].pos + offset.pos, centers[i].orientation * offset.orientation)
end

function Gen.logpdf(
    ::GaussianVMFWeightedMixture,
    val::Pose,
    centers::Array{Pose},
    weights::Array{<:Real},
    positionStdev::Real,
    orientationKappa::Real,
)
    N = length(centers)
    logpdfs = [
        log(w) + logpdf(gaussianVMF, val, c, positionStdev, orientationKappa) for (w,c) in zip(weights,centers)
    ]
    logsumexp(logpdfs) - log(N)
end
Gen.has_output_grad(::GaussianVMFWeightedMixture) = false
Gen.has_argument_grads(::GaussianVMFWeightedMixture) = (false, false, false)
function Gen.logpdf_grad(
    ::GaussianVMFWeightedMixture,
    center::Pose,
    weights::Array{Real},
    positionStdev::Real,
    orientationKappa::Real,
)
    (nothing, nothing, nothing, nothing)
end

export gaussianVMF_weighted_mixture


text:
function variablecost(m::JuMP.AbstractModel, devices::Array{PowerSystems.InterruptibleLoad,1}, device_formulation::Type{D}, system_formulation::Type{S}) where {D <: AbstractControllableLoadForm, S <: PM.AbstractPowerFormulation}

    p_cl = m[:p_cl]
    time_index = m[:p_cl].axes[2]
    name_index = m[:p_cl].axes[1]

    var_cost = AffExpr()

    for  (ix, name) in enumerate(name_index)
         c = gencost(m, p_cl[name,:], -1*devices[ix].sheddingcost)
        JuMP.add_to_expression!(var_cost,c)
    end

    return var_cost

end


text:
module DynamicFactorModels

#using GLM
using DataArrays
using DataFrames
using Distributions
#using Gadfly
# using DimensionalityReduction
#using GLMNet
#using Fred

# text represenation
import StatsBase.predict
#show()

export DynamicFactorModel, predict, calculate_factors,
    lag_vector, factor_model_DGP, normalize, pseudo_out_of_sample_forecasts,
    targeted_predictors,
    wild_bootstrap, residual_bootstrap,
    LR_test, LM_test, Wald_test

include("utils.jl")
include("DynamicFactorModel.jl")
include("targeted_predictors.jl")
include("bootstrap.jl")
include("chowtest.jl")

end # module


text:
################################################################################
# returns dense idenity matrix of eltype T, used in BFGS
identity(T::Type, n::Int) = Matrix(one(T) * I(n))
identity(x::AbstractVecOrMat) = identity(eltype(x), size(x, 1))

################################################################################
using DataStructures: CircularBuffer
# convenience for circular buffer
# does in place copy of x in front of buffer, important if cb is buffer of arrays
# used in LBFGS
function copyfirst!(cb::CircularBuffer, x)
    cb.first = (cb.first == 1 ? cb.capacity : cb.first - 1)
    if length(cb) < cb.capacity
        cb.length += 1
    end
    copy!(cb.buffer[cb.first], x)
end

# TODO: make negative shifts efficient
# TODO: if shifts > c.length ÷ 2,
# used in LBFGS
# Base.circshift
function _circshift!(cb::CircularBuffer, shifts::Int = 1)
    shifts = mod(shifts, cb.length) # since identity if shifts = cb.length
    if shifts > 0
        if cb.length < cb.capacity
            first, last = cb.first, cb.first + cb.length
            for i in 1:shifts
                src = mod1((i-1) + first, cb.capacity)
                dest = mod1(i + last, cb.capacity)
                copy!(cb.buffer[dest], cb.buffer[src])
            end
        end
        cb.first = mod1(cb.first + shifts, cb.capacity)
    end
    return cb
end

######################### Lazy Difference Vector ###############################
# IDEA: is it worth using the LazyArrays package?
# could be replaced by LazyVector(applied(-, x, y)), which is neat.
# lazy difference between two vectors, has no memory footprint
struct LazyDifference{T, U, V} <: AbstractVector{T}
    x::U
    y::V
    function LazyDifference(x, y)
        length(x) == length(y) || throw(DimensionMismatch("x and y do not have the same length: $(length(x)) and $(length(y))."))
        T = promote_type(eltype(x), eltype(y))
        new{T, typeof(x), typeof(y)}(x, y)
    end
end

difference(x::Number, y::Number) = x-y # avoid laziness for scalars
difference(x::AbstractVector, y::AbstractVector) = LazyDifference(x, y)
difference(x::Tuple, y::Tuple) = LazyDifference(x, y)
# this creates a type instability:
Base.size(d::LazyDifference) = (length(d.x),)
Base.getindex(d::LazyDifference, i::Integer) = d.x[i]-d.y[i]

######################### solves with eigen factorization, used for SaddleFreeNewton
function LinearAlgebra.ldiv!(y::AbstractVector, E::Eigen, x::AbstractVector)
     ldiv!!(y, E, copy(x))
end
# WARNING: uses x as a temporary
function ldiv!!(y::AbstractVector, E::Eigen, x::AbstractVector) # temporary
    mul!(y, E.vectors', x)
    @. x = y / E.values
    mul!(y, E.vectors, x)
end


text:
module PopPKBase

include("post.jl")

green = "#91BE33"
lightgreen = "#32A03E"
turqoise = "#18C8D0"
teal = "#16A0AC"
purple = "#6C62D0"
violet = "#A238BA"
yosemite = "#8A8CE3"
blue = "#178CCB"
navy = "#052049"
magenta = "#561038"
red = "#E61048"
orange = "#F48024"

end


text:
# testset for using them in loops
@testset "Loop Induction Variables" begin
	@show @__LINE__
    f(x) = cos(x) * log(x)
    function avxmax(v)
        max_x = -Inf
        @avx for i ∈ eachindex(v)
            x = f(i)
            max_x = max(max_x, x)
        end
        max_x
    end
    function avxextrema(v)
        max_x = -Inf
        min_x = Inf
        @avx for i ∈ eachindex(v)
            x = f(i)
            max_x = max(max_x, x)
            min_x = min(min_x, x)
        end
        min_x, max_x
    end

    v = 1:19
    minref, maxref = extrema(f, v)
    @test maxref ≈ avxmax(v)
    minavx, maxavx = avxextrema(v);
    @test minref ≈ minavx
    @test maxref ≈ maxavx
end


text:
using test
using Test

@testset "test.jl" begin
    # Write your tests here.
    x = 2
    y = 2
    @test test.sum_values(x,y) == 4
end


text:
using Gadfly
using Color
using DataFrames

Season = [div(j-1,8)+1 for j in 1:160];
Drops = [mod(j-1,8)+1 for j in 1:160];
Count = rand(1:20, 160) .+ 30;
df = DataFrame(Season=Season, Drops=Drops, Count=Count);

# Diferite variante

plot(df, x=:Season, y=:Count, color=:Drops, Geom.bar(position=:stack), Guide.yticks(ticks=[0:50:400]))

#=
plot(df, x=:Season, y=:Count, color=:Drops, Geom.bar(position=:stack), Guide.yticks(ticks=[0:50:400]), Scale.ContinuousColorScale(p -> RGB(0, 0, 0.2+0.7*(1-p))))
=#
#=
plot(df, x=:Season, y=:Count, color=:Drops, Geom.bar(position=:stack), Guide.yticks(ticks=[0:50:400]), Scale.ContinuousColorScale(p -> RGB(0.2+ 0.7*(1-p),0.2+0.7*p,0.2+0.7*(1-p))))
=#
#=
plot(df, x=:Season, y=:Count, color=:Drops, Geom.bar(position=:stack), Guide.yticks(ticks=[0:50:400]), Scale.ContinuousColorScale(p -> RGB(0.2+ 0.7*(1-p),0.2+0.7*p,0.2+0.7*(1-p))))
=#

max_stars_repo_name: tudor-berariu/phd

max_issues_repo_name: tudor-berariu/phd

max_forks_repo_name: tudor-berariu/phd


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Based on Models written by Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
# The complete documentation for these models needs to be updated !!!
type User_tray
	User_tray()=begin
		new(
			User_trayBasic(),
			dens_mass(),
			dens_mass(),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Temporary variable of Roffels liquid flow equation"
			)),
			[
				:(rhoL = _base_1.PP.LiquidDensity(_base_1.OutletL.T, _base_1.OutletL.P, _base_1.OutletL.z)),
				:(rhoV = _base_1.PP.VapourDensity(_base_1.InletV.T, _base_1.InletV.P, _base_1.InletV.z)),
			],
			[
				"Liquid Density","Vapour Density",
			],
			[:rhoL,:rhoV,:btemp,]
		)
	end
	_base_1::User_trayBasic
	rhoL::dens_mass
	rhoV::dens_mass
	btemp::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export User_tray
function setEquationFlow(in::User_tray)
	addEquation(1)
	addEquation(2)
end
function atributes(in::User_tray,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/Tray"
	fields[:Brief]="Description of a column tray."
	fields[:Info]="To be updated"
	drive!(fields,_)
	return fields
end
User_tray(_::Dict{Symbol,Any})=begin
	newModel=User_tray()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# # Wrapper blocks

abstract type WrapperBlock <: AbstractBlock end

wrapped(w::WrapperBlock) = w.block
function setwrapped(w::WrapperBlock, b)
    return Setfield.@set w.block = b
end
mockblock(w::WrapperBlock) = mockblock(wrapped(w))
checkblock(w::WrapperBlock, data) = checkblock(wrapped(w), data)

# If not overwritten, encodings are applied to the wrapped block

function encodedblock(enc::Encoding, wrapper::WrapperBlock)
    inner = encodedblock(enc, wrapped(wrapper))
    return isnothing(inner) ? nothing : setwrapped(wrapper, inner)
end
function decodedblock(enc::Encoding, wrapper::WrapperBlock)
    inner = decodedblock(enc, wrapped(wrapper))
    return isnothing(inner) ? nothing : setwrapped(wrapper, inner)
end
function encode(enc::Encoding, ctx, wrapper::WrapperBlock, data; kwargs...)
    return encode(enc, ctx, wrapped(wrapper), data; kwargs...)
end
function decode(enc::Encoding, ctx, wrapper::WrapperBlock, data; kwargs...)
    return decode(enc, ctx, wrapped(wrapper), data; kwargs...)
end

# ## Named

"""
    Named(name, block)

Wrapper `Block` to attach a name to a block. Can be used in conjunction
with [`Only`](#) to apply encodings to specific blocks only.
"""
struct Named{Name, B<:AbstractBlock} <: WrapperBlock
    block::B
end
Named(name::Symbol, block::B) where {B<:AbstractBlock} = Named{name, B}(block)

# the name is preserved through encodings and decodings
function encodedblock(enc::Encoding, named::Named{Name}) where Name
    outblock = encodedblock(enc, wrapped(named))
    return isnothing(outblock) ? nothing : Named(Name, outblock)
end

function decodedblock(enc::Encoding, named::Named{Name}) where Name
    outblock = decodedblock(enc, wrapped(named))
    return isnothing(outblock) ? nothing : Named(Name, outblock)
end

# ## Many

"""
    Many(block) <: WrapperBlock

`Many` indicates that you can variable number of data instances for
`block`. Consider a bounding box detection task where there may be any
number of targets in an image and this number varies for different
samples. The blocks `(Image{2}(), BoundingBox{2}()` imply that there is exactly
one bounding box for every image, which is not the case. Instead you
would want to use `(Image{2}(), Many(BoundingBox{2}())`.
"""
struct Many{B<:AbstractBlock} <: WrapperBlock
    block::B
end

FastAI.checkblock(many::Many, datas) = all(checkblock(wrapped(many), data) for data in datas)
FastAI.mockblock(many::Many) = [mockblock(wrapped(many)), mockblock(wrapped(many))]

function FastAI.encode(enc::Encoding, ctx, many::Many, datas)
    return map(datas) do data
        encode(enc, ctx, wrapped(many), data)
    end
end

function FastAI.decode(enc::Encoding, ctx, many::Many, datas)
    return map(datas) do data
        decode(enc, ctx, wrapped(many), data)
    end
end

# # Wrapper encodings

"""
    Only(name, encoding)
    Only(block, encoding)

Wrapper that conditionally applies `encoding` only if the block
equals `block` or is a `Named{name}`.
"""
struct Only{Name, E<:Encoding} <: StatefulEncoding
    encoding::E
end

function Only(name::Symbol, encoding::E) where E
    return Only{name, E}(encoding)
end

encodedblock(only::Only{Name}, block::Named{Name}) where Name = encodedblock(only.encoding, block)
decodedblock(only::Only{Name}, block::Named{Name}) where Name = decodedblock(only.encoding, block)

encodestate(only::Only, args...) = encodestate(only.encoding, args...)
decodestate(only::Only, args...) = decodestate(only.encoding, args...)

function encode(
        only::Only{Name},
        context,
        block::Named{Name},
        data;
        state=encodestate(only, context, block, data)) where Name
    return encode(only.encoding, context, block, data; state = state)
end

function decode(
        only::Only{Name},
        context,
        block::Named{Name},
        data;
        state=decodestate(only, context, block, data)) where Name
    return decode(only.encoding, context, block, data; state = state)
end


text:
## AbstractTrees interface

function AbstractTrees.children(tree::BT)
    tree.tasks
end

function AbstractTrees.printnode(io::IO, node::Sequence)
    if node.name != ""
        repr = "$(node.name) ->"
    else
        repr = "->"
    end
    print(io, repr)
end

function AbstractTrees.printnode(io::IO, node::Selector)
    if node.name != ""
        repr = "$(node.name) ?"
    else
        repr = "?"
    end
    print(io, repr)
end