text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

module ForestFireExample

export run_forest_fire_model!

using Agents
using Random
using Plots
using InteractiveDynamics
using CairoMakie

mutable struct Tree <: AbstractAgent
    id::Int
    pos::Dims{2}
    status::Symbol
end

function forest_fire(;
    density = 0.7,
    griddims = (100, 100)
)
    space = GridSpace(griddims; periodic = false, metric = :euclidean)
    forest = AgentBasedModel(Tree, space)
    # density is a probabilty of whether to add a tree
    for position in positions(forest)
        if rand() < density
            state = position[1] == 1 ? :burning : :green
            add_agent!(position, forest, state)
        end
    end
    return forest
end

function tree_step!(tree, forest)
    # the current tree is burning
    if tree.status == :burning
        # Find all green
        for neighbor in nearby_agents(tree, forest)
            if neighbor.status == :green
                neighbor.status = :burning
            end
        end
        tree.status = :burnt
    end
    return nothing
end

function burnt_percentage(m)
    burntTrees = count(t -> t.status == :burnt, allagents(m))
    totArea = length(positions(m))
    return burntTrees / totArea
end

function treecolor(a)
    if a.status == :burning
        return :red
    elseif a.status == :burnt
        return :darkred
    else
        return :green
    end
end

function run_forest_fire_model_old!()
    Random.seed!(2)
    forest = forest_fire()
    mdata = [burnt_percentage]
    df = init_model_dataframe(forest,mdata)
    anim = @animate for i in 0:50
        if i > 0
            step!(forest, tree_step!, 1)
        end
        collect_model_data!(df,forest,mdata,i)
        p1 = plotabm(forest; ac = treecolor, ms = 5, msw = 0)
        title!(p1,"step $(i)")
    end
    animated_gif = gif(anim, "forest.gif", fps=2)
    return (animated_gif, df)
end

function run_forest_fire_model!(file;frames=60,spf=1)
    Random.seed!(2)
    forest = forest_fire()
    mdata = [burnt_percentage]
    df = init_model_dataframe(forest, mdata)
    # Create 'Step N' as the Observable title
    s = Observable(0) # counter of current step
    t = lift(x -> "step = " * string(x), s)
    axiskwargs = (title = t, titlealign = :left) #title and position
    # create figure of the first frame to feed it into record
    fig, abmstepper = abm_plot(
        forest;
        ac = treecolor,
        axiskwargs = axiskwargs,
        ms = 5,
        msw = 0
    )
    # record the animation Makie style
    record(fig, file) do io
        for j in 1:frames-1
            recordframe!(io)
            # MUST USE THE SAME ABMSTEPPER! otherwise frame won't update
            step!(abmstepper, forest, tree_step!, dummystep, spf)
            collect_model_data!(df,forest,mdata,j)
            s[] += spf # update the observable step
            #s[] = s[] # Not sure why?
        end
        recordframe!(io) # last frame
    end
    isfile(file) || error("file not created")
    return (file,df)
end

#_,df = run_forest_fire_model!("forest_new.mp4")
end

ext: jl


text:
using ArgParse
using Random
using CUDAnative
using CuArrays
using Flux
using FileIO

include("DecisionStrain.jl")
using .DecisionStrain
using .DecisionStrain.Data
using .DecisionStrain.FluxFixes

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--gpu"
            help = "Select GPU to use."
            arg_type = Int
            default = 0
        "--dat", "-i"
            help = "path to directory with data set"
            arg_type = String
            default = "dat/FB15K"
        "--binary_dat"
            help = "Data set is stored in binary files."
            action = :store_true
        "--model"
            help = "select model. must be either DecisionStrain or ComplEx"
            arg_type = String
            default = "DecisionStrain"
        "--head_prediction"
            help = "Only used for --model DecisionStrain: train head prediction rather than tail prediction"
            action = :store_true
        "--aux"
            help = "path to auxiliary model"
            arg_type = String
        "--epochs", "-E"
            help = "number of training epochs"
            arg_type = Int
            default = 100
        "--batchsize", "-B"
            help = "minibatch size"
            arg_type = Int
            default = 500
        "--eval_batchsize"
            help = "minibatch size for evaluation"
            arg_type = Int
            default = 500
        "--hidden_dim", "-k"
            help = "dimension of embedding space"
            arg_type = Int
            default = 20
        "--lmbda"
            help = "prefactor of regularizer"
            arg_type = Float64
            default = 0.01
        "--lr"
            help = "learning rate"
            arg_type = Float64
            default = 1.0
        "--neg_ratio"
            help = "Number of negative samples per positive sample"
            arg_type = Int
            default = 10
        "--rngseed"
            help = """
                seed for random number generator; if not provided, a seed will
                automatically be generated from system entropy and written to
                the output file."""
            arg_type = UInt32
        "--population", "-n"
            help = "number of simultaneously trained decision trees"
            arg_type = Int
            default = 1
        "--out", "-o"
            help = "output directory; must not exist"
            arg_type = String
        "--continue", "-c"
            help = "provide file name to continue training from."
            arg_type = String
        "--eval_every"
            help = "evaluate model every n epochs on validation set (set to zero to evaluate only at end)"
            arg_type = Int
            default = 10
        "--save_every"
            help = "save checkpoing every n epochs (set to zero to save only at end)"
            arg_type = Int
            default = 100
        "--copy_subtrees"
            help = "copy branching vectors of winning subtrees"
            action = :store_true
        "--precision"
            help = """
                floating point precision;  allowed values are "half", "single",
                and "double".  Default is "single"."""
            arg_type = Union{Type{Float16}, Type{Float32}, Type{Float64}}
            default = Float32
    end

    return parse_args(s, as_symbols = true)
end

function ArgParse.parse_item(
        ::Type{Union{Type{Float16}, Type{Float32}, Type{Float64}}},
        s::AbstractString)
    if s == "half"
        return Float16
    elseif s == "single"
        return Float32
    elseif s == "double"
        return Float64
    else
        throw("""Unknown float type "$s".""")
    end
end

function dicttostr(dict::Dict)
    stream = IOBuffer()
    show(stream, MIME("text/plain"), dict)
    lines = split(String(take!(stream)), "\n")
    "Dict(\n" * join("$line,\n" for line in lines[2:end]) * ")"
end

function main()
    args = parse_commandline()

    device!(args[:gpu])

    # initialize random number generator
    if args[:rngseed] === nothing
        args[:rngseed] = rand(RandomDevice(), UInt32)
    end
    rng = MersenneTwister(args[:rngseed]);

    if args[:out] === nothing
        args[:out] = joinpath("out", repr(hash(args))[end - 5 : end])
    end

    # create output directory and open log file
    print("Output directory: $(args[:out])\n")
    mkdir(args[:out])
    try
        open(joinpath(args[:out], "log"), write = true) do log
            # TODO: print date, time, hostname, git commit, and $CUDA_VISIBLE_DEVICES
            print(log, """
            args = $(dicttostr(args))

            # Running on $(Threads.nthreads()) threads.

            """)

            dat = Dataset(args[:dat], binary_files = args[:binary_dat], log = log,
                          reverse = args[:head_prediction])
            if args[:continue] === nothing
                if args[:model] == "DecisionStrain"
                    model = DecisionStrainModel(rng, dat, args, log = log)
                elseif args[:model] == "ComplEx"
                    aux_model = nothing
                    if args[:aux] !== nothing
                        aux_model = gpu(load(args[:aux], "model"), args, aux_model=true)
                    end
                    model = ComplExModel(rng, dat, args, aux_model)
                else
                    throw("Unknown model $(args[:model]).")
                end
                print(log, "# Initialized new model.\n\n")
            else
                model = gpu(load(args[:continue], "model"), args)
                print(log, "# Loaded model from file $(args[:continue])\n\n")
            end
            flush(log)

            # Run the training loop in a separate function to allow the compiler to generate
            # specialized code based on floattype(model).
            training_loop!(
                model, dat, rng, args[:epochs], args[:batchsize], args[:eval_every],
                args[:save_every], args[:out], log = log)
        end
    catch e
        open(joinpath(args[:out], "err"), write = true) do err
            trace = catch_backtrace()
            showerror(err, e, trace)
            print(err, "\n")
            rethrow()
        end
    end
end

function training_loop!(
        model::AbstractKnowledgGraphModel,
        dat::Dataset,
        rng::AbstractRNG,
        nepochs::Integer,
        batchsize::Integer,
        eval_every::Integer,
        save_every::Integer,
        out_path::AbstractString;
        log::IO = stdout)

    trainlength = size(dat.train_unsorted, 2)
    num_batches = trainlength ÷ batchsize

    # print(log, "initial_performance = ", evaluate(model, dat.val, dat), "\n\n")
    print(log, """
    col_headers = ["epoch", "avg_duration_continuous", "avg_duration_discrete",
        "mrr", "hits_at", "log_likelihood_per_data_point", "eval_fail_cnt"]

    data = [
    """)
    flush(log)

    sum_duration_cont = 0.0
    sum_duration_disc = 0.0
    duration_disc = 0.0
    last_print = 0

    for epoch in 1 : nepochs
        # print_stats(model)
        _, duration_cont = @timed begin
            shuffled = CuArray(shuffle_train(rng, dat))
            for batchstart = 1 : batchsize : batchsize * num_batches
                train_minibatch!(
                    rng, model, shuffled, batchstart, batchsize, epoch, log = log)
            end
        @sync_cuda("continuous updates done")
        end
        sum_duration_cont += duration_cont

        if needs_discrete_optimization(model)
            _, duration_disc = @timed begin
                train_epoch!(rng, model, dat, epoch, batchsize, log = log)
                @sync_cuda("discrete updates done")
            end
            sum_duration_disc += duration_disc
        else
            duration_disc = 0.0
        end

        GC.gc()

        print(log, "# Epoch $epoch took $duration_cont + $duration_disc seconds.\n")
        flush(log)

        if save_every != 0 && epoch % save_every == 0
            save_checkpoint(model, out_path, epoch)
            GC.gc()
        end

        # if eval_every != 0 && epoch % eval_every == 0
        #     print(log, (
        #         epoch,
        #         sum_duration_cont / (epoch - last_print),
        #         sum_duration_disc / (epoch - last_print),
        #         evaluate(model, dat.val, dat)...
        #         ), ",\n")
        #     flush(log)
        #     sum_duration_cont = sum_duration_disc = 0.0
        #     last_print = epoch
        #     GC.gc()
        # end
    end

    if save_every == 0 || nepochs % save_every != 0
        save_checkpoint(model, out_path, nepochs)
        GC.gc()
    end

    # if eval_every == 0 || nepochs % eval_every != 0
    #     print(
    #         log,
    #         (nepochs, nothing, nothing, evaluate(model, dat.val, dat)...),
    #         ",\n"
    #     )
    # end
    print(log, "]\n")
    flush(log)
end

@noinline function save_checkpoint(model, directory, epoch)
    save(joinpath(directory, "checkpoint-$epoch.jld2"), "model", cpu(model))
end

@noinline function evaluate(
            model::AbstractKnowledgGraphModel,
            heldout::AbstractArray{<:Integer, 2},
            fulldat::Dataset;
            batchsize = 100,
            windows=[1, 3, 10])
    hits_cnt_at_t = zeros(Int, length(windows), Threads.nthreads())
    hits_cnt_at_h = zeros(Int, length(windows), Threads.nthreads())
    sum_rr_t = zeros(Float64, Threads.nthreads())
    sum_rr_h = zeros(Float64, Threads.nthreads())
    cnt = zeros(Int, Threads.nthreads())
    total_ll = 0.0

    for batchstart = 1 : batchsize : size(heldout, 2)
        batchend = min(batchstart + batchsize - 1, size(heldout, 2))
        batch = batchstart : batchend
        hs = heldout[1, batch]
        ts = heldout[2, batch]
        rs = heldout[3, batch]
        all_ll, predicted_lls, raw_ranks, sum_ll = log_likelihoods_and_raw_ranks(model, hs, ts, rs)
        total_ll += sum_ll

        Threads.@threads for i = 1 : length(batch)
            predicted_ll = predicted_lls[i]
            raw_rank = raw_ranks[i]
            other_correct_ts = view(fulldat.all_sorted, fulldat.all_ranges[(hs[i], rs[i])])
            other_correct_lls = all_ll[other_correct_ts, i]
            filtered_rank_t = raw_rank - sum(other_correct_lls .> predicted_ll)

            threadid = Threads.threadid()

            for (ii, w) = enumerate(windows)
                if filtered_rank_t <= w
                    hits_cnt_at_t[ii, threadid] += 1
                end
            end

            sum_rr_t[threadid] += 1 / filtered_rank_t
            cnt[threadid] += 1
        end
    end

    success_cnt = sum(cnt)
    hits_cnt_at_t = sum(hits_cnt_at_t, dims = 2)[:, 1]
    hits_at_t = Dict(w => hits / success_cnt for (w, hits) = zip(windows, hits_cnt_at_t))
    mrr_t = sum(sum_rr_t) / success_cnt
    return mrr_t, hits_at_t, total_ll / size(heldout, 2), size(heldout, 2) - success_cnt
end

main()

max_stars_repo_name: mandt-lab/adversarial-negative-sampling

max_issues_repo_name: mandt-lab/adversarial-negative-sampling

max_forks_repo_name: mandt-lab/adversarial-negative-sampling


text:
"""
```
label_components(OneComponent3D(), binary_volume, local_window)
```

Computes the connected components of a binary volume using 26-connectivity.

# Output

Returns a three-dimensional `Int` array that represents the label of every voxel.

# Arguments

## `binary_volume::AbstractArray`

The binary volume made up of values 0 (black) and 1 (white).

## `local_window::NTuple{3, Int}`

The size of the local window for each iteration.

# Example

Computes the labels for a randomly generated 20 x 20 x 20 array using 5 x 7 x 7 local window.

```julia
using ImageComponentAnalysis

binary_volume = Int.(rand(Bool, 20, 20, 20))

labels = label_components(OneComponent3D(), binary_volume, (5, 7, 7))
```

# Reference

[1] Q. Hu, G. Qian, and W.L. Nowinski, “Fast connected-component labelling in three-dimensional binary images based on iterative recursion”, Computer Vision and Image Understanding, vol. 99, no. 3, pp. 414–434, Sep. 2005. [doi:10.1016/j.cviu.2005.04.001](https://doi.org/10.1016/j.cviu.2005.04.001)
"""
function label_components(algorithm::OneComponent3D, binary_volume::AbstractArray, local_window::NTuple{3, Int})

	@inbounds begin
		height, width, length = size(binary_volume)
		queue = Queue{NTuple{3, Int}}()

		n₁, n₂, n₃ = map((b, l) -> ceil(Int, (max(1, min(b, l)-1)/2)), size(binary_volume), local_window)
		labels = padarray(Int, binary_volume, Fill(0, (n₁ + 1, n₂ + 1, n₃ + 1)))
		labelindex = 3

		for x = 1:height, y = 1:width, z = 1:length
			if (labels[x, y, z] == 1)
				enqueue!(queue, (x, y, z))
				while !isempty(queue)
					vᵢx, vᵢy, vᵢz = dequeue!(queue)
					label_subvolume!(algorithm, vᵢx, vᵢy, vᵢz, n₁, n₂, n₃, queue, labels, labelindex)
				end
				labelindex += 1
			end
		end

		labels = view(labels, 1:height, 1:width, 1:length)
		replace!(s -> s >= 3 ? s-2 : 0, labels)
	end
end

# Label 26-connected voxels in the subvolume with the labelindex.
function label_subvolume!(algorithm, vᵢx, vᵢy, vᵢz, n₁, n₂, n₃, queue, labels, labelindex)

	window_size = 0
	center = (vᵢx, vᵢy, vᵢz)

	# Set center voxel of subvolume as labelindex.
	labels[vᵢx, vᵢy, vᵢz] = labelindex

	# Label 3 x 3 x 3 neighbourhood of center voxel.
	for (x, y, z) = get_neighbours(algorithm, 1, center, (n₁, n₂, n₃))
		if (labels[x, y, z] == 1)
			labels[x, y, z] = -10
		end
	end

	# Grow window size to local window.
	while window_size != max(n₁, n₂, n₃) + 1
		for (vx, vy, vz) = get_neighbours(algorithm, window_size, center, (n₁, n₂, n₃))

			# Label object voxels with distance of window_size.
			if (labels[vx, vy, vz] == -10)
				if (abs(vx - vᵢx) == n₁ || abs(vy - vᵢy) == n₂ || abs(vz - vᵢz) == n₃)
					enqueue!(queue, (vx, vy, vz))
				end

				# Label 3 x 3 x 3 neighbourhood of current object voxel.
				for (x, y, z) = get_neighbours(algorithm, 1, (vx, vy, vz), (n₁, n₂, n₃), center)
					if (labels[x, y, z] == 1)
						labels[x, y, z] = -10
						if (abs(x - vᵢx) == n₁ || abs(y - vᵢy) == n₂ || abs(z - vᵢz) == n₃)
							enqueue!(queue, (x, y, z))
						end
					end
				end

				labels[vx, vy, vz] = labelindex
			end
		end
		window_size += 1
	end

	# Replace and queue all object voxels with label -10.
	for (x, y, z) = Iterators.product(vᵢx-n₁:vᵢx+n₁, vᵢy-n₂:vᵢy+n₂, vᵢz-n₃:vᵢz+n₃)
		if labels[x, y, z] == -10
			enqueue!(queue, (x, y, z))
			labels[x, y, z] = labelindex
		end
	end
end

# Get the neighbourhood of the center voxel.
function get_neighbours(algorithm::OneComponent3D, distance, center, n, truecenter = center)

	neighbour_ranges = map((c, n, t) -> max(t-n,c-distance):min(t+n,c+distance), center, n, truecenter)

	function select_neighbours((x, y, z))
		max(abs(x-center[1]), abs(y-center[2]), abs(z-center[3])) == distance
	end

	Iterators.filter(select_neighbours, Iterators.product(neighbour_ranges...))
end


text:
# Autogenerated wrapper script for Qt5Svg_jll for x86_64-apple-darwin
export libqt5svg

using Qt5Base_jll
JLLWrappers.@generate_wrapper_header("Qt5Svg")
JLLWrappers.@declare_library_product(libqt5svg, "@rpath/QtSvg")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Base_jll)
    JLLWrappers.@init_library_product(
        libqt5svg,
        "lib/QtSvg.framework/Versions/5/QtSvg",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
"""
function spg_standardize_cell(lattice::Array{Float64, 2},
                              positions::Array{Float64, 2},
                              types::Array{Int64, 1},
                              num_atom::Int64,
                              to_primitive::Int64,
                              no_idealize::Int64,
                              symprec::Float64)
    # transpose to colomn vector used by spglib
    lattice_ = copy(lattice)
    positions_ = copy(positions)

    allocN = 4
    positions_ = zeros(Float64, 3, num_atom*allocN)
    types_ = zeros(Int64, num_atom*allocN)

    positions_[:, 1:num_atom] = positions[:, :]
    types_[1:num_atom] = types

    types_ = Base.cconvert(Array{Int32, 1}, types_)
    num_atom = Base.cconvert(Int32, num_atom)
    to_primitive = Base.cconvert(Int32, to_primitive)
    no_idealize = Base.cconvert(Int32, no_idealize)

    num_primitive_atom =
    ccall( (:spg_standardize_cell, libsymspg), Int32,
           ( Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Int32, Int32, Int32, Float64 ),
           lattice_, positions_, types_, num_atom, to_primitive, no_idealize, symprec )

    positions = positions_[:, 1:num_primitive_atom]
    types = types_[1:num_primitive_atom]

    return lattice_, positions, Base.cconvert(Array{Int64, 1}, types), Base.cconvert(Int64,num_primitive_atom)
end

function spg_standardize_cell(lattice::Array{Float64, 2},
                              positions::Array{Float64, 2},
                              types::Array{Int64, 1},
                              to_primitive::Bool,
                              no_idealize::Bool,
                              symprec::Float64)

    to_primitive ? to_primitive_ = 1 : to_primitive_ = 0
    no_idealize ? no_idealize_ = 1 : no_idealize_ = 0

    num_atom = size(types)[1]

    return spg_standardize_cell(lattice, positions, types, num_atom, to_primitive_, no_idealize_, symprec)
end

function spg_find_primitive(lattice::Array{Float64, 2},
                            positions::Array{Float64, 2},
                            types::Array{Int64, 1},
                            symprec::Float64)

    return spg_standardize_cell(lattice, positions, types, true, false, symprec)
end

function spg_refine_cell(lattice::Array{Float64, 2},
                         positions::Array{Float64, 2},
                         types::Array{Int64, 1},
                         symprec::Float64)

    return spg_standardize_cell(lattice, positions, types, false, false, symprec)
end

"""
angle_tolerance
"""
function spgat_standardize_cell(lattice::Array{Float64, 2},
                              positions::Array{Float64, 2},
                              types::Array{Int64, 1},
                              num_atom::Int64,
                              to_primitive::Int64,
                              no_idealize::Int64,
                              symprec::Float64,
                              angle_tolerance::Float64)
    # transpose to colomn vector used by spglib
    lattice_ = copy(lattice)
    positions_ = copy(positions)

    allocN = 4
    positions_ = zeros(Float64, 3, num_atom*allocN)
    types_ = zeros(Int64, num_atom*allocN)

    positions_[:, 1:num_atom] = positions[:, :]
    types_[1:num_atom] = types

    types_ = Base.cconvert(Array{Int32, 1}, types_)
    num_atom = Base.cconvert(Int32, num_atom)
    to_primitive = Base.cconvert(Int32, to_primitive)
    no_idealize = Base.cconvert(Int32, no_idealize)

    num_primitive_atom =
    ccall( (:spgat_standardize_cell, libsymspg), Int32,
           ( Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Int32, Int32, Int32, Float64, Float64 ),
           lattice_, positions_, types_, num_atom, to_primitive, no_idealize, symprec, angle_tolerance )

    positions = positions_[:, 1:num_primitive_atom]
    types = types_[1:num_primitive_atom]

    return lattice_, positions, Base.cconvert(Array{Int64, 1}, types), Base.cconvert(Int64,num_primitive_atom)
end

function spgat_standardize_cell(lattice::Array{Float64, 2},
                              positions::Array{Float64, 2},
                              types::Array{Int64, 1},
                              to_primitive::Bool,
                              no_idealize::Bool,
                              symprec::Float64,
                              angle_tolerance::Float64)

    to_primitive ? to_primitive_ = 1 : to_primitive_ = 0
    no_idealize ? no_idealize_ = 1 : no_idealize_ = 0

    num_atom = size(types)[1]

    return spgat_standardize_cell(lattice, positions, types, num_atom, to_primitive_, no_idealize_, symprec, angle_tolerance)
end

function spgat_find_primitive(lattice::Array{Float64, 2},
                            positions::Array{Float64, 2},
                            types::Array{Int64, 1},
                            symprec::Float64,
                            angle_tolerance::Float64)

    return spgat_standardize_cell(lattice, positions, types, true, false, symprec, angle_tolerance)
end

function spgat_refine_cell(lattice::Array{Float64, 2},
                         positions::Array{Float64, 2},
                         types::Array{Int64, 1},
                         symprec::Float64,
                         angle_tolerance::Float64)

    return spgat_standardize_cell(lattice, positions, types, false, false, symprec, angle_tolerance)
end


text:
using FileIO,JLD2,MAT,YAML,Mmap,LoopVectorization
import FileIO: save

include("SpikeGLX.jl")
include("Imager.jl")

"""
Drop array dims according to `MATLAB` convention.

- scalar: scalar instead of 1x1 array
- rvector: vector instead of 1xN array
- cvector: vector instead of Nx1 array
"""
dropmatdim!(x;scalar = true, rvector = true, cvector = true) = x
function dropmatdim!(x::Dict;scalar = true, rvector = true, cvector = true)
    foreach(k->x[k]=dropmatdim!(x[k],scalar=scalar,rvector=rvector,cvector=cvector),keys(x))
    return x
 end
function dropmatdim!(x::Array;scalar = true, rvector = true, cvector = true)
    if ndims(x) == 2
        nr,nc = size(x)
        if nr==1 && nc==1 && scalar
            return dropmatdim!(x[1,1],scalar=scalar,rvector=rvector,cvector=cvector)
        elseif nr==1 && nc>1 && rvector
            x = dropdims(x,dims=1)
        elseif nr>1 && nc==1 && cvector
            x = dropdims(x,dims=2)
        end
    end
    if x isa Array{Any}
        foreach(i->x[i]=dropmatdim!(x[i],scalar=scalar,rvector=rvector,cvector=cvector),eachindex(x))
    end
    return x
end

"""
Read variables of a `MATLAB` MAT file into a Dictionary

1. f: MAT file path
2. vars: variable names in the varset to read

- varset: variable names to select in vars
- scalar: scalar instead of 1x1 matrix
- rvector: vector instead of 1xN matrix
- cvector: vector instead of Nx1 matrix
"""
function readmat(f::AbstractString,vars...;varset=["spike","lfp","digital","analog1k","image"],scalar = true, rvector = true, cvector = true)
    if isempty(vars)
        d=matread(f)
    else
        d=Dict()
        matopen(f) do file
            fvs = names(file)
            nonselectvars=setdiff(fvs,varset)
            selectvars = intersect(varset,vars)
            for i in setdiff(selectvars,fvs)
                @warn """Selected variable "$i" not found in "$f"."""
            end
            readvars=union(nonselectvars,intersect(fvs,selectvars))
            if !isempty(readvars)
                for v in readvars
                    d[v]=read(file,v)
                end
            end
        end
    end
    if scalar || rvector || cvector
        d = dropmatdim!(d,scalar=scalar,rvector=rvector,cvector=cvector)
    end
    return d
end

"Read `DataExport` Metadata MAT File"
readmeta(f::AbstractString) = DataFrame(readmat(f)["Tests"])

"Load images in directory"
function loadimageset(dir;name=[],n=typemax(Int),alpha=false)
    isdir(dir) || error("Invalid Directory")
    for (root, dirs, files) in walkdir(dir)
        if isempty(name)
            n = min(n,length(files))
            n<=0 && error("n <= 0")
            imgs=load.(joinpath.([root],string.(1:n).*splitext(files[1])[2]))
        else
            imgs=load.(joinpath.([root],string.(Int.(name)).*splitext(files[1])[2]))
        end
        if alpha
            imgs = map(i->coloralpha.(i),imgs)
        end
        return imgs
    end
end

"Default digital input channels"
const CondDCh=1
const MarkDCh=2
const StartDCh=3
const StopDCh=4
const Bits16DCh=5
"Get digital channel time and value."
function digitaldata(dataset::Dict,ch)
    chidx = findfirst(dataset["digital"]["channel"].==ch)
    if !isempty(chidx)
        return dataset["digital"]["time"][chidx],dataset["digital"]["data"][chidx]
    else
        return [],[]
    end
end
"Parse bit event time and value."
function digitalbit(dt,dv,bits...)
    n = length(dt)
    maxbit = maximum(bits)
    bn = length(bits)
    bt = [zeros(Float64,n+1) for _ in 1:bn]
    bv = [zeros(UInt8,n+1) for _ in 1:bn]
    j=ones(Int,bn)
    for i in 1:n
        ds = digits(UInt8,dv[i],2,maxbit)
        for b in 1:bn
            v=ds[bits[b]]
            if v != bv[b][j[b]]
                j[b]=j[b]+1
                bt[b][j[b]]=dt[i]
                bv[b][j[b]]=v
            end
        end
    end
    for b in 1:bn
        bt[b]=bt[b][2:j[b]]
        bv[b]=bv[b][2:j[b]]
    end
    if bn==1
        return bt[1],bv[1]
    end
    return bt,bv
end

"Read and Prepare Dataset in `MATLAB` MAT File"
prepare(f::AbstractString,vars...)=prepare!(readmat(f,vars...))
function prepare!(d::Dict)
    if haskey(d,"secondperunit")
        settimeunit(d["secondperunit"])
    end
    if haskey(d,"sourceformat")
        fmt = d["sourceformat"]
        if fmt=="Ripple"
            d=prepare_ripple!(d)
        elseif fmt=="OI"
            d=prepare_oi!(d)
        elseif fmt=="Experica"
            d=prepare_experica!(d)
        end
    end
    if haskey(d,"imecindex")
        imecspike=Dict()
        for i in d["imecindex"]
            s = "spike$(i)_kilosort"
            if haskey(d,s)
                imecspike[i] = unitspike_kilosort(d[s],syncindex=i,sortspike=true)
            end
        end
        d = mergeimecspike!(imecspike,d)
    end
    return d
end
"Prepare `Ripple` Data"
function prepare_ripple!(d::Dict)
    if haskey(d,"spike")
        st=d["spike"]["time"]
        su=d["spike"]["unitid"]
        for i in 1:length(st)
            if ndims(st[i])==0
                st[i]=[st[i]]
                su[i]=[su[i]]
            end
        end
    end
    return d
end
"Prepare `Optical Imaging` Block Data"
function prepare_oi!(d::Dict)
    return d
end
"Prepare `Experica` Experiment Data"
function prepare_experica!(d::Dict)
    return d
end

"Extract state time of condtest"
function statetime(ct::Dict;statetype::AbstractString="CONDSTATE",state::AbstractString="COND")
    if haskey(ct,statetype)
        filter!(l->!isempty(l),map(i->begin
        t = filter!(k->!isempty(k),map(j->haskey(j,state) ? j[state] : [],i))
        length(t)==1 ? t[1] : t
    end,ct[statetype]))
else
    []
end
end

function getparam(d::Dict{String,Any},name,fromobject="")
    for k in keys(d)
        startswith(k,name) && endswith(k,fromobject) && return d[k]
    end
    nothing
end

"Get `CondTestCond` DataFrame"
function condtestcond(ctcd::Dict)
    ctc = DataFrame(ctcd)
    return ctc
end
"Get `CondTest` DataFrame"
function condtest(ctd::Dict)
    ct = DataFrame(ctd)
    return ct
end
"Get `CondTest` and `CondTestCond` DataFrame"
function ctctc(ctd::Dict,ctcd::Dict)
    ct = condtest(ctd)
    ctc = condtestcond(ctcd)
    return ct,ctc
end
function ctctc(ex::Dict)
    return ctctc(ex["CondTest"],ex["CondTestCond"])
end

"Map time to data timeline, optionally add latency."
function maptodatatime(x,ex::Dict;addlatency=true)
    t=x*(1+ex["TimerDriftSpeed"])+ex["t0"]
    if addlatency
        t+=ex["Latency"]
    end
    return t
end

"Map time between ref and target, according to the timing diff between target and ref of the same sync signal"
function ref2sync(t,dataset::Dict,si="0")
    haskey(dataset,"syncdt") || return t
    tt = "syncdiff$si"
    if haskey(dataset,tt)
        ref2sync(t,dataset[tt],dataset["syncdt"])
    else
        @warn "No syncdiff for Sync$si, skip timing correction.";return t
    end
end
function sync2ref(t,dataset::Dict,si="0")
    haskey(dataset,"syncdt") || return t
    tt = "syncdiff$si"
    if haskey(dataset,tt)
        sync2ref(t,dataset[tt],dataset["syncdt"])
    else
        @warn "No syncdiff for Sync$si, skip timing correction.";return t
    end
end
ref2sync(t,syncdiff,syncdt=500) = t .+ @view syncdiff[clamp!(round.(Int,t./syncdt),1,length(syncdiff))]
sync2ref(t,syncdiff,syncdt=500) = t .- @view syncdiff[clamp!(round.(Int,t./syncdt),1,length(syncdiff))]

"Epochs of `Neuropixels` Channel Sample `x`, optionally gain corrected(voltage), line noise(60,120,180Hz) removed and bandpass filtered"
function epochsamplenp(x,fs,epochs,chs;meta=[],bandpass=[1,100])
    f = nothing
    if !isempty(meta)
        f = i -> gaincorrectnp(i,meta)
        if !isempty(bandpass)
            f = i -> hlpass(rmline!(gaincorrectnp(i,meta),fs),fs,high=bandpass[1],low=bandpass[2])
        end
    end
    epochsample(x,fs,epochs,chs,fun=f)
end

"""
Organize each spiking unit from `Kilosort` result.
"""
function unitspike_kilosort(data::Dict;syncindex="0",sortspike::Bool=true)
    # each unit
    unitid = data["clusterid"]
    unitgood = data["clustergood"].==1
    unitindex = [data["cluster"].==i for i in unitid]
    unitspike = map(i->data["time"][i],unitindex)
    unittemplate = map(i->data["template"][i],unitindex)
    unitamplitude = map(i->data["amplitude"][i],unitindex)

    chmap = data["chanmap"]
    chposition = data["channelposition"]
    unittemplatesindex = unique.(unittemplate)
    unitposition = vcat(map(i->mean(data["templatesposition"][i,:],dims=1),unittemplatesindex)...)
    unittemplateamplitude = map(i->mean(data["templatesamplitude"][i]),unittemplatesindex)
    # first found template feature as unit template feature
    unittemplatefeature = Dict(k=>data["templateswaveformfeature"][k][first(unittemplatesindex)] for k in keys(data["templateswaveformfeature"]))
    unitfeature = haskey(data,"clusterwaveformfeature") ? data["clusterwaveformfeature"] : unittemplatefeature

    sortspike && foreach(sort!,unitspike)
    return Dict("unitid"=>unitid,"unitgood"=>unitgood,"unitspike"=>unitspike,"chposition"=>chposition,"unitposition"=>unitposition,
                "unittemplateamplitude"=>unittemplateamplitude,"unitfeature"=>unitfeature,"isspikesorted"=>sortspike,"unitsync"=>fill(syncindex,size(unitspike)))
end

function mergeimecspike!(imecspike,d)
    n = length(imecspike)
    if n == 0
        @warn "No IMEC spike data"
    elseif n==1
        d["spike"] = first(values(imecspike))
    else
        @info "Merging spike from multiple IMEC, not implemented"
    end
    return d
end

"convert `Matlab` struct of array to `DataFrame`"
matdictarray2df(d) = DataFrame(Any[squeeze(v,2) for v in values(d)],[Symbol(k) for k in keys(d)])

"Regular expression to match `VLab` data file names"
function vlabregex(testtype;subject="[0-9]",maxch=1,cell="[a-z]",maxrepeat=3,format="mat")
    mr = lpad(maxrepeat,2,0)
    Regex("^$subject+[A-Za-z]+[1-$maxch]$cell$testtype[0-$(mr[1])][0-$(mr[2])]\\.$format")
end

"Get matched file names in path"
function getdatafile(testtype;subject="[0-9]",path="./data")
    matchfile(vlabregex(testtype,subject=subject),path=path)
end

"Regular Expression to match Optical Imaging `VDAQ` block file name"
function oifileregex(;subject="[A-Za-z0-9]",session="[A-Za-z0-9]",experimentid="[0-9]",format="mat")
    Regex("^$subject+_$session*_?E$experimentid+B[0-9]+[.]$format")
end

"Get Optical Imaging `VDAQ` matched files in directory"
function getoifile(;subject="[A-Za-z0-9]",session="[A-Za-z0-9]",experimentid="[0-9]",format="mat",dir="",adddir=false)
    matchfile(oifileregex(subject=subject,session=session,experimentid=experimentid,format=format),dir=dir,adddir=adddir)
end

"Regular Expression to match `Experica` data file name"
function expericafileregex(;subject="[A-Za-z0-9]",session="[A-Za-z0-9]",site="[A-Za-z0-9]",test="[A-Za-z0-9]",maxrepeat=10,format="mat")
    mr = lpad(maxrepeat,2,"0");d2=mr[1];d1=parse(Int,d2)>0 ? 9 : mr[2]
    Regex("^$subject+_$session*_?$site*_?$test+_[0-$d2]?[0-$d1][.]$format")
end

"Get matched `Experica` files in directory"
function getexpericafile(;subject="[A-Za-z0-9]",session="[A-Za-z0-9]",site="[A-Za-z0-9]",test="[A-Za-z0-9]",maxrepeat=10,format="mat",dir="",adddir=false)
    matchfile(expericafileregex(subject=subject,session=session,site=site,test=test,maxrepeat=maxrepeat,format=format),dir=dir,adddir=adddir)
end

"Get matched file names in directory, optionally join directory to get file path"
function matchfile(pattern::Regex;dir="",join::Bool=false)
    fs = filter!(f->occursin(pattern,f),readdir(dir))
    if join
        fs=joinpath.(dir,fs)
    end
    return fs
end

function save(filepath::Union{AbstractString,IO},cm::Dict)
    ext = splitext(filepath)[2]
    if ext == ".yaml"
        data = deepcopy(cm)
        foreach(n->data[n]=Dict(k=>cm isa String ? cm : vec.(cm) for (k,cm) in pairs(data[n])),keys(data))
        YAML.write_file(filepath,data)
    elseif ext == ".mat"
        data = deepcopy(cm)
        foreach(n->data[n]=Dict(k=>cm isa String ? cm : mapreduce(c->vec(c)',vcat,cm) for (k,cm) in pairs(data[n])),keys(data))
        matwrite(filepath,data)
    end
end


text:
const TYPE_CACHE = Dict{String,CLRObject}()

function gettypecached(typename)
    if haskey(TYPE_CACHE, typename)
        return TYPE_CACHE[typename]
    end
    ty = CLRBridge.GetType(typename)
    if !isnull(ty)
        TYPE_CACHE[typename] = ty
    end
    return ty
end

macro T_str(name)
    :(gettypecached($name))
end

function clrtypeof(x::CLRObject)
    CLRBridge.GetObjectType(gethandle(x))
end

function isclrtype(x::CLRObject, t::CLRObject)
    xt = clrtypeof(x)
    ret = CLRBridge.InvokeMember(gethandle(clrtypeof(xt)), "Equals", CLRBridge.BindingFlags.InvokeMethod,
    0, gethandle(xt), [gethandle(t)])
    return CLRBridge.GetBool(gethandle(ret))
end

function isassignable(totype::CLRObject, fromtype::CLRObject)
    ret = CLRBridge.InvokeMember(gethandle(T"System.Type"), "IsAssignableFrom",
    CLRBridge.BindingFlags.InvokeMethod, 0, gethandle(totype), [gethandle(fromtype)])
    return CLRBridge.GetBool(gethandle(ret))
end


text:
"""
The `DefUnPack` module only exports the `@defunpack` macro,
which defines a macro unpacking properties of an object.

__Simple Intended Usage:__
```julia
julia> using .DefUnPack

julia> struct Foo{A, B, C} a::A; b::B; c::C end

julia> @defunpack all_Foo Foo
@unpackall_Foo (macro with 1 method)

help?> @unpackall_Foo
@unpackall_Foo(obj) unpacks the fields (:a, :b, :c) of obj.

julia> @macroexpand @unpackall_Foo foo
:((a, b, c) = (foo.a, foo.b, foo.c))

julia> @unpackall_Foo Foo(1, 2.0, "three")
(1, 2.0, "three")

julia> a, b, c
(1, 2.0, "three")
```
"""
module DefUnPack

export @defunpack

"""
    @defunpack(name::Symbol, expr)

defines the macro named `Symbol(:unpack, name)`
which unpacks the properties specified by `expr` of an object.

Let `val` be the value of `expr`. Then the list of the unpacking properties is set to

* `val` if `val` is a tuple of symbols,
* `fieldnames(val)` if `val` is a type,
* `propertynames(val)` otherwise.

### Example

```julia
julia> @defunpack _cat_and_dog (:cat, :dog)
@unpack_cat_and_dog (macro with 1 method)

help?> @unpack_cat_and_dog
@unpack_cat_and_dog(obj) unpacks the fields (:cat, :dog) of obj.

julia> @macroexpand @unpack_cat_and_dog x
:((cat, dog) = (x.cat, x.dog))

julia> @unpack_cat_and_dog (dog = "bowwow", mouse = "squeak", cat = "meow")
("meow", "bowwow")

julia> cat, dog
("meow", "bowwow")

julia> struct Foo{A, B, C} a::A; b::B; c::C end

julia> @defunpack all_Foo Foo
@unpackall_Foo (macro with 1 method)

help?> @unpackall_Foo
@unpackall_Foo(obj) unpacks the fields (:a, :b, :c) of obj.

julia> @macroexpand @unpackall_Foo foo
:((a, b, c) = (foo.a, foo.b, foo.c))

julia> @unpackall_Foo Foo(1, 2.0, "three")
(1, 2.0, "three")

julia> a, b, c
(1, 2.0, "three")

julia> baz = (p = "one", q = 2.0, r = 3)

julia> @defunpack all_baz baz
@unpackall_baz (macro with 1 method)

help?> @unpackall_baz
@unpackall_baz(obj) unpacks the fields (:p, :q, :r) of obj.

julia> @macroexpand @unpackall_baz baz
:((p, q, r) = (baz.p, baz.q, baz.r))

julia> @unpackall_baz baz
("one", '2', 3)

julia> p, q, r
("one", '2', 3)
```

### Appendix

Instead of using `@defunpack`, we can also define the unpacking macro directly as follows:

```julia
struct Foo{A, B, C} a::A; b::B; c::C end

\"\"\"
`@unpackall_Foo(obj)` unpacks all fields of the object `obj` of type `Foo`.
\"\"\"
macro unpackall_Foo(obj)
    names = fieldnames(Foo)
    Expr(:(=),
        Expr(:tuple, names...),
        Expr(:tuple, (:(\$obj.\$name) for name in names)...)
    ) |> esc
end

@macroexpand @unpackall_Foo foo
```
```
:((a, b, c) = (foo.a, foo.b, foo.c))
```
"""
macro defunpack(name::Symbol, expr)
    macroname = Symbol(:unpack, name)
    atmacroname = Symbol('@', macroname)
    val = Core.eval(__module__, expr)
    names = val isa Tuple{Vararg{Symbol}} ? val :
            val isa Type ? fieldnames(val) : propertynames(val)
    docstr = """`$atmacroname(obj)` unpacks the properties `$names` of `obj`."""
    quote
        macro $macroname(obj)
            Expr(:(=),
                Expr(:tuple, $names...),
                Expr(:tuple, (:($obj.$name) for name in $names)...)
            ) |> esc
        end
        @doc $docstr $(:($atmacroname))
        $atmacroname
    end |> esc
end

unpackall_Foo = """
macro unpackall_Foo(obj)
    names = fieldnames(Foo)
    Expr(:(=),
        Expr(:tuple, names...),
        Expr(:tuple, (:(\$obj.\$name) for name in names)...)
    ) |> esc
end
"""

end


text:
using Test, LinearAlgebra
using CVChannel

@testset "./src/optimizations.jl" begin

maxEntState = 0.5 * [1 0 0 1 ; 0 0 0 0 ; 0 0 0 0 ; 1 0 0 1]
maxMixState = 0.25 * [1 0 0 0 ; 0 1 0 0 ; 0 0 1 0 ; 0 0 0 1]

max_ent_choi = Choi(2*maxEntState, 2, 2)
max_mix_choi = Choi(2*maxMixState, 2, 2)

@testset "eaCV" begin
    primal_opt_tuple = eaCV(max_ent_choi)
    @test length(primal_opt_tuple) == 2
    @test primal_opt_tuple[1] ≈ 4 atol=1e-6
    @test primal_opt_tuple[2] ≈ 2*maxEntState atol=1e-6

    dual_opt_tuple = eaCV(max_ent_choi, :dual)
    @test length(dual_opt_tuple) == 2
    @test dual_opt_tuple[1] ≈ 4 atol=1e-6
    @test dual_opt_tuple[2] ≈ [2 0;0 2]
end

@testset "eaCVPrimal" begin
    @test isapprox(eaCVPrimal(maxMixState,2,2)[1] ,1/2, atol=1e-6)
    @test isapprox(eaCVPrimal(maxEntState,2,2)[1], 2, atol = 1e-6)
    #We don't test second output on maxMix because the optimal optimizer set is too big
    @test isapprox(eaCVPrimal(maxEntState,2,2)[2], 2*maxEntState, atol = 1e-6)
end

@testset "eaCVDual" begin
    @test isapprox(eaCVDual(maxMixState,2,2)[1] ,1/2, atol=1e-6)
    @test isapprox(eaCVDual(maxEntState,2,2)[1], 2, atol = 1e-6)
    @test isapprox(eaCVDual(maxMixState,2,2)[2],1/4*[1 0 ; 0 1], atol = 1e-6)
    @test isapprox(eaCVDual(maxEntState,2,2)[2],[1 0 ; 0 1], atol = 1e-6)
end

@testset "pptCV" begin
    primal_opt_tuple = pptCV(max_ent_choi)
    @test length(primal_opt_tuple) == 2
    @test primal_opt_tuple[1] ≈ 2 atol=1e-6
    @test primal_opt_tuple[2] ≈ [2/3 0 0 1/3 ; 0 1/3 0 0 ; 0 0 1/3 0 ; 1/3 0 0 2/3] atol=1e-2

    dual_opt_tuple = pptCV(max_ent_choi, :dual)
    @test length(dual_opt_tuple) == 3
    @test dual_opt_tuple[1] ≈ 2 atol=1e-6
    @test dual_opt_tuple[2] ≈ [1 0;0 1] atol=1e-6
    @test dual_opt_tuple[3] ≈ [0 0 0 0;0 1 -1 0;0 -1 1 0;0 0 0 0] atol=1e-6
end

@testset "pptCVPrimal" begin
    @test isapprox(pptCVPrimal(maxMixState,2,2)[1] ,1/2, atol=1e-6)
    @test isapprox(pptCVPrimal(maxEntState,2,2)[1], 1, atol = 1e-6)
    #Again we don't test second output on maxMix
    #One can determine the optimizer is what I give below, but the numerics gets only close-ish, hence atol is big
    @test isapprox(pptCVPrimal(maxEntState,2,2)[2],[2/3 0 0 1/3 ; 0 1/3 0 0 ; 0 0 1/3 0 ; 1/3 0 0 2/3], atol=10^(-2.4))
end

@testset "pptCVDual" begin
    @test isapprox(pptCVDual(maxMixState,2,2)[1] ,1/2, atol=1e-6)
    @test isapprox(pptCVDual(maxEntState,2,2)[1], 1, atol = 1e-6)
    @test isapprox(pptCVDual(maxMixState,2,2)[2],1/4*[1 0 ; 0 1], atol = 1e-6)
    @test isapprox(pptCVDual(maxEntState,2,2)[2],1/2*[1 0 ; 0 1], atol = 1e-6)
end

@testset "pptCVMultiplicativity" begin
    @testset "Choi operator arguments" begin
        ident3_choi = Choi(vec(Matrix{Int}(I,3,3))*vec(Matrix{Int}(I,3,3))',3,3)
        v1 = pptCVMultiplicativity(max_ent_choi, ident3_choi)
        v2 = pptCVMultiplicativity(max_ent_choi, ident3_choi, singular_method = :dual)
        v3 = pptCVMultiplicativity(max_ent_choi, ident3_choi, parallel_method = :primal)
        @test_throws MethodError pptCVMultiplicativity(max_ent_choi, ident3_choi, singular_method = :cat)
        @test_throws TypeError pptCVMultiplicativity(max_ent_choi, ident3_choi, parallel_method = true)
        #Check single copy
        @test isapprox(v1[1],2, atol=1e-6)
        @test isapprox(v1[2],3, atol=1e-6)
        @test isapprox(v2[1],2, atol=1e-6)
        @test isapprox(v2[2],3, atol=1e-6)
        @test isapprox(v3[1],2, atol=1e-6)
        @test isapprox(v3[2],3, atol=1e-6)
        @test v1[1] != v2[1]
        @test v1[2] != v2[2]
        #Check parallel answer
        @test v1[3]==v2[3]
        @test v1[3]!=v3[3]
        @test v2[3]!=v3[3]
        @test isapprox(v1[3],v2[3], atol=1e-5)
        @test isapprox(v1[3],v3[3], atol=1e-5)
        @test isapprox(v1[3],6, atol=1e-5)
    end
    @testset "matrix and dims arguments" begin
        choi_ident3 = vec([1 0 0 ; 0 1 0 ; 0 0 1])*vec([1 0 0 ; 0 1 0 ; 0 0 1])'
        v1 = pptCVMultiplicativity(2*maxEntState,2,2,choi_ident3,3,3)
        v2 = pptCVMultiplicativity(2*maxEntState,2,2,choi_ident3,3,3,step1isdual = true)
        v3 = pptCVMultiplicativity(2*maxEntState,2,2,choi_ident3,3,3,step2isprimal = true)
        @test_throws TypeError pptCVMultiplicativity(2*maxEntState,2,2,choi_ident3,3,3,step2isprimal = "cat")
        @test_throws TypeError pptCVMultiplicativity(2*maxEntState,2,2,choi_ident3,3,3,step1isdual = 12)
        #Check single copy
        @test isapprox(v1[1],2, atol=1e-6)
        @test isapprox(v1[2],3, atol=1e-6)
        @test isapprox(v2[1],2, atol=1e-6)
        @test isapprox(v2[2],3, atol=1e-6)
        @test isapprox(v3[1],2, atol=1e-6)
        @test isapprox(v3[2],3, atol=1e-6)
        @test v1[1] != v2[1]
        @test v1[2] != v2[2]
        #Check parallel answer
        @test v1[3]==v2[3]
        @test v1[3]!=v3[3]
        @test v2[3]!=v3[3]
        @test isapprox(v1[3],v2[3], atol=1e-5)
        @test isapprox(v1[3],v3[3], atol=1e-5)
        @test isapprox(v1[3],6, atol=1e-5)
    end
end

@testset "WHIDLP" begin
    identChan(X) = X
    conditions = [[2,2,0],[3,2,0],[2,3,0.25]]
    for (d1f,d2f,λ) in conditions
        d1 = Int(d1f)
        d2 = Int(d2f)

        wern_choi = d1*wernerState(d1, λ)
        ident_choi = choi(identChan, d2, d2)
        kron_par_choi = kron(wern_choi,ident_choi)
        par_choi = permuteSubsystems(kron_par_choi, [1,3,2,4], [d1,d1,d2,d2])
        par_cv = pptCVPrimal(par_choi, d1*d2, d1*d2)
        par_cv_LP = WHIDLP(d1, d2, λ)
        @test isapprox(par_cv[1], par_cv_LP, atol=1e-5)
    end
end

@testset "generalWHLPConstraints" begin
    #Domain Error Checks
    @test_throws DomainError generalWHLPConstraints(3,3,ones(4))
    @test_throws DomainError generalWHLPConstraints(2,4,[1.,6])

    #Generating Constraints Checks
    n,d,λ = 1,2,1
    A,B,g,a = generalWHLPConstraints(n,d,λ*ones(n))
    @testset "n,d,λ = 1,2,1" begin
        @test A == [1 1 ; 1 -1]
        @test B == [1 0 ; 1 d]
        @test a == reshape([1.;(2λ-1)],:,1)
        @test g == [d 1]
    end

    λ = 0.8
    A,B,g,a = generalWHLPConstraints(n,d,λ*ones(n))
    @testset "n,d,λ = 1,2,0.8" begin
        @test A == [1 1 ; 1 -1]
        @test B == [1 0 ; 1 d]
        @test a == reshape([1.;(2λ-1)],:,1)
        @test g == [d 1]
    end

    d = 4
    A,B,g,a = generalWHLPConstraints(n,d,λ*ones(n))
    @testset "n,d,λ = 1,4,0.8" begin
        @test A == [1 1 ; 1 -1]
        @test B == [1 0 ; 1 d]
        @test a == reshape([1.;(2λ-1)],:,1)
        @test g == [d 1]
    end

    n,d,λ = 2,3,1
    A,B,g,a = generalWHLPConstraints(n,d,λ*ones(n))
    targA = [1 1 1 1 ; 1 -1 1 -1 ; 1 1 -1 -1 ; 1 -1 -1 1]
    targB = [1. 0 0 0 ; 1 d 0 0 ; 1 0 d 0 ; 1 d d d^2]
    targa = reshape([1. 1 1 1],:,1)
    targg = [d^2 d d 1]
    @testset "n,d,λ = 2,3,1" begin
        @test A == targA
        @test B == targB
        @test a == targa
        @test g == targg
    end
end

@testset "wernerHolevoCVPPT" begin
    n,d,λ = 1,3,0
    A,B,g,a = generalWHLPConstraints(n,d,λ*ones(n))
    lp_cv_ppt, opt = wernerHolevoCVPPT(n,d,A,B,g,a)
    orig_choi = 3*wernerState(3,0)
    test1 = pptCVDual(orig_choi,3,3)
    @test isapprox(test1[1], lp_cv_ppt, atol = 1e-6)

    λ = 0.2
    A,B,g,a = generalWHLPConstraints(n,d,λ*ones(n))
    lp_cv_ppt, opt = wernerHolevoCVPPT(n,d,A,B,g,a)
    choi_2 = 3*wernerState(3,λ)
    test2 = pptCVDual(choi_2,3,3)
    @test isapprox(test2[1], lp_cv_ppt, atol = 1e-6)

    n,λ = 2,0
    A,B,g,a = generalWHLPConstraints(n,d,λ*ones(n))
    lp_cv_ppt, opt = wernerHolevoCVPPT(n,d,A,B,g,a)
    kron_par_choi = kron(orig_choi,orig_choi)
    par_choi = permuteSubsystems(kron_par_choi,[1,3,2,4],[3,3,3,3])
    test3 = pptCVPrimal(par_choi,9,9)
    @test isapprox(test3[1], lp_cv_ppt, atol = 1e-5)
end

@testset "twoSymCVPrimal" begin
    #We can't really test a lot with this since we don't really know much
    #Here we verify that it gives answers we know it should give even over
    #the separable cone
    test1 = 3*wernerState(3,0)
    choi1 = Choi(test1, 3, 3)
    cv1 = pptCVPrimal(test1,3,3)
    cv1_two_sym = twoSymCVPrimal(test1,3,3)
    choi_cv1_two_sym = twoSymCVPrimal(choi1)
    @test isapprox(cv1[1], cv1_two_sym[1], atol=1e-6)
    @test isapprox(cv1[1], choi_cv1_two_sym[1], atol=1e-6)

    test2 = vec([1 0 0 ; 0 1 0 ; 0 0 1])*vec([1 0 0 ; 0 1 0 ; 0 0 1])'
    choi2 = Choi(test2, 3, 3)
    cv2 = pptCVPrimal(test2,3,3)
    cv2_two_sym = twoSymCVPrimal(test2,3,3)
    choi_cv2_two_sym = twoSymCVPrimal(choi2)
    @test isapprox(cv2[1], cv2_two_sym[1], atol=2e-6)
    @test isapprox(cv2[1], choi_cv2_two_sym[1], atol=2e-6)
end
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
using LinearAlgebra, DataFrames, CSV, Plots, Optim, Statistics

# ============================================================
# Ler dados
# ============================================================
data = CSV.read("MER_T07_01.csv", DataFrame);
YYYYMM = data[!,"YYYYMM"]
index = mod.(YYYYMM,100) .!= 13
Value = data[!,"Value"]
y = parse.(Float64,Value[index][1:588])
y_insample = y[1:120]
y_outofsample =y[121:132]
n = length(y_insample)

# ============================================================
# MEB com tendência linear local e sazonalidade estocástica
# ============================================================

function MEE(s::Int)
    if s == 6
        m = 2*s + 1
    else
        m = 2*s + 2
    end
    Z = [1 0 1 0 1 0 1 0 1 0 1 0 1];
    Z = Z[1:m]
    T = zeros(m, m);
    T[1,1] = 1;
    T[1,2] = 1;
    T[2,2] = 1;
    for j in 1:s
        if j == 6
            break
        else
            c = [[cos(pi*j/6), -sin(pi*j/6)] [sin(pi*j/6), cos(pi*j/6)]];
            T[2*j+1:2*j+2, 2*j+1:2*j+2] = c
        end
    end
    if s == 6
        T[m, m] = -1;
    end
    R = I;
    return Z, T, R, m
end
s = 3
Z, T, R, m = MEE(s)
# ============================================================
# Inicialização
# ============================================================
# Big Kappa
a = Vector{Vector{Real}}(undef, n+1);
a[1] = zeros(m);
P = Vector{Array{Real, 2}}(undef, n+1);
P[1] = diagm([1e6 for _ in 1:m]);

# Resultados armazenados para o smoother
F_aux = Vector{Any}(undef, n);
K_aux = Vector{Any}(undef, n);
V_aux = Vector{Any}(undef, n);

# ============================================================
# Estimação do Filtro de Kalman e cálculo da log verossimilhança
# ============================================================

function get_minus_log_likelihood(des)

    des_epsilon, des_eta, des_csi, des_omega = des
    var_epsilon = des_epsilon^2
    var_eta = des_eta^2
    var_csi = des_csi^2
    var_omega = des_omega^2
    # H
    H = var_epsilon;
    # Q
    Q = zeros(m,m);
    Q[1,1] = var_eta;
    Q[2,2] = var_csi;
    for j in 3:m
        Q[j,j]=var_omega;
    end
    # FK
    likelihood = 0.0;
    for t in 1:n

        F = Z'*P[t]*Z + H;
        K = T*P[t]*Z * ((F)^(-1));
        V = y_insample[t] - Z'*a[t];

        a[t+1] = T*a[t] + K*V;
        P[t+1] = T*P[t] * (T - K*Z')' + R*Q*R';

        if t > m
            likelihood += -log(2*pi)/2 + (-1/2)*log(F[1]) + (-1/2)*(V[1]^2)/F[1];
        end
        # Resultados armazenados para o smoother
        F_aux[t] = F;
        K_aux[t] = K;
        V_aux[t] = V;
    end
    return -likelihood/(n-m)
end

res = optimize(get_minus_log_likelihood, [sqrt(var(y_insample)), sqrt(var(y_insample)), 1.0, 1.0], LBFGS(),Optim.Options(show_trace=true)) #NelderMead

function AIC(ll,q,n,w)
    return (1/(n-q))*(-2*(n-q)*ll+2*(q+w))
end
ll = -res.minimum*(n-m)
@show ll
@show AIC(ll ,m,n,4)
# ============================================================
# Suavizador
# ============================================================

# Suavizador
r = Vector{Vector{Real}}(undef, n+1);
r[end] = zeros(m);
N = Vector{Array{Real, 2}}(undef, n+1);
N[end] = zeros(m,m);
alpha_hat = Vector{Vector{Real}}(undef, n);
V_true = Vector{Matrix{Real}}(undef, n);

for t in n:-1:1
    L = T - K_aux[t]*Z';
    r[t] = Z * ((F_aux[t])^(-1)) * V_aux[t] + L'*r[t+1];
    N[t] = Z * ((F_aux[t])^(-1)) * Z' + L' * N[t+1] * L;
end
for t in 1:n
    alpha_hat[t] = a[t] + P[t+1]*r[t];
    V_true[t] = P[t+1] - P[t+1]*N[t]*P[t+1];
end

# ============================================================
# Plot
# ============================================================

# Y
plot(y_insample, label = "Real", legend=:bottomright)
plot!([dot(Z,x) for x in a], label = "FK")
plot!([dot(Z,x) for x in alpha_hat], label = "Smoother")

savefig("figs\\filter_$(s).png")

psi_val = res.minimizer.^2;
println("Level: $(psi_val[2])")
println("Slope: $(psi_val[3])")
println("Seasonal: $(psi_val[4])")
println("Irregular: $(psi_val[1])")

# Seasonality
plot([dot(Z[3:end],x[3:end]) for x in alpha_hat], label = "Seasonality")
savefig("figs\\seasonality_$(s).png")
# Level
plot([x for x in y_insample], label = "Real")
plot!([dot(Z[1],x[1]) for x in alpha_hat], label = "Level")
savefig("figs\\level_$(s).png")

# Slope
plot([x[2] for x in alpha_hat], label = "Slope")
savefig("figs\\slope_$(s).png")


text:
# Gradient Boosted Decision Trees
module GBDecisionTree

using DecisionTree
using DataStructures

importall GradientBoost.GB
importall GradientBoost.LossFunctions

export GBDT,
       build_base_func

# Gradient boosted decision tree algorithm.
type GBDT <: GBAlgorithm
  loss_function::LossFunction
  sampling_rate::FloatingPoint
  learning_rate::FloatingPoint
  num_iterations::Int
  tree_options::Dict

  function GBDT(;loss_function=LeastSquares(),
    sampling_rate=0.6, learning_rate=0.1,
    num_iterations=100, tree_options=Dict())

    default_options = {
      :maxlabels => 5,
      :nsubfeatures => 0
    }
    options = merge(default_options, tree_options)
    new(loss_function, sampling_rate, learning_rate, num_iterations, options)
  end
end

function GB.build_base_func(
  gb::GBDT,
  instances,
  labels,
  prev_func_pred,
  psuedo)

  # Train learner
  model = build_tree(
    psuedo, instances,
    gb.tree_options[:maxlabels],
    gb.tree_options[:nsubfeatures]
  )
  psuedo_pred = apply_tree(model, instances)

  # Update regions (leaves)
  # NOTE(svs14): Trees are immutable,
  #              override leaves by having node-to-val mapping.
  inst_node_index = InstanceNodeIndex(model, instances)
  function val_func(node)
    inst_ind = inst_node_index.n2i[node]

    # If loss function is LeastSquares, we don't need need to change values.
    if typeof(gb.loss_function) <: LeastSquares
      val = node.majority
    else
      val = fit_best_constant(gb.loss_function,
        labels[inst_ind],
        psuedo[inst_ind],
        psuedo_pred[inst_ind],
        prev_func_pred[inst_ind]
      )
    end

    val
  end
  val_type = eltype(prev_func_pred)
  n2v = Dict{Leaf, val_type}()
  update_regions!(n2v, model, val_func)

  # Prediction function
  function pred(instances)
    num_instances = size(instances, 1)
    predictions = [
      n2v[instance_to_node(model, instances[i,:])]
      for i in 1:num_instances
    ]
    predictions
  end

  # Produce function that delegates prediction to model
  return (instances) -> pred(instances)
end

# DT Helper Functions

type InstanceNodeIndex
  i2n::Vector{Leaf}
  n2i::DefaultDict{Leaf, Vector{Int}}

  function InstanceNodeIndex(tree::Union(Leaf,Node), instances)
    num_instances = size(instances, 1)
    i2n = Array(Leaf, num_instances)
    n2i = DefaultDict(Leaf, Vector{Int}, () -> Int[])

    for i = 1:num_instances
      node = instance_to_node(tree, instances[i,:])
      i2n[i] = node
      push!(n2i[node], i)
    end

    new(i2n, n2i)
  end
end

# Returns respective node of instance.
function instance_to_node(tree::Node, instance)
  # Code adapted from DecisionTree.jl
  features = instance
  if tree.featval == nothing || features[tree.featid] < tree.featval
    return instance_to_node(tree.left, features)
  else
    return instance_to_node(tree.right, features)
  end
end
function instance_to_node(leaf::Leaf, instance)
  return leaf
end

# Update region by having updated leaf value encoded
# in a leaf-to-value mapping.
function update_regions!{T}(n2v::Dict{Leaf, T}, node::Node, val_func::Function)
  update_regions!(n2v, node.left, val_func)
  update_regions!(n2v, node.right, val_func)
end
function update_regions!{T}(n2v::Dict{Leaf, T}, leaf::Leaf, val_func::Function)
  n2v[leaf] = val_func(leaf)
end

# Loss function fits
function fit_best_constant(lf::LeastAbsoluteDeviation,
  labels, psuedo, psuedo_pred, prev_func_pred)

  values = labels .- prev_func_pred
  median(values)
end
function fit_best_constant(lf::BinomialDeviance,
  labels, psuedo, psuedo_pred, prev_func_pred)

  num = sum(psuedo)
  denom = sum((labels .- psuedo) .* (1 .- labels .+ psuedo))
  if denom == 0.0
    return 0.0
  else
    return num / denom
  end
end

end # module


text:
__precompile__(true)

module FastFilter

using Random
include("utils.jl")
include("binaryfusefilter.jl")
export BinaryFuseFilter

end


text:
#######################################################################
# Graceful exit
# and macro for printing information depending on verbosity
#
# ASSUMPTION: the package defines variable `VERBOSE`
#######################################################################

# String → Nothing
# Prints error message and terminates execution
exitErrWithMsg(msg :: String, err :: Exception) = begin
    @error(msg, err)
    exit(1)
end
exitErrWithMsg(msg :: String) = begin
    @error(msg)
    exit(1)
end

# If verbose, prints information `info`
macro status(info)
    :(if VERBOSE ; @info($(esc(info))) end)
end

splitlines(str) = split(str, r"\n|\r\n")


text:
module SparseMerkleTree

# Based on a slightly modified version of vbuterins new_bin_trie.py

using MLStyle
using SHA
using BitIntegers

import Base.getindex
import Base.setindex!

export new_tree, merkle_proof, verify_proof, compress_proof, decompress_proof
export Hash

# Assumes system is little endian
function from_bytes(b)::UInt256
    reinterpret(UInt256, reverse(b))[1]
end

# In system byte order
function to_bytes(b)
    reinterpret(UInt8, [b])
end

function to_bytes(b::Array{UInt8})
    b
end

struct Hash
    value::UInt256
end

function Base.:(==)(x::Hash, y::Hash)
    x.value == y.value
end

function Base.:+(x::Hash, y::Hash)
    bytes = map(to_bytes, (x.value, y.value))
    bytes = map(reverse, bytes)
    Hash(from_bytes(sha256(vcat(bytes[1], bytes[2]))))
end

zerohashes = [Hash(from_bytes(sha256(repeat(b"\0", 32))))]
for i in 1:256
    insert!(zerohashes, 1, zerohashes[1] + zerohashes[1])
end

NodeValue = Array{Array{UInt8, 1}, 1}  # Array of 2 elements with child hashes
LeafValue = Array{UInt8}

mutable struct SMT
    db::Dict{UInt256, Union{NodeValue, LeafValue}}
    root::Hash
end

function new_tree()
    smt = SMT(Dict(), Hash(0))
    db = smt.db
    h = sha256(repeat(b"\0", 32))
    for i in 1:256
        hh = vcat(h, h)
        newh = sha256(hh)
        db[from_bytes(newh)] = [h, h]
        h = newh
    end
    smt.root = Hash(from_bytes(h))
    return smt
end

function getindex(smt::SMT, k::Hash)
    x = smt.db[k.value]
    @match x begin
    nothing => nothing
    [a, b] => (from_bytes(a), from_bytes(b))
    end
end

function bitmask(n::UInt256)
    x = BitVector(undef, 256)
    n_chunks = reinterpret(UInt64, [n])
    for (i, ch) in enumerate(n_chunks)
        x.chunks[i] = ch
    end
    x
end

function getindex(smt::SMT, k::UInt256)
    s = bitmask(k)
    v = smt.root
    for b in s
        v = @match b begin
        false => Hash(smt[v][1])
        true  => Hash(smt[v][2])
        end
    end
    return v
end

struct Proof
    sidenodes
    zero_bitmask
end

function Base.:(==)(x::Proof, y::Proof)
    x.sidenodes == y.sidenodes && x.zero_bitmask == y.zero_bitmask
end

function merkle_proof(smt::SMT, k::UInt256)
    # Iterate the tree top down and compute sidenodes
    s = reverse(bitmask(k))
    cur = smt.root
    sidenodes = []
    for b in s
        children = smt[cur]
        cur, side = @match b begin
        true  => Hash.(reverse(children))
        false => Hash.(children)
        end
        append!(sidenodes, side.value)
    end
    return Proof(Hash.(sidenodes), 0)
end

function setindex!(smt::SMT, v::V, k::UInt256) where V
    sidenodes = merkle_proof(smt, k).sidenodes
    # Iterate the tree bottom up updating hashes
    cur = to_bytes(v)
    curhash = Hash(from_bytes(sha256(cur)))
    smt.db[curhash.value] = v
    for b in bitmask(k)
        side = last(sidenodes)
        newh, left, right = @match b begin
        true => (side + curhash, side, curhash)
        false => (curhash + side, curhash, side)
        end
        smt.db[newh.value] = [reverse(to_bytes(left.value)), reverse(to_bytes(right.value))]
        pop!(sidenodes)
        curhash = newh
    end
    smt.root = curhash
end

function verify_proof(proof::Proof, root, k, v)
    sidenodes = reverse(proof.sidenodes)
    # Iterate the tree bottom up updating hashes
    cur = to_bytes(v)
    curhash = Hash(from_bytes(sha256(cur)))
    for (i, b) in enumerate(bitmask(k))
        side = sidenodes[i]
        newh, left, right = @match b begin
        true => (side + curhash, side, curhash)
        false => (curhash + side, curhash, side)
        end
        curhash = newh
    end
    return root.value == curhash.value
end

function compress_proof(proof)
    sidenodes = proof.sidenodes
    mask = (sidenodes .== zerohashes[2:end])
    compressed = [s.value for (s,b) in zip(sidenodes, mask) if !b]
    Proof(Hash.(compressed), mask)
end

function decompress_proof(proof)
    sidenodes = proof.sidenodes
    mask = proof.zero_bitmask
    decompressed = []
    for (z, b) in zip(zerohashes[2:end], mask)
        d = if b
	    z.value
        else
            s = pop!(sidenodes)
            s.value
        end
        append!(decompressed, d)
    end
    Proof(Hash.(decompressed), 0)
end

end


text:
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Parse CoNLL-X (also known as CoNLL 2006) style dependency data.
#
# Data format description:
#
#   http://nextens.uvt.nl/depparse-wiki/DataFormat
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-05-01

module CoNLLX

export NOHEAD, NOVAL, Sentence, Sentences, Token, conllxparse

type Token
    # [2,|Sentence|+1] sentence-internal unique token identifier.
    id::Uint
    # Token form, the actual text.
    form::String
    # Lemma of the token form.
    lemma::String
    # Coarse part-of-speech.
    cpostag::String
    # Fine-grained part-of-speech.
    postag::String
    # Set of syntactic and/or morphological features.
    feats::String # TODO: Should be an ordered set, from DataStructures.jl?
    # [1,|Sentence|+1] id of the head of the token, 0 for no head.
    head::Uint
    # Dependency relation to the head.
    deprel::String
    # [1,|Sentence|+1] id of the projective head of the token, 0 for no phead.
    phead::Uint
    # Dependency relation to the projective head.
    pdeprel::String
end

const NOHEAD = 0
const NOVAL = "_"

import Base: show
function show(io::IO, t::Token)
    head, phead = map(i -> i == NOHEAD ? NOVAL : string(i - 1),
        (t.head, t.phead))

    print(io, string("$(t.id-1)\t$(t.form)\t$(t.lemma)\t$(t.cpostag)\t",
        "$(t.postag)\t$(t.feats)\t$head\t$(t.deprel)\t$phead\t",
        "$(t.pdeprel)"))
end

function useproj!(t::Token)
    t.head = t.phead
    t.deprel = t.pdeprel
end

function blind!(t::Token)
    t.head, t.phead = (NOHEAD, NOHEAD)
    t.deprel, t.pdeprel = (NOVAL, NOVAL)
end

typealias Sentence Vector{Token}

import Base: show
function show(io::IO, s::Sentence)
    print(io, join([string(t) for t in s], '\n'))
end

typealias Sentences Vector{Sentence}

import Base: show
function show(io::IO, ss::Sentences)
    print(io, join([string(s) for s in ss], "\n\n"))
    print(io, '\n')
end

# Not the friendliest of expressions, but it will do the trick.
const CONLLXREGEX = Regex(string("^([0-9]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t",
    "([^\t]+)\t([^\t]+)\t([0-9]+)\t([^\t]+)\t([0-9]+|_)\t([^\t]+)\n\$"))

# Parsing implemented using the iterator protocol.
type CoNLLXParse
    stream::IO
    useproj::Bool
    blind::Bool
    linenum::Uint
end

function conllxparse(stream::IO; blind=false, useproj=false)
    return CoNLLXParse(stream, useproj, blind, 0)
end

import Base: start
start(::CoNLLXParse) = nothing

import Base: next
function next(itr::CoNLLXParse, nada)
    sent = Token[]

    for line in eachline(itr.stream)
        itr.linenum += 1

        # Ignore comments (this is a format extension).
        if line[1] == '#'
            continue
        end

        # Empty lines are used as sentence separators.
        if line[1] == '\n'
            if !isempty(sent)
                return (sent, nothing)
            end
            continue
        end

        soup = match(CONLLXREGEX, line)
        if soup == nothing
            warn(string("failed to parse:$(itr.linenum) \"",
                rstrip(line, '\n'), "\""))
            continue
        end

        (id_str, form, lemma, cpostag, postag, feats, head_str, deprel,
            phead_str, pdeprel) = soup.captures
        id, head = map(int, (id_str, head_str))

        # Bump the ids internally to allow it to double as an index (necessary
        #   due to Julia being 1-indexed).
        id += 1
        head += 1

        # The projective head is optional.
        if phead_str != NOVAL
            # Convert and bump the id.
            phead = int(phead_str) + 1
        else
            phead = NOHEAD
        end

        tok = Token(id, form, lemma, cpostag, postag, feats, head, deprel,
            phead, pdeprel)

        # Override the potentially non-projective head/deprel with their
        #   projective counterparts.
        if itr.useproj
            useproj!(tok)
        end

        # Remove (blind) any gold annotations.
        if itr.blind
            blind!(tok)
        end

        push!(sent, tok)
    end

    # Yield what we have observed before EOF (could be an empty Sentence).
    return (sent, nothing)
end

import Base: done
done(itr::CoNLLXParse, nada) = eof(itr.stream)

end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
facts("ArrayQFunction construction") do
    @fact typeof(ArrayQFunction(zeros(5, 3))) => ArrayQFunction{Float64}
    @fact ArrayQFunction(5, 3) == ArrayQFunction(zeros(5, 3)) => true
    @fact ArrayQFunction(Float32, 5, 3) == ArrayQFunction(zeros(Float32, 5, 3)) => true
    @fact_throws MethodError ArrayQFunction(zeros(5))
    @fact_throws MethodError ArrayQFunction(zeros(5, 3, 2))
end

facts("ArrayQFunction methods") do
    qf = ArrayQFunction(Float64[1 2 3; 3 2 1; 2 2 2; 1 1 3; 2 3 1])
    @fact value(qf) => [3.0, 3.0, 2.0, 3.0, 3.0]
    @fact policy(qf) => [3, 1, 1, 3, 2]
    @fact policy(UInt8, qf) => [0x03, 0x01, 0x01, 0x03, 0x02]
    @fact value(qf, 2) => 3.0
    @fact policy(qf, 2) => 1
    @fact getvalue(qf, 1, 2) => 2.0
    @fact num_states(qf) => 5
    @fact num_actions(qf) => 3
    newq  = copy(qf)
    @fact newq == qf => true
    @fact newq === qf => false
    @fact isequal(newq, qf) => true
    @fact valuetype(qf) => Float64
end

facts("ArrayQFunction setvalue!") do
    qf = ArrayQFunction(Float64[1 2 3; 3 2 1; 2 2 2; 1 1 3; 2 3 1])
    setvalue!(qf, 4.0, 1, 2)
    @fact getvalue(qf, 1, 2) => 4.0
    setvalue!(qf, 3, 2, 1)
    @fact getvalue(qf, 2, 1) => 3.0
end

facts("VectorQFunction construction") do
    @fact typeof(VectorQFunction(zeros(5), zeros(Int, 5))) => VectorQFunction{Float64,Int}
    @fact VectorQFunction(5) == VectorQFunction(zeros(5), zeros(Int, 5)) => true
    @fact VectorQFunction(Float32, UInt8, 5) == VectorQFunction(zeros(Float32, 5), zeros(UInt8, 5)) => true
    @fact_throws DimensionMismatch VectorQFunction(zeros(5), zeros(Int, 6))
end

facts("VectorQFunction methods") do
    qf = VectorQFunction([3.0, 3.0, 2.0, 3.0, 3.0], [3, 1, 1, 3, 2])
    @fact value(qf) => [3.0, 3.0, 2.0, 3.0, 3.0]
    @fact policy(qf) => [3, 1, 1, 3, 2]
    @fact value(qf, 2) => 3.0
    @fact policy(qf, 2) => 1
    @fact_throws MethodError value(qf, 1, 2)
    @fact num_states(qf) => 5
    @fact_throws MethodError num_actions(qf)
    qf_new = copy(qf)
    @fact qf_new == qf => true
    @fact qf_new === qf => false
    @fact isequal(qf_new, qf) => true
    @fact valuetype(qf) => Float64
end

facts("VectorQFunction setvalue! assigned") do
    qf = VectorQFunction([3.0, 3.0, 2.0, 3.0, 3.0], [3, 1, 1, 3, 2])
    setvalue!(qf, 4.0, 1, 2)
    @fact value(qf, 1) => 4.0
    @fact policy(qf, 1) => 2
end

facts("VectorQFunction setvalue! not assigned") do
    qf = VectorQFunction([3.0, 3.0, 2.0, 3.0, 3.0], [3, 1, 1, 3, 2])
    setvalue!(qf, 2.0, 1, 2)
    @fact value(qf, 1) => 3.0
    @fact policy(qf, 1) => 3
end

facts("helper constructor") do
    @fact typeof(QFunction(zeros(5, 3))) => ArrayQFunction{Float64}
    @fact typeof(QFunction(zeros(5), zeros(Int, 5))) => VectorQFunction{Float64,Int}
end

facts("copy value! to vector") do
    qf = ArrayQFunction(rand(5, 3))
    val = Array(Float64, 5)
    value!(val, qf)
    @fact val == value(qf) => true
    @fact val === value(qf) => false
    @fact_throws BoundsError value!(Array(Float64, 4), qf)
    qf = VectorQFunction(rand(5), rand(Int, 5))
    val = Array(Float64, 5)
    value!(val, qf)
    @fact val == value(qf) => true
    @fact val === value(qf) => false
    @fact_throws BoundsError value!(Array(Float64, 4), qf)
end

facts("copy policy! to vector") do
    qf = ArrayQFunction(rand(5, 3))
    pol = Array(Int, 5)
    policy!(pol, qf)
    @fact pol == policy(qf) => true
    @fact pol === policy(qf) => false
    @fact_throws BoundsError policy!(Array(Int, 4), qf)
    qf = VectorQFunction(rand(5), rand(Int, 5))
    pol = Array(Int, 5)
    policy!(pol, qf)
    @fact pol == policy(qf) => true
    @fact pol === policy(qf) => false
    @fact_throws BoundsError policy!(Array(Int, 4), qf)
end


text:
function pose_to_matrix(pose::Pose)::Matrix{Float32}
    R = Matrix{Float32}(pose.orientation)
    mat = zeros(Float32,4,4)
    mat[end,end] = 1
    mat[1:3,1:3] = R
    mat[1:3,4] = pose.pos
    return mat
end

function I4(t)
    x = zeros(t, 4,4)
    x[1,1] = 1.0
    x[2,2] = 1.0
    x[3,3] = 1.0
    x[4,4] = 1.0
    x
end

function compute_projection_matrix(fx, fy, cx, cy, near, far, skew=0f0)
    proj = I4(Float32)
    proj[1, 1] = fx
    proj[2, 2] = fy
    proj[1, 2] = skew
    proj[1, 3] = -cx
    proj[2, 3] = -cy
    proj[3, 3] = near + far
    proj[3, 4] = near * far
    proj[4, 4] = 0.0f0
    proj[4, 3] = -1f0
    return proj
end

function compute_ortho_matrix(left, right, bottom, top, near, far)
    ortho = I4(Float32)
    ortho[1, 1] = 2f0 / (right-left)
    ortho[2, 2] = 2f0 / (top-bottom)
    ortho[3, 3] = - 2f0 / (far - near)
    ortho[1, 4] = - (right + left) / (right - left)
    ortho[2, 4] = - (top + bottom) / (top - bottom)
    ortho[3, 4] = - (far + near) / (far - near)
    return ortho
end

function get_perspective_matrix(width, height, fx, fy, cx, cy, near, far)
    # (height-cy) is used instead of cy because of the difference between
    # image coordinate systems between OpenCV and OpenGL. In the former,
    # the origin is at the top-left of the image while in the latter the
    # origin is at the bottom-left.
    proj_matrix = compute_projection_matrix(
            fx, fy, cx, (height-cy),
            near, far, 0.f0)
    ndc_matrix = compute_ortho_matrix(0, width, 0, height, near, far)
    ndc_matrix * proj_matrix
end

function CameraIntrinsics_from_fov_aspect(width, height, fov_y, aspect_ratio, near, far)
    # Camera principal point is the center of the image.
    cx, cy = width / 2.0, height / 2.0

    # Vertical field of view is given.
    fov_y = deg2rad(fov_y)
    # Convert field of view to distance to scale by aspect ratio and
    # convert back to radians to recover the horizontal field of view.
    fov_x = 2 * atan(aspect_ratio * tan(fov_y / 2.0))

    # Use the following relation to recover the focal length:
    #   FOV = 2 * atan( (0.5 * IMAGE_PLANE_SIZE) / FOCAL_LENGTH )
    fx = cx / tan(fov_x / 2.0)
    fy = cy / tan(fov_y / 2.0)

    CameraIntrinsics(width, height,
        fx, fy, cx, cy,
        near, far)
end

function scale_down_camera(camera, factor)
    camera_modified = CameraIntrinsics(
        width=round(Int,camera.width/factor), height=round(Int,camera.height/factor),
        fx=camera.fx/factor, fy=camera.fy/factor, cx=camera.cx/factor, cy=camera.cy/factor,
        near=camera.near, far=camera.far)
    camera_modified
end


text:
#
# Date created: 2022-03-08
# Author: aradclif
#
#
############################################################################################
# Re-work of vfindmax.jl

function findmaxbody(N, D)
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    i = ntuple(d -> Symbol(:i_, d), N)
    b = Expr(:ref, :B, ntuple(d -> D[d] == 1 ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> D[d] == 1 ? 1 : Symbol(:i_, d), N)...)
    n = Expr(:(=), :n, Expr(:call, :<, b, a))
    m = Expr(:(=), b, Expr(:if, :n, a, b))
    cₐ = Expr(:(=), c, Expr(:if, :n, Expr(:tuple, i...), c))
    push!(body.args, n)
    push!(body.args, m)
    push!(body.args, cₐ)
    body
end
function lvfindmaxbody(N, D)
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    i = ntuple(d -> Symbol(:i_, d), N)
    b = Expr(:ref, :B, ntuple(d -> D[d] == 1 ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> D[d] == 1 ? 1 : Symbol(:i_, d), N)...)
    j = Expr(:(=), :j, c)
    v = Expr(:(=), :v, b)
    n = Expr(:(=), :n, Expr(:call, :<, :v, a))
    # m = Expr(:(=), b, Expr(:if, :n, a, b))
    m = Expr(:(=), b, Expr(:call, :ifelse, :n, a, :v))
    sym = i[N]
    l = :($sym - 1)
    for d = (N - 1):-1:1
        sym = i[d]
        l = :($sym - 1 + dims[$d] * $l)
    end
    l = :($l + 1)
    # cₐ = Expr(:(=), c, Expr(:if, :n, l, c))
    cₐ = Expr(:(=), c, Expr(:call, :ifelse, :n, l, :j))
    push!(body.args, j)
    push!(body.args, v)
    push!(body.args, n)
    push!(body.args, m)
    push!(body.args, cₐ)
    body
    # #### A different approach
    # body = Expr(:block)
end
Meta.@dump(rand(Bool) ? 5 : 10)
Meta.@dump((1,2,3))
findmaxbody(N, Dᴮ)
lvfindmaxbody(N, Dᴮ)

# ex = :(j_2 - 1)
# for i = (N - 1):-1:1
#     sym = Symbol(:j_, i)
#     ex = :($sym - 1 + D[$i] * $ex)
# end
# ex = :($ex + 1)

function findmax_quote(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = fill(typemin(eltype(A)), Dᴮ)
    C = fill(ntuple(d -> 0, N), Dᴮ)
    ls = loopgen(A)
    body = findmaxbody(N, Dᴮ)
    push!(ls.args, body)
    ls
end

A = [1 2 3; 4 5 6]
N = ndims(A)
dims = (1,)
Dᴬ = size(A)
Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
B = fill(typemin(eltype(A)), Dᴮ)
C = fill(ntuple(d -> 0, N), Dᴮ)

ls = findmax_quote(A, dims)
eval(ls)
findmax(A, dims=dims)

function findmax_quote(N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    i = ntuple(d -> Symbol(:i_, d), N)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    n = Expr(:(=), :n, Expr(:call, :<, b, a))
    m = Expr(:(=), b, Expr(:if, :n, a, b))
    cₐ = Expr(:(=), c, Expr(:if, :n, Expr(:tuple, i...), c))
    push!(body.args, n)
    push!(body.args, m)
    push!(body.args, cₐ)
    push!(ls.args, body)
    ls
end

function lvfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = fill(typemin(T), Dᴮ)
    C = fill(ntuple(d -> 0, N), Dᴮ)
    Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
    _lvfindmax!(B, C, A, Dᴮ′)
    B, CartesianIndex.(C)
end

@generated function _lvfindmax!(B::AbstractArray{T, N}, C::AbstractArray{S, N},
                                A::AbstractArray{T, N}, dims::D) where {T, N, S, D}
    findmax_quote(N, D)
end

@timev lvfindmax(A, dims)

function findmax_quote2(N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    i = ntuple(d -> Symbol(:i_, d), N)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    # n = Expr(:(=), :n, Expr(:call, :<, b, a))
    # m = Expr(:(=), b, Expr(:call, :ifelse, :n, a, b))
    # sym = i[N]
    # l = :($sym - 1)
    # for d = (N - 1):-1:1
    #     sym = i[d]
    #     l = :($sym - 1 + dims[$d] * $l)
    # end
    # l = :($l + 1)
    # cₐ = Expr(:(=), c, Expr(:call, :ifelse, :n, l, c))
    # push!(body.args, n)
    # push!(body.args, m)
    # push!(body.args, cₐ)
    # push!(ls.args, body)
    j = Expr(:(=), :j, c)
    v = Expr(:(=), :v, b)
    n = Expr(:(=), :n, Expr(:call, :<, :v, a))
    # m = Expr(:(=), b, Expr(:if, :n, a, b))
    m = Expr(:(=), b, Expr(:call, :ifelse, :n, a, :v))
    sym = i[N]
    l = :($sym - 1)
    for d = (N - 1):-1:1
        sym = i[d]
        l = :($sym - 1 + static_dims[$d] * $l)
    end
    l = :($l + 1)
    # cₐ = Expr(:(=), c, Expr(:if, :n, l, c))
    cₐ = Expr(:(=), c, Expr(:call, :ifelse, :n, l, :j))
    push!(body.args, j)
    push!(body.args, v)
    push!(body.args, n)
    push!(body.args, m)
    push!(body.args, cₐ)
    push!(ls.args, body)
    return quote
        static_dims::Vector{Int} = collect(size(A))
        @turbo $ls
        return B, C
    end
end

function lvfindmax2(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = fill(typemin(T), Dᴮ)
    C = zeros(Int, Dᴮ)
    Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
    _lvfindmax2!(B, C, A, Dᴮ′)
    B, CartesianIndices(A)[C]
end

@generated function _lvfindmax2!(B::AbstractArray{T, N}, C::AbstractArray{S, N},
                                 A::AbstractArray{T, N}, dims::D) where {T, N, S, D}
    findmax_quote2(N, D)
end

AA = rand(100, 10, 10);

@timev lvfindmax2(A, dims)
Meta.parse("true ? 5 : 10")

A = rand(10,10);
dims = (1,)

lvfindmax2(A, dims) == lvfindmax(A, dims) == findmax(A, dims=dims)
ls = findmax_quote2(N, typeof(Dᴮ′))
Meta.@dump(ifelse(x > 0, log(x), inv(x)))

function findmaxturbo!(B, C, A::AbstractMatrix{T}, static_dims) where {T}
    v = zero(T)
    j = 0
    l = ((1 - 1) + static_dims[1] * (1 - 1)) + 1
    ii = 1
    @turbo for i_2 = axes(A, 2), i_1 = axes(A, 1)
        v = B[ii, i_2]
        j = C[ii, i_2]
        # l = ((i_1 - 1) + static_dims[1] * (i_2 - 1)) + 1
        l = i_1 - 1 + static_dims[1] * (i_2 - 1) + 1
        newm = A[i_1, i_2] > B[ii, i_2]
        B[ii, i_2] = newm ? A[i_1, i_2] : v
        C[ii, i_2] = newm ? l : j
    end
    B, C
end
function findmaxturbo(A::AbstractMatrix{T}) where {T}
    B = fill(typemin(eltype(A)), Base.setindex(size(A), 1, 1))
    C = ones(Int, Base.setindex(size(A), 1, 1))
    static_dims = collect(size(A))
    findmaxturbo!(B, C, A, static_dims)
    B, CartesianIndices(A)[C]
end
A = rand(2, 3)
static_dims = size(A)
@timev findmaxturbo6!(B, C, A, collect(size(A)))
findmaxturbo(A)
using LoopVectorization
LoopVectorization.check_args(B, C, A)


text:
@testset "text tests" begin
  @testset "text roundtrip" begin
    # Strings and symbols
    @test (s = "abc"; x = K(s); String(x) == s)
    @test (s = :abc; x = K(s); Symbol(x) == s && String(x) == String(s))
    @test (a = [:abc, :def]; x = K(a); Array(x) == a)
  end
  @testset "Scalar to string" begin
    @test string(K(42)) == "K(42)"
    @test string(K(:a)) == "a"
  end
end


text:
"These fonctions are based on the package https://github.com/cormullion/TickTock.jl and are used to compute the elapsed time during a search."

"""
    tic()
Start a timer.
"""
function tic()
    t0 = time_ns()
    task_local_storage(:TIMERS, (t0, get(task_local_storage(), :TIMERS, ())))
end

"""
    peektimer()
Return the elapsed seconds counted by the most recent timer, without stopping it.
"""
function peektimer()
    t1 = time_ns()
    timers = get(task_local_storage(), :TIMERS, ())
    if timers === ()
        error("Use `tic()` to start a timer.")
    end
    t0 = timers[1]::UInt64
    return (t1 - t0)/1e9
end

function toc()
    timers = get(task_local_storage(), :TIMERS, ())
    if timers === ()
        error("Use `tick()` to start a timer.")
    end
    t = peektimer()
    task_local_storage(:TIMERS, timers[2])
    return t
end


text:
# Some of the function does not exist in Julia/LinearAlgebra
export isunitary, isreflexive, iscommute, ishermitian

"""
    isunitary(op) -> Bool

check if this operator is a unitary operator.
"""
isunitary(op) = op' * op ≈ IMatrix(size(op, 1))
isunitary(op::Number) = op' * op ≈ one(op)

"""
    isreflexive(op) -> Bool

check if this operator is reflexive.
"""
isreflexive(op) = op * op ≈ IMatrix(size(op, 1))
isreflexive(op::Number) = op * op ≈ one(op)

"""
    ishermitian(op) -> Bool

check if this operator is hermitian.
"""
LinearAlgebra.ishermitian(op) = op' ≈ op

"""
    iscommute(ops...) -> Bool

check if operators are commute.
"""
function iscommute(ops...)
    n = length(ops)
    for i=1:n
        for j=i+1:n
            iscommute(ops[i], ops[j]) || return false
        end
    end
    true
end

iscommute(op1, op2) = op1*op2 ≈ op2*op1


text:
module libIIO_jl

using Preferences
using Libdl
# Regarding type convertions :
#   - Basic types as arguments need to be passed as the corresponding Julia type to the wrapper function.
#   - Basic types as return values keep the Ctype alias and are not converted to "native" Julia type.
#   Ex: If the C function returns a size_t, the wrapper function returns a Csize_t (or long long -> Clonglong)
# This is done to keep the maximum transparency as to what the C function needs, does, and returns
# When the C function needs a pointer, a reference to a native Julia type is used.

"""
Change libIIO Ldriver provider. Support "yggdrasil" to use shipped jll file or "local" to use custom installed library
set_provider("yggdrasil")
or
set_provider("local")
"""
# We export here to be sure it is seen @AdalmPluto level
function set_provider(new_provider::String)
    if !(new_provider in ("yggdrasil", "local"))
        throw(ArgumentError("Invalid provider: \"$(new_provider)\""))
    end
    # Set it in our runtime values, as well as saving it to disk
    @set_preferences!("provider" => new_provider)
    @info("New provider set; restart your Julia session for this change to take effect!")
end
function get_provider()
    return @load_preference("provider","yggdrasil")
end
const libiio_provider = get_provider()
@static  if libiio_provider == "yggdrasil"
    # --- Using Yggdrasil jll file
    using libiio_jll
    const libIIO = libiio_jll.libiio
    # --- Defning rootpath
    const libIIO_rootpath = libiio_jll.artifact_dir
end
@static if libiio_provider == "local"
    # --- Using local install, assuming it works
    libIIO_system_h = dlopen("libiio", false);
    const libIIO = dlpath(libIIO_system_h)
end

# needed for libIIO functions
const BUF_SIZE = 2^12; # same value as iio_common.h
const C_INT_MAX = 2^31 - 1;
# disable Julia errors and return the NULL/error code instead
NO_ASSERT = false;

# adds a udev rule needed for usb devices
# should be a volatile rule and will need to be added each boot
# --> Pain in the ass, we put that in the /etc/udev/rules.d folder
# but it makes it possible to delete the artifact without leftovers
function __init__()
    # This workflow is only valid for Linux machines, as there is no rules for MACOS
    if Sys.islinux()
        # Now the udev_rules is in another function
        if !isfile("/etc/udev/rules.d/53-adi-plutosdr-usb.rules")
            # print a warning here to be sure it is run when local provider is used
            if get_provider() == "default" || get_provider() == "yggdrasil"
                println("On Linux with Yggdrasil, specific udev-rules should be added with root privilege. Be sure to run AdalmPluto.set_udev_rules() ")
            end
        end
    end
end

function set_udev_rules()
    if Sys.islinux()
        if get_provider() == "default" || get_provider() == "yggdrasil"
            # if !isfile("/run/udev/rules.d/90-libiio.rules")
            println("\nAdding the udev rules to /run/udev/rules.d/53-adi-plutosdr-usb.rules.\nDirectory is write protected, password prompt does not come from Julia");
            rule = """SUBSYSTEM=="usb", ATTRS{idVendor}=="0456", ATTRS{idProduct}=="b674", MODE="0664", GROUP="plugdev"
            SUBSYSTEM=="usb", ATTRS{idVendor}=="2fa2", ATTRS{idProduct}=="5a32", MODE="0664", GROUP="plugdev"
            # SDR Device
            SUBSYSTEM=="usb", ATTRS{idVendor}=="0456", ATTRS{idProduct}=="b673", MODE="0664", GROUP="plugdev"
            SUBSYSTEM=="usb", ATTRS{idVendor}=="2fa2", ATTRS{idProduct}=="5a02", MODE="0664", GROUP="plugdev"
            # tell the ModemManager (part of the NetworkManager suite) that the device is not a modem,
            # and don't send AT commands to it
            SUBSYSTEM=="usb", ATTRS{idVendor}=="0456", ATTRS{idProduct}=="b673", ENV{ID_MM_DEVICE_IGNORE}="1"
            SUBSYSTEM=="usb", ATTRS{idVendor}=="2fa2", ATTRS{idProduct}=="5a02", ENV{ID_MM_DEVICE_IGNORE}="1"\n""";
            open("/tmp/53-adi-plutosdr-usb.rules", "w") do f
                write(f, rule);
                end
                # run(`sudo mkdir -p /run/udev/rules.d`);
                run(`sudo cp /tmp/53-adi-plutosdr-usb.rules /etc/udev/rules.d/53-adi-plutosdr-usb.rules`);
            # else
                # println("Udev rules is already present in /run/rules.d/. Nothing to do")
            # end
        else
            println("Local libiio lib is used, so Libiio documentation to set the udev rules (should be done when intalling the lib")
        end
    end
end

# helpers
include("helpers.jl");
export
    toggleNoAssertions
;

# structures exports
include("structures.jl");
export
    iio_context_info,
    iio_scan_block,
    iio_scan_context,
    iio_context,
    iio_device,
    iio_channel,
    iio_buffer,
    iio_data_format,
    iio_modifier,                           # enum
    iio_chan_type                           # enum
;

# scan exports
include("scan.jl");
export
    C_iio_context_info_get_description,
    C_iio_context_info_get_uri,
    C_iio_context_info_list_free,
    C_iio_create_scan_block,
    C_iio_create_scan_context,
    C_iio_scan_block_destroy,
    C_iio_scan_block_get_info,
    C_iio_scan_block_scan,
    C_iio_scan_context_destroy,
    C_iio_scan_context_get_info_list
;

# toplevel exports
include("toplevel.jl");
export
    C_iio_get_backend,
    C_iio_get_backends_count,
    C_iio_has_backend,
    C_iio_library_get_version,
    C_iio_strerror
;

# context exports
include("context.jl");
export
    C_iio_context_clone,
    C_iio_context_destroy,
    C_iio_context_find_device,
    C_iio_context_get_attr,
    C_iio_context_get_attr_value,
    C_iio_context_get_attrs_count,
    C_iio_context_get_description,
    C_iio_context_get_device,
    C_iio_context_get_devices_count,
    C_iio_context_get_name,
    C_iio_context_get_version,
    C_iio_context_get_xml,
    C_iio_context_set_timeout,
    C_iio_create_context_from_uri,
    C_iio_create_default_context,
    C_iio_create_local_context,
    C_iio_create_network_context,
    C_iio_create_xml_context,
        C_iio_create_xml_context_mem
    ;

    # device exports
    include("device.jl");
    export
        C_iio_device_attr_read,
        C_iio_device_attr_read_all,             # PLACEHOLDER : function pointer needed
        C_iio_device_attr_read_bool,
        C_iio_device_attr_read_double,
        C_iio_device_attr_read_longlong,
        C_iio_device_attr_write,
        C_iio_device_attr_write_all,            # PLACEHOLDER : function pointer needed
        C_iio_device_attr_write_bool,
        C_iio_device_attr_write_double,
        C_iio_device_attr_write_longlong,
        C_iio_device_attr_write_raw,
        C_iio_device_buffer_attr_read,
        C_iio_device_buffer_attr_read_all,      # PLACEHOLDER : function pointer needed
        C_iio_device_buffer_attr_read_bool,
        C_iio_device_buffer_attr_read_double,
        C_iio_device_buffer_attr_read_longlong,
        C_iio_device_buffer_attr_write,
        C_iio_device_buffer_attr_write_all,     # PLACEHOLDER : function pointer needed
        C_iio_device_buffer_attr_write_bool,
        C_iio_device_buffer_attr_write_double,
        C_iio_device_buffer_attr_write_longlong,
        C_iio_device_buffer_attr_write_raw,
        C_iio_device_find_attr,
        C_iio_device_find_buffer_attr,
        C_iio_device_find_channel,
        C_iio_device_get_attr,
        C_iio_device_get_attrs_count,
        C_iio_device_get_buffer_attr,
        C_iio_device_get_buffer_attrs_count,
        C_iio_device_get_channel,
        C_iio_device_get_channels_count,
        C_iio_device_get_context,
        C_iio_device_get_data,
        C_iio_device_get_id,
        C_iio_device_get_name,
        C_iio_device_get_trigger,
        C_iio_device_is_trigger,
        C_iio_device_set_data,
        C_iio_device_set_kernel_buffers_count,
        C_iio_device_set_trigger
;

# channel exports
include("channel.jl");
export
    C_iio_channel_attr_get_filename,
    C_iio_channel_attr_read,
    C_iio_channel_attr_read_all,            # PLACEHOLDER : function pointer needed
    C_iio_channel_attr_read_bool,
    C_iio_channel_attr_read_double,
    C_iio_channel_attr_read_longlong,
    C_iio_channel_attr_write,
    C_iio_channel_attr_write_all,           # PLACEHOLDER : function pointer needed
    C_iio_channel_attr_write_bool,
    C_iio_channel_attr_write_double,
    C_iio_channel_attr_write_longlong,
    C_iio_channel_attr_write_raw,
    C_iio_channel_disable,
    C_iio_channel_enable,
    C_iio_channel_find_attr,
    C_iio_channel_get_attr,
    C_iio_channel_get_attrs_count,
    C_iio_channel_get_data,
    C_iio_channel_get_device,
    C_iio_channel_get_id,
    C_iio_channel_get_modifier,
    C_iio_channel_get_name,
    C_iio_channel_get_type,
    C_iio_channel_is_enabled,
    C_iio_channel_is_output,
    C_iio_channel_is_scan_element,
    #  C_iio_channel_read,
    C_iio_channel_read!,
    C_iio_channel_read_raw,                 # PLACEHOLDER : data format
    C_iio_channel_set_data,
    C_iio_channel_write,                    # PLACEHOLDER : data format
    C_iio_channel_write_raw                 # PLACEHOLDER : data format
;

# buffer exports
include("buffer.jl");
export
    C_iio_buffer_cancel,
    C_iio_buffer_destroy,
    C_iio_buffer_end,
    C_iio_buffer_first,
    C_iio_buffer_foreach_sample,            # PLACEHOLDER : function pointer needed
    C_iio_buffer_get_data,
    C_iio_buffer_get_device,
    C_iio_buffer_get_poll_fd,
    C_iio_buffer_push,
    C_iio_buffer_push_partial,
    C_iio_buffer_refill,
    C_iio_buffer_set_blocking_mode,
    C_iio_buffer_set_data,
    C_iio_buffer_start,
    C_iio_buffer_step,
    C_iio_device_create_buffer
;

# debug exports
include("debug.jl");
export
    C_iio_device_get_sample_size,
    C_iio_device_identify_filename
;

end


text:
using LinearAlgebra

using Oceananigans.Architectures: array_type

function can_solve_single_tridiagonal_system(arch, N)
    ArrayType = array_type(arch)

    a = rand(N-1)
    b = 3 .+ rand(N) # +3 to ensure diagonal dominance.
    c = rand(N-1)
    f = rand(N)

    # Solve the system with backslash on the CPU to avoid scalar operations on the GPU.
    M = Tridiagonal(a, b, c)
    ϕ_correct = M \ f

    # Convert to CuArray if needed.
    a, b, c, f = ArrayType.([a, b, c, f])

    ϕ = reshape(zeros(N), (1, 1, N)) |> ArrayType

    grid = RectilinearGrid(arch, size=(1, 1, N), extent=(1, 1, 1))

    btsolver = BatchedTridiagonalSolver(arch, grid;
                                        lower_diagonal = a,
                                        diagonal = b,
                                        upper_diagonal = c)

    solve!(ϕ, btsolver, f)

    return Array(ϕ[:]) ≈ ϕ_correct
end

function can_solve_single_tridiagonal_system_with_functions(arch, N)
    ArrayType = array_type(arch)

    grid = RectilinearGrid(arch, size=(1, 1, N), extent=(1, 1, 1))

    a = rand(N-1)
    c = rand(N-1)

    @inline b(i, j, k, grid) = 3 .+ cos(2π * grid.zᵃᵃᶜ[k])  # +3 to ensure diagonal dominance.
    @inline f(i, j, k, grid) = sin(2π * grid.zᵃᵃᶜ[k])

    bₐ = [b(1, 1, k, grid) for k in 1:N]
    fₐ = [f(1, 1, k, grid) for k in 1:N]

    # Solve the system with backslash on the CPU to avoid scalar operations on the GPU.
    M = Tridiagonal(a, bₐ, c)
    ϕ_correct = M \ fₐ

    # Convert to CuArray if needed.
    a, c = ArrayType.((a, c))

    ϕ = reshape(zeros(N), (1, 1, N)) |> ArrayType

    btsolver = BatchedTridiagonalSolver(arch, grid;
                                        lower_diagonal = a,
                                        diagonal = b,
                                        upper_diagonal = c)

    solve!(ϕ, btsolver, f)

    return Array(ϕ[:]) ≈ ϕ_correct
end

function can_solve_batched_tridiagonal_system_with_3D_RHS(arch, Nx, Ny, Nz)
    ArrayType = array_type(arch)

    a = rand(Nz-1)
    b = 3 .+ rand(Nz) # +3 to ensure diagonal dominance.
    c = rand(Nz-1)
    f = rand(Nx, Ny, Nz)

    M = Tridiagonal(a, b, c)
    ϕ_correct = zeros(Nx, Ny, Nz)

    # Solve the systems with backslash on the CPU to avoid scalar operations on the GPU.
    for i = 1:Nx, j = 1:Ny
        ϕ_correct[i, j, :] .= M \ f[i, j, :]
    end

    # Convert to CuArray if needed.
    a, b, c, f = ArrayType.([a, b, c, f])

    grid = RectilinearGrid(arch, size=(Nx, Ny, Nz), extent=(1, 1, 1))
    btsolver = BatchedTridiagonalSolver(arch, grid;
                                        lower_diagonal = a,
                                        diagonal = b,
                                        upper_diagonal = c)

    ϕ = zeros(Nx, Ny, Nz) |> ArrayType

    solve!(ϕ, btsolver, f)

    return Array(ϕ) ≈ ϕ_correct
end

function can_solve_batched_tridiagonal_system_with_3D_functions(arch, Nx, Ny, Nz)
    ArrayType = array_type(arch)

    grid = RectilinearGrid(arch, size=(Nx, Ny, Nz), extent=(1, 1, 1))

    a = rand(Nz-1)
    c = rand(Nz-1)

    @inline b(i, j, k, grid) = 3 + grid.xᶜᵃᵃ[i] * grid.yᵃᶜᵃ[j] * cos(2π * grid.zᵃᵃᶜ[k])
    @inline f(i, j, k, grid) = (grid.xᶜᵃᵃ[i] + grid.yᵃᶜᵃ[j]) * sin(2π * grid.zᵃᵃᶜ[k])

    ϕ_correct = zeros(Nx, Ny, Nz)

    # Solve the system with backslash on the CPU to avoid scalar operations on the GPU.
    for i = 1:Nx, j = 1:Ny
        bₐ = [b(i, j, k, grid) for k in 1:Nz]
        M = Tridiagonal(a, bₐ, c)

        fₐ = [f(i, j, k, grid) for k in 1:Nz]
        ϕ_correct[i, j, :] .= M \ fₐ
    end

    # Convert to CuArray if needed.
    a, c = ArrayType.([a, c])

    btsolver = BatchedTridiagonalSolver(arch, grid;
                                        lower_diagonal = a,
                                        diagonal = b,
                                        upper_diagonal = c)

    ϕ = zeros(Nx, Ny, Nz) |> ArrayType
    solve!(ϕ, btsolver, f)

    return Array(ϕ) ≈ ϕ_correct
end

@testset "Batched tridiagonal solvers" begin
    @info "Testing BatchedTridiagonalSolver..."

    for arch in archs
        @testset "Batched tridiagonal solver [$arch]" begin
            for Nz in [8, 11, 18]
                @test can_solve_single_tridiagonal_system(arch, Nz)
                @test can_solve_single_tridiagonal_system_with_functions(arch, Nz)
            end

            for Nx in [3, 8], Ny in [5, 16], Nz in [8, 11]
                @test can_solve_batched_tridiagonal_system_with_3D_RHS(arch, Nx, Ny, Nz)
                @test can_solve_batched_tridiagonal_system_with_3D_functions(arch, Nx, Ny, Nz)
            end
        end
    end
end