text:
# -------------------------------------
# Results and related sub structures
# -------------------------------------

"""
	ResultTimes

Part of the Result object returned by the solver. ResultTimes contains timing results for certain parts of the algorithm:

Time Name  | Description
---  | :---
solver_time | Total time used to solve the problem
setup_time |  Setup time = graph\\_time + init\\_factor\\_time + scaling\\_time
scaling_time | Time to scale the problem data
graph_time | Time used to perform chordal decomposition
init\\_factor\\_time | Time used for initial factorisation of the system of linear equations
factor\\_update\\_time | Sum of times used to refactor the system of linear equations due to rho
iter_time | Time spent in iteration loop
proj_time | Time spent in projection functions
post_time | Time used for post processing
update_time | Time spent in the update! function of the accelerator
accelerate_time | Time spent in the accelerate! function of the accelerator

By default COSMO only measures `solver_time`, `setup_time` and `proj_time`. To measure the other times set `verbose_timing = true`.
"""
mutable struct ResultTimes{T <: AbstractFloat}
	solver_time::T
	setup_time::T
	scaling_time::T
	graph_time::T
	init_factor_time::T
	factor_update_time::T
	iter_time::T
	proj_time::T
	post_time::T
	update_time::T
	accelerate_time::T
end

ResultTimes{T}() where{T} = ResultTimes{T}(T(NaN), T(NaN), T(NaN), T(NaN), T(NaN), T(NaN), T(NaN), T(NaN), T(NaN), T(NaN), T(NaN))
ResultTimes(T::Type = DefaultFloat) = ResultTimes{T}()

function Base.show(io::IO, obj::ResultTimes)
  obj.iter_time != 0 ? verbose = true : verbose = false
  print(io,"Solver time:\t$(round.(obj.solver_time, digits = 4))s ($(round.(obj.solver_time * 1000, digits = 2))ms)\n",
"Setup time:\t$(round.(obj.setup_time, digits = 4))s ($(round.(obj.setup_time * 1000, digits = 2))ms)\n",
"Proj time:\t$(round.(obj.proj_time, digits = 4))s ($(round.(obj.proj_time * 1000, digits = 2))ms)\n")
  if verbose
    print(io,"Iter time:\t$(round.(obj.iter_time, digits = 4))s ($(round.(obj.iter_time * 1000, digits = 2))ms)\n",
    "Scaling time:\t$(round.(obj.scaling_time, digits = 4))s ($(round.(obj.scaling_time * 1000, digits = 2))ms)\n",
    "Graph time:\t$(round.(obj.graph_time, digits = 4))s ($(round.(obj.graph_time * 1000, digits = 2))ms)\n",
    "Initial Factor time:\t$(round.(obj.init_factor_time, digits = 4))s ($(round.(obj.init_factor_time * 1000, digits = 2))ms)\n",
    "Factor update time:\t$(round.(obj.factor_update_time, digits = 4))s ($(round.(obj.factor_update_time * 1000, digits = 2))ms)\n",
    "Post time:\t$(round.(obj.post_time, digits = 4))s ($(round.(obj.post_time * 1000, digits = 2))ms)\n",
    "Accelerator - update time:\t$(round.(obj.update_time, digits = 4))s ($(round.(obj.update_time * 1000, digits = 2))ms)\n",
    "Accelerator - accelerate time:\t$(round.(obj.accelerate_time, digits = 4))s ($(round.(obj.accelerate_time * 1000, digits = 2))ms)\n")
  end
end

"""
    ResultInfo{T <: AbstractFloat}

Object that contains further information about the primal residual, the dual residuals and the rho updates.
"""
struct ResultInfo{T <: AbstractFloat}
	r_prim::T
	r_dual::T
	max_norm_prim::T
	max_norm_dual::T
	rho_updates::Vector{T}
end

ResultInfo(rp, rd, ro, rho_updates) = ResultInfo{DefaultFloat}(rp, rd, rho_updates)

"""
    Result{T <: AbstractFloat}

Object returned by the COSMO solver after calling `optimize!(model)`. It has the following fields:

Fieldname | Type | Description
---  | :--- | :---
x | Vector{T}| Primal variable
y | Vector{T}| Dual variable
s | Vector{T}| (Primal) set variable
obj_val | T | Objective value
iter | Int | Total number of ADMM iterations (incl. safeguarding_iter)
safeguarding_iter | Int | Number of iterations due to safeguarding of accelerator
status | Symbol | Solution status
info | COSMO.ResultInfo | Struct with more information
times | COSMO.ResultTimes | Struct with several measured times
"""
struct Result{T <: AbstractFloat}
    x::Vector{T}
    y::Vector{T}
    s::Vector{T}
    obj_val::T
    iter::Int
    safeguarding_iter::Int
    status::Symbol
    info::ResultInfo
    times::ResultTimes

    function Result{T}() where {T <: AbstractFloat}
      return new(zeros(T, 1), zeros(T, 1), zeros(T, 1), zero(T), 0, 0, :Unsolved, ResultInfo{T}(0., 0., 0., 0., T[]), ResultTimes{T}())
    end

    function Result{T}(x, y, s, obj_val, iter, safeguarding_iter, status, info, times) where {T <: AbstractFloat}
      return new(x, y, s, obj_val, iter, safeguarding_iter, status, info, times)
    end

end

function Base.show(io::IO, obj::Result)
	print(io,">>> COSMO - Results\nStatus: ")
	if obj.status == :Solved
		result_color = :green
	else
		result_color = :red
	end
	printstyled("$(obj.status)\n", color = result_color)
	println("Iterations: $(obj.iter) (incl. $(obj.safeguarding_iter) safeguarding iterations)\nOptimal Objective: $(@sprintf("%.4g", obj.obj_val))\nRuntime: $(round.(obj.times.solver_time * 1000, digits = 2))ms\nSetup Time: $(round.(obj.times.setup_time * 1000, digits = 2))ms\n")
	!isnan(obj.times.iter_time) && print("Avg Iter Time: $(round.((obj.times.iter_time / obj.iter) * 1000, digits = 2))ms")
end

# -------------------------------------
# Problem scaling
# -------------------------------------

struct ScaleMatrices{Tf <: AbstractFloat}
	D::Union{UniformScaling{Bool}, Diagonal{Tf, Vector{Tf}} }
	Dinv::Union{UniformScaling{Bool}, Diagonal{Tf, Vector{Tf}} }
	E::Union{UniformScaling{Bool}, Diagonal{Tf, Vector{Tf}} }
	Einv::Union{UniformScaling{Bool}, Diagonal{Tf, Vector{Tf}} }
	c::Base.RefValue{Tf}
	cinv::Base.RefValue{Tf}
end

ScaleMatrices(args...) = ScaleMatrices{DefaultFloat}(args...)

ScaleMatrices{T}() where {T} = ScaleMatrices(I, I, I, I, Base.RefValue{T}(one(T)), Base.RefValue{T}(one(T)))

function ScaleMatrices{T}(m, n) where {T <: AbstractFloat}
	D    = Diagonal(ones(T, n))
	Dinv = Diagonal(ones(T, n))
	E    = Diagonal(ones(T, m))
	Einv = Diagonal(ones(T, m))
	c    = Base.RefValue{T}(one(T))
	cinv = Base.RefValue{T}(one(T))
	ScaleMatrices(D, Dinv, E, Einv, c, cinv)
end

# -------------------------------------
# Problem data
# -------------------------------------

mutable struct ProblemData{T <: AbstractFloat}
	P::AbstractMatrix{T}
	q::Vector{T}
	A::AbstractMatrix{T}
	b::Vector{T}
	C::CompositeConvexSet{T}
	model_size::Array{Integer,1}

	function ProblemData{T}() where {T <: AbstractFloat}
		return new(
			spzeros(T, 1, 1),             #P
			T[],                        #q
			spzeros(T, 1, 1),             #A
			T[],                        #b
			COSMO.CompositeConvexSet{T}([COSMO.ZeroSet{T}(1)]),     #C
			[0; 0])                 #model size
	end
end

ProblemData(args...) = ProblemData{DefaultFloat}(args...)

# ---------------------------
# Struct to hold clique and sparsity data for a constraint
# ---------------------------

mutable struct SparsityPattern
  sntree::SuperNodeTree
  ordering::Array{Int}
  reverse_ordering::Array{Int}
  row_range::UnitRange{Int} # the starting row of the psd cone in the original problem
  cone_ind::Int # this is the ind of the original psd cone in ws.p.C that is decomposed
  nz_ind_map::SparseVector{Int, Int} # maps a matrix entry k = svec(i, j) to the location of the entry in the sparse data structure

  # constructor for sparsity pattern
  function SparsityPattern(L::SparseMatrixCSC, N::Int, ordering::Array{Int, 1}, merge_strategy::Union{Type{<: AbstractMergeStrategy}, OptionsFactory{<: AbstractMergeStrategy}}, row_range::UnitRange{Int}, cone_ind::Int, nz_ind_map::SparseVector{Int, Int})

    merge_strategy = merge_strategy()
    sntree = SuperNodeTree(L, merge_strategy)

    # clique merging
    sntree.num > 1 && merge_cliques!(sntree)

    # reorder vertices in supernodes to have consecutive order
    # necessary for equal column structure for psd completion
    reorder_snd_consecutively!(sntree, ordering)

    # undo the reordering and sort
    # for iii = 1:sntree.num
    #   sep = get_sep(sntree, iii)
    #   snd = get_snd(sntree, iii)
    #   map!(v -> ordering[v], sep)
    #   map!(v -> ordering[v], snd)
    #   sort!(sep)
    #   sort!(snd)
    # end

    # for each clique determine the number of entries of the block represented by that clique
    calculate_block_dimensions!(sntree)#, merge_strategy)

    return new(sntree, ordering, invperm(ordering), row_range, cone_ind, nz_ind_map)
  end

  # For debugging
  function SparsityPattern(sntree::SuperNodeTree, ordering::Array{Int}, reverse_ordering::Array{Int}, row_range::UnitRange{Int}, cone_ind::Int)
    return new(sntree, ordering, reverse_ordering, row_range, cone_ind)
  end
end

# -------------------------------------
# Chordal Decomposition Information
# -------------------------------------
mutable struct ChordalInfo{T <: AbstractFloat}
  decompose::Bool # an internal flag to check if problem has been decomposed
  originalM::Int
  originalN::Int
  originalC::CompositeConvexSet{T}
  H::SparseMatrixCSC{T}
  sp_arr::Array{COSMO.SparsityPattern}
  psd_cones_ind::Array{Int} # stores the position of decomposable psd cones in the composite convex set
  num_psd_cones::Int # number of psd cones of original problem
  num_decomposable::Int #number of decomposable cones
  num_decom_psd_cones::Int #total number of psd cones after decomposition
  L::SparseMatrixCSC{T} #pre allocate memory for QDLDL
  cone_map::Dict{Int, Int} # map every cone in the decomposed problem to the equivalent or undecomposed cone in the original problem

  function ChordalInfo{T}(problem::COSMO.ProblemData{T}, settings::COSMO.Settings) where {T <: AbstractFloat}
    originalM = problem.model_size[1]
    originalN = problem.model_size[2]
    originalC = deepcopy(problem.C)
    num_psd_cones = length(findall(x -> typeof(x) <: Union{PsdConeTriangle{T}, PsdCone{T}} , problem.C.sets))
    # allocate sparsity pattern for each cone
    sp_arr = Array{COSMO.SparsityPattern}(undef, num_psd_cones)
    cone_map = Dict{Int, Int}()

    return new(settings.decompose, originalM, originalN, originalC, spzeros(1, 1), sp_arr, Int[], num_psd_cones, 0, 0, spzeros(1, 1), cone_map)
  end

	function ChordalInfo{T}() where{T}
		C = COSMO.CompositeConvexSet{T}([COSMO.ZeroSet{T}(1)])
		return new(false, 0, 0, C, spzeros(1, 1), COSMO.SparsityPattern[], [1])
	end

end

# -------------------------------------
# Structure of internal iterate variables
# -------------------------------------

struct Variables{T}
	w::Vector{T}
	w_prev::Vector{T}
	x::SubArray{T}
	s::SplitVector{T}
	μ::Vector{T}

	function Variables{T}(m::Int, n::Int, C::AbstractConvexSet{T}) where {T <: AbstractFloat}
		m == C.dim || throw(DimensionMismatch("set dimension is not m"))
		w = zeros(T, n + m)
		w_prev = zeros(T, n + m)
		x = view(w_prev, 1:n) #x is a view onto w_prev, to keep it in sync with s and μ
		s = SplitVector(zeros(T, m), C)
		μ = zeros(T, m)
		new(w, w_prev, x, s, μ)
	end
end

Variables(args...) = Variables{DefaultFloat}(args...)

mutable struct IterateHistory
	x_data::AbstractMatrix
	s_data::AbstractMatrix
	y_data::AbstractMatrix
	v_data::AbstractMatrix
	r_prim_data::AbstractVector
	r_dual_data::AbstractVector
	eta_data::AbstractArray
	alpha_data::AbstractArray
	cond_data::AbstractVector
	aa_fail_data::AbstractVector

	function IterateHistory(m, n, mem)
		new(zeros(n, 0), zeros(m, 0), zeros(m, 0), zeros(m + n, 0), Float64[], Float64[], zeros(mem, 0), zeros(mem + 1, 0), Float64[], Int64[])
	end
end

function update_iterate_history!(history::IterateHistory, x, s, y, v, r_prim, r_dual, eta::Vector{Float64}, cond::Float64)
	history.x_data = hcat(history.x_data, x)
	history.s_data = hcat(history.s_data, s)
	history.y_data = hcat(history.y_data, y)
	history.v_data = hcat(history.v_data, v)
	history.cond_data = push!(history.cond_data, cond)
	push!(history.r_prim_data, r_prim)
	push!(history.r_dual_data, r_dual)

	# alphas = compute_alphas(eta)
	# history.alpha_data = hcat(history.alpha_data, alphas)
	# history.eta_data = hcat(history.eta_data, eta)
end

struct UtilityVariables{T}
  vec_m::Vector{T}
  vec_n::Vector{T}
  vec_n2::Vector{T}

  function UtilityVariables{T}(m::Int, n::Int) where {T <: AbstractFloat}
    new(zeros(T, m), zeros(T, n), zeros(T, n))
  end
end

UtilityVariables(args...) = UtilityVariables{DefaultFloat}(args...)

# -------------------------------------
# a collection of state flags
# -------------------------------------

mutable struct States
	IS_ASSEMBLED::Bool # the workspace has been assembled with problem data
	IS_OPTIMIZED::Bool # the optimization function has been called on the model
	IS_CHORDAL_DECOMPOSED::Bool # the problem has been decomposed
	KKT_FACTORED::Bool # the KKT matrix has been factored
	IS_SCALED::Bool # the problem data has been scaled
	States() = new(false, false, false, false, false)
end

# -------------------------------------
# Top level container for all solver data
# -------------------------------------
"""
	Workspace{T <: AbstractFloat}()

Initializes an empty COSMO model that can be filled with problem data using `assemble!(model, P, q,constraints; [settings, x0, s0, y0])`.
"""
mutable struct Workspace{T}
	p::ProblemData{T}
	settings::Settings{T}
	sm::ScaleMatrices{T}
	ci::ChordalInfo{T}
	vars::Variables{T}
  	utility_vars::UtilityVariables{T}
	δx::Vector{T}
	δy::SplitVector{T}
	s_tl::Vector{T}
	ls::Vector{T}
	sol::Vector{T}
	x_tl::SubArray{T}
	ν::SubArray{T}
	ρ::T
	ρvec::Vector{T}
	kkt_solver::Union{AbstractKKTSolver,Nothing}
	states::States
	rho_updates::Vector{T} #keep track of the rho updates and the number of refactorisations
	rho_update_due::Bool # a flag that indicates that at the next possible iteration ρ should be updated
	infeasibility_check_due::Bool # a flag that indicates that at the next possible iteration we should check for infeasibility
	times::ResultTimes{Float64} #always 64 bit regardless of data type
	row_ranges::Array{UnitRange{Int}, 1} # store a set_ind -> row_range map
	safeguarding_iter::Int # count extra ADMM iterations due to bad quality accelerated steps
	accelerator::AbstractAccelerator
	accelerator_active::Bool
	activation_reason::AbstractActivationReason
	#constructor
	function Workspace{T}() where {T <: AbstractFloat}
		p = ProblemData{T}()
		sm = ScaleMatrices{T}()
		vars = Variables{T}(1, 1, p.C)
    	uvars = UtilityVariables{T}(1, 1)
		ci = ChordalInfo{T}()
		δx = zeros(T, 0)
		δy = SplitVector(zeros(T, 1), p.C)
		s_tl = zeros(T,0)
		ls = zeros(T, 0)
		sol = zeros(T, 1)
		x_tl = view(sol, 1:1)
		ν = view(sol, 1:1)
		return new(p, Settings{T}(), sm, ci, vars,  uvars, δx, δy, s_tl, ls, sol, x_tl, ν, zero(T), T[], nothing, States(), T[], false, false, ResultTimes(), [0:0], 0, EmptyAccelerator(), false, ImmediateActivation())
	end
end
Workspace(args...) = Workspace{DefaultFloat}(args...)

Base.show(io::IO, model::COSMO.Workspace{T}) where {T} = println(io, "A COSMO Model with Float precision: $(T)")

# Type alias facing the user
"""
	Model{T <: AbstractFloat}()

Initializes an empty COSMO model that can be filled with problem data using `assemble!(model, P, q,constraints; [settings, x0, s0, y0])`.
"""
const Model = Workspace;


text:
DECODER = GNSSDecoderState(PRN = 1)
PREV_29 = 0
PREV_30 = 0
BUFFER_TEST_TRUE = [
[
    BitArray([1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1]),
    BitArray([1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0]),
    BitArray([0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
    BitArray([0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1]),
    BitArray([0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0]),
    BitArray([1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1]),
    BitArray([0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0]),
    BitArray([0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1]),
    BitArray([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0]),
    BitArray([1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0])
    ],

[
    BitArray([1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1]),
    BitArray([1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0]),
    BitArray([0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0]),
    BitArray([0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1]),
    BitArray([1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1]),
    BitArray([1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0]),
    BitArray([1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0]),
    BitArray([0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1]),
    BitArray([1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1]),
    BitArray([1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0])
    ],

[
    BitArray([1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1]),
    BitArray([1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0]),
    BitArray([1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0]),
    BitArray([0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1]),
    BitArray([0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1]),
    BitArray([0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1]),
    BitArray([1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1]),
    BitArray([0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1]),
    BitArray([0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1]),
    BitArray([1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0])
    ],

[
    BitArray([1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1]),
    BitArray([1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0]),
    BitArray([0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1]),
    BitArray([0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1]),
    BitArray([0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1]),
    BitArray([0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1]),
    BitArray([0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0])
    ],

[
    BitArray([1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1]),
    BitArray([1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0]),
    BitArray([0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0]),
    BitArray([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0])
    ]
]

ext: jl


text:
__precompile__()
module DTWDT

include("functions.jl")
using .DTWDTfunctions

export dtwdt

"""
    dtwdt(u0::Union{Array{Float32,1},Array{Float64,1}}, u1::Union{Array{Float32,1},Array{Float64,1}}, dt::Float64; dtwnorm::String="L2"
        maxlag::Int64=80, b::Int64=1, direction::Int64=1)

returns minimum distance time lag and index in dist array, and dtw error between traces.

# Arguments
- `u0, u1::Union{Array{Float32,1},Array{Float64,1}}`: Time series.
- `dt::Float64`: time step (dt of u0 and u1 should be same)
- `dtwnorm::String`: norm to calculate distance; effect on the unit of dtw error. (L2 or L1)
-  Note: L2 is not squard, thus distance is calculated as (u1[i]-u0[j])^2
- `maxlag::Int64`: number of maxLag id to search the distance.
- `b::Int64t`: b value to controll in distance calculation algorithm (see Mikesell et al. 2015).
- `direction::Int64`: length of noise data window, in seconds, to cross-correlate.

# Outputs
- `stbarTime::Array{Float64,1}`: series of time shift at t.
- `stbar::Array{Int64,1}`: series of minimum distance index in distance array.
- `dist::Array{Int64,2}`: distance array.
- `dtwerror::Float64`: dtw error (distance) between two time series.

"""
function dtwdt(u0::Union{Array{Float32,1},Array{Float64,1}}, u1::Union{Array{Float32,1},Array{Float64,1}}, dt::Float64;
    dtwnorm::String="L2",    #norm to calculate distance; effect on the unit of dtw error
    maxLag::Int64=80,         #number of maxLag id to search the distance
    b::Int64=1,               #b value to controll in distance calculation algorithm (see Mikesell et al. 2015)
    direction::Int64=1,       #direction to accumulate errors (1=forward, -1=backward, 0=double to smooth)
    )

    #prepare datasize and time vector
    lvec   = (-maxLag:maxLag).*dt; # lag array for plotting below
    npts   = length(u0);            # number of samples
    if length(u0) != length(u1) error("u0 and u1 must be same length.") end
    tvec   = ( 0 : npts-1 ) .* dt; # make the time axis

    #compute distance between traces
    err = computeErrorFunction(u1, u0, npts, maxLag, norm=dtwnorm);

    #compute distance array and backtrack index
    if direction == 1 || direction == -1
        dist  = accumulateErrorFunction(direction, err, npts, maxLag, b); # forward accumulation to make distance function
        stbar = backtrackDistanceFunction( -1*direction, dist, err, -maxLag, b ); # find shifts

    elseif direction == 0
        #calculate double time to smooth distance array
        dist1 = accumulateErrorFunction( -1, err, npts, maxLag, b ); # forward accumulation to make distance function
        dist2 = accumulateErrorFunction( 1, err, npts, maxLag, b ); # backwward accumulation to make distance function
        dist  = dist1 .+ dist2 .- err; # add them and remove 'err' to not count twice (see Hale's paper)
        stbar = backtrackDistanceFunction( -1, dist, err, -maxLag, b );
    else
        error("direction must be +1, -1 or 0(smoothing).")
    end

    stbarTime = stbar .* dt;      # convert from samples to time
    tvec2     = tvec + stbarTime; # make the warped time axis

    #accumulate distance in distance array to calculate dtw error
    error = computeDTWerror( err, stbar, maxLag );

    return (stbarTime, stbar, dist, error)
end

end # module


text:
const ONEOSQPI(::Type{BigFloat}) = big"5.6418958354775628694807945156077258584405E-1"
const TWOOPI(::Type{BigFloat}) = big"6.3661977236758134307553505349005744813784E-1"
#const SQPIO2(::Type{BigFloat}) = big"1.253314137315500251207882642405522626503493370304969158314961788171146827303924"
#const SQ1O2PI(::Type{BigFloat}) = big"0.3989422804014326779399460599343818684758586311649346576659258296706579258993008"
#const SQ2OPI(::Type{BigFloat}) = big"0.7978845608028653558798921198687637369517172623298693153318516593413158517986017"
#const PIO4(::Type{BigFloat}) = big"0.7853981633974483096156608458198757210492923498437764552437361480769541015715495"

const PIO4(::Type{Float64}) = .78539816339744830962
const THPIO4(::Type{Float64}) = 2.35619449019234492885
const SQ2OPI(::Type{Float64}) = .79788456080286535588
const TWOOPI(::Type{Float64}) = 0.6366197723675814
const SQPIO2(::Type{Float64}) = 1.25331413731550025
const SQ1O2PI(::Type{Float64}) = 0.3989422804014327

const PIO4(::Type{Float32}) = 0.78539816339744830962f0
const TWOOPI(::Type{Float32}) = 0.636619772367581343075535f0
const THPIO4(::Type{Float32}) = 2.35619449019234492885f0


text:
## He gradient test
using Flux
include("src/helperFunctions.jl")
include("src/utils.jl")
using .MapInLayer

using Plots
using SmoothLivePlot
using Printf
gr(show = true)

## Plotting functions
function myPlotFun(tforPlot, tDomain, trueTemperature, tempData, initialEstAll, initialEst, finalEstAll, finalEst)
  plot(tforPlot, trueTemperature,
    label = "True",
    xmirror=true,
    legend=:bottomright)
  scatter!(tDomain, tempData, markersize = 5, label = "Data")
  plot!(tforPlot, initialEstAll, color = "green", label = "Original")
  scatter!(tDomain, initialEst, markersize = 5, color = "green", label = "")
  plot!(tforPlot, finalEstAll, color = "red", label = "Final")
  scatter!(tDomain, finalEst, markersize = 5, color = "red", label = "")
  plot!(xaxis = ("Time, Ma", :flip))
  plot!(yaxis = ("Temperature, C", :flip))
  plot!(ylims = (0, 110))
  #display(fig)
end

## Main functions
function trueTemperatureFun(t::Real)
  # Based on QTQt test temperature
  intc = Float32(100)
  midp = Float32(50)
  grad = Float32(8/5)
  Float32(intc - grad*abs(t - midp))
end
function formModel()
  # Set up neural network
  temperatureNet = Chain(
  Dense(1, 16, gelu),
  Dense(16, 32, gelu),
  Dense(32, 64),
  +)

  # NOTE: This function seems to be less effective than inlining
  initalizeWeights!(temperatureNet)

  temperatureMap = MapIn(temperatureNet)
end
function initalizeWeights!(net)
  # Randomize weights
  p, re = Flux.destructure(net)
  p[p.!=0.0] = Float32.(0.5*randn(sum(p.!=0.0)))
  p[p.==0.0] = (50.0 .+ Float32.(10.0*randn(sum(p.==0.0))))
  net = re(p)
  net = MapIn(net)
end
function makeData(tDomain, nTrain)
  tempData = trueTemperatureFun.(tDomain)
  #data = [(tDomain, trueTemperature)]
  data = Iterators.repeated((tDomain, tempData), nTrain)
  return data, tempData
end
function main()

  # True temperature for comparison
  tforPlot = collect(range32(0, 100, 0.1))
  trueTemperature = trueTemperatureFun.(tforPlot)

  # Sample true data
  tDomain = Float32[0.1, 5.3, 15.5, 25.5, 39.9, 45.3, 55.5, 65.5, 89.9, 99.5]
  data, tempData = makeData(tDomain, 200)

  # Model and initial prediction
  temperatureMap = formModel()
  initialEstAll = temperatureMap(tforPlot)
  initialEst = temperatureMap(tDomain)

  # Training requirements
  parameters = params(temperatureMap)
  loss(x, y) = Flux.mae(temperatureMap(x), y)
  opt = ADAM(0.1)
  #anim = Animation()
  evalcb() = begin
    println(loss(data.xs.x...))
    finalEstAll = temperatureMap(tforPlot)
    finalEst = temperatureMap(tDomain)
    # NOTE: Get SmoothLivePlot working with this
    display(myPlotFun(tforPlot, tDomain, trueTemperature, tempData, initialEstAll, initialEst, finalEstAll, finalEst))
    #frame(anim)
  end

  # Start simple training
  #Flux.train!(loss, parameters, data, opt, cb = Flux.throttle(evalcb, 0.1))
  #Flux.train!(loss, parameters, data, ADAM(0.01), cb = Flux.throttle(evalcb, 0.1))
  Flux.train!(loss, parameters, data, opt, cb = Flux.throttle(evalcb, 0.1))
  Flux.train!(loss, parameters, data, ADAM(0.01), cb = Flux.throttle(evalcb, 0.1))

  # Compare results
  finalEstAll = temperatureMap(tforPlot)
  finalEst = temperatureMap(tDomain)
  myPlotFun(tforPlot, tDomain, trueTemperature, tempData, initialEstAll, initialEst, finalEstAll, finalEst)

  #gif(anim, "anim_fps30b.gif", fps = 30)

end

main()

max_stars_repo_name: williamjsdavis/ThermoChroNN

max_issues_repo_name: williamjsdavis/ThermoChroNN

max_forks_repo_name: williamjsdavis/ThermoChroNN


text:
p = ccall(:strdup, Ptr{Cuchar}, (Ptr{Cuchar},), "Hello world")
@show unsafe_string(p) # "Hello world"
ccall(:free, Void, (Ptr{Cuchar},), p)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
println("Compile libraries on main processor..."); flush(stdout)
include("imports.jl")
include("scenarios.jl")
println("Finished loading libraries."); flush(stdout)

simdata = scenarios(1, seed=1, Ni=5000)
data = CDE.Model.Data(simdata[:yC], simdata[:yT])
K = 6
prior = CDE.Model.PriorAM(K, true, phi=Normal(-10, 5), data=data)
println(prior)

warmup(N) = let
  init = CDE.Model.StateAM(data, prior)
  tuner = CDE.Model.MCMC.MvTuner(CDE.Model.tovec(init))
  for n in 1:N
    println("n: $n")
    @time chain, init, summarystats, tuner = CDE.amfit(
      init, data, prior, nsamps=[1], nburn=200, tuner=tuner,
      nc=50*n, nt=50*n, temper=10)
    tuner.iter = 1
  end
  return init, tuner
end
init, tuner = warmup(0)

@time chain, laststate, summarystats, tuner = CDE.amfit(
  init, data, prior, thin=2, nsamps=[2000], nburn=5000,
  tuner=tuner, nc=500, nt=500, temper=1)

function convert_chain!(chain, beta)
  for d in chain[1]
    d[:psi] = CDE.Util.toaltskew.(d[:sigma], d[:phi])
    d[:omega] = CDE.Util.toaltscale.(d[:sigma], d[:phi]) .^ 2
    d[:beta] = beta
  end
end
convert_chain!(chain, true)

imgdir = "results/testam/img"; mkpath(imgdir)
CDE.Model.plotpostsummary([chain[1]], summarystats, simdata[:yC],
                          simdata[:yT], imgdir; digits=3, laststate=laststate,
                          bw_postpred=0.2, ygrid=ygrid, plotsize=(400,400),
                          density_legend_pos=:topleft, simdata=simdata)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module DiscretePDEs
    include("gmsh_interface.jl")
    include("geo.jl")
    include("gds.jl")
    include("electromagnetics.jl")
end


text:
#--------------------------------------------------------------
# MultilevelTT package
# example code for https://arxiv.org/abs/1802.09062
#--------------------------------------------------------------
# Code for 1D Poisson equation example,
# using soft thresholding solver
module Poisson

using TT, TTFEM, TTAlg

export quadraticRefSol, stsolve1D, stsolve2D

# exact solution of 1D test problem
function quadraticRefSol(L::Int64)
  t = Tensor(2, L, 3)
  x = .5;
  t[1] = reshape([0. 0. 1.; x^2 x 1.] , (1,2,3));
  t[1] /= 2*sqrt(2.)
  for i = 2:(L-1)
    x = 2.0^(-i);
    t[i] = reshape([ 1. 0. 0.; 0. 1. 0.; 0. 0. 1.;  1. 0. 0.;  2.*x 1. 0.; x^2 x 1.], (3,2,3))
    t[i] /= sqrt(2.)
  end
  x = 2.0^(-L);
  x0 = x
  x1 = x0 + x
  t[L] = reshape([ -1.; 2.*(1. - x0); x0*(2.-x0); -1.; 2.*(1. - x1); x1*(2.-x1)], (3, 2, 1))
  t[L] /= sqrt(2.)
  return t
end

# soft thresholding iteration with 1D BPX preconditioning,
# reduced-rank operator with intermediate recompression
function stsolve1D(L::Int64, ɛ::Float64)
  normest = 25.
  invnormest = .5

  C = ttBPXdnC(L, 1)
  F = ttBPXdnQLaplace(L, 1, 1)
  Ft = transpose(F)

  println("operator ranks: ", ranks(F)," / ", ranks(Ft))

  f = ttrhsdn(L)
  g = matvec(C, f)
  svdtrunc!(g, svd!(g), .005*ɛ);
  println("rhs ranks ", ranks(g))

  Θ = 0.5

  function residual(u::Tensor, δ::Float64)
    gδ = deepcopy(g)
    svdtrunc!(gδ, svd!(gδ), .1*δ);
    t = matvec(F, u)
    maxrank = [maximum(ranks(t))]
    svdtrunc!(t, svd!(t), .8*δ / sqrt(normest));
    push!(maxrank, maximum(ranks(t)))
    r = matvec(Ft, t)
    push!(maxrank, maximum(ranks(r)))
    add!(r, -1., g)
    push!(maxrank, maximum(ranks(r)))
    svdtrunc!(r, svd!(r), .1*δ);
    push!(maxrank, maximum(ranks(r)))
    return r, maxrank
  end

  u, cdata = inexstsolve(L, 1, residual, normest, invnormest, .995*ɛ, Θ, .95)

  return u, matvec(C, u), cdata
end

# soft thresholding iteration with 2D BPX preconditioning
function stsolve2D(L::Int64, ɛ::Float64)
  normest = 25.
  invnormest = .4

  C = ttBPXdnC(L, 2)

  Q11 = ttBPXdnQLaplace(L, 2, 1)
  Q12 = ttBPXdnQLaplace(L, 2, 2)
  Q11t = transpose(Q11)
  Q12t = transpose(Q12)

  f = ttkron(ttrhsdn(L), ttrhsdn(L))
  g = matvec(C, f)
  svdtrunc!(g, svd!(g), .005*ɛ);
  println("rhs ranks ", ranks(g))

  Θ = 0.5

  function residual(u::Tensor, δ::Float64)
	maxrank = Int64[]
    gδ = deepcopy(g)
    svdtrunc!(gδ, svd!(gδ), δ/7);

    s = matvec(Q11, u)
    svdtrunc!(s, svd!(s), δ / 7 / sqrt(normest));
    push!(maxrank, maximum(ranks(s)))

    r1 = matvec(Q11t, s)
	svdtrunc!(r1, svd!(r1), δ / 7);
    push!(maxrank, maximum(ranks(r1)))

	s = matvec(Q12, u)
    svdtrunc!(s, svd!(s), δ / 7 / sqrt(normest));
	push!(maxrank, maximum(ranks(s)))

	r2 = matvec(Q12t, s)
	svdtrunc!(r2, svd!(r2), δ / 7);
	push!(maxrank, maximum(ranks(r2)))

	add!(r1, r2)
	svdtrunc!(r1, svd!(r1), δ / 7);
	push!(maxrank, maximum(ranks(r1)))

    add!(r1, -1., g)
    svdtrunc!(r1, svd!(r1), δ / 7);
    push!(maxrank, maximum(ranks(r1)))
    return r1, maxrank
  end

  u, cdata = inexstsolve(L, 2, residual, normest, invnormest, .995*ɛ, Θ)

  return u, matvec(C, u), cdata
end

end


text:
# Statistical model for analyzing Radial Velocity observations
# Assumes existance of data in arrays times, obs, sigma_obs
# Assumes existance of functions: TODO

export loglikelihood, logprior_noise, is_valid_noise

num_noise_param() = 1

if(!isdefined(:PDMats))        using PDMats       end

function is_valid_noise(p::Vector)
  param = extract_noise_param(p)
  @assert length(param) == num_noise_param()
  sigma_cor = param[1]
  tau_evolve   = 50.0 # param[2]
  lambda_p     = 0.5  # param[3]
  tau_rotate   = 20.0 # param[4]
  const min_jitter = 0.0
  const max_jitter = 99.0
  #if !(min_jitter<=sigma_sq_cor<=max_jitter) || !(0<tau_evolve) || !(0<lambda_p) || !(0<tau_rotate)
  if !(min_jitter<=sigma_cor<=max_jitter)
     return false
  end
  return true
end

function kernel_quaesiperiodic(dt::Float64,param::Vector)
  sigma_cor = param[1]
  tau_evolve   = 50.0 # param[2]
  lambda_p     = 0.5  # param[3]
  tau_rotate   = 20.0 # param[4]
  sigma_cor*sigma_cor*exp(-0.5*((dt/tau_evolve)^2 + (sin(pi*dt/tau_rotate)/lambda_p)^2) )
end

function loglikelihood(p::Vector)
  num_pl = num_planets(p)
  if num_pl <= 1 return 0.0 end
  if !is_valid(p) return -Inf end  # prempt model evaluation
  # Set t, o, and so to point to global arrays with observational data, while enforcing types
  t::Array{Float64,1} = times
  o::Array{Float64,1} = obs
  so::Array{Float64,1} = sigma_obs
  @assert length(times) == length(obs) == length(sigma_obs)
  gp_param = extract_noise_param(p)
  model_minus_data = map(i->calc_model_rv(p,t[i])-o[i],1:length(t))
  Sigma = PDMat( [ kernel_quaesiperiodic(t[i]-t[j],gp_param) + (i==j ? so[i]^2 : 0.0) for i in 1:length(t), j in 1:length(t)] )

  # chisq = model_minus_data'*(Sigma*model_minus_data)
  chisq = invquad(Sigma, model_minus_data)
  log_likelihood = -0.5*( chi_sq + logdet(Sigma) + length(t)*log(2pi) )
end

function logprior_noise(p::Vector)
  logp = zero(eltype(p))
  if !is_valid_noise(p)
     logp -= Inf
     return logp
  end
  param = extract_noise_param(p)
  @assert length(param) == num_noise_param()
  sigma_cor = sqrt(param[1])
  # tau_evolve   = param[2]
  # lambda_p     = param[3]
  # tau_rotate   = param[4]
  const Jitter0 = 1.0
  const min_jitter = 0.0
  const max_jitter = 1000.0
  const jitter_prior_norm = -log(Jitter0*log1p(max_jitter/Jitter0))
  logp += -log1p(sigma_cor/Jitter0) + jitter_prior_norm
  # TODO: Add prior for other GP parameters if they become variable
  return logp::eltype(p)
end


text:
# RainFARM:  Stochastic downscaling following
# - *D'Onofrio et al. 2014, J of Hydrometeorology 15 , 830-843* and
# - *Rebora et. al 2006, JHM 7, 724*
# Includes orographic corrections
# Implementation in Julia language

# Copyright 2016 Jost von Hardenberg
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
# either express or implied.
# See the License for the specific language governing permissions
# and limitations under the License.

__precompile__()
module RainFARM
export agg, fft3d, initmetagauss, gaussianize, metagauss, smoothconv, smoothspec
export mergespec_spaceonly, downscale_spaceonly, lon_lat_fine, fitslopex
export read_netcdf2d, write_netcdf2d, rainfarm, interpola, smooth
export overwrite_netcdf2d, rfweights

using Interpolations, NetCDF, Compat
using Compat.Statistics, Compat.Printf, Compat.SparseArrays
if VERSION >= v"0.7.0-DEV.2005"
    using FFTW
end

include("rf/agg.jl")
include("rf/smoothconv.jl")
include("rf/smoothspec.jl")
include("rf/smooth.jl")
include("rf/fft3d.jl")
include("rf/initmetagauss.jl")
include("rf/gaussianize.jl")
include("rf/metagauss.jl")
include("rf/mergespec_spaceonly.jl")
include("rf/downscale_spaceonly.jl")
include("rf/lon_lat_fine.jl")
include("rf/fitslopex.jl")
include("rf/read_netcdf2d.jl")
include("rf/write_netcdf2d.jl")
include("rf/overwrite_netcdf2d.jl")
include("rf/interpola.jl")
include("rf/rainfarm.jl")
include("rf/rfweights.jl")

end


text:
# SPDX-License-Identifier: MIT

"""
    loggamma_dot

Wrapper around the log-gamma version to avoid conflicts with CUDA
"""
loggamma_dot(x::AbstractVector) = loggamma.(x)
loggamma_dot(x::CuArray) = lgamma.(x)

"""
    pdmat_logdet(M)

Log-determinant of a positive definite matrix.
"""
function pdmat_logdet(M::AbstractMatrix)
    U = cholesky(M).U
    sum(2 * log.(diag(U)))
end
pdmat_logdet(M::Diagonal) = sum(log.(diag(M)))

"""
    pdmat_inverse(M)

Inverse of a positive definite matrix.
"""
function pdmat_inverse(M::AbstractMatrix)
    U = cholesky(M).U
    I = typeof(M)(LinearAlgebra.I, size(M)...)
    U⁻ = I / U
    U⁻ * U⁻'
end
pdmat_inverse(M::Diagonal) = Diagonal(1 ./ diag(M))

"""
    vec_tril(M)

Returns the vectorized low-triangular part (diagonal not included) of
the matrix.

See also [`inv_vec_tril`](@ref), [`matrix`](@ref).
"""
function vec_tril(M)
    D, _ = size(M)
    vcat([diag(M, -i) for i in 1:D-1]...)
end

"""
    inv_vec_tril(v)

Returns a lower triangular matrix from a vectorized form. The diagonal
of the matrix is set to 0.

See also [`vec_tril`](@ref), [`matrix`](@ref)
"""
function inv_vec_tril(v)
    # Determine the dimension of the matrix from the length of the
    # vector. To do this, we solve the equation:
    #   x²/2 + x/2 - l = 0
    # The dimension is given by s (the positive solution) + 1.
    l = length(v)
    D = Int((-1 + sqrt(1 + 8*l))/2) + 1

    M = similar(v, D, D)
    fill!(M, 0)
    #M = zeros(eltype(v), D, D)
    offset = 0
    for i in 1:D-1
        M[diagind(M, -i)] = v[offset+1:offset+D-i]
        offset += D-i
    end
    M

    #LowerTriangular(M)
end

"""
    matrix(diagM, vec_trilM)

Returns a symmetrix matrix from the diagonal and the "tril" form of a
matrix.

See also [`vec_tril`](@ref), [`inv_vec_tril`](@ref)
"""
function matrix(diagM, vec_trilM)
    T = eltype(diagM)
    trilM = inv_vec_tril(vec_trilM)
    retval = trilM + trilM'
    retval[diagind(retval)] .= diagM
    retval
end

max_stars_repo_name: lucasondel/ExpFamilyDistributions-jl

max_issues_repo_name: BUTSpeechFIT/ExpFamilyDistributions

max_forks_repo_name: lucasondel/ExpFamilyDistributions


text:
"""
Concatenated-power-mean-embeddings (CPMean) structure for
document embedding using word vectors.

# References
  * [Rücklé et al. 2018 "Concatenated power mean word embeddings
     as universal cross-lingual sentence representations"](https://arxiv.org/abs/1803.01400a)
"""
struct CPMeanEmbedder{T,S} <: WordVectorsEmbedder{T,S}
    embeddings::EmbeddingsLibrary{T,S}
    powers::Vector{T}
    znorm::Bool
    config::NamedTuple
end

function CPMeanEmbedder(embeddings::EmbeddingsLibrary{T,S},
                        config;
                        powers::Vector{T}=T[-Inf, 0.0, 1.0, Inf],
                        znorm::Bool=true,
                        kwargs...
                      ) where {T<:AbstractFloat, S<:AbstractString}
    return CPMeanEmbedder(embeddings, powers, znorm, config)
end

# Dimensionality function
function dimensionality(embedder::CPMeanEmbedder)
    # Embedding corresponding to all powers are
    # concatenated vertically
    return length(embedder.powers) * size(embedder.embeddings)[1]
end

# Sentence embedding function - returns a `embedder.dim`×1 matrix
function sentences2vec(embedder::CPMeanEmbedder,
                       document_embedding::Vector{Matrix{T}};
                       kwargs...) where {T,S}
    if isempty(document_embedding)
        return zeros(T, dimensionality(embedder), 0)
    else
        return __cpmean(document_embedding,
                        embedder.powers,
                        embedder.znorm,
                        dimensionality(embedder))
    end
end

function __cpmean(document_embedding::Vector{Matrix{T}},
                  powers::Vector{T},
                  znorm::Bool,
                  dim::Int
                 ) where {T<:AbstractFloat}
    embs = hcat(document_embedding...)
    n = size(embs, 2)  # total number of embedded words in all sentences
    m = size(embs, 1)  # embedding dimensionality
    X = zeros(T, dim, 1)

    # Construct various power mean functions
    # (they all work on matrices where the columns
    # are embeddings)
    dict_f = Dict{T, Function}(-Inf => (A,p)->minimum(A, dims=2),
                               Inf => (A,p)->maximum(A, dims=2),
                               0.0 => (A,p)->begin
                                            p = prod(A, dims=2);
                                            sign.(p).*abs.(p).^(1/n);
                                        end
                              )
    default_f = (A,p)->(1/n .* sum(A.^p, dims=2)).^(1/p)

    # Build document embedding
    i = 1
    @inbounds @simd for p in powers
        f = get(dict_f, p, default_f)
        Xp = view(X, (i-1)*m+1 : i*m, 1:1)
        Xp .= f(embs, p)
        znorm && znormalize!(Xp)
        i+=1
    end
    return X
end

function znormalize!(X::AbstractMatrix{T}) where {T}
    μ, σ = mean(X), std(X)
    X .= (X.-μ)./(σ+eps(T))
end


text:
module TestDimensionalityReductionWrapper

include(joinpath("..", "fixture_learners.jl"))
using .FixtureLearners
fcp = FeatureClassification()

using FactCheck

importall Orchestra.Transformers.DimensionalityReductionWrapper

facts("DimensionalityReduction transformers") do
  context("PCA transforms features") do
    instances = [
      5 10;
      -5 0;
      0 5;
    ]
    labels = ["x"; "y"; "z"]
    options = {:center => false, :scale => false}
    pca = PCA(options)
    fit!(pca, instances, labels)
    transformed = transform!(pca, instances)

    @fact true => maximum(instances - transformed * pca.model.rotation') < 10e-4
  end
end

end # module


text:
is_optional(member::SpecStructMember) = member.name == :pNext || member.requirement ∈ [OPTIONAL, POINTER_OPTIONAL] || is_inferable_length(member)
is_optional(param::SpecFuncParam) = param.requirement ∈ [OPTIONAL, POINTER_OPTIONAL]

"""
Represent an integer that gives the start of a C pointer.
"""
function is_pointer_start(spec::Spec)
    params = children(parent_spec(spec))
    any(params) do param
        !isempty(param.arglen) &&
            spec.type == :UInt32 &&
            string(spec.name) == string("first", uppercasefirst(replace(string(param.name), r"Count$" => "")))
    end
end

is_semantic_ptr(type) = is_ptr(type) || type == :Cstring
needs_deps(spec::SpecStruct) = any(is_semantic_ptr, spec.members.type)
must_return_success_code(spec::SpecFunc) = length(spec.success_codes) > 1 && :VK_INCOMPLETE ∉ spec.success_codes
must_repeat_while_incomplete(spec::SpecFunc) = !must_return_success_code(spec) && :VK_INCOMPLETE ∈ spec.success_codes
is_data_with_retrievable_size(spec::SpecFuncParam) = is_data(spec) && len(spec).requirement == POINTER_REQUIRED
is_opaque_data(spec) = is_data(spec) && len(spec).requirement ≠ POINTER_REQUIRED
is_opaque_pointer(type) = is_ptr(type) && is_void(ptr_type(type))
is_opaque_pointer(spec::Spec) = is_opaque_pointer(spec.type)
is_opaque(spec) = is_opaque_data(spec) || is_opaque_pointer(spec)
is_implicit_return(spec::SpecFuncParam) =
    !is_opaque_data(spec) &&
    !spec.is_constant &&
    is_ptr(spec.type) &&
    !is_length(spec) &&
    spec.type ∉ extension_types &&
    ptr_type(spec.type) ∉ extension_types
has_implicit_return_parameters(spec::SpecFunc) = any(is_implicit_return, children(spec))
is_flag(type) = type in spec_flags.name
is_flag(spec::Union{SpecFuncParam,SpecStructMember}) = spec.type in spec_flags.name
is_flag_bitmask(type) = type ∈ getproperty.(filter(!isnothing, spec_flags.bitmask), :name)
is_fn_ptr(type) = startswith(string(type), "PFN_")
is_fn_ptr(spec::Spec) = is_fn_ptr(spec.type)

function is_hl(type)
    vktype = Symbol(:Vk, type)
    vktype in [spec_structs.name; spec_unions.name]
end

is_intermediate(type) = startswith(string(type), '_')

has_intermediate_type(::SpecHandle) = false
has_intermediate_type(spec::Union{SpecStruct,SpecUnion}) = true

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "STL"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/eschnett/STL.jl", "76f5f5eeaa0b789010b302ca5cf9d8551e3c83fb"),
    GitSource("https://github.com/eschnett/TestAbstractTypes.jl", "a23107bf47796db9d414c77801c6b3331f4950f0"),
    ArchiveSource("https://julialang-s3.julialang.org/bin/musl/x64/1.6/julia-1.6.5-musl-x86_64.tar.gz",
                  "e38eece6f9f20c7472caf3f8f74a99ad0880921c28e1301461fa7af919880383"),
]

# Bash recipe for building across all platforms
script = raw"""
export PATH="${PATH}:${WORKSPACE}/srcdir/julia-1.6.5/bin"

# Switch into source directory
cd STL.jl

# Create C++ wrapper code
julia --project=@. --eval '
    using Pkg

    oldpwd = pwd()
    cd("..")
    Pkg.generate("STL_jll")
    write("STL_jll/src/STL_jll.jl",
        \"\"\"
        module STL_jll
        # Just a placeholder library name
        const libSTL = "libSTL"
        export libSTL
        end
        \"\"\")
    cd(oldpwd)

    # Replace `STL_jll` which may not yet exist
    Pkg.rm("STL_jll")
    # Replace `TestAbstractTypes` which is not registered yet
    Pkg.rm("TestAbstractTypes")

    Pkg.develop(path="../STL_jll")
    Pkg.develop(path="../TestAbstractTypes.jl")

    # Pkg.build("STL")
    using STL
    STL.cxx_write_code!()
'

# Build and install C++ wrapper code
${CXX} -Drestrict=__restrict__ -std=c++17 -g -O2 -fPIC -shared -o ${libdir}/libSTL.${dlext} \
    StdMap.cxx StdSharedPtr.cxx StdString.cxx StdVector.cxx
install_license LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libSTL", :libSTL),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version = v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#=
Script to read in a circuit from a qasm file, convert it to a tensor network
and contract it.
=#

using PicoQuant
using ArgParse
using JSON

"""
    parse_commandline()

Parse command line options and return argument dictionary
"""
function parse_commandline(ARGS)
    s = ArgParseSettings()
        @add_arg_table! s begin
            "--qasm", "-q"
                help = "QASM file to load"
                required = true
            "--output", "-o"
                help = "Output file to write the contracted network structure to"
                arg_type = String
                required = true
            "--output_data", "-d"
                help = "Output file to write the contracted network data to"
                arg_type = String
                required = true
            "--decompose"
                help = "Decompose two qubit gates"
                action = :store_true
            "--transpile"
                help = "Transpile circuit so gates only operate between neighbouring gates"
                action = :store_true
            "--method", "-m"
                help = "The contraciton method to use (options, random, inorder, fullwf and mps, default fullwf)"
                arg_type = String
                default = "fullwf"
            "--indent"
                help = "Indent to use in output json file"
                action = :store_true
            "--input_config", "-c"
                help = "The starting configuration as a string of 0, 1, +, -. Default is all zeros"
                arg_type = String
                default = ""
            "--work_dir", "-w"
                help = "Folder to keep temporary files, if non given will use current folder"
                arg_type = String
                default = ""

        end
        return parse_args(ARGS, s)
end

function main(ARGS)
    parsed_args = parse_commandline(ARGS)

    work_dir = parsed_args["work_dir"] == "" ? "." : parsed_args["work_dir"]
    isdir(work_dir) || mkdir(work_dir)

    # prepend work files with work folder path
    qasm_filename = parsed_args["qasm"]
    filename_base = splitext(qasm_filename)[1]
    tl_filename = joinpath(work_dir, "$(filename_base).tl")
    data_filename = joinpath(work_dir, "$(filename_base).h5")

    # require a backend to save tensor data to
    circuit = load_qasm_as_circuit_from_file(qasm_filename)
    tng = convert_qiskit_circ_to_network(circuit, DSLBackend(tl_filename, data_filename),
                                         decompose=parsed_args["decompose"],
                                         transpile=parsed_args["transpile"])

    if parsed_args["input_config"] == ""
        add_input!(tng, "0"^tng.number_qubits)
    else
        @assert length(parsed_args["input_config"]) == tng.number_qubits
        add_input!(tng, parsed_args["input_config"])
    end

    method = parsed_args["method"]
    if method == "fullwf"
        full_wavefunction_contraction!(tng)
    elseif method == "mps"
        contract_mps_tensor_network_circuit!(tng)
    elseif method == "inorder"
        contract_network!(tng, inorder_contraction_plan(tng))
    elseif method == "random"
        contract_network!(tng, random_contraction_plan(tng))
    end

    execute_dsl_file(tl_filename, data_filename, parsed_args["output_data"])

    open(parsed_args["output"], "w") do io
        write(io, to_json(tng, parsed_args["indent"] ? 2 : 0))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end


text: Dict{String,Any}("with.dot" => Dict{String,Any}("value" => "2","type" => "integer"),"plain_table" => Dict{String,Any}("with.dot" => Dict{String,Any}("value" => "4","type" => "integer"),"plain" => Dict{String,Any}("value" => "3","type" => "integer")),"table" => Dict{String,Any}("withdot" => Dict{String,Any}("key.with.dots" => Dict{String,Any}("value" => "6","type" => "integer"),"plain" => Dict{String,Any}("value" => "5","type" => "integer"))),"plain" => Dict{String,Any}("value" => "1","type" => "integer"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
"""
	FormatHandler

Abstract type for format readers/writers

# Example

```julia
struct MyIO <: Formats.FormatHandler end
```
"""
abstract type FormatHandler end

struct FormatRegistry
	formats::Set{MIME}
	codings::Set{MIME}
	extensions::Dict{String,Vector{MIME}}
	signatures::Dict{Vector{UInt8},Vector{MIME}}
	readers::Dict{MIME,Vector{FormatHandler}}
	writers::Dict{MIME,Vector{FormatHandler}}
	favorite_readers::Dict{MIME,FormatHandler}
	favorite_writers::Dict{MIME,FormatHandler}
	global_favorite_readers::Vector{FormatHandler}
	global_favorite_writers::Vector{FormatHandler}
	decoders::Dict{MIME,Type}
	encoders::Dict{MIME,Type}
end

FormatRegistry() = FormatRegistry(
		Set{MIME}(),
		Set{MIME}(),
		Dict{String,Vector{MIME}}(),
		Dict{Vector{UInt8},Vector{MIME}}(),
		Dict{MIME,Vector{FormatHandler}}(),
		Dict{MIME,Vector{FormatHandler}}(),
		Dict{MIME,FormatHandler}(),
		Dict{MIME,FormatHandler}(),
		Vector{FormatHandler}(),
		Vector{FormatHandler}(),
		Dict{MIME,Type}(),
		Dict{MIME,Type}())

Base.copy(registry::FormatRegistry) = FormatRegistry(
		copy(registry.formats),
		copy(registry.codings),
		copy(registry.extensions),
		copy(registry.signatures),
		copy(registry.readers),
		copy(registry.writers),
		copy(registry.favorite_readers),
		copy(registry.favorite_writers),
		copy(registry.global_favorite_readers),
		copy(registry.global_favorite_writers),
		copy(registry.decoders),
		copy(registry.encoders))

function Base.empty!(registry::FormatRegistry)
	empty!(registry.formats)
	empty!(registry.codings)
	empty!(registry.extensions)
	empty!(registry.signatures)
	empty!(registry.readers)
	empty!(registry.writers)
	empty!(registry.favorite_readers)
	empty!(registry.favorite_writers)
	empty!(registry.global_favorite_readers)
	empty!(registry.global_favorite_writers)
	empty!(registry.decoders)
	empty!(registry.encoders)
end

function Base.merge!(registry::FormatRegistry, other::FormatRegistry)
	union!(registry.formats, other.formats)
	union!(registry.formats, other.codings)
	merge!(registry.extensions, other.extensions)
	merge!(registry.signatures, other.signatures)
	merge!(registry.readers, other.readers)
	merge!(registry.writers, other.writers)
	merge!(registry.favorite_readers, other.favorite_readers)
	merge!(registry.favorite_writers, other.favorite_writers)
	append!(registry.global_favorite_readers, other.global_favorite_readers)
	append!(registry.global_favorite_writers, other.global_favorite_writers)
	merge!(registry.decoders, other.decoders)
	merge!(registry.encoders, other.encoders)
end

const registry = FormatRegistry()

function newregistry(f::Function)
	tmp_registry = copy(registry)
	empty!(registry)
	f()
	merge!(registry, tmp_registry)
end

"""
	adddformat(name::AbstractString)

Register MIME type `name` as a data format.

# Example

```julia
Formats.addformat("image/png")
```
"""
function addformat(name::AbstractString)
	mime = MIME{Symbol(name)}()
	if mime in registry.codings
		error("$(name) is already registered as a data coding")
	else
		push!(registry.formats, mime)
	end
end

"""
	addcoding(name::AbstractString)

Register MIME type `name` as a data coding.

# Example

```julia
Formats.addcoding("application/gzip")
```
"""
function addcoding(name::AbstractString)
	mime = MIME{Symbol(name)}()
	if mime in registry.formats
		error("$(name) is already registered as a data format")
	else
		push!(registry.codings, mime)
	end
end

"""
	addextension(format::AbstractString, extension::AbstractString)

Register a filename `extension` for a specific `format`.

# Example

```julia
Formats.addextension("image/png", ".png")
```
"""
function addextension(format::AbstractString, extension::AbstractString)
	mime = MIME{Symbol(format)}()
	startswith(extension, ".") || error("invalid filename extension")
	extension = lowercase(extension)
	mimes = get!(registry.extensions, extension, [])
	if mime in mimes
		@debug("extension $(extension) already registered for $(format)")
	else
		push!(mimes, mime)
	end
	if length(mimes) > 1
		@warn("extension $(extension) registered for multiple formats")
	end
end

"""
	addsignature(format::AbstractString, signature::Vector{UInt8})
	addsignature(format::AbstractString, signature::AbstractString)

Register a (max 512-bytes) `signature` for a specific `format`.

# Example

```julia
Formats.addsignature("image/png", [0x89,0x50,0x4e,0x47,0x0d,0x0a,0x1a,0x0a])
```
"""
function addsignature(format::AbstractString, signature::Vector{UInt8})
	mime = MIME{Symbol(format)}()
	length(signature) > 0 || error("signature cannot be empty")
	length(signature) > 512 && error("maximum signature size is 512 bytes")
	mimes = get!(registry.signatures, signature, [])
	if mime in mimes
		@debug("signature already registered for $(format)")
	else
		push!(mimes, mime)
	end
	if length(mimes) > 1
		@warn("signature registered for multiple formats")
	end
end

addsignature(format::AbstractString, signature::AbstractString) =
		addsignature(format, Vector{UInt8}(signature))

"""
	addreader(format::AbstractString, reader::FormatHandler)

Register a `reader` for the specified `format`.

# Example

```julia
struct MyIO <: Formats.FormatHandler end

Formats.addreader("image/png", MyIO())
```
"""
function addreader(format::AbstractString, reader::FormatHandler)
	mime = MIME{Symbol(format)}()
	readers = get!(registry.readers, mime, [])
	if reader in readers
		error("reader $(reader) already registered for $(format)")
	else
		push!(readers, reader)
	end
	if length(readers) > 1
		@info("$(format) has multiple registered readers")
	end
end

"""
	addwriter(format::AbstractString, writer::FormatHandler)

Register a `writer` for the specified `format`.

# Example

```julia
struct MyIO <: Formats.FormatHandler end

Formats.addwriter("image/png", MyIO())
```
"""
function addwriter(format::AbstractString, writer::FormatHandler)
	mime = MIME{Symbol(format)}()
	writers = get!(registry.writers, mime, [])
	if writer in writers
		error("writer $(writer) already registered for $(format)")
	else
		push!(writers, writer)
	end
	if length(writers) > 1
		@info("$(format) has multiple registered writers")
	end
end

"""
	preferreader(reader::FormatHandler, [format::AbstractString])

Register `reader` as the preferred alternative for the specified `format`.

If `format` is omitted, `reader` becomes a globally preferred alternative: it
is used for all the formats it can read.
```
"""
function preferreader(reader::FormatHandler)
	if reader in registry.global_favorite_readers
		error("reader $(reader) is already globally preferred")
	end
	push!(registry.global_favorite_readers, reader)
end

function preferreader(reader::FormatHandler, format::AbstractString)
	mime = MIME{Symbol(format)}()
	readers = get(registry.readers, mime, [])
	if ! (reader in readers)
		error("reader $(reader) is not registered for $(format)")
	end
	if haskey(registry.favorite_readers, mime)
		@warn("replacing preferred reader for $(format)")
	end
	registry.favorite_readers[mime] = reader
end

"""
	preferwriter(writer::FormatHandler, [format::AbstractString])

Register `writer` as the preferred alternative for the specified `format`.

If `format` is omitted, `writer` becomes a globally preferred alternative: it
is used for all the formats it can write.
"""
function preferwriter(writer::FormatHandler)
	if writer in registry.global_favorite_writers
		error("writer $(writertype) is already globally preferred")
	end
	push!(registry.global_favorite_writers, writer)
end

function preferwriter(writer::FormatHandler, format::AbstractString)
	mime = MIME{Symbol(format)}()
	writers = get(registry.writers, mime, [])
	if ! (writer in writers)
		error("writer $(writer) is not registered for $(format)")
	end
	if haskey(registry.favorite_writers, mime)
		@warn("replacing preferred writer for $(format)")
	end
	registry.favorite_writers[mime] = writer
end

"""
	setdecoder(coding::AbstractString, decoder::Type)

Set the `decoder` to use for the specified `coding`.
"""
function setdecoder(coding::AbstractString, decoder::Type)
	mime = MIME{Symbol(coding)}()
	if haskey(registry.decoders, mime)
		@warn("replacing decoder for $(coding)")
	end
	registry.decoders[mime] = decoder
end

"""
	setencoder(coding::AbstractString, encoder::Type)

Set the `encoder` to use for the specified `coding`.
"""
function setencoder(coding::AbstractString, encoder::Type)
	mime = MIME{Symbol(coding)}()
	if haskey(registry.encoders, mime)
		@warn("replacing encoder for $(coding)")
	end
	registry.encoders[mime] = encoder
end

function resolvehandler(mime::MIME, thing::AbstractString,
		available::Dict{MIME,Vector{FormatHandler}},
		favorites::Dict{MIME,FormatHandler},
		global_favorites::Vector{FormatHandler})
	name = string(mime)
	handlers = get(available, mime) do
		error("no $(thing) registered for $(name)")
	end
	nfavorites = 0
	for favorite in global_favorites
		if favorite in handlers
			handler = favorite
			nfavorites += 1
		end
	end
	if nfavorites > 1
		error("multiple $(thing)s for $(name) are marked as favorites")
	elseif nfavorites == 0
		handler = get(favorites, mime, nothing)
		if handler == nothing
			handler = handlers[1]
			if length(handlers) > 1
				@warn("multiple applicable $(thing)s, using $(handler)")
			end
		else
			if ! (handler in handlers)
				error("$(thing) $(handler), preferred for $(name), " *
						"is not registered")
			end
		end
	else
		if haskey(favorites, mime)
			@warn("both a global and a format-specific preferred $(thing) " *
					"are registered for $(name); " *
					"using global preferred $(thing) $(handler)")
		end
	end
	handler
end

resolvereader(mime::MIME) =
		resolvehandler(mime, "reader", registry.readers,
		registry.favorite_readers, registry.global_favorite_readers)

resolvewriter(mime::MIME) =
		resolvehandler(mime, "writer", registry.writers,
		registry.favorite_writers, registry.global_favorite_writers)

function resolvedecoder(mime::MIME)
	get(registry.decoders, mime) do
		error("no decoder registered for format $(string(mime))")
	end
end

function resolveencoder(mime::MIME)
	get(registry.encoders, mime) do
		error("no encoder registered for format $(string(mime))")
	end
end


text:
using AlphaGo
using AlphaGo: go
using Flux
using BSON: @load

# Note:
# do Pkg.checkout("WebIO", "s/asset-registry")
# and Pkg.checkout("Mux", "s/asset-registry")
include("../src/interface.jl")
include("./helpers.jl")

set_all_params(9)

assets = "$(@__DIR__)/assets"
files=[
    "$(@__DIR__)/../lib/interface.js"
    "$(@__DIR__)/../lib/interface.css"
    "$(assets)/css/go.css"
    "$(assets)/css/wgo.player.css"
    "$(assets)/js/wgo.min.js"
    "$(assets)/js/go.js"]

gw = GameWindow(files; config=Dict("boardSize"=>9))
setup(gw) # default setup

# agz_nn = loadNeuralNet()
agz_nn = NeuralNet(;tower_height=1)
agz = MCTSPlayer(agz_nn, num_readouts = 64, two_player_mode = true)

playing = false
isPlaying() = playing
setPlaying(p=true) = (playing = p)

# game play
oninput(gw) do x
    println("js sent $x")
    is_done(agz) && ( return gw.msg[] = "Game over! " * agz.result_string )
    isPlaying() && return gw.msg[] = "computer's turn!"
    isIllegalMove(x) && return gw.msg[] = "Illegal Move"
    setPlaying()

    move = playMove(x) # user action
    println("user action...")
    updateBoard(x)   # this happens kinda late ( along with the next update )
    move = playMove()  # computer action
    println("computer action...")
    updateBoard(stringToAction(go.to_kgs(move), -1))

    setPlaying(false)
    is_done(agz) && ( return gw.msg[] = "Game over! " * agz.result_string )
end

function startGame()
    initialize_game!(agz)
    gw(dom"div#demo_wrapper"(
        dom"div#controls"(
            dom"h1"("Go game"),
            dom"button.pass"("Pass")
        ),
        dom"div#playground"()))
end

startGame()


text:
module Linda

using Printf

import MathOptInterface
const MOI = MathOptInterface

# package code goes here
include("env.jl")
include("column.jl")
include("oracle.jl")
include("master.jl")
include("colgen.jl")

end # module


text:
#!/usr/bin/env julia

"""
A colorscheme is an array of colors. To use the package:

```
using ColorSchemes, Colors
```

and, if you want image input and output:

```
using Images, FileIO
```

The names of the registered (built-in) colorschemes are listed in the `schemes` array.

To use one of the built-in colorschemes, use the symbol:

```
julia> ColorSchemes.picasso
```

or import it:

```
julia> import ColorSchemes.pigeon
julia> pigeon
```

## Functions:

```
extract(imfile, n=10, i=10, tolerance=0.01; kwargs...)
```
    - extract a new colorscheme from an image file, return a colorscheme

```
extract_weighted_colors(imfile, n=10, i=10, tolerance=0.01; shrink = 2.0)
```
    - return a colorscheme and weights for each entry

```
get(cscheme::Vector{C}, x) where {C <: Colorant}
```
    - return a single color from cscheme based on location of x (0 - 1) in cscheme

```
get(cscheme::Vector{C}, inData :: Array{Number, 2}, rangescale) where {C <: Colorant}
```
    - return an RGB image generated by applying the color scheme to the 2D input data.

The mapping between values and ColorScheme can be adjusted with rangescale, see help on `get()`.

```
colorscheme_weighted(cscheme::Vector{C}, weights, l = 50) where {C <: Colorant}
```
    - return a weighted colorscheme, given a colorscheme and an array of weights for each entry

```
colorscheme_to_image(cs::Vector{C}, nrows=50, tilewidth=5) where {C <: Colorant}
```
    - make an image of a scheme by repeating each color m times in h rows

```
compare_colors(color_a, color_b, field = :l)
```
    - compare colors, return true if the specified field of `color_a` is less than `color_b`.

```
colorscheme_to_text(cs::Vector{C}, schemename::String, file::String; comment="") where {C <: Colorant}
```
    - export a colorscheme to a text file

```
image_to_swatch(imagefilepath, n::Int64, destinationpath; nrows=50, tilewidth=5)
```
    - extract a colorscheme and save it as a swatch PNG in the destination file

```
schemes
```
    - an array of names of all the loaded schemes, as symbols

```
sortcolorscheme(colorscheme::Vector{C}, field = :l; kwargs...) where {C <: Colorant}
```
    - sort a colorscheme

"""
module ColorSchemes

using Images, Colors, Clustering, FileIO, Dates

import Base.get

export
colorscheme_to_image,
colorscheme_to_text,
colorscheme_weighted,
extract,
extract_weighted_colors,
get,
getinverse,
convert_to_scheme,
image_to_swatch,
schemes,
sortcolorscheme,
@reg

"""
`schemes` contains the names of all the available colorschemes.
"""
const schemes = Symbol[]

"""
    reg(vname, args)

Load a variable and some values, and add the symbol to the list of schemes in
`schemes`.
"""
macro reg(vname, args)
    quote
        $(esc(push!(schemes, vname)))
        $(esc(vname)) = $(args)
    end
end

# load the installed schemes
include(dirname(@__FILE__) * "/../data/allcolorschemes.jl")
include(dirname(@__FILE__) * "/../data/colorbrewerschemes.jl")
include(dirname(@__FILE__) * "/../data/matplotlib.jl")
include(dirname(@__FILE__) * "/../data/cmocean.jl")

# the `schemes` array now contains the names of the built-in ColorSchemes

# convert a value between oldmin/oldmax to equivalent value between newmin/newmax

remap(value, oldmin, oldmax, newmin, newmax) =
    ((value .- oldmin) ./ (oldmax .- oldmin)) .* (newmax .- newmin) .+ newmin

"""
    extract(imfile, n=10, i=10, tolerance=0.01; shrink=n)

`extract()` extracts the most common colors from an image from the image file `imfile`
by finding `n` dominant colors, using `i` iterations. You can (and probably should)
shrink larger images before running this function.

Returns a colorscheme (an array of colors)
"""
function extract(imfile, n=10, i=10, tolerance=0.01; kwargs...)
    return extract_weighted_colors(imfile, n, i, tolerance; kwargs...)[1] # throw away the weights
end

"""
    extract_weighted_colors(imfile, n=10, i=10, tolerance=0.01; shrink = 2)

Extract colors and weights of the clusters of colors in an image file.

Example:

    pal, wts = extract_weighted_colors(imfile, n, i, tolerance; shrink = 2)
"""
function extract_weighted_colors(imfile, n=10, i=10, tolerance=0.01; shrink = 2.0)
    img = load(imfile)
    # TODO this is wrong, isn't it?
    (!@isdefined img) && error("Can't load the image file \"$imfile\"")
    w, h = size(img)
    neww = round(Int, w/shrink)
    newh = round(Int, h/shrink)
    smaller_image = Images.imresize(img, (neww, newh))
    w, h = size(smaller_image)
    imdata = convert(Array{Float64}, channelview(smaller_image))
    d = reshape(imdata, 3, :)
    R = kmeans(d, n, maxiter=i, tol=tolerance)
    colscheme = RGB{Float64}[]
    for i in 1:3:length(R.centers)
        push!(colscheme, RGB(R.centers[i], R.centers[i+1], R.centers[i+2]))
    end
    return colscheme, R.cweights/sum(R.cweights)
end

"""
    colorscheme_weighted(colorscheme, weights, length)

Returns a new colorscheme of length `length` (default 50) where the proportion
of each color in `colorscheme` is represented by the associated weight of each entry.

Examples:

```
colorscheme_weighted(extract_weighted_colors("hokusai.jpg")...)
colorscheme_weighted(extract_weighted_colors("filename00000001.jpg")..., 500)
```
"""
function colorscheme_weighted(cscheme::Vector{C}, weights, l = 50) where {C <: Colorant}
    iweights = map(n -> convert(Integer, round(n * l)), weights)
    #   adjust highest or lowest so that length of result is exact
    while sum(iweights) < l
        val, ix = findmin(iweights)
        iweights[ix]=val+1
    end
    while sum(iweights) > l
        val,ix = findmax(iweights)
        iweights[ix]=val-1
    end
    a = Array{RGB{Float64}}(undef, 0)
    for n in 1:length(cscheme)
        a = vcat(a, repeat([cscheme[n]], iweights[n]))
    end
    return a
end

"""
    compare_colors(color_a, color_b, field = :l)

Compare two colors, using the Luv colorspace. `field` defaults to luminance `:l`
but could be `:u` or `:v`. Return true if the specified field of `color_a` is
less than `color_b`.
"""
function compare_colors(color_a, color_b, field = :l)
    if 1 < color_a.r < 255
        fac = 255
    else
        fac = 1
    end
    luv1 = convert(Luv, RGB(color_a.r/fac, color_a.g/fac, color_a.b/fac))
    luv2 = convert(Luv, RGB(color_b.r/fac, color_b.g/fac, color_b.b/fac))
    return getfield(luv1, field) < getfield(luv2, field)
end

"""
    sortcolorscheme(colorscheme, field; kwargs...)

Sort (non-destructively) a colorscheme using a field of the LUV colorspace.

The less than function is `lt = (x,y) -> compare_colors(x, y, field)`.

The default is to sort by the luminance field `:l` but could be by `:u` or `:v`.
"""
function sortcolorscheme(colorscheme::Vector{C}, field = :l; kwargs...) where {C <: Colorant}
    sort(colorscheme, lt = (x,y) -> compare_colors(x, y, field); kwargs...)
end

"""
    get(cscheme, x)

Find the nearest color in a colorscheme `cscheme` corresponding to a point `x` between 0 and 1.

Returns a single color.
"""
function get(cscheme::Vector{C}, x, rangescale) where {C<:Colorant}
    if rangescale==:clamp
        get(cscheme, x, (0.0, 1.0))
    elseif (rangescale==:extrema)
        get(cscheme, x, extrema(x))
    else
        error("rangescale ($rangescale) not supported, should be :clamp, :extrema or tuple (minVal, maxVal)")
    end
end

"""

    get(cscheme, inData :: Array{Number, 2}, rangescale=:clamp)
    get(cscheme, inData :: Array{Number, 2}, rangescale=(minVal, maxVal))

Return an RGB image generated by applying the colorscheme to the 2D input data.

If `rangescale` is `:clamp` the colorscheme is applied to values between 0.0-1.0, and values
outside this range get clamped to the ends of the colorscheme.

Else, if `rangescale` is `:extrema`, the colorscheme is applied to the range `minimum(indata)..maximum(indata)`.

# Examples

```
img = get(ColorSchemes.leonardo, rand(10,10))
save("testoutput.png", img)  # might need FileIO or similar

img2 = get(ColorSchemes.leonardo, 10.0*rand(10,10), :extrema)
img3 = get(ColorSchemes.leonardo, 10.0*rand(10,10), (1.0, 9.0))

# Also works with PerceptualColourMaps
using PerceptualColourMaps
img4 = get(PerceptualColourMaps.cmap("R1"), rand(10,10))
```
"""
function get(cscheme::Vector{C}, x, rangescale :: Tuple{Number, Number}=(0.0, 1.0)) where {C<:Colorant}
    x = clamp.(x, rangescale...)
    before_fp = remap(x, rangescale..., 1, length(cscheme))
    before = round.(Int, before_fp, RoundDown)
    after =  min.(before .+ 1, length(cscheme))
    # blend between the two colors adjacent to the point
    cpt = before_fp - before
    return weighted_color_mean.(1 .- cpt, cscheme[before], cscheme[after])
end

"""
    getinverse(cscheme, c, range=(0.0, 1.0))

Computes where the provided Color `c` would fit in `cscheme`.

This is the inverse of `get()` — it returns the value `x` in the provided `range`
for which `get(scheme, x)` would most closely match the provided Color `c`.

# Examples
```jldoctest
julia> getinverse(ColorSchemes.leonardo, RGB(1,0,0))
0.6248997995654847
julia> getinverse([RGB(0,0,0), RGB(1,1,1)], RGB(.5,.5,.5))
0.5432555858022048
julia> cs = range(RGB(0,0,0), stop=RGB(1,1,1), length=5);
julia> getinverse(cs, cs[3])
0.5
```
"""
function getinverse(cscheme::Vector{C}, c, rangescale :: Tuple{Number, Number}=(0.0, 1.0)) where {C<:Colorant}
    if length(cscheme) <= 1 ; throw(InexactError()) ; end
    cdiffs = map(c_i->colordiff(promote(c,c_i)...), cscheme)
    closest = argmin(cdiffs)
    left = right = 0;
    if closest == 1 ; left = closest; right = closest + 1;
    elseif closest == length(cscheme) ; left = closest - 1; right = closest;
    else
        next_closest = cdiffs[closest-1] < cdiffs[closest+1] ? closest-1 : closest+1
        left = min(closest, next_closest)
        right = max(closest, next_closest)
    end

    v = left
    if cdiffs[left] != cdiffs[right] ;  # Prevents divide by 0.
        v += ( cdiffs[left] / (cdiffs[left] + cdiffs[right]))
    end
    return ColorSchemes.remap(v, 1, length(cscheme), rangescale...)
end

"""
    convert_to_scheme(cscheme, img)

Converts img from its current color values to use only the colors defined in cscheme.

```
image = nonTransparentImg
convert_to_scheme(ColorSchemes.leonardo, image)
convert_to_scheme(ColorSchemes.Paired_12, image)
```
"""
convert_to_scheme(cscheme::Vector{C},img) where {C<:Colorant} =
    map(c->get(cscheme, getinverse(cscheme, c)), img)

"""
    colorscheme_to_text(cscheme, schemename, filename; comment="")

Write a colorscheme to a Julia file in a format suitable for `include`ing.

Example:

```
colorscheme_to_text(
    extract("/tmp/1920px-Great_Wave_off_Kanagawa2.jpg"),
        "hokusai_1",
        "/tmp/hok.jl",
        comment="from Hokusai's Great Wave")
```

To read a text file created thusly in and register it in `schemes`:

```
julia> include("/tmp/hok.jl")
julia> schemes[end]
:hokusai_1
julia> get(hokusai_1, .4)
RGB{Float64}(0.5787354153400166,0.49341844091747,0.22277034922842723)
```
"""
function colorscheme_to_text(cs::Vector{C}, schemename::String, file::String; comment="") where {C <: Colorant}
    fhandle = open(file, "w")
    write(fhandle, string("# ", comment, "\n"))
    write(fhandle, string("# created $(Dates.now())\n"))
    write(fhandle, string("@reg $schemename [\n", join(cs, ",\n"), " ]"))
    close(fhandle)
end

"""
    colorscheme_to_image(cs, nrows=50, tilewidth=5)

Make an image from a colorscheme by repeating the colors in a colorscheme.

Returns the image as an array.

Examples:

```
using FileIO

img = colorscheme_to_image(ColorSchemes.leonardo, 50, 200)
save("/tmp/cs_image.png", img)

save("/tmp/blackbody.png", colorscheme_to_image(ColorSchemes.blackbody, 10, 100))
```
"""
function colorscheme_to_image(cs::Vector{C}, nrows=50, tilewidth=5) where {C <: Colorant}
    ncols = tilewidth * length(cs)
    a = Array{RGB{Float64}}(undef, nrows, ncols)
    for row in 1:nrows
        for col in 1:ncols
            a[row, col] = cs[div(col-1, tilewidth) + 1]
        end
    end
    return a
end

"""
    image_to_swatch(imagefilepath, samples, destinationpath; nrows=50, tilewidth=5)

Extract a colorscheme from the image in `imagefilepath` to a swatch image PNG in
`destinationpath`. This just runs `sortcolorscheme()`, `colorscheme_to_image()`, and
`save()` in sequence.

Specify the number of colors. You can also specify the number of rows, and how many
times each color is repeated.

```
image_to_swatch("monalisa.jpg", 10, "/tmp/monalisaswatch.png")
```
"""
function image_to_swatch(imagefilepath, n::Int64, destinationpath; nrows=50, tilewidth=5)
    temp = sortcolorscheme(extract(imagefilepath, n))
    img = colorscheme_to_image(temp, nrows, tilewidth)
    save(destinationpath, img)
end

end


text:
using MLStyle
using Random: AbstractRNG
using KeywordCalls
using ConstructionBase
export @parameterized, @half

# A fold over ASTs. Example usage in `replace`
function foldast(leaf, branch; kwargs...)
    function go(ast)
        MLStyle.@match ast begin
            Expr(head, args...) => branch(head, map(go, args); kwargs...)
            x                   => leaf(x; kwargs...)
        end
    end

    return go
end

function replace(p, f, expr)
    leaf(s) = p(s) ? f(s) : s
    branch(head, newargs) = Expr(head, newargs...)
    foldast(leaf, branch)(expr)
end

# from https://thautwarm.github.io/MLStyle.jl/latest/tutorials/capture.html
function capture(template, ex, action)
    let template = Expr(:quote, template)
        quote
            @match $ex begin
                $template => $action
                _         => nothing
            end
        end
    end
end

macro capture(template, ex, action)
    capture(template, ex, action) |> esc
end

function _parameterized(__module__, expr)
    @capture ($op($μ($(p...)), $base)) expr begin
        @assert op ∈ [:<<, :≪, :≃]
        μbase = Symbol(:__, μ, :_base)

        μ = esc(μ)
        base = esc(base)

        # @gensym basename
        q = quote
            struct $μ{N,T} <: MeasureBase.ParameterizedMeasure{N}
                par::NamedTuple{N,T}
            end

            const $μbase = $base
            MeasureBase.basemeasure(::$μ) = $μbase
        end

        if !isempty(p)
            # e.g. Normal(μ,σ) = Normal((μ=μ, σ=σ))
            pnames = QuoteNode.(p)
            push!(q.args, :($μ($(p...)) = $μ(NamedTuple{($(pnames...),)}(($(p...),)))))
        end

        return q
    end

    @capture $μ($(p...)) expr begin
        μ = esc(μ)

        q = quote
            struct $μ{N,T} <: MeasureBase.ParameterizedMeasure{N}
                par::NamedTuple{N,T}
            end
        end

        if !isempty(p)
            # e.g. Normal(μ,σ) = Normal((μ=μ, σ=σ))
            pnames = QuoteNode.(p)
            push!(q.args, :($μ($(p...)) = $μ(NamedTuple{($(pnames...),)}(($(p...),)))))
        end

        return q
    end

    @capture $μ($(p...)) expr begin
        μ = esc(μ)

        q = quote
            struct $μ{N,T} <: MeasureBase.ParameterizedMeasure{N}
                par::NamedTuple{N,T}
            end
        end

        if !isempty(p)
            # e.g. Normal(μ,σ) = Normal((μ=μ, σ=σ))
            pnames = QuoteNode.(p)
            push!(q.args, :($μ($(p...)) = $μ(NamedTuple{($(pnames...),)}(($(p...),)))))
        end

        return q
    end
end

"""
    @parameterized <declaration>

The <declaration> gives a measure and its default parameters, and specifies
its relation to its base measure. For example,

    @parameterized Normal(μ,σ)

declares the `Normal` is a measure with default parameters `μ and σ`. The result is equivalent to
```
struct Normal{N,T} <: ParameterizedMeasure{N}
    par :: NamedTuple{N,T}
end

KeywordCalls.@kwstruct Normal(μ,σ)

Normal(μ,σ) = Normal((μ=μ, σ=σ))
```

See [KeywordCalls.jl](https://github.com/cscherrer/KeywordCalls.jl) for details on `@kwstruct`.
"""
macro parameterized(expr)
    _parameterized(__module__, expr)
end

"""
    @half dist([paramnames])

Starting from a symmetric univariate measure `dist ≪ Lebesgue(ℝ)`, create a new
measure `Halfdist ≪ Lebesgue(ℝ₊)`. For example,

    @half Normal()

creates `HalfNormal()`, and

    @half StudentT(ν)

creates `HalfStudentT(ν)`.
"""
macro half(ex)
    _half(__module__, ex)
end

function _half(__module__, ex)
    @match ex begin
        :($dist) => begin
            halfdist = esc(Symbol(:Half, dist))
            dist = esc(dist)
            quote
                $halfdist(args...) = half($dist(args...))
                $halfdist(; kwargs...) = half($dist(; kwargs...))
            end
        end
    end
end


text:
"""
SimulationModule (C) 2020, Martin Gurtner
Simulation of the dynamics for the dielectrophoretic, magnetophoretic, and acoustophoretic distributed manipulation platforms
"""
module SimulationModule

using ..ActuatorArrayModule
using ..ObjectAgentModule

using DifferentialEquations
using LinearAlgebra

export AgentCtrlDEP, AgentCtrlMAG, CircularTrajectoryGenerator, WaypointsTrajectoryGenerator, RegulatorP, control_law, generate_reference, simulate!, agent_pos

## Abstract types
abstract type AbstractTrajectoryGenerator end
abstract type AbstractRegulator end
abstract type AbstractAgentCtrl end

## Controllers
struct RegulatorP <: AbstractRegulator
    P::Real
    F_sat::Real
end

## Trajectory generators
struct CircularTrajectoryGenerator <: AbstractTrajectoryGenerator
    radius::Real
    freq::Real
    phase::Real
    offset::Tuple{T,T}  where T<:Real
end

mutable struct WaypointsTrajectoryGenerator <: AbstractTrajectoryGenerator
    points::Union{Array{Tuple{T, T}, 1}, Array{Tuple{T, T, T}, 1}} where T<:Real
    switch_dist::Real
    current_point_index
    repeat
    dir

    function WaypointsTrajectoryGenerator(points::Array{Tuple{T, T}, 1}, switch_dist, repeat=false) where T<:Real
        new(points, switch_dist, 1, repeat, 1)
    end
end

## Agents - structures storing the state of the manipulated objects (position and velocity)
mutable struct AgentCtrlMAG <: AbstractAgentCtrl
    state::Vector{Float64}
    ref_gen::AbstractTrajectoryGenerator
    reg::AbstractRegulator
end

mutable struct AgentCtrlDEP <: AbstractAgentCtrl
    state::Vector{Float64}
    ref_gen::AbstractTrajectoryGenerator
    reg::AbstractRegulator
end

## Functions
function control_law(reg::RegulatorP, agnt::AgentCtrlMAG, p_ref::Tuple{Float64, Float64})
    p_current = agnt.state[[1,3]];
    return (p_ref .- p_current) .* reg.P
end

function control_law(reg::RegulatorP, agnt::AgentCtrlDEP, p_ref::Tuple{Float64, Float64})
    p_current = agnt.state[[1,3]];

    # Compute the desired force in the 2D manipulation plane
    F_des_xy = (p_ref .- p_current) .* reg.P;

    # Saturate the desired force
    if norm(F_des_xy) > reg.F_sat
        F_des_xy = F_des_xy/norm(F_des_xy) .* reg.F_sat;
    end

    return (F_des_xy[1], F_des_xy[2], 3.34e-11)
end

function control_law(agnt::AbstractAgentCtrl, time)
    ref_pos = generate_reference(agnt.ref_gen, time, agent_pos(agnt));
    return control_law(agnt.reg, agnt, ref_pos);
end

function generate_reference(traj_gen::CircularTrajectoryGenerator, t, pos)
    x_traj = traj_gen.radius * sin(2*π*traj_gen.freq*t + traj_gen.phase) + traj_gen.offset[1];
    y_traj = traj_gen.radius * cos(2*π*traj_gen.freq*t + traj_gen.phase) + traj_gen.offset[2];

    x_traj, y_traj
end

function generate_reference(traj_gen::WaypointsTrajectoryGenerator, t, pos)
    if norm(traj_gen.points[traj_gen.current_point_index] .- pos) < traj_gen.switch_dist
        if (traj_gen.current_point_index + traj_gen.dir) > length(traj_gen.points) || (traj_gen.current_point_index + traj_gen.dir) < 1
            if traj_gen.repeat
                traj_gen.dir = -1*traj_gen.dir
            else
                if traj_gen.dir > 0
                    traj_gen.current_point_index = length(traj_gen.points)
                else
                    traj_gen.current_point_index = 1
                    traj_gen.dir = 0
                end
            end
        else
            traj_gen.current_point_index += traj_gen.dir
        end
    end

    traj_gen.points[traj_gen.current_point_index]
end

function agent_pos(agnt::AgentCtrlMAG)
    return agnt.state[1], agnt.state[3]
end

function agent_pos(agnt::AgentCtrlDEP)
    return agnt.state[1], agnt.state[3], agnt.state[5]
end

function ode_fun_MAG(du, u, p, t)
    # Calculate the approximate mass
    ballRadius = 10.0e-3
    mass = 4/3*π*ballRadius^3*8000    # 8000kg/m3 is the density of steel

    # Extract parameters
    actuatorArray, actuatorCommands = p

    # Calculate the force acting on the objects
    pos = (u[1], u[3])
    F = calcMAGForce(actuatorArray, pos, actuatorCommands)

    du[1] = u[2]
    du[2] = (-0.05*5*u[2] + F[1])/mass
    du[3] = u[4]
    du[4] = (-0.05*5*u[4] + F[2])/mass
end

function ode_fun_DEP(du, u, p, t)
    # Extract parameters
    actuatorArray, actuatorCommands = p

    # Calculate the force acting on the objects
    pos = (u[1], u[3], u[5])
    F = calcDEPForce(actuatorArray, pos, actuatorCommands)

    K = 5.4426e+05

    du[1] = K*F[1]
    du[2] = 0
    du[3] = K*F[2]
    du[4] = 0
    du[5] = 0
    du[6] = 0
end

function simulate!(agnt_ctrl::AgentCtrlMAG, aa::ActuatorArray, dt, actuatorCommands; solvopts...)
    sol = solve(ODEProblem(ode_fun_MAG, agnt_ctrl.state, (0, dt), (aa, actuatorCommands)); solvopts...)
    agnt_ctrl.state = sol.u[end]

    nothing
end

function simulate!(agnt_ctrl::AgentCtrlDEP, aa::ActuatorArray, dt, actuatorCommands; solvopts...)
    sol = solve(ODEProblem(ode_fun_DEP, agnt_ctrl.state, (0, dt), (aa, actuatorCommands)); solvopts...)
    agnt_ctrl.state = sol.u[end]

    nothing
end

end


text:
function DormandPrince853FieldIntegrator(arg0::Field, arg1::jdouble, arg2::jdouble, arg3::Vector{jdouble}, arg4::Vector{jdouble})
    return DormandPrince853FieldIntegrator((Field, jdouble, jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4)
end

function DormandPrince853FieldIntegrator(arg0::Field, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return DormandPrince853FieldIntegrator((Field, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function get_a(obj::DormandPrince853FieldIntegrator)
    return jcall(obj, "getA", Vector{Vector{RealFieldElement}}, ())
end

function get_b(obj::DormandPrince853FieldIntegrator)
    return jcall(obj, "getB", Vector{RealFieldElement}, ())
end

function get_c(obj::DormandPrince853FieldIntegrator)
    return jcall(obj, "getC", Vector{RealFieldElement}, ())
end

function get_order(obj::DormandPrince853FieldIntegrator)
    return jcall(obj, "getOrder", jint, ())
end


text:
#==== Plotting functions ======================================================#

using CSV
using DataFrames
using LaTeXStrings
using Query
using StatPlots

# function r(; name::AbstractString="samples", ks=0:Inf)

#     df = read(name*"txt")

# end

function ks(name::AbstractString="samples.txt")

    # We add the 'rows_for_type_detect' argument because there seems to be a bug
    # in the 'read' method: even with 'allowmissing=:all', a column whose first
    # missing value only appears after the type-detect rows will cause an error.

    df = CSV.read(name, rows_for_type_detect=typemax(Int))
    N = size(df,1)
    df = by(df, :k, d -> DataFrame(p=size(d,1)/N))
    @df df bar(:k, :p, xlabel="k", ylabel="P(k)", label="", size=(600,400))

end

function locs(name::AbstractString="samples.txt"; k::Int=1, ss=1:k)

    df = CSV.read(name, rows_for_type_detect=typemax(Int))
    df = @from r in df begin
         @where r.k == k
         @select r
         @collect DataFrame
    end
    ss = [Symbol("s$i") for i in ss]
    ps = vcat(
        [density(df[s], title=s, label="") for s in ss],
        [scatter(df[s], label="", ms=1) for s in ss]
    )
    m = length(ss)
    plot(ps..., layout=(2,m), size=(600m, 600))

end

function hts(name::AbstractString="samples.txt"; k::Int=1, hs=0:k)

    df = CSV.read(name, rows_for_type_detect=typemax(Int))
    df = @from r in df begin
         @where r.k == k
         @select r
         @collect DataFrame
    end
    hs = [Symbol("h$i") for i in hs]
    ps = vcat(
        [density(df[h], title=h, label="") for h in hs],
        [scatter(df[h], label="", ms=1) for h in hs]
    )
    m = length(hs)
    plot(ps..., layout=(2,m), size=(600m, 600))

end

# using Plots

# import Base.linspace

# #= I/O functions ==============================================================#

# function load(name::AbstractString)

#     ks = Vector{Int}()
#     open(name*".k") do f
#         ks = eval(parse(readline(f)))
#     end

#     K = maximum(find(ks)) - 1 # maximum sampled k.
#     S = [Matrix{Float64}(k+2,ks[k+1]) for k = 0:K] # step locations.
#     H = [Matrix{Float64}(k+1,ks[k+1]) for k = 0:K] # step heights.
#     js = zeros(Int, K+1) # track the samples per k.
#     open(name*".txt") do f
#         for (l,ln) in enumerate(eachline(f))
#             s,h = map(t -> eval(parse(t)), split(ln,';'))
#             k = length(s) - 2
#             j = (js[k+1] += 1)
#             S[k+1][:,j] = s
#             H[k+1][:,j] = h
#         end
#     end
#     return js,S,H

# end

# #= Statistical functions ======================================================#

# # Returns a range of linearly spaced elements over the data domain.
# linspace(n::Int, S) = linspace(0, S[end][end,end]-1, n)

# # Returns the corresponding step heights for the given (sorted) datums.
# function heights(xs,ss,hs)
#     is = similar(xs)
#     c = 1 # cursor.
#     for (i,s) in enumerate(@view(ss[2:end]))
#         j = c-1 + searchsortedfirst(@view(xs[c:end]), s)
#         is[c:j-1] = fill(hs[i],j-c)
#         c = j
#     end
#     return is
# end

# # Returns the sum of the step heights for given data, for a specific k.
# function hsum(xs,k,S,H)
#     ss,hs = S[k+1], H[k+1]
#     return (size(ss,2) == 0) ?
#         zeros(xs) : sum(i -> heights(xs,ss[:,i],hs[:,i]), indices(ss,2))
# end

# # Returns the posterior rate function at the given datums, for a specific k.
# function r(xs,k,S,H)
#     kn = size(S[k+1],2)
#     return (kn == 0) ? zeros(xs) : hsum(xs,k,S,H) / size(S[k+1],2)
# end

# # Returns the posterior rate function at the given datums.
# function r(xs,S,H)
#     K = length(S) - 1
#     ks = [size(S[k+1],2) for k = 0:K]
#     sum(k -> hsum(xs,k,S,H), 0:K) / sum(ks)
# end

# #= Plotting functions =========================================================#

# plot(xs, S, H) = Plots.plot(xs, r(xs,S,H), leg=false)
# plot(xs, k::Int, S, H) = Plots.plot(xs, r(xs,k,S,H), lab="k = $k")
# plot(xs, ks, S, H) = Plots.plot(xs, [r(xs,k,S,H) for k in ks],
#                                 lab=["k = $k" for k in ks])

max_stars_repo_name: kevdur/rjcp

max_issues_repo_name: kevdur/rjcp

max_forks_repo_name: kevdur/rjcp


text:
include("brent.jl")
include("goldensection.jl")


text:
import Cairo
import Cairo: CairoSurfaceImage
import Cairo: CairoSurface
import Cairo: CairoPattern

function image_surface(surface::CairoSurface)
    w, h = Int(surface.width), Int(surface.height)
    buffer = zeros(UInt32, w, h)
    img = Cairo.CairoImageSurface(buffer, Cairo.FORMAT_ARGB32, flipxy=false)
    cr = Cairo.CairoContext(img)
    Cairo.set_source_surface(cr, surface, 0, 0)
    Cairo.paint(cr)
    return img
end

function image_data(surface::CairoSurface)
    img = image_surface(surface)
    data = img.data
    Cairo.finish(img)
    Cairo.destroy(img)
    return data
end

function image_data(surface::CairoSurfaceImage)
    return surface.data
end

mutable struct Image <: Iterable
    surface::CairoSurfaceImage
    w::Int
    h::Int
    Image() = new()
end

function Image(surface::CairoSurface)
    this = Image()
    this.surface = image_surface(surface)
    this.w = surface.width
    this.h = surface.height
    return this
end

function Image(pixels::AbstractMatrix{UInt32})
    Image(Cairo.CairoImageSurface(pixels, Cairo.FORMAT_ARGB32, flipxy=false))
end

function Image(pixels::AbstractMatrix{ARGB32})
    Image(reinterpret(UInt32, pixels))
end

function Image(path::String)
    surface = Cairo.read_from_png(path)
    return Image(surface)
end

function drawimage(img::Image, x::Real = 0, y::Real = 0, w::Real = img.w, h::Real = img.h; centered = false)
    ctx = Luxor.get_current_cr()

    if centered
        x = x - w/2
        y = y - h/2
    end

    Cairo.image(ctx, img.surface, x, y, w, h)
end

function slice(image::Image, x::Int, y::Int, w::Int, h::Int)
    data = image.surface.data
    sx = (x + 1):w + x
    sy = (y + 1):h + y
    pixels = data[sx, sy]
    return Image(pixels)
end

function Base.getindex(image::Image, x::Int, y::Int, w::Int, h::Int)
    return slice(image, x, y, w, h)
end

function Base.replace!(this::Image, new::Image)
    this.surface = new.surface
    this.w = new.w
    this.h = new.h
end

Base.size(this::Image) = (this.w, this.h)

function rotl(this::Image)
    pixels = this.surface.data
    return Image(rotl90(pixels))
end

function rotl!(this::Image)
    new = rotl(this)
    replace!(this, new)
end

function rotr(this::Image)
    pixels = this.surface.data
    return Image(rotr90(pixels))
end

function rotr!(this::Image)
    new = rotr(this)
    replace!(this, new)
end

function Picture(image::Image, width = image.w, height = image.h)
    canvas = Canvas(width, height)
    ondraw(canvas) do c
        w, h = size(c)
        drawimage(image, 0, 0, w, h)
    end
    return canvas
end

Picture(image::AbstractString) = Picture(Image(image))

function Picture(image::AbstractString, width::Int, height::Int)
    return Picture(Image(image), width, height)
end


text:
import DataFrames: DataFrame
import StatsBase: countmap, weights, mean

# import Base: length
Base.length(x::Nothing) = 0
is_empty(x) = length(x) == 0
not_empty(x) = length(x) > 0

weighted_mean(x, w) = mean(x, weights(w))
weighted_sum(x, w) = sum(x, weights(w))

function nth(x, n)
    x = sort(x)
    x[n]
end

which_isna(x) = findall(x .== nothing)
which_notna(x) = findall(x .!= nothing)

"""
    match2(x, y)

# Examples
```julia
x = [1, 2, 3, 3, 4]
y = [0, 2, 2, 3, 4, 5, 6]
match2(x, y)
```

# Note: match2 only find the element in `y`
"""
function match2(x, y)
    # find x in y
    ind = indexin(x, y)
    I_x = which_notna(ind)
    I_y = something.(ind[I_x])
    # use `something` to suppress nothing `Union`
    DataFrame(value = x[I_x], I_x = I_x, I_y = I_y)
end

uniqueN(x) = length(unique(x))

# TODO: need to test
function CartesianIndex2Int(x, ind)
    # I = 1:prod(size(x))
    I = LinearIndices(x)
    I[ind]
end

table = countmap

"""
    duplicated(x::Vector{<:Real})

```julia
x = [1, 2, 3, 4, 1]
duplicated(x)
# [0, 0, 0, 0, 1]
```
"""
function duplicated(x::Vector)
    grps = table(x)
    grps = filter(x -> x[2] > 1, grps)

    n = length(x)
    res  = BitArray(undef, n)
    res .= false
    for (key, val) in grps
        k = 0
        for i = 1:n
            if x[i] == key
                k = k + 1
                if k >= 2; res[i] = true; end
                if k == val; break; end
            end
        end
    end
    res
end

seq_along(x) = 1:length(x)
seq_len(n) = 1:n

export table, which_isna, which_notna, match2, uniqueN, duplicated,
    is_empty, not_empty,
    mean, weighted_mean, weighted_sum,
    seq_along, seq_len


text:
# The simple
struct Newton{T1, Tlin, TSeq} <: QuasiNewton{T1}
   approx::T1
   linsolve::Tlin
   sequence::TSeq
end
# struct DefaultSequence end
Newton(;approx=Direct(), linsolve=\, sequence=NWI()) = Newton(approx, linsolve, sequence)


text:
"""
    gammaReg_log(topredict::String, varnames::Vector{Symbol}, train::DataFrame, test::DataFrame, RMSD::Bool=true)

Performs a gamma regression (log link) with variables from variables and data of D.

# Arguments

- `topredict::String`: The name of the column of the dataframe of interest that is to be predicted.
- `varnames::Vector{Symbol}`: The explanatory variables to include in the regression.
- `train::DataFrame`: The dataset used to train the model.
- `test::DataFrame`: The dataset used to test the model.
- `RMSD::Bool=true`: If true, returns the RMSD of the model. If false, returns the vector of predictions. Default = true.
"""
function gammaReg_log(topredict::String, varnames::Vector{Symbol}, train::DataFrame, test::DataFrame, RMSD::Bool=true)
    model = glm(Term(Symbol(topredict)) ~ sum(Term.(varnames)), train, Gamma(), LogLink())
    predictions = convert(Vector{Float64},GLM.predict(model, test))

    RMSD ? (return rmsd(predictions, test[:, Symbol(topredict)])) : return predictions
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
"""
    DistanceBasedCausalityTest{N}

The supertype of all abstract and composite types representing a causality
test based on some sort of distance computation.

The type parameter `N` indicates the number of returned elements when
applying the test.

Concrete subtypes are

- [`CrossMappingTest`](@ref)
- [`ConvergentCrossMappingTest`](@ref)
- [`JointDistanceDistributionTest`](@ref)
- [`JointDistanceDistributionTTest`](@ref)
- [`SMeasureTest`](@ref)

"""
abstract type DistanceBasedCausalityTest{N} <: CausalityTest end

export DistanceBasedCausalityTest