text:
using CUDA: CuArray, unsafe_free!
using Knet.CuArrays: cuarrays
using Knet.KnetArrays: cuallocator
using AutoGrad: Result, Node, Tape

# During the back pass we want to make pointers available as soon as we can to save memory
# without waiting for gc. This is risky as we have to make sure the pointers are not going to
# be used again.  We initialize a priority queue in gcnode_init() where CuArrays are mapped to
# the smallest index of the tape they are used.  If a CuArray is in a Param.value,
# Result.value, or Node{Param}.outgrad we do not want to free it, so we map those to 0. Only
# Node{Result}.outgrads can be safely freed, Result values are not safe because the user can
# assign them to global variables (as is done in RNNs).  Every time gcnode() is called after a
# back step, we update the queue with the pointers from the new outgrads of the node and its
# parents. The values in the queue are only decreased, never increased. Finally we free the
# pointers whose queue value is the current tape position. Note that tape indices go forward
# in time from parameters to loss but are processed backward in time from loss to parameters.

# The gcnode_queue maps CuArrays to the first index on tape they have a reference in. We use
# Base.Order.Reverse because we want to free the pointers with highest indices first.  Using
# ObjectId to allow fast hashing, hashing CuArrays directly is slow like Arrays.
const gcnode_queue = PriorityQueue{UInt,Int}(Base.Order.Reverse)

# To call unsafe_free gcnode_dict helps us find which CuArray belongs to an objectid.  Using
# WeakRef on CuArray values to allow garbage collection. Note that as a result sometimes
# WeakRef.value will metamorphose from CuArray to nothing if gc gets to it earlier.
const gcnode_dict = Dict{UInt,WeakRef}()

# We use node indices on the tape to use as values in the priority queue. Using
# ObjectId(::Node) for keys to allow gc after gcnode is done.
const gcnode_index = Dict{UInt,Int}()

# During the backward step parents of a node (who have lower indices) may have their
# outgrads modified, thus new references to CuArrays that we have already indexed may
# appear. We want to keep the smallest index for each CuArray.
function gcnode_setindex!(c::CuArray,v::Int)
    cid = objectid(c)
    get!(gcnode_dict, cid) do; WeakRef(c); end
    if v < get(gcnode_queue,cid,typemax(Int))
        gcnode_queue[cid] = v
    end
end

function gcnode_init(tape::Tape)
    empty!(gcnode_index)
    empty!(gcnode_queue)
    empty!(gcnode_dict)
    @inbounds for (i,n) in enumerate(tape.list)
        gcnode_index[objectid(n)] = i
        if n.Value isa Result
            for c in cuarrays(n.Value.value); gcnode_setindex!(c,0); end # pointers with index 0 will never get gc'ed
            for c in cuarrays(n.outgrad);     gcnode_setindex!(c,i); end # this only sets gcnode_queue[c] if it was not seen
        else # n.Value isa Param
            for c in cuarrays(n.Value.value); gcnode_setindex!(c,0); end
            for c in cuarrays(n.outgrad);     gcnode_setindex!(c,0); end
        end
    end
    # Mark this tape so we know when gcnode is called with a new tape and gcnode_init is needed
    tape.dict[gcnode_null] = gcnode_null_node
end

function gcnode(n::Node, tape::Tape)
    cuallocator[] || return knetgcnode(n,tape)
    if !haskey(tape.dict, gcnode_null)
        gcnode_init(tape)
    end
    ni = gcnode_index[objectid(n)]
    if n.Value isa Result
        for c in cuarrays(n.outgrad); gcnode_setindex!(c, ni); end
    end
    @inbounds for i in 1:length(n.parents);
        isassigned(n.parents, i) || continue
        parent = n.parents[i]
        if parent.Value isa Result
            pi = gcnode_index[objectid(parent)]
            for c in cuarrays(parent.outgrad); gcnode_setindex!(c, pi); end
        else
            for c in cuarrays(parent.outgrad); gcnode_setindex!(c,0); end # protect Params
        end
    end
    while !isempty(gcnode_queue) && peek(gcnode_queue)[2] >= ni  ## 5.62μs
        (cid,v) = dequeue_pair!(gcnode_queue)
        @assert v == ni
        c = gcnode_dict[cid].value
        if c isa CuArray ## c turns into nothing if gc'ed
            unsafe_free!(c)
        end
    end
    if n.Value isa Result
        n.Value, n.outgrad = gcnode_null, nothing
    end
end

const gcnode_null = Result{Nothing}(nothing,nothing,nothing,nothing)
const gcnode_null_node = Node(gcnode_null)


text:
mutable struct quadrotor_output_t <: LCMType
    timestamp::Int64
    position::SVector{3, Float64}
    orientation::SVector{4, Float64}
    twist::SVector{6, Float64}
    accelerometer::SVector{3, Float64}
    gyroscope::SVector{3, Float64}
end

@lcmtypesetup(quadrotor_output_t)

ext: jl


text:
using ACE1pack, Test, LazyArtifacts

test_train_set = joinpath(artifact"TiAl_tiny_dataset", "TiAl_tiny.xyz")

@info("Test constructing `data_params` and reading data")
params = data_params(fname = test_train_set, energy_key = "energy", force_key = "force", virial_key = "virial")
data = ACE1pack.read_data(params)


text:
using Pervii
using Documenter

DocMeta.setdocmeta!(Pervii, :DocTestSetup, :(using Pervii); recursive=true)

makedocs(;
    modules=[Pervii],
    authors="Nikita",
    repo="https://github.com/NikitaKutelev/Pervii.jl/blob/{commit}{path}#{line}",
    sitename="Pervii.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://NikitaKutelev.github.io/Pervii.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/NikitaKutelev/Pervii.jl",
)


text:
using RNLA, Test, LinearAlgebra, SparseArrays

@testset "rrpca_test1" begin
    m = 100
    k = 10
    n = 100
    p = 0.1

    L1 = rand(m, k) * rand(k, n)
    S1 = sprandn(m, n, p)
    A = L1 + S1 + rand(m, n)

    L2, S2, E2 = rrpca(A)

    @test norm(abs.(L1) - abs.(L2)) < sqrt(eps())
    @test norm(abs.(S1) - abs.(S2)) < sqrt(eps())
end


text:
@testset "arraymath" begin
    @testset "inv(::Matrix{$T})" for T in (Float64, ComplexF64)
        N = 3
        B = generate_well_conditioned_matrix(T, N)
        frule_test(inv, (B, randn(T, N, N)))
        rrule_test(inv, randn(T, N, N), (B, randn(T, N, N)))
    end

    @testset "*" begin
        @testset "Matrix-Matrix" begin
            dims = [3,4,5]
            for n in dims, m in dims, p in dims
                # don't need to test square case multiple times
                n > 3 && n == m == p && continue
                A = randn(m, n)
                B = randn(n, p)
                Ȳ = randn(m, p)
                rrule_test(*, Ȳ, (A, randn(m, n)), (B, randn(n, p)))
            end
        end
        @testset "Scalar-AbstractArray" begin
            for dims in ((3,), (5,4), (10,10), (2,3,4), (2,3,4,5))
                rrule_test(*, randn(dims), (1.5, 4.2), (randn(dims), randn(dims)))
                rrule_test(*, randn(dims), (randn(dims), randn(dims)), (1.5, 4.2))
            end
        end
    end

    @testset "$f" for f in (/, \)
        @testset "Matrix" begin
            for n in 3:5, m in 3:5
                A = randn(m, n)
                B = randn(m, n)
                Ȳ = randn(size(f(A, B)))
                rrule_test(f, Ȳ, (A, randn(m, n)), (B, randn(m, n)))
            end
        end
        @testset "Vector" begin
            x = randn(10)
            y = randn(10)
            ȳ = randn(size(f(x, y))...)
            rrule_test(f, ȳ, (x, randn(10)), (y, randn(10)))
        end
        if f == (\)
            @testset "Matrix $f Vector" begin
                X = randn(10, 4)
                y = randn(10)
                ȳ = randn(size(f(X, y))...)
                rrule_test(f, ȳ, (X, randn(size(X))), (y, randn(10)))
            end
            @testset "Vector $f Matrix" begin
                x = randn(10)
                Y = randn(10, 4)
                ȳ = randn(size(f(x, Y))...)
                rrule_test(f, ȳ, (x, randn(size(x))), (Y, randn(size(Y))))
            end
        end
    end
    @testset "/ and \\ Scalar-AbstractArray" begin
        A = randn(3, 4, 5)
        Ā = randn(3, 4, 5)
        Ȳ = randn(3, 4, 5)
        rrule_test(/, Ȳ, (A, Ā), (7.2, 2.3))
        rrule_test(\, Ȳ, (7.2, 2.3), (A, Ā))
    end
end


text:
macro consense(x, msgs...)
    msg = isempty(msgs) ? string(x, " has conflicting values") : :(string($(esc(msgs[1])), ": ", unique(v)))
    quote
        v = $(esc(x))
        @assert(all(broadcast(==, v, first(v))), $msg)
        first(v)
    end
end

macro adddevice(type, abstype, defaultformulation, name, axes, file)
    type = esc(type)
    abstype = esc(abstype)
    defaultformulation = esc(defaultformulation)
    quote
        struct $type{DF<:DeviceFormulation} <: $abstype{DF}
            model::EnergyModel
            class::Symbol
        end
        # Fall back to default formulation
        $type(model, class) = $type{$defaultformulation}(model, class)
        # Allow changing device formulation
        $type{DF}(d::$type) where DF = $type{DF}(d.model, d.class)
        $type{DF}(d::$type, ::Type{NewDF}) where {DF <: DeviceFormulation, NewDF <: DeviceFormulation} = $type{NewDF}(d)

        addcomponent($type, $(esc(name)), $(esc(axes)), $(esc(file)))
    end
end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
using Flux:onecold
using Test

@testset "onecold" begin
  a = [1, 2, 5, 3.]
  A = [1 20 5; 2 7 6; 3 9 10; 2 1 14]
  labels = ['A', 'B', 'C', 'D']

  @test onecold(a) == 3
  @test onecold(A) == [3, 1, 4]
  @test onecold(a, labels) == 'C'
  @test onecold(A, labels) == ['C', 'A', 'D']

  data = [:b, :a, :c]
  labels = [:a, :b, :c]
  hot = Flux.onehotbatch(data, labels)
  cold = onecold(hot, labels)

  @test cold == data
end

@testset "onehotbatch indexing" begin
  y = Flux.onehotbatch(ones(3), 1:10)
  @test y[:,1] isa Flux.OneHotVector
  @test y[:,:] isa Flux.OneHotMatrix
end


text:
"""
Explore Base.load_path() behaviour in diffferent situations.

Run this file as:

    julia init_me.jl
    julia --project=@. init_me.jl
    julia --project=. init_me.jl

Try running in a folder that has or does not have Project.toml file.
"""

MESSAGE = Dict(true => " (exists)", false => " (does not exist)")
exist(path)::Bool = isfile(path) || isdir(path)
printexist(path::String) = println("  ", path, MESSAGE[exist(path)])
printexist(nothing) = println("nothing")
printf2(s,n=8) = print("  " * s * ' '^(n-length(s)))

p_ = Base.JLOptions().project
project_option = (p_ != C_NULL) ? unsafe_string(p_) : "option not provided"
println("--project:\n  ", project_option)

println("load_path():")
for path in Base.load_path()
    printexist(path)
end

println("Base.active_project():")
printexist(Base.active_project())
println("Base.current_project():")
printexist(Base.current_project())

println("alias expansion with Base.load_path_expand():")
for alias in ["@.", "@", "@stdlib", "@v1.0"]
    printf2(alias)
    printexist(Base.load_path_expand(alias))
end

println("Base.DEPOT_PATH collects possible locations for registries:")
for p in Base.DEPOT_PATH
    printexist(p)
end


text:
__precompile__(true)

"""
Module `Approximations.jl` -- polygonal approximation of convex sets through
support vectors.
"""
module Approximations

using LazySets, Requires, LinearAlgebra, SparseArrays

export approximate,
       ballinf_approximation,
       box_approximation, interval_hull,
       decompose,
       overapproximate,
       box_approximation_symmetric, symmetric_interval_hull,
       BoxDirections,
       BoxDiagDirections,
       OctDirections,
       PolarDirections,
       SphericalDirections

const TOL(N::Type{Float64}) = eps(N)
const TOL(N::Type{Float32}) = eps(N)
const TOL(N::Type{Rational{INNER}}) where {INNER} = zero(N)

const DIR_EAST(N) = [one(N), zero(N)]
const DIR_NORTH(N) = [zero(N), one(N)]
const DIR_WEST(N) = [-one(N), zero(N)]
const DIR_SOUTH(N) = [zero(N), -one(N)]

include("iterative_refinement.jl")
include("box_approximations.jl")
include("template_directions.jl")
include("overapproximate.jl")
include("decompositions.jl")

end # module


text:
using Documenter, ChainMap

makedocs(
  modules = ChainMap
)

deploydocs(
  repo = "github.com/bramtayl/ChainMap.jl.git",
  deps = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
  osname = "osx"
)


text:
#this thing is useful for statmech and peps
#in principle we could let InfiniteMPS subtype from this thing
#but then we'd have to assert that numrows == 1 everywhere where it doesn't make sense ...

"
    2d extension of InfiniteMPS
"
struct MPSMultiline{A<:GenericMPSTensor,B<:MPSBondTensor}
    AL::PeriodicArray{A,2}
    AR::PeriodicArray{A,2}
    CR::PeriodicArray{B,2}
    AC::PeriodicArray{A,2}
end

Base.size(arr::MPSMultiline) = size(arr.AL)
Base.size(arr::MPSMultiline,i) = size(arr.AL,i)
Base.length(arr::MPSMultiline) = size(arr,1)
Base.eltype(arr::MPSMultiline) = eltype(arr.AL[1])
Base.lastindex(arr::MPSMultiline,i) = lastindex(arr.AL,i);
Base.similar(st::MPSMultiline) = MPSMultiline(similar(st.AL),similar(st.AR),similar(st.CR),similar(st.AC))
virtualspace(psi::MPSMultiline, a::Integer,b::Integer) = _firstspace(psi.AL[a,b])
function Base.convert(::Type{MPSMultiline},st::InfiniteMPS)
    convert(MPSMultiline,[st]);
end
function Base.convert(::Type{MPSMultiline},v::AbstractVector{T}) where T<:InfiniteMPS{A,B} where {A,B}
    ALs = PeriodicArray{A}(undef,length(v),length(v[1]));
    ARs = PeriodicArray{A}(undef,length(v),length(v[1]));
    CRs = PeriodicArray{B}(undef,length(v),length(v[1]));
    ACs = PeriodicArray{A}(undef,length(v),length(v[1]));

    for (i,row) in enumerate(v)
        ALs[i,:] = row.AL[:];
        ARs[i,:] = row.AR[:];
        CRs[i,:] = row.CR[:];
        ACs[i,:] = row.AC[:];
    end

    MPSMultiline(ALs,ARs,CRs,ACs);
end

function Base.convert(::Type{InfiniteMPS},st::MPSMultiline{A,B}) where {A,B}
    @assert size(st,1) == 1 #otherwise - how would we convert?
    convert(Vector,st)[1]
end

function Base.convert(::Type{Vector},st::MPSMultiline{A,B}) where {A,B}
    map(1:size(st,1)) do row
        InfiniteMPS{A,B}(   PeriodicArray(st.AL[row,:]),
                            PeriodicArray(st.AR[row,:]),
                            PeriodicArray(st.CR[row,:]),
                            PeriodicArray(st.AC[row,:]));
    end
end

function MPSMultiline(A::AbstractArray{T,2};tol = Defaults.tolgauge,maxiter = Defaults.maxiter,leftgauged=false) where T <: GenericMPSTensor

    ACs = PeriodicArray{T,2}(undef,size(A,1),size(A,2));
    ALs = PeriodicArray{T,2}(undef,size(A,1),size(A,2));
    ARs = PeriodicArray{T,2}(undef,size(A,1),size(A,2));;

    ctype = typeof(TensorMap(rand,eltype(A[1,1]),space(A[1,1],1),space(A[1,1],1)))
    Cs = PeriodicArray{ctype,2}(undef,size(A,1),size(A,2));

    @sync for row in 1:size(A,1)
        @Threads.spawn begin
            if !leftgauged
                tal,_,deltal= uniform_leftorth(PeriodicArray(A[row,:]); tol = tol, maxiter = maxiter)
            else
                tal = PeriodicArray(A[row,:]);
            end
            tar,tc,deltar = uniform_rightorth(tal; tol = tol, maxiter = maxiter)

            ALs[row,:] = tal[:];
            ARs[row,:] = tar[:];
            Cs[row,:] = circshift(tc[:],-1);

            for loc = 1:length(tal)
                ACs[row,loc] = ALs[row,loc]*Cs[row,loc]
            end
        end
    end

    return MPSMultiline(ALs,ARs,Cs,ACs)
end

l_RR(state::MPSMultiline,row,loc::Int=1) = @tensor toret[-1;-2]:=state.CR[row,loc-1][1,-2]*conj(state.CR[row,loc-1][1,-1])
l_RL(state::MPSMultiline,row,loc::Int=1) = state.CR[row,loc-1]
l_LR(state::MPSMultiline,row,loc::Int=1) = state.CR[row,loc-1]'
l_LL(state::MPSMultiline{A},row,loc::Int=1) where A= isomorphism(Matrix{eltype(A)}, space(state.AL[row,loc],1),space(state.AL[row,loc],1))

r_RR(state::MPSMultiline{A},row,loc::Int=length(state)) where A= isomorphism(Matrix{eltype(A)},domain(state.AR[row,loc]),domain(state.AR[row,loc]))
r_RL(state::MPSMultiline,row,loc::Int=length(state)) = state.CR[row,loc]'
r_LR(state::MPSMultiline,row,loc::Int=length(state)) = state.CR[row,loc]
r_LL(state::MPSMultiline,row,loc::Int=length(state))= @tensor toret[-1;-2]:=state.CR[row,loc][-1,1]*conj(state.CR[row,loc][-2,1])


text:
"""
	contour(cmd0::String="", arg1=[]; kwargs...)

Reads a table data and produces a raw contour plot by triangulation.

Full option list at [`contour`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html)

Parameters
----------

- $(GMT.opt_J)
- **A** : **annot** : -- Str or Number --       Flags = [-|[+]annot_int][labelinfo]

    *annot_int* is annotation interval in data units; it is ignored if contour levels are given in a file.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#a)
- $(GMT.opt_B)
- **C** : **cont** : **contours** : **levels** : -- Str or Number or GMTcpt --  Flags = [+]cont_int

    Contours contours to be drawn may be specified in one of three possible ways.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#c)
- **D** : **dump** : -- Str --

    Dump contours as data line segments; no plotting takes place.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#d)
- **E** : **index** : -- Str or Mx3 array --

    Give name of file with network information. Each record must contain triplets of node
    numbers for a triangle.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#e)
- **G** : **labels** : -- Str --

    Controls the placement of labels along the quoted lines.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#g)
- **I** : **colorize** : -- Bool or [] --

    Color the triangles using the color scale provided via **C**.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#i)
- $(GMT.opt_Jz)
- **L** : **mesh** : -- [] or Str or Number --

    Draw the underlying triangular mesh using the specified pen attributes (if not provided, use default pen)
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#l)
- **N** : **no_clip** : --- Bool or [] --

    Do NOT clip contours or image at the boundaries [Default will clip to fit inside region].
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#n)
- $(GMT.opt_P)
- **Q** : **cut** : -- Str or Number --         Flags = [cut[unit]][+z]]

    Do not draw contours with less than cut number of points.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#q)
- **S** : **skip** : -- Str or [] --            Flags = [p|t]

    Skip all input xyz points that fall outside the region.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#s)
- **T** : **ticks** : -- Str --                 Flags = [+|-][+a][+dgap[/length]][+l[labels]]

    Draw tick marks pointing in the downward direction every *gap* along the innermost closed contours.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#t)
- $(GMT.opt_R)
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** : **pen** : -- Str or Number --

    Sets the attributes for the particular line.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#w)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- **Z** : **scale** : -- Str --

    Use to subtract shift from the data and multiply the results by factor before contouring starts.
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/pscontour.html#z)
- $(GMT.opt_bi)
- $(GMT.opt_bo)
- $(GMT.opt_d)
- $(GMT.opt_di)
- $(GMT.opt_do)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function contour(cmd0::String="", arg1=[]; K=false, O=false, first=true, kwargs...)

    arg2 = []       # May will contain a CPT or a Mx3 indices array
    arg3 = []       # May will contain a Mx3 indices array
	N_args = isempty_(arg1) ? 0 : 1

	length(kwargs) == 0 && return monolitic("pscontour", cmd0, arg1)

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	cmd  = ""
    cmd, opt_B, opt_J, opt_R = parse_BJR(d, cmd, "", O, " -JX12c/0")
	cmd  = parse_UVXY(cmd, d)
	cmd, opt_bi = parse_bi(cmd, d)
	cmd, = parse_bo(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_di(cmd, d)
	cmd, = parse_do(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, opt_i = parse_i(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd, = parse_swap_xy(cmd, d)
	cmd  = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	# If file name sent in, read it and compute a tight -R if this was not provided
    cmd, arg1, opt_R, opt_i = read_data(d, cmd0, cmd, arg1, opt_R, opt_i, opt_bi, opt_di)

    cmd, arg1, arg2, N_args = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_args, arg1, arg2)

	cmd = add_opt(cmd, 'A', d, [:A :annot])
	if (!occursin(" -C", cmd))			# Otherwise ignore an eventual :cont because we already have it
		cmd = add_opt(cmd, 'C', d, [:cont :contours :levels])
	end

	for sym in [:E :index]
		if (haskey(d, sym))
			if (isa(d[sym], Array{Int}))
				cmd = string(cmd, " -E")
				# Now need to find the free slot where to store the indices array
				if (N_args == 0)		arg1 = d[sym]
				elseif (N_args == 1)	arg2 = d[sym]
				else					arg3 = d[sym]
				end
			else
				cmd = string(cmd, " -E", arg2str(d[sym]))
			end
			break
		end
	end

	cmd = add_opt(cmd, 'D', d, [:D :dump])
	cmd = add_opt(cmd, 'G', d, [:G :labels])
	cmd = add_opt(cmd, 'I', d, [:I :colorize])
	cmd = add_opt(cmd, 'L', d, [:L :mesh])
	cmd = add_opt(cmd, 'N', d, [:N :no_clip])
	cmd = add_opt(cmd, 'Q', d, [:Q :cut])
	cmd = add_opt(cmd, 'S', d, [:S :skip])
	cmd = add_opt(cmd, 'T', d, [:T :ticks])
	cmd = add_opt(cmd, 'W', d, [:W :pen])

	if (!occursin(" -W", cmd) && !occursin(" -I", cmd))		# Use default pen
		cmd = cmd * " -W"
	end

	cmd = finish_PS(d, cmd, output, K, O)
    return finish_PS_module(d, cmd, "-D", output, fname_ext, opt_T, K, "pscontour", arg1, arg2, arg3)
end

# ---------------------------------------------------------------------------------------------------
contour!(cmd0::String="", arg1=[]; K=true, O=true, first=false, kw...) =
	contour(cmd0, arg1; K=true, O=true, first=false, kw...)

contour(arg1, cmd0::String=""; K=false, O=false, first=true, kw...) =
	contour(cmd0, arg1; K=K, O=O, first=first, kw...)

contour!(arg1, cmd0::String=""; K=true, O=true, first=false, kw...) =
	contour(cmd0, arg1; K=true, O=true, first=false, kw...)

# ---------------------------------------------------------------------------------------------------
pscontour  = contour
pscontour! = contour!


text:
###############################
# kmedoids, Voronoi iteration #
###############################

function kmedoids_voronoi(D::Matrix{Float}, k::Int;
                    verbose::Bool = true,
                    epsi::Float = Float(1e-10),
                    avg_iters::Int = Int(1),
                    max_iters::Int = Int(1000),
                    calc_emp::Bool = true,
                    medoid_init::Vector{Int} = Int[],
                    init::Symbol = :per)

  n = size(D,1)

  # initialize optimization
  opt = Float[]
  queue = Array{Float, 1}()
  diff_ratio = Float(1)
  avg_objective = Float(0)
  medoid_idx = zeros(Int, k)

  max_val = maximum(D)
  if verbose
    println("initialization...")
  end
  if init == :rand
    init_idx = 1:n
  elseif init == :cen
    ~, init_idx = findmin(max(D,max_val*eye(n)),2)
    ~, init_idx = ind2sub(size(D), init_idx[:])
  elseif init == :per
    ~, init_idx = findmax(D,2)
    ~, init_idx = ind2sub(size(D), init_idx[:])
  else
    error("no such initialization")
  end
  medoid_idx = (length(medoid_init) == k) ? medoid_init : sample(init_idx, k, replace = false)
  ~,clust_idx = findmin(D[:, medoid_idx],2)

  if verbose
    @printf("%-20s %-20s %-20s\n", "iteration", "objective", "difference ratio");
  end
  local clust_idx
  local objective
  local min_clust_idx, min_medoid_idx, t
  min_obj = 1.
  mu = mean(D)
  for t in 0:max_iters

    if t > 0
      # find new medoids
      for c in 1:k
        cidx = find(clust_idx .== c)
        medoid_idx[c] = (length(cidx) > 0) ? cidx[indmin(mean(D[cidx, cidx],1))] : 0
      end
      # reassign empty medoids
      empt_ind = medoid_idx .== 0
      empt_num = sum(empt_ind)
      if empt_num > 0
        medoid_idx[empt_ind] = sample(init_idx, empt_num, replace = (length(init_idx) <= empt_num))
      end
    end

    # nearest cluster
    min_dist,clust_idx = findmin(D[:, medoid_idx],2)
    ~, clust_idx = ind2sub(size(D), clust_idx[:])
    # objective
    objective = mean(fixnan.(min_dist./mu))

    if objective < min_obj
      min_obj = objective
      min_clust_idx = clust_idx
      min_medoid_idx = medoid_idx
    end

    if t > 0
      # finalization
      Base.push!(opt, objective)

      if t > avg_iters
        diff_ratio = abs(objective - avg_objective)/fixdiv(avg_objective);
      end

      Base.push!(queue, objective/avg_iters)
      avg_objective = avg_objective + objective/avg_iters
      if t > avg_iters
        avg_objective = avg_objective - shift!(queue)
      end
    end

    if verbose
      if t > avg_iters
        @printf("%-20d %-20.6f %-20.6f\n", t, objective, diff_ratio)
      else
        @printf("%-20d %-20.6f %-20s\n", t, objective, "-")
      end
    end

    if (diff_ratio < epsi)
      break
    end

  end

  return min_clust_idx, min_medoid_idx, min_obj, t
end

#################
# kmedoids, PAM #
#################

function kmedoids_pam(D::Matrix{Float}, k::Int;
                    verbose::Bool = true,
                    avg_iters::Int = Int(10),
                    max_iters::Int = Int(1000),
                    min_obj::Float = eps(Float),
                    calc_emp::Bool = true,
                    medoid_init::Vector = Int[],
                    init::Symbol = :per)

  n = size(D,1)

  # initialize optimization
  opt = Float[]
  queue = Array{Float, 1}()
  diff_ratio = Float(1)
  avg_objective = Float(0)
  medoid_idx = zeros(Int, k)

  max_val = maximum(D)
  if length(medoid_init) != k
    if verbose
      println("initialization...")
    end
    if init == :rand
      init_idx = 1:n
    elseif init == :cen
      ~, init_idx = findmin(max(D,max_val*eye(n)),2)
      ~, init_idx = ind2sub(size(D), init_idx[:])
    elseif init == :per
      ~, init_idx = findmax(D,2)
      ~, init_idx = ind2sub(size(D), init_idx[:])
    else
      error("no such initialization")
    end
    medoid_idx = sample(init_idx, k, replace = false)
  else
    medoid_idx = deepcopy(medoid_init)
  end
  mu = mean(D)
  local clust_idx
  local objective
  new_min_dist = zeros(n)

  if verbose
    @printf("%-20s %-20s\n", "iteration", "objective");
  end

  total_samp = 0
  best_samp = 0
  switched = true
  t = 0
  while (switched) && (t<=max_iters)
    t += 1

    ###############
    # build phase #
    ###############

    # nearest mediod
    min_dist,clust_idx = findmin(D[:, medoid_idx],2)
    ~, clust_idx = ind2sub(size(D), clust_idx[:])
    # objective
    objective = fixnan(mean(min_dist)./mu)
    total_samp += 1
    best_samp = total_samp

    if verbose
      if t > avg_iters
        @printf("%-20d %-20.6f\n", t, objective)
      else
        @printf("%-20d %-20s\n", t, objective)
      end
    end

    if objective < min_obj
      break
    end

    ##############
    # swap phase #
    ##############

    # find new medoids
    switched = false
    for i in randperm(k)
      new_min_dist[:] = minimum(D[:, medoid_idx[1:k .!= i]],2)
      for j in setdiff(randperm(n), medoid_idx)
        new_objective = fixnan(mean(min.(new_min_dist, D[:, j]))./mu)
        total_samp += 1
        if new_objective < objective
          switched = true
          medoid_idx[i] = j
        end
        if switched
          break
        end
      end
      if switched
        break
      end
    end
  end

  return clust_idx, medoid_idx, objective, best_samp, total_samp
end

nothing


text:
export BrickletAmbientLightV3IlluminanceCallbackConfiguration
struct BrickletAmbientLightV3IlluminanceCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletAmbientLightV3Configuration
struct BrickletAmbientLightV3Configuration
    illuminance_range::Integer
    integration_time::Integer
end

export BrickletAmbientLightV3SPITFPErrorCount
struct BrickletAmbientLightV3SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletAmbientLightV3Identity
struct BrickletAmbientLightV3Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAmbientLightV3
"""
Measures ambient light up to 64000lux
"""
mutable struct BrickletAmbientLightV3 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAmbientLightV3(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_ambient_light_v3")
        deviceInternal = package.BrickletAmbientLightV3(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_illuminance
"""
    $(SIGNATURES)

Returns the illuminance of the ambient light sensor. The measurement range goes
up to about 100000lux, but above 64000lux the precision starts to drop.
The illuminance is given in lux/100, i.e. a value of 450000 means that an
illuminance of 4500lux is measured.

An illuminance of 0lux indicates an error condition where the sensor cannot
perform a reasonable measurement. This can happen with very dim or very bright
light conditions. In bright light conditions this might indicate that the sensor
is saturated and the configuration should be modified (:func:`Set Configuration`)
to better match the conditions.

If you want to get the value periodically, it is recommended to use the
:cb:`Illuminance` callback. You can set the callback configuration
with :func:`Set Illuminance Callback Configuration`.
"""
function get_illuminance(device::BrickletAmbientLightV3)
    return device.deviceInternal.get_illuminance()
end

export set_illuminance_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Illuminance` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Illuminance` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_illuminance_callback_configuration(device::BrickletAmbientLightV3, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_illuminance_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_illuminance_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Illuminance Callback Configuration`.
"""
function get_illuminance_callback_configuration(device::BrickletAmbientLightV3)
    return device.deviceInternal.get_illuminance_callback_configuration()
end

export set_configuration
"""
    $(SIGNATURES)

Sets the configuration. It is possible to configure an illuminance range
between 0-600lux and 0-64000lux and an integration time between 50ms and 400ms.

The unlimited illuminance range allows to measure up to about 100000lux, but
above 64000lux the precision starts to drop.

A smaller illuminance range increases the resolution of the data. A longer
integration time will result in less noise on the data.

If the actual measure illuminance is out-of-range then the current illuminance
range maximum +0.01lux is reported by :func:`Get Illuminance` and the
:cb:`Illuminance` callback. For example, 800001 for the 0-8000lux range.

With a long integration time the sensor might be saturated before the measured
value reaches the maximum of the selected illuminance range. In this case 0lux
is reported by :func:`Get Illuminance` and the :cb:`Illuminance` callback.

If the measurement is out-of-range or the sensor is saturated then you should
configure the next higher illuminance range. If the highest range is already
in use, then start to reduce the integration time.
"""
function set_configuration(device::BrickletAmbientLightV3, illuminance_range, integration_time)
    device.deviceInternal.set_configuration(illuminance_range, integration_time)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletAmbientLightV3)
    return device.deviceInternal.get_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletAmbientLightV3)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletAmbientLightV3, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletAmbientLightV3)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletAmbientLightV3, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletAmbientLightV3, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletAmbientLightV3, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletAmbientLightV3)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletAmbientLightV3)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletAmbientLightV3)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletAmbientLightV3, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletAmbientLightV3)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAmbientLightV3)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletAmbientLightV3, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
    LinearDecaySchedule
A schedule that linearly decreases a value from `start` to `stop` in `steps` steps.
if the value is greater or equal to `stop`, it stays constant.

# Constructor

`LinearDecaySchedule(;start, stop, steps)`
"""
@with_kw struct LinearDecaySchedule{R<:Real} <: Function
    start::R
    stop::R
    steps::Int
end

function (schedule::LinearDecaySchedule)(k)
    rate = (schedule.start - schedule.stop) / schedule.steps
    val = schedule.start - k*rate
    val = max(schedule.stop, val)
end

"""
    ExplorationPolicy <: Policy
An abstract type for exploration policies.
Sampling from an exploration policy is done using `action(exploration_policy, on_policy, k, state)`.
`k` is a value that is used to determine the exploration parameter. It is usually a training step in a TD-learning algorithm.
"""
abstract type ExplorationPolicy <: Policy end

"""
    loginfo(::ExplorationPolicy, k)
returns information about an exploration policy, e.g. epsilon for e-greedy or temperature for softmax.
It is expected to return a namedtuple (e.g. (temperature=0.5)). `k` is the current training step that is used to compute the exploration parameter.
"""
function loginfo end

"""
    EpsGreedyPolicy <: ExplorationPolicy

represents an epsilon greedy policy, sampling a random action with a probability `eps` or returning an action from a given policy otherwise.
The evolution of epsilon can be controlled using a schedule. This feature is useful for using those policies in reinforcement learning algorithms.

# Constructor:

`EpsGreedyPolicy(problem::Union{MDP, POMDP}, eps::Union{Function, Float64}; rng=Random.GLOBAL_RNG, schedule=ConstantSchedule)`

If a function is passed for `eps`, `eps(k)` is called to compute the value of epsilon when calling `action(exploration_policy, on_policy, k, s)`.

# Fields

- `eps::Function`
- `rng::AbstractRNG`
- `actions::A` an indexable list of action
"""
struct EpsGreedyPolicy{T<:Function, R<:AbstractRNG, A} <: ExplorationPolicy
    eps::T
    rng::R
    actions::A
end

function EpsGreedyPolicy(problem, eps::Function;
                         rng::AbstractRNG=Random.GLOBAL_RNG)
    return EpsGreedyPolicy(eps, rng, actions(problem))
end
function EpsGreedyPolicy(problem, eps::Real;
                         rng::AbstractRNG=Random.GLOBAL_RNG)
    return EpsGreedyPolicy(x->eps, rng, actions(problem))
end

function POMDPs.action(p::EpsGreedyPolicy, on_policy::Policy, k, s)
    if rand(p.rng) < p.eps(k)
        return rand(p.rng, p.actions)
    else
        return action(on_policy, s)
    end
end

loginfo(p::EpsGreedyPolicy, k) = (eps=p.eps(k),)

# softmax
"""
    SoftmaxPolicy <: ExplorationPolicy

represents a softmax policy, sampling a random action according to a softmax function.
The softmax function converts the action values of the on policy into probabilities that are used for sampling.
A temperature parameter or function can be used to make the resulting distribution more or less wide.

# Constructor

`SoftmaxPolicy(problem, temperature::Union{Function, Float64}; rng=Random.GLOBAL_RNG)`

If a function is passed for `temperature`, `temperature(k)` is called to compute the value of the temperature when calling `action(exploration_policy, on_policy, k, s)`

# Fields

- `temperature::Function`
- `rng::AbstractRNG`
- `actions::A` an indexable list of action

"""
struct SoftmaxPolicy{T<:Function, R<:AbstractRNG, A} <: ExplorationPolicy
    temperature::T
    rng::R
    actions::A
end

function SoftmaxPolicy(problem, temperature::Function;
                       rng::AbstractRNG=Random.GLOBAL_RNG)
    return SoftmaxPolicy(temperature, rng, actions(problem))
end
function SoftmaxPolicy(problem, temperature::Real;
                       rng::AbstractRNG=Random.GLOBAL_RNG)
    return SoftmaxPolicy(x->temperature, rng, actions(problem))
end

function POMDPs.action(p::SoftmaxPolicy, on_policy::Policy, k, s)
    vals = actionvalues(on_policy, s)
    vals ./= p.temperature(k)
    maxval = maximum(vals)
    exp_vals = exp.(vals .- maxval)
    exp_vals /= sum(exp_vals)
    return p.actions[sample(p.rng, Weights(exp_vals))]
end

loginfo(p::SoftmaxPolicy, k) = (temperature=p.temperature(k),)

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
module SplitDataset

using DiskArrays
using Glob
using HDF5
using JSON3
using Printf
using PaddedBlocks

struct H5SplitDataset{T,N} <: AbstractDiskArray{T, N}
    files::Array{String,1}
    blocks::Blocks{N}
    stackname::String
end

haschunks(a::H5SplitDataset) = DiskArrays.Chunked()
eachchunk(a::H5SplitDataset) = DiskArrays.GridChunks(a, a.blocks.block_size)
Base.size(a::H5SplitDataset) = a.blocks.full_size

"Collects filenames and block dimensions for a HDF5 dataset"
function file_blocks_names(path)
    return h5open(path, "r") do file
        return Blocks(
            Int64.(tuple(collect_attributes(file["stack_metadata"]["full_size"])...)[end:-1:1]),
            Int64.(tuple(collect_attributes(file["stack_metadata"]["block_size"])...)[end:-1:1]),
        ),
        collect_attributes(file["stack_metadata"]["files"])
    end
end

function stack_metadata_from_json(path)
    mddict = open(path, "r") do f
        return JSON3.read(f)
    end
    return Blocks(
        Int64.(Tuple(mddict["shape_full"]))[end:-1:1],
        Int64.(Tuple(mddict["shape_block"]))[end:-1:1],
    )
end

function H5SplitDataset(folder::String, prefix::String = "")
    jsonpath = joinpath(folder, "stack_metadata.json")
    files = sort(glob(prefix * "*.h5", folder))
    if isfile(jsonpath)
        blocks = stack_metadata_from_json(jsonpath)
        filenames = files
    else
        blocks, filenames = file_blocks_names(files[end])
    end
    n_dim = length(blocks.full_size)

    stack_name = h5open(files[1]) do f
        if "stack" in names(f)
            return "stack"
        else
            return "stack_$(n_dim)D"
        end
    end

    el1 = h5read(files[1], stack_name, tuple((1:1 for _ in 1:n_dim)...))
    return H5SplitDataset{eltype(el1),n_dim}(files, blocks, stack_name)
end

function H5SplitDataset(T, path, full_size::NTuple{N, Int}, block_size::NTuple{N, Int}) where {N}
    blocks = Blocks(Int64.(full_size), Int64.(block_size))
    files = [joinpath(path, @sprintf("%04d.h5", i)) for i in 0:(length(blocks)-1)]
    if !isdir(path)
        mkpath(path)
    end
    open(joinpath(path, "stack_metadata.json"), "w") do f
        return JSON3.write(
            f,
            Dict(
                "shape_full" => blocks.full_size[end:-1:1],
                "shape_block" => blocks.block_size[end:-1:1],
            ),
        )
    end
    return H5SplitDataset{T, length(full_size)}(files, blocks, "stack")
end

function H5SplitDataset(
    folder,
    a::AbstractArray{T,N},
    block_size::NTuple{N,Int},
) where {N,T}
    blocks = Blocks(size(a), block_size)
    mkpath(folder)
    files = String[]
    for (i_block, sl) in enumerate(slices(blocks))
        filename = @sprintf("%04d.h5", i_block)
        push!(files, filename)
        savepath = joinpath(folder, filename)
        el1 = h5write(savepath, "stack_$(N)D", a[sl...])
    end
    open(joinpath(folder, "stack_metadata.json"), "w") do f
        return JSON3.write(
            f,
            Dict(
                "shape_full" => blocks.full_size[end:-1:1],
                "shape_block" => blocks.block_size[end:-1:1],
            ),
        )
    end
    return H5SplitDataset(folder)
end

"Collects attributes form a group in a deepdish-saved HDF5 file"
function collect_attributes(hgroup)
    atrnms = names(attrs(hgroup))
    nums = sort([
        begin
            mtch = match(r"i(\d+)", nam)
            parse(mtch.captures[1])
        end for nam in filter(n -> ismatch(r"i(\d+)", n), atrnms)
    ])
    return [read(attrs(hgroup)["i$i"]) for i in nums]
end

function get_matching_slices(f_idx, file_limits, block_limits, block_size)
    source_slices = tuple((
        (ci == s ? si : 1):(ci == e ? ei : bs)
        for
        (ci, s, e, si, ei, bs) in zip(
            f_idx,
            file_limits[1, :],
            file_limits[2, :],
            block_limits[1, :],
            block_limits[2, :],
            block_size
        ))...)

    read_size = tuple((sl.stop - sl.start + 1 for sl in source_slices)...)

    rel_idx = f_idx .- tuple(file_limits[1, :]...) .+ 1
    target_slices = tuple((
        begin
            st = (st_idx == 1 ? 1 : (st_idx - 2) * bs + bs - first_idx + 2)
            st:(st+sz-1)
        end for
        (st_idx, bs, first_idx, sz) in
        zip(rel_idx, block_size, block_limits[1, :], read_size)
    )...)
    return source_slices, target_slices
end

function DiskArrays.readblock!(dset::H5SplitDataset{T, N}, target, i::AbstractUnitRange...) where {T,N}
    idx = i
    file_limits, block_limits = PaddedBlocks.blocks_to_take(dset.blocks, idx)
    s2i = LinearIndices(dset.blocks.blocks_per_dim)
    for f_idx in
        Iterators.product((s:e for (s, e) in zip(file_limits[1, :], file_limits[2, :]))...)
        i_file = s2i[f_idx...]
        source_slices, target_slices = get_matching_slices(f_idx, file_limits, block_limits, dset.blocks.block_size)
        ar = h5read(dset.files[i_file], dset.stackname, source_slices)
        target[target_slices...] .= ar

    end
end

function DiskArrays.writeblock!(dset::H5SplitDataset{T, N}, target, i::AbstractUnitRange...) where {T,N}
    idx = i
    file_limits, block_limits = PaddedBlocks.blocks_to_take(dset.blocks, idx)
    s2i = LinearIndices(dset.blocks.blocks_per_dim)
    for f_idx in
        Iterators.product((s:e for (s, e) in zip(file_limits[1, :], file_limits[2, :]))...)
        i_file = s2i[f_idx...]

        source_slices, target_slices = get_matching_slices(f_idx, file_limits, block_limits, dset.blocks.block_size)

        # TODO handle wrong existing datasets
        h5open(dset.files[i_file], "cw") do fid
            if !exists(fid, dset.stackname)
                h5dset = d_create(fid, dset.stackname, datatype(T),
                dataspace(dset.blocks.block_size...),
                "chunk", HDF5.heuristic_chunk(T, (dset.blocks.block_size)), "blosc", 3)
            else
                h5dset = fid[dset.stackname]
            end

            h5dset[source_slices...] = target[target_slices...]
        end

    end
end

export H5SplitDataset

end


text:
"""Module for the modelling of pollutants' dispersion in the atmosphere."""
module Plumes

using ArchGDAL
using ArgParse
using Dates

include(".\\Utils\\Functions.jl")
include(".\\Utils\\FunctionsDB.jl")

export run_plume

const agd = ArchGDAL

mutable struct Plume <: Functions.AbstractAnalysisObject
  # Parameters
    concentration::Float64      # Rate of chemical emission (m³/sec)
    x::Float64                  # Distance in wind direction, x coordinate (m)
    y::Float64                  # Y coordinate of target point (m)
    z::Float64                  # Height of target point, Z Coordinate (m)
    stability::String           # Atmosphere Pasquill class stability
    outdoor::String             # Outdoor class
    stack_height::Float64       # Height of the stack source of the pollutants
    stack_diameter::Float64     # Diameter of the stack source of the pollutants
    direction::Int64            # Main wind direction (°)
    wind_speed::Float64         # Wind speed in the main direction
    gas_velocity::Float64       # Gas velocity
    gas_temperature::Float64    # Absolute gas temperature
    temperature::Float64        # Absolute ambient air temperature
    max_domain::Float64         # Maximum concentration found during analysis
 # Computational results
    H::Float64                  # Real source height
    σy::Float64                 # Gaussian distribution horizontal standard deviation
    σz::Float64                 # Gaussian distribution vertical standard deviation
    g1::Float64                 # Gaussian horizontal distribution factor
    g2::Float64                 # Gaussian vertical distribution factor

    Plume(concentration,x,y,z,stability,outdoor,stack_height,stack_diameter,direction,wind_speed,gas_velocity,gas_temperature,temperature,max_domain)=new(concentration,x,y,z,stability,outdoor,stack_height,stack_diameter,direction,wind_speed,gas_velocity,gas_temperature,temperature,max_domain,0.0,0.0,0.0,0.0,0.0)
end

function calc_h!( p::Plume )
    try
        fb = 9.81 * ( (p.stack_diameter * p.gas_velocity) / 4.0 ) * ( ( p.gas_temperature / p.temperature ) / p.gas_temperature )
        Δh = 1.6 * fb^0.333333 * p.x^0.666667
        p.H = p.stack_height + Δh
    catch
        p.H = p.stack_height
    end
    return p.H
end

function calc_σ!( p::Plume )
    σy1, σy2, σyexp, σz1, σz2, σzexp = FunctionsDB.air_extract(p.stability, p.outdoor, ["sigmay1", "sigmay2", "sigmayexp", "sigmaz1", "sigmaz2", "sigmazexp"])[1, :]
    p.σy, p.σz = @. ( (σy1, σz1) * p.x ) / ( ( 1.0 + (σy2, σz2) * p.x )^(σyexp, σzexp) )
    return p.σy, p.σz
end

function calc_g!( p::Plume )
    p.g1 = ℯ^( ( -0.5 * p.y^2.0 ) / p.σy^2.0 )
    p.g2 = ℯ^( ( -0.5 * (p.z - p.stack_height)^2.0 ) / p.σz^2.0 ) + ℯ^( ( -0.5 * (p.z + p.stack_height)^2.0 ) / p.σz^2.0 )
    return p.g1, p.g2
end

function calc_C!( p::Plume )
    return (100.0p.concentration / 3600.0p.wind_speed) * ( (p.g1 * p.g2) / (2.0π * p.σy * p.σz) )
end

function Functions.compute_concentration!( p::Plume )
    if p.x > 0
        calc_σ!(p)
        calc_g!(p)
        calc_h!(p)
        cfinal = calc_C!(p)
        if cfinal > p.max_domain
            p.max_domain = cfinal
        end
        return cfinal
    end
    return 0.0
end

"""
    run_plume(; dem_file::String, source_file::String, stability::String, outdoor::String, concentration::Float64, tollerance::Int64=2, resolution::Int64,
                wind_direction::Int64, wind_speed::Float64, stack_height::Float64, stack_diameter::Float64=0.0, gas_velocity::Float64=0.0, gas_temperature::Float64=0.0,
                temperature::Float64=0.0, output_path::AbstractString=".\\plume_otput_model.tiff" )

Create and save as `output_path` a raster containing the results of model of dispersion of airborne pollutants.

# Arguments
- `dem_file::String`: path to the raster containing the height of the terrain in each cell.
- `source_file::String`: path to the shapefile containing the source point of the plume.
- `stability::String`: atmosphere Pasquill class stability.
- `outdoor::String`: type of environment, either `\"c\"` (country) or `\"u\"` (urban).
- `concentration::Float64`: rate of chemical emission.
- `tollerance::Int64=2`: value used to determine wether the concentration of pollutant in a cell is relevant.
    Specifically, a concentration value is considered relevant if its value is within "tollerance" orders of magnitute from the concentration on other cells.
- `resolution::Int64`: size of the cell in meters.
- `wind_direction::Int64`: angle of main direction of the wind in degrees.
- `wind_speed::Float64`: average wind speed in the main direction.
- `stack_height::Float64 `: height of the stack, or height of the source of the plume.
- `stack_diameter::Float64=0.0`: diameter of the stack emiting the plume.
- `gas_velocity::Float64=0.0`: gas velocity.
- `gas_temperature::Float64=0.0`: absolute temperature of the gas.
- `temperature::Float64=0.0`: absolute ambient air temperature.
- `output_path::String=".\\plume_otput_model.tiff"`: output file path.
"""
function run_plume(; dem_file::String, source_file::String, stability::String, outdoor::String, concentration::Float64, tollerance::Int64=2, resolution::Float64,
                     wind_direction::Int64, wind_speed::Float64, stack_height::Float64, stack_diameter::Float64=0.0, gas_velocity::Float64=0.0, gas_temperature::Float64=0.0,
                     temperature::Float64=0.0, output_path::AbstractString=".\\plume_otput_model.tiff" )

    if outdoor ∉ ["c", "o"]
        throw(DomainError(outdoor, "`outdoor` must either be `\"c\"` or `\"u\"`."))
    end

    geom = agd.getgeom(collect(agd.getlayer(agd.read(source_file), 0))[1])

    if agd.geomdim(geom) != 0
        throw(DomainError(source_file, "The shapefile must contain a point."))
    end

    dem = agd.read(dem_file)
    refsys = agd.getproj(dem)
    geotransform = agd.getgeotransform(dem)

    if refsys != agd.toWKT(agd.getspatialref(geom))
        throw(DomainError("The reference systems are not uniform. Aborting analysis."))
    end

    x_source = agd.getx(geom, 0)
    y_source = agd.gety(geom, 0)
    r_source, c_source = Functions.toIndexes(dem, x_source, y_source)
    plume = Plume(concentration, x_source, y_source, agd.getband(dem, 1)[r_source, c_source], stability, outdoor, stack_height, stack_diameter, wind_direction, wind_speed, gas_velocity, gas_temperature, temperature, 0.0)
    points = Functions.expand(r_source, c_source, concentration, tollerance, dem, plume)

    maxR = maximum( point -> point[1], points )
    minR = minimum( point -> point[1], points )
    maxC = maximum( point -> point[2], points )
    minC = minimum( point -> point[2], points )

    geotransform[[1, 4]] = Functions.toCoords(dem, minR, minC)
    noData = Float32(agd.getnodatavalue(agd.getband(dem, 1)))
    data = fill(noData, maxR-minR+1, maxC-minC+1)
    if isnothing(data) || isempty(data)
        throw(ErrorException("Analysis failed"))
    end
    @inbounds for r in minR:maxR, c in minC:maxC
        match = findfirst( p -> p[1] == r && p[2] == c, points )
        if !isnothing(match)
            data[r-minR+1, c-minC+1] = points[match][3]
        end
    end
    Functions.writeRaster(data, agd.getdriver("GTiff"), geotransform, refsys, noData, output_path)
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger
using .Numbers: is_palindrome

function compute(n::Int)::Int
    result = 0
    for i ∈ 1:n
        number = BigInteger(string(i))
        number += BigInteger(reverse(number.str))
        j = 1
        while j ≤ 50 && !is_palindrome(number)
            number += BigInteger(reverse(number.str))
            j += 1
        end
        if j > 50
            result += 1
        end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
mutable struct TreeLeaf
  key::Union{Symbol,AbstractString,Number}
  value::Any
end

mutable struct TreeNode
  key::Union{Symbol,AbstractString,Number}
  children::Vector{Union{TreeNode, TreeLeaf}}
end

mutable struct TreeList
  root::Union{TreeNode, Nothing}
  nodes::Vector{Union{TreeNode, TreeLeaf}}
end

function build_node!(tree_list::TreeList, cur_key::Union{Symbol,AbstractString,Number}, nestedDict::AbstractDict)
  cur_children = []
  for (cur_key, cur_value) in nestedDict
    tmp_node = build_node!(tree_list, cur_key, cur_value)
    push!(cur_children, tmp_node)
  end

  cur_node = TreeNode(cur_key, cur_children)
  push!(tree_list.nodes, cur_node)
  return cur_node
end

function build_node!(tree_list::TreeList, cur_key::Union{Symbol,AbstractString,Number}, cur_value)
  cur_node = TreeLeaf(cur_key, cur_value)
  push!(tree_list.nodes, cur_node)
  return cur_node
end

function build_tree(nested_dict::AbstractDict)
  tree_list = TreeList(nothing, [])
  tree_list.root = build_node!(tree_list, :root, nested_dict)

  return tree_list
end

function tree_to_dict(cur_node::TreeLeaf)
  return cur_node.value
end

function tree_to_dict(cur_node::TreeNode)
  cur_dict = OrderedDict()
  for sub_node in cur_node.children
    cur_dict[sub_node.key] = tree_to_dict(sub_node)
  end
  return cur_dict
end


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

JLLWrappers.@generate_main_file_header("LLVMLibUnwind")
JLLWrappers.@generate_main_file("LLVMLibUnwind", UUID("47c5dbc3-30ba-59ef-96a6-123e260183d9"))
end  # module LLVMLibUnwind_jll


text:
addprocs(3)
import POPT

@everywhere function widescaleRandomNoiseTL(x::Vector{Float64})
    n = length(x)
    expon = rand()*6 - 4
    displacement = randn(n)
    new_x = x + displacement*10^expon
end

@everywhere function f(x)
	sum(abs.(x))
end

dims = 50
start_x = 1000*randn(dims)
micro_steps = 100
macro_steps = 10
temps = [16, 4]

(m_12, m_22, b_12, b_22, s_12, s_22) = POPT.runblock(start_x, f, temps, widescaleRandomNoiseTL, micro_steps, macro_steps)
(m_1, m_2, b_1, b_2, s_1, s_2) = POPT.runnonblock(start_x, f, temps, widescaleRandomNoiseTL, micro_steps, macro_steps)
val_gt, x_gt, p_val_gt, swaps_made = PTChain.greedyTempering(start_x, f, temps, widescaleRandomNoiseTL, micro_steps, macro_steps)

println("---Blocking")

println("\t m_1: $(m_12) (actual value: $(f(b_12))")
println("\t m_2: $(m_22) (actual value: $(f(b_22))")
println("\t made $(s_12) swaps on 1 and $(s_22) swaps on 2")

println("---Non Blocking")
println("\t m_1: $(m_1) (actual value: $(f(b_1))")
println("\t m_2: $(m_2) (actual value: $(f(b_2))")
println("\t made $(s_1) swaps on 1 and $(s_2) swaps on 2")

println("---Running sequential")
start_x = 1000*randn(dims)
println("\t Sequential best was: $(val_gt)")
println("\t made $(swaps_made) swaps")

max_stars_repo_name: jacob-ruth/pt-optimization

max_issues_repo_name: jacob-ruth/pt-optimization

max_forks_repo_name: jacob-ruth/pt-optimization


text:
# Dataset Name
abstract type Name end

struct MandatoryMethod <: Exception
    method::String
    N::Type{<:Name}
end

function Base.showerror(io::IO, e::MandatoryMethod)
    print(io, "MandatoryMethod: method `$(e.method)(::Type{$(nameof(e.N))})` not defined")
end

## Mandatory methods
format(N::Type{<:Name}) = throw(MandatoryMethod("format", N))
task(N::Type{<:Name}) = throw(MandatoryMethod("problem", N))
source(N::Type{<:Name}) = throw(MandatoryMethod("source", N))
downloadlink(N::Type{<:Name}) = throw(MandatoryMethod("downloadlink", N))
preprocess(N::Type{<:Name}) = throw(MandatoryMethod("preprocess", N))
nattributes(N::Type{<:Name}) = throw(MandatoryMethod("nattributes", N))
ninstances(N::Type{<:Name}) = throw(MandatoryMethod("ninstances", N))

## Optional methods
checksum(::Type{<:Name}) = ""
datadepname(N::Type{<:Name})= string(nameof(N))
name(N::Type{<:Name})= string(nameof(N))
author(::Type{<:Name}) = String[]
licence(::Type{<:Name}) = ""
citation(::Type{<:Name}) = ""

## Registration
function description(N::Type{<:Name})
    return """
    Dataset: $(name(N))
    Author(s): $(join(author(N), ", "))
    Website: $(source(N))
    Licence: $(licence(N))
    Sample size: $(join(nattributes(N), "x"))
    Train/validation/test: $(join(ninstances(N), "/"))

    Please respect copyright and use data responsibly. For more details about copyright and license, visit the website mentioned above.
    """
end

function make_datadep(N::Type{<:Name})
    return DataDep(
        datadepname(N),
        description(N),
        downloadlink(N),
        checksum(N);
        post_fetch_method = preprocess(N),
    )
end

# Dataset Task
abstract type Task end
args(::Type{<:Task}; kwargs...) = (;)
postprocess(::Type{<:Name}, ::Type{<:Task}, data; kwargs...) = data

# Dataset Format
abstract type Format end
args(::Type{<:Format}; kwargs...) = (;)
postprocess(::Type{<:Name}, ::Type{<:Format}, data; kwargs...) = data

abstract type Split end

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
## NOTE: Some constants are defined in mxpr_type.jl
"""
    SYMATA_PACKAGE_PATH

Top level directory of the Symata package.
"""
const SYMATA_PACKAGE_PATH = joinpath(dirname(@__FILE__), "..")

"""
    SYMATA_LANG_TEST_PATH

Filesystem path to directory containing tests written in Symata.
"""
const SYMATA_LANG_TEST_PATH = joinpath(SYMATA_PACKAGE_PATH, "symata_test")

const NullMxpr = mxprcf(:Null)
const Null = :Null  # In Mma, Null is a Symbol. But, the Mma REPL prints nothing when encountering it (sometimes)

const ComplexInfinity = mxprcf(:DirectedInfinity)
const Infinity = mxprcf(:DirectedInfinity,1)
const MinusInfinity = mxprcf(:DirectedInfinity,-1)
setsymval(:Infinity, Infinity)
setsymval(:ComplexInfinity, ComplexInfinity)
mergesyms(Infinity,:nothing)
mergesyms(MinusInfinity,:nothing)
mergesyms(ComplexInfinity,:nothing)

function Base.convert(::Type{AbstractFloat}, mx::Mxpr{:DirectedInfinity})
    if length(mx) == 1
        margs(mx)[1] == 1 && return Inf
        margs(mx)[1] == -1 && return -Inf
    end
    symerror("Can't convert ", wrapout(mx), " to float.") # symerror already does wrapout
end

function AbstractFloat(mx::Mxpr{:DirectedInfinity})
    if length(mx) == 1
        margs(mx)[1] == 1 && return Inf
        margs(mx)[1] == -1 && return -Inf
    end
    symerror("Can't convert ", wrapout(mx), " to float.") # symerror already does wrapout
end

const Indeterminate = :Indeterminate
const Undefined = :Undefined

const I = complex(0,1)
const MinusI = complex(0,-1)

setsymval(:True, true)
setsymval(:False, false)  ## Took a long time for this bug to manifest. :False was true for a thousand commits.

for s in (:BigInt, :BigFloat, :Float64, :Int64, :Int, :DataType, :Char, :String)
    sq = QuoteNode(s)
    @eval setsymval($sq, $s)
end
setsymval(:DateTime, Dates.DateTime)

## 1/Sqrt(2) and -1/Sqrt(2)
## Making these constant does save a bit of time
## We could do this in an more organized way.
## 1/Sqrt(2)
const _moosq2 = mmul(-1,mpow(2,-1//2))
## -1/Sqrt(2)
const _oosq2 = mpow(2,-1//2)
## Pi^2
const Pisq = mmul(Pi,Pi)
## Log(2)
const Log2 = Log(2)
for s in (:_moosq2, :_oosq2 , :Pisq, :ComplexInfinity, :Infinity, :MinusInfinity, :NullMxpr, :Log2)
    @eval mergesyms($s, :nothing)
    @eval setfixed($s)
    @eval setcanon($s)
end


text:
module BenchCasing
using BenchmarkTools
using SymSpellChecker: get_opcodes, transfer_casing_for_similar_text

suite = BenchmarkGroup()
s1 = "qabxcd"
s2 = "abycdf"
suite["opcodes"] = @benchmarkable get_opcodes($s1, $s2)

text_w_casing = "Haaw is the weeather in New York?"
text_wo_casing = "how is the weather in new york?"
suite["transfer_casing"] = @benchmarkable transfer_casing_for_similar_text($text_w_casing, $text_wo_casing)

end

BenchCasing.suite


text:
using LinearAlgebra: normalize, ⋅, norm
export GraphPlot, graphplot, graphplot!

"""
    graphplot(graph::AbstractGraph)
    graphplot!(ax, graph::AbstractGraph)

Creates a plot of the network `graph`. Consists of multiple steps:
- Layout the nodes: see `layout` attribute. The node position is accessible from outside
  the plot object `p` as an observable using `p[:node_pos]`.
- plot edges as `edgeplot`-plot
- if `arrows_show` plot arrowheads as `scatter`-plot
- plot nodes as `scatter`-plot
- if `nlabels!=nothing` plot node labels as `text`-plot
- if `elabels!=nothing` plot edge labels as `text`-plot

The main attributes for the subplots are exposed as attributes for `graphplot`.
Additional attributes for the `scatter`, `edgeplot` and `text` plots can be provided
as a named tuples to `node_attr`, `edge_attr`, `nlabels_attr` and `elabels_attr`.

Most of the arguments can be either given as a vector of length of the
edges/nodes or as a single value. One might run into errors when changing the
underlying graph and therefore changing the number of Edges/Nodes.

## Attributes
### Main attributes
- `layout=Spring()`: function `AbstractGraph->Vector{Point}` determines the base layout
- `node_color=scatter_theme.color`
- `node_size=scatter_theme.markersize`
- `node_marker=scatter_theme.marker`
- `node_attr=(;)`: List of kw arguments which gets passed to the `scatter` command
- `edge_color=lineseg_theme.color`: Color for edges.
- `edge_width=lineseg_theme.linewidth`: Pass a vector with 2 width per edge to
  get pointy edges.
- `edge_attr=(;)`: List of kw arguments which gets passed to the `linesegments` command
- `arrow_show=Makie.automatic`: `Bool`, indicate edge directions with arrowheads?
  Defaults to `true` for `SimpleDiGraph` and `false` otherwise.
- `arrow_size=scatter_theme.markersize`: Size of arrowheads.
- `arrow_shift=0.5`: Shift arrow position from source (0) to dest (1) node.
- `arrow_attr=(;)`: List of kw arguments which gets passed to the `scatter` command

### Node labels
The position of each label is determined by the node position plus an offset in
data space.

- `nlabels=nothing`: `Vector{String}` with label for each node
- `nlabels_align=(:left, :bottom)`: Anchor of text field.
- `nlabels_distance=0.0`: Pixel distance from node in direction of align.
- `nlabels_color=labels_theme.color`
- `nlabels_offset=nothing`: `Point` or `Vector{Point}` (in data space)
- `nlabels_textsize=labels_theme.textsize`
- `nlabels_attr=(;)`: List of kw arguments which gets passed to the `text` command

### Edge labels
The base position of each label is determinded by `src + shift*(dst-src)`. The
additional `distance` parameter is given in pixels and shifts the text away from
the edge.

- `elabels=nothing`: `Vector{String}` with label for each edge
- `elabels_align=(:center, :bottom)`: Anchor of text field.
- `elabels_distance=0.0`: Pixel distance of anchor to edge.
- `elabels_shift=0.5`: Position between src and dst of edge.
- `elabels_opposite=Int[]`: List of edge indices, for which the label should be
  displayed on the opposite side
- `elabels_rotation=nothing`: Angle of text per label. If `nothing` this will be
  determined by the edge angle!
- `elabels_offset=nothing`: Additional offset in data space
- `elabels_color=labels_theme.color`
- `elabels_textsize=labels_theme.textsize`
- `elabels_attr=(;)`: List of kw arguments which gets passed to the `text` command

### self edges & curvy edges
- `edge_plottype=Makie.automatic()`: Either `automatic`, `:linesegments` or
  `:beziersegments`. `:beziersegments` are much slower for big graphs!
- `selfedge_size=Makie.automatic()`: Size of self-edge-loop (dict/vector possible).
- `selfedge_direction=Makie.automatic()`: Direction of self-edge-loop as `Point2` (dict/vector possible).
- `selfedge_width=Makie.automatic()`: Opening of selfloop in rad (dict/vector possible).
- `tangents=nothing`:

    Specify a pair of tangent vectors per edge (for src and dst). If `nothing`
    (or edge idx not in dict) draw a straight line.

- `tfactor=0.6`:

    Factor is used to calculate the bezier waypoints from the (normalized) tangents.
    Higher factor means bigger radius. Can be tuple per edge to specify different
    factor for src and dst.

- `waypoints=nothing`

    Specify waypoints for edges. This parameter should be given as a vector or
    dict. Waypoints will be crossed using natural cubic splines. The waypoints may
    or may not include the src/dst positions.

- `waypoint_radius=nothing`: If number (dict/vector possible) bent lines within radius of waypoints.
"""
@recipe(GraphPlot, graph) do scene
    # TODO: figure out this whole theme business
    scatter_theme = default_theme(scene, Scatter)
    lineseg_theme = default_theme(scene, LineSegments)
    labels_theme = default_theme(scene, Makie.Text)
    Attributes(
        layout = Spring(),
        # node attributes (Scatter)
        node_color = scatter_theme.color,
        node_size = scatter_theme.markersize,
        node_marker = scatter_theme.marker,
        node_attr = (;),
        # edge attributes (LineSegements)
        edge_color = lineseg_theme.color,
        edge_width = lineseg_theme.linewidth,
        edge_attr = (;),
        # arrow attributes (Scatter)
        arrow_show = automatic,
        arrow_size = scatter_theme.markersize,
        arrow_shift = 0.5,
        arrow_attr = (;),
        # node label attributes (Text)
        nlabels = nothing,
        nlabels_align = (:left, :bottom),
        nlabels_distance = 0.0,
        nlabels_color = labels_theme.color,
        nlabels_offset = nothing,
        nlabels_textsize = labels_theme.textsize,
        nlabels_attr = (;),
        # edge label attributes (Text)
        elabels = nothing,
        elabels_align = (:center, :bottom),
        elabels_distance = 0.0,
        elabels_shift = 0.5,
        elabels_opposite = Int[],
        elabels_rotation = nothing,
        elabels_offset = nothing,
        elabels_color = labels_theme.color,
        elabels_textsize = labels_theme.textsize,
        elabels_attr = (;),
        # self edge attributes
        edge_plottype = automatic,
        selfedge_size = automatic,
        selfedge_direction = automatic,
        selfedge_width = automatic,
        tangents=nothing,
        tfactor=0.6,
        waypoints=nothing,
        waypoint_radius=nothing,
    )
end

function Makie.plot!(gp::GraphPlot)
    graph = gp[:graph]

    # create initial vertex positions, will be updated on changes to graph or layout
    # make node_position-Observable available as named attribute from the outside
    gp[:node_pos] = @lift [Pointf(p) for p in ($(gp.layout))($graph)]

    node_pos = gp[:node_pos]

    # create array of pathes triggered by node_pos changes
    # in case of a graph change the node_position will change anyway
    edge_paths = lift(node_pos, gp.selfedge_size,
                      gp.selfedge_direction, gp.selfedge_width) do pos, s, d, w
        find_edge_paths(graph[], gp.attributes, pos)
    end

    sc = Makie.parent_scene(gp)

    # function which projects the point in px space
    to_px = lift(sc.px_area, sc.camera.projectionview) do pxa, pv
        # project should transform to 2d point in px space
        (point) -> project(sc, point)
    end
    # get angle in px space from path p at point t
    to_angle = @lift (path, p0, t) -> begin
        # TODO: maybe shorter tangent? For some perspectives this might give wrong angles in 3d
        p1 = p0 + tangent(path, t)
        tpx = $to_px(p1) - $to_px(p0)
        atan(tpx[2], tpx[1])
    end

    # plot edges
    edge_plot = edgeplot!(gp, edge_paths;
                          plottype=gp.edge_plottype,
                          color=gp.edge_color,
                          linewidth=gp.edge_width,
                          gp.edge_attr...)

    # plott arrow heads
    arrow_pos = @lift if !isempty($edge_paths)
        broadcast(interpolate, $edge_paths, $(gp.arrow_shift))
    else # if no edges return (empty) vector of points, broadcast yields Vector{Any} which can't be plotted
        Vector{eltype(node_pos[])}()
    end
    arrow_rot = @lift Billboard(broadcast($to_angle, edge_paths[], $arrow_pos, gp.arrow_shift[]))
    arrow_show = @lift $(gp.arrow_show) === automatic ? $graph isa SimpleDiGraph : $(gp.arrow_show)
    arrow_heads = scatter!(gp,
                           arrow_pos,
                           marker = '➤',
                           markersize = gp.arrow_size,
                           color = gp.edge_color,
                           rotations = arrow_rot,
                           strokewidth = 0.0,
                           markerspace = Pixel,
                           visible = arrow_show,
                           gp.arrow_attr...)

    # plot vertices
    vertex_plot = scatter!(gp, node_pos;
                           color=gp.node_color,
                           marker=gp.node_marker,
                           markersize=gp.node_size,
                           gp.node_attr...)

    # plot node labels
    if gp.nlabels[] !== nothing
        positions = @lift begin
            if $(gp.nlabels_offset) != nothing
                $node_pos .+ $(gp.nlabels_offset)
            else
                copy($node_pos)
            end
        end

        offset = @lift if $(gp.nlabels_align) isa Vector
            $(gp.nlabels_distance) .* align_to_dir.($(gp.nlabels_align))
        else
            $(gp.nlabels_distance) .* align_to_dir($(gp.nlabels_align))
        end

        nlabels_plot = text!(gp, gp.nlabels;
                             position=positions,
                             align=gp.nlabels_align,
                             color=gp.nlabels_color,
                             offset=offset,
                             textsize=gp.nlabels_textsize,
                             gp.nlabels_attr...)
    end

    # plot edge labels
    if gp.elabels[] !== nothing
        # positions: center point between nodes + offset + distance*normal + shift*edge direction
        positions = @lift begin
            pos = broadcast(interpolate, $edge_paths, $(gp.elabels_shift))

            if $(gp.elabels_offset) !== nothing
                pos .= pos .+ $(gp.elabels_offset)
            end
            pos
        end

        # rotations based on the edge_vec_px and opposite argument
        rotation = @lift begin
            if $(gp.elabels_rotation) isa Real
                # fix rotation to a single angle
                rot = $(gp.elabels_rotation)
            else
                # determine rotation for each position
                rot = broadcast($to_angle, edge_paths[], $positions, gp.elabels_shift[])

                for i in $(gp.elabels_opposite)
                    rot[i] += π
                end
                # if there are user provided rotation for some labels, use those
                if $(gp.elabels_rotation) isa Vector
                    for (i, α) in enumerate($(gp.elabels_rotation))
                        α !== nothing && (rot[i] = α)
                    end
                end
            end
            return rot
        end

        # calculate the offset in pixels in normal direction to the edge
        offsets = @lift begin
            tangent_px = broadcast(edge_paths[], $positions, gp.elabels_shift[]) do path, p0, t
                p1 = p0 + tangent(path, t)
                $to_px(p1) - $to_px(p0)
            end

            offsets = map(p -> Point(-p.data[2], p.data[1])/norm(p), tangent_px)
            offsets .= $(gp.elabels_distance) .* offsets
            for i in $(gp.elabels_opposite)
                offsets[i] = -1.0 * offsets[i] # flip offset if in opposite
            end
            offsets
        end

        elabels_plot = text!(gp, gp.elabels;
                             position=positions,
                             rotation=rotation,
                             offset=offsets,
                             align=gp.elabels_align,
                             color=gp.elabels_color,
                             textsize=gp.elabels_textsize,
                             gp.elabels_attr...)
    end

    return gp
end

"""
    find_edge_paths(g, attr, pos::AbstractVector{PT}) where {PT}

Returns an `AbstractPath` for each edge in the graph.
"""
function find_edge_paths(g, attr, pos::AbstractVector{PT}) where {PT}
    paths = Vector{AbstractPath{PT}}(undef, ne(g))

    ne(g) == 0 && return paths

    for (i, e) in enumerate(edges(g))
        if e.src == e.dst
            size = getattr(attr.selfedge_size, i)
            direction = getattr(attr.selfedge_direction, i)
            width = getattr(attr.selfedge_width, i)
            paths[i] = selfedge_path(g, pos, e.src, size, direction, width)
        else
            p1, p2 = pos[e.src], pos[e.dst]
            tangents = getattr(attr.tangents, i)
            tfactor = getattr(attr.tfactor, i)
            waypoints::Vector{PT} = getattr(attr.waypoints, i, PT[])
            if isnothing(waypoints) || isempty(waypoints)
                paths[i] = Path(p1, p2; tangents, tfactor)
            else
                # the waypoints may already include the endpoints
                waypoints[begin] == p1 && popfirst!(waypoints)
                waypoints[end] == p2 && pop!(waypoints)

                radius = getattr(attr.waypoint_radius, i, nothing)

                if isempty(waypoints) || radius === nothing || radius === :spline
                    paths[i] = Path(p1, waypoints..., p2; tangents, tfactor)
                elseif radius isa Real
                    paths[i] = Path(radius, p1, waypoints..., p2)
                else
                    throw(ArgumentError("Invalid radius $radius for edge $i!"))
                end
            end
        end
    end

    # try to narrow down the typ, i.e. just `Line`s
    T = mapreduce(typeof, promote_type, paths)
    return convert(Vector{T}, paths)
end

"""
    selfedge_path(g, pos, v, size, direction, width)

Return a Path for the
"""
function selfedge_path(g, pos::AbstractVector{<:Point2}, v, size, direction, width)
    vp = pos[v]
    # get the vectors to all the neighbors
    ndirs = [pos[n] - vp for n in neighbors(g, v) if n != v]

    # angle and maximum width of loop
    γ, Δ = 0.0, 0.0

    if direction === automatic && !isempty(ndirs)
        angles = SVector{length(ndirs)}(atan(p[2], p[1]) for p in ndirs)
        angles = sort(angles)

        for i in 1:length(angles)
            α = angles[i]
            β = get(angles, i+1, 2π + angles[1])
            if β-α > Δ
                Δ = β-α
                γ = (β+α) / 2
            end
        end

        # set width of selfloop
        Δ = min(.7*Δ, π/2)
    elseif direction === automatic && isempty(ndirs)
        γ = π/2
        Δ = π/2
    else
        @assert direction isa Point2 "Direction of selfedge should be 2 dim vector ($direction)"
        γ = atan(direction[2], direction[1])
        Δ = π/2
    end

    if width !== automatic
        Δ = width
    end

    # the size (max distance to v) of loop
    # if there are no neighbors set to 0.5. Else half dist to nearest neighbor
    if size === automatic
        size = isempty(ndirs) ? 0.5 : minimum(norm.(ndirs)) * 0.5
    end

    # the actual length of the tagent vectors, magic number from `CurveTo`
    l = Float32( size/(cos(Δ/2) * 2*0.375) )
    t1 = vp + l * Point2f(cos(γ-Δ/2), sin(γ-Δ/2))
    t2 = vp + l * Point2f(cos(γ+Δ/2), sin(γ+Δ/2))

    return BezierPath([MoveTo(vp),
                       CurveTo(t1, t2, vp)])
end

function selfedge_path(g, pos::AbstractVector{<:Point3}, v, size, direction, width)
    error("Self edges in 3D not yet supported")
end

"""
    edgeplot(paths::Vector{AbstractPath})
    edgeplot!(sc, paths::Vector{AbstractPath})

Recipe to draw the edges. Attribute `plottype` can be either

- `:linesegments`: Draw edges as `linesegments` (just straight lines)
- `:beziersegments`: Draw edges as `beziersegments` (separate `lines` plot for
  each edge. Much slower than `linesegments`!)
- `Makie.automatic()`: Only use `beziersegements` if there is at least one curvy line.
"""
@recipe(EdgePlot, paths) do scene
    Attributes(
        plottype = automatic
    )
end

function Makie.plot!(p::EdgePlot)
    N = length(p[:paths][])
    # remove plottype from attributes, otherwise attributs can't be passed on
    plottype = pop!(p.attributes, :plottype)[]
    if plottype === automatic
        justlines = all(isline, p[:paths][])

        if justlines
            plottype = :linesegments
        elseif N > 500
            @warn "Since there are a lot of edges ($N), they will be drawn as straight lines "*
                "even though they contain curvy edges. If you realy wan't to plot them as "*
                "bezier curves pass `edge_plottype=:beziersegments` explicitly. This will have "*
                "much worse performance!"
            plottype = :linesegments
        else
            plottype = :beziersegments
        end
    end

    if plottype === :linesegments
        PT = ptype(eltype(p[:paths][]))
        segs = Observable(Vector{PT}(undef, 2*N))

        update_segments!(segs, p[:paths][]) # first call to initialize
        on(p[:paths]) do paths # update if pathes change
            update_segments!(segs, paths)
        end
        linesegments!(p, segs; p.attributes...)
    elseif plottype === :beziersegments
        beziersegments!(p, p[:paths]; p.attributes...)
    else
        error("Unknown plottype $(repr(plottype)) for edge plot")
    end

    return p
end

function update_segments!(segs, paths)
    N = length(paths)
    if length(segs[]) != 2*N
        resize!(segs[], 2*N)
    end
    for (i, p) in enumerate(paths)
        segs[][2*i - 1] = interpolate(p, 0.0)
        segs[][2*i]     = interpolate(p, 1.0)
    end
    notify(segs)
end

"""
    beziersegments(paths::Vector{AbstractPath})
    beziersegments!(sc, paths::Vector{AbstractPath})

Recipe to draw bezier pathes. Each path will be descritized and ploted with a
separate `lines` plot. Scalar attributes will be used for all subplots. If you
provide vector attributs of same length als the pathes the i-th `lines` subplot
will see the i-th attribute.

TODO: Beziersegments plot won't work if the number of pathes changes.
"""
@recipe(BezierSegments, paths) do scene
    Attributes(default_theme(scene, LineSegments)...)
end

function Makie.plot!(p::BezierSegments)
    N = length(p[:paths][])
    PT = ptype(eltype(p[:paths][]))
    attr = p.attributes

    # set colorange automaticially if needed
    if attr[:color][] isa AbstractArray{<:Number}
        numbers = attr[:color][]
        colorrange = get(attr, :colorrange, nothing) |> to_value

        if colorrange === Makie.automatic
            attr[:colorrange] = extrema(numbers)
        end
    end

    # array which holds observables for the discretized values
    # this is needed so the `lines!` subplots will update
    disc = [Observable{Vector{PT}}() for i in 1:N]

    update_discretized!(disc, p[:paths][]) # first call to initialize
    on(p[:paths]) do paths # update if pathes change
        update_discretized!(disc, paths)
    end

    # plot all the lines
    for i in 1:N
        # each subplot will pick its specic arguments if there are vectors
        specific_attr = Attributes()

        for (k, v) in attr
            if k === :color && v[] isa AbstractVector{<:Number} && length(v[]) == N
                colormap = attr[:colormap][] |> to_colormap
                colorrange = attr[:colorrange][]
                specific_attr[k] = @lift Makie.interpolated_getindex(colormap,
                                                                     Float64($v[i]),
                                                                     colorrange)
            elseif v[] isa AbstractVector && length(v[]) == N
                specific_attr[k] = @lift $v[i]
            else
                specific_attr[k] = v
            end
        end
        lines!(p, disc[i]; specific_attr...)
    end

    return p
end

function update_discretized!(disc, pathes)
    for (i, p) in enumerate(pathes)
        disc[i][] = discretize(p)
    end
end


text:
via_df64(fn::F, data) where {F<:Function} = Float64.( fn( Double64.(data) ) )
dot_df64(fn::F, data) where {F<:Function} = Float64.( fn.( Double64.(data) ) )

via_df64(fn::F, data1, data2) where {F<:Function} = Float64.( fn(  Double64.(data1), Double64.(data2) ) )
dot_df64(fn::F, data1, data2) where {F<:Function} = Float64.( fn.( Double64.(data1), Double64.(data2) ) )

via_df64(fn1::F1, fn2::F2, data) where {F1<:Function, F2<:Function} = Float64.( fn1( fn2(  Double64.(data) ) ))
dot_df64(fn1::F1, fn2::F2, data) where {F1<:Function, F2<:Function} = Float64.( fn1.( fn2.( Double64.(data) ) ))

via_df64(fn1::F1, fn2::F2, data1, data2) where {F1<:Function, F2<:Function} = Float64.( fn1( fn2(  Double64.(data1), Double64.(data2) ) ))
dot_df64(fn1::F1, fn2::F2, data1, data2) where {F1<:Function, F2<:Function} = Float64.( fn1.( fn2.( Double64.(data1), Double64.(data2) ) ))

macro via_df64(fn, data)
    :( Float64.($(esc(fn)).(Double64.($(esc(data))) )) )
end

macro via_df64(fn, xdata, ydata)
    :( Float64.($(esc(fn)).(Double64.($(esc(xdata))), Double64.($(esc(ydata))) )) )
end

macro via_df64(fn, xdata, ydata, zdata)
    :( Float64.($(esc(fn)).(Double64.($(esc(xdata))), Double64.($(esc(ydata))), Double64.($(esc(zdata))) )) )
end

macro via_df32(fn, data)
    :( Float32.($(esc(fn)).(Double32.($(esc(data))) )) )
end

macro via_df32(fn, xdata, ydata)
    :( Float32.($(esc(fn)).(Double32.($(esc(xdata))), Double32.($(esc(ydata))) )) )
end

macro via_df32(fn, xdata, ydata, zdata)
    :( Float64.($(esc(fn)).(Double32.($(esc(xdata))), Double32.($(esc(ydata))), Double32.($(esc(zdata))) )) )
end


text:
__precompile__(true)

module Bicor

    export prepareColBicor;

    include("colwise.jl")

end # module


text:
### CONCRETE TYPE: DIRECT PROX EVALUATION
# prox! is computed using a Cholesky factorization of A'A + I/(lambda*gamma)
# or AA' + I/(lambda*gamma), according to which matrix is smaller.
# The factorization is cached and recomputed whenever gamma changes

using LinearAlgebra
using SparseArrays
using SuiteSparse

mutable struct LeastSquaresDirect{N, R <: Real, C <: RealOrComplex{R}, M <: AbstractMatrix{C}, V <: AbstractArray{C, N}, F <: Factorization} <: LeastSquares
    A::M # m-by-n
    b::V # m (by-p)
    lambda::R
    lambdaAtb::V
    gamma::Union{Nothing, R}
    shape::Symbol
    S::M
    res::Array{C, N} # m (by-p)
    q::Array{C, N} # n (by-p)
    fact::F
    function LeastSquaresDirect{N, R, C, M, V, F}(A::M, b::V, lambda::R) where {N, R <: Real, C <: RealOrComplex{R}, M <: AbstractMatrix{C}, V <: AbstractArray{C, N}, F <: Factorization}
        if size(A, 1) != size(b, 1)
            error("A and b have incompatible dimensions")
        end
        m, n = size(A)
        if m >= n
            S = lambda * (A' * A)
            shape = :Tall
        else
            S = lambda * (A * A')
            shape = :Fat
        end
        x_shape = infer_shape_of_x(A, b)
        new(A, b, lambda, lambda*(A'*b), nothing, shape, S, zero(b), zeros(C, x_shape))
    end
end

is_convex(f::LeastSquaresDirect) = f.lambda >= 0
is_concave(f::LeastSquaresDirect) = f.lambda <= 0

function LeastSquaresDirect(A::M, b::V, lambda::R) where {N, R <: Real, C <: Union{R, Complex{R}}, M <: DenseMatrix{C}, V <: AbstractArray{C, N}}
    LeastSquaresDirect{N, R, C, M, V, Cholesky{C, M}}(A, b, lambda)
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {N, R <: Real, C <: Union{R, Complex{R}}, I <: Integer, M <: SparseMatrixCSC{C, I}, V <: AbstractArray{C, N}}
    LeastSquaresDirect{N, R, C, M, V, SuiteSparse.CHOLMOD.Factor{C}}(A, b, lambda)
end

# Adjoint/Transpose versions
function LeastSquaresDirect(A::M, b::V, lambda::R) where {N, R <: Real, C <: Union{R, Complex{R}}, M <: TransposeOrAdjoint{<:DenseMatrix{C}}, V <: AbstractArray{C, N}}
    LeastSquaresDirect(copy(A), b, lambda)
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {N, R <: Real, C <: Union{R, Complex{R}}, I <: Integer, M <: TransposeOrAdjoint{<:SparseMatrixCSC{C, I}}, V <: AbstractArray{C, N}}
    LeastSquaresDirect(copy(A), b, lambda)
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {N, R <: Real, C <: Union{R, Complex{R}}, M <: AbstractMatrix{C}, V <: AbstractArray{C, N}}
    @warn "Could not infer type of Factorization for $M in LeastSquaresDirect, this type will be type-unstable"
    LeastSquaresDirect{N, R, C, M, V, Factorization}(A, b, lambda)
end

function (f::LeastSquaresDirect{N, R, C})(x::AbstractArray{C, N}) where {N, R, C}
    mul!(f.res, f.A, x)
    f.res .-= f.b
    return (f.lambda / 2) * norm(f.res, 2)^2
end

function prox!(y::AbstractArray{C, N}, f::LeastSquaresDirect{N, R, C, M, V, F}, x::AbstractArray{C, N}, gamma::R=R(1)) where {N, R, C, M, V, F}
    # if gamma different from f.gamma then call factor_step!
    if gamma != f.gamma
        factor_step!(f, gamma)
    end
    solve_step!(y, f, x, gamma)
    mul!(f.res, f.A, y)
    f.res .-= f.b
    return (f.lambda/2)*norm(f.res, 2)^2
end

function factor_step!(f::LeastSquaresDirect{N, R, C, M, V, F}, gamma::R) where {N, R, C, M, V, F}
    f.fact = cholesky(f.S + I/gamma)
    f.gamma = gamma
end

function factor_step!(f::LeastSquaresDirect{N, R, C, M, V, F}, gamma::R) where {N, R, C, M <: SparseMatrixCSC, V, F}
    f.fact = ldlt(f.S; shift = R(1)/gamma)
    f.gamma = gamma
end

function solve_step!(y::AbstractArray{C, N}, f::LeastSquaresDirect{N, R, C, M, V, F}, x::AbstractArray{C, N}, gamma::R) where {N, R, C, M, V, F <: Cholesky{C, M}}
    f.q .= f.lambdaAtb .+ x./gamma
    # two cases: (1) tall A, (2) fat A
    if f.shape == :Tall
        # y .= f.fact\f.q
        y .= f.q
        LAPACK.trtrs!('U', 'C', 'N', f.fact.factors, y)
        LAPACK.trtrs!('U', 'N', 'N', f.fact.factors, y)
    else # f.shape == :Fat
        # y .= gamma*(f.q - lambda*(f.A'*(f.fact\(f.A*f.q))))
        mul!(f.res, f.A, f.q)
        LAPACK.trtrs!('U', 'C', 'N', f.fact.factors, f.res)
        LAPACK.trtrs!('U', 'N', 'N', f.fact.factors, f.res)
        mul!(y, adjoint(f.A), f.res)
        y .*= -f.lambda
        y .+= f.q
        y .*= gamma
    end
end

function solve_step!(y::AbstractArray{C, N}, f::LeastSquaresDirect{N, R, C, M, V, F}, x::AbstractArray{C, N}, gamma::R) where {N, R, C, M, V, F}
    f.q .= f.lambdaAtb .+ x./gamma
    # two cases: (1) tall A, (2) fat A
    if f.shape == :Tall
        y .= f.fact\f.q
    else # f.shape == :Fat
        # y .= gamma*(f.q - lambda*(f.A'*(f.fact\(f.A*f.q))))
        mul!(f.res, f.A, f.q)
        f.res .= f.fact\f.res
        mul!(y, adjoint(f.A), f.res)
        y .*= -f.lambda
        y .+= f.q
        y .*= gamma
    end
end

function gradient!(y::AbstractArray{C, N}, f::LeastSquaresDirect{N, R, C, M, V, F}, x::AbstractArray{C, N}) where {N, R, C, M, V, F}
    mul!(f.res, f.A, x)
    f.res .-= f.b
    mul!(y, adjoint(f.A), f.res)
    y .*= f.lambda
    fy = (f.lambda/2)*dot(f.res, f.res)
end

function prox_naive(f::LeastSquaresDirect{N, R, C}, x::AbstractArray{C, N}, gamma::R=R(1)) where {N, R, C <: RealOrComplex{R}}
    lamgam = f.lambda*gamma
    y = (f.A'*f.A + I/lamgam)\(f.A' * f.b + x/lamgam)
    fy = (f.lambda/2)*norm(f.A*y-f.b)^2
    return y, fy
end


text:
@testset "detectors" begin
  Random.seed!(42)
  x = vcat(rand(8,2), 5*rand(1,2), 30*rand(1,2))
  x1 = vcat(rand(20,3), 5*rand(1,3), 30*rand(1,3))
  @test rxdetect(x, 0.8) == [false, false, false, false, false, false, false, false, true, true]
  @test hosvdc4detect(x, 4., 2; b=2) == [false, false, false, false, false, false, false, false, true, true]
  @test hosvdc4detect(x1, 3.9, 1; b=1) == vcat(fill(false, 20), [true, true])
  x2 = [[0. 0.]; [1. 2.]; [1. 1.]]
  @test hosvdc4detect(x2, .1, 2; b=2) == [false, false, false]
  @test rxdetect(x1, 0.9) == vcat(fill(false, 20), [true, true])
  ls = fill(true, 10)
  ls1 = [true, true, true, true, true, true, true, true, false, false]
  c = cumulants(x, 4)
  @test hosvdstep(x, ls, 4., 2, c[4])[1] == ls1
  @test hosvdstep(x, ls, 3., 1, c[4])[2] ≈ 1.15 atol= 0.1
end


text:
using Distributions
using LinearAlgebra

# Scene loader & Predictor parameters
scene_mode = "data";                                                                # "synthetic" or "data"
prediction_mode = "trajectron";                                                     # "gaussian" or "trajectron" or "oracle"
conf_file_name = "config.json"                                                      # trajectron config file name
test_data_name = "eth_test.pkl";                                                    # test data set name
test_scene_id = 0;                                                                  # test data id
start_time_idx = 905;                                                               # start time index in test data
incl_robot_node = false;                                                            # if true, robot node is created in trajectron
prediction_device = "cpu";                                                          # "cpu" or "cuda"
prediction_steps = 12;                                                              # number of steps to look ahead in the future
prediction_rng_seed = 1;                                                            # random seed for prediction (and stochastic transition for "synthetic" scenes)
deterministic = false;                                                              # if true, a single, deterministic sample is drawn regardless of random seed. (num_samples = 1 is needed)
num_samples = 30;                                                                   # number of trajectory samples (per ado agent)
use_robot_future = false;                                                           # if true, prediction is conditioned on robot's nominal control
# Cost Parameters
include("params_cost.jl")
# Control Parameters
include("params_control.jl")
# Ego initial state
ego_pos_init_vec = [5., 0.5] .+ [-5.263534, -5.314636];                             # initial ego position [x, y] [m]
ego_pos_goal_vec = [5., 8.9] .+ [-5.263534, -5.314636];                             # goal ego position [x, y] [m]
# Other parameters
pos_error_replan = 2.0;                                                             # position error for replanning target trajectory [m]
target_speed = 0.7;                                                                 # target speed [m/s]
sim_horizon = 10.0;                                                                 # simulation time horizon [s]


text:
__precompile__()

module Andes

using Conda
using PyCall

import SparseArrays

const py = PyNULL()

include("kvxopt.jl")

function __init__()
    copy!(py, pyimport_conda("andes", "andes", "conda-forge"))

    pytype_mapping(pyimport("kvxopt").spmatrix, SparseArrays.SparseMatrixCSC)
end

# --- export ---
export convert;

end # module