text:
export logP!

nplus1_impl(N::Integer) = :(Val{$(N+1)})
@generated nplus1(::Type{Val{N}}) where {N} = nplus1_impl(N)

@inline function logP!(grad::AbstractVector{T}, tmp::Vector{T}, θt::AbstractVector{T}, σ::T, prim::dcdp_primitives{FF,T}, isev::ItpSharedEV, uv::NTuple{2,T}, z::NTuple{NZ,Real}, d_obs::Integer, s_idx::Integer, itypidx::NTuple{NI,Real}, dograd::Bool=true) where {FF,T<:Real,NZ,NI}

  if dograd
    length(grad) == length(θt)+1 || throw(DimensionMismatch())
  end

  # unpack information about current state
  itype = getitype.(isev.itypes, itypidx)
  wp = prim.wp
  drng = actionspace(wp, s_idx)

  # information
  # NOTE: truncation of ψ1 can lead to errors in gradient!!!!!!
  ρ = _ρ(σ)
  ψ = s_idx <= end_ex0(wp) ? _ψ1clamp(uv..., ρ, prim) : _ψ2(uv...)

  # we'll reuse this
  ubV = view(tmp, drng.+1)

  @inbounds for d in drng
    sp_idx = sprime(wp, s_idx, d)
    ubV[d+1] = flow(FF, wp, s_idx, θt, σ, z, ψ, d, itype...) + prim.β * isev.EV[z..., ψ, sp_idx, itypidx...]
  end

  dograd ||  return ubV[d_obs+1] - logsumexp(ubV)

  # do this in two steps so that we don't accientally overwrite ubV[d_obs] with Pr(d_obs)
  logp = ubV[d_obs+1]
  logp -= logsumexp_and_softmax!(ubV)

  @inbounds for d in drng
    wt = d==d_obs ? one(T) - ubV[d+1] : -ubV[d+1]
    sp_idx = sprime(wp, s_idx, d)

    @inbounds for k in eachindex(θt) # NOTE: assumes 1-based linear indexing!!
      grad[k] += wt * (flowdθ(FF, wp, s_idx, θt, σ, z, ψ, k, d, itype...) + prim.β * isev.dEV[z..., ψ, k, sp_idx, itypidx...] )
    end

    if s_idx <= end_ex0(wp)
      dpsi = flowdψ(FF, wp, s_idx, θt, σ, z, ψ, d, itype...) + prim.β * gradient_d(nplus1(Val{NZ}), isev.EV, z..., ψ, sp_idx, itypidx...)::T
      dsig = flowdσ(FF, wp, s_idx, θt, σ, z, ψ, d, itype...) + prim.β * isev.dEVσ[z..., ψ, sp_idx, itypidx...]
      grad[end] += wt * (dpsi*_dψ1dθρ(uv..., ρ, σ) + dsig)
    end
  end

  return logp
end


text:
using Documenter
using DocumenterMarkdown
dir = mktempdir()
cp("Project.toml", joinpath(dir, "Project.toml"))
mkdir(joinpath(dir, "src/"))
cp("src/index.md", joinpath(dir, "src/README.md"))
makedocs(format = Markdown(), source = joinpath(dir, "src/"), build = joinpath(dir, "build/"), doctest=false)

s = read(joinpath(dir, "build/README.md"), String)
warn = "\n<!-- README.md is generated from docs/src/index.md. Please edit that file and rebuild with `cd docs/ && julia make_readme.jl`-->\n"
write("../README.md", warn*s)

ext: jl


text:
using Test
using SetProg

@testset "Variables" begin
    @testset "Ellipsoid" begin
        @testset "Dimension" begin
            model = Model()
            @variable(model, S, Ellipsoid())
            err = ErrorException("Missing dimension information, use Ellipsoid(dimension=...) or PolySet(dimension=...)")
            @test_throws err JuMP.optimize!(model)
        end
        for d in 1:3
            @test Ellipsoid(point=SetProg.InteriorPoint(ones(d))).dimension == d
        end
        @testset "Missing Point" begin
            model = Model()
            @variable(model, E, Ellipsoid(dimension=2))
            @objective(model, Max, nth_root(volume(E))) # Force dual space
            err = ArgumentError("Specify a point for nonsymmetric ellipsoid, e.g. `Ellipsoid(point=InteriorPoint([1.0, 0.0]))")
            @test_throws err JuMP.optimize!(model)
        end
    end
    @testset "PolySet" begin
        err = ErrorException("Degree of PolySet not specified, use PolySet(degree=..., ...)")
        @test_throws err PolySet(dimension=1)
        @test_throws ArgumentError("Degree of PolySet not even") PolySet(degree=1)
        @testset "Dimension" begin
            model = Model()
            @variable(model, S, PolySet(degree=2))
            err = ErrorException("Missing dimension information, use Ellipsoid(dimension=...) or PolySet(dimension=...)")
            @test_throws err JuMP.optimize!(model)
        end
        #@testset "Convex" begin
        #    model = Model()
        #    @variable(model, E, PolySet(degree=2, dimension=2))
        #    @objective(model, Max, nth_root(volume(E)))
        #    err = ErrorException("Cannot optimize volume of non-convex polynomial sublevel set. Use PolySet(convex=true, ...)")
        #    @test_throws err begin
        #        JuMP.optimize!(model)
        #    end
        #end
    end
end


text:
module UnstaggeredMesh1D2P
export Mesh

struct Mesh
  x
  p
  q
end

end


text:
# Practical Machine learning
# Association rule based learning - Apriori example
# Chapter 7

# Find k-freq-itemset in given transactions of items queried together
using StatsBase

# Find frequent itemsets from transactions
# @T: transaction list
# @minsupp: minimum support
function find_freq_itemset(T, minsupp)
    N = length(T)
    # Find itemset I from transaction list T
    I = Array(Int64,0)
    for t in T
        for i in t
            push!(I,i)
        end
    end
    I = Set(I)

    # Find freq-itemset when k = 1: Fₖ = {i | i ∈ I^σ({i}) ≥ N × minsupp}
    k = 1
    F = []
    push!(F,map(x->[x],filter(i->σ(i,T) >= N * minsupp, I))) # F₁
    while true
        Cₖ = gen_candidate(F[end]) # Generate candidate set Cₖ from Fₖ₋₁
        Fₖ = filter(c->σ(c,T) >= Nbumanzu * minsupp, Cₖ)
        if !isempty(Fₖ)
            push!(F,Fₖ) # Eliminate infrequent candidates, then set to Fₖ
        else break
        end
    end
    F
end

# Generate freq-itemset from a list of itemsets
# @x: list of itemsets
function gen_candidate(x)
    n = length(x)
    Cₖ = Array(Array{Int64,1},0)
    for a = 1:n, b = 1:n
        if a >= b;continue
        end
        is_candidate = true
        sort!(x[a]); sort!(x[b])
        for i in 1:length(x[1])-1
            if x[a][i] == x[b][i]; continue
            else is_candidate = false; break
            end
        end
        if is_candidate
            push!(Cₖ, sort!([ x[a][1:end-1], x[a][end], x[b][end] ]))
        end
    end
    Cₖ
end

# Generate rules from frequent itemsets
# @x: list of frequent itemsets
# @T: Transaction list
function gen_rules(x, T)
    if length(x) <= 1; return [] # F contains 1-itemsets only, hence no rules generated.
    end
    x = reduce(append!,x[2:end])
    R = Array(Rule,0)
    for f in x # f as each freq-f-itemset fₖ
        ap_genrules!(R,f,map(i->Array([i]),f),T) # H₁ itemset is same as f
    end
    R
end

function ap_genrules!(R, f, H, T)
    k, m = length(f), length(H[1])
    if k > m + 1
        H = gen_candidate(H)
        H_plus_1 = []
        for h in H
            p = setdiff(f,h)
            if conf(p, h, T) >= minconf
                push!(R, Rule(p,h))
                push!(H_plus_1, h)
            end
        end
        ap_genrules(R, f, H_plus_1, T)
    end
end


text:
struct DynamicsVals{T,N,A}
    fVal::Vector{SVector{N,T}}
    xMid::Vector{SVector{N,T}}
    ∇f::Vector{A}
end

function DynamicsVals(dyn_con::DynamicsConstraint)
	DynamicsVals(dyn_con.fVal, dyn_con.xMid, dyn_con.∇f)
end

struct ProblemInfo{T,N}
    model::AbstractModel
    obj::Objective
    conSet::ALConstraintSet{T}
    x0::SVector{N,T}
    xf::SVector{N,T}
end

function ProblemInfo(prob::Problem)
    n = size(prob)[1]
    ProblemInfo(prob.model, prob.obj, ALConstraintSet(prob), SVector{n}(prob.x0), SVector{n}(prob.xf))
end

"""
$(TYPEDEF)
Projected Newton Solver
Direct method developed by the REx Lab at Stanford University
Achieves machine-level constraint satisfaction by projecting onto the feasible subspace.
    It can also take a full Newton step by solving the KKT system.
This solver is to be used exlusively for solutions that are close to the optimal solution.
    It is intended to be used as a "solution polishing" method for augmented Lagrangian methods.
"""
struct ProjectedNewtonSolver{T,N,M,NM} <: ConstrainedSolver{T}
    # Problem Info
    prob::ProblemInfo{T,N}
    Z::Traj{N,M,T,KnotPoint{T,N,M,NM}}
    Z̄::Traj{N,M,T,KnotPoint{T,N,M,NM}}

    opts::SolverOptions{T}
    stats::SolverStats{T}
    P::Primals{T,N,M}
    P̄::Primals{T,N,M}

    H::SparseMatrixCSC{T,Int}
    g::Vector{T}
    # E::Vector{CostExpansion{T,N,N,M}}
    E::TO.CostExpansion{N,M,T}

    D::SparseMatrixCSC{T,Int}
    d::Vector{T}
    λ::Vector{T}

    dyn_vals::DynamicsVals{T}
    active_set::Vector{Bool}

    dyn_inds::Vector{SVector{N,Int}}
    con_inds::Vector{<:Vector}
end

function ProjectedNewtonSolver(prob::Problem,
        opts::SolverOptions=SolverOptions(), stats::SolverStats=SolverStats())
    Z = prob.Z  # grab trajectory before copy to keep associativity
    prob = copy(prob)  # don't modify original problem

    n,m,N = size(prob)
    NN = n*N + m*(N-1)

    # Add dynamics constraints
    TO.add_dynamics_constraints!(prob, integration(prob), 1)
    conSet = prob.constraints
    NP = sum(num_constraints(conSet))

    # Trajectory
    prob_info = ProblemInfo(prob)
    Z̄ = copy(prob.Z)

    # Create concatenated primal vars
    P = Primals(n,m,N)
    P̄ = Primals(n,m,N)

    # Allocate Cost Hessian & Gradient
    H = spzeros(NN,NN)
    g = zeros(NN)
    # E = [CostExpansion{Float64}(n,m) for k = 1:N]
    E = TO.CostExpansion(n,m,N)

    D = spzeros(NP,NN)
    d = zeros(NP)
    λ = zeros(NP)

    fVal = [@SVector zeros(n) for k = 1:N]
    xMid = [@SVector zeros(n) for k = 1:N-1]
    ∇F = [@SMatrix zeros(n,n+m+1) for k = 1:N]
    dyn_vals = DynamicsVals(fVal, xMid, ∇F)
    active_set = zeros(Bool,NP)

    con_inds = gen_con_inds(conSet)

    # Set constant pieces of the Jacobian
    xinds,uinds = P.xinds, P.uinds

    dyn_inds = SVector{n,Int}[]
    ProjectedNewtonSolver(prob_info, Z, Z̄, opts, stats,
        P, P̄, H, g, E, D, d, λ, dyn_vals, active_set, dyn_inds, con_inds)
end

primals(solver::ProjectedNewtonSolver) = solver.P.Z
primal_partition(solver::ProjectedNewtonSolver) = solver.P.xinds, solver.P.uinds

# AbstractSolver interface
Base.size(solver::ProjectedNewtonSolver{T,n,m}) where {T,n,m} = n,m,length(solver.Z)
TO.get_model(solver::ProjectedNewtonSolver) = solver.prob.model
TO.get_constraints(solver::ProjectedNewtonSolver) = solver.prob.conSet
TO.get_trajectory(solver::ProjectedNewtonSolver) = solver.Z
TO.get_objective(solver::ProjectedNewtonSolver) = solver.prob.obj
iterations(solver::ProjectedNewtonSolver) = solver.stats.iterations_pn
get_active_set(solver::ProjectedNewtonSolver) = solver.active_set
solvername(::Type{<:ProjectedNewtonSolver}) = :ProjectedNewton


text:
function load_experimental_data_set(path_to_data_file::String)::Array{Float64,2}

    try

        # load the experimental data -
        df = readdlm(path_to_data_file)

        # return -
        return df
    catch error
        throw(error)
    end
end

function load_configuration_dictionary(path_to_configuration_file)::Dict{String,Any}

    try

        # load the config TOML file -
        configuration_dictionary = TOML.parsefile(path_to_configuration_file)

        # return -
        return configuration_dictionary
    catch error
        throw(error)
    end
end


text:
using Documenter, GMT

makedocs(
	modules = [GMT],
	#format = :html,
	#format = Documenter.HTML(),
	#format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
	sitename = "GMT",
	#assets = ["assets/custom.css"],
	format = Documenter.HTML(assets = ["assets/custom.css"]),
	pages = Any[
		"Introduction"             => "usage.md",
		"Some examples"            => "examples.md",
		"Draw rectangles examples" => "rectangles.md",
		"Draw frames examples"     => "frames.md",
		"Map projections"          => "proj_examples.md",
		"Gallery"                  => [
			"AGU"                   => "gallery/tables.md",
			"Map projections"       => "gallery/mapprojs.md",
			"Historical collection" => "gallery/historic.md",
		],
		hide("gallery/scripts_agu/arrows_I.md"),
		hide("gallery/scripts_agu/arrows_II.md"),
		hide("gallery/scripts_agu/arrows_III.md"),
		hide("gallery/scripts_agu/bars_3D.md"),
		hide("gallery/scripts_agu/bars3_peaks.md"),
		hide("gallery/scripts_agu/bezier.md"),
		hide("gallery/scripts_agu/colored_bars.md"),
		hide("gallery/scripts_agu/compass.md"),
		hide("gallery/scripts_agu/decorated_I.md"),
		hide("gallery/scripts_agu/decorated_II.md"),
		hide("gallery/scripts_agu/error_bars.md"),
		hide("gallery/scripts_agu/flower.md"),
		hide("gallery/scripts_agu/histo_step.md"),
		hide("gallery/scripts_agu/matangles.md"),
		hide("gallery/scripts_agu/hello_I.md"),
		hide("gallery/scripts_agu/spiders.md"),
		hide("gallery/scripts_agu/snake.md"),
		hide("gallery/scripts_agu/solar.md"),
		hide("gallery/scripts_agu/scatter_cart.md"),
		hide("gallery/scripts_agu/scatter_polar.md"),
		hide("gallery/historic/ex01.md"),
		hide("gallery/historic/ex02.md"),
		hide("gallery/historic/ex03.md"),
		hide("gallery/historic/ex04.md"),
		hide("gallery/historic/ex05.md"),
		hide("gallery/historic/ex06.md"),
		hide("gallery/historic/ex07.md"),
		hide("gallery/historic/ex08.md"),
		hide("gallery/historic/ex09.md"),
		hide("gallery/historic/ex10.md"),
		hide("gallery/historic/ex11.md"),
		hide("gallery/historic/ex12.md"),
		hide("gallery/historic/ex13.md"),
		hide("gallery/historic/ex14.md"),
		hide("gallery/historic/ex15.md"),
		hide("gallery/historic/ex16.md"),
		hide("gallery/historic/ex17.md"),
		hide("gallery/historic/ex18.md"),
		hide("gallery/historic/ex19.md"),
		hide("gallery/historic/ex20.md"),
		hide("gallery/historic/ex21.md"),
		hide("gallery/historic/ex22.md"),
		hide("gallery/historic/ex23.md"),
		hide("gallery/historic/ex24.md"),
		hide("gallery/historic/ex25.md"),
		hide("gallery/historic/ex26.md"),
		hide("gallery/historic/ex27.md"),
		hide("gallery/historic/ex28.md"),
		hide("gallery/historic/ex29.md"),
		hide("gallery/historic/ex30.md"),
		hide("gallery/historic/ex32.md"),
		hide("gallery/historic/ex33.md"),
		hide("gallery/historic/ex34.md"),
		hide("gallery/historic/ex35.md"),
		hide("gallery/historic/ex36.md"),
		hide("gallery/historic/ex40.md"),
		hide("gallery/historic/ex41.md"),
		hide("gallery/historic/ex42.md"),
		hide("gallery/historic/ex43.md"),
		hide("gallery/historic/ex44.md"),
		hide("gallery/historic/ex45.md"),
		hide("gallery/historic/ex46.md"),
		hide("gallery/historic/ex48.md"),
		hide("pens.md"),
		hide("grdhisteq.md"),
		"Manual" => [
			"monolitic.md",
			"modules.md",
			"Common options"   => "common_opts.md",
			"General features" => [
				"arrows_control.md",
				"color.md",
				"decorated.md",
				"symbols.md",
			],
		],
		"Modules manuals" => [
			"arrows.md",
			"bar.md",
			"bar3.md",
			"coast.md",
			"colorbar.md",
			"contourf.md",
			"grdcontour.md",
			"grdimage.md",
			"grdgradient.md",
			"grdview.md",
			"histogram.md",
			"lines.md",
			"makecpt.md",
			"movie.md",
			"plot.md",
			"scatter.md",
			"scatter3.md",
			"solar.md",
			"subplot.md",
			"text.md",
		],
		"The GMT types"            => "types.md",
		"Index"                    => "index.md",
	],
)

deploydocs(
	repo   = "github.com/GenericMappingTools/GMT.jl.git",
)


text:
using PlutoUI
using FFTW, Statistics, LinearAlgebra

n = 1000
x = range(0, 1, length=n)
ς = Slider(0.5:0.5:2.0)  # noise level

μ = 3*x.*sin.(2pi*x) # periodic signal in time domain
#μ = 6*sqrt.(abs.(x .- 0.5)) # this one is difficult to estimate

# Model: Signal distorted by white noise
y = μ + ς*randn(n)

# Prior power decay of prior samples in frequency domain
α = Slider(0.5:0.5:2.0)
# α = 3.0 for type integrated Brownian motion
# α = 2.0 for type Brownian motion prior
# α = 1.0 for a very rough prior (pink noise)

H = (α - 1)/2 # equivalent Hurst exponent

# Fourier matrix
# ℱ = [exp(2pi*im*i*j/n) for i in fftshift(-n÷2:n÷2-1), j in fftshift(-n÷2:n÷2-1)]

# Define prior on unit interval in frequency domain
# by choosing power decay
function freq_decay(n, α)  #ifftshift(-n/2:n/2-1)
    fr = 1 ./ fftfreq(n)
    fr[1] = 2n # allow for intercept
    fr = abs.(fr) .^(α/2) # scale frequencies to "unit
    fr/norm(fr)*sqrt(n)
end

# Canonical posterior contraction rate
contract(n, H) =  n^(-H/(1 + 2H))

# Compute posterior mean
γ = ς^(-2) .+ (freq_decay(n, α)).^(-2) # posterior precision in frequency domain
μ̂ = real(ifft(γ.\(fft(ς^(-2)*y))))
m = μ̂ + sqrt(n)*real(ifft(sqrt.(γ).\randn(Complex{Float64},n)));
# note: `fft(x)/sqrt(length(x))` is isometric

# Compute L2 estimation error, compare with theoretical size
@show norm(μ̂ - μ)/sqrt(n),  ς*contract(n, H)

# We can even compute the posterior covariance

# Use circulant property of covariance in time domain, could be done with https://github.com/JuliaMatrices/ToeplitzMatrices.jl
σ = 1.96*sqrt.(abs.((ifft(diag(inv(Diagonal(γ)))))))[1]

# Equivalent to
# Γ = ℱ*Diagonal(freq_decay(n, α).^(-2))*ℱ'/n/n # ≈ ifft(Diagonal(n. + freq_decay(n, α).^(-2)))
# σ = 1.96sqrt.(real(diag(inv((I + Γ)))))

#=
# Plot
using Plots
p = scatter(x, y, markersize=0.5, label="obs")
plot!(x, μ̂, color=:blue, ribbon=σ,fillalpha=0.2, label="posterior mean")
plot!(x, μ, color=:green, label="truth")
z = sqrt(n)*real(ifft(randn(Complex{Float64},n).*abs.(freq_decay(n, α))));
plot!(x, z, label="prior sample")
plot!(x, m, label="posterior sample")

p
=#

max_stars_repo_name: mschauer/nonparbayes

max_issues_repo_name: mschauer/nonparbayes

max_forks_repo_name: mschauer/nonparbayes


text:
export gauss_line_fit

function gaussclossure(X::AbstractVector)
    function gaussloss(params)
        a, σ = params[1], params[2]
        x = a * [exp(-(t)^2 / (2 * σ^2)) for t in axes(X, 1)]
        return sum(abs2, x - X)
    end
    return gaussloss
end

function gauss_line_fit(X::AbstractVector; kwargs...)
    f = gaussclossure(X)
    sX = sum(X)
    σ² = sum(i^2 * X[i] for i in axes(X, 1)) / sX
    σ = sqrt(max(σ², oftype(σ², 1e-4)))  # σ >= 0.01 since we only have integer-scale resolution
    params = [sX / (σ * sqrt(2π)), σ]

    #result = optimize(f, params, Newton(), Optim.Options(; kwargs...); autodiff = :forward)
    result = optimize(f, params, BFGS(), Optim.Options(; kwargs...))
    Optim.converged(result) || @warn "Optimization failed to converge"
    return result
    #return Optim.minimizer(result)
end


text:
module Treaps

import Base: show, isempty, minimum, maximum, in

export Treap, TreapNode, show, isempty, add!, remove!, minimum, maximum, left, right, key, root

typealias PriorityT Float32

type TreapNode{K}
	priority::PriorityT
	key::K
	left::TreapNode{K}
	right::TreapNode{K}
	TreapNode(key, priority, left, right) = new(priority, key, left, right)
	TreapNode() = new(PriorityT(Inf))
end
show(io::IO, t::TreapNode) = show(io, "Key: $(t.key), Priority: $(t.priority)")
isempty(t::TreapNode) = t.priority == PriorityT(Inf)

key(t::TreapNode) = t.key
left(t::TreapNode) = t.left
right(t::TreapNode) = t.right

function add!{K}(t::TreapNode{K}, key::K)
	if isempty(t)
		t.key = key
		t.priority = rand(PriorityT)
		t.left = TreapNode{K}()
		t.right = TreapNode{K}()
		return t
	end

	if key < t.key
		t.left = add!(t.left, key)
		t.left.priority < t.priority ? rotate_right!(t) : t
	else
		@assert t.key < key "A treap may not contain duplicate keys: $key, $(t.key)"
		t.right = add!(t.right, key)
		t.right.priority < t.priority ? rotate_left!(t) : t
	end
end

function merge!{K}(left::TreapNode{K}, right::TreapNode{K})
	isempty(left)  && return right
	isempty(right) && return left
	if left.priority < right.priority
		result = left
		result.right = merge!(result.right, right)
	else
		result = right
		result.left = merge!(left, result.left)
	end
	result
end

function remove!{K}(t::TreapNode{K}, key::K)
	isempty(t) && throw(KeyError(key))
	if key == t.key
		merge!(t.left, t.right)
	elseif key < t.key
		t.left = remove!(t.left, key)
		t.left.priority < t.priority ? rotate_right!(t) : t
	else
		t.right = remove!(t.right, key)
		t.right.priority < t.priority ? rotate_left!(t) : t
	end
end

function in{K}(key::K, t::TreapNode{K})
	isempty(t) && return false
	if key == t.key
		true
	elseif key < t.key
		in(key, t.left)
	else
		in(key, t.right)
	end
end

function minimum(t::TreapNode)
	isempty(t) && error("An empty treap has no minimum.")
	while !isempty(t.left) t = t.left end
	t.key
end

function maximum(t::TreapNode)
	isempty(t) && error("An empty treap has no maximum.")
	while !isempty(t.right) t = t.right end
	t.key
end
function rotate_right!(root::TreapNode)
	@assert !isempty(root)
	newroot = root.left
	root.left, newroot.right = newroot.right, root
	newroot
end

function rotate_left!(root::TreapNode)
	@assert !isempty(root)
	newroot = root.right
	root.right, newroot.left = newroot.left, root
	newroot
end

type Treap{K}
	root::TreapNode{K}
	Treap() = new(TreapNode{K}())
end
add!{K}(t::Treap{K}, key::K) = t.root = add!(t.root, key)
remove!{K}(t::Treap{K}, key::K) = t.root = remove!(t.root, key)
show(io::IO, t::Treap) = show(io, t.root)
isempty(t::Treap) = isempty(t.root)
length(t::Treap) = length(t.root)
minimum(t::Treap) = minimum(t.root)
maximum(t::Treap) = maximum(t.root)
left(t::Treap) = left(t.root)
right(t::Treap) = right(t.root)
key(t::Treap) = key(t.root)
root(t::Treap) = t.root
in{K}(key::K, t::Treap{K}) = in(key, t.root)

include("nearest.jl")

end # module


text:
using AnotherParser
using Test

@testset "test StringLiteral" begin
    matched, v, i = recognize(StringLiteral("abc"),
                              "abcd") # ; index = 2)
    @test matched == true
    @test i == 4
    @test v == "abc"
end

include("SemVerBNF.jl")
include("test_note_BNFNode_location.jl")

include("../examples/BNF/test_bnf_grammar.jl")

@testset "AnotherParser.jl" begin
    # Write your tests here.
end


text:
"""
Constant magnetomotive force.

Parameters:
- `V_m`: [A] Magnetic potential difference
"""
function ConstantMagneticPotentialDifference(;name, V_m=1.0)
    port_p = PositiveMagneticPort()
    port_n = NegativeMagneticPort()
    @parameters V_m=V_m
    @variables Phi(t)
    eqs = [
        V_m ~ port_p.V_m - port_n.V_m
        Phi ~ port_p.Phi
        0 ~ port_p.Phi + port_n.Phi
    ]
    ODESystem(eqs, t, [Phi], [V_m], systems=[port_p, port_n], name=name)
end

"""
Source of constant magnetic flux.

Parameters:
- `Phi`: [Wb] Magnetic flux
"""
function ConstantMagneticFlux(;name, Phi=1.0)
    port_p = PositiveMagneticPort()
    port_n = NegativeMagneticPort()
    @parameters Phi=Phi
    @variables V_m(t)
    eqs = [
        V_m ~ port_p.V_m - port_n.V_m
        Phi ~ port_p.Phi
        0 ~ port_p.Phi + port_n.Phi
    ]
    ODESystem(eqs, t, [Phi], [V_m], systems=[port_p, port_n], name=name)
end


text:
function test_simulation()
    prob = LightDark1D()
    sim = HistoryRecorder(max_steps=100, capture_exception=true)

    solver = MCVISolver(MCVISimulator(), nothing, 1, 10, 8, 50, 100, 500, 10, LightDark1DLowerBound(sim.rng), LightDark1DUpperBound(sim.rng))
    println("Solving...")
    policy = solve(solver, prob)
    println("...Solved")
    up = updater(policy)
    hist = simulate(sim, prob, policy)
    println("Simulation ran for $(n_steps(hist)) steps.")
    io = IOBuffer()
    showerror(io, hist.exception)
    println("Was terminated because of: $(String(take!(io)))")
    println("Reward:", discounted_reward(hist))
    return true
end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

x = [1, 3, 2, 5]
@show x
@show length(x)

y = [4, 5, 6, 1]
z = x .+ y # elementwise operation

X = [1  2; 3 4]

X = reshape([1, 2, 3, 4], 2, 2)

X = permutedims(reshape([1, 2, 3, 4], 2, 2))

X = reshape([1,2,3,4], 2, 2) |> permutedims

x = 4
@show x^2
@show sqrt(x)

sqrt.([4, 9, 16])

using Statistics, StatsBase

x = randn(1_000) # 1_000 points iid from a N(0, 1)
μ = mean(x)
σ = std(x)
@show (μ, σ)

X = [1 2; 3 4; 5 6]
@show X[1, 2]
@show X[:, 1]
@show X[1, :]
@show X[[1, 2], [1, 2]]

size(X)

using CSV

using RDatasets

using DataFrames

auto = dataset("ISLR", "Auto")
first(auto, 3)

describe(auto, :mean, :median, :std)

names(auto)

mpg = auto.MPG
mpg = auto[:, 1]
mpg = auto[:, :MPG]
mpg |> mean

@show size(auto)
@show nrow(auto)
@show ncol(auto)

using PyPlot

figure(figsize=(8,6))
plot(mpg)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module TermInterface

using ..Util
using Base.Meta

# Interface for Expr. Implement these methods for your own type to
# Use it instead of Expr in egraphs!
gethead(e::Expr) = isexpr(e, :call) ? e.args[1] : e.head
getargs(e::Expr) = e.args[(isexpr(e, :call) ? 2 : 1):end]
istree(e::Expr) = true
getmetadata(e::Expr) = (iscall=isexpr(e, :call),)
preprocess(e::Expr) = cleanast(e)

# Fallback implementation for all other types
gethead(e) = e
getargs(e) = []
istree(a) = false
getmetadata(e) = nothing
preprocess(e) = e

export gethead
export getargs
export istree
export getmetadata
export preprocess
end


text:
# We output a message after loading each package to work around a
# ten-minute timeout limit on Travis. Travis assumes the tests have hung
# if the interval between printing something to stdio exceeds ten minutes.
using Images
@info "Finished loading Images package."
using TestImages
@info "Finished loading TestImages package."
using StaticArrays
@info "Finished loading StaticArrays package."
using OffsetArrays
@info "Finished loading OffsetArrays package."
using Random
@info "Finished loading Random package."
using CoordinateTransformations
@info "Finished loading CoordinateTransformations package."

function evaluate_error(dims, flow::Array{SVector{2, Float64}, 2},  Δ,  tol)
    error_count = 0
    maximum_error = 0
    for r = 1:first(dims)
        for c = 1:last(dims)
            δ = flow[r,c]
            if abs(δ[1] - Δ[1]) > tol || abs(δ[2] - Δ[2]) > tol
                    error_count += 1
            end
            error = sqrt( (δ[1] - Δ[1])^2 + (δ[2] - Δ[2])^2 )
            if error > maximum_error
                maximum_error = error
            end
        end
    end
    return error_count, maximum_error
end

@testset "Farneback" begin
    @testset "Polynomial Expansion" begin
        @info "Running Polynomial Expansion test."
        img = Gray{Float64}.(testimage("mandrill"))
        @info "Running Polynomial Expansion Convolution Implementation."
        @time A, B, C = polynomial_expansion(ConvolutionImplementation(), img, 6, 2)
        @info "Running Polynomial Expansion Matrix Implementation."
        @time P, Q, R = polynomial_expansion(MatrixImplementation(), img, 6, 2)
        @test sum(abs.(vec(A)-vec(P))) ≈  0 atol=1e-9
        @test sum(abs.(vec(B)-vec(Q))) ≈  0 atol=1e-9
        @test sum(abs.(vec(C)-vec(R))) ≈  0 atol=1e-9

    end

    @testset "Horizontal Motion" begin
        @info "Running Horizontal Motion test."
        maximum_percentage_error = 7.5
        img1 = Gray{Float64}.(testimage("mandrill"))

        tol = 0.3
        Δ = (0.0, 3.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        algorithm = Farneback(7, estimation_window = 19,
                                 σ_estimation_window = 6.0,
                                 expansion_window = 6,
                                 σ_expansion_window = 1.5)

        @time flow = optical_flow(img1, img2, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Horizontal Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4

        # Same as above, except that we pass in an initial displacement field.
        displacement = Array{SVector{2, Float64}, 2}(undef,size(img1))
        for i in eachindex(displacement)
                displacement[i] = SVector{2, Float64}(0.0, 0.0)
        end
        @time flow = optical_flow(img1, img2, displacement, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Horizontal Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4
    end

    @testset "Vertical Motion" begin
        @info "Running Vertical Motion test."
        maximum_percentage_error = 7.5
        img1 = Gray{Float64}.(testimage("mandrill"))

        tol = 0.3
        Δ = (3.0, 0.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        algorithm = Farneback(7, estimation_window = 19,
                                 σ_estimation_window = 6.0,
                                 expansion_window = 6,
                                 σ_expansion_window = 1.5)

        @time flow = optical_flow(img1, img2, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Vertical Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4

        # Same as above, except that we pass in an initial displacement field.
        displacement = Array{SVector{2, Float64}, 2}(undef,size(img1))
        for i in eachindex(displacement)
                displacement[i] = SVector{2, Float64}(0.0, 0.0)
        end
        @time flow = optical_flow(img1, img2, displacement, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Vertical Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4
    end

    @testset "Combined Motion" begin
        @info "Running Combined  Motion test."
        maximum_percentage_error = 13
        img1 = Gray{Float64}.(testimage("mandrill"))

        tol = 0.3
        Δ = (3.0, -1.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        algorithm = Farneback(7, estimation_window = 19,
                                 σ_estimation_window = 6.0,
                                 expansion_window = 6,
                                 σ_expansion_window = 1.5)

        @time flow = optical_flow(img1, img2, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Combined Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4

        # Same as above, except that we pass in an initial displacement field.
        displacement = Array{SVector{2, Float64}, 2}(undef,size(img1))
        for i in eachindex(displacement)
                displacement[i] = SVector{2, Float64}(0.0, 0.0)
        end
        @time flow = optical_flow(img1, img2, displacement, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Combined Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4

    end

end;


text:
"""
    MenuAdventures.Outfits

A sub-module that enables the player to wear clothes.

```jldoctest
julia> using MenuAdventures

julia> using MenuAdventures.Testing

julia> using MenuAdventures.Outfits

julia> import MenuAdventures: ever_possible

julia> @universe struct Universe <: AbstractUniverse
        end;

julia> @noun struct Room <: AbstractRoom
            already_lit::Bool = true
        end;

julia> @noun struct Person <: Noun
        end;

julia> @noun struct Clothes <: Noun
        end;

julia> ever_possible(::Dress, ::Inventory, ::Clothes) = true;

julia> ever_possible(::Take, ::Reachable, ::Clothes) = true;

julia> ever_possible(::Dress, ::Reachable, ::Person) = true;

julia> cd(joinpath(pkgdir(MenuAdventures), "test", "Outfits")) do
            check_choices() do interface
                you = Person(
                    "Brandon",
                    grammatical_person = second_person,
                    indefinite_article = "",
                )
                room = Room("room")
                coat = Clothes("coat")
                universe = Universe(you, interface = interface)
                universe[room, you] = Containing()
                universe[room, Clothes("coat")] = Containing()
                universe
            end
        end
true
```
"""
module Outfits

using MenuAdventures: Action, BE, get_parent, Inventory, Location, Reachable, Relationship, Sentence, subject_to_verb, Verb
import MenuAdventures: argument_domains, print_sentence, string_relationship_to, verb_for

"""
    Wearing()

A is `Wearing` B means B is worn by A.

A [`Relationship`](@ref)
"""
struct Wearing <: Relationship end
export Wearing

verb_for(::Wearing) = Verb("wear")

function string_relationship_to(thing_answer, ::Wearing, parent_thing)
    buffer = IOBuffer()
    print(buffer, thing_answer.text)
    print(buffer, ' ')
    print(buffer, "that ")
    show(buffer, parent_thing)
    print(buffer, ' ')
    print(buffer, subject_to_verb(parent_thing, BE))
    print(buffer, " wearing")
    String(take!(buffer))
end

"""
    Dress()

`Dress` someone [`Reachable`](@ref) in something from your [`Inventory`](@ref).

An [`Action`](@ref).
"""
struct Dress <: Action end

export Dress

function argument_domains(::Dress)
    Reachable(), Inventory()
end

function (::Dress)(universe, parent_thing, thing)
    universe[parent_thing, thing] = Wearing()
    return false
end

function print_sentence(io, ::Dress, parent_thing_answer, thing_answer)
    print(io, "Dress ")
    print(io, parent_thing_answer.text)
    print(io, " in ")
    print(io, thing_answer.text)
end

end


text:
using Libtask
using Test

include("ctask.jl")
include("tarray.jl")
include("tref.jl")

if get(ENV, "BENCHMARK", nothing) != nothing
    include("benchmarks.jl")
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Test, Tensors
# time = [...]
# strain = [...]
# stress = run_simulation(mat, time, strain)
# @test isapprox(stress, stress_expected)

#
# mat = Material(IdealPlastic, tuple())
# mat.properties.youngs_modulus = 200.0e3
# mat.properties.poissons_ratio = 0.3
# mat.properties.yield_stress = 100.0

parameters = IdealPlasticParameterState(youngs_modulus = 200.0e3,
                                        poissons_ratio = 0.3,
                                        yield_stress = 100.0)

dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64},1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
ddrivers = IdealPlasticDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
mat = IdealPlastic(parameters=parameters, ddrivers=ddrivers)
# mat.dtime = 0.25
# mat.dstrain .= 1.0e-3*[-0.3, -0.3, 1.0, 0.0, 0.0, 0.0]*mat.dtime
integrate_material!(mat)
update_material!(mat)
@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [50.0, 0.0, 0.0, 0.0, 0.0, 0.0]))

# mat.time += mat.dtime
# mat.strain .+= mat.dstrain
# mat.stress .+= mat.dstress
mat.ddrivers = ddrivers
integrate_material!(mat)
update_material!(mat)
@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [100.0, 0.0, 0.0, 0.0, 0.0, 0.0]))
@test isapprox(mat.variables.cumeq, 0.0; atol=1.0e-12)

# mat.time += mat.dtime
# mat.strain .+= mat.dstrain
# mat.stress .+= mat.dstress
# mat.dstrain[:] .= 1.0e-3*[-0.5, -0.5, 1.0, 0.0, 0.0, 0.0]*mat.dtime
dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0]; offdiagscale=2.0)
ddrivers = IdealPlasticDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
mat.ddrivers = ddrivers
integrate_material!(mat)
update_material!(mat)
@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [100.0, 0.0, 0.0, 0.0, 0.0, 0.0]); atol=1.0e-12)
@test isapprox(mat.variables.cumeq, 0.25*1.0e-3)

# mat.time += mat.dtime
# mat.strain .+= mat.dstrain
# mat.stress .+= mat.dstress
# mat.dstrain[:] .= -1.0e-3*[-0.3, -0.3, 1.0, 0.0, 0.0, 0.0]*mat.dtime
dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, -1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
ddrivers = IdealPlasticDriverState(time = 0.25, strain = 0.25*dstrain_dtime)
mat.ddrivers = ddrivers
integrate_material!(mat)
update_material!(mat)
@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [50.0, 0.0, 0.0, 0.0, 0.0, 0.0]); atol=1.0e-12)

# mat.time += mat.dtime
# mat.strain .+= mat.dstrain
# mat.stress .+= mat.dstress
# mat.dtime = 1.0
# m1 = 1.0e-3*[-0.3, -0.3, 1.0, 0.0, 0.0, 0.0]
# m2 = 1.0e-3*[-0.5, -0.5, 1.0, 0.0, 0.0, 0.0]
# mat.dstrain[:] .= -m1*0.75 - m2*0.25
dstrain_dtime = (-0.75*fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0)
                -0.25*fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0]; offdiagscale=2.0))
ddrivers = IdealPlasticDriverState(time = 1.0, strain = dstrain_dtime)
mat.ddrivers = ddrivers
integrate_material!(mat)
integrate_material!(mat)
update_material!(mat)
@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [-100.0, 0.0, 0.0, 0.0, 0.0, 0.0]))
# @test isapprox(mat.properties.dplastic_multiplier, 0.25*1.0e-3)


text:
 using LightGraphs
using MetaGraphs
using Dates
using CSV
using MathOptInterface
const MOI = MathOptInterface

#=
utility.jl

this script implements various helper functions for LPBasic module usage

=#

export IsEdgeOnLine

"""
    IsEdgeOnLine(instance::ProblemInstance, edge::LightGraphs.SimpleGraphs.SimpleEdge, lineIndex::Int64; numerical_returns = false)
Check if the a given *edge* is contained in a the line indexed by *lineIndex* in the given LPBasic problem *instance*.

# Returns:
- true/false
- 1/0, if numerical_returns == true

# Implementation:
- Has O(1) time complexity, by using an internal dictionary pairing each edge to its respective lines
"""
function IsEdgeOnLine(instance::ProblemInstance, edge::LightGraphs.SimpleGraphs.SimpleEdge, lineIndex::Int64; numerical_returns = false)
    line_coverage = get_prop(instance.network, edge, :line_coverage)
    if haskey(line_coverage, lineIndex)
        if numerical_returns return 1 else return true end
    else
        if numerical_returns return 0 else return false end
    end

end

function redirect_to_files(dofunc, outfile, errfile)
    open(outfile, "w") do out
        open(errfile, "w") do err
            redirect_stdout(out) do
                redirect_stderr(err) do
                    dofunc()
                end
            end
        end
    end
end

#=calculates great circle distance using the haversine formula https://en.wikipedia.org/wiki/Haversine_formula
this will be used to calculate point to point distances given latitudes and longitudes
supposes the earth to be perfectly round, but the associated error is very small(given our application)
=#
"""
    GreatCircleDistance(latitude1::Number, longitude1::Number, latitude2::Number, longitude2::Number)::Number
Calculate the great circle distance between points (latitude1, longitude1) and (latitude2, longitude2) on the earth's surface.

# Implementation
Uses the Haversine formula(https://en.wikipedia.org/wiki/Haversine_formula) supposing the Earth is perfectly round
"""
function GreatCircleDistance(latitude1::Number, longitude1::Number, latitude2::Number, longitude2::Number)::Number
    ϕ1, ϕ2 = deg2rad(latitude1), deg2rad(latitude2)
    λ1, λ2 = deg2rad(longitude1), deg2rad(longitude2)
    r = 6371e3 #earth radius in meters

    d = 2*r * asin(sqrt(
                sin((ϕ2 - ϕ1) / 2)^2 +
                cos(ϕ1)*cos(ϕ2)*(sin((λ2-λ1)/2) ^ 2)
    ))

    return d
end

"""
    RunSuite(gtfs_path::String, outpath::String = "./")
Run through the optimization pipeline reading data from gtfs_path and outputting to outpath

Can be used as a sort of example script for the module's usage.
"""
function RunSuite(gtfs_path::String, outpath::String = "./")

    instance = nothing
    try
        instance = ReadNatureFeed(gtfs_path)
    catch err
        #stacktrace()
    end

    try
        instance = ReadGTFS(gtfs_path, transportModeCode = 3)
    catch err
        println(err)
        #stacktrace()
    end

    if instance === nothing
        error("Unable to open given gtfs")
    end

    #save initial directory
    initial_dir = pwd()
    #create new directory for results
    path = "$(outpath)_$(Dates.minute(now()))"
    if !isdir(path)
        working_dir = mkdir(path)
    else
        working_dir = path
    end
    cd(working_dir)

    #aux functions
    ms  = (penalty) -> (instance) -> ConstructModel(instance, useSlack = true, slack_penalty = penalty)
    mb = (instance) -> ConstructModel(instance, useBinaryVariant = true)
    mbs = (penalty) -> (instance) -> ConstructModel(instance, useSlack = true, slack_penalty = penalty, useBinaryVariant = true)
    mm = (min_freq) -> (instance) -> ConstructModel(instance, minFrequencyMultiplier = min_freq)
    mmm = (min_freq, max_freq) -> (instance) -> ConstructModel(instance, minFrequencyMultiplier = min_freq, maxFrequencyMultiplier = max_freq)
    mbmm = (min_freq, max_freq) -> (instance) -> ConstructModel(instance, minFrequencyMultiplier = min_freq, maxFrequencyMultiplier = max_freq, useBinaryVariant = true)
    mnmm = (min_freq, max_freq) -> (instance) -> ConstructModel(instance, minFrequencyMultiplier = min_freq, maxFrequencyMultiplier = max_freq, useNonIntegerDecisionVariables = true)
    mmmslack = (min_freq, max_freq, penalty) -> (instance) -> ConstructModel(instance, minFrequencyMultiplier = min_freq, maxFrequencyMultiplier = max_freq, useSlack = true, slack_penalty = penalty, useNonIntegerDecisionVariables = true)
    #modeling_functions = [mm(0.1), mm(0.3), mm(0.5), mm(0.7), mm(0.8), mm(0.9), mm(1.0), mm(1.1), mm(1.2)]
    #names = ["def-01", "def-03", "def-05", "def-07", "def-08", "def-09", "def-10", "def-11", "def-12"]
    #binary = [false, false, false, false, false, false, false, false, false]

    #default
    #modeling_functions = [mmm(0.1, 20.0), mmm(0.2, 20.0), mmm(0.3, 20.0), mmm(0.4, 20.0), mmm(0.5, 20.0), mmm(0.6, 20.0), mmm(0.7, 20.0), mmm(0.8, 20.0), mmm(0.9, 20.0), mmm(1.0, 20.0), mmm(1.1, 20.0), mmm(1.2, 20.0), mmm(1.3, 20.0), mmm(1.4, 20.0), mmm(1.5, 20.0)]
    #names = ["def-01", "def-02", "def-03", "def-04", "def-05", "def-06", "def-07", "def-08", "def-09", "def-10", "def-11", "def-12", "def-13", "def-14", "def-15"]
    #binary = [false, false, false, false, false, false, false, false, false, false, false, false, false, false, false]

    #non-integer
    #modeling_functions = [mnmm(0.1, 20.0), mnmm(0.2, 20.0), mnmm(0.3, 20.0), mnmm(0.4, 20.0), mnmm(0.5, 20.0), mnmm(0.6, 20.0), mnmm(0.7, 20.0), mnmm(0.8, 20.0), mnmm(0.9, 20.0), mnmm(1.0, 20.0), mnmm(1.1, 20.0), mnmm(1.2, 20.0), mnmm(1.3, 20.0), mnmm(1.4, 20.0), mnmm(1.5, 20.0)]
    #names = ["ni-01", "ni-02", "ni-03", "ni-04", "ni-05", "ni-06", "ni-07", "ni-08", "ni-09", "ni-10", "ni-11", "ni-12", "ni-13", "ni-14", "ni-15"]
    #binary = [false, false, false, false, false, false, false, false, false, false, false, false, false, false, false]

    #slack-100
    penalty = 200.0
    modeling_functions = [mmmslack(0.1, 20.0, penalty), mmmslack(0.2, 20.0, penalty), mmmslack(0.3, 20.0,penalty), mmmslack(0.4, 20.0, penalty), mmmslack(0.5, 20.0, penalty), mmmslack(0.6, 20.0, penalty), mmmslack(0.7, 20.0, penalty), mmmslack(0.8, 20.0, penalty), mmmslack(0.9, 20.0, penalty), mmmslack(1.0, 20.0, penalty), mmmslack(1.1, 20.0, penalty), mmmslack(1.2, 20.0, penalty), mmmslack(1.3, 20.0, penalty), mmmslack(1.4, 20.0, penalty), mmmslack(1.5, 20.0, penalty)]
    names = ["s$(penalty)-01", "s$(penalty)-02", "s$(penalty)-03", "s$(penalty)-04", "s$(penalty)-05", "s$(penalty)-06", "s$(penalty)-07", "s$(penalty)-08", "s$(penalty)-09", "s$(penalty)-10", "s$(penalty)-11", "s$(penalty)-12", "s$(penalty)-13", "s$(penalty)-14", "s$(penalty)-15"]
    binary = [false, false, false, false, false, false, false, false, false, false, false, false, false, false, false]

    #plot instance, without solution
    PlotInstance(instance, out_path = "instance_plot.png", should_preview = false)

    #new dataframes:
    df1 = DataFrame()
    df_upper = DataFrame()
    df_lower = DataFrame()
    df_rest = DataFrame()
    df_dual = DataFrame()
    df_stats = nothing

    println("length: $(length(instance.lines))")
    yb = Array{Bool, 1}()
    for i in 1:length(instance.lines)
        push!(yb, true)
    end
    total_possible_cost = CalculateCost(instance,yb)

    i = 1
    #redirects stdout and stderr
    redirect_to_files("stdout.log", "stderr.log") do
        for f in modeling_functions
            model = f(instance)
            name = names[i]

            println("Starting $name model")

                println("Optimization model number $(i): $(name)")

                time = @elapsed begin optimize!(model, with_optimizer(Cbc.Optimizer)) end

				try

					#dump log_string into file
					open("$(name)-out.txt", "w") do f
						println(f, name)
						print(f, model)
						println(f, termination_status(model))
						if termination_status(model) == MOI.OPTIMAL
                            println(f, primal_status(model))
                            println(f, "objective value", objective_value(model))
						else
                            println(f, "primal_status, NaN")
                            println(f, "objective value", "NaN")
						end
					end
				catch e

				end

                #new dataframe for results
                vars = JuMP.all_variables(model)
                qtd_lines = length(instance.lines)
                n_edges = ne(instance.network)

                y = Array{Int64, 1}()
                yb = Array{Bool, 1}()

                l_slack = Array{Int64, 1}()
                u_slack = Array{Int64, 1}()
                if termination_status(model) == MOI.OPTIMAL
                    for i = 1:length(vars)
                        if i <= qtd_lines
                            push!(y, round(Int64, value(vars[i])))
                            push!(yb, value(vars[i]) > 0)
                        elseif i <= qtd_lines + n_edges
                            push!(l_slack, round(Int64, value((vars[i]))))
                        else #if i <= qtd_lines + n_edges + n_edges ==? length(vars)
                            push!(u_slack, round(Int64, value((vars[i]))))

                        end

                    end
                else
                    for i = 1:length(vars)
                        if i <= qtd_lines
                            push!(y, 0)
                            push!(yb, false)
                        elseif i <= qtd_lines + n_edges
                            push!(l_slack, 0)
                        else #if i <= qtd_lines + n_edges + n_edges ==? length(vars)
                            push!(u_slack, 0)
                        end

                    end
                end

                col_name = Symbol(name)
                df1[!, col_name] = y

                if(length(l_slack) > 0)

                    df_lower[!,col_name] = l_slack
                    df_upper[!,col_name] = u_slack
                else
                    df_lower[!,col_name] = [0 for i in 1:n_edges]
                    df_upper[!,col_name] = [0 for i in 1:n_edges]
                end

                if binary[i]
                    stats = BasicStats(instance, yb)
                else
                    stats = BasicStats(instance, y)
                end

				obj_value = nothing
                if termination_status(model) == MOI.OPTIMAL
					obj_value = JuMP.objective_value(model)
				else
					obj_value = Inf
				end

                stats[:PenaltyCost] = obj_value - stats[:LineCost]
                stats[:ObjectiveValue] = obj_value
                stats[:ElapsedTime] = time
                stats[:PercOfMaxCost] = obj_value / total_possible_cost
                #total possible cost, for for flexible optimization, is a good upper bound, but not really the maximum cost anymore
                stats[:MaxCost] = total_possible_cost

                if(df_stats === nothing)
                    df_stats = DataFrame()
                    keys = [key for (key, _) in stats]
                    df_stats[!,:model] = Any[]
                    for (key,_) in stats
                        df_stats[!, key] = Any[]
                    end
                end

                values = Array{Any,1}()
                pushfirst!(values, name)
                for (_,value) in stats
                    push!(values, value)
                end

                push!(df_stats, values)

                #a dataframe for this model, so its easier to interpret

                df_results = DataFrame()
                df_restrictions = DataFrame()
                df_results[!,:y] = y
                df_results[!,:yb] = yb
                constraint_type = Array{String, 1}()
                constraint_type2 = Array{String, 1}()
                constraint_slack = Array{Float64,1}()
                shadow_price_arr = Array{Float64,1}()
                dual_arr = Array{Float64,1}()
                for ctr_type in list_of_constraint_types(model)
                    for ctr in all_constraints(model, ctr_type[1], ctr_type[2])
                        push!(constraint_type, string(ctr_type[1]))
                        push!(constraint_type2, string(ctr_type[2]))
                        if has_duals(model)
                            push!(shadow_price_arr, shadow_price(ctr))
                            push!(dual_arr, dual(ctr))
                        else
                            push!(shadow_price_arr, NaN)
                            push!(dual_arr, NaN)
                        end
                    end
                end
                df_restrictions[!,:constraint_type] = constraint_type
                df_restrictions[!,:constraint_type2] = constraint_type2
                df_restrictions[!,:dual] = dual_arr
                df_restrictions[!,:shadow_price] = shadow_price_arr
                #shadow_price

                CSV.write("$(name)-results.csv", df_results)
                CSV.write("$(name)-restriction.csv", df_restrictions)

                PlotSolution(instance, yb, out_path = name * "_solution.png", should_preview = false)

            i += 1

        end
    end

    println("Writing data to csv")
    #write dataframes
    CSV.write("solutions.csv", df1)
    CSV.write("slack_upper.csv", df_upper)
    CSV.write("slack_lower.csv", df_lower)
    CSV.write("stats.csv", df_stats)

    #go back to initial dir
    cd(initial_dir)

    println("Finished")

end


text:
export findallmax,
       discount_rewards,
       discount_rewards!,
       CachedSampleAvg,
       SampleAvg,
       CachedSum,
       huber_loss,
       drop_last_dim

"""
    huber_loss(labels, predictions;δ = 1.0)

See [huber_loss](https://en.m.wikipedia.org/wiki/Huber_loss)
and the [implementation](https://github.com/tensorflow/tensorflow/blob/r1.12/tensorflow/python/ops/losses/losses_impl.py#L394-L469) in TensorFlow.

!!! warning
    The return is not reduced!
"""
function huber_loss(labels, predictions;δ = 1.0f0)
    abs_error = abs.(predictions .- labels)
    quadratic = min.(abs_error, 1.0f0)  # quadratic = min.(abs_error, δ)
    linear = abs_error .- quadratic
    0.5f0 .* quadratic .* quadratic .+ linear  # 0.5f0 .* quadratic .* quadratic .+ δ .* linear
end

"""
    findallmax(A::AbstractArray)

Like `findmax`, but all the indices of the maximum value are returned.

!!! warning
    All elements of value `NaN` in `A` will be ignored, unless all elements are `NaN`.
    In that case, the returned maximum value will be `NaN` and the returned indices will be `collect(1:length(A))`

#Examples
```julia-repl
julia> findallmax([-Inf, -Inf, -Inf])
(-Inf, [1, 2, 3])

julia> findallmax([Inf, Inf, Inf])
(Inf, [1, 2, 3])

julia> findallmax([Inf, 0, Inf])
(Inf, [1, 3])

julia> findallmax([0,1,2,1,2,1,0])
(2, [3, 5])
```
"""
function findallmax(A)
    maxval = typemin(eltype(A))
    idxs = Int[]
    for (i, x) in enumerate(A)
        if !isnan(x)
            if x > maxval
                maxval = x
                empty!(idxs)
                push!(idxs, i)
            elseif x == maxval
                push!(idxs, i)
            end
        end
    end
    if length(idxs) == 0
        NaN, collect(1:length(A))
    else
        maxval, idxs
    end
end

function discount_rewards!(new_rewards, rewards, γ)
    new_rewards[end] = rewards[end]
    for i = (length(rewards) - 1):-1:1
        new_rewards[i] = rewards[i] + new_rewards[i+1] * γ
    end
    new_rewards
end

discount_rewards(rewards, γ) = discount_rewards!(similar(rewards), rewards, γ)

discount_rewards_reduced(rewards, γ) = foldr((r, g) -> r + γ * g, rewards)

mutable struct SampleAvg
    t::Int
    avg::Float64
    SampleAvg() = new(0, 0)
end

function (s::SampleAvg)(x)
    s.t += 1
    s.avg += (x - s.avg) / s.t
    s.avg
end

struct CachedSampleAvg
    cache::Dict{Any,SampleAvg}
    CachedSampleAvg() = new(Dict())
end

function (c::CachedSampleAvg)(k, x)
    if !haskey(c.cache, k)
        c.cache[k] = SampleAvg()
    end
    c.cache[k](x)
end

struct CachedSum
    cache::Dict{Any,Float64}
    CachedSum() = new(Dict{Any,Float64}())
end

function (c::CachedSum)(k, x)
    c.cache[k] = get!(c.cache, k, 0.0) + x
end

drop_last_dim(x) = dropdims(x; dims=ndims(x))


text:
"""
    draw_gif_result(time, x, y, ψ, shape, file_path, [, interval, fps]) -> gif

Draw the gif animation from simulation result.
"""
function draw_gif_result(time, x, y, ψ, shape, file_path; interval = 10, fps = 100)
    anim = @animate for i = 1:interval:length(time)
        plot(
            x,
            y,
            label = "",
            xlabel = "x [m]",
            ylabel = "y [m]",
            linestyle = :dot,
            aspect_ratio = :equal,
        )
        ship = square(x[i], y[i], shape, ψ[i])
        plot!(Shape(ship[1], ship[2]), label = "")
        scatter!(
            [x[i]],
            [y[i]],
            seriestype = :scatter,
            title = "time = $(time[i])",
            label = "",
        )
    end
    gif(anim, file_path, fps = fps)
end

function rotate_pos(pos, angle)
    rotate_matrix = [cos(angle) -sin(angle); sin(angle) cos(angle)]
    rotate_matrix * pos
end

function square(center_x, center_y, shape, angle)
    square_xy = [
        rotate_pos([shape[1], shape[2]] / 2, angle) + [center_x, center_y],
        rotate_pos([-shape[1], shape[2]] / 2, angle) + [center_x, center_y],
        rotate_pos([-shape[1], -shape[2]] / 2, angle) + [center_x, center_y],
        rotate_pos([shape[1], -shape[2]] / 2, angle) + [center_x, center_y],
    ]
    xy = hcat(square_xy...)
    [xy[1, :], xy[2, :]]
end

max_stars_repo_name: mitsuyukiLab/JuliaShipMMG

max_issues_repo_name: mitsuyukiLab/JuliaShipMMG


text:
"""
    find_best(s1::AbstractString, iter, dist::PreMetric)

`find_best` returns the element of the iterator `iter` that has the highest similarity score with `s1` according to the distance `dist`.
The function is optimized for `Levenshtein` and `DamerauLevenshtein` distances (potentially modified by `Partial`, `TokenSort`, `TokenSet`, or `TokenMax`)
"""
function find_best(s1::AbstractString, iter_s2, dist::PreMetric; min_score = 0.0)
    best_s2 = nothing
    for s2 in iter_s2
        score = compare(s1, s2, dist; min_score = min_score)
        if score >= min_score
            score == 1.0 && return s2
            best_s2 = s2
            min_score = score
        end
    end
    return best_s2
end

"""
    find_all(s1::AbstractString, iter, dist::PreMetric; min_score = 0.8)
`find_all` returns the vector with all the elements of `iter` that have a similarity score higher than `min_score` according to the distance `dist`.
The function is optimized for `Levenshtein` and `DamerauLevenshtein` distances (potentially modified by `Partial`, `TokenSort`, `TokenSet`, or `TokenMax`)
"""
function find_all(s::AbstractString, iter, dist::PreMetric; min_score = 0.8)
    filter(x -> compare(s, x, dist; min_score = min_score) >= min_score, iter)
end


text:
import Base: show, ==, push!, append!

## nice-valued ranges for histograms
function histrange{T}(v::AbstractArray{T}, n::Integer, closed::Symbol=:right)
    F = float(T)
    nv = length(v)
    if nv == 0 && n < 0
        throw(ArgumentError("number of bins must be ≥ 0 for an empty array, got $n"))
    elseif nv > 0 && n < 1
        throw(ArgumentError("number of bins must be ≥ 1 for a non-empty array, got $n"))
    elseif nv == 0
        return zero(F):zero(F)
    end

    lo, hi = extrema(v)
    histrange(F(lo), F(hi), n, closed)
end

function histrange{F}(lo::F, hi::F, n::Integer, closed::Symbol=:right)
    if hi == lo
        start = F(hi)
        step = one(F)
        divisor = one(F)
        len = one(F)
    else
        bw = (F(hi) - F(lo)) / n
        lbw = log10(bw)
        if lbw >= 0
            step = exp10(floor(lbw))
            r = bw / step
            if r <= 1.1
                nothing
            elseif r <= 2.2
                step *= 2
            elseif r <= 5.5
                step *= 5
            else
                step *= 10
            end
            divisor = one(F)
            start = step*floor(lo/step)
            len = ceil((hi - start)/step)
        else
            divisor = exp10(-floor(lbw))
            r = bw * divisor
            if r <= 1.1
                nothing
            elseif r <= 2.2
                divisor /= 2
            elseif r <= 5.5
                divisor /= 5
            else
                divisor /= 10
            end
            step = one(F)
            start = floor(lo*divisor)
            len = ceil(hi*divisor - start)
        end
    end
    # fix up endpoints
    if closed == :right #(,]
        while lo <= start/divisor
            start -= step
        end
        while (start + (len-1)*step)/divisor < hi
            len += one(F)
        end
    else
        while lo < start/divisor
            start -= step
        end
        while (start + (len-1)*step)/divisor <= hi
            len += one(F)
        end
    end
    FloatRange(start,step,len,divisor)
end

histrange{N}(vs::NTuple{N,AbstractVector},nbins::NTuple{N,Integer},closed::Symbol) =
    map((v,n) -> histrange(v,n,closed),vs,nbins)
histrange{N}(vs::NTuple{N,AbstractVector},nbins::Integer,closed::Symbol) =
    map(v -> histrange(v,nbins,closed),vs)

## histograms ##
function sturges(n)  # Sturges' formula
    n==0 && return one(n)
    ceil(Integer, log2(n))+1
end

abstract AbstractHistogram{T<:Real,N,E}

# N-dimensional histogram object
type Histogram{T<:Real,N,E} <: AbstractHistogram{T,N,E}
    edges::E
    weights::Array{T,N}
    closed::Symbol
    function Histogram(edges::NTuple{N,AbstractArray},weights::Array{T,N},closed::Symbol)
        closed == :right || closed == :left || error("closed must :left or :right")
        map(x -> length(x)-1,edges) == size(weights) || error("Histogram edge vectors must be 1 longer than corresponding weight dimensions")
        new(edges,weights,closed)
    end
end

Histogram{T,N}(edges::NTuple{N,AbstractVector},weights::AbstractArray{T,N},closed::Symbol=:right) =
    Histogram{T,N,typeof(edges)}(edges,weights,closed)

Histogram{T,N}(edges::NTuple{N,AbstractVector},::Type{T},closed::Symbol=:right) =
    Histogram(edges,zeros(T,map(x -> length(x)-1,edges)...),closed)

Histogram{N}(edges::NTuple{N,AbstractVector},closed::Symbol=:right) =
    Histogram(edges,Int,closed)

function show(io::IO, h::AbstractHistogram)
    println(io, typeof(h))
    println(io,"edges:")
    for e in h.edges
        println(io,"  ",e)
    end
    println(io,"weights: ",h.weights)
    print(io,"closed: ",h.closed)
end

(==)(h1::Histogram,h2::Histogram) = (==)(h1.edges,h2.edges) && (==)(h1.weights,h2.weights) && (==)(h1.closed,h2.closed)

# 1-dimensional
Histogram{T}(edge::AbstractVector, weights::AbstractVector{T}, closed::Symbol=:right) =
    Histogram((edge,), weights, closed)

Histogram{T}(edge::AbstractVector, ::Type{T}, closed::Symbol=:right) =
    Histogram(edge, zeros(T,length(edge)-1), closed)

Histogram(edge::AbstractVector,closed::Symbol=:right) =
    Histogram(edge, Int, closed)

function push!{T,E}(h::Histogram{T,1,E}, x::Real,w::Real)
    i = if h.closed == :right
        searchsortedfirst(h.edges[1], x) - 1
    else
        searchsortedlast(h.edges[1], x)
    end
    if 1 <= i <= length(h.weights)
        @inbounds h.weights[i] += w
    end
    h
end
push!{T,E}(h::AbstractHistogram{T,1,E}, x::Real) = push!(h,x,one(T))

function append!{T}(h::AbstractHistogram{T,1}, v::AbstractVector)
    for x in v
        push!(h,x)
    end
    h
end
function append!{T}(h::AbstractHistogram{T,1}, v::AbstractVector,wv::WeightVec)
    for (x,w) in zip(v,wv.values)
        push!(h,x,w)
    end
    h
end

fit(::Type{Histogram},v::AbstractVector, edg::AbstractVector; closed::Symbol=:right) =
    append!(Histogram(edg,closed), v)
fit(::Type{Histogram},v::AbstractVector; closed::Symbol=:right, nbins=sturges(length(v))) =
    fit(Histogram, v, histrange(v,nbins,closed); closed=closed)

fit{W}(::Type{Histogram},v::AbstractVector, wv::WeightVec{W}, edg::AbstractVector; closed::Symbol=:right) =
    append!(Histogram(edg,W,closed), v, wv)
fit(::Type{Histogram},v::AbstractVector, wv::WeightVec; closed::Symbol=:right, nbins=sturges(length(v))) =
    fit(Histogram, v, wv, histrange(v,nbins,closed); closed=closed)

# N-dimensional
function push!{T,N}(h::Histogram{T,N},xs::NTuple{N,Real},w::Real)
    is = if h.closed == :right
        map((edge, x) -> searchsortedfirst(edge,x) - 1, h.edges, xs)
    else
        map(searchsortedlast, h.edges, xs)
    end
    try
        h.weights[is...] += w
    catch e
        isa(e,BoundsError) || rethrow(e)
    end
    h
end
push!{T,N}(h::AbstractHistogram{T,N},xs::NTuple{N,Real}) = push!(h,xs,one(T))

function append!{T,N}(h::AbstractHistogram{T,N}, vs::NTuple{N,AbstractVector})
    for xs in zip(vs...)
        push!(h,xs)
    end
    h
end
function append!{T,N}(h::AbstractHistogram{T,N}, vs::NTuple{N,AbstractVector},wv::WeightVec)
    for (xs,w) in zip(zip(vs...),wv.values)
        push!(h,xs,w)
    end
    h
end

fit{N}(::Type{Histogram}, vs::NTuple{N,AbstractVector}, edges::NTuple{N,AbstractVector}; closed::Symbol=:right) =
    append!(Histogram(edges,closed), vs)
fit{N}(::Type{Histogram}, vs::NTuple{N,AbstractVector}; closed::Symbol=:right, nbins=sturges(length(vs[1]))) =
    fit(Histogram, vs, histrange(vs,nbins,closed); closed=closed)

fit{N,W}(::Type{Histogram}, vs::NTuple{N,AbstractVector}, wv::WeightVec{W}, edges::NTuple{N,AbstractVector}; closed::Symbol=:right) =
    append!(Histogram(edges,W,closed), vs, wv)
fit{N}(::Type{Histogram},vs::NTuple{N,AbstractVector}, wv::WeightVec; closed::Symbol=:right, nbins=sturges(length(vs[1]))) =
    fit(Histogram, vs, wv, histrange(vs,nbins,closed); closed=closed)


text:
using PowerModels

testcases = [
	Dict(
	"file" => "../data/case9.m",
 	"name" => "case9-1",
 	"attack_budget" => 3,
 	"inactive_indices" => [],
 	"protected_indices" => []
	),
	Dict(
	"file" => "../data/case9.m",
	"name" => "case9-2",
 	"attack_budget" => 3,
 	"inactive_indices" => [1,2,3,4,5,6,7,8,9],
 	"protected_indices" => []
	),
	Dict(
	"file" => "../data/case9.m",
	"name" => "case9-3",
 	"attack_budget" => 3,
 	"inactive_indices" => [],
 	"protected_indices" => [1,2,3,4,5,6,7,8,9]
	),
	Dict(
	"file" => "../data/case9.m",
	"name" => "case9-4",
 	"attack_budget" => 3,
 	"inactive_indices" => [1,4],
 	"protected_indices" => []
	),
	Dict(
	"file" => "../data/case9.m",
	"name" => "case9-5",
 	"attack_budget" => 3,
 	"inactive_indices" => [],
 	"protected_indices" => [2,3]
	),
	Dict(
	"file" => "../data/case9.m",
	"name" => "case9-6",
 	"attack_budget" => 3,
 	"inactive_indices" => [1,4,7],
 	"protected_indices" => [2,3,5,6,8,9]
	),
	Dict(
	"file" => "../data/case9.m",
	"name" => "case9-7",
 	"attack_budget" => 3,
 	"inactive_indices" => [1],
 	"protected_indices" => [2]
	)
]


text:
@testset "app" begin

@testset "exit" begin
	app = MiniApp()
	put!(app.fromGUI, :exit=>[])
	process_thread(app.jg, app.fromGUI, app.toGUI; verbosityLevel=2)
	didExit = false
	while true
		@test isready(app.toGUI)
		msg = take!(app.toGUI)
		msg.first == :exited && (didExit=true; break)
	end
	@test didExit
end

@testset "exit2" begin
	app = MiniApp()
	init(app)
	lastSchedulerTime = Ref{UInt64}(0)
	put!(app.fromGUI, :exit=>[])
	@test !runall(app)
end

@testset "data.tsv" begin
	filepath = joinpath(@__DIR__, "data/data.tsv")

	varIds = [Symbol("V$i") for i=1:40]
	sampleIds = [string('S',lpad(i,2,'0')) for i=1:20]
	groupAnnot = repeat(["A","B"],inner=10)
	timeAnnot = vcat(0.01:0.01:0.1, 0.01:0.01:0.1)

	app = MiniApp()
	init(app)

	# Open file
	setvalue(app, Any["samplefilepath", filepath])
	setvalue(app, Any["lastsampleannot", "Time"])
	put!(app.fromGUI, :loadsample=>[])
	@test runall(app)
	dfSample = app.jg.scheduler.jobs[app.jg.loadSampleID].result
	@test propertynames(dfSample) == vcat([:SampleId, :Group, :Time], varIds)
	@test dfSample.SampleId == sampleIds
	@test dfSample.Group == groupAnnot
	@test dfSample.Time == timeAnnot
	@test size(dfSample)==(20,3+40)

	X = Matrix(Matrix{Float64}(dfSample[:,4:end])')
	normalizemean!(X)

	# PMA (groups)
	setvalue(app, Any["samplesimplexmethod", "SA"])
	setvalue(app, Any["sampleannot", "Group"])
	put!(app.fromGUI, :dimreduction=>[])
	@test runall(app)
	reduced = app.jg.scheduler.jobs[app.jg.dimreductionID].result
	@test reduced.sa == dfSample[!,1:3]
	@test reduced.va[!,1] == varIds
	factorizationcmp(pma(X, groupsimplices(groupAnnot); nsv=10), reduced.F)

	# PMA (time series)
	setvalue(app, Any["samplesimplexmethod", "Time"])
	setvalue(app, Any["sampleannot", "Group"])
	setvalue(app, Any["timeannot", "Time"])
	put!(app.fromGUI, :dimreduction=>[])
	@test runall(app)
	reduced = app.jg.scheduler.jobs[app.jg.dimreductionID].result
	@test reduced.sa == dfSample[!,1:3]
	@test reduced.va[!,1] == varIds
	#factorizationcmp(pma(X, timeseriessimplices(timeAnnot, groupby=groupAnnot); nsv=10), reduced.F)
	G = timeseriessimplices(timeAnnot, groupby=groupAnnot)
	factorizationcmp(pma(X, G; nsv=10), reduced.F)

	# PMA (NN)
	setvalue(app, Any["samplesimplexmethod", "NN"])
	setvalue(app, Any["knearestneighbors", "2"])
	setvalue(app, Any["distnearestneighbors", "0.5"])
	put!(app.fromGUI, :dimreduction=>[])
	@test runall(app)
	reduced = app.jg.scheduler.jobs[app.jg.dimreductionID].result
	@test reduced.sa == dfSample[!,1:3]
	@test reduced.va[!,1] == varIds
	factorizationcmp(pma(X, neighborsimplices(X,k=2,r=0.5,dim=50), nsv=10), reduced.F)

	# PMA (NN withing groups)
	setvalue(app, Any["samplesimplexmethod", "NNSA"])
	setvalue(app, Any["sampleannot", "Group"])
	setvalue(app, Any["knearestneighbors", "2"])
	setvalue(app, Any["distnearestneighbors", "0.5"])
	put!(app.fromGUI, :dimreduction=>[])
	@test runall(app)
	reduced = app.jg.scheduler.jobs[app.jg.dimreductionID].result
	@test reduced.sa == dfSample[!,1:3]
	@test reduced.va[!,1] == varIds
	F = pma(X, neighborsimplices(X,k=2,r=0.5,dim=50,groupby=groupAnnot), nsv=10)
	factorizationcmp(F, reduced.F)

	# Exports
	mktempdir() do temppath
		for (dim,order,mode) in zip((1,3,7,10), ("Abs","Descending","Ascending","Original"), ("Samples","Variables","Samples","Variables"))
			setvalue(app, Any["exportsingledim", "$dim"])
			setvalue(app, Any["exportsinglesort", order])
			setvalue(app, Any["exportmode", mode])
			filepath = joinpath(temppath,"PMA$(dim)_$(mode)_$order.tsv")
			setvalue(app, Any["exportsinglepath", filepath])
			put!(app.fromGUI, :exportsingle=>[])
			@test runall(app)
			result = CSV.read(filepath, DataFrame; delim='\t', ntasks=1)
			colName = Symbol(:PMA,dim)
			resultPMA = result[:,colName]

			if mode=="Samples"
				idCol = :SampleId
				ids = sampleIds
				PMA = reduced.F.V[:,dim]
			else
				idCol = :VariableId
				ids = string.(varIds)
				PMA = reduced.F.U[:,dim]
			end

			perm = collect(1:length(PMA))
			if order=="Abs"
				perm = sortperm(PMA, by=abs, rev=true)
			elseif order=="Descending"
				perm = sortperm(PMA, rev=true)
			elseif order=="Ascending"
				perm = sortperm(PMA)
			end

			@test propertynames(result)==[idCol,colName]
			@test result[!,idCol] == ids[perm]
			@test resultPMA ≈ PMA[perm]
		end

		for (dim,mode) in zip((4,5), ("Samples","Variables"))
			setvalue(app, Any["exportmultipledim", "$dim"])
			setvalue(app, Any["exportmode", mode])
			filepath = joinpath(temppath,"PMA$(dim)_$(mode).csv")
			setvalue(app, Any["exportmultiplepath", filepath])
			put!(app.fromGUI, :exportmultiple=>[])
			@test runall(app)
			result = CSV.read(filepath, DataFrame; delim=',', ntasks=false)

			colNames = Symbol.(:PMA,1:dim)

			if mode=="Samples"
				@test propertynames(result)==vcat(:SampleId, colNames)
				@test result.SampleId == sampleIds
				@test Matrix(result[:,colNames]) ≈ reduced.F.V[:,1:dim]
			else
				@test propertynames(result)==vcat(:VariableId, colNames)
				@test result.VariableId == string.(varIds)
				@test Matrix(result[:,colNames]) ≈ reduced.F.U[:,1:dim]
			end
		end
	end

	# Exit
	put!(app.fromGUI, :exit=>[])
	@test !runall(app)
end

@testset "MissingValueReconstruction" begin
	reconstructed = Float64[0.06 0.26 0.0 0.0 0.28 0.03 0.04 0.12 0.05 0.0 0.02 0.31 0.0 0.85 0.41 0.06 0.04 0.92 0.32 0.12; 0.01 0.13 0.48 0.64 0.23 0.01 0.0 0.18 0.03 0.14 0.26 0.0 0.09 0.27 0.1 0.0 0.48 0.06 0.07 0.56; 0.0 0.22 0.02 0.0 0.02 0.06 0.05 0.59 0.64 0.81 0.98 0.24 0.08 0.24 0.41 0.0 0.65 0.15 0.11 0.0; 0.01 0.01 0.31 0.23 0.02 0.0 0.01 0.15 0.03 0.64 0.21 0.16 0.0 0.02 0.12 0.84 0.06 0.0 0.0 0.2; 0.03 0.55 0.0 0.06 0.01 0.01 0.16 0.01 0.42 0.73 0.01 0.15 0.27 0.02 0.02 0.37 0.16 0.19 0.04 0.0; 0.36 0.97 0.54 0.98 0.03 0.52 0.1 0.06 0.14 0.0 0.46 0.38 0.29 0.93 0.06 0.07 0.05 0.45 0.0 0.01; 0.12 0.97 0.3 0.07 0.11 0.28 0.31 0.07 0.0 0.16 0.0 0.29 0.33 0.01 0.01 0.02 0.17 0.0 0.67 0.45; 0.0 0.08 0.49 0.34 0.07 0.01 0.03 0.1 0.05 0.16 0.05 0.36 0.0 0.38 0.41 0.0 0.09 0.67 0.61 1.0; 0.01 0.3 0.19 0.0 0.0 0.0 0.32 0.72 0.0 0.04 0.7 1.0 0.55 0.12 0.03 0.19 0.52 0.04 0.04 0.8; 0.12 0.0 0.04 0.02 0.19 0.09 0.01 0.66 0.25 0.21 0.01 0.02 0.38 0.62 0.3 0.99 0.38 0.1 0.54 0.07; 0.08 0.0 0.16 0.01 0.69 0.23 0.81 0.29 0.13 0.4 0.04 0.15 0.17 0.64 0.03 0.32 0.78 0.28 0.12 0.36; 0.89 0.02 0.0 0.02 0.26 0.06 0.01 0.01 0.48 0.58 0.25 0.29 0.51 0.31 0.22 0.0 0.27 0.01 0.04 0.05; 0.55 0.04 0.47 0.0 0.24 0.07 0.64 0.01 0.02 0.0 0.42 0.0 0.36 0.22 0.06 0.95 0.66 0.27 0.09 0.05; 0.46 0.19 0.1 0.01 0.13 0.15 0.06 0.84 0.02 0.0 0.09 0.62 0.01 0.0 0.01 0.0 0.01 0.33 0.47 0.58; 0.04 0.28 0.27 0.69 0.0 0.78 0.12 0.01 0.13 0.0 0.0 0.09 0.11 0.53 0.0 0.29 0.39 0.67 0.42 0.29; 0.01 0.09 0.0 0.49 0.18 0.0 0.64 0.67 0.32 0.37 0.08 0.0 0.04 0.17 0.0 0.06 0.01 0.0 0.43 0.22; 0.27 0.02 0.05 0.88 0.02 0.16 0.09 0.0 0.46 0.05 0.18 0.11 0.02 0.37 0.88 0.47 0.92 0.41 0.96 0.0; 0.0 0.99 0.98 0.31 0.05 0.02 0.06 0.03 0.34 0.37 0.17 0.0 0.33 0.78 0.47 0.75 0.34 0.04 0.69 0.0; 0.69 0.34 0.15 0.12 0.11 0.0 0.91 0.35 0.05 0.0 0.26 0.01 0.62 0.82 0.48 0.4 0.08 0.0 0.32 0.0; 0.08 0.01 0.85 0.82 0.0 0.0 0.99 0.0 0.12 0.05 0.04 0.53 0.0 0.89 0.26 0.01 0.17 0.39 0.05 0.0; 0.01 0.01 0.72 0.2647368421052631 0.0 0.6 0.01 0.08 0.18 0.42 0.48 0.02 0.34 0.59 0.02 0.6 0.02 0.76 0.03 0.14; 0.17555555555555552 0.13 0.47 0.01 0.08 0.26 0.25 0.17555555555555552 0.0 0.0 0.01 0.82 0.0 0.0 0.03 0.0 0.0 0.3 0.07 0.73; 0.0 0.3664705882352941 0.39 0.99 0.03 0.5 0.88 0.13 0.8 0.0 0.0 0.0 0.3664705882352941 0.3664705882352941 0.43 0.7 0.72 0.22 0.04 0.4; 0.0 0.5 0.5 0.97 0.02 0.08 0.35 0.238125 0.02 0.26 0.02 0.238125 0.1 0.0 0.0 0.92 0.0 0.238125 0.238125 0.07; 0.09 0.0 0.09733333333333334 0.09733333333333334 0.56 0.0 0.09733333333333334 0.09733333333333334 0.0 0.14 0.01 0.0 0.02 0.08 0.12 0.0 0.09733333333333334 0.09 0.1 0.25; 0.34500000000000003 0.34500000000000003 0.13 0.63 0.34500000000000003 0.03 0.0 0.06 0.56 0.55 0.15 0.34500000000000003 0.87 0.9 0.79 0.04 0.34500000000000003 0.34500000000000003 0.09 0.03; 0.85 0.79 0.32769230769230767 0.32769230769230767 0.35 0.32769230769230767 0.32769230769230767 0.21 0.03 0.01 0.36 0.44 0.52 0.01 0.32769230769230767 0.66 0.32769230769230767 0.03 0.0 0.32769230769230767; 0.86 0.24499999999999997 0.24499999999999997 0.24 0.01 0.24499999999999997 0.84 0.24499999999999997 0.24499999999999997 0.15 0.24499999999999997 0.24499999999999997 0.24499999999999997 0.1 0.13 0.3 0.03 0.04 0.24 0.0; 0.3427272727272727 0.92 0.3427272727272727 0.3427272727272727 0.09 0.3427272727272727 0.9 0.0 0.3427272727272727 0.3427272727272727 0.02 0.3427272727272727 0.63 0.05 0.19 0.09 0.0 0.88 0.3427272727272727 0.3427272727272727; 0.6 0.271 0.271 0.271 0.271 0.271 0.16 0.271 0.271 0.271 0.271 0.271 0.02 0.4 0.04 0.0 0.18 0.34 0.01 0.96; 0.0 0.1688888888888889 0.24 0.0 0.1688888888888889 0.2 0.02 0.1688888888888889 0.33 0.1688888888888889 0.1688888888888889 0.1688888888888889 0.17 0.1688888888888889 0.06 0.1688888888888889 0.1688888888888889 0.1688888888888889 0.1688888888888889 0.5; 0.26375000000000004 0.26375000000000004 0.84 0.26375000000000004 0.26375000000000004 0.42 0.26375000000000004 0.26375000000000004 0.26375000000000004 0.26375000000000004 0.54 0.26375000000000004 0.26375000000000004 0.0 0.13 0.01 0.04 0.26375000000000004 0.26375000000000004 0.13; 0.21142857142857144 0.21142857142857144 0.21142857142857144 0.21142857142857144 0.27 0.03 0.21142857142857144 0.11 0.21142857142857144 0.5 0.21142857142857144 0.21142857142857144 0.21142857142857144 0.21142857142857144 0.21142857142857144 0.21142857142857144 0.21142857142857144 0.15 0.36 0.06; 0.22833333333333336 0.22833333333333336 0.22833333333333336 0.22833333333333336 0.1 0.22833333333333336 0.22833333333333336 0.01 0.02 0.22833333333333336 0.01 0.22833333333333336 0.22833333333333336 0.22833333333333336 0.22833333333333336 0.64 0.22833333333333336 0.59 0.22833333333333336 0.22833333333333336; 0.08800000000000001 0.08800000000000001 0.08800000000000001 0.04 0.08800000000000001 0.08800000000000001 0.08 0.15 0.17 0.08800000000000001 0.08800000000000001 0.08800000000000001 0.08800000000000001 0.08800000000000001 0.08800000000000001 0.0 0.08800000000000001 0.08800000000000001 0.08800000000000001 0.08800000000000001; 0.52 0.62 0.30250000000000005 0.07 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.0 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.30250000000000005 0.30250000000000005; 0.043333333333333335 0.08 0.043333333333333335 0.043333333333333335 0.043333333333333335 0.043333333333333335 0.0 0.043333333333333335 0.043333333333333335 0.043333333333333335 0.043333333333333335 0.043333333333333335 0.043333333333333335 0.05 0.043333333333333335 0.043333333333333335 0.043333333333333335 0.043333333333333335 0.043333333333333335 0.043333333333333335; 0.34 0.34 0.34 0.34 0.34 0.34 0.34 0.68 0.34 0.0 0.34 0.34 0.34 0.34 0.34 0.34 0.34 0.34 0.34 0.34; 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
	logReconstructed = Float64[-0.6666 -0.2688 -0.811 -0.811 -0.2345 -0.737 -0.7131 -0.5353 -0.6897 -0.811 -0.7612 -0.1844 -0.811 0.5059 -0.02915 -0.6666 -0.7131 0.5753 -0.1681 -0.5353; -0.7859 -0.5146 0.07039 0.275 -0.3219 -0.7859 -0.811 -0.415 -0.737 -0.4941 -0.2688 -0.811 -0.5995 -0.2515 -0.5778 -0.811 0.07039 -0.6666 -0.6439 0.1763; -0.811 -0.3401 -0.7612 -0.811 -0.7612 -0.6666 -0.6897 0.2141 0.275 0.4647 0.6323 -0.304 -0.6215 -0.304 -0.02915 -0.811 0.2869 -0.4739 -0.5564 -0.811; -0.7859 -0.7859 -0.1844 -0.3219 -0.7612 -0.811 -0.7859 -0.4739 -0.737 0.275 -0.3585 -0.454 -0.811 -0.7612 -0.5353 0.4957 -0.6666 -0.811 -0.811 -0.3771; -0.737 0.1635 -0.811 -0.6666 -0.7859 -0.7859 -0.454 -0.7859 -0.0145 0.3785 -0.7859 -0.4739 -0.2515 -0.7612 -0.7612 -0.08927 -0.454 -0.3959 -0.7131 -0.811; -0.1047 0.6229 0.1506 0.6323 -0.737 0.1243 -0.5778 -0.6666 -0.4941 -0.811 0.04264 -0.074 -0.2176 0.585 -0.6666 -0.6439 -0.6897 0.02857 -0.811 -0.7859; -0.5353 0.6229 -0.2009 -0.6439 -0.5564 -0.2345 -0.1844 -0.6439 -0.811 -0.454 -0.811 -0.2176 -0.152 -0.7859 -0.7859 -0.7612 -0.4344 -0.811 0.3103 0.02857; -0.811 -0.6215 0.08406 -0.1361 -0.6439 -0.7859 -0.737 -0.5778 -0.6897 -0.454 -0.6897 -0.1047 -0.811 -0.074 -0.02915 -0.811 -0.5995 0.3103 0.2388 0.6508; -0.7859 -0.2009 -0.3959 -0.811 -0.811 -0.811 -0.1681 0.3674 -0.811 -0.7131 0.3448 0.6508 0.1635 -0.5353 -0.737 -0.3959 0.1243 -0.7131 -0.7131 0.4542; -0.5353 -0.811 -0.7131 -0.7612 -0.3959 -0.5995 -0.7859 0.2987 -0.2863 -0.3585 -0.7859 -0.7612 -0.074 0.251 -0.2009 0.6415 -0.074 -0.5778 0.1506 -0.6439; -0.6215 -0.811 -0.454 -0.7859 0.3334 -0.3219 0.4647 -0.2176 -0.5146 -0.04394 -0.7131 -0.4739 -0.4344 0.275 -0.737 -0.1681 0.433 -0.2345 -0.5353 -0.1047; 0.546 -0.7612 -0.811 -0.7612 -0.2688 -0.6666 -0.7859 -0.7859 0.07039 0.2016 -0.2863 -0.2176 0.111 -0.1844 -0.3401 -0.811 -0.2515 -0.7859 -0.7131 -0.6897; 0.1635 -0.7131 0.05658 -0.811 -0.304 -0.6439 0.275 -0.7859 -0.7612 -0.811 -0.0145 -0.811 -0.1047 -0.3401 -0.6666 0.6041 0.2987 -0.2515 -0.5995 -0.6897; 0.04264 -0.3959 -0.5778 -0.7859 -0.5146 -0.4739 -0.6666 0.4957 -0.7612 -0.811 -0.5995 0.251 -0.7859 -0.811 -0.7859 -0.811 -0.7859 -0.152 0.05658 0.2016; -0.7131 -0.2345 -0.2515 0.3334 -0.811 0.433 -0.5353 -0.7859 -0.5146 -0.811 -0.811 -0.5995 -0.5564 0.1375 -0.811 -0.2176 -0.05889 0.3103 -0.0145 -0.2176; -0.7859 -0.5995 -0.811 0.08406 -0.415 -0.811 0.275 0.3103 -0.1681 -0.08927 -0.6215 -0.811 -0.7131 -0.4344 -0.811 -0.6666 -0.7859 -0.811 0.0 -0.3401; -0.2515 -0.7612 -0.6897 0.5361 -0.7612 -0.454 -0.5995 -0.811 0.04264 -0.6897 -0.415 -0.5564 -0.7612 -0.08927 0.5361 0.05658 0.5753 -0.02915 0.6135 -0.811; -0.811 0.6415 0.6323 -0.1844 -0.6897 -0.7612 -0.6666 -0.737 -0.1361 -0.08927 -0.4344 -0.811 -0.152 0.433 0.05658 0.4005 -0.1361 -0.7131 0.3334 -0.811; 0.3334 -0.1361 -0.4739 -0.5353 -0.5564 -0.811 0.5656 -0.1203 -0.6897 -0.811 -0.2688 -0.7859 0.251 0.4751 0.07039 -0.04394 -0.6215 -0.811 -0.1681 -0.811; -0.6215 -0.7859 0.5059 0.4751 -0.811 -0.811 0.6415 -0.811 -0.5353 -0.6897 -0.7131 0.1375 -0.811 0.546 -0.2688 -0.7859 -0.4344 -0.05889 -0.6897 -0.811; -0.7859 -0.7859 0.3674 -0.3344 -0.811 0.2265 -0.7859 -0.6215 -0.415 -0.0145 0.07039 -0.7612 -0.1361 0.2141 -0.7612 0.2265 -0.7612 0.4114 -0.737 -0.4941; -0.4889 -0.5146 0.05658 -0.7859 -0.6215 -0.2688 -0.2863 -0.4889 -0.811 -0.811 -0.7859 0.4751 -0.811 -0.811 -0.737 -0.811 -0.811 -0.2009 -0.6439 0.3785; -0.811 -0.1888 -0.05889 0.6415 -0.737 0.09761 0.5361 -0.5146 0.4542 -0.811 -0.811 -0.811 -0.1888 -0.1888 0.0 0.3448 0.3674 -0.3401 -0.7131 -0.04394; -0.811 0.09761 0.09761 0.6229 -0.7612 -0.6215 -0.1203 -0.3979 -0.7612 -0.2688 -0.7612 -0.3979 -0.5778 -0.811 -0.811 0.5753 -0.811 -0.3979 -0.3979 -0.6439; -0.5995 -0.811 -0.6093 -0.6093 0.1763 -0.811 -0.6093 -0.6093 -0.811 -0.4941 -0.7859 -0.811 -0.7612 -0.6215 -0.5353 -0.811 -0.6093 -0.5995 -0.5778 -0.2863; -0.2232 -0.2232 -0.5146 0.263 -0.2232 -0.737 -0.811 -0.6666 0.1763 0.1635 -0.4739 -0.2232 0.5261 0.5558 0.4436 -0.7131 -0.2232 -0.2232 -0.5995 -0.737; 0.5059 0.4436 -0.2349 -0.2349 -0.1203 -0.2349 -0.2349 -0.3585 -0.737 -0.7859 -0.1047 0.01436 0.1243 -0.7859 -0.2349 0.2987 -0.2349 -0.737 -0.811 -0.2349; 0.516 -0.3675 -0.3675 -0.304 -0.7859 -0.3675 0.4957 -0.3675 -0.3675 -0.4739 -0.3675 -0.3675 -0.3675 -0.5778 -0.5146 -0.2009 -0.737 -0.7131 -0.304 -0.811; -0.2489 0.5753 -0.2489 -0.2489 -0.5995 -0.2489 0.5558 -0.811 -0.2489 -0.2489 -0.7612 -0.2489 0.263 -0.6897 -0.3959 -0.5995 -0.811 0.5361 -0.2489 -0.2489; 0.2265 -0.328 -0.328 -0.328 -0.328 -0.328 -0.454 -0.328 -0.328 -0.328 -0.328 -0.328 -0.7612 -0.04394 -0.7131 -0.811 -0.415 -0.1361 -0.7859 0.6135; -0.811 -0.4688 -0.304 -0.811 -0.4688 -0.3771 -0.7612 -0.4688 -0.152 -0.4688 -0.4688 -0.4688 -0.4344 -0.4688 -0.6666 -0.4688 -0.4688 -0.4688 -0.4688 0.09761; -0.3384 -0.3384 0.4957 -0.3384 -0.3384 -0.0145 -0.3384 -0.3384 -0.3384 -0.3384 0.1506 -0.3384 -0.3384 -0.811 -0.5146 -0.7859 -0.7131 -0.3384 -0.3384 -0.5146; -0.3846 -0.3846 -0.3846 -0.3846 -0.2515 -0.737 -0.3846 -0.5564 -0.3846 0.09761 -0.3846 -0.3846 -0.3846 -0.3846 -0.3846 -0.3846 -0.3846 -0.4739 -0.1047 -0.6666; -0.4036 -0.4036 -0.4036 -0.4036 -0.5778 -0.4036 -0.4036 -0.7859 -0.7612 -0.4036 -0.7859 -0.4036 -0.4036 -0.4036 -0.4036 0.275 -0.4036 0.2141 -0.4036 -0.4036; -0.6108 -0.6108 -0.6108 -0.7131 -0.6108 -0.6108 -0.6215 -0.4739 -0.4344 -0.6108 -0.6108 -0.6108 -0.6108 -0.6108 -0.6108 -0.811 -0.6108 -0.6108 -0.6108 -0.6108; 0.1243 0.251 -0.2699 -0.6439 -0.2699 -0.2699 -0.2699 -0.2699 -0.2699 -0.2699 -0.2699 -0.2699 -0.2699 -0.2699 -0.811 -0.2699 -0.2699 -0.2699 -0.2699 -0.2699; -0.7074 -0.6215 -0.7074 -0.7074 -0.7074 -0.7074 -0.811 -0.7074 -0.7074 -0.7074 -0.7074 -0.7074 -0.7074 -0.6897 -0.7074 -0.7074 -0.7074 -0.7074 -0.7074 -0.7074; -0.2445 -0.2445 -0.2445 -0.2445 -0.2445 -0.2445 -0.2445 0.3219 -0.2445 -0.811 -0.2445 -0.2445 -0.2445 -0.2445 -0.2445 -0.2445 -0.2445 -0.2445 -0.2445 -0.2445; -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131 -0.7131; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]

	filepath = joinpath(@__DIR__, "data/data_missing.tsv")

	varIds = [Symbol("V$i") for i=1:40]
	sampleIds = [string('S',lpad(i,2,'0')) for i=1:20]
	groupAnnot = repeat(["A","B"],inner=10)
	timeAnnot = vcat(0.01:0.01:0.1, 0.01:0.01:0.1)

	app = MiniApp()
	init(app)

	# Open file
	setvalue(app, Any["samplefilepath", filepath])
	setvalue(app, Any["lastsampleannot", "Time"])
	setvalue(app, Any["normalizemethod", "None"])
	put!(app.fromGUI, :normalize=>[])
	@test runall(app)
	sampleData = app.jg.scheduler.jobs[app.jg.normalizeID].result
	@test sampleData.data ≈ reconstructed

	setvalue(app, Any["logtransform", "true"])
	setvalue(app, Any["logtransformoffset", "0.57"])
	put!(app.fromGUI, :normalize=>[])
	@test runall(app)
	sampleData = app.jg.scheduler.jobs[app.jg.normalizeID].result
	@test sampleData.data ≈ logReconstructed atol=1e-3
end

end


text:
# ----------------------------------------------------------------------
# type
struct SimpleLockFile
    path::String
end

SimpleLockFile() = SimpleLockFile(tempname())

lock_path(slf::SimpleLockFile) = slf.path

# ----------------------------------------------------------------------
# Base
import Base.isfile
isfile(slf::SimpleLockFile) = isfile(lock_path(slf))

import Base.rm
rm(slf::SimpleLockFile; kwargs...) = rm(lock_path(slf); kwargs...)

import Base.basename
basename(slf::SimpleLockFile) = basename(lock_path(slf))

import Base.dirname
dirname(slf::SimpleLockFile) = dirname(dirname(lock_path(slf)))

import Base.mkpath
mkpath(slf::SimpleLockFile) = mkpath(dirname(lock_path(slf)))

# ----------------------------------------------------------------------
# id
rand_lkid(n = 10) = randstring(n)

# ----------------------------------------------------------------------
# write

const _LOCK_FILE_SEP = ","

_validate_id(lock_id::String) =
    contains(lock_id, _LOCK_FILE_SEP) &&
        error("Separator '", _LOCK_FILE_SEP, "' found in the lock id")

const _LOCK_DFT_TIME_OUT = 0.0
const _LOCK_DFT_WAIT_TIME = 1.0
const _LOCK_DFT_VALID_TIME = 30.0
const _LOCK_DFT_CHECK_TIME = 0.1

function _write_lock_file(lf::String;
        lkid::String = rand_lkid(),
        vtime::Float64 = _LOCK_DFT_VALID_TIME
    )
    _validate_id(lkid)
    mkpath(dirname(lf))
    ttag = time() + vtime
    write(lf, string(lkid, _LOCK_FILE_SEP, ttag))
    return (lkid, ttag)
end

write_lock_file(slf::SimpleLockFile; kwargs...) =
    _write_lock_file(lock_path(slf); kwargs...)

# ----------------------------------------------------------------------
# read
function _tryparse(type::Type{T}, str::AbstractString, dft::T) where {T}
    val = tryparse(type, str)
    return isnothing(val) ? dft : val
end

function _read_lock_file(lf::String)
    try
        !isfile(lf) && return ("", -1.0)
        txt = read(lf, String)
        spt = split(txt, _LOCK_FILE_SEP)
        length(spt) != 2 && return ("", -1.0)
        lkid, ttag_str = spt
        ttag = _tryparse(Float64, ttag_str, -1.0)
        return (lkid, ttag)
    catch ignore; end
    return ("", -1.0)
end
read_lock_file(slf::SimpleLockFile) = _read_lock_file(lock_path(slf))

# ----------------------------------------------------------------------
# has lock

_is_valid_ttag(ttag) = ttag > time()

function _is_locked(lf::String, lkid::String)

    !isfile(lf) && return false

    # read
    curr_lid, ttag = _read_lock_file(lf)

    # del if invalid
    if !_is_valid_ttag(ttag)
        _force_unlock(lf)
        return false
    end

    # test
    return lkid == curr_lid
end

"""
    is_locked(slf::SimpleLockFile, lkid::String)::Bool

Check if the lock is valid and is owned by `lkid`
"""
is_locked(slf::SimpleLockFile, lkid::String) = _is_locked(lock_path(slf), lkid)

# ----------------------------------------------------------------------
# release
function _unlock(lf::String, lkid::String)
    !isfile(lf) && return false
    !_is_locked(lf, lkid) && return false
    _force_unlock(lf)
    return true
end

import Base.unlock
"""
    unlock(slf::SimpleLockFile, lkid::String)

Releases ownership of the lock (if `lkid` is valid).
"""
unlock(slf::SimpleLockFile, lkid::String) = _unlock(lock_path(slf), lkid)

# ----------------------------------------------------------------------
# force_unlock
function _force_unlock(lf::String)
    rm(lf; force = true)
    return nothing
end

"""
    force_unlock(slf::SimpleLockFile)

Releases ownership of the lock (even if it is valid).
"""
force_unlock(slf::SimpleLockFile) = _force_unlock(lock_path(slf))

# ----------------------------------------------------------------------
# acquire_lock

function _acquire(lf::String, lkid::String = rand_lkid();
        vtime = _LOCK_DFT_VALID_TIME
    )
    if isfile(lf)
        curr_lid, ttag = _read_lock_file(lf)

        # check if is taken
        if _is_valid_ttag(ttag)
            return (curr_lid == lkid) ?
                (curr_lid, ttag) : # is mine
                ("", ttag) # is taken
        else
            # del if invalid
            _force_unlock(lf)
        end
    end
    return _write_lock_file(lf; lkid, vtime)
end

function _acquire_lock(lf::String, lkid::String = rand_lkid();
        vtime = _LOCK_DFT_VALID_TIME,
        wtime = _LOCK_DFT_WAIT_TIME,
        tout = _LOCK_DFT_TIME_OUT,
        ctime = _LOCK_DFT_CHECK_TIME,
        force = false
    )
    ctime = abs(ctime)

    if tout > 0.0
        t0 = time()

        # try to _acquire till tout
        while true
            lkid0, ttag = _acquire(lf, lkid; vtime)
            if !isempty(lkid0)
                sleep(ctime) # wait before confirm
                _is_locked(lf, lkid) && return (lkid, ttag)
            end

            if (time() - t0) > tout
                if force
                    _force_unlock(lf)
                    return _acquire(lf, lkid; vtime)
                else
                    return ("", ttag)
                end
            end
            sleep(wtime)
        end
    else
        force && _force_unlock(lf)
        return _acquire(lf, lkid; vtime)
    end
end

acquire_lock(slf::SimpleLockFile, lkid::String = rand_lkid(); kwargs...) = _acquire_lock(lock_path(slf), lkid; kwargs...)

# ----------------------------------------------------------------------
# Base.lock

import Base.lock
"""
    lock(f::Function, slf::SimpleLockFile, lkid::String = rand_lkid();
        vtime = $(_LOCK_DFT_VALID_TIME),
        wtime = $(_LOCK_DFT_WAIT_TIME),
        tout = $(_LOCK_DFT_TIME_OUT),
        ctime = $(_LOCK_DFT_CHECK_TIME),
        force = false
    )

Acquire the lock, execute `f()` with the lock held, and release the lock when f returns.
If the lock is already locked by a different `lkid`, the method waits for it to become available and then executes `f()`.
During waiting, it will sleep `wtime` seconds before re-attempting to acquire the lock again and again till `tout` seconds elapsed.
Once acquired the lock the method will wait `ctime` and double check that it is still secure (to reduce races).
The acquisition will take at least `ctime` seconds and if the double check fails it will keep attempting it till `tout`.
If acquired, the lock will be consider valid for `vtime` seconds (no other process/thread asking politely for the lock should get it).
(_**WARNING**_: If `f()` execution time is greater than `vtime`, the lock can be acquired legally by other owner before it finished).
If `force = true` it will forcefully acquire the lock anyway (after `tout`) and then executes `f()` (this is a deadlock free configuration).
This method is not fully secure from racing, but it must be ok for slow applications.
Returns `true` if the locking process was successful (`f()` finished and the lock is still valid).
"""

function lock(
        f::Function, slf::SimpleLockFile, lkid::String = rand_lkid();
        vtime = _LOCK_DFT_VALID_TIME,
        wtime = _LOCK_DFT_WAIT_TIME,
        tout = _LOCK_DFT_TIME_OUT,
        ctime = _LOCK_DFT_CHECK_TIME,
        force = false
    )

    lf = lock_path(slf)
    ok_flag = false
    try
        _acquire_lock(lf, lkid; force, vtime, wtime, tout, ctime)
        ok_flag = _is_locked(lf, lkid)
        ok_flag && f()
    finally
        ok_flag = _is_locked(lf, lkid)
        _unlock(lf, lkid)
    end

    return ok_flag
end


text:
#!/usr/bin/julia

include("Common.jl")
include("SizeModel.jl")
include("IntervalModel.jl")
include("Analysis.jl")

#Analysis.run()

max_stars_repo_name: efficient/msls-eval

max_issues_repo_name: theopengroup/msls-eval

max_forks_repo_name: theopengroup/msls-eval


text:
using MDTable
using Test
using DataFrames: DataFrame, size, names

@testset "Read Tests" begin
    let loaded = readMDTable("file.md")
        @test loaded isa DataFrame
        @test size(loaded) == (2, 4)
    end

    let err = nothing
        try
            readMDTable("fileErr.md")
        catch err
        end

        @test err isa Exception
        @test sprint(showerror, err) == "The number of values in a line must be the same as the number of headers of the file"
    end

    let noHeader = readMDTable("fileNoHeader.md", header=false)
        @test noHeader isa DataFrame
        @test names(noHeader) == ["x1", "x2", "x3", "x4"]
    end
end


text:
    using PyPlot

    x = admmst.x
    mydata = skyst.mydata
    sky = skyst.sky
    nu = psfst.nu
    nu0 = psfst.nu0

    #############################################################################
    #############################################################################
    ######################################  sky / dirty / rec_obj / error
    figure(1)
    clf()
    title("multiplot")
    e = 1
    for z in [1 5 10 15]
        subplot(4,4,e)
        axis("off")
        a = trunc(nu[z]/1e9,3)
        imshow(sky[:,:,z],vmin=0,vmax=1.3)
        title(string(L"$\nu \, = \,$","$a ", L"\, GHz"),fontsize=14)
        e += 1
    end
    e = 5
    for z in [1 5 10 15]
        subplot(4,4,e)
        axis("off")
        imshow(mydata[:,:,z]./maximum(mydata).*1.2,vmin=0,vmax=1.3)
        e += 1
    end
    e = 9
    for z in [1 5 10 15]
        subplot(4,4,e)
        axis("off")
        imshow(x[:,:,z],vmin=0,vmax=1.3)
        e += 1
    end

    subplots_adjust(bottom=0.1, right=0.8, top=0.9)
    cax = axes([0.85, 0.3, 0.025, 0.6])
    colorbar(cax=cax)

    e = 13
    for z in [1 5 10 15]
        subplot(4,4,e)
        axis("off")
        imshow(((sky[:,:,z]-x[:,:,z])./sum(sky[:,:,z])),vmin = 0)
        e += 1
    end

    subplots_adjust(bottom=0.1, right=0.8, top=0.9)
    cax = axes([0.85, 0.1, 0.025, 0.17])
    colorbar(cax=cax)


text:
# ---
# title: 944. Delete Columns to Make Sorted
# id: problem944
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Greedy
# link: <https://leetcode.com/problems/delete-columns-to-make-sorted/description/>
# hidden: true
# ---
#
# We are given an array `A` of `N` lowercase letter strings, all of the same
# length.
#
# Now, we may choose any set of deletion indices, and for each string, we delete
# all the characters in those indices.
#
# For example, if we have an array `A = ["abcdef","uvwxyz"]` and deletion
# indices `{0, 2, 3}`, then the final array after deletions is `["bef", "vyz"]`,
# and the remaining columns of `A` are `["b","v"]`, `["e","y"]`, and
# `["f","z"]`.  (Formally, the `c`-th column is `[A[0][c], A[1][c], ...,
# A[A.length-1][c]]`).
#
# Suppose we chose a set of deletion indices `D` such that after deletions, each
# remaining column in A is in **non-decreasing** sorted order.
#
# Return the minimum possible value of `D.length`.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = ["cba","daf","ghi"]
#     Output: 1
#     Explanation:
#     After choosing D = {1}, each column ["c","d","g"] and ["a","f","i"] are in non-decreasing sorted order.
#     If we chose D = {}, then a column ["b","a","h"] would not be in non-decreasing sorted order.
#
#
# **Example 2:**
#
#
#
#     Input: A = ["a","b"]
#     Output: 0
#     Explanation: D = {}
#
#
# **Example 3:**
#
#
#
#     Input: A = ["zyx","wvu","tsr"]
#     Output: 3
#     Explanation: D = {0, 1, 2}
#
#
#
#
# **Constraints:**
#
#   * `1 <= A.length <= 100`
#   * `1 <= A[i].length <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
abstract type AbstractColors <: Abstractstr end
@enum Colors::String begin
    RED = "red"
    GREEN = "green"
    BLUE = "blue"
end
function show()
    color_map = Dict(Colors.RED => "1", Colors.GREEN => "2", Colors.BLUE => "3")
    a = Colors.GREEN
    if a == Colors.GREEN
        println("green")
    else
        println("Not green")
    end
    println(length(color_map))
end

if abspath(PROGRAM_FILE) == @__FILE__
    show()
end


text:
struct Q
    x::Array{Float64, 1}
end

x = rand(1000);
q = Q(x)

function fx(x::Array{Float64})
    for i in eachindex(x)
        x[i] = exp(-2x[i])
    end
end

function fq(q)
    for i in eachindex(q.x)
        q.x[i] = exp(-2q.x[i])
    end
end

fx1(x) = @time fx(x);
fx1(x);

fq1(q) = @time fq(q);
fq1(q);

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev