text:
using Koma, Plots, LaTeXStrings
pgfplotsx()
rf_wf = [
    #Waveform am_sg_100_100_0
    0, 413, 859, 1337, 1847, 2389, 2962, 3568, 4204, 4870,
    5565, 6288, 7038, 7814, 8613, 9435, 10277, 11138, 12015, 12906,
    13809, 14721, 15640, 16563, 17487, 18409, 19327, 20237, 21136, 22023,
    22892, 23742, 24570, 25372, 26146, 26889, 27598, 28271, 28905, 29497,
    30046, 30549, 31005, 31411, 31767, 32070, 32320, 32515, 32655, 32739,
    32767, 32739, 32655, 32515, 32320, 32070, 31767, 31411, 31005, 30549,
    30046, 29497, 28905, 28271, 27598, 26889, 26146, 25372, 24570, 23742,
    22892, 22023, 21136, 20237, 19327, 18409, 17487, 16563, 15640, 14721,
    13809, 12906, 12015, 11138, 10277, 9435, 8613, 7814, 7038, 6288,
    5565, 4870, 4204, 3568, 2962, 2389, 1847, 1337, 859, 413,
    0
] ./ (2^16/2 - 1) #signed int 16 bits
Trf = 0.486e-3
Arf = 21e-6 * rf_wf
Gx1 = 12.89e-3
ζ11 = 0.29e-3
ζ21 = 0.37e-3
Tg1 = Trf
delay = Tg1 + ζ11 - Trf
ex1 = Sequence([Grad(0,0); Grad(0,0); Grad(Gx1,Tg1,ζ11,ζ21,0);;],[RF(Arf,Trf,-5000,delay);;],[ADC(0,0)])
ex2 = Sequence([Grad(0,0); Grad(0,0); Grad(Gx1,Tg1,ζ11,ζ21,0);;],[RF(Arf,Trf,    0,delay);;],[ADC(0,0)])
ex3 = Sequence([Grad(0,0); Grad(0,0); Grad(Gx1,Tg1,ζ11,ζ21,0);;],[RF(Arf,Trf, 5000,delay);;],[ADC(0,0)])
Gx2 = -.83e-3
ζ12 = 0.026e-3
ζ22 = 0.28e-3
Tg2 = 5.75e-3
ref = Sequence([Grad(0,0); Grad(0,0); Grad(Gx2,Tg2,ζ12,ζ22,0);;])
#Different slices
seq1 = ex1 + ref
seq2 = ex2 + ref
seq3 = ex3 + ref
#Focus Time
t_focus_aux = (abs(Gx2)*(2Tg2+ζ12+ζ22)-ζ21*Gx1)/(2Gx1)
t_focus = (Tg1+ζ11-t_focus_aux)
#SIM
z = range(-2e-2,2e-2,200)
M1 = simulate_slice_profile(seq1;z)
M2 = simulate_slice_profile(seq2;z)
M3 = simulate_slice_profile(seq3;z)
#Plots
plot_seq(seq2; width=480, height=350, slider=false)
f = γ*Gx1*z*1e-3
plot(f,  abs.(M1.xy),label=L"|M_{xy}| @ \Delta f = -5\, \mathrm{kHz}",
        xlabel=L"\mathrm{kHz}",line=3,title="α=$(Koma.get_flip_angle(seq2.RF[1]))")
plot!(f, abs.(M2.xy),label=L"|M_{xy}| @ \Delta f = +0\, \mathrm{kHz}",line=3)
plot!(f, abs.(M3.xy),label=L"|M_{xy}| @ \Delta f = +5\, \mathrm{kHz}",line=3)

#GIF
# anim = @animate for tf = range(0,7.5e-3,60)
#     simParams = Dict("return_Mag"=>true, "end_sim_at"=>tf, "gpu"=>false)
#     M = simulate(phantom, seq2, sys; simParams)
#     #Plots
#     plot(xx*1e2, abs.(M.xy),label="Mxy",line=4,ylim=(-1,1),title="Simulation ended at $(round(simParams["end_sim_at"]*1e3,digits=3)) ms")
#     plot!(xx*1e2, abs.(M.z),label="Mz",line=:dash)
#     plot!(xx*1e2, real.(M.xy),label="Mx",line=:dash)
#     plot!(xx*1e2, imag.(M.xy),label="My",line=:dash)
# end
# gif(anim, "AnimSliceExc.gif", fps = 10)


text:
"""
    solve(L,Γ;<keywordarguments>)

Find eigenvalues of linear operator family `L` inside contour `Γ`.

# Arguments
- `L::LinearOperatorFamily`: Linear operator family for which eigenvalues are to be found
- `Γ::List`: List of complex points defining the contour.

# Keywordarguments (optional)
- `Δl::Int=1`: Increment to the search space dimension in Beyn's algorithm
- `N::Int=16`: Number of integration points per edge of `Γ`.
- `tol::Float=1E-8`: Tolerance for local eigenvalue solver.
- `eigvals::Dict=Dict()`: Dictionary of known eigenvalues and corresponding eigenvectors.
- `maxcycles::Int=1`: Number of local correction cycles to global eigenvalue estimate.
- `nev::Int=1`: Number of tested eigenvalues for finding the best local match to the global estimates.
- `max_outer_cycles::Int=1`: Number of outer cycles, i.e., incremental increases to the search space dimension by `Δl`.
- `atol_σ::Float=1E-12`: absolute tolerance for terminating the inner cycle if the maximum singular value is `σ<atol_σ``
- `rtol_σ::Float=1E-8`: relative tolerance for terminating the inner cycle if the the ratio of the initial and the current maximum singular value is `σ/σ0<rtol_σ`
- `loglevel::Int=0`: loglevel value between 0 (no output) and 2 (max output) to control the detail of printed output.

# Returns
- `eigvals::Dict`:List of computed eigenvalues and corresponding eigenvectors.

# Notes
The algorithm is experimental. It startes with a low-dimensional Beyn integral.
From which estimates to the eigenvalues are computed. These estimates are used
to analytically correct Beyn's integral. This step requires no additional
numerical integration and is therefore relatively quick. New local estimates are
computed and if the local solver then found new eigenvalues these are used to
further correct Beyn's integral. In an optional outer loop the entire procedure
can be repeated, after increasing the search space dimension.

See also: [`householder`](@ref), [`beyn`](@ref)
"""
function solve(L::LinearOperatorFamily,Γ;Δl=1,N=16, tol=1E-8,eigvals=Dict(),maxcycles=1,nev=1,max_outer_cycles=1,atol_σ=1E-12,rtol_σ=1E-8,loglevel=0)
    #header
    d=size(L.terms[1].coeff)[1]#system dimension
    A=[]
    l=Δl
    σ0,σ,σmax=0.0,0.0,0.0
    #outer loop
    while l<=max_outer_cycles*Δl
        V=zeros(ComplexF64,d,Δl)
        for ll=1:Δl
            #TODO: randomization
            V[(l-Δl)+ll,ll]=1.0+0.0im
        end
        A=append!(A,[compute_moment_matrices(L,Γ,V,K=1,N=N,output=true),])
        #σmax=max(σmax,maximum(abs.(A[1][:,:,1])))
        if l>Δl #reinitialize σmax in new iteration of outer loop
            Ω,P,Σ=moments2eigs(A,return_σ=true)
            #println(maximum(Σ))
            σmax,σ0,σ=max(σmax,maximum(Σ)),maximum(Σ),0.0
        end
        #update moment matrix with known solutions
        for (ω, val) in eigvals
            w=wn(ω,Γ)
            for ll=1:Δl
                moment=-2*pi*1.0im*w*val[1].v*val[1].v_adj[ll]'
                A[l÷Δl][:,ll,1]+=moment
                A[l÷Δl][:,ll,2]+=ω*moment
            end
        end

        number_of_interior_eigvals=0
        for val in values(eigvals)
            if val[2]
                number_of_interior_eigvals+=1
            end
        end
        neigval=number_of_interior_eigvals
        cycle=0
        while cycle<maxcycles# inner cycle
            cycle+=1
            if loglevel>=1
                flush(stdout)
                println("####cycle$cycle####")
            end
            Ω,P,Σ=moments2eigs(A,return_σ=true)
            #println(maximum(Σ))
            σmax,σ0,σ=max(σmax,σ),σ,maximum(Σ)

            #compute accurate local solutions
            for idx=1:length(Ω)
                ω=Ω[idx]
                if loglevel>=2
                    flush(stdout)
                    println("omega:$ω")
                end
                #remove known vectors from first Krylov vector
                 v0=P[:,idx]
                v0/=sqrt(v0'v0)
                for (key,val) in eigvals
                    v=val[1].v
                    v/=sqrt(v'v)
                    h=v'v0
                    v0-=h*v
                    v0/=sqrt(v0'v0)
                end
                # #TODO: initilaize adjoint vector
                #
                #
                # sol,nn,flag=householder(L,ω,maxiter=10,tol=tol,output=false,order=3,nev=nev,v0=v0)
                sol,nn,flag=mehrmann(L,ω,maxiter=10,tol=tol,output=false,x0=v0,v=v0)
                ω=sol.params[sol.eigval]
                if flag>=0
                    is_new_value=true
                    for val in keys(eigvals)
                        if abs(ω-val)<10*tol
                            is_new_value=false
                            #break
                        end
                    end
                else
                    is_new_value=false
                end
                if loglevel>=2
                    println("conv:$ω flag:$flag new: $is_new_value")
                end
                #update moment matrix with local solutions
                if is_new_value && inpoly(ω,Γ)
                    #R=5*tol
                    #dφ=2pi/(3+1)
                    #φ=0:dφ:2pi-dφ
                    w=wn(ω,Γ)
                    #C=ω.+R*exp.(-1im.*φ*w) #this circle is winded opposite to Γ
                    #A+=compute_moment_matrices(L,C;l=l,N=N,output=true,random=false)
                    #moment2=compute_moment_matrices(L,C;l=l,N=N,output=true,random=false)
                    for idx=1:length(A)
                        for ll=1:Δl
                            moment=-2*pi*1.0im*w*sol.v*sol.v_adj[(idx-1)*Δl+ll]'
                            A[idx][:,ll,1]+=moment
                            A[idx][:,ll,2]+=ω*moment
                        end
                    end
                    eigvals[ω]=[sol,true]
                elseif is_new_value
                    #handle outside values
                    eigvals[ω]=[sol,false]
                end
            end

            number_of_interior_eigvals=0
            for val in values(eigvals)
                if val[2]
                    number_of_interior_eigvals+=1
                end
            end

            if neigval==number_of_interior_eigvals
                break
            else
                neigval=number_of_interior_eigvals
            end
        end
        # σ0,σ=σ,0.0
        # for idx =1:length(A)
        #     σ=max(σ,maximum(abs.(A[idx][:,:,1])))
        # end
        # σmax=max(σmax,σ)
        if loglevel>=1
            flush(stdout)
            println("maximum σmax:$σmax")
            println("final σ:$σ")
            println("relative σ/σ0:$(σ/σ0)")
            println("relative σ/σmax:$(σ/σmax)")
            println("cycles:$cycle,maxcycles:$(maxcycles<=cycle)")
        end
        #stopping criteria
        if σ/σmax<rtol_σ || σ<atol_σ
            if loglevel>=1
                println("finished after $l integrations")
            end
            break
        else
            l+=Δl
            println("you should run more outer cycles!")
        end

    end #of outer loop
    return eigvals
end

ext: jl


text:
using Distributed
using Workspace
using Test
using Configurations
using StructTypes
using JSON

d = Dict(
    "name" => "test_workspace",
    "uuid" => "94f78cb8-c85a-11eb-0dc6-f7f6bd9facf1",
    "compiler" => Dict(
        "project" => Dict(
            "deps" => Dict(
                "JSON" => "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
            ),
            "compat" => Dict(
                "JSON" => "0.20.0",
            )
        )
    )
)

s = JSON.json(d)
ins = from_dict(WorkspaceInstance, JSON.parse(s))
table = WorkspaceTable()
current_nproc = nprocs()
add_workspace_process!(table, ins)
eval_in_workspace(table, ins, quote
    println("aaaaa")
end)

@test nprocs() == current_nproc + 1
rm_workspace!(table, ins)
@test nprocs() == current_nproc

f = Distributed.remotecall(2) do ch
end
fetch(f)


text:
immutable PointwiseLayerState{P<:AbstractPointwise} <:
        AbstractScatteredLayerState
    blobs::Vector{Mocha.Blob}
    layer::PointwiseLayer{P}
end

function PointwiseLayerState(
        backend::Mocha.CPUBackend,
        layer::PointwiseLayer,
        inputs::Vector{Mocha.Blob})
    blobs = Vector{Mocha.Blob}(length(inputs))
    for idblob in eachindex(inputs)
        blobs[idblob] = layer.ρ(inputs[idblob])
    end
    return PointwiseLayerState(blobs, layer)
end

function Mocha.setup(
        backend::Mocha.CPUBackend,
        layer::PointwiseLayer,
        inputs::Vector{Mocha.Blob},
        diffs::Vector{Mocha.Blob})
    return PointwiseLayerState(backend, layer, inputs)
end

function forward(
        backend::Mocha.CPUBackend,
        layerstate::PointwiseLayerState,
        inputs::Vector{Mocha.Blob})
    for id in eachindex(inputs)
        map!(layerstate.ρ, layerstate.blobs[id], inputs[id])
    end
end

function forward!(
        backend::Mocha.CPUBackend,
        state::PointwiseLayerState,
        ρ::AbstractPointwise,
        inputs::Vector)
    @inbounds for idblob in eachindex(inputs)
        map!(ρ, state.blobs[idblob], inputs[idblob])
    end
end


text:
macro left_associative(f::Symbol, x::Symbol, y::Symbol, z::Symbol)
    return quote
        $(esc(f))($(esc(f))($(esc(x)), $(esc(y))), $(esc(z))),
        $(esc(f))($(esc(x)), $(esc(f))($(esc(y)), $(esc(z))))
    end
end

macro right_associative(f::Symbol, x::Symbol, y::Symbol, z::Symbol)
    return quote
        $(esc(f))($(esc(x)), $(esc(f))($(esc(y)), $(esc(z)))),
        $(esc(f))($(esc(f))($(esc(x)), $(esc(y))), $(esc(z)))
    end
end

"""
f is commutative if f(x, y) = f(y, x)
"""
macro commutative(f::Symbol, x::Symbol, y::Symbol)
    return quote
        $(esc(f))($(esc(x)), $(esc(y))),
        $(esc(f))($(esc(y)), $(esc(x)))
    end
end

"""
A function f is left distributive over g if f(x, g(y, z)) = g(f(x, y), f(x, z))
"""
macro distributive(f::Symbol, g::Symbol, x::Symbol, y::Symbol, z::Symbol)
    return quote
        $(esc(f))($(esc(x)), $(esc(g))($(esc(y)), $(esc(z)))),
        $(esc(g))($(esc(f))($(esc(x)), $(esc(y))), $(esc(f))($(esc(x)), $(esc(z))))
    end
end

"""
A function f is left factorisable over g iff f(g(x, y), g(x, z)) = g(x, f(y, z))
"""
macro factorisable(f::Symbol, g::Symbol, x::Symbol, y::Symbol, z::Symbol)
    return quote
        $(esc(f))($(esc(g))($(esc(x)), $(esc(y))), $(esc(g))($(esc(x)), $(esc(z)))),
        $(esc(g))($(esc(x)), $(esc(f))($(esc(y)), $(esc(z))))
    end
end

"""
f is idempotent if f(x, x) = x
"""
macro idempotent(f::Symbol, x::Symbol)
    return quote
        $(esc(f))($(esc(x)), $(esc(x))), $(esc(x))
    end
end

"""
A function f is invloute if f(f(x)) = x
"""
macro involute(f::Symbol, x::Symbol)
    return quote
        $(esc(f))($(esc(f))($(esc(x)))), $(esc(x))
    end
end

macro re_involute(f::Symbol, x::Symbol)
    return quote
        $(esc(x)), $(esc(f))($(esc(f))($(esc(x))))
    end
end


text:
export dimtree, children, count_leaves, dims, display, height, isequal, ==, is_leaf, is_left, is_right
export left_child_length, lvl, nodes_on_lvl, node2ind, non, parent, positions, show, sibling, structure, subnodes, subtree

"""
    dimtree(leaves::AbstractVector{<:Integer}[,internal_nodes::AbstractVector{<:Integer}])
    dimtree(N[,treetype])

Dimension tree. Create from:
- a vector of leaves (and a vector of internal nodes),
- an order of a tensor N and a type treetype of a dimtree. Default: treetype="balanced".
"""
mutable struct dimtree
  leaves::AbstractVector{<:Integer}
  internal_nodes::AbstractVector{<:Integer}
  function dimtree(leaves::AbstractVector{<:Integer},internal_nodes::AbstractVector{<:Integer})
    @assert(length(leaves)==length(internal_nodes)+1,"Incorrect nodes.")
    nodes=sort([leaves;internal_nodes])
    @assert(nodes==collect(1:maximum(nodes)),"Missing nodes.")
    new(leaves,internal_nodes)
  end
end
#dimtree(leaves::AbstractVector{<:Integer},internal_nodes::AbstractVector{<:Integer})=dimtree{D1,D2}(leaves,internal_nodes)
function dimtree(leaves::AbstractVector{<:Integer})
  N=length(leaves)
  nn=2*N-1
  internal_nodes=setdiff(collect(1:nn),leaves)
  @assert(sort([leaves;internal_nodes])==collect(1:nn),"Wrong input.")
  dimtree(leaves,internal_nodes)
end

function dimtree(N::Integer,treetype="balanced")
  @assert(treetype=="balanced","Only balanced trees are supported.")
  C=zeros(Int,2*N-1,2)
  dims=Array{Any}(undef,2*N-1)
  dims[1]=collect(1:N)
  nn=1
  i=1
  while i<=nn
    if length(dims[i]) == 1
      C[i,:]=[0 0]
    else
      ind_left=nn+1
      ind_right=nn+2
      nn+=2
      C[i,:]=[ind_left,ind_right]
      dims[ind_left]=dims[i][1:ceil(Int,end/2)]
      dims[ind_right]=dims[i][ceil(Int,end/2)+1:end]
    end
    i+=1
  end
  L=findall((in)([0,0]),C[:,1])
  leaves=zeros(Int,N)
  [leaves[dims[l]].=l for l in L]
  dimtree(leaves)
end

"""
    children(T[,node])

Matrix of children for each node of a dimtree T.
"""
function children(T::dimtree)
  l=length(T.leaves)+length(T.internal_nodes)
  C=zeros(Int,l,2)
  C[1,:]=[2 3]
  next_node=4
  for n=2:l
    if n in T.internal_nodes
      C[n,:]=[next_node next_node+1]
      next_node+=2
    end
  end
  C
end
function children(T::dimtree,node::Integer)
    children(T)[node,:]
end

"""
    count_leaves(T,node)

Number of leaves under a node of a dimtree T.
"""
function count_leaves(T::dimtree,node::Int)
  @assert(0<node<=non(T),"Input exceeds number of nodes.")
  nodes=subnodes(T,node)
  l=0
  for n in nodes
    if is_leaf(T,n)
      l+=1
    end
  end
  l
end

"""
    dims(T,node)

Content of a node of a dimtree T.
"""
function dims(T::dimtree)
    dims_mat=Array{Vector}(undef,non(T))
    [dims_mat[T.leaves[n]]=[n] for n=1:length(T.leaves)]
    for n in reverse(T.internal_nodes)
        c=children(T,n)
        dims_mat[n]=[dims_mat[c[1]]...,dims_mat[c[2]]...]
    end
    dims_mat
end
dims(T::dimtree,node::Integer)=dims(T)[node]

"""
    display(T)

Display a dimtree T.
"""
function display(T::dimtree)
  println("Dimensional tree:" )
  if T.leaves==[1]
    println("   1   ")
    return
  end
  t,tl,tr=structure(T)
  h=height(T)
  initial_blank_len=zeros(Int,h)
  a=2
  for l=2:h
    initial_blank_len[l]=initial_blank_len[l-1]+a
    a=a*2
  end
  blank_len=2*initial_blank_len[1:end-1].+1
  l=h;k=h-1;
  [print(" ") for i=1:initial_blank_len[l]]; l-=1
  print(minimum(t[1]),"-",maximum(t[1]))
  println()
  [print(" ") for i=1:initial_blank_len[l]]; l-=1
  length(tl[1]) == 1 ? print(" ",minimum(tl[1])," ") : print(minimum(tl[1]),"-",maximum(tl[1]))
  [print(" ") for i=1:blank_len[k]]; k-=1
  length(tr[1]) == 1 ? print(" ",minimum(tr[1])," ") : print(minimum(tr[1]),"-",maximum(tr[1]))
  println()
  if l>0
  [print(" ") for i=1:initial_blank_len[l]]; l-=1
  m=2
  nodes=collect(1:non(T)) #nodes
  nodepos=positions(T) #node position
  nnmbr=4 #node counter
  level_old=2
  bingo=0
  for p=4:2^h-1  #loop over all positions in a binary tree
    if bingo == 1 || nnmbr>non(T)
      bingo=0
      continue
    end
    level=lvl(T,nodes[nnmbr])
     if level!=level_old && p==2^level
       println()
       [print(" ") for i=1:initial_blank_len[l]];  l-=1
       level_old=level;
       level+=1
       k-=1
     end
    #println("p = $p, nodepos[$nnmbr] = $(nodepos[nnmbr])")
    if nodepos[nnmbr]==p
      length(tl[m]) == 1 ? print(" ",tl[m][1]," ") : print(tl[m][1],"-",maximum(tl[m]))
      [print(" ") for i=1:blank_len[k]]
      length(tr[m]) == 1 ? print(" ",tr[m][1]," ") : print(tr[m][1],"-",maximum(tr[m]))
      [print(" ") for i=1:blank_len[k]]
      m+=1
      nnmbr+=2
      bingo=1
    else
      #print(" ")
      [print(" ") for i=1:3+blank_len[k]]
      #print("")
    end
  end
  end
  println()
end

"""
    height(T)

Height of a dimtree T.
"""
function height(T::dimtree)
    maximum(lvl(T))+1
end

"""
    isequal(T1,T2)

True if two dimensional trees are equal, false otherwise.
"""
function isequal(T1::dimtree,T2::dimtree)
  if T1.leaves==T2.leaves
    return true
  else
    return false
  end
end
==(T1::dimtree,T2::dimtree)=isequal(T1,T2)

"""
    is_leaf(T,node)

True if a node is a leaf in a dimtree T. False otherwise.
"""
function is_leaf(T::dimtree,node::Integer)
   node in T.leaves ? true : false
end

"""
    is_left(T,node)

True if a node is a left node of dimtree T. False otherwise.
"""
function is_left(T::dimtree,node::Integer)
    ind=findall(children(T)[:,1].==node)
    if length(ind)>0
        return true
    else
        return false
    end
end

"""
    is_right(T,node)

True if a node is a right node of dimtree T. False otherwise.
"""
function is_right(T::dimtree,node::Integer)
    ind=findall(children(T)[:,2].==node)
    if length(ind)>0
        return true
    else
        return false
    end
end

"""
    left_child_length(T[,node])

Length of the left child of a node of a dimtree T. Default: node=1.
"""
function left_child_length(T::dimtree,node::Int)
    if node in T.leaves
        return 0
    else
        return count_leaves(T,children(T)[node,1])
    end
end
left_child_length(T::dimtree)=left_child_length(T,1)

"""
    lvl(T[,node])

Vector of levels for each node of a dimtree T.
"""
function lvl(T::dimtree)
    nn=non(T)
    L=zeros(Int,nn)
    for n in collect(2:nn)
        L[n]=L[parent(T,n)]+1
    end
    L
end
lvl(T::dimtree,node::Integer)=lvl(T)[node]
function lvl(T::dimtree,nodes::AbstractVector{<:Integer})
    [lvl(T)[n] for n in nodes]
end

"""
    nodes_on_lvl(T,l)

Nodes on a level l in a dimtree T.
"""
function nodes_on_lvl(T::dimtree,l::Integer)
    findall(lvl(T).==l)
end

"""
    node2ind(T,nodes)

Convert node numbers to transfer tensor or frames indices in a dimtree T.
"""
function node2ind(T::dimtree,nodes::AbstractVector{<:Integer})
    ind=zeros(Int,length(nodes))
    k=1
    for n in nodes
        if is_leaf(T,n)
            ind[k]=findall((in)(n),T.leaves)[1]
        else
            ind[k]=findall((in)(n),T.internal_nodes)[1]
        end
        k+=1
    end
    ind
end
node2ind(T::dimtree,nodes::Integer)=node2ind(T,[nodes])[1]
node2ind(T::dimtree)=node2ind(T,collect(1:non(T)))

"""
    non(T)

Number of nodes of a dimtree T.
"""
function non(T::dimtree)
  maximum([T.leaves;T.internal_nodes])
end

"""
    parent(T[,node])

Vectors of parents for each node of a dimtree T.
"""
function parent(T::dimtree)
  nn=non(T)
  C=children(T)
  P=zeros(Int,nn) #parent[node]
  P[1]=0
  for n=2:nn
    P[n]=CartesianIndices(C)[findall(C.==n)][1][1]
  end
  P
end
parent(T::dimtree,node::Integer)=parent(T)[node]
parent(T::dimtree,nodes::AbstractVector{<:Integer})=[parent(T)[n] for n in nodes]

"""
    positions(T)

Positions of nodes of a dimtree T in a full binary tree.
"""
function positions(T::dimtree)
    h=height(T)
    pos_nmbr=2^h-1
    position=zeros(Int,non(T))
    position[1]=1;
    i=2;
    maxit=100;
    for l=1:h
        nds=nodes_on_lvl(T,l)
        if length(nds)==2^l
            for n in nds
                position[n]=i
                i+=1
            end
        else
            for n in nds
                if position[n]!=0
                    continue
                end
                p=parent(T,n)
                while i!=position[p]*2 && i<maxit
                    i+=2
                end
                position[n]=i
                i+=1
                position[n+1]=i
                i+=1
                n+=1
            end
        end
    end
    position
end

function Base.show(T::dimtree)
    display(T)
end

"""
    sibling(T,node)

Sibling of a node in a dimtree T.
"""
function sibling(T::dimtree,node::Integer)
  C=children(T)
  ind=findall(C.==node)[1]
  if ind[2]==1
      C[ind[1]...,2]
  else
      C[ind[1]...,1]
  end
end

"""
    structure(T)

For each internal node of a given dimtree, returns its representation t, its left child and right child representations tl and tr.
"""
function structure(T::dimtree)
  t=dims(T)[T.internal_nodes]
  tl=dims(T)[children(T)[T.internal_nodes,1]]
  tr=dims(T)[children(T)[T.internal_nodes,2]]
  t,tl,tr
end

"""
    subnodes(T,node)

Nodes in the subtree of a node in a dimtree T.
"""
function subnodes(T::dimtree,node::Integer)
    sub_ind=[node]
    n=1
    for i=node:non(T)
        if !is_leaf(T,sub_ind[n])
            sub_ind=[sub_ind; children(T,sub_ind[n])]
        end
        if n!=length(sub_ind)
            n+=1
        end
    end
    sub_ind
end

"""
    subtree(T,node)

Dimensional tree with the same structure as a subtree of a node in a dimtree T. See also: subnodes.
"""
function subtree(T::dimtree,node::Integer)
  if is_leaf(T,node)
    return dimtree([1],Int[])
  end
  nodes=subnodes(T,node)
  #L=T.leaves[findin(T.leaves,nodes)]
  L=T.leaves[findall((in)(nodes),T.leaves)]
  #I=T.internal_nodes[findin(T.internal_nodes,nodes)]
  I=T.internal_nodes[findall((in)(nodes),T.internal_nodes)]
  dimtree(sortperm(L).+length(I))
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using Base.Iterators: partition
using Random
using BSON
using Statistics
using MLDatasets
using Flux
using Flux: onehotbatch, onecold, crossentropy
using Flux.Data: DataLoader
using Plots
using ImageInspector
import RDatasets: dataset

plot(1:10)

iris = dataset("datasets", "iris")

T = Float32
X_train, y_train = MLDatasets.MNIST.traindata(T)
X_train = reshape(X_train, size(X_train,1), size(X_train,2), 1, size(X_train,3))
y_train = onehotbatch(y_train, 0:9)
MLDatasets.MNIST.testdata(T)
MLDatasets.CIFAR10.traindata(T)

imageplot(1 .- X_train, 1:3; nrows=1, size=(800,480))

m_aux = Chain(
    Conv((2,2), 1=>16, relu),
    MaxPool((2,2)),
    Conv((2,2), 16=>8, relu),
    MaxPool((2,2)),
    flatten,
    Dense(288, size(y_train,1)),
    softmax,
)

L_aux(X, y) = crossentropy(m_aux(X), y)

batches_aux = DataLoader((X_train, y_train); batchsize=64, shuffle = true)

gradient(() -> L_aux(X_train[:,:,:,1:10], y_train[:,1:10]), params(m_aux))

onecold(m_aux(X_train[:,:,:,1:10]))

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
using NatureGas
using Test
#=
gas=[
"CO2"    0.006  0.005  0.015  0.016  0.076  0.011;
"N2"     0.003  0.031  0.010  0.100  0.057  0.117;
"H2"     0.00   0.00   0.00   0.095  0.00   0.00;
"CO"     0.00   0.00   0.00   0.010  0.00   0.00;
"CH4"    0.965  0.907  0.859  0.735  0.812  0.826;
"C2H6"   0.018  0.0450 0.085  0.033  0.043  0.035;
"C3H8"   0.0045 0.0084 0.023  0.0074 0.009  0.0075;
"iC4H10" 0.0010 0.0010 0.0035 0.0012 0.0015 0.0012;
"nC4H10" 0.0010 0.0015 0.0035 0.0012 0.0015 0.0012;
"iC5H12" 0.0005 0.0003 0.0005 0.0004 0.00   0.0004;
"nC5H12" 0.0003 0.0004 0.0005 0.0004 0.00   0.0004;
"nC6H14" 0.0007 0.0004 0.00   0.0002 0.00   0.0002;
"nC7H16" 0.00   0.00   0.00   0.0001 0.00   0.000 ;
"nC8H18" 0.00   0.00   0.00   0.0001 0.00   0.00]
=#
case=[
60  -3.15 0.84053 0.83348 0.79380 0.88550 0.82609 0.85380;
60   6.85 0.86199 0.85596 0.82206 0.90144 0.84969 0.87370;
60  16.85 0.88006 0.87484 0.84544 0.91501 0.86944 0.89052;
60  36.85 0.90867 0.90466 0.88183 0.93674 0.90052 0.91723;
60  56.85 0.93011 0.92696 0.90868 0.95318 0.92368 0.93730;
120 -3.15 0.72133 0.71044 0.64145 0.81024 0.69540 0.75074;
120  6.85 0.76025 0.75066 0.68971 0.83782 0.73780 0.78586;
120 16.85 0.79317 0.78475 0.73123 0.86137 0.77369 0.81569;
120 36.85 0.84515 0.83863 0.79697 0.89913 0.83022 0.86311;
120 56.85 0.88383 0.87870 0.84553 0.92766 0.87211 0.89862]
atol=0.005
@testset "Gas 1" begin
    for i=1:10
        @test isapprox(ng_zfactor(case[i,1]/10,case[i,2]+273.15,CO2=0.006,N2=0.003,CH4=0.965,C2H6=0.018,
        C3H8=0.0045,iC4H10=0.001,nC4H10=0.001,iC5H12=0.0005,nC5H12=0.0003,nC6H14=0.0007),case[i,3];rtol=atol)
    end
end

@testset "Gas 2" begin
    for i=1:10
        @test isapprox(ng_zfactor(case[i,1]/10,case[i,2]+273.15,CO2=0.005,N2=0.031,CH4=0.907,C2H6=0.045,
        C3H8=0.0084,iC4H10=0.001,nC4H10=0.0015,iC5H12=0.0003,nC5H12=0.0004,nC6H14=0.0004),case[i,4];rtol=atol)
    end
end

@testset "Gas 3" begin
    for i=1:10
        @test isapprox(ng_zfactor(case[i,1]/10,case[i,2]+273.15,CO2=0.015,N2=0.01,CH4=0.859,C2H6=0.085,
        C3H8=0.023,iC4H10=0.0035,nC4H10=0.0035,iC5H12=0.0005,nC5H12=0.0005,nC6H14=0.0000),case[i,5];rtol=atol)
    end
end

@testset "Gas 4" begin
    for i=1:10
        @test isapprox(ng_zfactor(case[i,1]/10,case[i,2]+273.15,CO2=0.016,N2=0.1,H2=0.095,CO=0.01,CH4=0.735,C2H6=0.033,
        C3H8=0.0074,iC4H10=0.0012,nC4H10=0.0012,iC5H12=0.0004,nC5H12=0.0004,nC6H14=0.0002,nC7H16=0.0001,nC8H18=0.0001),case[i,6];rtol=atol)
    end
end

@testset "Gas 5" begin
    for i=1:10
        @test_skip isapprox(ng_zfactor(case[i,1]/10,case[i,2]+273.15,CO2=0.076,N2=0.057,CH4=0.812,C2H6=0.043,
        C3H8=0.009,iC4H10=0.0015,nC4H10=0.0015),case[i,7];rtol=atol)
    end
end

@testset "Gas 6" begin
    for i=1:10
             @test isapprox(ng_zfactor(case[i,1]/10,case[i,2]+273.15,CO2=0.011,N2=0.117,CH4=0.826,C2H6=0.035,
            C3H8=0.0075,iC4H10=0.0012,nC4H10=0.0012,iC5H12=0.0004,nC5H12=0.0004,nC6H14=0.0002,nC7H16=0.0001),case[i,8];rtol=atol)
    end
end

max_stars_repo_name: hzgzh/NG

max_issues_repo_name: hzgzh/NG

max_forks_repo_name: hzgzh/NG


text:
function findneighbor(square, ci, dir)
    if dir == :L
        res = findlast(x -> x > 0, square[1:ci[1] - 1, ci[2]])
        res !== nothing && return CartesianIndex(res, ci[2])
    elseif dir == :R
        res = findfirst(x -> x > 0, square[ci[1] + 1:end, ci[2]])
        res !== nothing && return CartesianIndex(ci[1] + res, ci[2])
    elseif dir == :D
        res = findlast(x -> x > 0, square[ci[1], 1:ci[2] - 1])
        res !== nothing && return CartesianIndex(ci[1], res)
    elseif dir == :U
        res = findfirst(x -> x > 0, square[ci[1], ci[2] + 1:end])
        res !== nothing && return CartesianIndex(ci[1], ci[2] + res)
    end

    return nothing
end

const dirs = (:L, :R, :D, :U)

function process!(case, square)
    total_interest = 0

    while true
        next_square = copy(square)
        total_interest += sum(square)
        tocheck = findall(x -> x > 0, square)
        length(tocheck) > 0 || break

        for ci in tocheck
            neighbor_skill = 0
            neighbors = 0

            for dir in dirs
                cn = findneighbor(square, ci, dir)

                if cn !== nothing
                    neighbors += 1
                    neighbor_skill += square[cn]
                end
            end

            if neighbors > 0 && neighbor_skill / neighbors > square[ci]
                next_square[ci] = 0
            end
        end

        square == next_square && break
        square = next_square
    end

    println("Case #$(case): $(total_interest)")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        R, C = parse.(Int, stdin |> readline |> split)
        square = Matrix{Int}(undef, R, C)
        for i in 1:R
            square[i, :] = parse.(Int, stdin |> readline |> split)
        end

        process!(case, square)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
@testset "Simulation" begin
    @testset "LogNormal" begin
        Random.seed!(123)
        y = simulate_GAS_1_12(LogNormal, 0.0, 123)

        v = [0.1, 0.1]
        gas = Model([1, 12], [1, 12], LogNormal, 0.0)
        gas.ω = v
        gas.A[1]  = convert(Matrix{Float64}, Diagonal(3*v))
        gas.A[12] = convert(Matrix{Float64}, Diagonal(-3*v))
        gas.B[1]  = convert(Matrix{Float64}, Diagonal(3*v))
        gas.B[12] = convert(Matrix{Float64}, Diagonal(-3*v))

        scenarios_observations, scenarios_params = simulate(y, gas, 50, 1000)
        @test maximum(scenarios_observations) < 1e5
        @test minimum(scenarios_observations) > 0

        forec = forecast(y, gas, 50)
        @test all([forec.observation_quantiles[i, 1] .<
                   forec.observation_quantiles[i, 2] .<
                   forec.observation_quantiles[i, 3] for i in 1:50])
        @test all([forec.parameter_quantiles[i, 1, 1] <=
                   forec.parameter_quantiles[i, 1, 2] <=
                   forec.parameter_quantiles[i, 1, 3] for i in 1:50])
        @test all([forec.parameter_quantiles[i, 2, 1] <=
                   forec.parameter_quantiles[i, 2, 2] <=
                   forec.parameter_quantiles[i, 2, 3] for i in 1:50])
        @test maximum(forec.observation_quantiles) < 1e3
        @test minimum(forec.observation_quantiles) > 1e-4

        @test minimum(forec.parameter_scenarios[:, 2, :]) > 0.0
    end
end


text:
function get_celestial_bodies(obj::DataContext)
    return jcall(obj, "getCelestialBodies", CelestialBodies, ())
end

function get_default(::Type{DataContext})
    return jcall(DataContext, "getDefault", LazyLoadedDataContext, ())
end

function get_frames(obj::DataContext)
    return jcall(obj, "getFrames", Frames, ())
end

function get_geo_magnetic_fields(obj::DataContext)
    return jcall(obj, "getGeoMagneticFields", GeoMagneticFields, ())
end

function get_gravity_fields(obj::DataContext)
    return jcall(obj, "getGravityFields", GravityFields, ())
end

function get_time_scales(obj::DataContext)
    return jcall(obj, "getTimeScales", TimeScales, ())
end

function set_default(::Type{DataContext}, arg0::LazyLoadedDataContext)
    return jcall(DataContext, "setDefault", void, (LazyLoadedDataContext,), arg0)
end


text:
@with_kw struct Problem{T <: AbstractFloat}
    A::SparseMatrixCSC{T, Int64}
    G::SparseMatrixCSC{T, Int64}
    c::Vector{T}
    b::Vector{T}
    h::Vector{T}
    C::CompositeCone{T}

    @assert all(isfinite.(A.nzval))
    @assert all(isfinite.(G.nzval))
    @assert all(isfinite.(b))
    @assert all(isfinite.(c))
    @assert all(isfinite.(h))
    @assert size(A) == (length(b), length(c))
    @assert size(G) == (length(h), length(c))
    @assert dim(C) == length(h)
end

Problem(A, G, c, b, h, C::NonnegCone) = Problem(A, G, c, b, h, CompositeCone(C))
Problem(A, G, c, b, h, C::QuadCone) = Problem(A, G, c, b, h, CompositeCone(C))

get_num_vars(prob::Problem) = length(prob.c)
get_cone_size(prob::Problem) = length(prob.h)
get_num_constrs(prob::Problem) = length(prob.b)
get_prob_sizes(prob::Problem) = (get_num_vars(prob), get_num_constrs(prob), get_cone_size(prob))

function LinearProgram(A::AbstractMatrix{T}, G::AbstractMatrix{T}, c::Vector{T},
                       b::Vector{T}, h::Vector{T}) where T <: AbstractFloat
    return Problem(A, G, c, b, h, NonnegCone{T}(length(h)))
end

function QuadraticProgram(A::SparseMatrixCSC{T, Int64}, G::SparseMatrixCSC{T, Int64},
                          P::SparseMatrixCSC{T, Int64}, q::Vector{T}, b::Vector{T},
                          h::Vector{T}) where T <: AbstractFloat
    m, n = size(A)
    p, n = size(G)
    idx = diag(P) .> 0
    Pe = P[idx, idx]
    ne = sum(idx)
    ldl_fact = QDLDL.qdldl(Pe; reg=1e-14, is_psd=true)
    L = (ldl_fact.L + sparse(Diagonal(ones(ne)))) * sqrt.(ldl_fact.D)
    D = ldl_fact.D
    perm = ldl_fact.P
    Ue = sparse(L')[:, sortperm(perm)] / sqrt(2)
    U = spzeros(ne, n)
    U[:, idx] = Ue

    At = [spzeros(m) A]
    bt = b
    ct = [1; q]
    Gt = [spzeros(p) G ; -1 spzeros(1,n); spzeros(ne) -2*U; 1 spzeros(1,n)]
    ht = [h; 1; zeros(ne); 1]
    if p > 0
        Ct = CompositeCone([NonnegCone{T}(p), QuadCone{T}(ne+2)])
    else
        Ct = CompositeCone([QuadCone{T}(ne+2)])
    end
    return Problem(At, Gt, ct, bt, ht, Ct)
end


text:
println("Testing NaN detection.")

# Creating a NaN via computation.
tree = cos(exp(exp(exp(exp(Node("x1"))))))
X = randn(MersenneTwister(0), Float32, 1, 100) * 100.0f0
output, flag = eval_tree_array(tree, X, options)
@test !flag

# Creating a NaN/Inf via division by constant zero.
tree = cos(Node("x1") / 0.0f0)
output, flag = eval_tree_array(tree, X, options)
@test !flag

# Having a NaN/Inf constants:
tree = cos(Node("x1") + Inf)
output, flag = eval_tree_array(tree, X, options)
@test !flag
tree = cos(Node("x1") + NaN)
output, flag = eval_tree_array(tree, X, options)
@test !flag

println("Passed.")


text:
include("../Tests/Tests.jl")

using Test
using JuMP
using GLPK
optimizer = GLPK.Optimizer()
bridged = MOI.Bridges.full_bridge_optimizer(optimizer, Float64)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5, query=false)
exclude = ["plus_minus"]
Tests.linear_test(bridged, config, exclude)


text:
const NanosPerNanosecond  = 1
const NanosPerMicrosecond = 1000
const NanosPerMillisecond = 1000*1000
const NanosPerSecond      = 1000*1000*1000
const NanosPerMinute      = 1000*1000*1000*60
const NanosPerHour        = 1000*1000*1000*60*60

const NanosPer = (NanosPerHour, NanosPerMinute, NanosPerSecond, NanosPerMillisecond, NanosPerMicrosecond, NanosPerNanosecond)
(3600000000000, 60000000000, 1000000000, 1000000, 1000, 1)

const TimeInParts = NamedTuple{(:hours, :minutes, :seconds, :millis, :micros, :nanos)}

struct TimeParts <: Function
end

TimeParts(x::Vector) = TimeInParts(x...)

function timeparts(nanos::Int128, parts::Vector{Int128})
    ns = nanos

    for (idx, unit_in_nanos) in enumerate(NanosPer)
        units = fld(ns, unit_in_nanos)
        ns -= units * unit_in_nanos
        parts[idx] = units
    end

    parts
end

function timeparts(nanos::Int128)
    ns = nanos
    result = zeros(Int128, length(NanosPer))

    for (idx, unit_in_nanos) in enumerate(NanosPer)
        units = fld(ns, unit_in_nanos)
        ns -= units * unit_in_nanos
        result[idx] = units
    end

    result
end

#=

julia> tm = Time(now()); tmv = tm.instant.value;

julia> string(tm), tm, tmv
("21:05:38.181", Time(21, 5, 38, 181), 75938181000000)

julia> timeparts(tmv)
6-element Vector{Int128}:
  21
   5
  38
 181
   0
   0

julia> forparts = zeros(Int128,6);

julia> timeparts(tmv, forparts)
6-element Vector{Int128}:
  21
   5
  38
 181
   0
   0

julia> TimeInParts(timeparts(tmv, forparts))
(hours = 21, minutes = 5, seconds = 38, millis = 181, micros = 0, nanos = 0)
=#

function timeparts(nanos::S) where {S<:Signed}
    ns = nanos%Int128
    hr = fld(ns, NanosPerHour)
    ns -= hr * NanosPerHour
    mn = fld(ns, NanosPerMinute)
    ns -= mn * NanosPerMinute
    sc = fld(ns, NanosPerSecond)
    ns -= sc * NanosPerSecond
    ms = fld(ns, NanosPerMillisecond)
    ns -= ms * NanosPerMillisecond
    us = fld(ns, NanosPerMicrosecond)
    ns -= us * NanosPerMicrosecond
    (hr, mn, sc, ms, us, ns)
end


text:
using LinearAlgebra
using ExtendableSparse

mutable struct RadialBasis{F,Q,X,Y,L,U,C,S,D} <: AbstractSurrogate
    phi::F
    dim_poly::Q
    x::X
    y::Y
    lb::L
    ub::U
    coeff::C
    scale_factor::S
    sparse::D
end

mutable struct RadialFunction{Q,P}
    q::Q # degree of polynomial
    phi::P
end

linearRadial = RadialFunction(0,z->norm(z))

cubicRadial = RadialFunction(1,z->norm(z)^3)
multiquadricRadial = RadialFunction(1,z->sqrt(norm(z)^2+1))

thinplateRadial = RadialFunction(2, z->begin
    result = norm(z)^2 * log(norm(z))
    ifelse(iszero(z), zero(result), result)
end)

"""
    RadialBasis(x,y,lb::Number,ub::Number; rad::RadialFunction = linearRadial,scale::Real=1.0)

Constructor for RadialBasis surrogate.
"""
function RadialBasis(x, y, lb::Number, ub::Number; rad::RadialFunction=linearRadial, scale_factor::Real=1.0, sparse = false)
    q = rad.q
    phi = rad.phi
    coeff = _calc_coeffs(x, y, lb, ub, phi, q,scale_factor, sparse)
    return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor,sparse)
end

"""
RadialBasis(x,y,lb,ub,rad::RadialFunction, scale_factor::Float = 1.0)

Constructor for RadialBasis surrogate
"""
function RadialBasis(x, y, lb, ub; rad::RadialFunction = linearRadial, scale_factor::Real=1.0, sparse = false)
    q = rad.q
    phi = rad.phi
    coeff = _calc_coeffs(x, y, lb, ub, phi, q, scale_factor, sparse)
    return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor, sparse)
end

function _calc_coeffs(x, y, lb, ub, phi, q, scale_factor, sparse)
    nd = length(first(x))
    num_poly_terms = binomial(q + nd, q)

    D = _construct_rbf_interp_matrix(x, first(x), lb, ub, phi, q, scale_factor, sparse)
    Y = _construct_rbf_y_matrix(y, first(y), length(y) + num_poly_terms)

    coeff = D \ Y
    return coeff
end

function _construct_rbf_interp_matrix(x, x_el::Number, lb, ub, phi, q, scale_factor, sparse)
    n = length(x)

    num_poly_terms = binomial(q + 1, q)
    m = n + num_poly_terms

    if sparse
        D = ExtendableSparseMatrix{eltype(x_el),Int}(m,m)
    else
        D = zeros(eltype(x_el), m, m)
    end
    @inbounds for i = 1:n
        for j = 1:n
            D[i,j] = phi( (x[i] .- x[j]) ./ scale_factor )
        end
        if i <= n
            for k = 1:num_poly_terms
                    D[i,n+k] = _scaled_chebyshev(x[i], k-1, lb, ub)
            end
        end
    end
    D_sym = Symmetric(D, :U)
    return D_sym
end

function _construct_rbf_interp_matrix(x, x_el, lb, ub, phi, q, scale_factor,sparse)
    n = length(x)
    nd = length(x_el)

    num_poly_terms = binomial(q + nd, q)
    m = n + num_poly_terms

    if sparse
        D = ExtendableSparseMatrix{eltype(x_el),Int}(m,m)
    else
        D = zeros(eltype(x_el), m, m)
    end
    @inbounds for i = 1:n
        for j = 1:n
            D[i,j] = phi( (x[i] .- x[j]) ./ scale_factor)
        end
        if i < n + 1
            for k = 1:num_poly_terms
                D[i,n+k] = multivar_poly_basis(x[i], k-1, nd, q)
            end
        end
    end
    D_sym = Symmetric(D, :U)
    return D_sym
end

_construct_rbf_y_matrix(y, y_el::Number, m) = [i <= length(y) ? y[i] : zero(y_el) for i = 1:m]
_construct_rbf_y_matrix(y, y_el, m) = [i <= length(y) ? y[i][j] : zero(first(y_el)) for i=1:m, j=1:length(y_el)]

using Zygote: @nograd

function _make_combination(n, d, ix)
    exponents_combinations = [
        e
        for e
        in collect(
            Iterators.product(
                Iterators.repeated(0:n, d)...
            )
        )[:]
        if sum(e) <= n
    ]

    return exponents_combinations[ix + 1]
end
# TODO: Is this correct? Do we ever want to differentiate w.r.t n, d, or ix?
# By using @nograd we force the gradient to be 1 for n, d, ix
@nograd _make_combination

"""
    multivar_poly_basis(x, ix, d, n)

Evaluates in `x` the `ix`-th element of the multivariate polynomial basis of maximum
degree `n` and `d` dimensions.

Time complexity: `(n+1)^d.`

# Example
For n=2, d=2 the multivariate polynomial basis is
````
1,
x,y
x^2,y^2,xy
````
Therefore the 3rd (ix=3) element is `y` .
Therefore when x=(13,43) and ix=3 this function will return 43.
"""
multivar_poly_basis(x, ix, d, n) = prod(
    a^d
    for (a, d)
    in zip(x, _make_combination(n, d, ix))
)

"""
Calculates current estimate of value 'val' with respect to the RadialBasis object.
"""
function (rad::RadialBasis)(val)
    approx = _approx_rbf(val, rad)
    return _match_container(approx, first(rad.y))
end

function _approx_rbf(val::Number, rad)
    n = length(rad.x)
    q = rad.dim_poly
    num_poly_terms = binomial(q + 1, q)
    lb = rad.lb
    ub = rad.ub
    approx = zero(rad.coeff[1, :])
    for i = 1:n
        approx += rad.coeff[i, :] * rad.phi( (val .- rad.x[i]) / rad.scale_factor)
    end
    for k = 1:num_poly_terms
        approx += rad.coeff[n+k, :] * _scaled_chebyshev(val, k-1, lb, ub)
    end
    return approx
end
function _approx_rbf(val, rad)
    n = length(rad.x)
    d = length(rad.x[1])
    q = rad.dim_poly
    num_poly_terms = binomial(q + d, q)
    lb = rad.lb
    ub = rad.ub
    sum_half_diameter = sum((ub[k]-lb[k])/2 for k = 1:d)
    mean_half_diameter = sum_half_diameter/d
    central_point = _center_bounds(first(rad.x), lb, ub)

    approx = zero(rad.coeff[1, :])
    @views approx += sum( rad.coeff[i, :] * rad.phi( (val .- rad.x[i]) ./rad.scale_factor) for i = 1:n)
    for k = 1:num_poly_terms
        @views approx += rad.coeff[n+k, :] .* multivar_poly_basis(val, k-1, d, q)
    end
    return approx
end

_scaled_chebyshev(x, k, lb, ub) = cos(k*acos(-1 + 2*(x-lb)/(ub-lb)))
_center_bounds(x::Tuple, lb, ub) = ntuple(i -> (ub[i] - lb[i])/2, length(x))
_center_bounds(x, lb, ub) = (ub .- lb) ./ 2

"""
    add_point!(rad::RadialBasis,new_x,new_y)

Add new samples x and y and update the coefficients. Return the new object radial.
"""
function add_point!(rad::RadialBasis,new_x,new_y)
    if (length(new_x) == 1 && length(new_x[1]) == 1) || ( length(new_x) > 1 && length(new_x[1]) == 1 && length(rad.lb)>1)
        push!(rad.x,new_x)
        push!(rad.y,new_y)
    else
        append!(rad.x,new_x)
        append!(rad.y,new_y)
    end
    rad.coeff = _calc_coeffs(rad.x,rad.y,rad.lb,rad.ub,rad.phi,rad.dim_poly, rad.scale_factor, rad.sparse)
    nothing
end


text:
using BigCombinatorics
using Test

@testset "Fibonacci" begin
    @test Fibonacci(10) == Fibonacci(9) + Fibonacci(8)
    n = 5000
    a = Fibonacci(n)
    b = Fibonacci(n + 1)
    @test gcd(a, b) == 1
end

@testset "Factorial" begin
    @test Factorial(10) == big(factorial(10))
    @test FallingFactorial(10, 3) == 10 * 9 * 8
    @test FallingFactorial(10, 10) == Factorial(10)
    @test FallingFactorial(10, 12) == 0
    @test RisingFactorial(10, 3) == 10 * 11 * 12
    @test DoubleFactorial(9) == 9 * 7 * 5 * 3
    @test DoubleFactorial(10) == 10 * 8 * 6 * 4 * 2
    @test HyperFactorial(5) == prod(big(k)^big(k) for k = 1:5)

    @test Derangements(1) == 0
    @test Derangements(12) == 176214841
end

@testset "Binomial" begin
    @test sum(Binomial(10, k) for k = 0:10) == 2^10
    @test Catalan(12) == 208012
    @test MultiChoose(10, 1) == 10
    @test MultiChoose(10, 0) == 1
    @test MultiChoose(10, 10) == 92378

    @test Multinomial(5, 5, 5) == div(Factorial(15), Factorial(5)^3)
    @test Multinomial([6, 6, 6, 6]) == Multinomial(6, 6, 6, 6)
end

@testset "Partitions" begin
    @test Bell(10) == 115975
    @test sum(Stirling2(10, k) for k = 0:10) == Bell(10)
    @test Stirling1(10, 10) == 1
    @test Stirling1(10, 0) == 0
    @test sum(Stirling1(10, k) for k = 0:10) == 0

    @test IntPartitions(10) == 42
    @test IntPartitionsDistinct(10) == 10

    n = 20
    k = 5
    a = Stirling2(n, k)
    b = sum(Binomial(n - 1, j) * Stirling2(j, k - 1) for j = (k-1):(n-1))
    @test a == b

    c = sum((-1)^j * Binomial(k, j) * (k - j)^n for j = 0:k) ÷ Factorial(k)
    @test a == c
end

@testset "Euler[ian]" begin
    @test Euler(12) == 2702765
    @test sum(Eulerian(10, k) for k = 1:10) == Factorial(10)
end

@testset "Misc" begin
    BigCombinatorics.cache_clear()
    @test PowerSum(10, 3) == sum(k^3 for k = 1:10)
    @test PowerSum(100, 1) == sum(1:100)

    n = 10
    # Test Menage recursion by Touchard's formula
    T =
        2 *
        Factorial(n) *
        sum((-1)^k * 2n * Binomial(2n - k, k) * Factorial(n - k) ÷ (2n - k) for k = 0:n)
    @test T == Menage(n)
end


text:
using Plots
using DataFrames
using CSV

frameRate = 1/24
global timeInSec=0.0-frameRate
global secsToPlot = 3

function getDataFromFile()
    global timeInSec += frameRate
    df = CSV.read("data/"*string(timeInSec)*".csv")
    stats = convert(Matrix{Float64}, df)
    numOfBodies = size(stats)[1]
    xs = [stats[1,2]]
    ys = [stats[1,3]]
    zs = [stats[1,4]]
    ms = [stats[1,1]]
    for i in 2:numOfBodies
        push!(xs,stats[i,2])
        push!(ys,stats[i,3])
        push!(zs,stats[i,4])
        push!(ms,stats[i,1])
    end
    rm("data/"*string(timeInSec)*".csv")
    return xs, ys, zs, ms
end

function getMarkerSize(ms::Array{Float64,1})
    size = []
    for i in 1:length(ms)
        dot = log(ms[i]/M)
        push!(size, dot)
    end
    return size
end
function plotting()
    # initialize a 3D plot with 1 empty series
    global timeInSec = nothing
    global timeInSec = 0.0-frameRate
    xs, ys, zs, ms = getDataFromFile()
    len = length(xs)
    minX = minimum(xs)
    minY = minimum(ys)
    minZ = minimum(zs)
    maxX = maximum(xs)
    maxY = maximum(ys)
    maxZ = maximum(zs)
    plt = Plots.plot([xs[1]],[ys[1]], [zs[1]],
        seriestype=:scatter, marker= (:yellow, stroke(3, 0.2, :yellow)), markersize = 0, camera = (30, 30),background_color = :black,
        xlim = (minX - 100, maxX + 100),
        ylim = (minY - 100, maxY + 100),
        zlim = (minZ - 100, maxZ + 100)
    )
    size = getMarkerSize(ms[2:len])
    scatter!(xs[2:len],ys[2:len], zs[2:len],marker= (:white, stroke(3, 0.2, :white)), markersize = size, camera = (30,30),background_color = :black)

    display(plt)

    # build an animated gif by pushing new points to the plot, saving every 10th frame
    @gif for i=2:secsToPlot*24+1
        xs, ys, zs, ms = getDataFromFile()
        len = length(xs)

        plot([xs[1]],[ys[1]], [zs[1]],
            seriestype=:scatter, marker= (:yellow, stroke(3, 0.2, :yellow)), markersize = 0, camera = (30,30),background_color = :black,
            xlim = (minX - 1e2, maxX + 1e2),
            ylim = (minY - 1e2, maxY + 1e2),
            zlim = (minZ - 1e2, maxZ + 1e2)
        )

        size = getMarkerSize(ms[2:len])
        scatter!(xs[2:len],ys[2:len], zs[2:len],marker= (:white, stroke(3, 0.2, :white)), markersize  = size, camera = (30,30),background_color = :black)
    end

end

max_stars_repo_name: KeKeJin/Barnes-Hut-Galaxy-Simulator

max_issues_repo_name: KeKeJin/Barnes-Hut-Galaxy-Simulator

max_forks_repo_name: KeKeJin/Barnes-Hut-Galaxy-Simulator


text:
h1("JuLIP vs ASE Calculator Tests")
using Test, JuLIP.Testing, ASE, PyCall, JuLIP, LinearAlgebra

##

@info("These tests to compare JuLIP vs ASE implementations of some potentials")

h3("Compare JuLIP vs ASE: EMT")
# JuLIP's EMT implementation
at = set_pbc!( bulk(:Cu, cubic=true) * (2,2,2), (true,false,false) )
rattle!(at, 0.1)
emt = EMT()

@info("Test JuLIP vs ASE EMT implementation")
pyemt = ASE.Models.EMT()
print("   energy: ")
println(@test abs(energy(emt, at) - energy(pyemt, at)) < 1e-10)
print("   forces: ")
println(@test norm(forces(pyemt, at) - forces(emt, at), Inf) < 1e-10)

##

h3("Compare JuLIP vs ASE: EMT - Multi-species")
# JuLIP's EMT implementation
at = set_pbc!( bulk(:Cu, cubic=true) * (2,2,2), (true,false,false) )
at.Z[5:10] .= atomic_number(:Al)
@show unique(chemical_symbols(at))
rattle!(at, 0.1)
emt = EMT()

@info("Test JuLIP vs ASE EMT implementation")
pyemt = ASE.Models.EMT()
print("   energy: ")
println(@test abs(energy(emt, at) - energy(pyemt, at)) < 1e-10)
print("   forces: ")
println(@test norm(forces(pyemt, at) - forces(emt, at), Inf) < 1e-10)

##

h3("Compare JuLIP vs ASE: EAM")

# pot_file = joinpath(dirname(pathof(JuLIP)), "..", "data", "w_eam4.fs")
pot_file = @__DIR__() * "/w_eam4.fs"
url = "https://www.ctcms.nist.gov/potentials/Download/2013--Marinica-M-C-Ventelon-L-Gilbert-M-R-et-al--W-4/1/w_eam4.fs"
if !isfile(pot_file)
   @info("Download an Example EAM potential")
   run(`curl -o $pot_file $url`)
end

@info("Generate the ASE potential")
eam4_ase = ASE.Models.EAM(pot_file)

@info("Generate low-, med-, high-accuracy JuLIP potential")
eam4_jl1 = EAM(pot_file)
eam4_jl2 = EAM(pot_file; s = 1e-4)
eam4_jl3 = EAM(pot_file; s = 1e-6)

at1 = rattle!(bulk(:W, cubic=true) * 3, 0.1)
at2 = deleteat!(bulk(:W, cubic=true) * 3, 1)
at1_ase = ASEAtoms(at1)
at2_ase = ASEAtoms(at2)

for (i, (at, at_ase)) in enumerate(zip([at1, at2], [at1_ase, at2_ase]))
   @info("Test $i")
   err_low = (energy(eam4_ase, at_ase) - energy(eam4_jl1, at)) / length(at)
   err_med = (energy(eam4_ase, at_ase) - energy(eam4_jl2, at)) / length(at)
   err_hi = (energy(eam4_ase, at_ase) - energy(eam4_jl3, at)) / length(at)
   println("      Low Accuracy energy error:", err_low, "; ",
         (@test abs(err_low) < 0.03))
   println("   Medium Accuracy energy error:", err_med, "; ",
         (@test abs(err_med) < 0.006))
   println("     High Accuracy energy error:", err_hi, "; ",
         (@test abs(err_hi) < 0.0005))
end

##

h3("Compare JuLIP vs Quippy: SW")
if ASE.Models.has_quippy()
   swq = ASE.Models.QuippySW()
   sw = StillingerWeber()
   at0 = bulk(:Si) * (4, 4, 4)
   for k = 1:30
      at = rattle!(deepcopy(at0), 0.1)
      at_ase = ASEAtoms(at)
      err = (energy(swq, at_ase) - energy(sw, at))/length(at)
      print_tf(@test abs(err) < 1e-6)
   end
   println()
else
   @info("quippy was not found")
end


text:
"""
    JSFun

Wraps around a Julia function keeping track of its name and javascript
representation.
"""
mutable struct JSFun
    f::Function
    s::JSString
    name::Symbol
end
(jsf::JSFun)(x...) = jsf.f(x...)

"""
    @jsf

Make a function usable on a board. The function can still be used as a normal
Julia function.
"""
macro jsf(ex)
    if !(ex isa Expr) || ex.head ∉ (:function, :(=))
        error("@jsf only works with function definitions.")
    end
    fname = ex.args[1].args[1]
    ex.head = :function
    jsf = JSFun(identity, eval(:(@js($ex))), fname)
    ex.args[1].args[1] = Symbol("__jsf_$(fname)__")
    jsf.f = eval(ex)
    return esc(:($fname = $jsf))
end

function Base.show(io::IO, jsf::JSFun)
    println(io, "JSFun: $(jsf.name)")
    return nothing
end

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

"""
$SIGNATURES

Internal function to recover the javascript corresponding to a JSFun.
"""
function str(f::JSFun, b=nothing)
    isdef(f, b) && return ""
    rs = replace(f.s.s, r"\^\((.*?,.*?)\)" => s"Math.pow(\1)")
    return rs * ";"
end

"""
$SIGNATURES

Internal function to get separate strings corresponding to a JSFun to
facilitate the use of functions in plots etc. If a board is given, we first
check whether the given function has already been defined or not.
"""
function strf(f::JSFun, n::String="FN", b=nothing)
    fn   = String(f.name)
    jsfn = JSString(fn)
    jsn  = JSString("INSERT_$n")
    fs   = str(f, b)
    fss  = js"function(t){return $jsn(t);}"
    return fs, fss, "INSERT_$n" => fn
end
function strf(x::Union{Real,AbstractArray{<:Real}}, ::String="", b=nothing)
    fss = JSString("$x")
    return "", fss, nothing
end

"""
$SIGNATURES

Internal function to replace a function insertion.
"""
replacefn(s::AbstractString, ::Nothing) = s
replacefn(s::AbstractString, p::Pair) = replace(s, p)
function replacefn(s::AbstractString, ps...)
    for p in ps
        s = replacefn(s, p)
    end
    return s
end


text:
# canonical data version: 2.0.0

using Test

include("allergies.jl")

#canonical data
@testset "testing for single allergies" begin
    @testset "testing for eggs allergy" begin
        @test allergic_to(0, "eggs") == false
        @test allergic_to(1, "eggs") == true
        @test allergic_to(3, "eggs") == true
        @test allergic_to(2, "eggs") == false
        @test allergic_to(255, "eggs") == true
    end

    @testset "testing for peanuts allergy" begin
        @test allergic_to(0, "peanuts") == false
        @test allergic_to(2, "peanuts") == true
        @test allergic_to(7, "peanuts") == true
        @test allergic_to(5, "peanuts") == false
        @test allergic_to(255, "peanuts") == true
    end

    @testset "testing for shellfish allergy" begin
        @test allergic_to(0, "shellfish") == false
        @test allergic_to(4, "shellfish") == true
        @test allergic_to(14, "shellfish") == true
        @test allergic_to(10, "shellfish") == false
        @test allergic_to(255, "shellfish") == true
    end

    @testset "testing for strawberries allergy" begin
        @test allergic_to(0, "strawberries") == false
        @test allergic_to(8, "strawberries") == true
        @test allergic_to(28, "strawberries") == true
        @test allergic_to(20, "strawberries") == false
        @test allergic_to(255, "strawberries") == true
    end

    @testset "testing for tomatoes allergy" begin
        @test allergic_to(0, "tomatoes") == false
        @test allergic_to(16, "tomatoes") == true
        @test allergic_to(56, "tomatoes") == true
        @test allergic_to(40, "tomatoes") == false
        @test allergic_to(255, "tomatoes") == true
    end

    @testset "testing for chocolate allergy" begin
        @test allergic_to(0, "chocolate") == false
        @test allergic_to(32, "chocolate") == true
        @test allergic_to(112, "chocolate") == true
        @test allergic_to(80, "chocolate") == false
        @test allergic_to(255, "chocolate") == true
    end

    @testset "testing for pollen allergy" begin
        @test allergic_to(0, "pollen") == false
        @test allergic_to(64, "pollen") == true
        @test allergic_to(224, "pollen") == true
        @test allergic_to(160, "pollen") == false
        @test allergic_to(255, "pollen") == true
    end

    @testset "testing for cats allergy" begin
        @test allergic_to(0, "cats") == false
        @test allergic_to(128, "cats") == true
        @test allergic_to(192, "cats") == true
        @test allergic_to(64, "cats") == false
        @test allergic_to(255, "cats") == true
    end
end

@testset "testing the allergy_list" begin
        @testset "testing single allergies" begin
        @test allergy_list(0) == Set([])
        @test allergy_list(1) == Set(["eggs"])
        @test allergy_list(2) == Set(["peanuts"])
        @test allergy_list(8) == Set(["strawberries"])
    end

    @testset "testing for more allergies" begin
        @test allergy_list(3) == Set(["eggs", "peanuts"])
        @test allergy_list(5) == Set(["eggs", "shellfish"])
        @test allergy_list(248) == Set(["strawberries", "tomatoes", "chocolate", "pollen", "cats"])
        @test allergy_list(255) == Set(["eggs", "peanuts", "shellfish", "strawberries", "tomatoes", "chocolate", "pollen", "cats"])
        @test allergy_list(509) == Set(["eggs", "shellfish", "strawberries", "tomatoes", "chocolate", "pollen", "cats"])
    end
end


text:
# This file is a part of InverseFunctions.jl, licensed under the MIT License (MIT).

"""
    inverse(f)

Return the inverse of function `f`.

`inverse` supports mapped and broadcasted functions (via `Base.Fix1`) and
function composition (requires Julia >= 1.6).

# Examples

```jldoctest
julia> foo(x) = inv(exp(-x) + 1);

julia> inv_foo(y) = log(y / (1 - y));

julia> InverseFunctions.inverse(::typeof(foo)) = inv_foo;

julia> InverseFunctions.inverse(::typeof(inv_foo)) = foo;

julia> x = 4.2;

julia> inverse(foo)(foo(x)) ≈ x
true

julia> inverse(inverse(foo)) === foo
true

julia> broadcast_foo = Base.Fix1(broadcast, foo);

julia> X = rand(10);

julia> inverse(broadcast_foo)(broadcast_foo(X)) ≈ X
true

julia> bar = log ∘ foo;

julia> VERSION < v"1.6" || inverse(bar)(bar(x)) ≈ x
true
```

# Implementation

Implementations of `inverse(::typeof(f))` have to satisfy

* `inverse(f)(f(x)) ≈ x` for all `x` in the domain of `f`, and
* `inverse(inverse(f))` is defined and `inverse(inverse(f))(x) ≈ f(x)` for all `x` in the domain of `f`.

You can check your implementation with [`InverseFunctions.test_inverse`](@ref).
"""
inverse(f)
export inverse

"""
    struct NoInverse{F}

An instance `NoInverse(f)` signifies that `inverse(f)` is not defined.
"""
struct NoInverse{F}
    f::F
end
export NoInverse

(f::NoInverse)(x) = error("inverse of ", f.f, " is not defined")

inverse(f) = NoInverse(f)

inverse(f::NoInverse) = f.f

inverse(::typeof(inverse)) = inverse

@static if VERSION >= v"1.6"
    function inverse(f::Base.ComposedFunction)
        inv_inner = inverse(f.inner)
        inv_outer = inverse(f.outer)
        if inv_inner isa NoInverse || inv_outer isa NoInverse
            NoInverse(f)
        else
            Base.ComposedFunction(inv_inner, inv_outer)
        end
    end
end

function inverse(mapped_f::Base.Fix1{<:Union{typeof(map),typeof(broadcast)}})
    inv_f_kernel = inverse(mapped_f.x)
    if inv_f_kernel isa NoInverse
        NoInverse(mapped_f)
    else
        Base.Fix1(mapped_f.f, inverse(mapped_f.x))
    end
end

inverse(::typeof(identity)) = identity
inverse(::typeof(inv)) = inv
inverse(::typeof(adjoint)) = adjoint
inverse(::typeof(transpose)) = transpose
inverse(::typeof(conj)) = conj

inverse(::typeof(!)) = !
inverse(::typeof(+)) = +
inverse(::typeof(-)) = -

inverse(f::Base.Fix1{typeof(+)}) = Base.Fix2(-, f.x)
inverse(f::Base.Fix2{typeof(+)}) = Base.Fix2(-, f.x)
inverse(f::Base.Fix1{typeof(-)}) = Base.Fix1(-, f.x)
inverse(f::Base.Fix2{typeof(-)}) = Base.Fix1(+, f.x)
inverse(f::Base.Fix1{typeof(*)}) = iszero(f.x) ? throw(DomainError(f.x, "Cannot invert multiplication by zero")) : Base.Fix1(\, f.x)
inverse(f::Base.Fix2{typeof(*)}) = iszero(f.x) ? throw(DomainError(f.x, "Cannot invert multiplication by zero")) : Base.Fix2(/, f.x)
inverse(f::Base.Fix1{typeof(/)}) = Base.Fix2(\, f.x)
inverse(f::Base.Fix2{typeof(/)}) = Base.Fix2(*, f.x)
inverse(f::Base.Fix1{typeof(\)}) = Base.Fix1(*, f.x)
inverse(f::Base.Fix2{typeof(\)}) = Base.Fix1(/, f.x)

inverse(::typeof(deg2rad)) = rad2deg
inverse(::typeof(rad2deg)) = deg2rad

inverse(::typeof(exp)) = log
inverse(::typeof(log)) = exp

inverse(::typeof(exp2)) = log2
inverse(::typeof(log2)) = exp2

inverse(::typeof(exp10)) = log10
inverse(::typeof(log10)) = exp10

inverse(::typeof(expm1)) = log1p
inverse(::typeof(log1p)) = expm1

inverse(::typeof(sqrt)) = square
inverse(::typeof(square)) = sqrt

inverse(::typeof(cbrt)) = Base.Fix2(^, 3)
inverse(f::Base.Fix2{typeof(^)}) = iszero(f.x) ? throw(DomainError(f.x, "Cannot invert x^$(f.x)")) : Base.Fix2(invpow2, f.x)
inverse(f::Base.Fix2{typeof(invpow2)}) = Base.Fix2(^, f.x)
inverse(f::Base.Fix1{typeof(^)}) = Base.Fix1(invpow1, f.x)
inverse(f::Base.Fix1{typeof(invpow1)}) = Base.Fix1(^, f.x)
inverse(f::Base.Fix1{typeof(log)}) = Base.Fix1(invlog1, f.x)
inverse(f::Base.Fix1{typeof(invlog1)}) = Base.Fix1(log, f.x)
inverse(f::Base.Fix2{typeof(log)}) = Base.Fix2(invlog2, f.x)
inverse(f::Base.Fix2{typeof(invlog2)}) = Base.Fix2(log, f.x)


text:
using Test
using Pkg

@testset "inference" begin
    f() = Pkg.Types.STDLIBS_BY_VERSION
    @inferred f()
    f() = Pkg.Types.UNREGISTERED_STDLIBS
    @inferred f()
end

@testset "hashing" begin
    @test hash(Pkg.Types.Project()) == hash(Pkg.Types.Project())
    @test hash(Pkg.Types.VersionBound()) == hash(Pkg.Types.VersionBound())
    @test hash(Pkg.Resolve.Fixed(VersionNumber(0,1,0))) == hash(Pkg.Resolve.Fixed(VersionNumber(0,1,0)))

    hash(Pkg.Types.VersionSpec()) # hash isn't stable
    hash(Pkg.Types.PackageEntry()) # hash isn't stable because the internal `repo` field is a mutable struct
end


text:
module gRPC

using Sockets
using ProtoBuf
using HTTP2
using HTTP2.Session

import Base: close, run
import Sockets: TCPServer
import ProtoBuf: call_method
import HTTP2: Headers
import HTTP2.Session: HTTPConnection, EvtGoaway, EvtRecvHeaders, EvtRecvData

include("common.jl")
include("server.jl")
include("client.jl")

export gRPCController
export gRPCClient, close, stub
export gRPCServer, close, run

end # module


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: March 2022
# Copyright: Georgia Institute of Technology, 2022

export loss_supervised, amortized_MAP_loss!

function loss_supervised(
    Net::NetworkConditionalHINT,
    X::AbstractArray{Float32,4},
    Y::AbstractArray{Float32,4};
    grad::Bool = true,
)

    Zx, Zy, logdet = Net.forward(X, Y)
    if CUDA.functional()
        CUDA.reclaim()
    end
    z_size = size(Zx)

    f = sum(logpdf(0.0f0, 1.0f0, Zx))
    f = f + sum(logpdf(0.0f0, 1.0f0, Zy))
    f = f + logdet * z_size[4]

    if grad
        ΔZx = -gradlogpdf(0.0f0, 1.0f0, Zx) / z_size[4]
        ΔZy = -gradlogpdf(0.0f0, 1.0f0, Zy) / z_size[4]

        ΔX, ΔY = Net.backward(ΔZx, ΔZy, Zx, Zy)[1:2]
        if CUDA.functional()
            CUDA.reclaim()
        end
        GC.gc()

        return -f / z_size[4], ΔX, ΔY
    else
        return -f / z_size[4]
    end
end

function amortized_MAP_loss!(
    ΔZx::AbstractArray{Float32,4},
    Net::ReverseNetwork,
    Zx::AbstractArray{Float32,4},
    Zy::AbstractArray{Float32,4},
    Y_obs::Array{Float32,1},
    J::JUDI.judiJacobian{Float32,Float32},
    sigma::Float32,
    AN_x_rev::InvertibleNetworks.Reverse,
)

    X, Y = Net.forward(Zx, Zy)
    if CUDA.functional()
        CUDA.reclaim()
    end

    X = AN_x_rev.forward(wavelet_unsqueeze(X |> cpu))
    Zx = wavelet_unsqueeze(Zx |> cpu)

    shape = size(X)

    residual = J * convert(Array{Float32,1}, vec(X[:, :, 1, 1]')) - Y_obs
    f = -sum(logpdf(0.0f0, sigma, residual)) - sum(logpdf(0.0f0, 1.0f0, Zx))

    ΔX = reshape(
        reshape(-adjoint(J) * gradlogpdf(0.0f0, sigma, residual), shape[1], shape[2])',
        shape[1],
        shape[2],
        1,
        1,
    )

    ΔX = wavelet_squeeze(AN_x_rev.backward(ΔX, X)[1]) |> gpu

    X = wavelet_squeeze(AN_x_rev.inverse(X)) |> gpu
    Zx = wavelet_squeeze(Zx) |> gpu

    ΔZx .= Net.backward(ΔX, 0.0f0 * ΔX, X, Y)[1] - gradlogpdf(0.0f0, 1.0f0, Zx)
    if CUDA.functional()
        CUDA.reclaim()
    end

    clear_grad!(Net)
    GC.gc()

    return f / shape[4], X
end


text:
################################################################################
#
#  Small groups
#
################################################################################

#export number_of_small_groups, small_groups_limit
export small_group
#CF: as of now, nothing works here and the exports cause errors in Oscar

mutable struct SmallGroupDB
  path::String
  max_order::Int
  db::Vector{Vector{NamedTuple{(:name, :gens, :rels, :nontrivrels, :orderdis,
                                :ordersubdis, :isabelian, :iscyclic,
                                :issolvable, :isnilpotent, :autorder, :aut_gens,
                                :nchars, :dims, :schur, :galrep, :fields, :mod),
                               Tuple{String,Array{Perm{Int64},1},
                                     Vector{Vector{Int64}},
                                     Vector{Vector{Int64}},
                                     Vector{Tuple{Int64,Int64}},
                                     Vector{Tuple{Int64,Int64}},
                                     Bool,Bool,Bool,Bool,BigInt,
                                     Vector{Vector{Vector{Int64}}},Int64,
                                     Vector{Int64},Vector{Int64}, Vector{Int},
                                     Vector{Vector{Rational{BigInt}}},
                                     Vector{Vector{Vector{Vector{Rational{BigInt}}}}}}}}}

  function SmallGroupDB(path::String)
    db = Hecke.eval(Meta.parse(Base.read(path, String)))
    max_order = length(db)
    return new(path, max_order, db)
  end
end

# TODO: Write a parser for the data

function show(io::IO, L::SmallGroupDB)
  print(io, "Database of small groups (order limit = ", L.max_order, ")")
end

const default_small_group_db = [joinpath(pkgdir, "data/small_groups_extended"), joinpath(pkgdir, "data/small_groups_default")]

function small_group_database()
  for pa in default_small_group_db
    if isfile(pa)
      return SmallGroupDB(pa)
    end
  end
  throw(error("No database for small groups found"))
end

const DefaultSmallGroupDB = small_group_database()

isfrom_db(G::GrpGen) = G.isfromdb

function small_group(i, j; DB = DefaultSmallGroupDB)
  data = DB.db[i][j]
  #i < 1 || i > 63 && error("Group order ($i) must be between 1 and $small_groups_limit")
  #j < 1 || j > number_of_small_groups(i) && error("Index ($j) must be between 1 and $(number_of_small_groups(i))")
  P = SymmetricGroup(i)
  Gens = [P()]
  for p in data.gens
    push!(Gens, P(p))
  end
  G, PtoG, GtoP = generic_group(closure(Gens, *, P()), *)

  # The small groups have generators, for which we know
  # generators for the relations.
  # This is used in the homomorphism computation.
  # To make this work, the generators must much exactly
  # with the generators in the database.
  G.gens = Vector{Int}(undef, length(data.gens))
  for k in 1:length(data.gens)
    p = data.gens[k]
    G.gens[k] = PtoG[P(p)][]
  end

  G.issolvable = data.issolvable
  G.isnilpotent = data.isnilpotent
  G.isfromdb = true
  G.small_group_id = (i, j)
  return G
end


text:
using RDatasets, Gadfly

set_default_plot_size(9inch, 9inch)

df1 = dataset("ggplot2", "diamonds")
df2 = dataset("car", "SLID")

# colorful_hist, histogram_density, histogram_explicit_bins
p1 = plot(df1, x=:Price, color=:Cut, Geom.histogram)
# p2: see issue #880
p2 = plot(df1, x=:Price, color=:Cut, Geom.histogram(bincount=30, density=true))
p3 = plot(df1, x=:Price, Geom.histogram(bincount=30))

# stacked_continuous_histogram, stacked_discrete_histogram, stacked_discrete_histogram_horizontal
p4 = plot(df2, x=:Wages, color=:Language, Geom.histogram)
p5 = plot(df2, x=:Wages, color=:Language, Geom.histogram(position=:stack))
p6 = plot(df2, color=:Language, y=:Wages, Geom.histogram(position=:stack, orientation=:horizontal))

# dodged_discrete_histogram, dodged_discrete_histogram_horizontal
p7 = plot(df2, x=:Wages, color=:Language, Geom.histogram(position=:dodge))
p8 = plot(df2, y=:Wages, color=:Language, Geom.histogram(position=:dodge, orientation=:horizontal))
p9 = plot()

gridstack([p1 p2 p3; p4 p5 p6; p7 p8 p9])


text:
export
ruleSPTransitionOutNPP,
ruleSPTransitionIn1PNP,
ruleSPTransitionOutNCP,
ruleSPTransitionIn1CNP,
ruleVBTransitionOut,
ruleVBTransitionIn1,
ruleVBTransitionA,
ruleSVBTransitionOutVCD,
ruleSVBTransitionIn1CVD,
ruleSVBTransitionADV,
ruleMTransitionCCD,
ruleMTransitionCCN

# Note that vanishing probabilities are softened to prevent singularities
function ruleSPTransitionOutNPP(msg_out::Nothing,
                                msg_in1::Message{PointMass, Multivariate},
                                msg_a::Message{PointMass, MatrixVariate})

    a = clamp.(msg_a.dist.params[:m]*msg_in1.dist.params[:m], tiny, Inf)

    Message(Univariate, Categorical, p=a./sum(a))
end

function ruleSPTransitionIn1PNP(msg_out::Message{PointMass, Multivariate},
                                msg_in1::Nothing,
                                msg_a::Message{PointMass, MatrixVariate})

    a = clamp.(msg_a.dist.params[:m]'*msg_out.dist.params[:m], tiny, Inf)

    Message(Univariate, Categorical, p=a./sum(a))
end

function ruleSPTransitionOutNCP(msg_out::Nothing,
                                msg_in1::Message{Categorical, Univariate},
                                msg_a::Message{PointMass, MatrixVariate})

    a = clamp.(msg_a.dist.params[:m]*msg_in1.dist.params[:p], tiny, Inf)

    Message(Univariate, Categorical, p=a./sum(a))
end

function ruleSPTransitionIn1CNP(msg_out::Message{Categorical, Univariate},
                                msg_in1::Nothing,
                                msg_a::Message{PointMass, MatrixVariate})

    a = clamp.(msg_a.dist.params[:m]'*msg_out.dist.params[:p], tiny, Inf)

    Message(Univariate, Categorical, p=a./sum(a))
end

function ruleVBTransitionOut(   dist_out::Nothing,
                                dist_in1::ProbabilityDistribution,
                                dist_a::ProbabilityDistribution{MatrixVariate})

    a = clamp.(exp.(unsafeLogMean(dist_a)*unsafeMeanVector(dist_in1)), tiny, Inf)

    Message(Univariate, Categorical, p=a./sum(a))
end

function ruleVBTransitionIn1(   dist_out::ProbabilityDistribution,
                                dist_in1::Nothing,
                                dist_a::ProbabilityDistribution{MatrixVariate})

    a = clamp.(exp.(unsafeLogMean(dist_a)'*unsafeMeanVector(dist_out)), tiny, Inf)

    Message(Univariate, Categorical, p=a./sum(a))
end

ruleVBTransitionA(  dist_out::ProbabilityDistribution,
                    dist_in1::ProbabilityDistribution,
                    dist_a::Nothing) =
    Message(MatrixVariate, Dirichlet, a=unsafeMeanVector(dist_out)*unsafeMeanVector(dist_in1)' .+ 1)

function ruleSVBTransitionOutVCD(   dist_out::Nothing,
                                    msg_in1::Message{Categorical, Univariate},
                                    dist_a::ProbabilityDistribution{MatrixVariate})

    a = clamp.(exp.(unsafeLogMean(dist_a))*msg_in1.dist.params[:p], tiny, Inf)

    Message(Univariate, Categorical, p=a./sum(a))
end

function ruleSVBTransitionIn1CVD(   msg_out::Message{Categorical, Univariate},
                                    dist_in1::Nothing,
                                    dist_a::ProbabilityDistribution{MatrixVariate})

    a = clamp.(exp.(unsafeLogMean(dist_a))'*msg_out.dist.params[:p], tiny, Inf)

    Message(Univariate, Categorical, p=a./sum(a))
end

ruleSVBTransitionADV(   dist_out_in1::ProbabilityDistribution{Multivariate, Contingency},
                        dist_a::Nothing)=
    Message(MatrixVariate, Dirichlet, a=dist_out_in1.params[:p] .+ 1)

function ruleMTransitionCCD(msg_out::Message{Categorical, Univariate},
                            msg_in1::Message{Categorical, Univariate},
                            dist_a::ProbabilityDistribution{MatrixVariate})

    B = Diagonal(msg_out.dist.params[:p])*clamp.(exp.(unsafeLogMean(dist_a)), tiny, huge)*Diagonal(msg_in1.dist.params[:p])

    ProbabilityDistribution(Multivariate, Contingency, p=B./sum(B))
end

function ruleMTransitionCCN(msg_out::Message{Categorical, Univariate},
                            msg_in1::Message{Categorical, Univariate},
                            msg_a::Message{PointMass, MatrixVariate})

    B = Diagonal(msg_out.dist.params[:p])*msg_a.dist.params[:m]*Diagonal(msg_in1.dist.params[:p])

    ProbabilityDistribution(Multivariate, Contingency, p=B./sum(B))
end


text:
#Demo 2: Data patterns
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
tvsbit = cons(:a, labels = set(yaxis="Value", xaxis="Bit Position"))
vvst = cons(:a, labels = set(yaxis="Amplitude (V)", xaxis="Time (s)"))
color1 = cons(:a, line = set(color=:red))
color2 = cons(:a, line = set(color=:blue))

#==Input data
===============================================================================#
tbit = 1e-9 #Bit period
osr = 20 #samples per bit
nsamples = 20

#==Computations
===============================================================================#
seq = 1.0*prbs(reglen=5, seed=1, nsamples=nsamples)
#seq = [1,0,1,1,1,0,0,0]
nsamples = length(seq)
t = DataTime(0:(tbit/osr):(nsamples*tbit))
Π = pulse(t, Pole(3/tbit,:rad), npw=Index(osr))
pat = pattern(seq, Π, nbit=Index(osr))

#==Generate plot
===============================================================================#
p1 = push!(cons(:plot, tvsbit, title="PRBS Sequence"),
	cons(:wfrm, DataF1(collect(1:length(seq)), seq), color2),
)
p2 = push!(cons(:plot, vvst, title="PRBS Pattern"),
	cons(:wfrm, DataF1(Π), color1, label="Pulse"),
	cons(:wfrm, DataF1(pat), color2, label="Pattern"),
)

pcoll = push!(cons(:plotcoll, title="Generating Patterns"), p1, p2)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
# function tailiterate(X)
#     Q,L = ql(X)
#     I = Eye(2)
#     Z = Zeros(2,2)
#     [I Z; Z Z]*X + [Z Z; I Z]*L*
#                 [Z I Z; Z Z I; Z Z Z]
# end

# doesn't normalize last column
function _qlfactUnblocked!(A::AbstractMatrix{T}) where {T}
    require_one_based_indexing(A)
    m, n = size(A)
    τ = zeros(T, min(m,n))
    for k = min(m,n):-1:2
        ν = k+n-min(m,n)
        x = view(A, k:-1:1, ν)
        τk = reflector!(x)
        τ[k] = τk
        reflectorApply!(x, τk, view(A, k:-1:1, 1:ν-1))
    end
    QL(A, τ)
end

function tailiterate!(X::AbstractMatrix{T}) where T
    c,a,b = X[1,:]
    h = zero(T)
    for _=1:10_000_000
        QL = ql!(X)
        if h == X[1,3]
            return QL
        end
        h = X[1,3]
        X[2,:] .= (zero(T), X[1,1], X[1,2]);
        X[1,:] .= (c,a,b);
    end
    error("Did not converge")
end
tailiterate(c,a,b) = tailiterate!([c a b; 0 c a])

struct ContinuousSpectrumError <: Exception end

function qltail(Z::Number, A::Number, B::Number)
    T = promote_type(eltype(Z),eltype(A),eltype(B))
    ñ1 = (A + sqrt(A^2-4B*Z))/2
    ñ2 = (A - sqrt(A^2-4B*Z))/2
    ñ = abs(ñ1) > abs(ñ2) ? ñ1 : ñ2
    (n,σ) = (abs(ñ),sign(ñ))
    if n^2 < abs2(B)
        throw(ContinuousSpectrumError())
    end

    e = sqrt(n^2 - abs2(B))
    d = σ*e*Z/n

    Q =
    ql!([Z A B; 0 d e])
end

ql(A::SymTriPertToeplitz{T}; kwds...) where T = ql_hessenberg!(BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2))); kwds...)
ql(A::SymTridiagonal{T}; kwds...) where T = ql_hessenberg!(BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2))); kwds...)
ql(A::TriPertToeplitz{T}; kwds...) where T = ql_hessenberg!(BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2))); kwds...)
ql_hessenberg(A::InfBandedMatrix{T}; kwds...) where T = ql_hessenberg!(BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2))); kwds...)

toeptail(B::BandedMatrix{T}) where T =
    _BandedMatrix(B.data.args[end].args[1][1:end-B.u]*Ones{T}(1,∞), size(B,1), B.l-B.u, B.u)

# asymptotics of A[:,j:end] as j -> ∞
rightasymptotics(d::Hcat) = last(d.args)
rightasymptotics(d::Vcat) = Vcat(rightasymptotics.(d.args)...)
rightasymptotics(d) = d

function ql_hessenberg!(B::InfBandedMatrix{TT}; kwds...) where TT
    l,u = bandwidths(B)
    @assert u == 1
    T = toeptail(B)
    # Tail QL
    F∞ = ql_hessenberg(T; kwds...)
    Q∞, L∞ = F∞

    # populate finite data and do ql!
    data = bandeddata(B).args[1]
    B̃ = _BandedMatrix(data, size(data,2), l,u)
    B̃[end,end] = L∞[1,1]
    B̃[end:end,end-l+1:end-1] = adjoint(Q∞)[1:1,1:l-1]*T[l:2(l-1),1:l-1]
    F = ql!(B̃)

    # fill in data with L∞
    B̃ = _BandedMatrix(B̃.data, size(data,2)+l, l,u)
    B̃[size(data,2)+1:end,end-l+1:end] = adjoint(Q∞)[2:l+1,1:l+1] * T[l:2l,1:l]

    # combine finite and infinite data
    H = Hcat(B̃.data, rightasymptotics(F∞.factors.data))
    QLHessenberg(_BandedMatrix(H, ∞, l, 1), Vcat( LowerHessenbergQ(F.Q).q, F∞.q))
end

getindex(Q::QLPackedQ{T,<:InfBandedMatrix{T}}, i::Integer, j::Integer) where T =
    (Q'*Vcat(Zeros{T}(i-1), one(T), Zeros{T}(∞)))[j]'

getL(Q::QL, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = LowerTriangular(Q.factors)
getL(Q::QLHessenberg, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = LowerTriangular(Q.factors)

# number of structural non-zeros in axis k
nzzeros(A::AbstractArray, k) = size(A,k)
nzzeros(::Zeros, k) = 0
nzzeros(B::Vcat, k) = sum(size.(B.args[1:end-1],k))
nzzeros(B::CachedArray, k) = max(B.datasize[k], nzzeros(B.array,k))
function nzzeros(B::AbstractMatrix, k)
    l,u = bandwidths(B)
    k == 1 ? size(B,2) + l : size(B,1) + u
end

function lmul!(A::QLPackedQ{<:Any,<:InfBandedMatrix}, B::AbstractVecOrMat)
    require_one_based_indexing(B)
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    l,u = bandwidths(Afactors)
    D = Afactors.data
    for k = 1:∞
        ν = k
        allzero = k > nzzeros(B,1) ? true : false
        for j = 1:nB
            vBj = B[k,j]
            for i = max(1,ν-u):k-1
                Bij = B[i,j]
                if !iszero(Bij)
                    allzero = false
                    vBj += conj(D[i-ν+u+1,ν])*Bij
                end
            end
            vBj = A.τ[k]*vBj
            B[k,j] -= vBj
            for i = max(1,ν-u):k-1
                B[i,j] -= D[i-ν+u+1,ν]*vBj
            end
        end
        allzero && break
    end
    B
end

function lmul!(adjA::Adjoint{<:Any,<:QLPackedQ{<:Any,<:InfBandedMatrix}}, B::AbstractVecOrMat)
    require_one_based_indexing(B)
    A = adjA.parent
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    l,u = bandwidths(Afactors)
    D = Afactors.data
    @inbounds begin
        for k = nzzeros(B,1)+u:-1:1
            ν = k
            for j = 1:nB
                vBj = B[k,j]
                for i = max(1,ν-u):k-1
                    vBj += conj(D[i-ν+u+1,ν])*B[i,j]
                end
                vBj = conj(A.τ[k])*vBj
                B[k,j] -= vBj
                for i = max(1,ν-u):k-1
                    B[i,j] -= D[i-ν+u+1,ν]*vBj
                end
            end
        end
    end
    B
end

function (*)(A::QLPackedQ{T,<:InfBandedMatrix}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    lmul!(A, cache(convert(AbstractVector{TS},x)))
end

function (*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    lmul!(A, cache(convert(AbstractVector{TS},x)))
end

function blocktailiterate(c,a,b, d=c, e=a)
    z = zero(c)
    for _=1:1_000_000
        X = [c a b; z d e]
        F = ql!(X)
        d̃,ẽ = F.L[1:2,1:2], F.L[1:2,3:4]

        d̃,ẽ = QLPackedQ(F.factors[1:2,3:4],F.τ[1:2])*d̃,QLPackedQ(F.factors[1:2,3:4],F.τ[1:2])*ẽ  # undo last rotation
        if ≈(d̃, d; atol=1E-10) && ≈(ẽ, e; atol=1E-10)
            X[1:2,1:2] = d̃; X[1:2,3:4] = ẽ
            return PseudoBlockArray(X,fill(2,2), fill(2,3)), F.τ[3:end]
        end
        d,e = d̃,ẽ
    end
    error("Did not converge")
end

###
# BlockTridiagonal
####

function _ql(A::BlockTriPertToeplitz, d, e)
    N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
    c,a,b = A[Block(N+1,N)],A[Block(N,N)],A[Block(N-1,N)]
    P,τ = blocktailiterate(c,a,b,d,e)
    B = BlockBandedMatrix(A,(2,1))

    BB = _BlockBandedMatrix(B.data.args[1], (fill(2,N+2), fill(2,N)), (2,1))
    BB[Block(N),Block.(N-1:N)] .= P[Block(1), Block.(1:2)]
    F = ql!(view(BB, Block.(1:N), Block.(1:N)))
    BB[Block(N+1),Block.(N-1:N)] .= P[Block(2), Block.(1:2)]
    BB[Block(N+2),Block(N)] .= P[Block(2), Block.(1)]

    QL(_BlockSkylineMatrix(Vcat(BB.data, mortar(Fill(vec(Vcat(P[Block(1,3)], P[Block(2,3)], P[Block(2,2)], P[Block(2,1)])),∞))),B.block_sizes),
            Vcat(F.τ,mortar(Fill(τ,∞)))), P[Block(1,1)], P[Block(1,2)]
end

ql(A::BlockTriPertToeplitz) = _ql(A, A[Block(2,3)], A[Block(3,3)])[1]

ql(A::Adjoint{T,BlockTriPertToeplitz{T}}) where T = ql(BlockTridiagonal(A))

const InfBlockBandedMatrix{T} = BlockSkylineMatrix{T,<:Vcat{T,1,<:Tuple{Vector{T},<:BlockArray{T,1,<:Fill{<:Any,1,Tuple{OneToInf{Int64}}}}}}}

function lmul!(adjA::Adjoint{<:Any,<:QLPackedQ{<:Any,<:InfBlockBandedMatrix}}, B::AbstractVector)
    require_one_based_indexing(B)
    A = adjA.parent
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    l,u = blockbandwidths(Afactors)
    # todo: generalize
    l = 2l+1
    u = 2u+1
    @inbounds begin
        for k = nzzeros(B,1)+u:-1:1
            ν = k
            for j = 1:nB
                vBj = B[k,j]
                for i = max(1,ν-u):k-1
                    vBj += conj(Afactors[i,ν])*B[i,j]
                end
                vBj = conj(A.τ[k])*vBj
                B[k,j] -= vBj
                for i = max(1,ν-u):k-1
                    B[i,j] -= Afactors[i,ν]*vBj
                end
            end
        end
    end
    B
end

getindex(Q::QLPackedQ{T,<:InfBlockBandedMatrix{T}}, i::Integer, j::Integer) where T =
    (Q'*Vcat(Zeros{T}(i-1), one(T), Zeros{T}(∞)))[j]'

function (*)(A::QLPackedQ{T,<:InfBlockBandedMatrix}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    lmul!(A, cache(convert(AbstractVector{TS},x)))
end

function (*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBlockBandedMatrix}}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    lmul!(A, cache(convert(AbstractVector{TS},x)))
end

ldiv!(F::QLProduct, b::AbstractVector) = ldiv!(F.L, lmul!(F.Q',b))

function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',BandedColumns{PertConstRows}},<:PaddedLayout})
    A,b = M.A,M.B
    require_one_based_indexing(A, b)
    n = size(A, 2)
    if !(n == length(b))
        throw(DimensionMismatch("second dimension of left hand side A, $n, and length of right hand side b, $(length(b)), must be equal"))
    end
    data = triangulardata(A)
    nz = nzzeros(b,1)
    @inbounds for j in 1:n
        iszero(data[j,j]) && throw(SingularException(j))
        bj = b[j] = data[j,j] \ b[j]
        allzero = j > nz && iszero(bj)
        for i in (j+1:n) ∩ colsupport(data,j)
            b[i] -= data[i,j] * bj
            allzero = allzero && iszero(b[i])
        end
        allzero && break
    end
    b
end


text:
#export JULIA_NUM_THREADS=4 (before starting julia)

"""
parallelPeelAndBound(
    model::SequencingModel,
    maxWidth::U;

    orderingType::OrderingOption=edge,
    bestKnownSolution::Union{Array{U},Nothing}=nothing,
    bestKnownValue::Union{T,Nothing}=nothing,
    peelSetting::PeelSetting=peelf,
    rrbSetting::Union{RRBSetting,Nothing}=nothing,
    rrbPacket=nothing,
    timeLimit::Union{Nothing, W}=nothing,
    timeElapsed::X=0.0,
    loggingOn::Bool=true,
    debugOn::Bool=false,
    fileName::Union{Nothing,String}=nothing
)where{T<:Real,U<:Int, W<:Real, X<:Real}

Solve a sequencing problem to optimality using a peel and bound scheme
Returns bestKnownSolution (Vector{Int}), bestKnownValue (Real), bound (Real)

# Arguments
- `model::SequencingModel`: The sequencing problem
- `maxWidth::Int`: The max width allowed in a decision diagram layer
- `orderingType::OrderingOption`:Which heuristic to use to sort the edges, defaults to 'edge'
- `bestKnownSolution::Union{Vector{Int},Nothing}`: An optional parameter that gives the DD a starting solution
- `bestKnownValue::Union{T,Nothing}`: An optional parameter, The value of the best known feasible solution
- `peelSetting::PeelSetting`:Which heuristic to use to decide where to start the peel from, defaults to 'peelf'
- `rrbSetting::Union{RRBSetting,Nothing}`: An optional parameter that chooses an rough relaxed bounding method
- `rrbPacket`: Optional data to be passed to the rrbFunction
- `timeLimit::Union{Nothing, W}`: An optional parameter that adds a time limit
- `timeElapsed::X`: An optional parameter denoting how much time passed already
- `loggingOn::Bool`: An optional parameter that can be used to turn off logging to the console and local file
- `debugOn::Bool`: An optional parameter that can be used to turn on debug statements in other functions,
- `fileName::Union{Nothing,String}`: The file to write results to
"""
function parallelPeelAndBound(
    model::SequencingModel,
    maxWidth::U;
    orderingType::OrderingOption=edge,
    bestKnownSolution::Union{Array{U},Nothing}=nothing,
    bestKnownValue::Union{T,Nothing}=nothing,
    peelSetting::PeelSetting=peelf,
    rrbSetting::Union{RRBSetting,Nothing}=nothing,
    rrbPacket=nothing,
    timeLimit::Union{Nothing, W}=nothing,
    timeElapsed::X=0.0,
    loggingOn::Bool=true,
    debugOn::Bool=false,
    fileName::Union{Nothing,String}=nothing
)where{T<:Real,U<:Int, W<:Real, X<:Real}
    startTime = time()
    bestBound = nothing
    timedOut() = isTimedOut(timeLimit,startTime, timeElapsed)
    cleanRestrictedDD(node) = restrictedSequencingDD(model, maxWidth,rootNode = node, debugOn=debugOn,bestKnownSolution=bestKnownSolution, bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket)
    cleanRelaxedDD(dd) = relaxedSequencingDD(model, maxWidth,splitOrder, dd = dd, debugOn=debugOn,bestKnownSolution=bestKnownSolution, bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket, timeLimit=timeLimit, timeElapsed=timeElapsed, startTime=startTime)
    splitOrder = findSplitOrdering(model,orderingType)

    #get initial bounds
    if isnothing(bestKnownSolution)
        bestKnownSolution, bestKnownValue, exact = restrictedSequencingDD(model,2048,debugOn=debugOn)
        if exact
            if loggingOn
                finishPeelAndBound(bestKnownSolution, bestKnownValue,startTime=startTime, timeElapsed=timeElapsed)
            end
            return bestKnownSolution, bestKnownValue, bestKnownValue
        end
    end
    if timedOut()
        return bestKnownSolution, bestKnownValue, 0
    end
    #improve initial bounds using rrb
    if !isnothing(rrbSetting)
        bestKnownSolution, bestKnownValue, exact = cleanRestrictedDD(nothing)
        #if solution is found finish
        if exact
            if loggingOn
                finishPeelAndBound(bestKnownSolution, bestKnownValue,startTime=startTime, timeElapsed=timeElapsed)
            end
            return bestKnownSolution, bestKnownValue, bestKnownValue
        end
        if timedOut()
            return bestKnownSolution, bestKnownValue, 0
        end
    end
    if loggingOn
        logSequencingSolution(bestKnownSolution, bestKnownValue,firstLine="Initial Solution:")
        logRuntime(startTime, timeElapsed)
        println()
    end
    #generate initial relaxed DD
    relaxedDD, relaxedBound, relaxedSolution, bestKnownSolution, bestKnownValue, exact = cleanRelaxedDD(nothing)
    if exact || relaxedBound == bestKnownValue
        if loggingOn
            finishPeelAndBound(bestKnownSolution, bestKnownValue,startTime=startTime, timeElapsed=timeElapsed)
        end
        return bestKnownSolution, bestKnownValue, bestKnownValue
    else
        bestBound = relaxedBound
        if loggingOn
            logOptimalityGap(bestBound, bestKnownValue)
            logRuntime(startTime, timeElapsed)
            println()
        end
    end

    #time limit block
    if timedOut()
        return bestKnownSolution, bestKnownValue, bestBound
    end

    #initialize queue
    ddType = Vector{Vector{RelaxedSequencingNode}}
    queue = Vector{Tuple{ddType, typeof(bestKnownValue), Vector{Int}, Bool}}()
    push!(queue, (relaxedDD, bestBound, relaxedSolution, false))

    queueLock = ReentrantLock()
    countLock = ReentrantLock()
    solutionLock = ReentrantLock()
    runningJobs = Ref(0)
    bestKnownValue = Ref(bestKnownValue)
    #process queue
    while !isempty(queue) || runningJobs[] > 0
        while runningJobs[] >= Threads.nthreads()
            sleep(.001)
        end
        currentBlock = nothing
        begin
            lock(queueLock)
            lock(countLock)
            try
                if !isempty(queue)
                    currentBlock = pop!(queue)
                    runningJobs[] += 1
                    if loggingOn
                        println("Running Jobs: ", runningJobs[], " | Queue: ", length(queue), " | Local Bound: ", currentBlock[2])
                        logRuntime(startTime, timeElapsed)
                    end
                end
            finally
                unlock(queueLock)
                unlock(countLock)
            end
        end
        if !isnothing(currentBlock)
            Threads.@spawn asyncProcessNode!(model, maxWidth, queue, queueLock, solutionLock, countLock, runningJobs,currentBlock, copy(bestKnownSolution), bestKnownValue[], bestKnownSolution, bestKnownValue, splitOrder, peelSetting, rrbSetting, rrbPacket, timeLimit, timeElapsed, startTime, debugOn, loggingOn)
        end #process node
    end

    #log solution
    if loggingOn
        finishPeelAndBound(bestKnownSolution, bestKnownValue[], startTime=startTime, timeElapsed=timeElapsed)
    end
    return bestKnownSolution, bestKnownValue[], bestKnownValue[]
end

"""
    asyncProcessNode!(model::SequencingModel, maxWidth::Int, queue::Vector{Tuple{Vector{Vector{RelaxedSequencingNode}}, U, Vector{Int}, Bool}}, queueLock::ReentrantLock, solutionLock::ReentrantLock,countLock::ReentrantLock, jobCount::Ref{Int}, currentBlock::Tuple{Vector{Vector{RelaxedSequencingNode}}, U, Vector{Int}, Bool}, bestKnownSolution::Union{Vector{Int},Nothing}, bestKnownValue::Union{T,Nothing}, solutionRef::Union{Vector{Int},Nothing}, valueRef::Ref{T}, splitOrder::Vector{Int}, peelSetting::PeelSetting, rrbSetting::Union{RRBSetting,Nothing}, rrbPacket, timeLimit::Union{Nothing, W}, timeElapsed::X, startTime::Y, debugOn::Bool, loggingOn::Bool)where{T<:Real, U<:Real, W<:Real, X<:Real, Y<:Real}

An asynchronous call to the process the next node
Returns bestKnownSolution (Vector{Int}), bestKnownValue (Real), jobs added (Int)

# Arguments
- `model::SequencingModel`: The sequencing problem
- `maxWidth::Int`: The max width allowed in a decision diagram layer
- `queue`::Vector{Tuple{Vector{Vector{RelaxedSequencingNode}}, U, Vector{Int}, Bool}}: The node queue
- `queueLock::ReentrantLock`: The queue lock for parallel processing
- `solutionLock::ReentrantLock`:
- `countLock::::ReentrantLock`:
- `jobCount::Ref{Int}`:
- `currentBlock::Tuple{Vector{Vector{RelaxedSequencingNode}}, U, Vector{Int}, Bool}`: the dd to be processed
- `bestKnownSolution::Union{Vector{Int},Nothing}`: An optional parameter that gives the DD a starting solution
- `bestKnownValue::Union{T,Nothing}`: An optional parameter, The value of the best known feasible solution
- `solutionRef::Union{Vector{Int},Nothing}`:
- `valueRef::Ref{T}`:
- `splitOrder::Vector{Int}`: A heuristic ordering of nodes by importance
- `peelSetting::PeelSetting`:Which heuristic to use to decide where to start the peel from, defaults to 'peelf'
- `rrbSetting::Union{RRBSetting,Nothing}`: An optional parameter that chooses an rough relaxed bounding method
- `rrbPacket`: Optional data to be passed to the rrbFunction
- `timeLimit::Union{Nothing, W}`: An optional parameter that adds a time limit
- `timeElapsed::X`: An optional parameter denoting how much time passed already
- `startTime::Y`: An optional parameter denoting when the countdown started if there is a timer
- `loggingOn::Bool`: An optional parameter that can be used to turn off logging to the console and local file
- `debugOn::Bool`: An optional parameter that can be used to turn on debug statements in other functions
"""
function asyncProcessNode!(model::SequencingModel, maxWidth::Int, queue::Vector{Tuple{Vector{Vector{RelaxedSequencingNode}}, U, Vector{Int}, Bool}}, queueLock::ReentrantLock, solutionLock::ReentrantLock,countLock::ReentrantLock, jobCount::Ref{Int}, currentBlock::Tuple{Vector{Vector{RelaxedSequencingNode}}, U, Vector{Int}, Bool}, bestKnownSolution::Union{Vector{Int},Nothing}, bestKnownValue::Union{T,Nothing}, solutionRef::Union{Vector{Int},Nothing}, valueRef::Ref{T}, splitOrder::Vector{Int}, peelSetting::PeelSetting, rrbSetting::Union{RRBSetting,Nothing}, rrbPacket, timeLimit::Union{Nothing, W}, timeElapsed::X, startTime::Y, debugOn::Bool, loggingOn::Bool)where{T<:Real, U<:Real, W<:Real, X<:Real, Y<:Real}
    if peelSetting == peelf
        frontierNode, fNodeIndex = getFrontierNode(currentBlock[1], currentBlock[3])
    else
        frontierNode, fNodeIndex = getLastExactNode(currentBlock[1], currentBlock[3])
    end
    filterArcs!(frontierNode.arcsOut,model;bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket)
    bestKnownSolution, bestKnownValue, restrictedExact = restrictedSequencingDD(model, maxWidth,rootNode = frontierNode, debugOn=debugOn,bestKnownSolution=bestKnownSolution, bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket)

    #time limit block
    if isTimedOut(timeLimit,startTime, timeElapsed)
        return bestKnownSolution, bestKnownValue
    end

    #if node is closed
    if restrictedExact
        #delete the node
        ddInfo, bestPeeledSolution, bestPeeledValue = peelAndRemove!(model,currentBlock[1],frontierNode,fNodeIndex,bestKnownValue,rrbSetting,rrbPacket,loggingOn,debugOn)

        #if local bound
        if !ddInfo[4] && isBetterSolutionValue(model,ddInfo[2], bestKnownValue)
            begin
                lock(queueLock)
                try
                    addDDToQueue!(queue, ddInfo, model)
                finally
                    unlock(queueLock)
                end
            end
        elseif isBetterSolutionValue(model, bestPeeledValue, bestKnownValue)
            bestKnownValue = bestPeeledValue
            bestKnownSolution = bestPeeledSolution
        end
    else #if node is not closed
        #peel the dd
        ddInfo, newDDInfo, bestPeeledSolution, bestPeeledValue = peelDD!(model, maxWidth, currentBlock[1], currentBlock[3], frontierNode, fNodeIndex, splitOrder,bestKnownValue,rrbSetting,rrbPacket,loggingOn,debugOn)
        if isBetterSolutionValue(model, bestPeeledValue, bestKnownValue)
            bestKnownValue = bestPeeledValue
            bestKnownSolution = bestPeeledSolution
        end
        #add the DDs to the queue
        #if both nodes need to be added
        if !ddInfo[4] && isBetterSolutionValue(model,ddInfo[2], bestKnownValue) && !newDDInfo[4] && isBetterSolutionValue(model,newDDInfo[2], bestKnownValue)
            begin
                lock(queueLock)
                try
                    addDDToQueue!(queue, ddInfo, model)
                    addDDToQueue!(queue, newDDInfo, model)
                finally
                    unlock(queueLock)
                end
            end
        elseif !ddInfo[4] && isBetterSolutionValue(model,ddInfo[2], bestKnownValue)
            begin
                lock(queueLock)
                try
                    addDDToQueue!(queue, ddInfo, model)
                finally
                    unlock(queueLock)
                end
            end
        elseif !newDDInfo[4] && isBetterSolutionValue(model,newDDInfo[2], bestKnownValue)
            begin
                lock(queueLock)
                try
                    addDDToQueue!(queue, newDDInfo, model)
                finally
                    unlock(queueLock)
                end
            end
        end
    end
    #update solution value
    if isBetterSolutionValue(model, bestKnownValue, valueRef[])
        begin
            lock(solutionLock)
            try
                valueRef[] = bestKnownValue
                solutionRef = bestKnownSolution
                if loggingOn
                    logSequencingSolution(bestKnownSolution, bestKnownValue,firstLine="Improved Solution:")
                    logRuntime(startTime, timeElapsed)
                    println()
                end
            finally
                unlock(solutionLock)
            end
        end
    end
    #update job count
    begin
        lock(countLock)
        try
            jobCount[] -=1
        finally
            unlock(countLock)
        end
    end
    return bestKnownSolution, bestKnownValue
    #=if !isempty(queue)
        if loggingOn && isBetterSolutionValue(model, bestKnownValue, currentValue)
            logSequencingSolution(bestKnownSolution, bestKnownValue,firstLine="Improved Solution:")
            logOptimalityGap(last(queue)[2], bestKnownValue)
            logRuntime(startTime, timeElapsed)
        elseif loggingOn && isBetterSolutionValue(model, bestBound, last(queue)[2])
            println("Running Jobs: $running_jobs")
            logOptimalityGap(last(queue)[2], bestKnownValue)
            logRuntime(startTime, timeElapsed)
        end
        begin
            lock(boundLock)
            try
                bestBound = last(queue)[2]
            finally
                unlock(boundLock)
            end
        end
    end

    #time limit block
    isTimedOut(timeLimit,startTime, timeElapsed)
        return bestKnownSolution, bestKnownValue, bestBound
    end=#
end