text:
# # TEST 13H: square plate under harmonic loading. Parallel execution.

# Source code: [`TEST13H_par_tut.jl`](TEST13H_par_tut.jl)

# ## Description

# Harmonic forced vibration problem is solved for a homogeneous square plate,
# simply-supported on the circumference. This is the TEST 13H from the Abaqus v
# 6.12 Benchmarks manual. The test is recommended by the National Agency for
# Finite Element Methods and Standards (U.K.): Test 13 from NAFEMS “Selected
# Benchmarks for Forced Vibration,” R0016, March 1993.

# The plate is discretized with hexahedral solid elements. The simple support
# condition is approximated by distributed rollers on the boundary. Because
# only the out of plane displacements are prevented, the structure has three
# rigid body modes in the plane of the plate.

# Homogeneous square plate, simply-supported on the circumference from the test
# 13 from NAFEMS “Selected Benchmarks for Forced Vibration,” R0016, March 1993.
# The nonzero benchmark frequencies are (in hertz): 2.377, 5.961, 5.961, 9.483,
# 12.133, 12.133, 15.468, 15.468 [Hz].

# The magnitude of the displacement for the fundamental frequency (2.377 Hz) is
# 45.42mm according to the reference solution.

# ![](test13h_real_imag.png)

# The harmonic response loop is processed with multiple threads. The algorithm
# is embarrassingly parallel (i. e. no communication is required). Hence the
# parallel execution is particularly simple.

# ## References

# [1] NAFEMS “Selected Benchmarks for Forced Vibration,” R0016, March 1993.

# ## Goals

# - Show how to generate hexahedral mesh, mirroring and merging together parts.
# - Execute transient simulation by the trapezoidal-rule time stepping of [1].

##
# ## Definitions

# Bring in required support.
using FinEtools
using FinEtoolsDeforLinear
using LinearAlgebra
using Arpack

# The input parameters come from [1].
E = 200*phun("GPa");# Young's modulus
nu = 0.3;# Poisson ratio
rho = 8000*phun("kg*m^-3");# mass density
qmagn = 100.0*phun("Pa")
L = 10.0*phun("m"); # side of the square plate
t = 0.05*phun("m"); # thickness of the square plate
nL = 16; nt = 4;
tolerance = t/nt/100;
frequencies = vcat(linearspace(0.0,2.377,15), linearspace(2.377,15.0,40))

# Compute the parameters of Rayleigh damping. For the two selected
# frequencies we have the relationship between the damping ratio and
# the Rayleigh parameters
#
# $\xi_m=a_0/\omega_m+a_1\omega_m$
# where $m=1,2$.  Solving for the Rayleigh parameters $a_0,a_1$ yields:
zeta1 = 0.02; zeta2 = 0.02;
o1 = 2*pi*2.377;  o2 = 2*pi*15.468;
a0 = 2*(o1*o2)/(o2^2-o1^2)*(o2*zeta1-o1*zeta2);# a0
a1 = 2*(o1*o2)/(o2^2-o1^2)*(-1/o2*zeta1+1/o1*zeta2);# a1

##
# ## Discrete model

# Generate the finite element domain as a block.
fens,fes  = H8block(L, L, t, nL, nL, nt)

# Create the geometry field.
geom = NodalField(fens.xyz)
# Create the displacement field. Note that it holds complex numbers.
u = NodalField(zeros(FCplxFlt, size(fens.xyz,1), 3)) # displacement field

# In order to apply the essential boundary conditions we need to select the
# nodes along the side faces of the plate and support them in the Z direction.
nl = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
setebc!(u, nl, true, 3)
nl = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
setebc!(u, nl, true, 3)
nl = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
setebc!(u, nl, true, 3)
nl = selectnode(fens, box=[-Inf Inf L L -Inf Inf], inflate=tolerance)
setebc!(u, nl, true, 3)
# Those boundary conditions can now be applied to the displacement field,...
applyebc!(u)
# ... and the degrees of freedom can be numbered.
numberdofs!(u)
println("nfreedofs = $(u.nfreedofs)")

# The model is three-dimensional.
MR = DeforModelRed3D
material = MatDeforElastIso(MR, rho, E, nu, 0.0)

# Given how relatively thin the plate is we choose an effective element: the
# mean-strain hexahedral element which is quite tolerant of the high aspect
# ratio.
femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)), material)
# These elements require  to know the geometry before anything else can be
# computed using them in a finite element machine. Hence we first need to
# associate the geometry with the FEMM.
femm = associategeometry!(femm, geom)
# Now we can calculate the stiffness matrix and the mass matrix: both evaluated
# with the high-order Gauss rule.
K = stiffness(femm, geom, u)
femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
M = mass(femm, geom, u)
# The damping matrix is a linear combination of the mass matrix and the
# stiffness matrix (Rayleigh model).
C = a0*M + a1*K

# Find the boundary finite elements on top of the plate. The uniform distributed
# loading will be applied to these elements.
bdryfes = meshboundary(fes)
# Those facing up (in the positive Z direction) will be chosen:
topbfl = selectelem(fens, bdryfes, facing=true, direction=[0.0 0.0 1.0])
# A base finite element model machine will be created to evaluate the loading.
# The force intensity is created as driven by a function, but the function
# really only just fills the buffer with the constant loading vector.
function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
    forceout .=  [0.0, 0.0, qmagn]
    return forceout
end
fi = ForceIntensity(FFlt, 3, pfun);

# The loading vector is lumped from the distributed uniform loading by
# integrating on the boundary. Hence, the dimension of the integration domain
# is 2.
el1femm =  FEMMBase(IntegDomain(subset(bdryfes,topbfl), GaussRule(2,2)))
F = distribloads(el1femm, geom, u, fi, 2);

##
# ## Sweep through the frequencies

# Sweep through the frequencies and calculate the complex displacement vector
# for each  of the frequencies from the complex balance equations of the
# structure.

# The entire solution will be stored  in this array:
U1 = zeros(FCplxFlt, u.nfreedofs, length(frequencies))

# It is best to prevent the BLAS library from using threads concurrently with
# our own use of threads. The threads might easily become oversubscribed, with
# attendant slowdown.
LinearAlgebra.BLAS.set_num_threads(1)

# We utilize all the threads with which Julia was started. We can select the
# number of threads to use by running the executable as `julia -t n`, where `n`
# is the number of threads.
using Base.Threads
print("Number of threads: $(nthreads())\n")
print("Sweeping through $(length(frequencies)) frequencies\n")
t0 = time()
Threads.@threads for k in 1:length(frequencies)
    f = frequencies[k];
    omega = 2*pi*f;
    U1[:, k] = (-omega^2*M + 1im*omega*C + K)\F;
    print(".")
end
print("\nTime = $(time()-t0)\n")

# On Windows the scaling is not great, which is not Julia's fault, but rather
# the operating system's failing. Linux usually gives much greater speedups.

# Find the midpoint of the plate bottom surface.  For this purpose the number of
# elements along the edge of the plate needs to be divisible by two.
midpoint = selectnode(fens, box=[L/2 L/2 L/2 L/2 0 0], inflate=tolerance);
# Check that we found that node.
@assert midpoint != []
# Extract the displacement component in the vertical direction (Z).
midpointdof = u.dofnums[midpoint, 3]

##
# ## Plot the results

using Gnuplot

# Plot the amplitude of the FRF.
umidAmpl = abs.(U1[midpointdof, :])/phun("MM")
@gp "set terminal windows 0 " :-
@gp  :- vec(frequencies) vec(umidAmpl) "lw 2 lc rgb 'red' with lines title 'Displacement of the corner' " :-
@gp  :- "set xlabel 'Frequency [Hz]'" :-
@gp  :- "set ylabel 'Midpoint  displacement amplitude [mm]'"

# Plot the FRF real and imaginary components.
umidReal = real.(U1[midpointdof, :])/phun("MM")
umidImag = imag.(U1[midpointdof, :])/phun("MM")
@gp  "set terminal windows 1 "  :-
@gp  :- vec(frequencies) vec(umidReal) "lw 2 lc rgb 'red' with lines title 'Real' "  :-
@gp  :- vec(frequencies) vec(umidImag) "lw 2 lc rgb 'blue' with lines title 'Imaginary' "  :-
@gp  :- "set xlabel 'Frequency [Hz]'" :-
@gp  :- "set ylabel 'Midpoint  displacement FRF [mm]'"

# # Plot the shift of the FRF.
umidPhase = atan.(umidImag, umidReal)/pi*180
@gp  "set terminal windows 2 "  :-
@gp  :- vec(frequencies) vec(umidPhase) "lw 2 lc rgb 'red' with lines title 'Phase shift' "  :-
@gp  :- "set xlabel 'Frequency [Hz]'" :-
@gp  :- "set ylabel 'Phase shift [deg]'"

true


text:
@testset "sim" begin
    import StableRNGs
    rng = StableRNGs.StableRNG(1)
    stopping_time = fixed_time(3.2)
    n_sims = 100_000
    initial_condition = (X=zeros(n_sims),φ=ones(Int,n_sims))
    sims = simulate(model,stopping_time,initial_condition,rng)
    f(x,i) = cdf(sims)(x,i)
    @test sum(f.(10.0,phases(model)))≈1.0
    p_3_2 = ([1.0 0 0] * exp(model.T*3))[:]
    @test f.(10.0,phases(model))≈p_3_2 rtol=2e-3
end

ext: jl

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
"""
    Support for indicator constraints with a set constraint as the right hand side
"""
function JuMP._build_indicator_constraint(
    _error::Function,
    variable::JuMP.AbstractVariableRef,
    constraint::JuMP.VectorConstraint,
    ::Type{MOI.IndicatorSet{A}},
) where {A}
    S = typeof(constraint.set)
    F = typeof(JuMP.moi_function(constraint))
    set = CS.IndicatorSet{A,F,S}(constraint.set, 1 + length(constraint.func))
    if constraint.func isa Vector{VariableRef}
        vov = JuMP.VariableRef[variable]
    else
        vov = JuMP.AffExpr[variable]
    end
    append!(vov, constraint.func)
    return JuMP.VectorConstraint(vov, set)
end


text:
#   This file is part of Fatou.jl. It is licensed under the MIT license
#   Copyright (C) 2019 Michael Reed

import ColorSchemes

AbstractPlotting.plot(K::FilledSet;bare=false,args...) = AbstractPlotting.heatmap(K;bare=bare,args...)

for plt ∈ (:contour,:contourf,:heatmap)
    plt! = Symbol(plt,:!)
    @eval function AbstractPlotting.$plt(K::FilledSet;bare::Bool=false,args...)
        scene, layout = AbstractPlotting.layoutscene()
        ax = layout[1,1] = if K.meta.newt
            AbstractPlotting.Axis(scene,title=String(K),ylabel="Fatou set: z ↦ z-m×f(z)/f'(z)")
        else
            AbstractPlotting.Axis(scene,title=String(K))
        end
        plt = AbstractPlotting.$plt!(ax,K;args...)
        layout[1,2] = AbstractPlotting.Colorbar(scene,plt;width=30)
        return scene
    end
    @eval function AbstractPlotting.$plt!(ax,K::FilledSet;bare::Bool=false,args...)
        !haskey(args,:colormap) && !isempty(K.meta.cmap) && (return AbstractPlotting.$plt!(ax,K;bare=bare,colormap=Symbol(K.meta.cmap),args...))
        r1,r2 = ranges(K)
        Z = $(plt≠:heatmap ? :(transpose(K.meta.iter ? K.iter : K.mix)) : :(reverse(transpose(K.meta.iter ? K.iter : K.mix),dims=2)))
        if bare
            AbstractPlotting.$plt!(ax,r1,r2,Z;args...)
        else
            AbstractPlotting.$plt!(ax,r1,r2,Z;axis=(title=String(K),),args...)
        end
    end
end # center!(scene)

function AbstractPlotting.surface(K::FilledSet;bare::Bool=false,args...)
    !haskey(args,:colormap) && !isempty(K.meta.cmap) && (return AbstractPlotting.surface(K;bare=bare,colormap=Symbol(K.meta.cmap),args...))
    r1,r2 = ranges(K)
    Z = reverse(reverse(transpose(K.meta.iter ? K.iter : K.mix),dims=1),dims=2)
    if bare
        AbstractPlotting.surface(reverse(r1),r2,Z;(xreversed=true,)args...)
    else
        AbstractPlotting.surface(reverse(r1),r2,Z;axis=(title=String(K),xreversed=true),args...)
    end
end

function AbstractPlotting.arrows(K::FilledSet;bare::Bool=false,args...)
    r1,r2 = ranges(K)
    Z = transpose(K.set.Ω)
    AbstractPlotting.arrows(r1,r2,real.(Z),imag.(Z);args...)
end

function orbit(E,f::Function,bi::Matrix{Float64},orb::Int=0,depth::Int=1,incr::Int=384)
    x,N,N2,orbit,bis = real_orb(E,f,bi,orb,depth,incr)
    # setup plot
    funs,funt = ~(orb == 0) ? (["ϕ(x₀₋ₙ)"],", IC: x₀ = $(bis[3]), n∈0:$orb") : ([],"")
    scene, layout = AbstractPlotting.layoutscene()
    ax = layout[1,1] = AbstractPlotting.Axis(scene,title="x ↦ $E$funt")
    # plot background lines
    plt = [AbstractPlotting.lines!(ax,x[:],N[:,1];linestyle=:dash),
    AbstractPlotting.lines!(ax,x[:],N[:,2]),
    # plot orbit cobweb path
    AbstractPlotting.lines!(ax,orbit[:,1],orbit[:,2];color=:red)]
    # plot f^2,f^3,f^4,...
    for h ∈ 3:depth+1
        push!(plt,AbstractPlotting.lines!(ax,x,N[:,h];linewidth=1,color=ColorSchemes.tab10[h>5 ? h-1 : h-2]))
    end
    if orb≠0
        ran = range(bi[1],stop=bi[2],length=length(N2))
        plt = [plt...,[
            AbstractPlotting.lines!(ax,ran,N2[:];color=:gray,linestyle=:dot,linewidth=1),
            AbstractPlotting.scatter!(ax,ran,N2[:];color=:gray,marker=:x)]]
    end
    # trim graph
    d=1.07
    AbstractPlotting.xlims!(ax,bi[1],bi[2])
    AbstractPlotting.ylims!(ax,minimum([d*minimum(N[:,2]),0]),maximum([d*maximum(N[:,2]),0]))
    # set legend
    layout[1,2] = Legend(scene,plt,vcat(["y=x","ϕ(x)","(xₙ,ϕ(xₙ))"],["ϕ^$x(x)" for x ∈ 2:depth],funs))
    return scene
end


text:
const Maybe{T} = Union{Nothing, T}

abstract type ComoniconExpr end

Base.@kwdef struct Description <: ComoniconExpr
    brief::Union{Nothing, String} = nothing
    content::Union{Nothing, String} = nothing
end

Base.convert(::Type{Description}, ::Nothing) = Description()
Base.convert(::Type{Description}, x::String) = Description(x)
Base.convert(::Type{Description}, x::AbstractString) = Description(String(x))

Description(::Nothing) = Description(nothing, nothing)
function Description(text::String)
    return Description(brief(text), text)
end

Base.@kwdef struct Argument <: ComoniconExpr
    name::String
    type = Any
    vararg::Bool = false
    require::Bool = true
    # this is only for docs, since Julia
    # function will handle the actual default
    # value
    default::Maybe{String} = nothing
    description::Description = Description()
    line::Maybe{LineNumberNode} = nothing

    function Argument(name, type, vararg, require, default, description, line)
        require = vararg ? false : require # force require=false for vararg
        new(name, type, vararg, require, default, description, line)
    end
end

Base.@kwdef struct Option <: ComoniconExpr
    sym::Symbol
    name::String = replace(string(sym), '_'=>'-')
    hint::Maybe{String} = nothing
    type = Any
    short::Bool = false
    description::Description = Description()
    line::Maybe{LineNumberNode} = nothing
end

Base.@kwdef struct Flag <: ComoniconExpr
    sym::Symbol
    name::String = replace(string(sym), '_'=>'-')
    short::Bool = false
    description::Description = Description()
    line::Maybe{LineNumberNode} = nothing
end

Base.@kwdef struct NodeCommand <: ComoniconExpr
    name::String
    subcmds::Dict{String, Any}
    description::Description = Description()
    line::Maybe{LineNumberNode} = nothing

    function NodeCommand(name, subcmds, description, line)
        !isempty(subcmds) || error("list of subcommands should not be empty")
        new(name, subcmds, description, line)
    end
end

Base.@kwdef struct LeafCommand <: ComoniconExpr
    fn::Any
    name::String
    args::Vector{Argument} = Argument[]
    nrequire::Int = count(x->x.require, args)
    vararg::Maybe{Argument} = nothing
    flags::Dict{String, Flag} = Dict{String, Flag}()
    options::Dict{String, Option} = Dict{String, Option}()
    description::Description = Description()
    line::Maybe{LineNumberNode} = nothing

    function LeafCommand(fn, name, arg, nrequire, vararg,
            flags, options, description, line)

        isnothing(vararg) || vararg.vararg == true ||
            error("expect vararg $(vararg.name) " *
                "to have property vararg=true")
        new(fn, name, arg, nrequire, vararg,
            flags, options, description, line)
    end
end

Base.@kwdef struct Entry <: ComoniconExpr
    root::Union{NodeCommand, LeafCommand}
    version::Maybe{VersionNumber} = nothing
    line::Maybe{LineNumberNode} = nothing
end


text:
function psim_vrf(sample_function::Function, operands; diagnostics::Union(Vector{ASCIIString}, Nothing)=nothing)
  zvsamples = pmap(sample_function, operands)
  means = mean_of_zvmean(zvsamples)
  vars = var_of_zvmean(zvsamples, means=means)
  vrfs = vrf(vars)

  results = (ASCIIString=>Any)["means"=>means, "vars"=>vars, "vrfs"=>vrfs]
  if diagnostics!= nothing
    for d in diagnostics
      merge!(results, [d=>collect_diagnostic(zvsamples; diagnostic=d)])
    end
  end

  results
end

function psim_serialmc_vrf(tasks::Array{MCMCTask, 1};
  transform::Union(Function, Nothing)=nothing, rate::Bool=true,lzv::Bool=true, qzv::Bool=true)
  sample_function(t::MCMCTask) = sample_serialmc_zvmean(t::MCMCTask; transform=transform, rate=rate, lzv=lzv, qzv=qzv)
  diagnostics = rate == true ? ["rate"] : nothing
  psim_vrf(sample_function, tasks; diagnostics=diagnostics)
end

function psim_serialmc_vrf(t::MCMCTask; nchains::Int64=100,
  transform::Union(Function, Nothing)=nothing, rate::Bool=true,lzv::Bool=true, qzv::Bool=true)
  sample_function(x::MCMCTask) = sample_serialmc_zvmean(x::MCMCTask; transform=transform, rate=rate, lzv=lzv, qzv=qzv)
  diagnostics = rate == true ? ["rate"] : nothing
  tasks = MCMCTask[t for i = 1:nchains]
  psim_vrf(sample_function, tasks; diagnostics=diagnostics)
end

function psim_rand_vrf{T<:Distribution}(dists::Array{T, 1}; nsamples::Int64=100,
  transform::Union(Function, Nothing)=nothing, rate::Bool=true,lzv::Bool=true, qzv::Bool=true)
  sample_function(d::Distribution) = sample_rand_zvmean(d, nsamples; transform=transform, lzv=lzv, qzv=qzv)
  psim_vrf(sample_function, dists)
end

function psim_rand_vrf(d::Distribution; nsets::Int64=100, nsamples::Int64=100,
  transform::Union(Function, Nothing)=nothing, rate::Bool=true,lzv::Bool=true, qzv::Bool=true)
  sample_function(x::Distribution) = sample_rand_zvmean(x, nsamples; transform=transform, lzv=lzv, qzv=qzv)
  dists = typeof(d)[d for i = 1:nsets]
  psim_vrf(sample_function, dists)
end


text:
include("spawner.jl")
include("become_pingpong.jl")
include("functional_style.jl")

function run_suite(libtype::Type{<:ActorLib})
    run_functional_style_test(libtype())
    run_becomepingpong_test(libtype())
    run_spawnertest(libtype())
end


text:
using ResponseManneville
using ValidatedNumerics

D = ApproxInducedLSV(0.5, 10)
B = C2Basis(1024)
Q = DiscretizedOperator(B, D)

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
abstract type AbstractTrieKey end

struct TrieKey{E <: Tuple} <: AbstractTrieKey
    elements::E
    TrieKey(k::Tuple) = new{typeof(k)}(k)
    TrieKey(ks...) = TrieKey(ks)
end
Base.Tuple(k::TrieKey) = k.elements

Base.first(k::TrieKey) = first(Tuple(k))
Base.tail(k::TrieKey) = TrieKey(Base.tail(Tuple(k)))
_static(s::TrieKey) = StaticTrieKey(_static(s.elements))

struct StaticTrieKey{E} <: AbstractTrieKey
    StaticTrieKey(k::Tuple) = new{k}()
    StaticTrieKey(ks...) = StaticTrieKey(ks)
end
Base.Tuple(::StaticTrieKey{k}) where {k} = k

Base.first(::StaticTrieKey{k}) where {k} = first(k)
Base.tail(::StaticTrieKey{k}) where {k} = StaticTrieKey{Base.tail(k)}()

Base.size(k::AbstractTrieKey) = size(Tuple(k))
Base.length(k::AbstractTrieKey) = length(Tuple(k))
Base.eltype(k::AbstractTrieKey) = eltype(Tuple(k))
function Base.iterate(k::AbstractTrieKey, i::Integer=1)
    i <= length(k) && return Tuple(k)[i]
    return nothing
end

struct TrieMultiKey{K}
    keys::K
end

import JSON
Base.show(io::IO, key::TrieMultiKey) = print(io, JSON.json(key, 2))
Base.println(key::TrieMultiKey) = println(JSON.json(key, 2))
Base.print(key::TrieMultiKey) = print(JSON.json(key, 2))

export TrieKey, StaticTrieKey, TrieMultiKey

# Missing functionality
# - subsume
# - turn to graph for visualisation
# - allow tuple of indices keys - change to cartesian indices


text:
#=
  Test set 1: Construct a spine from 1 table (influenza cases) and link events where possible.
=#

println("\nTEST SET 1")

println("Write LinkageConfig")
lc1 = SpineBasedRecordLinkage.config.LinkageConfig(joinpath("config", "construct_spine_from_influenza_cases.toml"))
SpineBasedRecordLinkage.config.write_config(joinpath("output", "test_linkage_config.toml"), lc1)

println("Construct spine from influenza cases")
outdir1a = run_linkage(joinpath("config", "construct_spine_from_influenza_cases.toml"))
spine    = DataFrame(CSV.File(joinpath(outdir1a, "output", "spine.tsv"); delim='\t'))

@test size(spine, 1) == 4

println("Link influenza cases to events")
cp(joinpath(outdir1a, "output", "spine.tsv"), joinpath(outdir, "spine.tsv"); force=true)
outdir1b = run_linkage(joinpath("config", "link_influenza_cases_to_events.toml"))

println("Reporting")
outfile = joinpath(outdir, "linkage_report.csv")
summarise_linkage_run(outdir1b, outfile)
result  = table_to_set_of_dicts(outfile)

@test length(result) == 11
@test in(Dict(:tablename => "LINKAGE RUNS",            :status => outdir1b,                    :nrecords => missing), result)
@test in(Dict(:tablename => "spine",                   :status => "existent",                  :nrecords => 4), result)
@test in(Dict(:tablename => "emergency_presentations", :status => "linked with criteria ID 1", :nrecords => 2), result)
@test in(Dict(:tablename => "emergency_presentations", :status => "linked with criteria ID 2", :nrecords => 1), result)
@test in(Dict(:tablename => "emergency_presentations", :status => "unlinked",                  :nrecords => 2), result)
@test in(Dict(:tablename => "hospital_admissions",     :status => "linked with criteria ID 3", :nrecords => 3), result)
@test in(Dict(:tablename => "hospital_admissions",     :status => "unlinked",                  :nrecords => 2), result)
@test in(Dict(:tablename => "influenza_cases",         :status => "linked with criteria ID 4", :nrecords => 4), result)
@test in(Dict(:tablename => "influenza_cases",         :status => "linked with criteria ID 5", :nrecords => 2), result)
@test in(Dict(:tablename => "influenza_cases",         :status => "linked with criteria ID 6", :nrecords => 1), result)
@test in(Dict(:tablename => "influenza_cases",         :status => "linked with criteria ID 7", :nrecords => 1), result)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function content(blob::GitBlob)
    return ccall((:git_blob_rawcontent, :libgit2), Ptr{Void}, (Ptr{Void},), blob.ptr)
end

function Base.length(blob::GitBlob)
    return ccall((:git_blob_rawsize, :libgit2), Int64, (Ptr{Void},), blob.ptr)
end

function lookup(repo::GitRepo, oid::Oid)
    blob_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_blob_lookup, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ref{Oid}),
                   blob_ptr_ptr, repo.ptr, Ref(oid))
    return GitBlob(blob_ptr_ptr[])
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
"""
    `RBF_gauss(R;A=1.0,sigma=0.2,deriv=0)`

Kernel function for Gaussian RBF Interpolation. Returns the magnitude of evaluated
function.

The optional argument `deriv` is for the order of the derivative desired, meaning
0 for base, 1 for first derivative (gradient) and 2 for second derivative (Laplacian).

# ARGUMENTS
* `R::Float64` : Magnitude of radius vector -> R = |X_j - X_i|
* `A=1.0`      : Normalization value
* `sigma=0.2`  : Gaussian spreading
* `deriv=0`    : Order of desired derivative

# OUTPUTS
* `res`        : Calculated magnitude
"""
function RBF_gauss(R::Float64;
                   A=1.0,sigma=0.2,deriv=0)

    if deriv == 0
        # const1 = 1/(2*pi)^(3/2)
        # res = const1*exp(-R^2/2)

        res = SVPM.zeta_gauserf(R)
        #res = A*exp(-(R^2)/(2*sigma^2))
    elseif deriv == 1
        res = -((A*R)/(sigma^2))*exp(-(R^2)/(2*sigma^2))
    elseif deriv == 2
        res = ((A/sigma)^2)*((R/sigma)^2 - 3)*exp(-(R^2)/(2*sigma^2))
    end

    return res
end


text:
function interpolate_with_tspan(itp, tspan, t)
    t0 = tspan[1]
    t1 = tspan[end]
    N = length(tspan)
    if t < t0
        t_i = 1
    elseif t > t1
        t_i = N
    else
        t_i = (t-t0)/t1 * (N-1) + 1
    end
    return itp[t_i]
end

function M_aabb_W_bath_intp(t, s, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, W0, W0_bath, tspan, rho_0_int_t_itp, W_0_bath_t_itp, _ = p

    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)
    rho_t = interpolate_with_tspan(rho_0_int_t_itp, tspan, t)
    W_bath_s = interpolate_with_tspan(W_0_bath_t_itp, tspan, s)

    elLen = aggCore.molCount+1
    M_aabb_ = zeros(ComplexF64, elLen, elLen)
    for b=1:elLen
        rho_mantis = zeros(Float64, elLen, elLen)
        rho_mantis[b, b] = 1.

        tmp1[:, :] = ad(rho_mantis, W_bath_s, aggCore, aggTools)
        tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
        tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data
        M_tr = trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

        for a=1:elLen
            if rho_t[a, a] != 0.
                M_aabb_[a, b] = M_tr[a, a] / rho_t[a, a]
            end
        end
    end
    return M_aabb_
end

function K_aabb_W_bath_intp(t, p, tmp1, tmp2; rtol=1e-12, atol=1e-12)
    aggCore, aggTools, aggOperators, W0, W0_bath, tspan, rho_0_int_t_itp, W_0_bath_t_itp, elementtype = p
    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    M_aabb_integrated, err = QuadGK.quadgk(
        s -> M_aabb_W_bath_intp(t, s, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = rtol,
        atol = atol,
    )
    # K_const_ab_ = K_const_ab(t, p, tmp1, tmp2, Ham_II_t)
    # K_ab_ = -elementtype(im) * K_const_ab_ - M_aabb_integrated
    K_ab_ = - M_aabb_integrated
    return K_ab_
end

function K_aabb_W_bath_intp(t, W0, W0_bath, agg::Aggregate; rtol=1e-12, atol=1e-12)
    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))
    return K_aabb_W_bath_intp(t, p, tmp1, tmp2; rtol=rtol, atol=atol)
end

function M_abcd_W_bath_intp(t, s, p, tmp1, tmp2, Ham_II_t)
    aggCore, aggTools, aggOperators, W0, W0_bath, tspan, rho_0_int_t_itp, W_0_bath_t_itp, _ = p

    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham_II_s = getInteractionHamIPicture(Ham_0, Ham_I, s)
    rho_t = interpolate_with_tspan(rho_0_int_t_itp, tspan, t)
    W_bath_s = interpolate_with_tspan(W_0_bath_t_itp, tspan, s)

    elLen = aggCore.molCount+1
    M_abcd_ = zeros(ComplexF64, elLen, elLen, elLen, elLen)
    for c=1:elLen, d=1:elLen
        rho_mantis = zeros(Float64, elLen, elLen)
        rho_mantis[c, d] = 1.

        tmp1[:, :] = ad(rho_mantis, W_bath_s, aggCore, aggTools)
        tmp2[:, :] = Ham_II_s.data * tmp1 - tmp1 * Ham_II_s.data
        tmp1[:, :] = Ham_II_t.data * tmp2 - tmp2 * Ham_II_t.data
        M_tr = trace_bath(tmp1, aggCore, aggTools; vib_basis=aggOperators.vib_basis)

        for a=1:elLen, b=1:elLen
            if rho_t[a, b] != 0.
                M_abcd_[a, b, c, d] = M_tr[a, b] / rho_t[a, b]
            end
        end
    end
    return M_abcd_
end

function K_abcd_W_bath_intp(t, p, tmp1, tmp2; rtol=1e-12, atol=1e-12)
    aggCore, aggTools, aggOperators, W0, W0_bath, tspan, rho_0_int_t_itp, W_0_bath_t_itp, elementtype = p
    Ham_II_t = getInteractionHamIPicture(aggOperators.Ham_0, aggOperators.Ham_I, t)
    M_aabb_integrated, err = QuadGK.quadgk(
        s -> M_abcd_W_bath_intp(t, s, p, tmp1, tmp2, Ham_II_t),
        0,
        t,
        rtol = rtol,
        atol = atol,
    )
    # K_const_ab_ = K_const_ab(t, p, tmp1, tmp2, Ham_II_t)
    # K_ab_ = -elementtype(im) * K_const_ab_ - M_aabb_integrated
    K_ab_ = - M_aabb_integrated
    return K_ab_
end

function K_abcd_W_bath_intp(t, W0, W0_bath, agg::Aggregate; rtol=1e-12, atol=1e-12)
    tmp1 = copy(W0.data)
    tmp2 = copy(W0.data)
    p = (agg.core, agg.tools, agg.operators, W0, W0_bath, eltype(W0))
    return K_abcd_W_bath_intp(t, p, tmp1, tmp2; rtol=rtol, atol=atol)
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
###
### Reassemble: structural information -> system
###

function pantelides_reassemble(state::TearingState, var_eq_matching)
    fullvars = state.fullvars
    @unpack var_to_diff, eq_to_diff = state.structure
    sys = state.sys
    # Step 1: write derivative equations
    in_eqs = equations(sys)
    out_eqs = Vector{Any}(undef, nv(eq_to_diff))
    fill!(out_eqs, nothing)
    out_eqs[1:length(in_eqs)] .= in_eqs

    out_vars = Vector{Any}(undef, nv(var_to_diff))
    fill!(out_vars, nothing)
    out_vars[1:length(fullvars)] .= fullvars

    D = Differential(get_iv(sys))

    for (varidx, diff) in edges(var_to_diff)
        # fullvars[diff] = D(fullvars[var])
        vi = out_vars[varidx]
        @assert vi!==nothing "Something went wrong on reconstructing states from variable association list"
        # `fullvars[i]` needs to be not a `D(...)`, because we want the DAE to be
        # first-order.
        if isdifferential(vi)
            vi = out_vars[varidx] = diff2term(vi)
        end
        out_vars[diff] = D(vi)
    end

    d_dict = Dict(zip(fullvars, 1:length(fullvars)))
    lhss = Set{Any}([x.lhs for x in in_eqs if isdiffeq(x)])
    for (eqidx, diff) in edges(eq_to_diff)
        # LHS variable is looked up from var_to_diff
        # the var_to_diff[i]-th variable is the differentiated version of var at i
        eq = out_eqs[eqidx]
        lhs = if !(eq.lhs isa Symbolic)
            0
        elseif isdiffeq(eq)
            # look up the variable that represents D(lhs)
            lhsarg1 = arguments(eq.lhs)[1]
            @assert !(lhsarg1 isa Differential) "The equation $eq is not first order"
            i = get(d_dict, lhsarg1, nothing)
            if i === nothing
                D(eq.lhs)
            else
                # remove clashing equations
                lhs = Num(nothing)
            end
        else
            D(eq.lhs)
        end
        rhs = ModelingToolkit.expand_derivatives(D(eq.rhs))
        substitution_dict = Dict(x.lhs => x.rhs
                                 for x in out_eqs if x !== nothing && x.lhs isa Symbolic)
        sub_rhs = substitute(rhs, substitution_dict)
        out_eqs[diff] = lhs ~ sub_rhs
    end

    final_vars = unique(filter(x -> !(operation(x) isa Differential), fullvars))
    final_eqs = map(identity,
                    filter(x -> value(x.lhs) !== nothing,
                           out_eqs[sort(filter(x -> x !== unassigned, var_eq_matching))]))

    @set! sys.eqs = final_eqs
    @set! sys.states = final_vars
    return sys
end

"""
    pantelides!(state::TransformationState; kwargs...)

Perform Pantelides algorithm.
"""
function pantelides!(state::TransformationState; maxiters = 8000)
    @unpack graph, var_to_diff, eq_to_diff = state.structure
    neqs = nsrcs(graph)
    nvars = nv(var_to_diff)
    vcolor = falses(nvars)
    ecolor = falses(neqs)
    var_eq_matching = Matching(nvars)
    neqs′ = neqs
    for k in 1:neqs′
        eq′ = k
        isempty(𝑠neighbors(graph, eq′)) && continue
        pathfound = false
        # In practice, `maxiters=8000` should never be reached, otherwise, the
        # index would be on the order of thousands.
        for iii in 1:maxiters
            # run matching on (dx, y) variables
            #
            # the derivatives and algebraic variables are zeros in the variable
            # association list
            varwhitelist = var_to_diff .== nothing
            resize!(vcolor, nvars)
            fill!(vcolor, false)
            resize!(ecolor, neqs)
            fill!(ecolor, false)
            pathfound = construct_augmenting_path!(var_eq_matching, graph, eq′,
                                                   v -> varwhitelist[v], vcolor, ecolor)
            pathfound && break # terminating condition
            for var in eachindex(vcolor)
                vcolor[var] || continue
                # introduce a new variable
                nvars += 1
                add_vertex!(graph, DST)
                # the new variable is the derivative of `var`

                add_edge!(var_to_diff, var, add_vertex!(var_to_diff))
                push!(var_eq_matching, unassigned)
                var_derivative!(state, var)
            end

            for eq in eachindex(ecolor)
                ecolor[eq] || continue
                # introduce a new equation
                neqs += 1
                add_vertex!(graph, SRC)
                # the new equation is created by differentiating `eq`
                eq_diff = add_vertex!(eq_to_diff)
                add_edge!(eq_to_diff, eq, eq_diff)
                eq_derivative!(state, eq)
            end

            for var in eachindex(vcolor)
                vcolor[var] || continue
                # the newly introduced `var`s and `eq`s have the inherits
                # assignment
                var_eq_matching[var_to_diff[var]] = eq_to_diff[var_eq_matching[var]]
            end
            eq′ = eq_to_diff[eq′]
        end # for _ in 1:maxiters
        pathfound ||
            error("maxiters=$maxiters reached! File a bug report if your system has a reasonable index (<100), and you are using the default `maxiters`. Try to increase the maxiters by `pantelides(sys::ODESystem; maxiters=1_000_000)` if your system has an incredibly high index and it is truly extremely large.")
    end # for k in 1:neqs′
    return var_eq_matching
end

"""
    dae_index_lowering(sys::ODESystem; kwargs...) -> ODESystem

Perform the Pantelides algorithm to transform a higher index DAE to an index 1
DAE. `kwargs` are forwarded to [`pantelides!`](@ref). End users are encouraged to call [`structural_simplify`](@ref)
instead, which calls this function internally.
"""
function dae_index_lowering(sys::ODESystem; kwargs...)
    state = TearingState(sys)
    var_eq_matching = pantelides!(state; kwargs...)
    return invalidate_cache!(pantelides_reassemble(state, var_eq_matching))
end


text:
using Oscar
using Combinatorics
export makeSimplicial
export makeSmooth

function printn(A)
    print(A)
    print("\n")
end

"""
    slicematrix(::AbstractMatrix{<:Number})

    Take a two-dimensional matrix and output a list of its row vectors.

# Examples
```jldoctest makeSmoothWithDependencies
julia> A=[1 2; 3 4]

julia> slicematrix(A)
[[ 1 ,  2 ], [ 3 ,  4 ]]
"""

function slicematrix(A::AbstractMatrix{<:Number})
    return [A[i, :] for i in 1:size(A,1)]
end

"""
    rowMinors(::AbstractMatrix{<:Number},::Union{AbstractSet,AbstractVector})

    Identical to slicematrix, except only returns row vectors indexed by a set S.

# Examples
```jldoctest makeSmoothWithDependencies
julia> A=[1 2 3;4 5 6; 7 8 9]

julia> S=Set([1,3])

julia> rowMinors(A,S)
2×3 LinearAlgebra.Transpose{Int64,Array{Int64,2}}:
 1  2  3
 7  8  9
"""

function rowMinors(A::AbstractMatrix{<:Number},S::Union{AbstractSet,AbstractVector})
    outList=[]
    slices=slicematrix(A)
    for i in 1:size(slices,1)
        if i in S
            append!(outList,[slices[i]])
        end
    end
    return transpose(hcat(outList...))
end

"""
    convertIncidenceMatrix(::Polymake.IncidenceMatrixAllocated{Polymake.NonSymmetric})

    Takes a Polymake incidence matrix (e.g., the output of X.MAXIMAL_CONES for a toric variety X) and outputs a list of vectors,
    with each vector recording the indices marked on a given row of the incidence matrix.

# Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0;1 1; 0 1],INPUT_CONES=[[0,1],[1,2]])

julia> M=X.MAXIMAL_CONES

julia> convertIncidenceMatrix(M)
[[ 1 ,  2 ], [ 2 ,  3 ]]
"""

function convertIncidenceMatrix(A::Polymake.IncidenceMatrixAllocated{Polymake.NonSymmetric})
    A=Array(A)
    dim1=size(A,1)
    dim2=size(A,2)
    out=[]
    for i in 1:dim1
        members=[]
        for j in 1:dim2
            if A[i,j]==true
                append!(members,j)
            end
        end
        append!(out,[members])
    end
    return convert.(Array{Int64, 1}, out)
end

"""

    findBarycenter(::Union{AbstractSet,AbstractVector},::Polymake.BigObjectAllocated)

    Takes a normal toric variety X and a set s corresponding to a subset of rays of X, and outputs a polymake vector
    corresponding to the barycenter of those rays.

# Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0;1 1; 0 1],INPUT_CONES=[[0,1],[1,2]])

julia> s=[1,2]

julia> findBarycenter(s,X)
pm::Matrix<pm::Integer>
2 1

"""

function findBarycenter(s::Union{AbstractSet,AbstractVector},X::Polymake.BigObjectAllocated)
    rays = rowMinors(Array(X.RAYS), s)
    dim=size(rays,2)
    bary=zeros(Polymake.Rational,dim,1)
    for i in 1:size(rays,1)
        bary+=rays[i,:]
    end
    bary=Polymake.common.primitive(transpose(bary))
    return bary
end

"""

    coneConvert(::abstractVector{Int64},::abstractMatrix{Int64})

    Takes a matrix where the columns represent rays, and a list of indices, and forms a Polymake cone object generated by the rays corresponding to those indices.

# Examples
```jldoctest StackyFan

julia> typeof(coneConvert([1, 2, 4],[1 0 0; 0 1 0; 0 0 1; 1 1 1]))
Polymake.BigObjectAllocated

"""

function coneConvert(cone::Array{Int64,1},rayMatrix::Array)
    coneRays=rowMinors(rayMatrix,cone)
    C=Polymake.polytope.Cone(RAYS=coneRays)
    return C
end

"""

    getCones(X::Polymake.BigObjectAllocated)

    Returns all the cones of a fan X as a list of lists, with each interior list containing the indices of the rays generating a given cone.

# Examples
```jldoctest StackyFan
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 1 0; 1 0 1; 1 1 ],INPUT_CONES=[[0,1,2,3]])

julia> getCones(X)
[[ 0 ,  1 ,  2 ,  3 ], [ 0 ,  1 ], [ 0 ,  2 ], [ 2 ,  3 ], [ 1 ,  3 ], [ 0 ], [ 1 ], [ 2 ], [ 3 ]]

"""

function getCones(X::Polymake.BigObjectAllocated)
    lattice=X.HASSE_DIAGRAM
    faces=@Polymake.convert_to Array{Set{Int}} lattice.FACES
    out=[]
    for i in 2:(size(faces,1)-1)
        push!(out,Array(@Polymake.convert_to Array{Int} faces[i]))
    end
    return out
end

function findFaceContainingRay(C::Polymake.BigObjectAllocated,v::Array{Int64, 1})
    rayMatrix=Array(C.RAYS)
    faces=convertIncidenceMatrix(C.RAYS_IN_FACETS)
    if faces==[[]]
        return nothing
    end
    for face in faces
        faceCone=coneConvert(face,rayMatrix)
        if Polymake.polytope.contains(faceCone, v)
            return faceCone
        end
    end
    return nothing
end

function findMinimalCone(X::Polymake.BigObjectAllocated,v::Array{Int64, 1})
    rayMatrix=Array(X.RAYS)
    cones=convertIncidenceMatrix(X.MAXIMAL_CONES)
    #print(cones)
    startCone=nothing
    for cone in cones
        #print(cone)
        polyCone=coneConvert(cone,rayMatrix)
        if Polymake.polytope.contains(polyCone, v)
            startCone=polyCone
        end
    end
    if startCone==nothing
        error("The given ray is not contained in any cone of the fan.")
    end
    currentCone=startCone
    while(true)
        nextCone=findFaceContainingRay(currentCone,v)
        if nextCone==nothing
            break
        end
        currentCone=nextCone
    end
    currentRays=slicematrix(Array(currentCone.RAYS))
    fanRays=slicematrix(rayMatrix)
    indices=findall(x->x in currentRays,fanRays)
    return indices
end

function starSubdivision(X::Polymake.BigObjectAllocated, v::Array{Int64, 1})
    minimalCone = findMinimalCone(X, v)
    s = [i - 1 for i in minimalCone]
    v = transpose(v)
    return toric_blowup(s, X, v)
end

"""

    toric_blowup(::Union{AbstractSet,AbstractVector},::Polymake.BigObjectAllocated,::AbstractVector)

    Takes a normal toric variety X, a set s corresponding to a subset of rays of X, and a (optional) polymake vector,
    v, blow up X at v. If v is not provided, blow up X at the barycenter of s.

# Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2,3]])

julia>toric_blowup([0, 1, 2, 3], X, nothing)
INPUT_CONES
{0 1 4}
{0 2 4}
{2 3 4}
{1 3 4}
INPUT_RAYS
1 0 0
1 1 0
1 0 1
1 1 1
2 1 1

julia>toric_blowup([0], X, nothing)
INPUT_CONES
{0 2 3}
{0 1 3}
INPUT_RAYS
1 0 0
1 1 0
1 0 1
1 1 1

"""

function toric_blowup(s, X, v)
    s = [i + 1 for i in s]
    # If v is not provided, take v as the barycenter
    if v==nothing
        v=findBarycenter(s,X)
    end
    coneList = convertIncidenceMatrix(X.MAXIMAL_CONES)
    # Extracting the indices of all the cones in X that contains the set of rays s
    starIndex = findall((t) -> all(((i) -> i in t).(s)), coneList)
    star = [coneList[i] for i in starIndex]
    rayMatrix = X.RAYS

    lattice = X.HASSE_DIAGRAM
    faces = @Polymake.convert_to Array{Set{Int}} lattice.FACES

    # Get all the subcones of X that is contained in one of the cones in star and has the same rank
    clStar = []
    # Iterate over star
    for t in star
        # Get the rank of t
        c = rank(Array(rowMinors(rayMatrix, t))) - 1
        # Get all the subcones of X with rank c
        rank_c_subcone_indices = @Polymake.convert_to Array{Int} Polymake.graph.nodes_of_rank(lattice,c)
        rank_c_subcones = [faces[i + 1] for i in rank_c_subcone_indices]
        # Iterate over rank_c_subcones, and put the cones that is contained in t into clStar
        for cone in rank_c_subcones
            new_cone = [i+1 for i in cone]
            if all((i -> i in t).(new_cone))
                push!(clStar, new_cone)
            end
        end
    end
    # Remove duplicates
    clStar = unique(clStar)

    n = size(rayMatrix, 1) + 1
    # Filter out the cones in star from conelist
    coneList = filter(x -> !(x in star), coneList)

    if length(s) == 1
        # If s consists of a single ray, find all the cones in clStar that does not contain s
        newCones = []
        for t in clStar
            if !(s[1] in t)
                push!(newCones, sort(push!(t, s[1])))
            end
        end
        # return newCones plus coneList
        finalCones = [[i - 1 for i in cone] for cone in append!(coneList, newCones)]
        return Polymake.fulton.NormalToricVariety(INPUT_RAYS = Array(X.RAYS), INPUT_CONES = finalCones)
    end
    newCones = []
    for t in clStar
        # Find all the cones in clStar that does not contain at least one ray in s
        # QUESTION: Why seperate this from the one element case? Any won't work with one element list?
        if any(((i) -> !(i in t)).(s))
            push!(newCones, push!(t, n))
        end
    end
    # return newCones plus coneList
    finalRays = vcat((X.RAYS),v)
    finalCones = [[i - 1 for i in cone] for cone in append!(coneList, newCones)]
    return Polymake.fulton.NormalToricVariety(INPUT_RAYS = finalRays, INPUT_CONES = finalCones)
end

"""

    convertBool(::AbstractVector)

    Takes a column vector of boolean values and converts it to a vector of indices marked 'true'.

#Examples
```jldoctest makeSmoothWithDependencies
julia> B=[true true false true]

julia> convertBool(transpose(B))
[0, 1, 3]
"""

function convertBool(B::AbstractVector)
    out=[]
    for i in 1:size(B,1)
        if B[i]==true
           append!(out,i-1)
        end
    end
    return out
end

"""

    getConeRank(::AbstractMatrix,::AbstractVector)

    Takes a matrix and a vector containing indices corresponding to rows of a matrix,
    and calculates the rank of the matrix consisting only of those rows.

#Examples
```jldoctest makeSmoothWithDependencies
julia> v=[1,2]

julia> M=[0 1; 1 1; 1 0]

julia> getConeRank(v,M)
2
"""

function getConeRank(coneRayIndices::AbstractVector, rayMatrix::AbstractMatrix)
    coneRays = rowMinors(rayMatrix,coneRayIndices)
    return rank(Matrix(coneRays))
end

"""
    getDimension(::Polymake.BigObjectAllocated)

    Returns the ambient dimension of a normal toric variety.

#Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 2 0;0 0 1;0 1 0; 1 1 1],INPUT_CONES=[[0,1,2],[0,2,3,4]])

julia> getDimension(X)
3

"""

function getDimension(X)
    return size(X.RAYS, 2)
end

"""
    getConeFaces(::Polymake.BigObjectAllocated,::AbstractVector,::AbstractMatrix)

    Takes a fan, its ray matrix, and a vector corresponding to one of its cones, and returns a list of maximal strict faces of that cone.

#Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0; 1 1 0; 1 0 1; 1 1 1],INPUT_CONES=[[0,1,2,3]])

julia> getConeFaces(X,[1,2,3,4],Array(X.RAYS))
[[ 1 ,  2 ], [ 1 ,  3 ], [ 3 ,  4 ], [ 2 ,  4 ]]
"""

function getConeFaces(fan::Polymake.BigObjectAllocated,cone::AbstractVector,rayMatrix::AbstractMatrix)
    lattice = fan.HASSE_DIAGRAM
    faces = @Polymake.convert_to Array{Set{Int}} lattice.FACES
    cone_faces=[]
    c = rank(Array(rowMinors(rayMatrix, cone))) - 1
    rank_c_subcone_indices = @Polymake.convert_to Array{Int} Polymake.graph.nodes_of_rank(lattice,c)
    rank_c_subcones = [faces[i + 1] for i in rank_c_subcone_indices]
    for subcone in rank_c_subcones
        new_cone = [i+1 for i in subcone]
        if all((i -> i in cone).(new_cone))
            push!(cone_faces, new_cone)
        end
    end
    return cone_faces
end

"""
    makeSimplicial(::Polymake.BigObjectAllocated)

    Takes in a normal toric variety and returns a simplicial toric variety
    by subdividing (blowing up) the non-simplicial maximal cones.

#Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2,3]])

julia> X.SIMPLICIAL
false

julia> makeSimplicial(X).SIMPLICIAL
true

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0 0;0 1 0 0;0 0 1 0;1 -1 1 0; 1 0 -2 0],INPUT_CONES=[[0,1,2,3],[0,4]])

julia> X.SIMPLICIAL
false

julia> makeSimplicial(X).SIMPLICIAL
true
"""

function makeSimplicial(X::Polymake.BigObjectAllocated)
    Y = copy(X)
    while (true)
        # If the initial toric variety is simplicial, the program terminates and returns it.
        if Y.SIMPLICIAL==true
            break
        end
        #Maximal cones and ray matrix
        coneList = convertIncidenceMatrix(Y.MAXIMAL_CONES)
        rayMatrix = Y.RAYS
        badCone = nothing
        for i in 1:size(coneList,1)
            cone = coneList[i]
            if (getConeRank(cone, rayMatrix) != size(cone)[1])
                badCone = cone
            end
        end
        if (badCone == nothing)
            # All cones are linearly independent
            break
        else
            # Find the first ray that is contained in more than one orbit
            # and subdivide at that ray, using toricBlowup

            # Get faces (need to replace this)
            edges = getConeFaces(Y,badCone,rayMatrix)
            # Find the first ray that is contained in more than one orbit
            i = 1
            while count(r->(badCone[i] in r), edges) == 1
                i += 1
            end
            # Subdivide at the cone containing just that ray
            badCone=[i-1 for i in badCone]
            Y = toric_blowup(badCone, Y,nothing)
            #Y = toric_blowup([badCone[i]], Y,nothing)
        end
        # Repeat this process until there are no more bad cones
    end
    return Y
end

"""
    makeSmooth(::Polymake.BigObjectAllocated)

    Takes in a normal toric variety X and output a new smooth toric variety by iteratively blowing up.
    In the language of fans, these blowups are achieved by subdividing non-smooth cones.

#Examples
```jldoctest makeSmoothWithDependencies
julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[4 -1; 0 1],INPUT_CONES=[[0, 1]])

julia> X.SMOOTH_FAN
false

julia> makeSmooth(X).SMOOTH_FAN
true

julia> makeSmooth(X).INPUT_RAYS
pm::Matrix<pm::Rational>
1 -1/4
0 1
1 0

julia>X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;1 1 0;1 0 1;1 1 1],INPUT_CONES=[[0,1,2,3]])

julia>X.SMOOTH_FAN
false

julia>makeSmooth(X).SMOOTH_FAN
true

julia>makeSmooth(X).INPUT_RAYS
pm::Matrix<pm::Rational>
1 0 0
1 1 0
1 0 1
1 1 1
2 1 1

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 0;0 1 0;0 0 1;0 -1 -1; -1 0 -1; -2 -1 0],INPUT_CONES=[[0,1,2],[0,1,3],[1,3,4],[1,2,4],[2,4,5],[0,2,5],[0,3,5],[3,4,5]])

julia> X.SMOOTH_FAN
false

julia> makeSmooth(X).SMOOTH_FAN
true

"""

function makeSmooth(X::Polymake.BigObjectAllocated)
    Y  = copy(X)
    while(true)
        coneList = convertIncidenceMatrix(Y.MAXIMAL_CONES)
        rayMatrix = Array(Y.RAYS)
        k = 1
        # Iterate through the coneList, getting the index of the first cone not smooth
        for coneSet in coneList
            # Getting the number of rays in coneSet
            S=size(coneSet)[1]
            coneRays=rowMinors(rayMatrix,coneSet)
            # Checking whether this cone is smooth
            smoothCheck=Polymake.fan.check_fan_objects(Polymake.polytope.Cone(RAYS=coneRays)).SMOOTH_FAN
            if !smoothCheck
                # If the cone is not simplicial or not smooth, we have found the cone that we need to make smooth
                break
            else
                k+=1
            end
        end
        # At this point, all the cones are smooth. The program terminates.
        if k == size(coneList,1)+1
            break
        end

        # Get the cone that we found to be not smooth
        sigma=coneList[k]
        sigmaRays=slicematrix(rowMinors(rayMatrix,sigma))
        tau=0; tauRays=0; tauCone=0
        # Iterate over the subcones of sigma, finding tau, the smallest one that is not smooth
        for subset in collect(powerset(sigma))
            if size(subset,1) > 1
                S=size(subset)[1]
                subsetRays=rowMinors(rayMatrix,subset)
                subsetCone=Polymake.polytope.Cone(RAYS=subsetRays)
                smoothCheck=Polymake.fan.check_fan_objects(subsetCone).SMOOTH_FAN
                if !smoothCheck
                    tau=subset
                    tauRays=subsetRays
                    tauCone=subsetCone
                    break
                end
            end
        end

        # Getting the Hilbert Basis of tau
        H=slicematrix(Matrix(tauCone.HILBERT_BASIS_GENERATORS[1]))
        rayIndex=0
        # Iterate over the Hilbert Basis, finding the first ray that is not the generator of sigma
        for i in 1:size(H,1)
            if !(H[i] in sigmaRays)
                rayIndex=i
                break
            end
        end
        if rayIndex==0
            # Every Hilbert Basis of tau is a generator of sigma. Make Y simplicial is sufficient to make sigma smooth
            Y=makeSimplicial(Y)
        else
            # blowupRay is not a generator of sigma, blow up tau at blowupRay
            blowupRay=H[rayIndex]
            tau=[i-1 for i in tau]
            Y=toric_blowup(tau,Y,transpose(blowupRay))
        end
    end
    return Y
end

max_stars_repo_name: nbenabla/destackification

max_issues_repo_name: nbenabla/destackification

max_forks_repo_name: nbenabla/destackification


text:
export TRMA57_abs

function TRMA57_abs(nlp 		:: AbstractNLPModel,
              		nlpstop 	:: NLPStopping;
					kwargs...
               		)

	return TRARC(nlp, nlpstop; TR = TrustRegion(10.0),
				  c = Combi(hessian_sparse, PDataMA57{eltype(nlp.meta.x0)}, solve_modelTRDiagAbs, preprocessMA57, decreaseFact, Tparam{eltype(nlp.meta.x0)}()),
				  kwargs...
				  )
end


text:
abstract AbstractGeometry

import Base: size

include("getNodeProperty.jl")
include("cuboidRoom.jl")
#include("LShapedRoom.jl")

size(f::AbstractGeometry) = (f.Nx,f.Ny,f.Nz)

function Base.show(io::IO, f::AbstractGeometry)

	println(io, "geometry     : ", fun_name(f))
	println(io, "dimensions   : ", fun_dim(f))
	println(io, "samples      : ", fun_Nxyz(f))
	println(io, "ξ            : ", fun_ξ(f))
	println(io, "spatial step : ", fun_X(f))
	println(io, "sampling frq.: ", fun_Fs(f))

end

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


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

JLLWrappers.@generate_main_file_header("STL")
JLLWrappers.@generate_main_file("STL", UUID("69cf8a62-f9ae-5d65-8d00-30585a4771ef"))
end  # module STL_jll


text:
using PointProcessGeneration
using Test

@testset "PointProcessGeneration.jl" begin
    # Write your tests here.

    dim = 2
    n = 100
    ns = 10
    poisson(n; dim)

    colored_clusters(ns, n)

    normal(n)

    normal_clusters(ns, n; dim)
end


text:
# -*- coding: utf-8 -*-
import Revise
import GLRenderer as GL
import ThreeDP3 as T
import Images as I
import MiniGSG as S
import NearestNeighbors as NN
import PoseComposition: Pose, IDENTITY_POSE, IDENTITY_ORN
using Gen
import Rotations as R
import DataStructures as DS
import MeshCatViz

import MeshCatViz
MeshCatViz.setup_visualizer()

YCB_DIR = "/home/nishadg/mcs/ThreeDVision.jl/data/ycbv2"
world_scaling_factor = 100.0
id_to_cloud, id_to_shift, id_to_box  = T.load_ycbv_models_adjusted(YCB_DIR, world_scaling_factor);
all_ids = sort(collect(keys(id_to_cloud)));

IDX =400
@show T.get_ycb_scene_frame_id_from_idx(YCB_DIR,IDX)
gt_poses, ids, rgb_image, gt_depth_image, cam_pose, original_camera = T.load_ycbv_scene_adjusted(
    YCB_DIR, IDX, world_scaling_factor, id_to_shift
);
camera = T.scale_down_camera(original_camera, 6)
@show camera
img = I.colorview(I.Gray, gt_depth_image ./ maximum(gt_depth_image))
I.colorview(I.RGB, permutedims(Float64.(rgb_image)./255.0, (3,1,2)))

# +
resolution = 0.5

renderer = GL.setup_renderer(camera, GL.DepthMode())
for id in all_ids
    cloud = id_to_cloud[id]
    v,n,f = GL.mesh_from_voxelized_cloud(GL.voxelize(cloud, resolution), resolution)
    GL.load_object!(renderer, v, f)
end
depth_image = GL.gl_render(
    renderer, ids, gt_poses, T.IDENTITY_POSE)
depth_image[depth_image .== 50000.0] .= 200.0

function get_cloud_from_ids_and_poses(ids, poses, camera_pose)
    depth_image = GL.gl_render(renderer, ids, poses, camera_pose)
    cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(depth_image, camera))
    cloud = cloud[:,cloud[3,:] .< (camera.far - 1.0)]
    if size(cloud)[2] == 0
        cloud = zeros(3,1)
    end
    cloud
end

img = I.colorview(I.Gray, depth_image ./ maximum(depth_image))
# -

obs_cloud = GL.flatten_point_cloud(GL.depth_image_to_point_cloud(gt_depth_image, original_camera))
obs_cloud = GL.move_points_to_frame_b(obs_cloud, cam_pose)
obs_cloud =obs_cloud[:, 1.0 .< obs_cloud[3, :] .< 40.0];
obs_cloud =obs_cloud[:, -25.0 .< obs_cloud[1, :] .< 25.0];
obs_cloud =obs_cloud[:, -25.0 .< obs_cloud[2, :] .< 25.0];
voxelized_obs_cloud = T.voxelize(obs_cloud, resolution)
MeshCatViz.viz(voxelized_obs_cloud ./ 100.0)

get_cloud_func(id,p) = T.move_points_to_frame_b(T.voxelize(get_cloud_from_ids_and_poses(
                                [id], [p], cam_pose),resolution), cam_pose)

blobs = T.get_entities_from_assignment(voxelized_obs_cloud,
    T.dbscan_cluster(voxelized_obs_cloud, radius=0.6))

@show ids
obs_blob = blobs[4]
MeshCatViz.reset_visualizer()
MeshCatViz.viz(obs_blob ./ 100.0)

# +
parent_box = S.Box(40.0,40.0,0.1)
parent_pose = Pose(0.0, 0.0, 0.0)
parent_face = :top

c1_tree = NN.KDTree(obs_blob)
id_to_face_cloud_and_score =  DS.DefaultDict{Any,Vector{Any}}(()->[])

for id in all_ids
    print(id," ")
    for child_face in S.BOX_SURFACE_IDS
        child_box = id_to_box[id]
        (x,y,ang) = T.icp_project_to_planar_contact(
            obs_blob, p -> get_cloud_func(id,p),
            parent_pose, parent_box, parent_face,
            child_box, child_face,
            0.0, 0.0, 0.0;
            c1_tree= c1_tree, outer_iterations=10
        )

        contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang))
        p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)
        c = get_cloud_func(id,p)
        score = Gen.logpdf(
            T.uniform_mixture_from_template,
            obs_blob, c, 0.001, 0.5*2,(-1000.0, 1000.0, -1000.0,1000.0,-1000.0,1000.0)
        )

        push!(id_to_face_cloud_and_score[id], (face=child_face, cloud=c, score=score, pose=p, contact_params=(x,y,ang)))
    end
#     println(size(id_to_face_cloud_and_score[id]))
    sort!(id_to_face_cloud_and_score[id],by=x->x.score,rev=true)
end

# -

MeshCatViz.reset_visualizer()
MeshCatViz.viz(obs_blob ./100.0; color=I.colorant"black", channel_name=:obs_cloud)

# +
scores = [id_to_face_cloud_and_score[id][1].score for id in all_ids]
perm = sortperm(-scores)
best_score_id = perm[1]
@show best_score_id
c = id_to_face_cloud_and_score[best_score_id][1].cloud

MeshCatViz.reset_visualizer()
MeshCatViz.viz(obs_blob ./100.0; color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(c ./100.0; color=I.colorant"red", channel_name=:gen_cloud)
# -

for id in all_ids
    c = id_to_face_cloud_and_score[id][1].cloud
    MeshCatViz.reset_visualizer()
    MeshCatViz.viz(obs_blob ./100.0; color=I.colorant"black", channel_name=:obs_cloud)
    MeshCatViz.viz(c ./100.0; color=I.colorant"red", channel_name=:gen_cloud)
    sleep(0.1)
end

# +
c = id_to_face_cloud_and_score[6][1].cloud

MeshCatViz.reset_visualizer()
MeshCatViz.viz(obs_blob ./100.0; color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(c ./100.0; color=I.colorant"red", channel_name=:gen_cloud)

# +

id_and_id_to_face_cloud_and_score =  DS.DefaultDict{Any,Vector{Any}}(()->[])

for id_1 in all_ids
    print(id_1," ")

    data = id_to_face_cloud_and_score[id_1][1]

    unexplained_idxs = T.get_unexplained_points(obs_blob, data.cloud)
    unexplained_points = obs_blob[:, unexplained_idxs]

    c1_tree = NN.KDTree(unexplained_points)

    parent_pose = data.pose
    parent_face = T.opposite_face(data.face)
    parent_box = id_to_box[id_1]

    for id in all_ids
        for child_face in S.BOX_SURFACE_IDS
            child_box = id_to_box[id]

            (x,y,ang) = T.icp_project_to_planar_contact(
                unexplained_points, p -> get_cloud_func(id,p), parent_pose, parent_box, parent_face,
                child_box, child_face,
                0.0,0.0,0.0;
                c1_tree= c1_tree
            )

            contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang))
            p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)
            c = get_cloud_func(id,p)
            score = Gen.logpdf(
                T.uniform_mixture_from_template,
                unexplained_points, c, 0.01, 0.5*2,(-1000.0, 1000.0, -1000.0,1000.0,-1000.0,1000.0)
            )

            push!(id_and_id_to_face_cloud_and_score[(id_1,id)],
                (face=child_face, cloud=c, score=score, pose=p, contact_params=(x,y,ang)))
        end
        sort!(id_and_id_to_face_cloud_and_score[(id_1,id)],by=x->x.score,rev=true)
    end
end

# -

id_and_id_to_score = Dict()
for id_1 in all_ids
    for id in all_ids
        pose_1 = id_to_face_cloud_and_score[id_1][1].pose
        pose_2 = id_and_id_to_face_cloud_and_score[(id_1,id)][1].pose

        c = T.move_points_to_frame_b(T.voxelize(get_cloud_from_ids_and_poses(
                                [id_1,id], [pose_1, pose_2], cam_pose), 0.5), cam_pose)

        score = Gen.logpdf(
            T.uniform_mixture_from_template,
            obs_blob, c, 0.01, 0.5*2,(-1000.0, 1000.0, -1000.0,1000.0,-1000.0,1000.0)
        )
        id_and_id_to_score[(id_1,id)] = score
    end
end

# +
(id_1,id) = argmax(id_and_id_to_score)
@show (id_1,id)
pose_1 = id_to_face_cloud_and_score[id_1][1].pose
pose_2 = id_and_id_to_face_cloud_and_score[(id_1,id)][1].pose

c = T.move_points_to_frame_b(T.voxelize(get_cloud_from_ids_and_poses(
                        [id_1,id], [pose_1, pose_2], cam_pose),0.5), cam_pose)
MeshCatViz.viz(obs_blob ./100.0; color=I.colorant"black", channel_name=:obs_cloud)
MeshCatViz.viz(c ./100.0; color=I.colorant"red", channel_name=:gen_cloud)
# -

# +
renderer_texture = GL.setup_renderer(original_camera, GL.TextureMode())
obj_paths = T.load_ycb_model_obj_file_paths(YCB_DIR)
texture_paths = T.load_ycb_model_texture_file_paths(YCB_DIR)
for id in all_ids
    v,n,f,t = renderer_texture.gl_instance.load_obj_parameters(
        obj_paths[id]
    )
    v = v * world_scaling_factor
    v .-= id_to_shift[id]'

    GL.load_object!(renderer_texture, v, n, f, t,
        texture_paths[id]
    )
end

# -

@show ids

img = I.colorview(I.RGB, permutedims(Float64.(rgb_image)./255.0, (3,1,2)))

# +
parent_box = S.Box(40.0,40.0,0.1)
parent_pose = Pose(0.0, 0.0, 0.0)
parent_face = :top
id = 6
child_box = id_to_box[id]

x = 17.0
y = -6.0
ang = Gen.uniform(0.0, 2*π)
child_face = :top

contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang))
p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)

renderer_texture.gl_instance.lightpos = cam_pose.pos
@time rgb_image, depth_image = GL.gl_render(
    renderer_texture, [id], [p], cam_pose)
overlay_img = I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))
mask = (overlay_img .!= overlay_img[1,1])
new_img = copy(img)
alpha = 0.1
new_img[mask] .= alpha*new_img[mask] .+ (1-alpha).*overlay_img[mask]
new_img

# +
parent_box = S.Box(40.0,40.0,0.1)
parent_pose = Pose(0.0, 0.0, 0.0)
parent_face = :top
id = 6
child_box = id_to_box[id]

param_sweep = [
    (x,y,0.0)
    for x in -20.0:2.5:13.0 for y in -10.0:2.5:10.0
]
child_face = :top
@show size(param_sweep)
imgs = [
    let
        contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang))
        p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)

        renderer_texture.gl_instance.lightpos = cam_pose.pos
        rgb_image, depth_image = GL.gl_render(
            renderer_texture, [id], [p], cam_pose)
        overlay_img = I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))
        mask = (overlay_img .!= overlay_img[1,1])
        new_img = copy(img)
        alpha = 0.1
        new_img[mask] .= alpha*new_img[mask] .+ (1-alpha).*overlay_img[mask]
        new_img
    end
    for (x,y,ang) in param_sweep
];
# -

import FileIO
gif = cat(imgs...;dims=3);
FileIO.save("test.gif", gif)

# +
new_img = copy(img)
function highlight_bounding_box(img_in, is,js, w, h)
    img = copy(img_in)
    col = I.colorant"red"
    border_width = 5
    for i in is:is+border_width
        for j in js:js+border_width

        img[i,j:j+w] .= col
        img[i+h,j:j+w] .= col
        img[i:i+h,j] .= col
        img[i:i+h,j+w] .= col
        end
    end
    img
end

param_sweep = [
    [
    highlight_bounding_box(new_img, x,y,z,z)
    for x in 1:50:(480-z-5) for y in 1:50:(640-z-5)
    ]
    for z in 100:50:200
];
param_sweep = vcat(param_sweep...)

@show size(param_sweep)
# -

import FileIO
gif = cat(param_sweep...;dims=3);
FileIO.save("test2.gif", gif)

Gen.categorical(ones(6)./6.0)

imgs = [
    let
        x = Gen.uniform(-20.0,13.0)
        y = Gen.uniform(-10.0,10.0)
        ang = Gen.uniform(0,2*pi)
        child_face = S.BOX_SURFACE_IDS[Gen.categorical(ones(6)./6.0)]

        contact = S.ShapeContact(parent_face, Real[], child_face, Real[], S.PlanarContact(x,y,ang))
        p = parent_pose * S.getRelativePoseFromContact(parent_box, child_box, contact)

        renderer_texture.gl_instance.lightpos = cam_pose.pos
        rgb_image, depth_image = GL.gl_render(
            renderer_texture, [id], [p], cam_pose)
        overlay_img = I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))
        mask = (overlay_img .!= overlay_img[1,1])
        new_img = copy(img)
        alpha = 0.1
        new_img[mask] .= alpha*new_img[mask] .+ (1-alpha).*overlay_img[mask]
        new_img
    end
    for _ in 1:50
];

import FileIO
gif = cat(imgs...;dims=3);
FileIO.save("test2.gif", gif)

imgs = [
    let
        p = T.uniformPose(-30.0, 33.0, -20.0,20.0, -10.0, 20.0)
        renderer_texture.gl_instance.lightpos = cam_pose.pos
        rgb_image, depth_image = GL.gl_render(
            renderer_texture, [id], [p], cam_pose)
        overlay_img = I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))
        mask = (overlay_img .!= overlay_img[1,1])
        new_img = copy(img)
        alpha = 0.1
        new_img[mask] .= alpha*new_img[mask] .+ (1-alpha).*overlay_img[mask]
        new_img
    end
    for _ in 1:50
];

import FileIO
gif = cat(imgs...;dims=3);
FileIO.save("test2.gif", gif)


text:
#!/usr/bin/env julia

import Pkg
#Pkg.add("Plots")
#Pkg.add("GR")
#Pkg.add("PlotlyBase")

using Plots

x = 1:10; y = rand(10, 2) # 2 columns means two lines
println(x)
println(y)
#plotly() # Set the backend to Plotly
# This plots into the web browser via Plotly
#gr()
plotly()

plt = plot(x, y, title = "This is Plotted using Plotly")
display(plt)

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
mutable struct NormAttributes
    affine::Bool
    track_stats::Bool
    active::Union{Bool,Nothing}
end

Flux.hasaffine(attr::NormAttributes) = attr.affine

function get_stats(::Val{true}, ::Val{false}, l, x::AbstractArray{T,N}, reduce_dims) where {T,N}
    # testmode with tracked stats
    stats_shape = ntuple(i -> i == N - 1 ? size(x, N - 1) : 1, N)
    return reshape(l.μ, stats_shape), reshape(l.σ², stats_shape)
end

function get_stats(::Val{false}, active, l, x, reduce_dims)
    # trainmode or testmode without tracked stats
    μ = mean(x; dims=reduce_dims)
    diff = x .- μ
    return μ, mean(abs2, diff; dims=reduce_dims)
end

function get_stats(::Val{true}, active::Val{true}, l, x::AbstractArray{T,N}, reduce_dims) where {T,N}
    # trainmode with tracked stats
    μ, σ² = get_stats(Val(false), active, l, x, reduce_dims)
    Zygote.ignore() do
        mtm = l.momentum
        m = prod(size(x)[reduce_dims])  # needed for computing corrected var
        μnew = vec(N ∈ reduce_dims ? μ : mean(μ; dims=N))
        σ²new = vec(N ∈ reduce_dims ? σ² : mean(σ²; dims=N))
        l.μ .= (1 - mtm) .* l.μ .+ mtm .* μnew
        return l.σ² .= (1 - mtm) .* l.σ² .+ mtm .* (m / (m - one(eltype(l.σ²)))) .* σ²new
    end
    return μ, σ²
end

function norm_forward(l, x::AbstractArray{T,N}, reduce_dims, affine_shape) where {T,N}
    μ, σ² = get_stats(Val(l.attrs.track_stats), Val(_isactive(l)), l, x, reduce_dims)
    if hasaffine(l)
        γ = reshape(l.γ, affine_shape)
        β = reshape(l.β, affine_shape)
        return l.λ.(_norm_forward(μ, σ², x, γ, β, l.ϵ))
    else
        return l.λ.(_norm_forward(μ, σ², x, l.ϵ))
    end
end

_norm_forward(μ, σ², x, ϵ) = (x .- μ) ./ sqrt.(σ² .+ ϵ)

_norm_forward(μ, σ², x, γ, β, ϵ) = γ .* (x .- μ) ./ sqrt.(σ² .+ ϵ) .+ β

Zygote.@adjoint function _norm_forward(μ, σ², x, γ, β, ϵ)
    N = ndims(x)

    σ²ϵ = σ² .+ ϵ
    inv_deno = 1 ./ sqrt.(σ²ϵ)
    res_1 = (x .- μ) .* inv_deno
    res_2 = γ .* res_1
    res = res_2 .+ β

    function norm_backward(Δ)
        reduce_dims_affine = filter(i -> isone(size(β, i)), 1:N)
        reduce_dims_stats = filter(i -> isone(size(σ², i)), 1:N)

        Δx = inv_deno .* Δ
        Δμ = -sum(Δx; dims=reduce_dims_stats)
        Δσ² = sum(-eltype(x)(0.5) .* res_2 .* Δ ./ σ²ϵ; dims=reduce_dims_stats)
        Δγ = sum(res_1 .* Δ; dims=reduce_dims_affine)
        Δβ = sum(Δ; dims=reduce_dims_affine)

        return (Δμ, Δσ², Δx, Δγ, Δβ, nothing)
    end

    return res, norm_backward
end

Zygote.@adjoint function _norm_forward(μ, σ², x, ϵ)
    N = ndims(x)

    σ²ϵ = σ² .+ ϵ
    inv_deno = 1 ./ sqrt.(σ²ϵ)
    res = (x .- μ) .* inv_deno

    function norm_backward(Δ)
        reduce_dims_stats = filter(i -> isone(size(σ², i)), 1:N)

        Δx = inv_deno .* Δ
        Δμ = -sum(Δx; dims=reduce_dims_stats)
        Δσ² = sum(-eltype(x)(0.5) .* res .* Δ ./ σ²ϵ; dims=reduce_dims_stats)

        return (Δμ, Δσ², Δx, nothing)
    end

    return res, norm_backward
end


text:
"""
A two dimensional nonlinear oscillator introduced in [1].
Note that the system is constructed by converse HJB method [2].
# References
[1] J. A. Primbs, “Nonlinear Optimal Control: A Receding Horizon Approach,”
California Institute of Technology, Pasadena, California, 1999.
[2] J. Doyle, J. A. Primbs, B. Shapiro, and V. Nevistic,
“Nonlinear Games: Examples and Counterexamples,”
in Proceedings of 35th IEEE Conference on Decision and Control, Kobe, Japan, 1996,
vol. 4, pp. 3915–3920. doi: 10.1109/CDC.1996.577292.
"""
struct TwoDimensionalNonlinearOscillator <: AbstractEnv
end

function State(env::TwoDimensionalNonlinearOscillator)
    return function (x1, x2)
        ComponentArray(x1=x1, x2=x2)
    end
end

function Dynamics!(env::TwoDimensionalNonlinearOscillator)
    @Loggable function dynamics!(dx, x, p, t; u)
        @unpack x1, x2 = x
        @log state = x
        @log input = u
        dx.x1 = x2
        dx.x2 = (
                 -x1 * (π/2 + atan(5*x1))
                 - (5*x1^2)/(2*(1+25*x1^2))
                 + 4*x2
                 + 3*u[1]
                )
    end
end

function RunningCost(env::TwoDimensionalNonlinearOscillator)
    return function (x, u)
        @unpack x1, x2 = x
        r = x2^2 + u[1]^2
    end
end

function OptimalValue(env::TwoDimensionalNonlinearOscillator)
    return function (x)
        @unpack x1, x2 = x
        V_star = x1^2 * (π/2 + atan(5*x1)) + x2^2
    end
end

function OptimalControl(env::TwoDimensionalNonlinearOscillator)
    return function (x)
        @unpack x1, x2 = x
        u_star = [-3*x2]  # make it vector
    end
end


text:
"""
   compute the frequency amplitude spectra for a vector signal, if fhigh_end is not
specified, the default is nyquensit frequency.
"""
function amplitude_spectra(f::Vector{Tv}, dt;
                           nf=0, fhigh_end=0.0) where {Tv<:AbstractFloat, Ti<:Int64}

    # length of signal
    nt = length(f)

    # determine the number of zeros need to be padded
    if nf == 0 || nf < nt
       nf = nextpow(2, nt)
    end

    #  make sure nf is an even number
    if mod(nf, 2) != 0
       nf = nf + 1
    end

    # padding zeros
    f_pad = vcat(f, zeros(Tv, nf-nt))

    # frequency sampling interval
    fmax = 1.0 / dt
    df   = fmax/ nf

    # taking normalized transform
    F = 1.0 / sqrt(nf) * fft(f_pad)

    # only output amplitude spectra
    if fhigh_end == 0.0
       nw = floor(Int64, nf/2) + 1
    else
       nw = floor(Int64, fhigh_end/df) + 1
    end

    # make the frequency axis
    faxis = collect(0 : nw-1) * df

    # make amplitude spectrum
    amp= abs.(F[1:nw])

    return faxis, amp

end

"""
   compute the mean time-frequency spectrum when input is a matrix, the amplitude spectra
of each column is averaged.
"""
function amplitude_spectra(d::Matrix{Tv}, dt;
         nf=0, fhigh_end=0.0) where {Tv<:AbstractFloat, Ti<:Int64}

    # size of input
    (nt, ntrace) = size(d)

    # determine the number of zeros need to be padded
    if nf == 0 || nf < nt
       nf = nextpow(2, nt)
    end

    #  make sure nf is an even number
    if mod(nf, 2) != 0
       nf = nf + 1
    end

    # padding zeros to input
    d_pad = vcat(d, zeros(Tv, nf-nt, ntrace))

    # frequency sampling interval
    fmax = 1.0 / dt
    df   = fmax/ nf

    # taking normalized transform
    F = 1.0 / sqrt(nf) * fft(d_pad, 1)

    # only output amplitude spectra
    if fhigh_end == 0.0
       nw = floor(Int64, nf/2) + 1
    else
       nw = floor(Int64, fhigh_end/df) + 1
    end
    faxis = collect(0:nw-1) * df

    # smooth the amplitude spectra
    amp = abs.(F[1:nw, :])
    amp = sum(amp; dims=2) / ntrace

    return faxis, amp
end

# """
#    compute time derivative in frequency domain
# """
# function spectra_time_derivative(f::Vector{Tv}, dt) where {Tv<:AbstractFloat}
#
#     nt = length(f)
#
#     # padding zeros to input
#     nf = nextpow(2, nt)   # nf is always an even number
#     fp = vcat(f, zeros(Tv, nf-nt))
#
#     # frequency sampling interval
#     fmax = 1.0 / dt
#     df   = fmax/ nf
#
#     # transform to frequency domain
#     F = fft(fp)
#
#     # index of nyquensit frequency
#     iw_upper = floor(Int64, nf/2)
#
#     # multiplication
#     F[1] = zero(Complex{Tv})
#     for i = 2 : iw_upper
#         omega = 2 * pi * (i-1) * df
#         F[i]  = 1im * omega * F[i]
#
#         # conjugcy of real signal
#         F[nf-i+2] = conj(F[i])
#     end
#
#     # inverse fft, take the real part and truncate to original length
#     return real(ifft(F))[1:nt]
#
# end
#
# """
#    band pass filtering along the first dimension
# """
# function band_pass_filter(din::Array{Tv}, dt::Tv,
#          f1::Tv, f2::Tv, f3::Tv, f4::Tv) where {Tv<:AbstractFloat}
#
#     # get the data format
#     data_format = eltype(din)
#
#     # dimensions of input data
#     dims = collect(size(din))
#
#     # padding zeros
#     nt  = dims[1]
#     nf  = 4 * nextpow2(nt)
#     dims_pad = copy(dims)
#     dims_pad[1] = nf - nt
#     pad = zeros(data_format, dims_pad...)
#     dp  = vcat(din, pad)
#     dp  = fft(dp, 1)
#
#     # the index of Nyquensit frequency
#     nyqind = floor(Int64, nf/2)+1
#
#     # the index of corner frequency
#     i1 = floor(Int64, f1*dt*nf) + 1
#     i2 = floor(Int64, f2*dt*nf) + 1
#     i3 = floor(Int64, f3*dt*nf) + 1
#     i4 = floor(Int64, f4*dt*nf) + 1
#
#     # amplitude taper in frequency domain
#     up  = collect(1:(i2-i1)) / (i2-i1)
#     down= collect((i4-i3):-1:1) / (i4-i3)
#     aux = vcat(zeros(i1), up, ones(i3-i2), down, zeros(nyqind-i4))
#     aux2= flipdim(aux[2:nyqind-1], 1)
#
#     # complex band pass filter
#     c   = 0.
#     amp = vcat(aux, aux2);
#     pha = (pi/180.) * vcat(0., c*ones(nyqind-2), 0., -c*ones(nyqind-2))
#     bpfilter = amp .* exp.(im*pha)
#
#     # applying filter trace by trace
#     if length(dims) == 1
#        dout  = zeros(Complex{data_format}, nf)
#        dout  = bpfilter .* dp
#        dout  = ifft(dout, 1)
#        dfilt = real(dout[1:nt])
#
#     elseif length(dims) >= 2
#        nx   = prod(dims[2:end])
#        dout = zeros(Complex{data_format}, nf, nx)
#        dp   = reshape(dp, nf, nx)
#        for ix = 1: nx
#            dout[:, ix] = bpfilter .* dp[:, ix]
#        end
#        dout  = ifft(dout, 1)
#        dfilt = real(dout[1:nt, :])
#        dfilt = reshape(dfilt, dims...)
#     end
#
#     return dfilt
# end
#
# # upsampling of time series via sinc interpolation
# function sinc_resample1D(d::Vector{Tv}, order::Ti) where {Tv<:AbstractFloat, Ti<:Integer}
#
#     N    = length(d)
# 	  Npad = order*N - 1
#     nf   = 2*Npad
#     nw   = convert(Int,floor(nf/2)) + 1
# 	  dd   = zeros(eltype(d), nf)
#     dd[1:order:Npad] = d[1:1:end]
# 	  D = fft(dd)
#     nyq = convert(Int,floor(nw/order)) + 1
#
# 	  for iw = nyq : nw
# 		    D[iw] *= 0
# 	  end
#
# 	  # symmetries
# 	  for iw=nw+1:nf
# 		    D[iw] = conj(D[nf-iw+2])
# 	  end
#
#     dd = real(ifft(D,1))
# 	  return dd[1:Npad]*order
#
# end


text:
@testset "Approximators" begin

    @testset "TabularApproximator" begin
        A = TabularVApproximator(; n_state = 2, opt = InvDecay(1.0))

        @test A(1) == 0.0
        @test A(2) == 0.0

        update!(A, 2 => A(2) - 3.0)
        @test A(2) == 1.5
        update!(A, 2 => A(2) - 6.0)
        @test A(2) == 3.0
    end

    @testset "NeuralNetworkApproximator" begin
        NN = NeuralNetworkApproximator(; model = Dense(2, 3), optimizer = Descent())

        q_values = NN(rand(2))
        @test size(q_values) == (3,)

        gs = gradient(params(NN)) do
            sum(NN(rand(2, 5)))
        end

        old_params = deepcopy(collect(params(NN).params))
        update!(NN, gs)
        new_params = collect(params(NN).params)

        @test old_params != new_params
    end

    @testset "ActorCritic" begin
        ac_cpu = ActorCritic(
            actor = NeuralNetworkApproximator(model = Dense(3, 2)),
            critic = NeuralNetworkApproximator(model = Dense(3, 1)),
        )

        ac = ac_cpu |> gpu

        # make sure optimizer is not changed
        @test ac_cpu.optimizer === ac.optimizer

        D = ac.actor.model |> gpu |> device
        @test D === device(ac) === device(ac.actor) == device(ac.critic)

        A = send_to_device(D, rand(3))
        ac.actor(A)
        ac.critic(A)
    end
end


text:
#=

This file implements rt_kernel!, which performs the core RT routines (elemental, doubling, interaction)

=#

# Perform the Core RT routines (elemental, doubling, interaction)
function rt_kernel!(pol_type, SFI, added_layer, composite_layer, computed_layer_properties, m, quad_points, I_static, architecture, qp_μN, iz)

    @unpack τ_λ, ϖ_λ, τ, ϖ, Z⁺⁺, Z⁻⁺, dτ_max, dτ, ndoubl, dτ_λ, expk, scatter, τ_sum, scattering_interface = computed_layer_properties

    # If there is scattering, perform the elemental and doubling steps
    if scatter
        @timeit "elemental" elemental!(pol_type, SFI, τ_sum, dτ_λ, dτ, ϖ_λ, ϖ, Z⁺⁺, Z⁻⁺, m, ndoubl, scatter, quad_points,  added_layer,  I_static, architecture)
        @timeit "doubling"   doubling!(pol_type, SFI, expk, ndoubl, added_layer, I_static, architecture)
    else # This might not work yet on GPU!
        # If not, there is no reflectance. Assign r/t appropriately
        added_layer.r⁻⁺[:] .= 0;
        added_layer.r⁺⁻[:] .= 0;
        added_layer.J₀⁻[:] .= 0;
        temp = Array(exp.(-τ_λ./qp_μN'))
        #added_layer.t⁺⁺, added_layer.t⁻⁻ = (Diagonal(exp(-τ_λ / qp_μN)), Diagonal(exp(-τ_λ / qp_μN)))
        for iλ = 1:length(τ_λ)
            added_layer.t⁺⁺[:,:,iλ] = Diagonal(temp[iλ,:]);
            added_layer.t⁻⁻[:,:,iλ] = Diagonal(temp[iλ,:]);
        end
    end

    # @assert !any(isnan.(added_layer.t⁺⁺))

    # If this TOA, just copy the added layer into the composite layer
    if (iz == 1)
        composite_layer.T⁺⁺[:], composite_layer.T⁻⁻[:] = (added_layer.t⁺⁺, added_layer.t⁻⁻)
        composite_layer.R⁻⁺[:], composite_layer.R⁺⁻[:] = (added_layer.r⁻⁺, added_layer.r⁺⁻)
        composite_layer.J₀⁺[:], composite_layer.J₀⁻[:] = (added_layer.J₀⁺, added_layer.J₀⁻ )

    # If this is not the TOA, perform the interaction step
    else
        @timeit "interaction" interaction!(scattering_interface, SFI, composite_layer, added_layer, I_static)
    end
end


text:
__precompile__(true)

module Display

export @display
macro display(args...)
	blk = Expr(:block)
	push!(blk.args, :(io = IOBuffer()))
	for arg in args
		push!(blk.args, :(show(io, "text/plain", $(esc(arg)))))
	end
	push!(blk.args, :(println(String(take!(io)))))
	blk
end

end


text:
using DataStructures
using Plots; pyplot(fmt="png")
PyPlot.matplotlib.rc("mathtext", fontset="cm");
PyPlot.matplotlib.rc("text", usetex=true);
PyPlot.matplotlib.rc("font", family="serif", size=12);
using LaTeXStrings
using LsqFit
using Statistics
using JLD

raw_data_dir  = "/home/user/Downloads/raw_data";
data_dir = "/home/user/thesis/data/scaling";
t_0_dir   = joinpath(raw_data_dir, "cluster_size_distribution/t_0");
t_1_dir   = joinpath(raw_data_dir, "cluster_size_distribution/t_1");
t_0_files = joinpath.(t_0_dir, readdir(t_0_dir));
t_1_files = joinpath.(t_1_dir, readdir(t_1_dir));

savepath = "/home/user/thesis/latex/images"
save_bool = true
dpi = 144
if save_bool
	dpi = 300
end

# %%

file_dict = Dict(
	"t_0" => joinpath.(t_0_dir, readdir(t_0_dir)),
	"t_1" => joinpath.(t_1_dir, readdir(t_1_dir))
);

function distribution_files_to_dict(file_dict, Ns)
	distributions = Dict(
		t => Dict(N => Dict() for N in Ns) for t in keys(file_dict)
	)
	N_counts = Dict(
		t => Dict() for t in keys(file_dict)
	)

	for (t, file_list) in file_dict
		N_counts[t] = Dict(
			N => sum([1 for file_name in file_list if occursin(string("_", Int(log2(N)), "_"), file_name)])
			for N in Ns
		)

		for file in file_list
			a             = maximum(findfirst("acceptance_", file)) + 1
			b             = minimum(findfirst("_t_", file)) - 1
			N             = 2^parse(Int, file[a:b])
			a             = maximum(findfirst("_seed_", file)) + 1
			b             = minimum(findfirst(".csv", file)) - 1
			seed          = parse(Int, file[a:b])
			distribution  = Dict()

			for (i, line) in enumerate(eachline(file))
				i ∈ (1, 2, 3) ? continue : nothing

				s, N_s = parse.(Int, split(line, ","))
				distribution[s] = N_s / N
			end

			distributions[t][N][seed] = distribution
		end
	end

	distribution_avgs = Dict(
		t => Dict(N => DefaultDict{Int, Float64}(0.) for N in Ns) for t in keys(file_dict)
	)
	distribution_arrays = Dict(
		t => Dict(N => DefaultDict{Int, Array}([]) for N in Ns) for t in keys(file_dict)
	)
	for (t, N_dicts) in distributions
		for (N, seed_dicts) in N_dicts
			for (seed, distribution) in seed_dicts
				for (s, n_s) in distribution
					distribution_avgs[t][N][s] += n_s / N_counts[t][N]
					push!(distribution_arrays[t][N][s], n_s)
				end
			end
		end
	end

	# # CHECK
	# distribution_avgs_v2 = Dict(
	# 	t => Dict(N => DefaultDict{Int, Float64}(0.) for N in Ns) for t in keys(file_dict)
	# )
	# for (t, N_dicts) in distributions
	#
	# end

	distribution_avgs_ = Dict()
	for (t, N_dicts) in distribution_avgs
		t_dict = Dict()
		for (N, distribution) in N_dicts
			t_dict[N] = Dict(distribution)
		end
		distribution_avgs_[t] = t_dict
	end

	return distributions, distribution_avgs_
end

function distribution_file_to_dict(file)
	N = 0
	distribution = Dict()

	for (i, ln) in enumerate(eachline(file))
		if i == 1
			N = parse(Int, ln[7:end])
			continue
		end
		if i == 2
			continue
		end
		s, N_s = parse.(Int, split(ln, ","))
		distribution[s] = N_s / N
	end

	return (distribution, N)
end

# %% Save and load data

Ns = 2 .^ collect(15:24);
distributions, distribution_avgs = distribution_files_to_dict(file_dict, Ns);

save_file = joinpath(data_dir, "cluster_size_distribution_avgs.jld")
save(save_file, "distribution_avgs", distribution_avgs)

max_stars_repo_name: cameronperot/explosive-percolation

max_issues_repo_name: cameronperot/explosive-percolation

max_forks_repo_name: cameronperot/explosive-percolation


text:
# This file was generated, do not modify it. # hide
@load KNNClassifier pkg="NearestNeighbors"
@load MultinomialClassifier pkg="MLJLinearModels";

KnnPipe = @pipeline(Standardizer(), KNNClassifier())
MnPipe = @pipeline(Standardizer(), MultinomialClassifier());


text:
"""
Distributed AdaOja's algorithm for streaming/online PCA with MPI

Author: Wayne Wang
Email: wayneyw@umich.edu
Last modified: 08/26/2021
"""

function oja_mpi(node_id::Int, num_nodes::Int, comm, input_dir::AbstractString;
    k::Int = 5)
    """
    AdaOja's algorithm with MPI
    """
    # initialization
    Xt = data_par_loader(node_id, input_dir) #data stream of dimension d_b × N_t
    d_par, N_t = size(Xt)
    V_par, grad_par, α = init(d_par, k)
    # save intial estimate from each partition
    save(string(input_dir, "Vhat_init_par_", node_id + 1, ".jld2"), Dict("Vhat_par" => V_par))

    # run streaming algorithm on each rank
    Xt = data_par_loader(node_id, input_dir)
    X_par = zeros((d_par, 1))

    # start a timer
    MPI.Barrier(comm) #make sure every process reaches this point
    tstart = MPI.Wtime()

    for i = 1:N_t
        X_par[:, 1] .= Xt[:, i] #making sure the (d_par,1) deos not collapse to (d_par,)
        # compute gradient
        ojagrad_mpi!(grad_par, X_par, V_par, node_id, num_nodes, comm)
        # update learning rates and estimate
        ojaupdate_mpi!(V_par, α, grad_par, node_id, num_nodes, comm)
        # QR orthogonalization
        cholesky_tsqr!(V_par, node_id, num_nodes, comm; k = k)
    end

    # stop the timer
    MPI.Barrier(comm) #make sure every process reaches this point
    tend = MPI.Wtime() - tstart
    # if node_id == 0 #output timing result only on the master node
    #     # print(string(tend, "\n"))
    #     open("time_oja_mpi_scaling.txt", "a") do f
    #         write(f, string(num_nodes, "\t", tend, "\n"))
    #     end
    # end

    # save final estimate from each partition
    # save(string(input_dir, "Vhat_par_", node_id + 1, ".jld2"), Dict("Vhat_par" => V_par))

    return nothing
end

function ojagrad_mpi!(grad_par::AbstractArray, X_par::AbstractArray, V_par::AbstractArray,
    node_id::Int, num_nodes::Int, comm;
    N::Int = 1, master::Int = 0,
    tag_from_master::Int = 1,
    tag_from_worker::Int = 2)
    """
    Compute X × (X^T × V) in the Oja's update using MPI
    """
    # each node (including master) get a partition of
    # X (data) and V (current estimate)
    @assert size(X_par, 1) == size(V_par, 1)

    # compute X^T × V partitions and collect via allreduce, i.e.,
    # X^T_V = X^T_V_par_1 + X^T_V_par_2 + ⋯
    XT_V_par = copy(X_par') * V_par
    XT_V = MPI.Allreduce(XT_V_par, +, comm)

    # compute X × (X^T × V)
    mul!(grad_par, X_par, XT_V, 1.0/N, 0.0)

    return nothing
end

function ojaupdate_mpi!(V_par::AbstractArray, α::AbstractArray,
    grad_par::AbstractArray, node_id::Int, num_nodes::Int, comm;
    master::Int = 0,
    tag_from_master::Int = 1,
    tag_from_worker::Int = 2)
    """
    Oja × Adagrad
    """
    # compute local sum of squared grad and reduce
    grad_norm_sq = MPI.Allreduce(sum(abs2, grad_par, dims = 1)[:], +, comm) # ∈ k × 1
    α .+= grad_norm_sq # α ∈ k × 1

    # update partitions of the current estimate
    V_par .+= grad_par ./ reshape(sqrt.(α), (1, size(α, 1)))

    return nothing
end

function direct_tsqr!(V_par::AbstractArray, node_id::Int, num_nodes::Int, comm;
    k::Int = 5, master::Int = 0, tag_from_master::Int = 1, tag_from_worker::Int = 2)
    """
    Implementation of the DirectTSQR algorithm (Benson, Cleich, and Demmel, 2013)
    """
    # local qr on each partition
    QR1 = qr!(V_par)
    Q1_par = Array(QR1.Q)
    R1_par = Array(QR1.R)

    # allocate intermediate data
    Q2_par = zeros((k, k))

    # master tasks
    if node_id == master
        # gather all R1 partitions on master
        R1 = zeros((k * num_nodes, k))
        R1[1:k, :] .= R1_par

        rreqs = MPI.Request[]
        for i = 1:(num_nodes - 1)
            rreq = MPI.Irecv!(view(R1, (k * i + 1):k * (i + 1), :), i, tag_from_worker, comm)
            push!(rreqs, rreq)
        end

        # perform QR for R1 on master node
        MPI.Waitall!(rreqs) #make sure received the full R1 matrix
        QR2 = qr!(R1)
        Q2 = Array(QR2.Q)

        # send Q2 partitions to workers
        Q2_par .= Q2[1:k, :]
        for i = 1:(num_nodes - 1)
            MPI.Isend(view(Q2, (k * i + 1):k * (i + 1), :), i, tag_from_master, comm)
        end
    end

    # worker tasks
    if node_id != master
        # send R1 partitions to the master
        MPI.Send(R1_par, master, tag_from_worker, comm)

        # receive Q2 partitions from master
        MPI.Recv!(Q2_par, master, tag_from_master, comm)
    end

    # form final partitions of the Q matrix inplace to V_par
    V_par .= Q1_par * Q2_par

    return nothing
end

function direct_tsqr_collective!(V_par::AbstractArray, node_id::Int, num_nodes::Int, comm;
    k::Int = 5, master::Int = 0, tag_from_master::Int = 1, tag_from_worker::Int = 2)
    """
    Implementation of the DirectTSQR algorithm (Benson, Cleich, and Demmel, 2013)
     * this version uses MPI collective communications only *
    """
    # local qr on each partition
    QR1 = qr!(V_par)
    Q1_par = Array(QR1.Q)
    R1_par = Array(QR1.R) #∈ k × k

    # allocate intermediate data
    Q2_par_1d = zeros(k^2)
    Q2_par = zeros(k, k)

    # gather all R1 partitions on master
    R1_1d = MPI.Gather(R1_par, master, comm) #∈ num_nodes*k*k

    if node_id == master
        ## convert to 2D Array
        R1 = zeros(num_nodes * k, k) #∈ num_nodes*k × k
        for i = 1:num_nodes
            R1[(i - 1) * k + 1:i * k, :] .= reshape(R1_1d[(i - 1) * k^2 + 1:(i * k^2)], (k, k))
        end

        ## perform QR for R1 on master node
        QR2 = qr!(R1)
        Q2 = copy(Array(QR2.Q)')[:]
    else
        Q2 = Float64[]
    end

    # scatter partitions of Q2
    MPI.Scatter!(MPI.UBuffer(Q2, k^2), Q2_par_1d, master, comm)

    # convert Q2_par to 2D Arrays
    for _ = 1:num_nodes
        Q2_par .= copy(reshape(Q2_par_1d, (k, k))')
    end

    # form final partitions of the Q matrix inplace to V_par
    V_par .= Q1_par * Q2_par

    return nothing
end

function cholesky_tsqr!(V_par::AbstractArray, node_id::Int, num_nodes::Int, comm;
    k::Int = 5, master::Int = 0, tag_from_master::Int = 1, tag_from_worker::Int = 2)
    """
    Implementation of the DirectTSQR algorithm (Benson, Cleich, and Demmel, 2013)
    """
    # compute V^T × V partitions and collect via reduce, i.e.,
    VT_V_par = copy(V_par') * V_par
    VT_V = MPI.Reduce(VT_V_par, +, master, comm) #∈ k × k

    # master tasks
    if node_id == master
        ## cholesky decomposition
        ### could add a small diagonal matrix to VT_V to ensure positiveness!
        chol = cholesky!(VT_V)
        ## compute R inverse
        Rinv = inv(chol.U)
    else
        Rinv = nothing
    end

    # broadcast to all ranks
    Rinv = MPI.bcast(Rinv, master, comm)

    # form final partitions of the Q matrix inplace to V_par
    V_par .= V_par * Rinv

    return nothing
end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
using CLIMA
using CLIMA.Mesh.Topologies: StackedCubedSphereTopology, cubedshellwarp, grid1d
using CLIMA.Mesh.Grids: DiscontinuousSpectralElementGrid
using CLIMA.Mesh.Filters
using CLIMA.DGmethods: DGModel, init_ode_state, VerticalDirection
using CLIMA.DGmethods.NumericalFluxes: Rusanov, CentralNumericalFluxGradient,
                                       CentralNumericalFluxDiffusive
using CLIMA.ODESolvers
using CLIMA.GeneralizedMinimalResidualSolver
using CLIMA.ColumnwiseLUSolver: ManyColumnLU
using CLIMA.VTK: writevtk, writepvtu
using CLIMA.GenericCallbacks: EveryXWallTimeSeconds, EveryXSimulationSteps
using CLIMA.PlanetParameters: planet_radius, day
using CLIMA.MoistThermodynamics: air_density, soundspeed_air, internal_energy
using CLIMA.Atmos: AtmosModel, SphericalOrientation,
                   DryModel, NoPrecipitation, NoRadiation, NoFluxBC,
                   ConstantViscosityWithDivergence,
                   vars_state, vars_aux,
                   Gravity, HydrostaticState, IsothermalProfile,
                   AtmosAcousticGravityLinearModel, AtmosLESConfiguration
using CLIMA.VariableTemplates: flattenednames

using MPI, Logging, StaticArrays, LinearAlgebra, Printf, Dates, Test

const output_vtk = false

function main()
  CLIMA.init()
  ArrayType = CLIMA.array_type()

  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = Dict("DEBUG" => Logging.Debug,
                  "WARN"  => Logging.Warn,
                  "ERROR" => Logging.Error,
                  "INFO"  => Logging.Info)[ll]

  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  polynomialorder = 5
  numelem_horz = 10
  numelem_vert = 5

  timeend = 60 * 60
  # timeend = 33 * 60 * 60 # Full simulation
  outputtime = 60 * 60

  expected_result = Dict()
  expected_result[Float32] = 9.5064378310656000e+13
  expected_result[Float64] = 9.5073452847081828e+13

  for FT in (Float32, Float64)
    result = run(mpicomm, polynomialorder, numelem_horz, numelem_vert,
                 timeend, outputtime, ArrayType, FT)
    @test result ≈ expected_result[FT]
  end
end

function run(mpicomm, polynomialorder, numelem_horz, numelem_vert,
             timeend, outputtime, ArrayType, FT)

  setup = AcousticWaveSetup{FT}()

  vert_range = grid1d(FT(planet_radius), FT(planet_radius + setup.domain_height), nelem = numelem_vert)
  topology = StackedCubedSphereTopology(mpicomm, numelem_horz, vert_range)

  grid = DiscontinuousSpectralElementGrid(topology,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = polynomialorder,
                                          meshwarp = cubedshellwarp)

  model = AtmosModel{FT}(AtmosLESConfiguration;
                         orientation=SphericalOrientation(),
                           ref_state=HydrostaticState(IsothermalProfile(setup.T_ref), FT(0)),
                          turbulence=ConstantViscosityWithDivergence(FT(0)),
                            moisture=DryModel(),
                              source=Gravity(),
                          init_state=setup)
  linearmodel = AtmosAcousticGravityLinearModel(model)

  dg = DGModel(model, grid, Rusanov(),
               CentralNumericalFluxDiffusive(), CentralNumericalFluxGradient())

  lineardg = DGModel(linearmodel, grid, Rusanov(),
                     CentralNumericalFluxDiffusive(),
                     CentralNumericalFluxGradient();
                     direction=VerticalDirection(),
                     auxstate=dg.auxstate)

  # determine the time step
  element_size = (setup.domain_height / numelem_vert)
  acoustic_speed = soundspeed_air(FT(setup.T_ref))
  dt_factor = 445
  dt = dt_factor * element_size / acoustic_speed / polynomialorder ^ 2
  # Adjust the time step so we exactly hit 1 hour for VTK output
  dt = 60 * 60 / ceil(60 * 60 / dt)
  nsteps = ceil(Int, timeend / dt)

  Q = init_ode_state(dg, FT(0))

  linearsolver = ManyColumnLU()

  odesolver = ARK2GiraldoKellyConstantinescu(dg, lineardg, linearsolver, Q;
                                             dt = dt, t0 = 0, split_nonlinear_linear=false)

  filterorder = 18
  filter = ExponentialFilter(grid, 0, filterorder)
  cbfilter = EveryXSimulationSteps(1) do
    Filters.apply!(Q, 1:size(Q, 2), grid, filter, VerticalDirection())
    nothing
  end

  eng0 = norm(Q)
  @info @sprintf """Starting
                    ArrayType       = %s
                    FT              = %s
                    polynomialorder = %d
                    numelem_horz    = %d
                    numelem_vert    = %d
                    dt              = %.16e
                    norm(Q₀)        = %.16e
                    """ "$ArrayType" "$FT" polynomialorder numelem_horz numelem_vert dt eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      runtime = Dates.format(convert(DateTime, now() - starttime[]), dateformat"HH:MM:SS")
      @info @sprintf """Update
                        simtime = %.16e
                        runtime = %s
                        norm(Q) = %.16e
                        """ gettime(odesolver) runtime energy
    end
  end
  callbacks = (cbinfo, cbfilter)

  if output_vtk
    # create vtk dir
    vtkdir = "vtk_acousticwave" *
      "_poly$(polynomialorder)_horz$(numelem_horz)_vert$(numelem_vert)" *
      "_dt$(dt_factor)x_$(ArrayType)_$(FT)"
    mkpath(vtkdir)

    vtkstep = 0
    # output initial step
    do_output(mpicomm, vtkdir, vtkstep, dg, Q, model)

    # setup the output callback
    cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
      vtkstep += 1
      Qe = init_ode_state(dg, gettime(odesolver))
      do_output(mpicomm, vtkdir, vtkstep, dg, Q, model)
    end
    callbacks = (callbacks..., cbvtk)
  end

  solve!(Q, odesolver; numberofsteps=nsteps, adjustfinalstep=false, callbacks=callbacks)

  # final statistics
  engf = norm(Q)
  @info @sprintf """Finished
  norm(Q)                 = %.16e
  norm(Q) / norm(Q₀)      = %.16e
  norm(Q) - norm(Q₀)      = %.16e
  """ engf engf/eng0 engf-eng0
  engf
end

Base.@kwdef struct AcousticWaveSetup{FT}
  domain_height::FT = 10e3
  T_ref::FT = 300
  α::FT = 3
  γ::FT = 100
  nv::Int = 1
end

function (setup::AcousticWaveSetup)(bl, state, aux, coords, t)
  # callable to set initial conditions
  FT = eltype(state)

  r = norm(coords, 2)
  @inbounds λ = atan(coords[2], coords[1])
  @inbounds φ = asin(coords[3] / r)
  h = r - FT(planet_radius)

  β = min(FT(1), setup.α * acos(cos(φ) * cos(λ)))
  f = (1 + cos(FT(π) * β)) / 2
  g = sin(setup.nv * FT(π) * h / setup.domain_height)
  Δp = setup.γ * f * g
  p = aux.ref_state.p + Δp

  state.ρ = air_density(setup.T_ref, p)
  state.ρu = SVector{3, FT}(0, 0, 0)
  state.ρe = state.ρ * (internal_energy(setup.T_ref) + aux.orientation.Φ)
  nothing
end

function do_output(mpicomm, vtkdir, vtkstep, dg, Q, model, testname = "acousticwave")
  ## name of the file that this MPI rank will write
  filename = @sprintf("%s/%s_mpirank%04d_step%04d",
                      vtkdir, testname, MPI.Comm_rank(mpicomm), vtkstep)

  statenames = flattenednames(vars_state(model, eltype(Q)))
  auxnames = flattenednames(vars_aux(model, eltype(Q)))
  writevtk(filename, Q, dg, statenames, dg.auxstate, auxnames)

  ## Generate the pvtu file for these vtk files
  if MPI.Comm_rank(mpicomm) == 0
    ## name of the pvtu file
    pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

    ## name of each of the ranks vtk files
    prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
      @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
    end

    writepvtu(pvtuprefix, prefixes, (statenames..., auxnames...))

    @info "Done writing VTK: $pvtuprefix"
  end
end

main()

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
export DPM, VCM, HDP

# define models

mutable struct DPM <: ModelType
  H::ConjugatePostDistribution
  α::Float64

  DPM(H::ConjugatePostDistribution; α = 1.0) = new(H, α)
end

mutable struct VCM <: ModelType
  α::Float64

  VCM(;α = 1.0) = new(α)
end

mutable struct HDP <: ModelType
  H::ConjugatePostDistribution
  α::Float64
  γ::Float64

  HDP(H::ConjugatePostDistribution; α = 0.1, γ = 5.0) = new(H, α, γ)
end