text:
const Parent = Tuple{ENode,Int64} # parent enodes and eclasses
mutable struct EClassData
    id::Int64
    nodes::Vector{ENode}
    parents::Vector{Parent}
end

# Interface for indexing EClassData
Base.getindex(a::EClassData, i) = a.nodes[i]
Base.setindex!(a::EClassData, v, i) = setindex!(a.nodes, v, i)
Base.firstindex(a::EClassData) = firstindex(a.nodes)
Base.lastindex(a::EClassData) = lastindex(a.nodes)

# Interface for iterating EClassData
Base.iterate(a::EClassData) = iterate(a.nodes)
Base.iterate(a::EClassData, state) = iterate(a.nodes, state)

# Showing
Base.show(io::IO, a::EClassData) = Base.show(io, a.nodes)

EClassData(id::Int64) = EClassData(id, [], Parent[])

function addparent!(a::EClassData, parent::Parent)
    if parent ∉ a.parents
        push!(a.parents, parent)
    end
end

function Base.union(from::EClassData, to::EClassData)
    EClassData(to.id, from.nodes ∪ to.nodes, from.parents ∪ to.parents)
end

function Base.union!(to::EClassData, from::EClassData)
    union!(to.nodes, from.nodes)
    union!(to.parents, from.parents)
    return to
end


text:
using PlotlyJS

function logistic(r,x)
    r*x*(1.0-x)
end

function logisticmap(;x0=0.21,r=4,n=10)
    results = [(x0,logistic(r,x0))];
    for i in 1:(n-1)
      xold = results[end][2]
      xnew = logistic(r,xold)
      append!(results,[(xold,xold)])
      append!(results,[(xold,xnew)])
  end
  results
end

simres = logisticmap(x0=0.2,r=4,n=20)

simres2 = logisticmap(x0=0.2001,r=4,n=20)

xvals = [x for (x, y) in simres]
yvals = [y for (x, y) in simres]

xvals2 = [x for (x, y) in simres2]
yvals2 = [y for (x, y) in simres2]

trace = scatter(;x=xvals, y=yvals, mode="lines", name = "x0=0.2")

trace2 = scatter(;x=xvals2, y=yvals2, mode="lines", name = "x0=0.2001")

plot([ trace, trace2])

savefig("../../images/doublelines.png")

ext: jl


text:
"""
	read_layerResponse(filein)
Read `layerResponse` results. See GravityEffect.jl package for details

**Input**
* filein: input file
* par: parameter to be read: "results","zones","layers","exclude","nanheight","outfile","def_density"

**Example**
```
fileinput = "test/input/layerResponse.txt";
response = read_layerResponse(fileinput,"results");
layers = read_layerResponse(fileinput,"layers");
```
"""
function read_layerResponse(filein::String,par::String="results")
	if par=="results"
		d = readdlm(filein,comments=true,comment_char='%');
		out = DataFrame(layer=d[:,1],start=d[:,2],stop=d[:,3],total=d[:,4]);
		for i in 1:size(d,2)-4
			out[Symbol("zone"*string(i))] = d[:,i+4]
		end
		return out
	else
		row = ""
		open(filein,"r") do fid
			while !occursin(lowercase(par),lowercase(row))
				row = eof(fid) ? par : readline(fid)
			end
		end
		return occursin("->",row) ? eval(Meta.parse(split(row,"->")[2])) : []
	end
end

"""
	write_layerResponse(sensor,layers,zones,exclude,nanheight,outfile,def_density,outdata)
Write `layerResponse` results. See GravityEffect.jl package for details

"""
function write_layerResponse(sensor,layers,zones,exclude,nanheight,outfile,def_density,outdata)
	open(outfile,"w") do fid
		@printf(fid,"%% layerResponse.jl SETTINGS:\n");
		# Write settings
		println(fid,"% sensor->$sensor");
		println(fid,"% layers->$layers");
		println(fid,"% zones->$zones");
		println(fid,"% exclude->$exclude");
		println(fid,"% nanheight->$nanheight");
		println(fid,"% outfile->$outfile");
		println(fid,"% def_density->$def_density");
		@printf(fid,"%% Computation date: %s\n",Dates.format(now(),"dd/mm/yyyy HH:MM:SS"))
		@printf(fid,"%% RESULTS (gravity effect in nm/s^2, depth in m): \n");
		@printf(fid,"%%  Nr start  stop    total ");
		for i in 1:length(zones[:radius])
			@printf(fid," %8s%i","zone",i);
		end
		@printf(fid,"\n");
		for i in 1:size(outdata,1)
			@printf(fid,"%4i %6.3f %6.3f %9.5f",i,
					outdata[:start][i],outdata[:stop][i],outdata[:total][i]*1e+9);
			for j in 1:length(zones[:radius])
				@printf(fid," %9.5f",outdata[Symbol("zone"*string(j))][i]*1e+9);
			end
			@printf(fid,"\n");
		end
	end
end


text:
using GraphViz
fname = "graph"
dot = open(readall,"$fname.dot")
gviz = Graph(dot)
GraphViz.layout!(gviz,engine="dot")
open("$fname.svg","w") do f
       GraphViz.writemime(f, MIME"image/svg+xml"(),gviz)
end #do
open("$fname.png","w") do f
       GraphViz.writemime(f, MIME"image/png"(),gviz)
end #do


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ImplicitFreeSurface
using Statistics
using IterativeSolvers
using LinearAlgebra, SparseArrays
using Oceananigans.Solvers: constructors, unpack_constructors

function geostrophic_adjustment_simulation(free_surface, topology)

    Lh = 100kilometers
    Lz = 400meters

    grid = RectilinearGrid(size = (64, 3, 1),
                                  x = (0, Lh), y = (0, Lh), z = (-Lz, 0),
                                  topology = topology)

    coriolis = FPlane(f=1e-4)

    model = HydrostaticFreeSurfaceModel(grid = grid,
                                        coriolis = coriolis,
                                        free_surface = free_surface)

    gaussian(x, L) = exp(-x^2 / 2L^2)

    U = 0.1 # geostrophic velocity
    L = grid.Lx / 40 # gaussian width
    x₀ = grid.Lx / 4 # gaussian center

    vᵍ(x, y, z) = - U * (x - x₀) / L * gaussian(x - x₀, L)

    g = model.free_surface.gravitational_acceleration

    η₀ = coriolis.f * U * L / g # geostrohpic free surface amplitude

    ηᵍ(x) = η₀ * gaussian(x - x₀, L)

    ηⁱ(x, y) = 2 * ηᵍ(x)

    set!(model, v=vᵍ, η=ηⁱ)
    gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed
    wave_propagation_time_scale = model.grid.Δxᶜᵃᵃ / gravity_wave_speed
    simulation = Simulation(model, Δt=2wave_propagation_time_scale, stop_iteration=10000)

    return simulation
end

function run_and_analyze(simulation)
    η = simulation.model.free_surface.η
    u, v, w = simulation.model.velocities
    Δt = simulation.Δt

    ηx = Field(∂x(η))
    compute!(ηx)

    u₀  = interior(u)[:, 1, 1]
    v₀  = interior(v)[:, 1, 1]
    η₀  = interior(η)[:, 1, 1]
    ηx₀ = interior(ηx)[:, 1, 1]

    solver_method = string(simulation.mode.free_surface.solver_method)

    simulation.output_writers[:fields] = JLD2OutputWriter(simulation.model, (η, ηx, u, v, w),
                                                          schedule = TimeInterval(Δt),
                                                          prefix = "solution_$(solver_method)")

    run!(simulation)

    compute!(ηx)

    u₁ = interior(u)[:, 1, 1]
    v₁ = interior(v)[:, 1, 1]
    η₁ = interior(η)[:, 1, 1]
    ηx₁ = interior(ηx)[:, 1, 1]

    @show mean(η₀)
    @show mean(η₁)

    Δη = η₁ .- η₀

    return (; η₀, η₁, Δη, ηx₀, ηx₁, u₀, u₁, v₀, v₁)
end

# fft_based_free_surface = ImplicitFreeSurface()
pcg_free_surface = ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient);
matrix_free_surface = ImplicitFreeSurface(solver_method=:MatrixIterativeSolver);

topology_types = [(Bounded, Periodic, Bounded), (Periodic, Periodic, Bounded)]

free_surfaces = [pcg_free_surface, matrix_free_surface];
simulations = [geostrophic_adjustment_simulation(free_surface, topology_type) for free_surface in free_surfaces, topology_type in topology_types];
data = [run_and_analyze(sim) for sim in simulations];

using GLMakie
using JLD2

file1 = jldopen("solution_PreconditionedConjugateGradient.jld2")
file2 = jldopen("solution_MatrixIterativeSolver.jld2")

grid = file1["serialized/grid"]

x  = grid.xᶜᵃᵃ[1:grid.Nx]
xf = grid.xᶠᵃᵃ[1:grid.Nx+1]
y  = grid.yᵃᶜᵃ[1:grid.Ny]

iterations = parse.(Int, keys(file1["timeseries/t"]))
iterations = iterations[1:200]

iter = Node(0)

mid = Int(floor(grid.Ny/2))
η0 = file1["timeseries/1/0"][:, mid, 1]
η1 = @lift(Array(file1["timeseries/1/" * string($iter)])[:, mid, 1])
η2 = @lift(Array(file2["timeseries/1/" * string($iter)])[:, mid, 1])
u1 = @lift(Array(file1["timeseries/3/" * string($iter)])[:, mid, 1])
u2 = @lift(Array(file2["timeseries/3/" * string($iter)])[:, mid, 1])

fig = Figure(resolution=(1000, 500))
plot(fig[1,1] , x, η0, color = :green)
plot!(fig[1,1], x, η1, color = :red)
plot!(fig[1,1], x, η2, color = :blue)
plot(fig[1,2], xf, u1, color = :red)
plot!(fig[1,2],xf, u2, color = :blue)
ylims!(-5e-5, 5e-5)
GLMakie.record(fig, "free_surface_bounded.mp4", iterations, framerate=12) do i
    @info "Plotting iteration $i of $(iterations[end])..."
    iter[] = i
end

mid = Int(floor(grid.Ny/2))
η3 = @lift(Array(file3["timeseries/1/" * string($iter)])[:, mid, 1])
η4 = @lift(Array(file4["timeseries/1/" * string($iter)])[:, mid, 1])
u3 = @lift(Array(file3["timeseries/3/" * string($iter)])[:, mid, 1])
u4 = @lift(Array(file4["timeseries/3/" * string($iter)])[:, mid, 1])

fig = Figure(resolution=(1000, 500))
plot(fig[1,1] , x, η0, color = :green)
plot!(fig[1,1], x, η3, color = :red)
plot!(fig[1,1], x, η4, color = :blue)
plot(fig[1,2],  x, u3, color = :red )
plot!(fig[1,2], x, u4, color = :blue)
ylims!(-5e-5, 5e-5)
GLMakie.record(fig, "free_surface_periodic.mp4", iterations, framerate=12) do i
    @info "Plotting iteration $i of $(iterations[end])..."
    iter[] = i
end


text:
module HMCUtilities

export free,
    constrain,
    constrain_with_pushlogpdf,
    constrain_with_pushlogpdf_grad

using ForwardDiff
using DiffResults
using Zygote
using Einsum

include("constraint.jl")
include("advanced_hmc.jl")

end # module


text:
# run this file with julia -p 4 to use 4 cores in the comparison
using JLD2, Distributed
@everywhere begin
using ReinforcementLearningEnvironmentClassicControl, Flux,
ReinforcementLearning
getenv() = (CartPole(), 4, 2)
# getenv() = (MountainCar(maxsteps = 10^4), 2, 3)
function setup(learner, env, preprocessor = NoPreprocessor())
    cb = EvaluateGreedy(callback = EvaluationPerEpisode(TotalReward(),
                                                        returnmean = true),
                        stoppingcriterion = ConstantNumberEpisodes(200),
                        every = Episode(200))
    RLSetup(learner, env, ConstantNumberEpisodes(2000),
            callbacks = [cb], preprocessor = preprocessor)
end
function acpg(i)
    env, ns, na = getenv()
    learner = ActorCriticPolicyGradient(na = na, ns = ns, α = .02,
                                        αcritic = 0.0, nsteps = 25)
    setup(learner, env)
end
function dqn(i)
    env, ns, na = getenv()
    learner = DQN(Chain(Dense(ns, 48, relu), Dense(48, 24, relu), Dense(24, na)),
                  updateevery = 1, updatetargetevery = 100,
                  startlearningat = 50, minibatchsize = 32,
                  doubledqn = false, replaysize = 10^3,
                  opttype = x -> ADAM(x, .0005))
    setup(learner, env)
end
function tilingsarsa(i)
    env, ns, na = getenv()
    if typeof(env) <: CartPole
        high = [2.6, 4, .24, 3.4]
        low = -high
        nbins = 20 * ones(Int64, 4)
        initvalue = 200
    elseif typeof(env) <: MountainCar
        high = [.5, .07]
        low = [-1.2, -.07]
        nbins = [16, 16]
        initvalue = 0
    end
    p0 = StateAggregator(low, high, nbins)
    preprocessor = TilingStateAggregator(p0, 8)
    learner = Sarsa(ns = 8*prod(nbins), na = na, λ = 0, α = .2/8,
                    initvalue = initvalue)
    s = setup(learner, env, preprocessor)
    s.policy.ϵ = 0
    s
end
rlsetupcreators = Dict("linear ACPG" => acpg, "DQN" => dqn,
                       "tiling Sarsa" => tilingsarsa)
end # everywhere
@time res = compare(rlsetupcreators, 10, verbose = true)

@save tempname() * ".jld2" res

a = plotcomparison(res);
a["xlabel"] = "epochs";
a["ylabel"] = "average episode length greedy policy"
a


text:
recurs_sym_type(ex::Any) =
  (ex==None || typeof(ex)==Symbol || length(ex.args)==1) ? eval(ex) : Expr(ex.head, ex.args[1], recurs_sym_type(ex.args[2]))
macro c(ret_type, func, arg_types, lib)
  local _arg_types = Expr(:tuple, [recurs_sym_type(a) for a in arg_types.args]...)
  local _ret_type = recurs_sym_type(ret_type)
  local _args_in = Any[ symbol(string('a',x)) for x in 1:length(_arg_types.args) ]
  local _lib = eval(lib)
  quote
    $(esc(func))($(_args_in...)) = ccall( ($(string(func)), $(Expr(:quote, _lib)) ), $_ret_type, $_arg_types, $(_args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  real_t = recurs_sym_type(real_t)
  quote
    typealias $fake_t $real_t
  end
end

# header: /usr/local/include/ida/ida_band.h
@ctypedef IDADlsDenseJacFn Ptr{:Void}
@ctypedef IDADlsBandJacFn Ptr{:Void}
@c Int32 IDADlsSetDenseJacFn (Ptr{:None},:IDADlsDenseJacFn) shlib
@c Int32 IDADlsSetBandJacFn (Ptr{:None},:IDADlsBandJacFn) shlib
@c Int32 IDADlsGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 IDADlsGetNumJacEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDADlsGetNumResEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDADlsGetLastFlag (Ptr{:None},Ptr{:Int64}) shlib
@c Ptr{:Uint8} IDADlsGetReturnFlagName (:Int64,) shlib
@c Int32 IDABand (Ptr{:None},:Int64,:Int64,:Int64) shlib

# header: /usr/local/include/ida/ida_bbdpre.h
@ctypedef IDABBDLocalFn Ptr{:Void}
@ctypedef IDABBDCommFn Ptr{:Void}
@c Int32 IDABBDPrecInit (Ptr{:None},:Int64,:Int64,:Int64,:Int64,:Int64,:realtype,:IDABBDLocalFn,:IDABBDCommFn) shlib
@c Int32 IDABBDPrecReInit (Ptr{:None},:Int64,:Int64,:realtype) shlib
@c Int32 IDABBDPrecGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 IDABBDPrecGetNumGfnEvals (Ptr{:None},Ptr{:Int64}) shlib

# header: /usr/local/include/ida/ida_dense.h
@c Int32 IDADense (Ptr{:None},:Int64) shlib

# header: /usr/local/include/ida/ida_direct.h

# header: /usr/local/include/ida/ida.h
@ctypedef IDAResFn Ptr{:Void}
@ctypedef IDARootFn Ptr{:Void}
@ctypedef IDAEwtFn Ptr{:Void}
@ctypedef IDAErrHandlerFn Ptr{:Void}
@c Ptr{:None} IDACreate () shlib
@c Int32 IDASetErrHandlerFn (Ptr{:None},:IDAErrHandlerFn,Ptr{:None}) shlib
@c Int32 IDASetErrFile (Ptr{:None},Ptr{:FILE}) shlib
@c Int32 IDASetUserData (Ptr{:None},Ptr{:None}) shlib
@c Int32 IDASetMaxOrd (Ptr{:None},:Int32) shlib
@c Int32 IDASetMaxNumSteps (Ptr{:None},:Int64) shlib
@c Int32 IDASetInitStep (Ptr{:None},:realtype) shlib
@c Int32 IDASetMaxStep (Ptr{:None},:realtype) shlib
@c Int32 IDASetStopTime (Ptr{:None},:realtype) shlib
@c Int32 IDASetNonlinConvCoef (Ptr{:None},:realtype) shlib
@c Int32 IDASetMaxErrTestFails (Ptr{:None},:Int32) shlib
@c Int32 IDASetMaxNonlinIters (Ptr{:None},:Int32) shlib
@c Int32 IDASetMaxConvFails (Ptr{:None},:Int32) shlib
@c Int32 IDASetSuppressAlg (Ptr{:None},:Int32) shlib
@c Int32 IDASetId (Ptr{:None},:N_Vector) shlib
@c Int32 IDASetConstraints (Ptr{:None},:N_Vector) shlib
@c Int32 IDASetRootDirection (Ptr{:None},Ptr{:Int32}) shlib
@c Int32 IDASetNoInactiveRootWarn (Ptr{:None},) shlib
@c Int32 IDAInit (Ptr{:None},:IDAResFn,:realtype,:N_Vector,:N_Vector) shlib
@c Int32 IDAReInit (Ptr{:None},:realtype,:N_Vector,:N_Vector) shlib
@c Int32 IDASStolerances (Ptr{:None},:realtype,:realtype) shlib
@c Int32 IDASVtolerances (Ptr{:None},:realtype,:N_Vector) shlib
@c Int32 IDAWFtolerances (Ptr{:None},:IDAEwtFn) shlib
@c Int32 IDASetNonlinConvCoefIC (Ptr{:None},:realtype) shlib
@c Int32 IDASetMaxNumStepsIC (Ptr{:None},:Int32) shlib
@c Int32 IDASetMaxNumJacsIC (Ptr{:None},:Int32) shlib
@c Int32 IDASetMaxNumItersIC (Ptr{:None},:Int32) shlib
@c Int32 IDASetLineSearchOffIC (Ptr{:None},:Int32) shlib
@c Int32 IDASetStepToleranceIC (Ptr{:None},:realtype) shlib
@c Int32 IDARootInit (Ptr{:None},:Int32,:IDARootFn) shlib
@c Int32 IDACalcIC (Ptr{:None},:Int32,:realtype) shlib
@c Int32 IDASolve (Ptr{:None},:realtype,Ptr{:realtype},:N_Vector,:N_Vector,:Int32) shlib
@c Int32 IDAGetDky (Ptr{:None},:realtype,:Int32,:N_Vector) shlib
@c Int32 IDAGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 IDAGetNumSteps (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDAGetNumResEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDAGetNumLinSolvSetups (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDAGetNumErrTestFails (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDAGetNumBacktrackOps (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDAGetConsistentIC (Ptr{:None},:N_Vector,:N_Vector) shlib
@c Int32 IDAGetLastOrder (Ptr{:None},Ptr{:Int32}) shlib
@c Int32 IDAGetCurrentOrder (Ptr{:None},Ptr{:Int32}) shlib
@c Int32 IDAGetActualInitStep (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 IDAGetLastStep (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 IDAGetCurrentStep (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 IDAGetCurrentTime (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 IDAGetTolScaleFactor (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 IDAGetErrWeights (Ptr{:None},:N_Vector) shlib
@c Int32 IDAGetEstLocalErrors (Ptr{:None},:N_Vector) shlib
@c Int32 IDAGetNumGEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDAGetRootInfo (Ptr{:None},Ptr{:Int32}) shlib
@c Int32 IDAGetIntegratorStats (Ptr{:None},Ptr{:Int64},Ptr{:Int64},Ptr{:Int64},Ptr{:Int64},Ptr{:Int32},Ptr{:Int32},Ptr{:realtype},Ptr{:realtype},Ptr{:realtype},Ptr{:realtype}) shlib
@c Int32 IDAGetNumNonlinSolvIters (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDAGetNumNonlinSolvConvFails (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDAGetNonlinSolvStats (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Ptr{:Uint8} IDAGetReturnFlagName (:Int64,) shlib
@c None IDAFree (Ptr{Ptr{:None}},) shlib

# header: /usr/local/include/ida/ida_impl.h
@ctypedef IDAMem Ptr{:Void}
@c Int32 IDAEwtSet (:N_Vector,:N_Vector,Ptr{:None}) shlib
@c None IDAProcessError (:IDAMem,:Int32,Ptr{:Uint8},Ptr{:Uint8},Ptr{:Uint8}) shlib
@c None IDAErrHandler (:Int32,Ptr{:Uint8},Ptr{:Uint8},Ptr{:Uint8},Ptr{:None}) shlib

# header: /usr/local/include/ida/ida_spbcgs.h
@ctypedef IDASpilsPrecSetupFn Ptr{:Void}
@ctypedef IDASpilsPrecSolveFn Ptr{:Void}
@ctypedef IDASpilsJacTimesVecFn Ptr{:Void}
@c Int32 IDASpilsSetPreconditioner (Ptr{:None},:IDASpilsPrecSetupFn,:IDASpilsPrecSolveFn) shlib
@c Int32 IDASpilsSetJacTimesVecFn (Ptr{:None},:IDASpilsJacTimesVecFn) shlib
@c Int32 IDASpilsSetGSType (Ptr{:None},:Int32) shlib
@c Int32 IDASpilsSetMaxRestarts (Ptr{:None},:Int32) shlib
@c Int32 IDASpilsSetMaxl (Ptr{:None},:Int32) shlib
@c Int32 IDASpilsSetEpsLin (Ptr{:None},:realtype) shlib
@c Int32 IDASpilsSetIncrementFactor (Ptr{:None},:realtype) shlib
@c Int32 IDASpilsGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 IDASpilsGetNumPrecEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDASpilsGetNumPrecSolves (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDASpilsGetNumLinIters (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDASpilsGetNumConvFails (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDASpilsGetNumJtimesEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDASpilsGetNumResEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 IDASpilsGetLastFlag (Ptr{:None},Ptr{:Int64}) shlib
@c Ptr{:Uint8} IDASpilsGetReturnFlagName (:Int64,) shlib
@c Int32 IDASpbcg (Ptr{:None},:Int32) shlib

# header: /usr/local/include/ida/ida_spgmr.h
@c Int32 IDASpgmr (Ptr{:None},:Int32) shlib

# header: /usr/local/include/ida/ida_spils.h

# header: /usr/local/include/ida/ida_sptfqmr.h
@c Int32 IDASptfqmr (Ptr{:None},:Int32) shlib


text:
# module DFSPlanner
#
# using Parameters
# using LightGraphs
# using DataStructures
# using JLD2, FileIO
# using CRCBS
#
# using ..TaskGraphs

export
    DFSRoutePlanner

"""
    DFSRoutePlanner

Prioritized Depth-First Search route planner.
"""
@with_kw struct DFSRoutePlanner{C}
    logger::SolverLogger{C} = SolverLogger{C}()
end
# function construct_heuristic_model(solver::DFSRoutePlanner,env_graph;
#         ph = PerfectHeuristic(get_dist_matrix(env_graph)),
#         kwargs...)
#     construct_composite_heuristic(ph,ph,NullHeuristic())
# end
# function construct_cost_model(solver::DFSRoutePlanner,
#         schedule, cache, problem_spec, env_graph;
#         extra_T=400, primary_objective=SumOfMakeSpans(), kwargs...)
#     cost_model = construct_composite_cost_model(
#         typeof(primary_objective)(schedule,cache),
#         FullCostModel(maximum,TravelTime()),
#         FullCostModel(maximum,TravelDistance())
#         )
#     heuristic_model = construct_heuristic_model(solver,env_graph;kwargs...)
#     # ph = PerfectHeuristic(get_dist_matrix(env_graph))
#     # heuristic_model = construct_composite_heuristic(ph,ph,NullHeuristic())
#     cost_model, heuristic_model
# end

export
    # sorted_actions,
    get_conflict_idx,
    select_action_dfs!,
    get_next_node_matching_agent_id,
    update_envs!,
    prioritized_dfs_search

@with_kw struct DFS_SearchState
    pickup_i::Int   = 0 # if i < start, iterate from current action (inclusive). Defines where to "pick back up"
    reset_i::Int    = 0 # if i > reset_i, iterate through all possible actions. Defines the point below which to "reset" the action vector
end
function update_search_state(s,i)
    s = DFS_SearchState(s,reset_i = max(s.reset_i,s.pickup_i+1))
    if s.pickup_i > i # reset pickup_i once it has been exceeded
        s = DFS_SearchState(s, pickup_i = 0)
    end
    if s.reset_i < i
       s = DFS_SearchState(s,reset_i = i)
    end
    s
end
# function sorted_actions(env,s)
#     f = (s,a,sp)->add_heuristic_cost(env,get_transition_cost(env,s,a,sp),get_heuristic_cost(env,sp))
#     sort(
#         collect(get_possible_actions(env,s)),
#         by=a->f(s,a,get_next_state(env,s,a))
#     )
# end

"""
    get_conflict_idx(envs,states,actions,i,ordering,idxs)

Check if the planned action for agent with priority `i` conflicts with the
    planned action of any agent with higher priority.
"""
function get_conflict_idx(envs,states,actions,i,ordering,idxs)
    idx = ordering[i]
    env = envs[idx]
    s = states[idx]
    a = actions[idx]
    pi = PathNode(s,a,get_next_state(env,s,a))
    for (j,idx) in enumerate(ordering[1:max(1,i-1)])
        if j == i
            continue
        end
        env = envs[idx]
        s = states[idx]
        a = actions[idx]
        pj = PathNode(s,a,get_next_state(env,s,a))
        if detect_state_conflict(pi,pj) || detect_action_conflict(pi,pj)
            return j
        end
    end
    return -1
end
function update_envs!(solver,search_env,envs,paths)
    cache = get_cache(search_env)
    schedule = get_schedule(search_env)
    cbs_node = initialize_root_node(search_env)
    # update_planning_cache!(solver,search_env)        cache.tF[v] = get_final_state(path).t

    # update cache times
    up_to_date = true
    for i in 1:length(envs)
        env = envs[i]
        path = paths[i]
        v = get_vtx(schedule,env.node_id)
        s = get_final_state(path)
        t_arrival = max(cache.tF[v], s.t + get_distance(get_graph(search_env).dist_function,s.vtx,env.goal.vtx))
        if is_goal(envs[i],s)
            if t_arrival > cache.tF[v] && env.goal.vtx != -1
                @log_info(-1,verbosity(solver),"DFS update_envs!(): extending tF[v] from $(cache.tF[v]) to $t_arrival in ",string(get_node(env)),", s = ",string(s))
                cache.tF[v] = t_arrival
                up_to_date = false
            end
        end
    end
    # reset cache
    if !up_to_date
        t0,tF,slack,local_slack = process_schedule(schedule;t0=cache.t0,tF=cache.tF)
        cache.t0            .= t0
        cache.tF            .= tF
        cache.slack         .= slack
        cache.local_slack   .= local_slack
        # rebuild envs to reset goal time
        for i in 1:length(envs)
            env = envs[i]
            path = paths[i]
            v = get_vtx(schedule,env.node_id)
            envs[i] = build_env(solver,search_env,cbs_node,v)
        end
    end
    # mark finished envs as complete
    for i in 1:length(envs)
        env = envs[i]
        path = paths[i]
        if is_goal(envs[i],get_final_state(path))
            # update schedule and cache
            v = get_vtx(schedule,env.node_id)
            if !(v in cache.closed_set)
                update_planning_cache!(solver,search_env,v,path)
                @assert v in cache.closed_set
                @assert i == env.agent_idx
            end
        end
    end
    update_planning_cache!(solver,search_env)
    i = 0
    while i < length(envs)
        i += 1
        env = envs[i]
        path = paths[i]
        v = get_vtx(schedule,env.node_id)
        if is_goal(envs[i],get_final_state(path))
            # update schedule and cache
            if !(v in cache.closed_set)
                update_planning_cache!(solver,search_env,v,path)
                update_planning_cache!(solver,search_env)
            end
            @assert v in cache.closed_set
            @assert i == env.agent_idx
            for v2 in outneighbors(schedule,v)
                if get_path_spec(schedule,v).agent_id == i
                    if !(v2 in cache.active_set)
                        @log_info(-1,solver.l2_verbosity,"node ",string(get_node_from_vtx(schedule,v2))," not yet in active set")
                        @log_info(-1,solver.l2_verbosity,"inneighbors of ",string(get_node_from_vtx(schedule,v2))," are ",map(v3->string(get_node_from_vtx(schedule,v3)),neighborhood(schedule,v2,3,dir=:in))...)
                    end
                end
            end
            # swap out env for new env
            node_id = get_next_node_matching_agent_id(schedule,cache,env.agent_idx)
            @assert node_id != env.node_id
            if get_vtx(schedule,node_id) in cache.active_set
                envs[i] = build_env(solver,search_env,cbs_node,get_vtx(schedule,node_id))
                # i = 0
                i -= 1
            else
                node_string = string(get_node_from_id(get_schedule(search_env),node_id))
                @log_info(4,verbosity(solver),"cannot update environment for agent $i because next node ",node_string," not in cache.active_set")
            end
        end
    end
    update_cost_model!(search_env)
    envs,paths
end
function select_ordering(solver,search_env,envs)
    schedule = get_schedule(search_env)
    cache = get_cache(search_env)
    ordering = sort(
        collect(1:search_env.num_agents),
        by = i->(
            ~isa(envs[i].schedule_node,Union{COLLECT,DEPOSIT}),
            ~isa(envs[i].schedule_node,CARRY),
            minimum(cache.slack[get_vtx(schedule,envs[i].node_id)])
            )
        )
    ordering
end
function select_action_dfs!(solver,envs,states,actions,i,ordering,idxs,search_state=SearchState())
    search_state = update_search_state(search_state,i)
    if i <= 0
        return false
    elseif i > length(states)
        for (env,s,a) in zip(envs,states,actions)
            if a != CRCBS.wait(env,s) || length(get_possible_actions(env,s)) == 1 || s.t < env.goal.t
                return true
            end
        end
        @log_info(4,verbosity(solver),"action vector $(map(a->string(a),actions)) not eligible")
        return false
    elseif !(ordering[i] in idxs)
        return select_action_dfs!(solver,envs,states,actions,i+1,ordering,idxs,search_state)
    else
        # idx = env.ordering_map[i]zs
        j = 0
        idx = ordering[i]
        env = envs[idx]
        s = states[idx]
        for ai in sorted_actions(env,s)
            a = actions[idx]
            c = get_transition_cost(env,s,ai)
            c0 = get_transition_cost(env,s,a)
            if (i >= search_state.reset_i) || (i < search_state.pickup_i && a == ai) || ((c >= c0 || is_valid(env,a)) && a != ai)
                actions[idx] = ai
                @log_info(5,verbosity(solver),"$(repeat(" ",i))i = $i, trying a=",string(ai)," from s = ",string(s),"for env ",string(get_node(env)), " with env.goal = ",string(env.goal))
                k = get_conflict_idx(envs,states,actions,i,ordering,idxs)
                # @assert k < i "should only check for conflicts with 1:$i, but found conflict with $k"
                if k <= 0
                    if select_action_dfs!(solver,envs,states,actions,i+1,ordering,idxs,search_state)
                        return true
                    end
                elseif !(ordering[k] in idxs)
                    # i = i - 1
                    # break
                    # return false
                else
                    @log_info(5,verbosity(solver),"--- conflict betweeh $i and $k")
                    j = max(k,j)
                end
            end
        end
        if j <= 0
            return false
        end
        # if j > 0
        search_state = DFS_SearchState(pickup_i=j,reset_i=0)
        # search_states[j] = search_states[i]
        # i = j
        # break
        # end
        return select_action_dfs!(solver,envs,states,actions,j,ordering,idxs,search_state) #
    end
    # end
end
function prioritized_dfs_search(solver,search_env,envs,paths;
        t0 = max(0, minimum(map(path->length(path),paths))),
        search_state = DFS_SearchState()
        )
     tip_times = map(path->length(path),paths)
     t = t0
     states     = map(path->get_s(get_path_node(path,t+1)), paths)
     actions    = map(path->get_a(get_path_node(path,t+1)), paths)
     while true
         increment_iteration_count(solver)
         enforce_iteration_limit!(solver)
         update_envs!(solver,search_env,envs,paths)
         if length(get_cache(search_env).active_set) == 0
             return envs, paths, true
         end
         @assert all(map(p->length(p) >= t,paths))
         ordering   = select_ordering(solver,search_env,envs)
         idxs       = Set(findall(tip_times .<= t))
         if select_action_dfs!(solver,envs,states,actions,1,ordering,idxs,search_state)
             # step forward in search
             for idx in idxs
                 env = envs[idx]
                 path = paths[idx]
                 s = states[idx]
                 a = actions[idx]
                 sp = get_next_state(env,s,a)
                 push!(path,PathNode(s,a,sp))
                 set_cost!(path, accumulate_cost(env,get_cost(path),get_transition_cost(env,s,a,sp)))
             end
             t += 1
             search_state = DFS_SearchState()
             @log_info(3,verbosity(solver),"stepping forward, t = $t")
         else
             # step backward in search
             all(tip_times .> t) == 0 ? break : nothing
             for (idx,path) in enumerate(paths)
                 if tip_times[idx] < t
                     pop!(path.path_nodes)
                     @assert length(path) == t-1
                 end
             end
             t -= 1
             idxs    = Set(findall(tip_times .<= t))
              # start new search where previous search left off
             search_state = DFS_SearchState(pickup_i = maximum(idxs))
             @log_info(-1,verbosity(solver),"stepping backward, t = $t")
         end
         states     = map(path->get_s(get_path_node(path,t+1)), paths)
         actions    = map(path->get_a(get_path_node(path,t+1)), paths)
     end
     return envs, paths, false
end

"""
    Iterate over agents in order of priority, allowing them to fill in a
    reservation table for the vtxs they would like to occupy at the next time
    step(s).
"""
function CRCBS.solve!(
    solver::DFSRoutePlanner,
    mapf::P;kwargs...) where {P<:PC_MAPF}

    search_env = mapf.env
    update_planning_cache!(solver,search_env) # NOTE to get rid of all nodes that don't need planning but are in the active set

    route_plan = deepcopy(get_route_plan(search_env))
    paths = get_paths(route_plan)
    # initialize envs
    envs = Vector{PCCBS.LowLevelEnv}([PCCBS.LowLevelEnv() for p in paths])
    cbs_node = initialize_root_node(search_env)
    for i in 1:search_env.num_agents
        node_id = get_next_node_matching_agent_id(get_schedule(search_env),get_cache(search_env),i)
        envs[i] = build_env(solver,search_env,cbs_node,get_vtx(get_schedule(search_env),node_id))
    end

    envs, paths, status = prioritized_dfs_search(solver,search_env,envs,paths;
        max_iters = solver.cbs_model.max_iters
    )
    if validate(get_schedule(search_env),convert_to_vertex_lists(route_plan),get_cache(search_env).t0,get_cache(search_env).tF)
        @log_info(0,verbosity(solver),"DFS: Succeeded in finding a valid route plan!")
    else
        # throw(SolverCBSMaxOutException("ERROR in DFS! Failed to find a valid route plan!"))
        # DUMP
        filename = joinpath(DEBUG_PATH,string("DFS_demo",get_debug_file_id(),".jld2"))
        mkpath(DEBUG_PATH)
        for env in envs
            v = get_vtx(get_schedule(search_env),env.node_id)
            @log_info(-1,verbosity(solver),"node ",string(get_node(env))," t0=$(get_t0(search_env,v)), tF=$(get_tF(search_env,v)), closed=$(v in get_cache(search_env).closed_set),")
        end
        robot_paths = convert_to_vertex_lists(route_plan)
        object_paths, object_intervals, object_ids, path_idxs = get_object_paths(route_plan,search_env)
        @log_info(-1,verbosity(solver),"Dumping DFS route plan to $filename")
        @save filename robot_paths object_paths object_intervals object_ids path_idxs

        throw(AssertionError("ERROR in DFS! Failed to find a valid route plan!"))
    end
    # for (path,base_path,env) in zip(paths,get_route_plan(search_env).paths,envs)
    #     c = base_path.cost
    #     for p in path.path_nodes[length(base_path)+1:end]
    #         c = accumulate_cost(env,c,get_transition_cost(env,p.s,p.a,p.sp))
    #     end
    #     path.cost = c
    # end
    cost = aggregate_costs(get_cost_model(search_env),map(p->get_cost(p),paths))

    return route_plan, get_cache(search_env), cost
end

# end


text:
using Images, File
import DataDeps: unpack_cmd, unpack

#TODO: not returning extracted path
#TODO: performance is slow and sometimes progress is not visible
#TODO: tests pending
function unpack(f; filepath=pwd(), keep_originals=false)
    run(unpack_cmd(f, filepath, last(splitext(f)), last(splitext(first(splitext(f))))))
    rm("pax_global_header"; force=true)
    !keep_originals && rm(f)
end

#=load_image
read an image file and load it to memory

It is a wrapper over the Images.load functions=#
function load_image(path::AbstractString)
	if isfile(path)
		return load(path)
	end
end

#=
image2array
convert image to it's array representation

it is a wrapper over the Images.channelview function=#
function image2array(image)
	convert(Array{Float32}, channelview(image))
end


text:
using ExprOptimization
using Test, Random

let
    grammar = @grammar begin
        R = |(1:3)
        R = R + R
    end

    function loss(node::RuleNode, grammar::Grammar)
        Core.eval(node, grammar)
    end

    Random.seed!(0)
    p = MonteCarlo(20, 5)
    res = optimize(p, grammar, :R, loss)
    @test res.expr == 1
    @test Core.eval(res.tree, grammar) == 1
    @test res.loss == 1

    p = MonteCarlo(20, 5; track_method=MonteCarlos.TopKTracking(3))
    res = optimize(p, grammar, :R, loss)
    res.alg_result[:top_k]
end


text:
using CodecZlib
using Statistics
export compress, decompress

struct CompressedMask{A}
  size::Tuple{Int,Int}
  axes::A
  data::Array{UInt8}
end

# NOTE: this somewhat odd duplication of types is to deal with legacy code. I
# want the newly stored data to use only builtin types of julia, which save to
# jld2 files in a simple format. This helps ensure it can be readily loaded in
# with future versions of julia and in other languages

struct CMask{T,F}
  size::Tuple{Int,Int}
  data::Array{UInt8}
  times_s::T
  freqs_Hz::F
end

function CMask(x::CompressedMask)
  @assert x.axes[1] isa Axis{:time}
  @assert x.axes[2] isa Axis{:freq}
  @assert unit(eltype(x.axes[1])) == s
  @assert unit(eltype(x.axes[2])) == Hz

  CMask(x.size,x.data,ustrip.(x.axes[1].val),
                         ustrip.(x.axes[2].val))
end

function CompressedMask(x::CMask)
  CompressedMask(x.size, (Axis{:time}(s .* x.times_s),
                          Axis{:freq}(Hz .* x.freqs_Hz)), x.data)
end

compress_axes(x::MetaUnion{AxisArray}) = AxisArrays.axes(x)
compress_axes(x) = nothing

function compress(x::AbstractMatrix)
  quantized = Array{UInt8}(undef,size(x))
  @. quantized = floor(UInt8,max(0,x/$maximum(x))*typemax(UInt8))
  CMask(CompressedMask(size(x), compress_axes(x),
                                        transcode(ZlibCompressor,vec(quantized))))
end

withaxes(x,axes) = AxisArray(x, axes...)
withaxes(x,::Nothing) = x

decompress(x::CMask) = decompress(CompressedMask(x))
function decompress(x::CompressedMask)
  mask = transcode(ZlibDecompressor,x.data) ./ typemax(UInt8)
  withaxes(reshape(mask,x.size...), x.axes)
end

function CorticalSpecrtralTemporalResponses.audiospect(x::CompressedMask,settings)
  settings = read_settings(settings)
  audiospect(decompress(x);settings.freqs.analyze...)
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
using Symbolics

nx = 6
nu = 1
@variables δx[1:nx] δu[1:nu] Bf[1:nx*nu]

function cost(δx, δu, Bf)
    B = reshape(Bf, nx, nu)
    return 1/2 * (δx - B*δu)' * (δx - B*δu)
end

c = cost(δx, δu, Bf)
dc = Symbolics.gradient(c, Bf)
ddc = Symbolics.hessian(c, Bf)

c_fct = eval(Symbolics.build_function([c], δx, δu, Bf)[1])
dc_fct = eval(Symbolics.build_function(dc, δx, δu, Bf)[1])
ddc_fct = eval(Symbolics.build_function(ddc, δx, δu, Bf)[1])

δx0 = rand(nx)
δu0 = rand(nu)
Bf0 = rand(nx*nu)

c_fct(δx0, δu0, Bf0)
dc_fct(δx0, δu0, Bf0)
ddc_fct(δx0, δu0, Bf0)

function leastsquares(δx, δu)
    N = length(δx)
    nx = length(δx[1])
    nu = length(δu[1])
    Bf = zeros(nx*nu)
    e = 0.0
    grad = zeros(nx*nu)
    hess = zeros(nx*nu, nx*nu)
    for i = 1:N
        e += c_fct(δx[i], δu[i], Bf)[1]
        grad += dc_fct(δx[i], δu[i], Bf)
        hess += ddc_fct(δx[i], δu[i], Bf)
    end
    Bf = Bf - hess \ grad
    B = reshape(Bf, (nx, nu))
    return B
end


text:
using MomentClosure
using MomentClosure: define_μ, define_M
using Test

@testset "creating symbolic moments" begin

    μ = define_μ(3, 3)
    @test μ[(0,0,0)] == 1
    @test string(μ[(2,1,0)]) == "μ₂₁₀(t)"
    @test string(μ[(0,0,3)]) == "μ₀₀₃(t)"
    M = define_M(3, 3)
    @test M[(1,0,0)] == 0
    @test string(M[1,0,2]) == "M₁₀₂(t)"

end


text:
using XSV

# Run on a small sample first
xsv_string, xsv_data = genxsv(num_rows=5, num_cols=2)
xsv_string_h = join(["C"*string(i) for i=1:length(xsv_data[1])], ",") * "\n" * xsv_string
xsv = readxsv(xsv_string)
for row in iterxsv(xsv_string) end
for row in iterxsvh(xsv_string_h) end
for row in iterxsvht(xsv_string_h, row_type_name="XsvTestRow") end

xsv_string, xsv_data = genxsv(num_rows=20_000, num_cols=20)
println("readxsv:")
@time xsv = readxsv(xsv_string)
println("\niterxsv:")
@time for row in iterxsv(xsv_string) end

xsv_string_h = join(["C"*string(i) for i=1:length(xsv_data[1])], ",") * "\n" * xsv_string
println("\niterxsvh:")
@time for row in iterxsvh(xsv_string_h) end

println("\niterxsvht:")
@time for row in iterxsvht(xsv_string_h) end


text:
__precompile__()

module DrakeVisualizer

using PyLCM
using GeometryTypes
import GeometryTypes: origin, radius
import Meshing
import PyCall: pyimport, PyObject, PyNULL, PyVector
import Rotations: Rotation, Quat
import CoordinateTransformations: Transformation, transform_deriv, IdentityTransformation, AbstractAffineMap
import ColorTypes: RGBA, Colorant, red, green, blue, alpha
import StaticArrays: SVector
import Base: convert, length

export GeometryData,
        Link,
        Robot,
        Visualizer,
        HyperRectangle,
        HyperEllipsoid,
        HyperCylinder,
        HyperSphere,
        HyperCube,
        contour_mesh,
        draw,
        reload_model

const drakevis = PyNULL()

destructure{T,N}(A::Array{T,N}) = reinterpret(eltype(T), A, (size(T)..., size(A)...))

# GeometryTypes doesn't define an Ellipsoid type yet, so we'll make one ourselves!
type HyperEllipsoid{N, T} <: GeometryPrimitive{N, T}
    center::Point{N, T}
    radii::Vec{N, T}
end

origin{N, T}(geometry::HyperEllipsoid{N, T}) = geometry.center
radii{N, T}(geometry::HyperEllipsoid{N, T}) = geometry.radii

type HyperCylinder{N, T} <: GeometryTypes.GeometryPrimitive{N, T}
    length::T # along last axis
    radius::T
    # origin is at center
end

length(geometry::HyperCylinder) = geometry.length
radius(geometry::HyperCylinder) = geometry.radius

type GeometryData{TransformType <: Transformation, GeometryType <: AbstractGeometry}
    geometry::GeometryType
    transform::TransformType
    color::RGBA{Float64}
end

GeometryData{TransformType <: Transformation, GeometryType <: AbstractGeometry}(geometry::GeometryType,
    transform::TransformType=IdentityTransformation(),
    color=RGBA{Float64}(1., 0, 0, 0.5)) = GeometryData{TransformType, GeometryType}(geometry, transform, convert(RGBA{Float64}, color))

function contour_mesh(f::Function,
                      bounds::HyperRectangle,
                      isosurface_value::Number=0.0,
                      resolution::Real=maximum(widths(bounds))/20.0)
    # Extract an isosurface from a vector-input function in 3 dimensions and
    # return it as a mesh. This function mostly exists to work around bugs in
    # Meshing.jl: specifically its weird handling of non-zero isosurfaces
    # and its incorrect mesh point scaling:
    # https://github.com/JuliaGeometry/GeometryTypes.jl/issues/49
    window_width = widths(bounds)
    sdf = SignedDistanceField(x -> f(SVector{3, Float64}(x[1], x[2], x[3])) - isosurface_value,
                              bounds,
                              resolution)

    # We've already accounted for the isosurface value in the construction of
    # the SDF, so we set the iso value here to 0.
    mesh = HomogenousMesh(sdf, 0.0)

    # Rescale the mesh points and then construct a new mesh using the rescaled
    # points and the original faces.
    lb = minimum(bounds)
    rescaled_points = Point{3,Float64}[Vec(v-1) ./ (Vec(size(sdf))-1) .* window_width .+ lb for v in vertices(mesh)]
    HomogenousMesh(rescaled_points, mesh.faces)
end

contour_mesh(f::Function,
             lower_bound::Vec,
             upper_bound::Vec,
             isosurface_value::Number=0.0,
             resolution::Real=maximum(upper_bound - lower_bound)/20.0) = (
    contour_mesh(f,
                 HyperRectangle(lower_bound, upper_bound - lower_bound),
                 isosurface_value,
                 resolution))

contour_mesh(f::Function,
             lower_bound::AbstractVector,
             upper_bound::AbstractVector,
             isosurface_value::Number=0.0,
             resolution::Real=maximum(upper_bound - lower_bound)/20.0) = (
    contour_mesh(f,
                 Vec{3, Float64}(lower_bound[1], lower_bound[2], lower_bound[3]),
                 Vec{3, Float64}(upper_bound[1], upper_bound[2], upper_bound[3]),
                 isosurface_value,
                 resolution))

type Link
    geometry_data::Vector{GeometryData}
    name::String
end
Link{T <: GeometryData}(geometry_data::Vector{T}) = Link(geometry_data, "link")

type Robot
    links::Vector{Link}
end

convert(::Type{Link}, geom::GeometryData) = Link([geom])
convert(::Type{Link}, geometry::AbstractGeometry) = Link(GeometryData(geometry))
convert(::Type{Robot}, link::Link) = Robot([link])
convert(::Type{Robot}, links::AbstractArray{Link}) = Robot(convert(Vector{Link}, links))
convert(::Type{Robot}, geom::GeometryData) = convert(Robot, convert(Link, geom))
convert(::Type{Robot}, geometry::AbstractGeometry) = convert(Robot, GeometryData(geometry))

to_lcm(q::Quat) = SVector{4, Float64}(q.w, q.x, q.y, q.z)
to_lcm_quaternion(matrix::AbstractMatrix) = to_lcm(Quat(matrix))
to_lcm_quaternion(matrix::UniformScaling) = to_lcm(Quat(1.0, 0, 0, 0))
to_lcm_quaternion(transform::AbstractAffineMap) = to_lcm_quaternion(transform_deriv(transform, SVector{3, Float64}(0,0,0)))
to_lcm_quaternion(transform::IdentityTransformation) = to_lcm_quaternion(Quat(1.0, 0, 0, 0))

to_lcm(color::Colorant) = Float64[red(color); green(color); blue(color); alpha(color)]

center(geometry::HyperRectangle) = minimum(geometry) + 0.5 * widths(geometry)
center(geometry::HyperCube) = minimum(geometry) + 0.5 * widths(geometry)
center(geometry::HyperSphere) = origin(geometry)
center(geometry::HyperEllipsoid) = origin(geometry)

function fill_geometry_data!(msg::PyObject, geometry::AbstractMesh, transform::Transformation)
    msg[:type] = msg[:MESH]
    msg[:position] = transform(SVector{3, Float64}(0, 0, 0))
    msg[:quaternion] = to_lcm_quaternion(transform)
    msg[:string_data] = ""
    float_data = Float64[length(vertices(geometry));
        length(faces(geometry));
        vec(destructure(vertices(geometry)));
        vec(destructure(map(f -> convert(Face{3, Int, -1}, f), faces(geometry))))]
    msg[:float_data] = float_data
    msg[:num_float_data] = length(float_data)
end

function fill_geometry_data!(msg::PyObject, geometry::HyperRectangle, transform::Transformation)
    msg[:type] = msg[:BOX]
    msg[:position] = SVector{3, Float64}(transform(center(geometry))...)
    msg[:quaternion] = to_lcm_quaternion(transform)
    msg[:string_data] = ""
    float_data = convert(Vector, widths(geometry))
    msg[:float_data] = float_data
    msg[:num_float_data] = length(float_data)

end

function fill_geometry_data!(msg::PyObject, geometry::HyperCube, transform::Transformation)
    msg[:type] = msg[:BOX]
    msg[:position] = SVector{3, Float64}(transform(center(geometry))...)
    msg[:quaternion] = to_lcm_quaternion(transform)
    msg[:string_data] = ""
    float_data = convert(Vector, widths(geometry))
    msg[:float_data] = float_data
    msg[:num_float_data] = length(float_data)
end

function fill_geometry_data!(msg::PyObject, geometry::HyperSphere, transform::Transformation)
    msg[:type] = msg[:SPHERE]
    msg[:position] = SVector{3, Float64}(transform(center(geometry))...)
    msg[:quaternion] = to_lcm_quaternion(transform)
    msg[:string_data] = ""
    msg[:float_data] = [radius(geometry)]
    msg[:num_float_data] = 1
end

function fill_geometry_data!(msg::PyObject, geometry::HyperEllipsoid, transform::Transformation)
    msg[:type] = msg[:ELLIPSOID]
    msg[:position] = SVector{3, Float64}(transform(center(geometry))...)
    msg[:quaternion] = to_lcm_quaternion(transform)
    msg[:string_data] = ""
    msg[:float_data] = convert(Vector, radii(geometry))
    msg[:num_float_data] = 3
end

function fill_geometry_data!(msg::PyObject, geometry::HyperCylinder{3}, transform::Transformation)
    msg[:type] = msg[:CYLINDER]
    msg[:position] = transform(SVector{3, Float64}(0, 0, 0))
    msg[:quaternion] = to_lcm_quaternion(transform)
    msg[:string_data] = ""
    msg[:float_data] = [radius(geometry); length(geometry)]
    msg[:num_float_data] = 2
end

function to_lcm{T, GeomType}(geometry_data::GeometryData{T, GeomType})
    msg = drakevis[:lcmt_viewer_geometry_data]()
    msg[:color] = to_lcm(geometry_data.color)

    fill_geometry_data!(msg, geometry_data.geometry, geometry_data.transform)
    msg
end

function to_lcm(link::Link, robot_id_number::Real)
    msg = drakevis[:lcmt_viewer_link_data]()
    msg[:name] = link.name
    msg[:robot_num] = robot_id_number
    msg[:num_geom] = length(link.geometry_data)
    for geometry_data in link.geometry_data
        push!(msg["geom"], to_lcm(geometry_data))
    end
    msg
end

function to_lcm(robot::Robot, robot_id_number::Real)
    msg = drakevis[:lcmt_viewer_load_robot]()
    msg[:num_links] = length(robot.links)
    for link in robot.links
        push!(msg["link"], to_lcm(link, robot_id_number))
    end
    msg
end

type Visualizer
    robot::Robot
    robot_id_number::Int
    lcm::LCM

    function Visualizer(robot::Robot, robot_id_number::Integer=1, lcm::LCM=LCM(); load_now::Bool=true)
        vis = new(robot, robot_id_number, lcm)
        if load_now
            reload_model(vis)
        end
        vis
    end

end

Visualizer(data, robot_id_number::Integer=1, lcm::LCM=LCM()) = Visualizer(convert(Robot, data), robot_id_number, lcm)

function Visualizer(robots::Vector{Robot}, lcm::LCM=LCM())
    visualizers = Vector{Visualizer}()
    for (i, robot) in enumerate(robots)
        push!(visualizers, Visualizer(robot, i, lcm; load_now=false))
    end
    reload_model(visualizers)
    visualizers
end

function reload_model(vis::Visualizer)
    msg = to_lcm(vis.robot, vis.robot_id_number)
    publish(vis.lcm, "DRAKE_VIEWER_LOAD_ROBOT", msg)
end

function reload_model(visualizers::AbstractArray{Visualizer})
    msg = drakevis[:lcmt_viewer_load_robot]()
    msg[:num_links] = 0
    for vis in visualizers
        msg[:num_links] = msg[:num_links] + length(vis.robot.links)
        for link in vis.robot.links
            push!(msg["link"], to_lcm(link, vis.robot_id_number))
        end
    end
    publish(visualizers[1].lcm, "DRAKE_VIEWER_LOAD_ROBOT", msg)
end

function new_window()
    installed_visualizer_path = joinpath(Pkg.dir("DrakeVisualizer"), "deps", "usr", "bin", "drake-visualizer")
    if isfile(installed_visualizer_path)
        # If we built drake-visualizer, then use it
        (stream, proc) = open(`$installed_visualizer_path`)
    else
        # Otherwise let the system try to find it
        (stream, proc) = open(`drake-visualizer`)
    end
    proc
end

function draw{T <: Transformation}(model::Visualizer, link_origins::Vector{T})
    msg = drakevis[:lcmt_viewer_draw]()
    msg[:timestamp] = convert(Int64, time_ns())
    msg[:num_links] = length(link_origins)
    for (i, origin) in enumerate(link_origins)
        push!(msg["link_name"], model.robot.links[i].name)
        push!(msg["robot_num"], model.robot_id_number)
        push!(msg["position"], convert(SVector{3, Float64}, origin(SVector{3, Float64}(0,0,0))))
        # push!(msg["position"], origin.offset)
        push!(msg["quaternion"], to_lcm_quaternion(origin))
        # push!(msg["quaternion"], to_lcm(qrotation(rotationparameters(origin.scalefwd))))
    end
    publish(model.lcm, "DRAKE_VIEWER_DRAW", msg)
end

function __init__()
    lcmtypes_path = abspath(joinpath(Pkg.dir("DrakeVisualizer"), "src", "lcmtypes"))
    println("adding: $(lcmtypes_path) to the python path")
    unshift!(PyVector(pyimport("sys")["path"]), lcmtypes_path)
    copy!(drakevis, pyimport("drakevis"))
end

end


text:
#=
  Sum of 2 dice.

  I asked how one should observe a sum of the two dice here:
  https://discourse.julialang.org/t/turing-observing-sum-of-two-random-dice/50486

=#

using Turing, StatsPlots, DataFrames
include("../jl_utils.jl")

@model function two_dice(ss)
    d1 ~ DiscreteUniform(1, 6)
    d2 ~ DiscreteUniform(1, 6)

    ss ~ Dirac(d1 + d2)

    return ss
end

# Note: Here we observe c=5 from outside
model = two_dice(4)

# chains = sample(model, MH(), 10000) # This don't work anymore! c is 1..5!
num_chains = 4 # 4
# This works: i.e. forces d1 and d2 to be 2
chains = sample(model, MH(), MCMCThreads(), 10000, num_chains)

# Both d1 and d2 have means 3.5
# chains = sample(model, IS(), MCMCThreads(), 1000, num_chains) # Nope

# Both d1 and d2 have means 3.5
# chains = sample(model, PG(20), 1000) # Nope

# # Both d1 and d2 have means 3.5
# chains = sample(model, SMC(1000), 1000) # Nope

# chains = sample(model, NUTS(), 1000) # Error
# chains = sample(model, HMC(0.1, 5), 1000) # Error
# chains = sample(model, HMCDA(0.15, 0.65), 1000) # ERROR

display(chains)
# display(plot(chains))
# display(histogram(chains))
# display(gelmandiag(chains))

show_var_dist_pct(chains, :d1)
show_var_dist_pct(chains, :d2)

println("\nhash d1+d2")
display(sort(make_hash([chains[:d1][i]+chains[:d2][i] for i in 1:length(chains)])))

gen = generated_quantities(model, chains)
println("After removing the 0s from gen")
# gen2 = [g for g in gen if g != 0]
show_var_dist_pct(gen, 120)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
mutable struct NVector{T,S} <: AbstractVector{T}
    data::NTuple{S, T}
end
Base.size(xs::NVector) = (length(xs.data),)
Base.getindex(xs::NVector, i::Integer) = xs.data[i]
Base.setindex!(xs::NVector, v, i::Integer)  = xs.data = ntuple(j -> j == i ? v : xs.data[j], length(xs.data))

# Maps.
include("core/address_map.jl")

# Selections.
include("core/selections.jl")

# Visitor.
include("core/visitor.jl")

# Traces.
include("core/traces.jl")
Trace() = DynamicTrace()

# Learnables and gradients.
include("core/learnables.jl")


text:
# This file is a part of SimilaritySearch.jl

export L1Distance, L2Distance, SqL2Distance, LInftyDistance, LpDistance
import Distances: evaluate

"""
    L1Distance()

The manhattan distance or ``L_1`` is defined as

```math
L_1(u, v) = \\sum_i{|u_i - v_i|}
```
"""
struct L1Distance <: PreMetric end

"""
    L2Distance()

The euclidean distance or ``L_2`` is defined as

```math
L_2(u, v) = \\sqrt{\\sum_i{(u_i - v_i)^2}}
```
"""
struct L2Distance <: PreMetric end

"""
    SqL2Distance()

The squared euclidean distance is defined as

```math
L_2(u, v) = \\sum_i{(u_i - v_i)^2}
```

It avoids the computation of the square root and should be used
whenever you are able do it.
"""
struct SqL2Distance <: PreMetric end

"""
    LInftyDistance()

The Chebyshev or ``L_{\\infty}`` distance is defined as

```math
L_{\\infty}(u, v) = \\max_i{\\left| u_i - v_i \\right|}
```

"""
struct LInftyDistance <: PreMetric end

"""
    LpDistance(p)
    LpDistance(p, pinv)

The general Minkowski distance ``L_p`` distance is defined as

```math
L_p(u, v) = \\left|\\sum_i{(u_i - v_i)^p}\\right|^{1/p}
```

Where ``p_{inv} = 1/p``. Note that you can specify unrelated `p` and `pinv` if you need an specific behaviour.
"""
struct LpDistance <: PreMetric
    p::Float32
    pinv::Float32
end

LpDistance(p::Real) = LpDistance(p, 1/p)

"""
    evaluate(::L1Distance, a, b)

Computes the Manhattan's distance between `a` and `b`
"""
function evaluate(::L1Distance, a, b)
    d = zero(eltype(a))

    @fastmath @inbounds @simd for i in eachindex(a)
	    m = a[i] - b[i]
        d += ifelse(m > 0, m, -m)
    end

    d
end

"""
    evaluate(::L2Distance, a, b)

Computes the Euclidean's distance betweem `a` and `b`
"""
function evaluate(::L2Distance, a, b)
    d = zero(eltype(a))

    @fastmath @inbounds @simd for i in eachindex(a)
        d += (a[i] - b[i])^2 #m * m
    end

    sqrt(d)
end

"""
    evaluate(::SqL2Distance, a, b)

Computes the squared Euclidean's distance between `a` and `b`
"""
function evaluate(::SqL2Distance, a, b)
    d = zero(eltype(a))

    @fastmath @inbounds @simd for i in eachindex(a)
        d += (a[i] - b[i])^2
    end

    d
end

"""
    evaluate(::LInftyDistance, a, b)

Computes the maximum distance or Chebyshev's distance
"""
function evaluate(::LInftyDistance, a, b)
    d = zero(eltype(a))

    @inbounds @simd for i in eachindex(a)
        m = abs(a[i] - b[i])
        d = max(d, m)
    end

    d
end

"""
    evaluate(lp::LpDistance, a, b)

Computes generic Minkowski's distance
"""
function evaluate(lp::LpDistance, a, b)
    d = zero(eltype(a))

    @inbounds @simd for i in eachindex(a)
        m = abs(a[i] - b[i])
        d += m ^ lp.p
    end

    d ^ lp.pinv
end


text:
using BenchmarkTools, JLD

# Define a parent BenchmarkGroup to contain our suite
#const suite = BenchmarkGroup()
#suite["psMethods"] = BenchmarkGroup(["integrationScheme","Nck"])

paramspath = joinpath(dirname(@__FILE__), "params.jld")

loadparams!(suite, BenchmarkTools.load(paramspath, "suite"), :evals, :samples);
results = run(suite; verbose = true) # run all benchmarks
BenchmarkTools.save("pr.jld", "results", results)   # save results to JLD file

master = load("master.jld", "results")
pr = load("pr.jld", "results")
regs = regressions(judge(minimum(pr), minimum(master))) # a BenchmarkGroup containing the regressions
pairs = leaves(regs) # an array of (ID, `TrialJudgement`) pairs

m1=median(suite["psMethods"])


text:
module SharedDistributedCode

using Distributed
using Test
using XUnit
using XUnit: create_deep_copy, has_wrapped_exception

function do_work(jobs, results) # define work function everywhere
    try
        tls = task_local_storage()
        if Main.GLOBAL_HAS_XUNIT_STATE
            xs = create_deep_copy(Main.GLOBAL_XUNIT_STATE)
            empty!(xs.test_suites_stack)
            empty!(xs.stack)
            empty!(xs.seen)
            tls[:__XUNIT_STATE__] = xs
        end
        attempt_cnt = 1
        # we wait until either tests get scheduled (indicated by `Main.__SCHEDULED_DISTRIBUTED_TESTS__`)
        # or at least tried to get scheduled (indicated by `Main.__RUN_DISTRIBUTED_TESTS_CALLED__`)
        while !isdefined(Main, :__SCHEDULED_DISTRIBUTED_TESTS__) &&
              !isdefined(Main, :__RUN_DISTRIBUTED_TESTS_CALLED__)
            @info "Main.__SCHEDULED_DISTRIBUTED_TESTS__ IS NOT still defined on process $(myid())."
            @info "Sleeping for a second (attempt #$(attempt_cnt))."
            sleep(1)
            attempt_cnt += 1
        end
        # If the test scheduling code encounters an error, then `Main.__RUN_DISTRIBUTED_TESTS_CALLED__`
        # is defined but `Main.__SCHEDULED_DISTRIBUTED_TESTS__` won't be available.
        @assert isdefined(Main, :__SCHEDULED_DISTRIBUTED_TESTS__) "Main.__SCHEDULED_DISTRIBUTED_TESTS__ IS NOT defined on process $(myid())"
        scheduled_tests = Main.__SCHEDULED_DISTRIBUTED_TESTS__

        task_count = 0
        while true
            (scheduled_tests_index, scheduled_test_name) = take!(jobs)
            scheduled_tests_index < 1 && break
            if scheduled_tests_index > length(scheduled_tests)
                throw("scheduled_tests_index ($scheduled_tests_index) is outside of bound for scheduled_tests (with $(length(scheduled_tests)) elements)")
            end

            task_count += 1

            println("Process $(myid()) is handling task #$(task_count) (which is $(scheduled_tests_index)/$(length(scheduled_tests)))")

            st = scheduled_tests[scheduled_tests_index]
            try
                if get_description(st) != scheduled_test_name
                    println("A critical test scheduling error:")
                    println("  get_description(st) (\"$(get_description(st))\") != scheduled_test_name (\"$(scheduled_test_name)\")")
                    println("List of scheduled tests on worker #$(myid()):")
                    for (i, tst) in enumerate(scheduled_tests)
                        println("$i ==> $(get_description(tst))")
                    end
                    put!(jobs, (scheduled_tests_index, scheduled_test_name))
                    break
                end

                if XUnit.TESTSET_PRINT_ENABLE[]
                    path = XUnit._get_path(vcat(st.parent_testsets, [st.target_testcase]))
                    std_io = IOBuffer()
                    print(std_io, "~> Running ")
                    printstyled(std_io, path; bold=true)
                    println(std_io, string(" test-case (on pid=", myid(), ")..."))
                    seekstart(std_io)
                    # thread-safe print
                    print(read(std_io, String))
                end

                XUnit.run_single_testcase(st.parent_testsets, st.target_testcase)

                put!(results, (
                    scheduled_tests_index,
                    XUnit.DistributedAsyncTestMessage(st.target_testcase),
                    myid(),
                ))
            catch err
                has_wrapped_exception(err, InterruptException) && rethrow()

                println("A critical error occued while running '$scheduled_test_name': ", err)
                for s in stacktrace(catch_backtrace())
                    println(s)
                end

                # something in the test block threw an error. Count that as an
                # error in this test set
                ts = st.target_testcase.testset_report.reporting_test_set[]
                XUnit.record(
                    ts,
                    XUnit.Error(
                        :nontest_error,
                        Expr(:tuple),
                        err,
                        Base.catch_stack(),
                        st.target_testcase.source
                    )
                )

                put!(results, (
                    scheduled_tests_index,
                    XUnit.DistributedAsyncTestMessage(st.target_testcase),
                    myid(),
                ))
            end
        end
        println("Process $(myid()) is done with handling task after running $(task_count) tasks (out of $(length(scheduled_tests)))")
    catch err
        has_wrapped_exception(err, InterruptException) && rethrow()

        println("A critical error occurred in XUnit while running tests: ", err)
        for s in stacktrace(catch_backtrace())
            println(s)
        end

        # cancel all remaining tests by putting an empty result with index `0`
        put!(results, (0, nothing, myid()))
    end
end

end


text:
using DataFrames, StatsBase

mutable struct problem
    data::Array{Float64}
    Sigma::Array{Float64}
end

# return obejctive and solution given a solution vector
function evaluate(solution, prob)
	return (transpose(solution)*prob.Sigma*solution)[1], solution
end

# return objective and solution given a support vector (of 0s and 1s)
function evalSupport(support, prob)
	xVal = LinearAlgebra.eig(Qhat(prob.Sigma,support))[2][:,size(prob.Sigma,1)]
	return evaluate(xVal, prob)
end

# return objective and solution given the indices of a support
function evalIndices(support, prob)
	y = round.(Int64,zeros(size(prob.Sigma,1)))
	y[support[support.!=0]] = 1
	xVal = LinearAlgebra.eig(Qhat(prob.Sigma,y))[2][:,size(prob.Sigma,1)]
	return evaluate(xVal, prob)
end

# routine for creating a new problem out of a large data set by selecting
# a subset of dimensions out of all the variables in the problem
function shrinkProblem(prob, dimensions)
	finished = false
	problem2 = problem(prob.data, prob.Sigma);
	while !finished
		problem2 = problem(prob.data, prob.Sigma);
		r = problem2.data
		if size(r)[2]>dimensions
			selection = sample(1:size(r)[2] ,dimensions, replace=false)
			r = r[:, selection]
		end
		r = r .- Transpose(mean.(eachcol(r)))
		problem2.data = r
		problem2.Sigma = r'*r/size(r)[1];
		# finished = isposdef(problem2.Sigma) # removing so we can go higher than 250
		finished = true
	end
	return problem2
end

# zeroes out rows and columns from mat that correspond to support = 0
function Qhat(mat, support)
	yDel = (support.==0)
	Q = copy(mat)
	Q[yDel , : ] = 0
	Q[ : , yDel] = 0
	return Q
end

#Faster computation of maximum eigenvalue and eigenvector
function myeigmax(SigOrA, thisStart, y, highDim, refineCap)
	yKeep = .!(y.==0)
	m=0
	if highDim
		thisA = copy(SigOrA[:, yKeep])
		m, n = size(thisA)
		Q = thisA*thisA'/(m-1)
		b = thisStart[yKeep]
		b = thisA*(b / LinearAlgebra.LinearAlgebra.norm(b))
		normb = 0
		newnorm = LinearAlgebra.norm(b)
	else
		Q = copy(SigOrA[yKeep,yKeep])
		beta0 = thisStart[yKeep]
		normb = LinearAlgebra.norm(beta0)
		b = Q*beta0
		newnorm = LinearAlgebra.norm(b)
	end

	cycles = 0
	while abs(normb - newnorm) > .0001*normb && cycles < 100
		normb = newnorm
		b = Q*(b / normb)
		newnorm = LinearAlgebra.norm(b)
		if refineCap>0
			if newnorm > refineCap
				return refineCap, copy(b/normb)
			end
		end
		cycles = cycles + 1
	end

	if highDim
		b = thisA'*b
		b = thisA'*(thisA*b/LinearAlgebra.norm(b))/(m-1)
		normb = LinearAlgebra.norm(b)
	end

	b = b/normb
	expandedb = zeros(length(y))
	expandedb[yKeep]=b

	return normb, expandedb
end

#Creates a sparse version of a vector by killing smallest components and scaling up
function Hk(origlist, sparsity, support)
	list = real(copy(origlist))
	ksparse = zeros(length(list))
	indicesToKeep = (support.==1)
	dummyvalue = minimum(list)-1
	list[(support.>-1)] .= dummyvalue

	newIndices = selectperm2(list, sparsity-sum(indicesToKeep))
	indicesToKeep[newIndices].=true

	ksparse[indicesToKeep]=origlist[indicesToKeep]

	ksparse = ksparse / LinearAlgebra.norm(ksparse[indicesToKeep])

	return ksparse
end

# selects *indices* of x that correspond to the k largest values of x
# for use in Hk utility
function selectperm2(x,k)
    if k > 1
        kk = 1:k
    else
        kk = 1
    end
    z = collect(1:length(x))
    return partialsort!(z,1:k,by = (i)->abs(x[i]), rev = true)
end

# selects out of x the k largest values
# for use in B+B algorithm
function selectsorted(x,k)
    if k > 1
        kk = 1:k
    else
        kk = 1
    end
    z = collect(1:length(x))
    return partialsort!(x,1:k,by = (i)->(i), rev = true)
end

# partitions a vector of length x into n parts
yourpart(x,n) = Any[Any[x[i:n:length(x)]] for i in 1:n]

# faster matrix vector multiplication when the vector is sparse
function sparseTimes(mat, vec)
    ind = abs.(vec).>.000001
    prod = mat[:,ind]*vec[ind]
    return prod
end

# Yuan and zhang algorithm with random restarts
function subset(prob, k; timeLimit = 7200, support = zeros(1), countdown = 100)
	n = size(prob.Sigma,1)
	lambdas, betas, =Arpack.eigs(prob.Sigma, nev=1, which=:LR)
	beta0=betas[:,1]
	if length(support)==1
		support = zeros(n).-1
	end

	#Starts with local search starting at first eigenvector
	bestBeta = eigSubset(prob, support, k, beta0)
	bestObj, ~ = evaluate(bestBeta, prob)
	start = time()
	#Tries a thousand random other starting points
	margin = 1
	while countdown > 0 && time()-start < timeLimit
		beta = rand(size(prob.Sigma,1))
		beta = beta / LinearAlgebra.norm(beta0)
		beta = eigSubset(prob, support, k, beta)
		obj, ~ = evaluate(beta, prob)
		if obj > bestObj
			bestObj = obj
			bestBeta = copy(beta)
			countdown = 100
		end
		countdown = countdown - 1
	end

	return bestObj, bestBeta
end

#Applies Yuan and Zhang heuristic starting from first eigenvector
function eigSubset(prob, support, k, beta0)
	beta = Hk(beta0, k, support)
	for i=1:100
		beta = Hk(prob.Sigma*beta, k, support)
	end
	return beta
end

# Creates a new file at filename with specified header
function refreshFile(filename, header)
	outfile = open(filename, "w");
	close(outfile);
	outfile = open(filename, "a");
	write(outfile,header);
	close(outfile);
end


text:
include("containers.jl")
include("math.jl")

# Enable debug timings `@timed timer() "name" stuff...`.
# This allows us to disable timings completely by executing
# `TimerOutputs.disable_debug_timings(Trixi)`
# and to enable them again by executing
# `TimerOutputs.enable_debug_timings(Trixi)`
timeit_debug_enabled() = true

# Store main timer for global timing of functions
const main_timer = TimerOutput()

# Always call timer() to hide implementation details
timer() = main_timer

"""
    PerformanceCounter

A `PerformanceCounter` be used to track the runtime performance of some calls.
Add a new runtime measurement via `put!(counter, runtime)` and get the averaged
runtime of all measurements added so far via `take!(counter)`, resetting the
`counter`.
"""
mutable struct PerformanceCounter
  ncalls_since_readout::Int
  runtime::Float64
end

PerformanceCounter() = PerformanceCounter(0, 0.0)

function Base.take!(counter::PerformanceCounter)
  time_per_call = counter.runtime / counter.ncalls_since_readout
  counter.ncalls_since_readout = 0
  counter.runtime = 0.0
  return time_per_call
end

function Base.put!(counter::PerformanceCounter, runtime::Real)
  counter.ncalls_since_readout += 1
  counter.runtime += runtime
end

"""
    examples_dir()

Return the directory where the example files provided with Trixi.jl are located. If Trixi is
installed as a regular package (with `]add Trixi`), these files are read-only and should *not* be
modified. To find out which files are available, use, e.g., `readdir`:

# Examples
```julia
julia> readdir(examples_dir())
5-element Array{String,1}:
 "1d"
 "2d"
 "3d"
 "README.md"
 "paper-self-gravitating-gas-dynamics"
```
"""
examples_dir() = joinpath(pathof(Trixi) |> dirname |> dirname, "examples")

"""
    get_examples()

Return a list of all example elixirs that are provided by Trixi.
"""
function get_examples()
  examples = String[]
  for (root, dirs, files) in walkdir(examples_dir())
    for f in files
      if startswith(f, "elixir_") && endswith(f, ".jl")
        push!(examples, joinpath(root, f))
      end
    end
  end

  return examples
end

"""
    default_example()

Return the path to an example elixir that can be used to quickly see Trixi in action on a
[`TreeMesh`]@(ref).
"""
default_example() = joinpath(examples_dir(), "2d", "elixir_advection_basic.jl")

"""
    default_example_unstructured()

Return the path to an example elixir that can be used to quickly see Trixi in action on an
[`UnstructuredQuadMesh`]@(ref). This simulation is run on the example curved, unstructured mesh
given in the Trixi documentation regarding unstructured meshes.
"""
default_example_unstructured() = joinpath(examples_dir(), "2d", "elixir_euler_unstructured_quad_basic.jl")

# Print informative message at startup
function print_startup_message()
  s = """

    ████████╗██████╗ ██╗██╗  ██╗██╗
    ╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
       ██║   ██████╔╝██║ ╚███╔╝ ██║
       ██║   ██╔══██╗██║ ██╔██╗ ██║
       ██║   ██║  ██║██║██╔╝ ██╗██║
       ╚═╝   ╚═╝  ╚═╝╚═╝╚═╝  ╚═╝╚═╝
    """
  mpi_println(s)
end

"""
    get_name(x)

Returns a name of `x` ready for pretty printing.
By default, return `string(y)` if `x isa Val{y}` and return `string(x)` otherwise.

# Examples

```jldoctest
julia> Trixi.get_name("test")
"test"

julia> Trixi.get_name(Val(:test))
"test"
```
"""
get_name(x) = string(x)
get_name(::Val{x}) where x = string(x)

"""
    @threaded for ... end

Semantically the same as `Threads.@threads` when iterating over a `AbstractUnitRange`
but without guarantee that the underlying implementation uses `Threads.@threads`
or works for more general `for` loops.
In particular, there may be an additional check whether only one thread is used
to reduce the overhead of serial execution or the underlying threading capabilities
might be provided by other packages such as [Polyester.jl](https://github.com/JuliaSIMD/Polyester.jl).

!!! warn
    This macro does not necessarily work for general `for` loops. For example,
    it does not necessarily support general iterables such as `eachline(filename)`.

Some discussion can be found at https://discourse.julialang.org/t/overhead-of-threads-threads/53964
and https://discourse.julialang.org/t/threads-threads-with-one-thread-how-to-remove-the-overhead/58435.
"""
macro threaded(expr)
  # Use `esc(quote ... end)` for nested macro calls as suggested in
  # https://github.com/JuliaLang/julia/issues/23221
  #
  # The following code is a simple version using only `Threads.@threads` from the
  # standard library with an additional check whether only a single thread is used
  # to reduce some overhead (and allocations) for serial execution.
  #
  # return esc(quote
  #   let
  #     if Threads.nthreads() == 1
  #       $(expr)
  #     else
  #       Threads.@threads $(expr)
  #     end
  #   end
  # end)
  #
  # However, the code below using `@batch` from Polyester.jl is more efficient,
  # since this packages provides threads with less overhead. Since it is written
  # by Chris Elrod, the author of LoopVectorization.jl, we expect this package
  # to provide the most efficient and useful implementation of threads (as we use
  # them) available in Julia.
  # !!! danger "Heisenbug"
  #     Look at the comments for `wrap_array` when considering to change this macro.

  return esc(quote @batch $(expr) end)
end

#     @timed timer() "some label" expression
#
# Basically the same as a special case of `@timeit_debug` from
# [TimerOutputs.jl](https://github.com/KristofferC/TimerOutputs.jl),
# but without `try ... finally ... end` block. Thus, it's not exception-safe,
# but it also avoids some related performance problems. Since we do not use
# exception handling in Trixi, that's not really an issue.
macro timed(timer_output, label, expr)
  timeit_block = quote
    if timeit_debug_enabled()
      local to = $(esc(timer_output))
      local enabled = to.enabled
      if enabled
        local accumulated_data = $(TimerOutputs.push!)(to, $(esc(label)))
      end
      local b₀ = $(TimerOutputs.gc_bytes)()
      local t₀ = $(TimerOutputs.time_ns)()
    end
    local val = $(esc(expr))
    if timeit_debug_enabled() && enabled
      $(TimerOutputs.do_accumulate!)(accumulated_data, t₀, b₀)
      $(TimerOutputs.pop!)(to)
    end
    val
  end
end


text:
using SuperResolution
using Test

@testset "fctSR.jl" begin
    @test true
# # function dim(A)
# # 	d1=length(A[:,1])
# # 	d2=length(A[1,:])
# # 	return d1,d2
# # end
#     @test dim(zeros(Int64,3,3))==(3,3)
#     @test dim(zeros(Int64,3,3))[1]==3
#     @test dim(zeros(Int64,3,3))[2]==3
# # function eye(d::Int64)
# # 	E=zeros(Float64,d,d)
# # 	map(i->E[i,i]=1,1:dim(E)[1])
# # 	return E
# # end
#     @test eye(3)*eye(3)==eye(3)
# # function ctrans(X::Array{Float64,2})
# # 	convert(Array{Float64,2},transpose(X))
# # end
# # function ctrans(X::Array{Array{Float64,2},1})
# # 	convert(Array{Array,1},map(i->ctrans(𝛍[1][i])[:],collect(1:length(𝛍[1]))))
# # end
#     @test ctrans([[1 2];[3 4]])==[[1 3];[2 4]]
# # function dim2clean(MU_star::Array{Array,1})
# # 	MU_star = dim2adj(MU_star)
# # 	N=maximum(size(MU_star))
# # 	𝛅=Int64(3)
# # 	for i=1:N
# # 		if isnan(MU_star[i])
# # 			for j=1:𝛅
# # 				MU_star[i,j]=-10.0
# # 			end
# # 		end
# # 	end
# # 	MU_star = dim2adj(MU_star)
# # 	MU_star = MU_star[MU_star .> dim2adj([100.0,100.0,100.0])]
# # end
# # function dim2clean(MU_star::Array{Array,1},𝛆::Array{Float64,1})
# # 	MU_star = dim2adj(MU_star)
# # 	N=maximum(size(MU_star))
# # 	𝛅=Int64(3)
# # 	for i=1:N
# # 		if isnan(MU_star[i])
# # 			for j=1:𝛅
# # 				MU_star[i,j]=-10.0
# # 			end
# # 		end
# # 	end
# # 	MU_star = dim2adj(MU_star)
# # 	MU_star = MU_star[MU_star .>= dim2adj(𝛆)]
# # end
# # function dim2crop(X::Array{Float64,2},xmin::Float64,xmax::Float64,ymin::Float64,ymax::Float64)
# # 	Xog = X
# # 	X=X[X[:,1].>=xmin,:]
# # 	X=X[X[:,1].<=xmax,:]
# # 	X=X[X[:,2].>=ymin,:]
# # 	X=X[X[:,2].<=ymax,:]
# # 	val=(Xog[:,1] .∈ [X[:,1]]) .& (Xog[:,2] .∈ [X[:,2]])
# # 	return[X,val]
# # end
# #
# # function dim2adj(X::Array{Float64,2})
# # 	J = size(X,1)
# # 	# if size(X,1)<size(X,2)
# # 	# 	J = size(X,2)
# # 	# 	X = X'
# # 	# end
# # 	MU_0 = Array{Array,1}(undef,J)
# # 	for j=1:J
# # 		MU_0[j] = X[j,:]
# # 	end
# # 	return MU_0
# # end
# # function dim2adj(X::Array{Float64,1})
# # 	MU_0 = Array{Array,1}(undef,1)
# # 	MU_0[1] = X[:,1]
# # 	return MU_0
# # end
# # function dim2adj(X::Array{Array,1},δ=Int64(2))
# # 	if δ != Int64(2)
# # 		Y_raw = Array{Float64,2}(undef,maximum(size(X)),size(X[1],1))
# # 		for i=1:maximum(size(X))
# # 			for j=1:size(X[1],1)
# # 				Y_raw[i,j] = X[i][j]
# # 			end
# # 		end
# # 		return Y_raw
# # 	end
# # 	Y_raw = Array{Float64,2}(undef,maximum(size(X)),2)
# # 	for i=1:maximum(size(X))
# # 		Y_raw[i,1] = X[i][1]
# # 		Y_raw[i,2] = X[i][2]
# # 	end
# # 	return Y_raw
# # end
# #
# # ####Generate MixtureModel off given data
# # function dim2mix(X::Array{Array,1},S=eye(size(X,2)))
# # 	MM = MixtureModel(map(x->MvNormal(x,S),X))
# # 	return MM
# # end
# # function dim2mix(X::Array{Float64,2},S=eye(size(X,2)))
# # 	X = dim2adj(X)
# # 	MM = MixtureModel(map(x->MvNormal(x,S),X))
# # 	return MM
# # end
# # function dim2mix(X::Array{Float64,2},w::Array{Float64,1},S=eye(size(X,2)))
# # 	#print("[dim2mix] X=",X,"\n")
# # 	X = dim2adj(X)
# # 	#print("[dim2mix..dim2adj(X)] X=",X,"\n")
# # 	#print("[dim2mix] S=",S,"\n")
# # 	MM = MixtureModel(map(x->MvNormal(x,S),X))
# # 	return MM
# # end
# # function dim2mix(X::Array{Array,1},Σ::Array{Array{Float64,2}})
# # 	MM = MixtureModel(map(i->MvNormal(X[i],Σ[i]),1:length(X)))
# # 	return MM
# # end
# # function dim2mix(X::Array{Array,1},Σ::Array{Float64,2},p::Array{Float64,1})
# # 	MM = MixtureModel(map(x->MvNormal(x,Σ),X),p)
# # 	return MM
# # end
# # function dim2mix(X::Array{Float64,2},Σ::Array{Float64,2})
# # 	X = dim2adj(X)
# # 	MM = MixtureModel(map(x->MvNormal(x,Σ),X))
# # 	return MM
# # end
# # ####Evaluate bounds of data and generate approriate range
# # ####in bounds for package KernelDensity
# # function dim2bounds(X::Array{Float64,2})
# # 	M = 0*Matrix{Float64}(undef,2,2)
# # 	M[1,1] = minimum(X[:,1])
# # 	M[2,1] = maximum(X[:,1])
# # 	M[1,2] = minimum(X[:,2])
# # 	M[2,2] = maximum(X[:,2])
# # 	return M
# # end
# # function dim2range(strt::Float64,stp::Float64,len::Int64)
# # 	Omega = range(strt,stop=stp,length=Int64(len))
# # 	return Omega
# # end
# # ####Generate appropriate range in bounds for data X_raw in dimension d
# # ####The range has to be generated for each dimension
# # function dim2range(X_raw::Array{Float64,2},d::Int64)
# # 	if size(X_raw,1)<size(X_raw,2)
# # 		X_raw=ctrans(X_raw)
# # 	end
# # 	len = Int64(floor(sqrt(maximum(size(X_raw)))))
# # 	M = dim2bounds(X_raw)
# # 	Omega = range(M[1,d],stop=M[2,d],length=Int64(len))
# # 	return Omega
# # end
# # function dim2range(X_raw::Array{Float64,2})
# # 	if size(X_raw,1)<size(X_raw,2)
# # 		X_raw=ctrans(X_raw)
# # 	end
# # 	len = Int64(floor(sqrt(maximum(size(X_raw)))))
# # 	M = dim2bounds(X_raw)
# # 	Omega1 = range(M[1,1],stop=M[2,1],length=Int64(len))
# # 	Omega2 = range(M[1,2],stop=M[2,2],length=Int64(len))
# # 	delta = [Omega1[end]-Omega1[1];Omega2[end]-Omega2[1]]
# # 	return [Omega1,Omega2,delta]
# # end
# # #MLE geom dis given FrNm
# # function geomMLE(FrNm::Array{Int64,1})
# # 	β=collect(1e-8:1e-2:1.0)
# # 	βpost=[]
# # 	#Likelihood for β parameter of Geometric Distribution, based on Framenmbers of Observations
# # 	FrNm = FrRef(FrNm)
# # 	for j=1:length(β)
# # 		push!(βpost,sum([FrNm[i]*log(1-β[j])+log(β[j])] for i in 1:length(FrNm))[1])
# # 	end
# # 	#MLE for β Parameter
# # 	βhat=β[collect(1:length(β))[βpost.==maximum(βpost)]][1]
# # 	#Geometric Distribution with MLE for Parameter
# # 	G=Geometric(βhat)
# # 	return G
# # end
# # function normalMLE(FrNm::Array{Int64,1})
# # 	μMLE=1/length(FrNm) * sum(FrNm)
# # 	σ=1.0
# # 	nMLE=Normal(μMLE,σ)
# # 	# σMLE=sqrt(1/length(FrNm) * sum( (FrNm .- μMLE).^2 ) )
# # 	# if sum(FrNm .== μMLE) == length(FrNm) σMLE = 1.0 end #avoid 0 value for σMLE (occurs for all FrNm==μMLE)
# # 	#nMLE=Normal(μMLE,σMLE)
# # 	return nMLE
# # end
# # ##ref FrNm of obs to first frame
# # function FrRef(FrNm::Array{Int64,1})
# # 	FrNm = sort(FrNm) .- minimum(FrNm)
# # end
# # #
# # function FrRef(FrNm::Int64)
# # 	return FrNm
# # end
# # #compute data based on cluster corresponding assignment in AcollElem
# # function Feval(AcollElem::Array{Int64,1})#,𝔽sim::Array{Int64,1})
# # 	uniqueLengthElem=length(unique(AcollElem))
# # 	uniqueAcollElem=unique(AcollElem)
# # 	#get framenumbers for each cluster
# # 	𝔽evalElem=Array{Array,1}(undef,uniqueLengthElem)
# # 	𝔽evalElem=map(i->𝔽sim[AcollElem .== unique(AcollElem)[i]],1:uniqueLengthElem)
# # 	#eval indicies of observations (to evaluate locations in 𝕐sim)
# # 	𝔽evalElemIdx=Array{Array,1}(undef,uniqueLengthElem)
# # 	𝔽evalElemIdx=map(i->findall(x->x∈𝔽evalElem[i],𝔽sim),1:uniqueLengthElem)
# # 	#get observations corresponding to clusters
# # 	𝕐evalElem=Array{Array,1}(undef,uniqueLengthElem)
# # 	𝕐evalElem=map(i->𝕐sim[𝔽evalElemIdx[i],:],1:uniqueLengthElem)
# # 	#FrmNm f.e. Cl ID||Indices of obs f.e. Cl ID||Loc of obs f.e. Cl ID||number of unique Cl||ID of unique Cl
# # 	return 𝔽evalElem,𝔽evalElemIdx,𝕐evalElem,uniqueLengthElem,uniqueAcollElem
# # end
# # #Kullback Leibler
# # function mu_star(mu1::Array{Float64,1},mu2::Array{Float64,1},W::Array{Float64,1})
# # 	mu_star = (W[1]*mu1+W[2]*mu2)/(W[1]+W[2])
# # 	return mu_star
# # end
# # function s_star(mu1::Array{Float64,1},mu2::Array{Float64,1},S1::Array{Float64,2},S2::Array{Float64,2},W::Array{Float64,1})
# # 	Σ_star = (W[1]*S1 + W[2]*S2)/(W[1]+W[2]) + (W[1]*W[2]*(mu1-mu2)*(mu1-mu2)')/((W[1]+W[2])^2)
# # 	return Σ_star
# # end
# # function dim2KLD(mu1::Array{Float64,1},mu2::Array{Float64,1},S1::Array{Float64,2},S2::Array{Float64,2},W::Array{Float64,1})
# # 	S=[S1,S2]
# # 	Σ_star = s_star(mu1,mu2,S[1],S[2],[W[1],W[2]])
# # 	𝛅=Int64(3)
# # 	KLD = (sum([W[i]/2.0 * (tr(Σ_star^-1 * S[i]) + log(det(Σ_star * S[i]^-1))-𝛅)] for i in 1:2)[1] + (W[1]*W[2]/(2.0*(W[1]+W[2])^2) * (mu1-mu2)'*Σ_star^-1*(mu1-mu2)))
# # 	if isnan(KLD) KLD=0.0 end
# # 	return KLD
# # end
# # #Extrema of KLD matrix
# # function dim2KLDex(KLDfull::Array{Float64,2})
# # 	𝛕=1e-5
# # 	EX = 0*Array{Int64,2}(undef,3,2)
# # 	K = maximum(size(KLDfull))
# # 	#column1 minimum KLD indices [vectorindex, matrixindices[j,k]]
# # 	EX[2,1] = findall(KLDfull .== minimum(KLDfull[findall((KLDfull .> 𝛕))]))[1][1] #column index
# # 	EX[3,1] = findall(KLDfull .== minimum(KLDfull[findall((KLDfull .> 𝛕))]))[1][2] #row index
# # 	EX[1,1] = Int64(floor((EX[3,1]-1)*K+EX[2,1])) #single indexing
# # 	EXmin = KLDfull[EX[1,1]]
# # 	#column2 maximum KLD indices [vectorindex, matrixindices[j,k]]
# # 	EX[2,2] = findall(KLDfull .== maximum(KLDfull[findall((KLDfull .> 𝛕))]))[1][1]
# # 	EX[3,2] = findall(KLDfull .== maximum(KLDfull[findall((KLDfull .> 𝛕))]))[1][2]#Int64(EX[1,2]-K*(EX[3,2]-1))
# # 	EX[1,2] = Int64(floor((EX[3,2]-1)*K+EX[2,2]))
# # 	EXmean=mean(collect(KLDfull[:]))
# # 	#print("Idx min: ",EX[1,1],"= ","[",EX[2,1],",",EX[3,1],"]","min value: ", EXmin,"(mean value: ",EXmean,")\n")
# # 	return EX
# # end
# # ##Weighted Kullback Leibler Divergence
# # function dim2WKLDfull(MU::Array{Array,1},S::Array{Array{Float64,2},1},W::Array{Float64,1},A=NaN,Atmp=NaN)
# # 	#entries are merged during each collapse, therefore K represents the current number of clusters
# # 	K=maximum(size(MU))
# # 	#if no A or Atmp is submitted, it's set as NaN (though A and Atmp are by default submitted during dim2WKLDfull computation)
# # 	if isnan(A[1]) A=collect(1:K) end
# # 	if isnan(Atmp[1]) Atmp=collect(1:K) end
# # 	KLD = 1e3*ones(K,K)
# # 	global MU_star = map(i->MU[i],1:length(MU)) #(clustered-)data based on 𝕐sim
# # 	global A_star = map(i->A[i],1:length(A)) #cluster ID
# # 	#S_star = S
# # 	W_star = map(i->W[i],1:length(W)) #weights
# # 	A_coll = map(i->Atmp[i],1:length(Atmp)) #A_coll stores the cluster IDs which are updated for every collapse
# # 	𝛅 = Int64(3) #dimension
# # 	S_star= repeat([1.0*eye(𝛅)],length(S)) #COV Matrices
# # 	FrNm = Feval(A_coll)[1] #evaluate frames f.e. clust ID in A_coll
# # 	#compute normal distribution with μMLE and σMLE for each current cluster
# # 	#NormDistClustID = map(i->normalMLE(FrNm[i]),1:length(FrNm))
# # 	GeomDistClustID = map(i->geomMLE(FrRef(FrNm[i])),1:length(FrNm))
# # 	#KL distances
# # 	K = maximum(size(MU_star))
# # 	#print("K=",K,"\n")
# # 	#initialize "distance" matrix
# # 	KLD = 1e-15*ones(K,K)
# # 	#compute KL-distances for each and every pair of points
# # 	for k=1:K
# # 		for j=1:K
# # 			if j==k
# # 				KLD[j,k]=1e5
# # 			else
# # 				#KLD[j,k] = (2-sum(pdf.(NormDistClustID[k],FrNm[j]))) * dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0])
# # 				KLD[j,k] = (1/exp(sum(pdf.(GeomDistClustID[k],FrNm[j])))) * dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0])
# # 			end
# # 		end
# # 	end
# # 	#dim2KLDex computes extrema and corresponding indices in the distance matrix
# # 	KLD[map(i->isnan(KLD[i]),1:length(KLD[:]))].= maximum(KLD[map(i->!isnan(KLD[i]),1:length(KLD[:]))])
# # 	EX = dim2KLDex(KLD)
# # 	#min1,min2 is the pair of values with minimal KLD and get merged
# # 	min1=EX[2,1]
# # 	min2=EX[3,1]
# # 	global MU_star[EX[2,1]] = mu_star(MU_star[EX[2,1]],MU_star[EX[3,1]],[W_star[EX[2,1]],W_star[EX[3,1]]])
# # 	S_star[EX[2,1]] = s_star(MU_star[EX[2,1]],MU_star[EX[3,1]],S_star[EX[2,1]],S_star[EX[3,1]],[W_star[EX[2,1]],W_star[EX[3,1]]])
# # 	W_star[EX[2,1]] += W_star[EX[3,1]]
# # 	S_star[EX[3,1]] = NaN.*S_star[EX[3,1]]
# # 	global MU_star[EX[3,1]] =NaN.*MU_star[EX[3,1]]
# # 	W_star[EX[3,1]] = NaN.*W_star[EX[3,1]]
# # 	A_coll[A_coll .== A_star[EX[3,1]]].= A_star[EX[2,1]]
# # 	#global A_star[EX[3,1]] = 0*A_star[EX[3,1]]
# # 	global A_star[EX[3,1]] = 0*A_star[EX[3,1]]
# # 	global A_star=A_star[map(i->A_star[i]!=0,1:length(A_star))]
# # 	S_star=S_star[map(i->!isnan(S_star[i][1]),1:length(S_star))]
# # 	global MU_star=MU_star[map(i->!isnan(MU_star[i][1]),1:length(MU_star))]
# # 	W_star=W_star[map(i->!isnan(W_star[i]),1:length(W_star))]
# # 	S_star=repeat([1.0*eye(𝛅)],length(S_star))
# # 	KLD[EX[1,1]]=0.0
# # 	#print("END OF while loop dim2WKLDfull\n")
# # 	return MU_star,W_star,S_star,A_star,A_coll,KLD
# # end
end

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
import Base: sin, cos, sincos, tan, sinh, cosh, tanh, asin, acos, atan, asinh, acosh,
             atanh, sec, csc, cot, asec, acsc, acot, sech, csch, coth, asech, acsch, acoth,
             sinpi, cospi, sinc, cosc, cosd, cotd, cscd, secd, sind, tand, acosd, acotd,
             acscd, asecd, asind, atand,
             cbrt,
             log, log2, log10, log1p,
             exp, frexp, expm1, exp2, exponent

import Base: prevfloat, floatmax, floatmin, maxintfloat, significand

float_one2one = (:sin, :cos, :tan, :sinh, :cosh, :tanh, :asin, :acos, :atan, :asinh,
                 :acosh, :atanh, :sec, :csc, :cot, :asec, :acsc, :acot, :sech, :csch,
                 :coth, :asech, :acsch, :acoth, :sinpi, :cospi, :sinc, :cosc, :cosd, :cotd,
                 :cscd, :secd, :sind, :tand, :acosd, :acotd, :acscd, :asecd, :asind,
                 :atand, :cbrt, :exp, :frexp, :expm1, :exp2, :significand)

function deffloat(::Type{ST}, ::Type{BT}) where {ST,BT<:BaseFloat}
    for f in float_one2one
        @eval begin
            $(f)(::$ST{X}) where X = $ST{$(f)(X::$BT)}()
        end
    end

    @eval begin
        Base.rem2pi(::$ST{X}, r::RoundingMode) where X = $ST{rem2pi(X::$BT, r)}()

        @generated function Base.sqrt(::$ST{X}) where X
            x = $ST{sqrt(X::$BT)::$BT}()
            :($x)
        end
#=
        @generated function Base.rat(::$ST{X}) where X
            n, d = Base.rat(X::$BT)
            :($(SInt{n}()), $(SInt{d}()))
        end
=#
        @generated function log(::$ST{X}) where X
            x = $ST{log(X::$BT)::$BT}()
            :($x)
        end

        @generated function log2(::$ST{X}) where X
            x = $ST{log2(X::$BT)::$BT}()
            :($x)
        end

        @generated function log10(::$ST{X}) where X
            x = $ST{log10(X::$BT)::$BT}()
            :($x)
        end

        @generated function log1p(::$ST{X}) where X
            x = $ST{log1p(X::$BT)::$BT}()
            :($x)
        end

        exponent(::$ST{X}) where X = SInt{exponent(X::$BT)}()

        function sincos(::$ST{V}) where V
            x, y = sincos(V::$BT)
            return $ST{x}(), $ST{y}()
        end

        maxintfloat(::Type{<:$ST}) = $ST{maxintfloat($BT)}()
        maxintfloat(::$ST{X}) where X = $ST{maxintfloat(X::$BT)}()

        (/)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = $ST{(/)(V1::$BT, V2::$BT)}()

        Base.fma(::$ST{V1}, ::$ST{V2}, ::$ST{V3}) where {V1,V2,V3} = $ST{fma(V1::$BT, V2::$BT, V3::$BT)}()

        prevfloat(x::$ST{V}) where V = $ST{prevfloat(V::$BT)}()
        prevfloat(x::$ST{V}, n::Integer) where V = $ST{prevfloat(V::$BT, n)}()

        floatmax(x::$ST) = $ST{floatmax($BT)}()
        floatmax(::Type{$ST}) = $ST{floatmax($BT)}()

        floatmin(x::$ST) = $ST{floatmin($BT)}()
        floatmin(::Type{$ST}) = $ST{floatmin($BT)}()

        function Base.mul12(x::$ST, y::$ST)
            Base.@_inline_meta
            h = x * y
            ifelse(iszero(h) | !isfinite(h), (h, h), Base.canonicalize2(h, fma(x, y, -h)))
        end

        function Base.div12(x::$ST, y::$ST)
            xs, xe = frexp(x)
            ys, ye = frexp(y)
            r = xs / ys
            rh, rl = canonicalize2(r, -fma(r, ys, -xs)/ys)
            ifelse(iszero(r) | !isfinite(r), (r, r), (ldexp(rh, xe-ye), ldexp(rl, xe-ye)))
        end
    end
end

struct SFloat64{V} <: AbstractFloat
    function SFloat64{V}() where V
        !(V isa Float64) && throw(ArgumentError("SFloat64 only supports static Float64 storage, got $(typeof(V))"))
        new{V}()
    end
end

struct SFloat32{V} <: AbstractFloat
    function SFloat32{V}() where V
        !(V isa Float32) && throw(ArgumentError("SFloat32 only supports static Float32 storage, got $(typeof(V))"))
        new{V}()
    end
end

struct SFloat16{V} <: AbstractFloat
    function SFloat16{V}() where V
        !(V isa Float16) && throw(ArgumentError("SFloat16 only supports static Float16 storage, got $(typeof(V))"))
        new{V}()
    end
end

const SFloat{V} = Union{SFloat16{V},SFloat32{V},SFloat64{V}}

SFloat(x::Float16) = SFloat16{x}()
SFloat(x::Float32) = SFloat32{x}()
SFloat(x::Float64) = SFloat64{x}()

function SFloat(val::Val{V}) where V
    if V isa Float16
        SFloat16{V}()
    elseif V isa Float32
        SFloat32{V}()
    else
        SFloat64{V}()
    end
end

Base.show(io::IO, v::SFloat) = print(io, "SFloat($(values(v)))")

const AbstractFloat64 = Union{SFloat64,Float64}
const AbstractFloat32 = Union{SFloat32,Float32}
const AbstractFloat16 = Union{SFloat16,Float16}

SF2BF = Dict(SFloat64 => Float64,
             SFloat32 => Float32,
             SFloat16 => Float16)

BF2SF = Dict(Float64 => SFloat64,
             Float32 => SFloat32,
             Float16 => SFloat16)

for (S,B) in SF2BF
    defbasics(S,B)
    defbool(S,B)
    deffloat(S,B)
end

SFloat64One, SFloat64Zero = defmath(SFloat64, Float64)
SFloat32One, SFloat32Zero = defmath(SFloat32, Float32)
SFloat16One, SFloat16Zero = defmath(SFloat16, Float16)

function Base.decompose(::SFloat{X}) where X
    s, e, d = Base.decompose(X)
    return SInt{s}(), SInt{e}(), SInt{d}()
end

for (S,B) in SI2BI
    @eval begin
        Base.trunc(::Type{$B}, x::SFloat{X}) where X = $S{trunc($B, X)}()
    end
end

Base.trunc(::Type{Integer}, x::SFloat{X}) where X = SInt{trunc(Integer, X)}()


text:
function mcmc(setup_mcmc::MCMCSetup, schedule::MCMCSchedule, setups::Vector{T}) where T <: ModelSetup
    workspaces = [create_workspace(setup) for setup in setups]
    wss, lls = [ws[1] for ws in workspaces], [ws[2] for ws in workspaces]
    θ = workspaces[1][3]
    ws_mcmc = create_workspace(setup_mcmc, schedule, θ)
    # adpt = adaptation_object(setup, ws) # adaptive auxiliary law currently not supported

    aux = nothing
    for step in schedule
        step.save && [save_imputed!(ws) for ws in wss]
        for i in step.idx
            if step.param_updt || typeof(ws_mcmc.updates[i]) <: Imputation
                wss = [next(ws, ws_mcmc.updates[i]) for ws in wss]
                lls, acc, θ = update!(ws_mcmc.updates[i], wss, θ, lls, step, aux)
                aux = aux_params(ws_mcmc.updates[i], aux)
                #TODO: make it nicer, currently we need a little dance here
                dance = typeof(ws_mcmc.updates[i]) <: Imputation
                dance && begin
                    update!(ws_mcmc, acc[1], θ, i)
                    for k in 2:length(setups)
                        register_accpt!(ws_mcmc.updates[i], acc[k])
                    end
                end
                !dance && update!(ws_mcmc, acc, θ, i)
                step.verbose && print("\n")
            end
        end
        step.verbose && print("-----------------------------------------------",
                              "------\n")
        step.readjust && readjust!(ws_mcmc, step.iter)
        step.fuse && fuse!(ws_mcmc, schedule)
        # ll = adaptation!(ws, adpt, step.iter, ll) adaptation currently not supported
    end
    [display_summary(ws, ws_mcmc) for ws in wss]
    wss, ws_mcmc
end

function update!(updt::Imputation{NoBlocking}, wss::Vector{<:Workspace{OS}}, θ,
                 lls, step, ::Any, headstart=false) where OS
    ρ, K = updt.ρs[1], length(wss)
    accpts = repeat([false], K)

    for k in 1:K
        ws = wss[k]
        # sample proposal starting point
        zᵒ, yᵒ = proposal_start_pt(ws, ws.P[1], ρ)

        # sample proposal path
        m = length(ws.WWᵒ)
        success = sample_segments!(1:m, ws, yᵒ, ρ, Val{headstart}())
        llᵒ = success ? (logpdf(ws.x0_prior, yᵒ) +
                         path_log_likhd(OS(), ws.XXᵒ, ws.P, 1:m, ws.fpt) +
                         lobslikelihood(ws.P[1], yᵒ)) : -Inf

        print_info(step, value(lls[k]), value(llᵒ), "impute")

        if accept_sample(llᵒ-lls[k], step.verbose)
            swap!(ws.XX, ws.XXᵒ, ws.WW, ws.WWᵒ, 1:m)
            set!(ws.z, zᵒ)
            lls[k] = llᵒ
            accpts[k] = true
        end
    end
    lls, accpts, θ
end

function update!(pu::ParamUpdate{MetropolisHastingsUpdt},
                 wss::Vector{<:Workspace{OS}}, θ, lls, step,
                 blocking::NoBlocking) where OS
    K = length(wss)

    #NOTE for now, let's sample parameter corresponding to a single path, this
    # will need to be changed for Mixed effect models
    θᵒ = rand(pu.t_kernel, θ, pu.updt_coord)               # sample new parameter
    (logpdf(pu.priors, θᵒ) === -Inf) && (return lls, false, θ)

    llᵒs = copy(lls)
    zᵒs = [copy(ws.z.val) for ws in wss]
    for k in 1:K
        ws = wss[k]
        WW, Pᵒ, P, XXᵒ, XX, fpt = ws.WW, ws.Pᵒ, ws.P, ws.XXᵒ, ws.XX, ws.fpt
        m = length(WW)
        update_laws!(Pᵒ, θᵒ)
        pu.aux.recompute_ODEs && solve_back_rec!(blocking, pu.aux.solver, Pᵒ) # compute (H, Hν, c)

        # find white noise which for a given θᵒ gives a correct starting point
        y = XX[1].yy[1]
        zᵒs[k] = inv_start_pt(y, ws.x0_prior, Pᵒ[1])

        success = find_path_from_wiener!(XXᵒ, y, WW, Pᵒ, 1:m)

        llᵒs[k] = success ? (logpdf(ws.x0_prior, y) +
                             path_log_likhd(OS(), XXᵒ, Pᵒ, 1:m, fpt) +
                             lobslikelihood(Pᵒ[1], y)) : -Inf
    end
    ll, llᵒ = sum(lls), sum(llᵒs)
    print_info(step, value(ll), value(llᵒ))

    llr = ( llᵒ - ll + prior_kernel_contrib(pu, θ, θᵒ))

    # Accept / reject
    if accept_sample(llr, step.verbose)
        for k in 1:K
            ws = wss[k]
            WW, Pᵒ, P, XXᵒ, XX, fpt = ws.WW, ws.Pᵒ, ws.P, ws.XXᵒ, ws.XX, ws.fpt
            m = length(WW)
            swap!(XX, XXᵒ, P, Pᵒ, 1:m)
            set!(ws.z, zᵒs[k])
        end
        return llᵒs, true, θᵒ
    else
        return lls, false, θ
    end
end

function update!(pu::ParamUpdate{ConjugateUpdt}, wss::Vector{<:Workspace{OS}},
                 θ, lls, step, blocking::NoBlocking) where OS
    K = length(wss)

    θᵒ = conjugate_draw(θ, [ws.XX for ws in wss], wss[1].P[1].Target, pu.priors[1], pu.updt_coord)

    total_ll_old = sum(lls)
    for k in 1:K
        ws = wss[k]
        WW, P, XX, fpt = ws.WW, ws.P, ws.XX, ws.fpt
        m = length(WW)

        update_laws!(P, θᵒ)
        pu.aux.recompute_ODEs && solve_back_rec!(blocking, pu.aux.solver, P) # compute (H, Hν, c)

        for i in 1:m    # compute wiener path WW that generates XX
            inv_solve!(EulerMaruyamaBounded(), XX[i], WW[i], P[i])
        end
        # compute white noise that generates starting point
        y = XX[1].yy[1]
        z = inv_start_pt(y, ws.x0_prior, P[1])

        llᵒ = logpdf(ws.x0_prior, y)
        llᵒ += path_log_likhd(OS(), XX, P, 1:m, fpt; skipFPT=true)
        llᵒ += lobslikelihood(P[1], y)
        lls[k] = llᵒ
        set!(ws.z, z)
    end
    print_info(step, sum(total_ll_old), sum(lls))
    return lls, true, θᵒ
end

function conjugate_draw(θ, XX::Vector{<:Vector}, PT, prior, updtIdx)
    μ = mustart(updtIdx)
    𝓦 = μ*μ'
    ϑ = SVector(thetaex(updtIdx, θ))
    for X in XX
        μ, 𝓦 = _conjugate_draw(ϑ, μ, 𝓦, X, PT, updtIdx)
    end
    Σ = inv(𝓦 + inv(Matrix(prior.Σ)))
    Σ = (Σ + Σ')/2 # eliminates numerical inconsistencies
    μ_post = Σ * (μ + Vector(prior.Σ\prior.μ))
    ϑ = rand(Gaussian(μ_post, Σ))
    move_to_proper_place(ϑ, θ, updtIdx)     # align so that dimensions agree
end