text:
"""
Convolution for CTMC
"""

function _dger!(x::Array{Tv,1}, y::Array{Tv,1}, A::Matrix{Tv}) where {Tv}
    m, n = size(A)
    @inbounds for j = 1:n
        for i = 1:m
            A[i,j] += x[i] * y[j]
        end
    end
    nothing
end

function _dger!(x::Array{Tv,1}, y::Array{Tv,1}, A::SparseCSR{Tv,Ti}) where {Ti,Tv}
    m, n = size(A)
    @inbounds for i = 1:m
        for z = A.rowptr[i]:A.rowptr[i+1]-1
            j = A.colind[z]
            A.val[z] += x[i] * y[j]
        end
    end
    nothing
end

function _dger!(x::Array{Tv,1}, y::Array{Tv,1}, A::SparseCSC{Tv,Ti}) where {Ti,Tv}
    m, n = size(A)
    @inbounds for j = 1:n
        for z = A.colptr[j]:A.colptr[j+1]-1
            i = A.rowind[z]
            A.val[z] += x[i] * y[j]
        end
    end
    nothing
end

function _dger!(x::Array{Tv,1}, y::Array{Tv,1}, A::SparseCOO{Tv,Ti}) where {Ti,Tv}
    @inbounds for z = 1:nnz(A)
        i = A.rowind[z]
        j = A.colind[z]
        A.val[z] += x[i] * y[j]
    end
    nothing
end

function _dger!(x::Array{Tv,1}, y::Array{Tv,1}, A::SparseMatrixCSC{Tv,Ti}) where {Ti,Tv}
    m, n = size(A)
    @inbounds for j = 1:n
        for z = A.colptr[j]:A.colptr[j+1]-1
            i = A.rowval[z]
            A.nzval[z] += x[i] * y[j]
        end
    end
    nothing
end

"""
convunifstep!(trQ, trH, P, poi, range, weight, qv_weight, x, y, z, H)

! Description: convolution integral operation for matrix exp form;
!
!          |t
! trH(H) = | exp(trQ(Q)*s) * x * y' * exp(trQ(Q)*(t-s)) ds
!          |0
!
!        and
!
!        z = exp(trQ(Q)*t) * x
!
!        t is involved in the Poisson probability vector.
!        qv is an uniformed parameter
!        return value is z

Parameters:
- trQ: transpose operator
- trH: transpose operator
- P: The uniformed matrix
- poi: Poisson p.m.f.
- range: domain of Poisson p.m.f
- weight: The normalizing constant for Poisson p.m.f.
- qv_weight: The normalizing constant for Poisson c.c.d.f.
- x: Vector (inout). x may be changed after executing
- y: Vector (out). y should be zero before executing
- z: Vector (out). cy should be zero before executing
- H: Array (out). H should be zero before executing
Return value: nothing
"""

function convunifstep!(trQ::Symbol, trH::Symbol,
    P::AbstractMatrix{Tv},
    poi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv, qv_weight::Tv,
    x::Array{Tv,N}, y::Array{Tv,N}, z::Array{Tv,N},
    H::AbstractMatrix{Tv})::Nothing where {Ti,Tv,N}
    _convunifstep!(Val(trQ), Val(trH), P, poi, range, weight, qv_weight, x, y, z, H)
end

@origin (vc => left, poi => left) function _convunifstep!(::Val{:N}, ::Val{:N},
    P::AbstractMatrix{Tv},
    poi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv, qv_weight::Tv,
    x::Array{Tv,N}, y::Array{Tv,N}, z::Array{Tv,N},
    H::AbstractMatrix{Tv})::Nothing where {Ti,Tv,N}
    @inbounds begin
        left, right = range
        Pdash = P'
        vc = Vector{Vector{Tv}}(undef, right - left + 1)
        vc[right] = zero(x)
        @axpy(poi[right], y, vc[right])
        for l = right-1:-1:left+1
            vc[l] = Pdash * vc[l+1]
            @axpy(poi[l], y, vc[l])
        end

        @axpy(poi[left], x, z)
        _dger!(x, vc[left+1], H)
        for l = left+1:right-1
            x .= P * x
            @axpy(poi[l], x, z)
            _dger!(x, vc[l+1], H)
        end
        @scal(1/weight, z)
        @scal(1/qv_weight, H)
        nothing
    end
end

@origin (vc => left, poi => left) function _convunifstep!(::Val{:T}, ::Val{:N},
    P::AbstractMatrix{Tv},
    poi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv, qv_weight::Tv,
    x::Array{Tv,N}, y::Array{Tv,N}, z::Array{Tv,N},
    H::AbstractMatrix{Tv})::Nothing where {Ti,Tv,N}
    @inbounds begin
        left, right = range
        Pdash = P'
        vc = Vector{Vector{Tv}}(undef, right - left + 1)
        vc[right] = zero(x)
        @axpy(poi[right], y, vc[right])
        for l = right-1:-1:left+1
            vc[l] = P * vc[l+1]
            @axpy(poi[l], y, vc[l])
        end

        @axpy(poi[left], x, z)
        _dger!(x, vc[left+1], H)
        for l = left+1:right-1
            x .= Pdash * x
            @axpy(poi[l], x, z)
            _dger!(x, vc[l+1], H)
        end
        @scal(1/weight, z)
        @scal(1/qv_weight, H)
        nothing
    end
end

@origin (vc => left, poi => left) function _convunifstep!(::Val{:N}, ::Val{:T},
    P::AbstractMatrix{Tv},
    poi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv, qv_weight::Tv,
    x::Array{Tv,N}, y::Array{Tv,N}, z::Array{Tv,N},
    H::AbstractMatrix{Tv})::Nothing where {Ti,Tv,N}
    @inbounds begin
        left, right = range
        Pdash = P'
        vc = Vector{Vector{Tv}}(undef, right - left + 1)
        vc[right] = zero(x)
        @axpy(poi[right], y, vc[right])
        for l = right-1:-1:left+1
            vc[l] = Pdash * vc[l+1]
            @axpy(poi[l], y, vc[l])
        end

        @axpy(poi[left], x, z)
        _dger!(vc[left+1], x, H)
        for l = left+1:right-1
            x .= P * x
            @axpy(poi[l], x, z)
            _dger!(vc[l+1], x, H)
        end
        @scal(1/weight, z)
        @scal(1/qv_weight, H)
        nothing
    end
end

@origin (vc => left, poi => left) function _convunifstep!(::Val{:T}, ::Val{:T},
    P::AbstractMatrix{Tv},
    poi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv, qv_weight::Tv,
    x::Array{Tv,N}, y::Array{Tv,N}, z::Array{Tv,N},
    H::AbstractMatrix{Tv})::Nothing where {Ti,Tv,N}
    @inbounds begin
        left, right = range
        Pdash = P'
        vc = Vector{Vector{Tv}}(undef, right - left + 1)
        vc[right] = zero(x)
        @axpy(poi[right], y, vc[right])
        for l = right-1:-1:left+1
            vc[l] = P * vc[l+1]
            @axpy(poi[l], y, vc[l])
        end

        @axpy(poi[left], x, z)
        _dger!(vc[left+1], x, H)
        for l = left+1:right-1
            x .= Pdash * x
            @axpy(poi[l], x, z)
            _dger!(vc[l+1], x, H)
        end
        @scal(1/weight, z)
        @scal(1/qv_weight, H)
        nothing
    end
end


text:
module Day6

export get_inputs, get_solution1, get_solution2

## Input getting
function get_inputs()
    test_input1 = test_input2 = read_inputs("test_input1.txt")
    test_output1 = 11
    test_output2 = 6
    data = read_inputs("input.txt")
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

read_inputs(f_name) = parse_declaration.(split(read(joinpath(@__DIR__, f_name), String), r"\n\n|\r\n\r\n"))

parse_declaration(str) = split(str, r"\n|\r\n", keepempty=false)

## Solution functions
# Part 1
get_solution1(data) = sum(length.(reduce.(union, data)))

# Part 2
get_solution2(data) = sum(length.(reduce.(intersect, data)))

end

ext: jl


text:
"""
    binary_pivot_search(array::Array{Int64,1}, target::Int64)

The idea is to find the pivot point for finding the target in an unsorted array. For this
reason, the array has to be divided into two subarrays; a binary search is performed on
the subarrays.

# Arguments
- `array::Array{Int64,1}`: Unsorted array of integers
- `target::Int64`: Target-value to find the position

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [2, 5, 4, 7, 2, 8, 9, 3, 10, 2]
julia> target = 3
julia> ClassicAlgorithmsCollections.binary_pivot_search(arr, target)
8
```
"""
function binary_pivot_search(array::Array{Int64,1}, target::Int64)
    left = 1
    right = length(array)

    while left <= right
        # Generating the middle element position
        mid = fld((left + right), 2)

        # Check if target is hit
        if array[mid] == target
            return mid
            # Check if array is sorted from left to mid
        elseif array[left] <= array[mid]
            # Check if key lies in the left or right half
            if target >= array[left] && target <= array[mid]
                right = mid - 1
            else
                left = mid + 1
            end
            # Check if array is sorted from mid to the right
        elseif array[right] >= array[mid]
            # If array from left to mid is not sorted, array from mid to right is sorted
            if target >= array[mid] && target <= array[right]
                left = mid + 1
            else
                right = mid - 1
            end
        end
    end
    return nothing
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -12.75118444054646 <= q <= 9.039664285972906)

                     add_NL_constraint(m, :(log(1 + exp(0.715796138346545 + -0.1516695982843861*log(1 + exp(-0.317423658088571 + -0.8843064794993589*log(1 + exp(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3]))) + 0.6219431368805965*log(1 + exp(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3]))) + -0.24757066957523488*log(1 + exp(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))))) + -0.0820632819191216*log(1 + exp(-0.8960501996212993 + -0.7682728024607406*log(1 + exp(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3]))) + -0.3378399703859083*log(1 + exp(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3]))) + 0.4153867033969405*log(1 + exp(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))))) + -0.0571486344530765*log(1 + exp(0.054021626349418206 + -0.05819989275612114*log(1 + exp(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3]))) + 0.8875139231924525*log(1 + exp(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3]))) + 0.11280184637312685*log(1 + exp(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))))))) + log(1 + exp(-0.8464220616292248 + -0.9331432460947604*log(1 + exp(-0.317423658088571 + -0.8843064794993589*log(1 + exp(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3]))) + 0.6219431368805965*log(1 + exp(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3]))) + -0.24757066957523488*log(1 + exp(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))))) + 0.425513343659472*log(1 + exp(-0.8960501996212993 + -0.7682728024607406*log(1 + exp(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3]))) + -0.3378399703859083*log(1 + exp(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3]))) + 0.4153867033969405*log(1 + exp(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))))) + 0.1421259709107252*log(1 + exp(0.054021626349418206 + -0.05819989275612114*log(1 + exp(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3]))) + 0.8875139231924525*log(1 + exp(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3]))) + 0.11280184637312685*log(1 + exp(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))))))) + log(1 + exp(-0.5205460270098605 + 0.13812073353255405*log(1 + exp(-0.317423658088571 + -0.8843064794993589*log(1 + exp(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3]))) + 0.6219431368805965*log(1 + exp(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3]))) + -0.24757066957523488*log(1 + exp(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))))) + 0.7749706249828683*log(1 + exp(-0.8960501996212993 + -0.7682728024607406*log(1 + exp(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3]))) + -0.3378399703859083*log(1 + exp(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3]))) + 0.4153867033969405*log(1 + exp(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))))) + 0.5007112790542019*log(1 + exp(0.054021626349418206 + -0.05819989275612114*log(1 + exp(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3]))) + 0.8875139231924525*log(1 + exp(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3]))) + 0.11280184637312685*log(1 + exp(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
### A Pluto.jl notebook ###
# v0.12.17

using Markdown
using InteractiveUtils

# ╔═╡ a10105d6-4ae0-11eb-0192-2f51f6c7cf2b
begin
	using Flux
	using Gadfly
end

# ╔═╡ fb607214-4ae0-11eb-1642-277331455b29
begin #household functions and variables
	model(X) = W*X .+ b
	loss(W,b,X,y) = sum(((W.-X).^2).+((b.-y).^2))#sum((model(x) .- y).^2)

	function init(n::Int)
	W = rand(n)
	b = zeros(n)
	X, y = rand(n), rand(n)
		return W,b,X,y
	end
	W,b,X,y=init(100)
end

# ╔═╡ b17b71c6-4ae0-11eb-2693-7dc486447d0e
begin
	d = gradient(() -> loss(W,b,X,y), params(W, b))
	dW,db = d[W],d[b]

	W .-= 0.01 * dW
	b .-= 0.01 * db

	l=loss(W,b,X,y)

	plot(layer(x=X,y=y),layer(x=W,y=b,Geom.line))
end

# ╔═╡ 237e0690-4ae5-11eb-3c50-c1be36fa96ef
l

# ╔═╡ Cell order:
# ╠═a10105d6-4ae0-11eb-0192-2f51f6c7cf2b
# ╠═fb607214-4ae0-11eb-1642-277331455b29
# ╠═237e0690-4ae5-11eb-3c50-c1be36fa96ef
# ╠═b17b71c6-4ae0-11eb-2693-7dc486447d0e

max_stars_repo_name: punkmath/Misc

max_issues_repo_name: punkmath/Misc

max_forks_repo_name: punkmath/Misc


text:
function view_profile(; C = false, kwargs...)
    d = Dict()

    if VERSION >= v"1.8.0-DEV.460"
        threads = ["all", 1:Threads.nthreads()...]
    else
        threads = ["all"]
    end
    data = Profile.fetch()

    if isempty(data)
        Profile.warning_empty()
        return
    end

    lidict = Profile.getdict(unique(data))
    data_u64 = convert(Vector{UInt64}, data)
    for thread in threads
        graph = stackframetree(data_u64, lidict; thread = thread, kwargs...)
        d[thread] = dicttree(Dict(
            :func => "root",
            :file => "",
            :path => "",
            :line => 0,
            :count => graph.count,
            :flags => 0x0,
            :children => []
        ), graph; C = C, kwargs...)
    end

    JSONRPC.send(conn_endpoint[], repl_showprofileresult_notification_type, (; trace = d))
end

function stackframetree(data_u64, lidict; thread = nothing, combine = true, recur = :off)
    root = combine ? Profile.StackFrameTree{Profile.StackFrame}() : Profile.StackFrameTree{UInt64}()
    if VERSION >= v"1.8.0-DEV.460"
        thread = thread == "all" ? (1:Threads.nthreads()) : thread
        root, _ = Profile.tree!(root, data_u64, lidict, true, recur, thread)
    else
        root = Profile.tree!(root, data_u64, lidict, true, recur)
    end
    if !isempty(root.down)
        root.count = sum(pr -> pr.second.count, root.down)
    end

    return root
end

# https://github.com/timholy/FlameGraphs.jl/blob/master/src/graph.jl
const runtime_dispatch = UInt8(2^0)
const gc_event         = UInt8(2^1)
const repl             = UInt8(2^2)
const compilation      = UInt8(2^3)
const task_event       = UInt8(2^4)
# const              = UInt8(2^5)
# const              = UInt8(2^6)
# const              = UInt8(2^7)
# const              = UInt8(2^8)

function status(sf::Profile.StackFrame)
    st = UInt8(0)
    if sf.from_c && (sf.func === :jl_invoke || sf.func === :jl_apply_generic || sf.func === :ijl_apply_generic)
        st |= runtime_dispatch
    end
    if sf.from_c && startswith(String(sf.func), "jl_gc_")
        st |= gc_event
    end
    if !sf.from_c && sf.func === :eval_user_input && endswith(String(sf.file), "REPL.jl")
        st |= repl
    end
    if !sf.from_c && occursin("./compiler/", String(sf.file))
        st |= compilation
    end
    if !sf.from_c && occursin("task.jl", String(sf.file))
        st |= task_event
    end
    return st
end

function status(node::Profile.StackFrameTree, C::Bool)
    st = status(node.frame)
    C && return st
    # If we're suppressing C frames, check all C-frame children
    for child in values(node.down)
        child.frame.from_c || continue
        st |= status(child, C)
    end
    return st
end

function add_child(graph, node, C::Bool)
    name = string(node.frame.file)
    func = String(node.frame.func)

    if func == ""
        func = "unknown"
    end

    d = Dict(
        :func => func,
        :file => basename(name),
        :path => fullpath(name),
        :line => node.frame.line,
        :count => node.count,
        :flags => status(node, C),
        :children => []
    )
    push!(graph[:children], d)

    return d
end

function dicttree(graph, node::Profile.StackFrameTree; C = false)
    for child_node in sort!(collect(values(node.down)); rev = true, by = node -> node.count)
        # child not a hidden frame
        if C || !child_node.frame.from_c
            child = add_child(graph, child_node, C)
            dicttree(child, child_node; C = C)
        else
            dicttree(graph, child_node)
        end
    end

    return graph
end

"""
    @profview f(args...) [C = false]

Clear the Profile buffer, profile `f(args...)`, and view the result graphically.

The default of `C = false` will only show Julia frames in the profile graph.
"""
macro profview(ex, args...)
    return quote
        Profile.clear()
        Profile.@profile $(esc(ex))
        view_profile(; $(esc.(args)...))
    end
end

max_stars_repo_name: TanookiToad/julia-vscode

max_issues_repo_name: TanookiToad/julia-vscode

max_forks_repo_name: TanookiToad/julia-vscode


text:
# Benchmark Test for YaoArrayRegister
# NOTE: we only test the standard interface here:
# 1. locations are all tuples
#
# forwarded bindings should be tested in YaoBlocks via blocks.

using PkgBenchmark, BenchmarkTools
using YaoArrayRegister, BitBasis, Random, YaoBase, StaticArrays, LuxurySparse
using LinearAlgebra, SparseArrays

bench(n, U, loc::Tuple) =
    @benchmarkable instruct!(st, $U, $loc) setup = (st = statevec(rand_state($n)))
bench(n, U, loc::Tuple, control_locs::Tuple, control_bits::Tuple) =
    @benchmarkable instruct!(st, $U, $loc, $control_locs, $control_bits) setup = (st = statevec(rand_state($n)))
bench(n, B::Int, U, loc::Tuple) =
    @benchmarkable instruct!(st, $U, $loc) setup = (st = statevec(rand_state($n, nbatch = $B)))
bench(n, B::Int, U, loc::Tuple, control_locs::Tuple, control_bits::Tuple) =
    @benchmarkable instruct!(st, $U, $loc, $control_locs, $control_bits) setup = (st = statevec(rand_state(
        $n,
        nbatch = $B,
    )))

const SUITE = BenchmarkGroup()
SUITE["specialized"] = BenchmarkGroup()

@info "generating benchmark for specialized operators for single qubits"
# Specialized Gate Instruction
SUITE["specialized"]["single qubit"] = BenchmarkGroup()
## single qubit benchmark
for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 1:4:25
    SUITE["specialized"]["single qubit"][string(U), n] = bench(n, Val(U), (1,))
end

SUITE["specialized"]["single control"] = BenchmarkGroup()
for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 4:4:25
    SUITE["specialized"]["single control"][string(U), n, (2,), (1,)] = bench(
        n,
        Val(U),
        (1,),
        (2,),
        (1,),
    )
end

SUITE["specialized"]["multi control"] = BenchmarkGroup()

for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 4:4:25
    control_locs = Tuple(2:n)
    control_bits = ntuple(x -> 1, n - 1)
    SUITE["specialized"]["multi control"][string(U), n, 2:n, control_locs] = bench(
        n,
        Val(U),
        (1,),
        control_locs,
        control_bits,
    )
end

SUITE["specialized"]["multi qubit"] = BenchmarkGroup()
const location_sparsity = 0.4
@info "generating benchmark for specialized operators for multi qubits"
## multi qubit benchmark
for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 4:4:25
    perms = randperm(n)[1:ceil(Int, location_sparsity * n)]
    SUITE["specialized"]["multi qubit"][string(U), n] = bench(n, Val(U), Tuple(perms))
end

SUITE["specialized"]["multi qubit multi control"] = BenchmarkGroup()
SUITE["specialized"]["single qubit multi control"] = BenchmarkGroup()

const control_rate = 0.3
for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 4:4:25
    num_controls = ceil(Int, n * control_rate)
    perms = randperm(n)
    control_locs = Tuple(perms[1:num_controls])
    control_bits = ntuple(x -> rand(0:1), num_controls)
    perms = perms[num_controls+1:num_controls+round(Int, location_sparsity * n)]

    SUITE["specialized"]["multi qubit multi control"][string(U), n, num_controls] = bench(
        n,
        Val(U),
        Tuple(perms),
        control_locs,
        control_bits,
    )
    SUITE["specialized"]["single qubit multi control"][string(U), n, num_controls] = bench(
        n,
        Val(U),
        (perms[1],),
        control_locs,
        control_bits,
    )
end

for n in 4:4:25
    SUITE["specialized"]["multi qubit"]["SWAP", n] = bench(n, Val(:SWAP), (1, 2))
    SUITE["specialized"]["multi qubit"]["SWAP", "random", n] = bench(
        n,
        Val(:SWAP),
        Tuple(randperm(n)[1:2]),
    )
end

# General Instructions (matrices based)
SUITE["matrices"] = BenchmarkGroup()
SUITE["matrices"]["contiguous"] = BenchmarkGroup()
SUITE["matrices"]["contiguous"]["ordered"] = BenchmarkGroup()
SUITE["matrices"]["contiguous"]["random"] = BenchmarkGroup()

SUITE["matrices"]["in-contiguous"] = BenchmarkGroup()
SUITE["matrices"]["in-contiguous"]["ordered"] = BenchmarkGroup()
SUITE["matrices"]["in-contiguous"]["random"] = BenchmarkGroup()

SUITE["matrices"]["single qubit"] = BenchmarkGroup()
SUITE["matrices"]["single qubit"]["ordered"] = BenchmarkGroup()
SUITE["matrices"]["single qubit"]["random"] = BenchmarkGroup()

## General Matrix Instruction
function matrices(::Type{T}, N) where {T}
    list = Any[
        rand_unitary(T, N), # dense matrices
        # SparseMatrixCSC
        sprand_hermitian(T, N, 0.1),
        # PermMatrix
        pmrand(T, N),
        Diagonal(rand(T, N)),
    ]
    if N < 100
        # StaticArrays
        push!(list, @SArray(rand(T, N, N)))
        push!(list, @MArray(rand(T, N, N)))
    end
    return list
end

# default test type is ComplexF64
matrices(N) = matrices(ComplexF64, N)

@info "generating benchmark for contiguous matrices locs"
### contiguous
for n in 1:2:10, T in [ComplexF64], U in matrices(T, 1 << n)
    # contiguous ordered address
    SUITE["matrices"]["contiguous"]["ordered"][n, string(T), string(typeof(U))] = bench(
        n,
        U,
        Tuple(1:n),
    )
    # contiguous random address
    SUITE["matrices"]["contiguous"]["random"][n, string(T), string(typeof(U))] = bench(
        n,
        U,
        Tuple(randperm(n)),
    )
end

@info "generating benchmark for in-contiguous matrices locs"
### in-contiguous
for m in 1:3, T in [ComplexF64], U in matrices(T, 1 << m)
    n = 10
    N = 1 << n
    # in-contiguous ordered address
    SUITE["matrices"]["in-contiguous"]["ordered"][m, string(T), string(typeof(U))] = bench(
        n,
        U,
        Tuple(sort(randperm(n)[1:m])),
    )
    # in-contiguous random address
    SUITE["matrices"]["in-contiguous"]["random"][m, string(T), string(typeof(U))] = bench(
        n,
        U,
        Tuple(randperm(n)[1:m]),
    )
end

@info "generating benchmark for single qubit matrices"
### single qubit
for T in [ComplexF64], U in matrices(T, 2), n in 1:4:25
    SUITE["matrices"]["single qubit"]["ordered"][string(T), string(typeof(U)), n] = bench(
        n,
        U,
        (rand(1:n),),
    )
    SUITE["matrices"]["single qubit"]["random"][string(T), string(typeof(U)), n] = bench(
        n,
        U,
        (rand(1:n),),
    )
end

SUITE["matrices"]["controlled"] = BenchmarkGroup()
SUITE["matrices"]["controlled"]["ordered"] = BenchmarkGroup()
SUITE["matrices"]["controlled"]["random"] = BenchmarkGroup()

test_bench(n, U, loc, control_locs, control_bits) =
    instruct!(statevec(rand_state(n)), U, loc, control_locs, control_bits)
for T in [ComplexF64], n in 1:2:10, m in 3:2:6, U in matrices(T, 1 << n)
    SUITE["matrices"]["controlled"]["ordered"][n, m, string(T), string(typeof(U))] = bench(
        n + m,
        U,
        Tuple(1:n),
        Tuple(n+1:n+m),
        ntuple(x -> 1, m),
    )

    perms = randperm(n + m)
    SUITE["matrices"]["controlled"]["random"][n, m, string(T), string(typeof(U))] = bench(
        n + m,
        U,
        Tuple(perms[1:n]),
        Tuple(perms[n+1:n+m]),
        ntuple(x -> 1, m),
    )
end

# batched
SUITE["batched specialized"] = BenchmarkGroup()

@info "generating benchmark for batched specialized operators for single qubits"
# Specialized Gate Instruction
SUITE["batched specialized"]["single qubit"] = BenchmarkGroup()
## single qubit benchmark
for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 5:4:15, B in 10:20:60
    SUITE["batched specialized"]["single qubit"][string(U), n, B] = bench(n, B, Val(U), (1,))
end

SUITE["batched specialized"]["single control"] = BenchmarkGroup()
for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 5:4:15, B in 10:20:60
    SUITE["batched specialized"]["single control"][string(U), n, B, (2,), (1,)] = bench(
        n,
        B,
        Val(U),
        (1,),
        (2,),
        (1,),
    )
end

SUITE["batched specialized"]["multi control"] = BenchmarkGroup()

for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 5:4:15, B in 10:20:60
    control_locs = Tuple(2:n)
    control_bits = ntuple(x -> 1, n - 1)
    SUITE["batched specialized"]["multi control"][string(U), n, B, 2:n, control_locs] = bench(
        n,
        B,
        Val(U),
        (1,),
        control_locs,
        control_bits,
    )
end

SUITE["batched specialized"]["multi qubit"] = BenchmarkGroup()
const location_sparsity = 0.4
@info "generating benchmark for specialized operators for multi qubits"
## multi qubit benchmark
for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 5:4:15, B in 10:20:60
    perms = randperm(n)[1:ceil(Int, location_sparsity * n)]
    SUITE["batched specialized"]["multi qubit"][string(U), n, B] = bench(n, B, Val(U), Tuple(perms))
end

SUITE["batched specialized"]["multi qubit multi control"] = BenchmarkGroup()
SUITE["batched specialized"]["single qubit multi control"] = BenchmarkGroup()

const control_rate = 0.3
for U in YaoArrayRegister.SPECIALIZATION_LIST, n in 5:4:15, B in 10:20:60
    num_controls = ceil(Int, n * control_rate)
    perms = randperm(n)
    control_locs = Tuple(perms[1:num_controls])
    control_bits = ntuple(x -> rand(0:1), num_controls)
    perms = perms[num_controls+1:num_controls+round(Int, location_sparsity * n)]

    SUITE["batched specialized"]["multi qubit multi control"][string(U), n, B, num_controls] = bench(
        n,
        B,
        Val(U),
        Tuple(perms),
        control_locs,
        control_bits,
    )
    SUITE["batched specialized"]["single qubit multi control"][string(U), n, B, num_controls] = bench(
        n,
        B,
        Val(U),
        (perms[1],),
        control_locs,
        control_bits,
    )
end

for n in 5:4:15, B in 10:20:60
    SUITE["batched specialized"]["multi qubit"]["SWAP", n, B] = bench(n, B, Val(:SWAP), (1, 2))
    SUITE["batched specialized"]["multi qubit"]["SWAP", "random", n, B] = bench(
        n,
        B,
        Val(:SWAP),
        Tuple(randperm(n)[1:2]),
    )
end

# General Instructions (matrices based)
SUITE["batched matrices"] = BenchmarkGroup()
SUITE["batched matrices"]["contiguous"] = BenchmarkGroup()
SUITE["batched matrices"]["contiguous"]["ordered"] = BenchmarkGroup()
SUITE["batched matrices"]["contiguous"]["random"] = BenchmarkGroup()

SUITE["batched matrices"]["in-contiguous"] = BenchmarkGroup()
SUITE["batched matrices"]["in-contiguous"]["ordered"] = BenchmarkGroup()
SUITE["batched matrices"]["in-contiguous"]["random"] = BenchmarkGroup()

SUITE["batched matrices"]["single qubit"] = BenchmarkGroup()
SUITE["batched matrices"]["single qubit"]["ordered"] = BenchmarkGroup()
SUITE["batched matrices"]["single qubit"]["random"] = BenchmarkGroup()

## General Matrix Instruction
function matrices(::Type{T}, N) where {T}
    list = Any[
        rand_unitary(T, N), # dense matrices
        # SparseMatrixCSC
        sprand_hermitian(T, N, 0.1),
        # PermMatrix
        pmrand(T, N),
        Diagonal(rand(T, N)),
    ]
    if N < 100
        # StaticArrays
        push!(list, @SArray(rand(T, N, N)))
        push!(list, @MArray(rand(T, N, N)))
    end
    return list
end

# default test type is ComplexF64
matrices(N) = matrices(ComplexF64, N)

@info "generating benchmark for batched contiguous matrices locs"
### contiguous
for n in 1:2:5, T in [ComplexF64], U in matrices(T, 1 << n), B in 10:20:60
    # contiguous ordered address
    SUITE["batched matrices"]["contiguous"]["ordered"][n, B, string(T), string(typeof(U))] = bench(
        n,
        B,
        U,
        Tuple(1:n),
    )
    # contiguous random address
    SUITE["batched matrices"]["contiguous"]["random"][n, B, string(T), string(typeof(U))] = bench(
        n,
        B,
        U,
        Tuple(randperm(n)),
    )
end

@info "generating benchmark for in-contiguous matrices locs"
### in-contiguous
for m in 1:3, T in [ComplexF64], U in matrices(T, 1 << m), B in 10:20:60
    n = 10
    N = 1 << n
    # in-contiguous ordered address
    SUITE["batched matrices"]["in-contiguous"]["ordered"][m, B, string(T), string(typeof(U))] = bench(
        n,
        B,
        U,
        Tuple(sort(randperm(n)[1:m])),
    )
    # in-contiguous random address
    SUITE["batched matrices"]["in-contiguous"]["random"][m, B, string(T), string(typeof(U))] = bench(
        n,
        B,
        U,
        Tuple(randperm(n)[1:m]),
    )
end

@info "generating benchmark for single qubit matrices"
### single qubit
for T in [ComplexF64], U in matrices(T, 2), n in 1:4:25, B in 10:20:60
    SUITE["batched matrices"]["single qubit"]["ordered"][string(T), B, string(typeof(U)), n] = bench(
        n,
        B,
        U,
        (rand(1:n),),
    )
    SUITE["batched matrices"]["single qubit"]["random"][string(T), B, string(typeof(U)), n] = bench(
        n,
        B,
        U,
        (rand(1:n),),
    )
end

SUITE["batched matrices"]["controlled"] = BenchmarkGroup()
SUITE["batched matrices"]["controlled"]["ordered"] = BenchmarkGroup()
SUITE["batched matrices"]["controlled"]["random"] = BenchmarkGroup()

test_bench(n, U, loc, control_locs, control_bits) =
    instruct!(statevec(rand_state(n)), U, loc, control_locs, control_bits)
for T in [ComplexF64], n in 1:2:5, m in 3:2:6, U in matrices(T, 1 << n), B in 10:20:60
    SUITE["batched matrices"]["controlled"]["ordered"][n, B, m, string(T), string(typeof(U))] = bench(
        n + m,
        B,
        U,
        Tuple(1:n),
        Tuple(n+1:n+m),
        ntuple(x -> 1, m),
    )

    perms = randperm(n + m)
    SUITE["batched matrices"]["controlled"]["random"][n, B, m, string(T), string(typeof(U))] = bench(
        n + m,
        B,
        U,
        Tuple(perms[1:n]),
        Tuple(perms[n+1:n+m]),
        ntuple(x -> 1, m),
    )
end


text:
@testset "Sparse logistic small ($T)" for T in [Float32, Float64]
    using ProximalOperators
    using ProximalAlgorithms
    using LinearAlgebra

    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]

    m, n = size(A)

    R = real(T)

    f = Translate(LogisticLoss(ones(R, m), R(1)), -b)
    lam = R(0.1)
    g = NormL1(lam)

    x_star = T[0, 0, 2.114635341704963e-01, 0, 2.845881348733116e+00]

    TOL = R(1e-6)

    # Nonfast/Adaptive

    x0 = zeros(T, n)
    x0_backup = copy(x0)
    solver = ProximalAlgorithms.ForwardBackward(tol = TOL, adaptive = true, fast = false)
    x, it = solver(x0, f = f, A = A, g = g)
    @test eltype(x) == T
    @test norm(x - x_star, Inf) <= 1e-4
    @test it < 1100
    @test x0 == x0_backup

    # Fast/Adaptive

    x0 = zeros(T, n)
    x0_backup = copy(x0)
    solver = ProximalAlgorithms.ForwardBackward(tol = TOL, adaptive = true, fast = true)
    x, it = solver(x0, f = f, A = A, g = g)
    @test eltype(x) == T
    @test norm(x - x_star, Inf) <= 1e-4
    @test it < 500
    @test x0 == x0_backup

    # ZeroFPR/Adaptive

    x0 = zeros(T, n)
    x0_backup = copy(x0)
    solver = ProximalAlgorithms.ZeroFPR(adaptive = true, tol = TOL)
    x, it = solver(x0, f = f, A = A, g = g)
    @test eltype(x) == T
    @test norm(x - x_star, Inf) <= 1e-4
    @test it < 25
    @test x0 == x0_backup

    # PANOC/Adaptive

    x0 = zeros(T, n)
    x0_backup = copy(x0)
    solver = ProximalAlgorithms.PANOC(adaptive = true, tol = TOL)
    x, it = solver(x0, f = f, A = A, g = g)
    @test eltype(x) == T
    @test norm(x - x_star, Inf) <= 1e-4
    @test it < 50
    @test x0 == x0_backup

end


text:
"""
`process_streaks`
"""

function process_streaks(df::DataFrames.AbstractDataFrame)
    streak_table = combine(groupby(df, :Trial)) do dd
        dt = DataFrame(
        Num_pokes = size(dd,1),
        Num_Rewards = length(findall(skipmissing(dd[!,:Reward]))),
        Start_Reward = dd[1,:Reward],
        Last_Reward = findlast(dd[!,:Reward]).== nothing ? 0 : findlast(dd[!,:Reward]),
        Prev_Reward = findlast(dd[!,:Reward]).== nothing ? 0 : findprev(dd[!,:Reward], findlast(dd[!,:Reward])-1),
        Poking_duration = (dd[end,:PokeOut]-dd[1,:PokeIn]),
        Trial_duration = (dd[end,:ROI_Out]-dd[1,:ROI_In]),
        Start_poking = (dd[1,:PokeIn]),
        Stop_poking = (dd[end,:PokeOut]),
        Start_trial = (dd[1,:ROI_In]),
        Stop_trial = (dd[end,:ROI_Out]),
        Pre_Interpoke = any(map(!,ismissing.(dd[!,:Pre_Interpoke]))) ? maximum(skipmissing(dd[!,:Pre_Interpoke])) : missing,
        Post_Interpoke = any(map(!,ismissing.(dd[!,:Post_Interpoke]))) ? maximum(skipmissing(dd[!,:Post_Interpoke])) : missing,
        #PokeSequence = [SVector{nrow(dd),Union{Bool,Missing}}(dd[!,:Reward])],
        Protocol = dd[1,:Protocol],
        Side = dd[1,:Side],
        ReverseTrial = dd[1,:ReverseTrial],
        Stim = dd[1,:Stim],
        Stim_Day = dd[1,:Stim_Day],
        AverageRewRate = mean(dd.AverageRewRate),
        Leaving_NextPrew = dd[end,:Next_Prew]
        )
        return dt
    end
    streak_table[!,:Prev_Reward] = [x .== nothing ? 0 : x for x in streak_table[:,:Prev_Reward]]
    streak_table[!,:AfterLast] = streak_table[!,:Num_pokes] .- streak_table[!,:Last_Reward];
    streak_table[!,:BeforeLast] = streak_table[!,:Last_Reward] .- streak_table[!,:Prev_Reward].-1;
    prov = lead(streak_table[!,:Start_poking],default = 0.0) .- streak_table[!,:Stop_poking];
    streak_table[!,:Poking_Travel_to] = [x.< 0 ? 0 : x for x in prov]
    prov = lead(streak_table[!,:Start_trial],default = 0.0) .- streak_table[!,:Stop_trial];
    streak_table[!,:Trial_Travel_to] = [x.< 0 ? 0 : x for x in prov]
    # streak_table[!,:Travel_to] = streak_table.Poking_Travel_to + streak_table.Trial_Travel_to
    streak_table[!,:ROI_Leaving_Time] = streak_table.Stop_trial - streak_table.Stop_poking
    return streak_table
end


text:
# A helper function to let our `method` keyword correctly understands `Degree` inputs.
@inline wrap_BSpline(itp::Interpolations.InterpolationType) = itp
@inline wrap_BSpline(degree::Interpolations.Degree) = BSpline(degree)

# The default values used by extrapolation for off-domain points
const FillType = Union{Number,Colorant,Flat,Periodic,Reflect}
const FloatLike{T<:AbstractFloat} = Union{T,AbstractGray{T}}
const FloatColorant{T<:AbstractFloat} = Colorant{T}
@inline _default_fillvalue(::Type{T}) where {T<:FloatLike} = convert(T, NaN)
@inline _default_fillvalue(::Type{T}) where {T<:FloatColorant} = _nan(T)
@inline _default_fillvalue(::Type{T}) where {T} = zero(T)

@inline _make_compatible(T, fill) = fill
@inline _make_compatible(::Type{T}, fill::Number) where {T} = T(fill)

Interpolations.tweight(A::AbstractArray{C}) where C<:Colorant{T} where T = T

const MethodType = Union{Interpolations.Degree, Interpolations.InterpolationType}

function box_extrapolation(
        parent::AbstractArray;
        fillvalue::FillType = _default_fillvalue(eltype(parent)),
        method=Linear(),
        kwargs...)
    T = typeof(zero(Interpolations.tweight(parent)) * zero(eltype(parent)))
    itp = maybe_skip_prefilter(T, parent, method)
    extrapolate(itp, _make_compatible(T, fillvalue))
end
box_extrapolation(etp::AbstractExtrapolation) = etp
box_extrapolation(itp::AbstractInterpolation{T}; fillvalue=_default_fillvalue(T)) where T =
    extrapolate(itp, _make_compatible(T, fillvalue))

@inline function maybe_skip_prefilter(::Type{T}, A::AbstractArray, degree::D) where {T, D<:Union{Linear, Constant}}
    axs = axes(A)
    return Interpolations.BSplineInterpolation{T,ndims(A),typeof(A),BSpline{D},typeof(axs)}(A, axs, BSpline(degree))
end
@inline function maybe_skip_prefilter(::Type{T}, A::AbstractArray, method::BSpline{D}) where {T, D<:Union{Linear, Constant}}
    maybe_skip_prefilter(T, A, construct_interpolation_type(D))
end

@inline function maybe_skip_prefilter(::Type{T}, A::AbstractArray, method::MethodType) where T
    return interpolate(A, wrap_BSpline(method))
end

# This is type-piracy, but necessary if we want Interpolations to be
# independent of OffsetArrays.
function AxisAlgorithms.A_ldiv_B_md!(dest::OffsetArray, F, src::OffsetArray, dim::Integer, b::AbstractVector)
    indsdim = axes(parent(src), dim)
    indsF = axes(F)[2]
    if indsF == indsdim
        AxisAlgorithms.A_ldiv_B_md!(parent(dest), F, parent(src), dim, b)
        return dest
    end
    throw(DimensionMismatch("indices $(axes(parent(src))) do not match $(axes(F))"))
end


text:
const DIMS = (:, 0x2, 3, 4)
const ITRS = (
    DIMS...,
    1:1,
    1:2,
    1:3,
    [true, true, true],
    not(1),
    FI(iseven),
    [false, true, true],
    [true, false, true],
    [true, true, false],
    [false, false, true],
    [false, true, false],
    [true, false, false],
    # [false, false, false],
)

struct WarpedArray{T,N,P} <: ResizingTools.AbstractRDArray{T,N}
    parent::P
    WarpedArray(A::P) where {T,N,P<:AbstractArray{T,N}} = new{T,N,P}(A)
end
warpsimplerdarray(A::AbstractArray) = WarpedArray(SimpleRDArray(A))

Base.parent(A::WarpedArray) = A.parent
ArrayInterface.parent_type(::Type{<:WarpedArray{T,N,P}}) where {T,N,P} = P
Base.IndexStyle(::Type{A}) where {A<:WarpedArray} =
    IndexStyle(ArrayInterface.parent_type(A))
Base.@propagate_inbounds Base.getindex(A::WarpedArray, Is...) = parent(A)[Is...]
Base.@propagate_inbounds Base.setindex!(A::WarpedArray, v, Is...) = parent(A)[Is...] = v
ResizingTools.getsize(A::WarpedArray) = ResizingTools.getsize(parent(A))
ResizingTools.getsize(A::WarpedArray, d::Int) = ResizingTools.getsize(parent(A), d)

_geninds(A::AbstractArray, d::Integer, I) =
    ntuple(i -> i == d ? _to_inds(A, d, I) : axes(A, i), Val(ndims(A)))
_geninds(A::AbstractArray, Is::Tuple) = ntuple(i -> _to_inds(A, i, Is[i]), Val(ndims(A)))

_to_inds(A::AbstractArray, d::Integer, itr) = ResizingTools.to_index(A, axes(A, d), itr)
_to_inds(A::AbstractArray, d::Integer, ::Colon) = axes(A, d)
_to_inds(A::AbstractArray, d::Integer, n::Integer) = Base.OneTo(min(n, size(A, d)))

_to_sinds(inds) = map(ind -> Base.OneTo(length(ind)), inds)

_to_size(A::AbstractArray, Is::Tuple) = ntuple(i -> _to_len(A, i, Is[i]), Val(ndims(A)))
_to_len(::AbstractArray, ::Integer, itr) = eltype(itr) <: Bool ? sum(itr) : length(itr)
_to_len(A::AbstractArray, d::Integer, I::FunctionIndices.AbstractFunctionIndex) =
    length(_to_inds(A, d, I))
_to_len(::AbstractArray, ::Integer, n::Integer) = n
_to_len(A::AbstractArray, d::Integer, ::Colon) = size(A, d)

struct NLoop{N,T} <: AbstractArray{Any,N}
    itr::T
end
NLoop{N}(tp::Tuple) where {N} = NLoop{N,typeof(tp)}(tp)
Base.IndexStyle(::Type{<:NLoop}) = IndexCartesian()
Base.size(A::NLoop{N}) where {N} = ntuple(_ -> length(A.itr), Val(N))
Base.getindex(A::NLoop{N}, is::Vararg{Int,N}) where {N} = ntuple(i -> A.itr[is[i]], Val(N))

function test_resize(f, g, A::AbstractArray{T,N}, itrs, dnums=1:N; pre=identity) where {T,N}
    @testset "resize!($(typeof(A)), $_Is)" for _Is in NLoop{N}(itrs)
        Is = pre(_Is)
        tA = g(f(A))
        fA = f(A)
        resize!(tA, Is)
        inds = _geninds(fA, Is)
        @test tA[_to_sinds(inds)...] == fA[inds...]
        @test size(tA) == _to_size(fA, Is)
    end
    @testset "resize!($(typeof(A)), $d, $I)" for d in dnums, I in itrs
        tA = g(f(A))
        fA = f(A)
        resize!(tA, d, I)
        inds = _geninds(fA, d, I)
        @test tA[_to_sinds(inds)...] == fA[inds...]
        @test size(tA, d) == _to_len(fA, d, I)
    end
end
test_resize(g, A::AbstractArray{T,N}, itrs, dnums=1:N; pre=identity) where {T,N} =
    test_resize(identity, g, A, itrs, dnums; pre=pre)

@testset "resize!: $f" for f in (SimpleRDArray, warpsimplerdarray)
    for A in (V, M, T)
        test_resize(f, A, ITRS)
    end
end

adjIS(I) = length(I) == 1 ? (1, I[1]) : I
@testset "resize!: adjoint" begin
    test_resize(adjoint, SimpleRDArray, V, ITRS, 2; pre=adjIS)
    test_resize(adjoint, SimpleRDArray, M, ITRS)
end

@testset "sizehint!: $f" for f in (SimpleRDArray, warpsimplerdarray)
    @testset "Vector" begin
        @test begin
            tV = f(V)
            @allocated resize!(tV, (4,))
        end > begin
            tV = f(V)
            sizehint!(tV, (4,))
            @allocated resize!(tV, (4,))
        end
        @test begin
            tV = f(V)
            sizehint!(parent(tV), 4)
            @allocated resize!(tV, (4,))
        end == begin
            tV = f(V)
            sizehint!(tV, (4,))
            @allocated resize!(tV, (4,))
        end
    end
    @testset "Matrix" begin
        @test begin
            tM = f(M)
            @allocated resize!(tM, (4, 4))
        end > begin
            tM = f(M)
            sizehint!(tM, (4, 4))
            @allocated resize!(tM, (4, 4))
        end
        @test begin
            tM = f(M)
            sizehint!(parent(tM), 16)
            @allocated resize!(tM, (4, 4))
        end == begin
            tM = f(M)
            sizehint!(tM, (4, 4))
            @allocated resize!(tM, (4, 4))
        end
    end
    @testset "3-rank array" begin
        @test begin
            tT = f(T)
            @allocated resize!(tT, (4, 4, 4))
        end > begin
            tT = f(T)
            sizehint!(tT, (4, 4, 4))
            @allocated resize!(tT, (4, 4, 4))
        end
        @test begin
            tT = f(T)
            sizehint!(parent(tT), 64)
            @allocated resize!(tT, (4, 4, 4))
        end == begin
            tT = f(T)
            sizehint!(tT, (4, 4, 4))
            @allocated resize!(tT, (4, 4, 4))
        end
    end
end

@testset "resize! with both Integer and AbstractArray #4" begin
    tM = SimpleRDArray(reshape(1:4, 2, 2))
    @test resize!(tM, (2, 1:2)) == reshape(1:4, 2, 2)
    @test resize!(tM, (1:2, 1)) == reshape(1:2, 2, 1)
end

@testset "resize!(::Vector, ::Any)" begin
    @test resize!(collect(1:4), (3,)) == 1:3
    @test resize!(collect(1:4), (5,))[1:4] == 1:4
    @test resize!(collect(1:4), (:,)) == 1:4
    @test resize!(collect(1:4), (1:3,)) == 1:3
    @test resize!(collect(1:4), (Bool[1, 0, 1, 1],)) == [1, 3, 4]
end


text:
# Unit 02 - Nonlinear Classification, Linear regression, Collaborative Filtering

# Lecture 5

# Hinge loss computation

x = [1 0 1; 1 1 1; 1 1 -1; -1 1 1]
y = [2,2.7,-0.7, 2]
θ = [0,1,2]
z = y - (θ' * x')'

hinge(z) = (z >= (1 - eps()) ) ? 0 : (1 - z)
l = mean(map(z -> hinge(z), z))

squared_error(z) = z^2 / 2
l = mean(map(z -> squared_error(z), z))

# Lecture 6

# radial basis kernel
using LinearAlgebra
radial_kernel(x,xᵖ) = exp(-1/2 * norm(x-xᵖ)^2)

x = [1,0,0]

xᵖ = [0,1,0]

k = radial_kernel(x,xᵖ)

max_stars_repo_path: Unit 02 - Nonlinear Classification, Linear regression, Collaborative Filtering/Unit 02 - Nonlinear Classification, Linear regression, Collaborative Filtering.jl


text:
using Random
using Distributions
using MultiAgentPathFinding
using MultiAgentAllocationTransit
using BenchmarkTools
using Statistics
using JSON

rng = MersenneTwister(2345)

const params_file = ARGS[1]
city_params = parse_city_params(params_file)

const out_file_pref = ARGS[2]

lat_dist = Uniform(city_params.lat_start, city_params.lat_end)
lon_dist = Uniform(city_params.lon_start, city_params.lon_end)

# N_agents = N_depots
const N_DEPOTS = parse(Int64, ARGS[3])
const depots = [LatLonCoords((lat = rand(rng, lat_dist), lon = rand(rng, lon_dist))) for i = 1:N_DEPOTS]
const N_AGENTS = N_DEPOTS

# const N_SITE_VALS = [50, 100, 200, 500, 1000, 5000]
const N_SITE_VALS = [5000]

for N_SITES in N_SITE_VALS

    sites = [LatLonCoords((lat = rand(rng, lat_dist), lon = rand(rng, lon_dist))) for i = 1:N_SITES]

    depot_sites = vcat(depots, sites)
    cost_fn(i, j) = MultiAgentAllocationTransit.distance_lat_lon_euclidean(depot_sites[i], depot_sites[j])

    b = @benchmarkable task_allocation($N_DEPOTS, $N_SITES, $N_AGENTS, $depot_sites, $cost_fn)

    t = run(b, seconds=100000, samples=5) # Do once to trigger...something
    # t = run(b)

    @show N_SITES
    @show mean(t.times)
    @show dump(t)

    allocation_results = Dict("n_depots"=>N_DEPOTS, "n_agents"=>N_AGENTS, "site_stats"=>Dict())
    allocation_results["site_stats"][N_SITES] = Dict("mean" => mean(t.times),
                                                     "median" => median(t.times),
                                                     "std" => std(t.times),
                                                     "samples" => t.params.samples)

    out_file = string(out_file_pref,"_",N_SITES,"sites.json")
    open(out_file, "w") do f
        JSON.print(f, allocation_results, 2)
    end
end


text:
using Distributed
using Logging
using Printf
using Dates

function reduce_results(res1, results...)
    new_res = typeof(res1)()
    for k in keys(res1)
        new_res[k] = vcat(res1[k], (r[k] for r in results)...)
    end
    new_res
end

const MAX_BATCH_SIZE = 1

function run_parallel(systemfn, algorithm, num_samples; compile_args = (;))
    batches = Int[]

    N = num_samples
    batch_size = clamp(floor(num_samples / nworkers()), 1, MAX_BATCH_SIZE)

    while num_samples > 0
        batch = min(num_samples, batch_size)
        push!(batches, batch)
        num_samples -= batch
    end

    @everywhere begin
        system = $systemfn()
        global compiled_system = PathWeightSampling.compile(system; $compile_args...)
    end

    channel = RemoteChannel(() -> Channel())
    @info "Starting Parallel computation"
    result = @sync begin
        # display progress
        progress = 0
        @async while true
            val = take!(channel)
            if val === nothing
                break
            end
            (hostname, elapsed_time, batch_size) = val
            progress += batch_size
            done = @sprintf "%i/%i" progress N
            percent_done = @sprintf "%6.2f %%" (progress / N * 100)
            time = now()
            @info "Finished batch" hostname time elapsed_time batch_size done percent_done
            flush(stderr)
        end

        @sync begin
            result = pmap(batches) do batch
                time_stats = @timed result = _mi_inner(Main.compiled_system, algorithm, batch, false)
                elapsed_time = time_stats.time
                hostname = gethostname()
                put!(channel, (hostname, elapsed_time, batch))
                yield()
                result
            end
            put!(channel, nothing)
            vcat(result...)
        end
    end
    @info "Finished Parallel computation"
    result
end


text:
###############
# H-OBCA: Hierarchical Optimization-based Collision Avoidance - a path planner for autonomous parking
# Copyright (C) 2018
# Alexander LINIGER [liniger@control.ee.ethz.ch; Automatic Control Lab, ETH Zurich]
# Xiaojing ZHANG [xiaojing.zhang@berkeley.edu; MPC Lab, UC Berkeley]
# Atsushi SAKAI [atsushisakai@global.komatsu; Komatsu Ltd / MPC Lab]
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
###############
# The paper describing the theory can be found here:
# 	X. Zhang, A. Liniger and F. Borrelli; "Optimization-Based Collision Avoidance"; Technical Report, 2017, [https://arxiv.org/abs/1711.03449]
#   X. Zhang, A. Liniger, A. Sakai and F. Borrelli; "Autonomous  Parking  using  Optimization-Based  Collision  Avoidance"; Technical Report, 2018 [add URL]
###############

###############
# computes collision-free trajectory by appropriately reformulating the distance function
###############

function ParkingSignedDist(x0,xF,N,Ts,L,ego,XYbounds,nOb,vOb, A, b,fixTime,xWS,uWS)

	# desired safety distance
	dmin = 0.05			# anything bigger than 0, e.g. 0.05

	##############################
	# Define JuMP file
	##############################
	# Define IPOPT as solver and well as solver settings
	##############################
	# seems to work best
	m = Model(solver=IpoptSolver(hessian_approximation="exact",mumps_pivtol=1e-6,alpha_for_y="min",recalc_y="yes",
	                             mumps_mem_percent=6000,max_iter=200,tol=1e-5, print_level=0,
	                             min_hessian_perturbation=1e-12,jacobian_regularization_value=1e-7))#,nlp_scaling_method="none"

	##############################
	# defining optimization variables
	##############################
	#state
	@variable(m, x[1:4,1:(N+1)])
	#scaling on sampling time
	if fixTime == 0
		@variable(m, timeScale[1:N+1])
	end
	# timeScale = ones(1,N+1)
	#control
	@variable(m, u[1:2,1:(N)])
	# lagrange multipliers for dual dist function
	@variable(m, l[1:sum(vOb),1:(N+1)])	# dual multiplier associated with obstacleShape
	@variable(m, n[1:nOb*4,1:(N+1)])	# dual multiplier associated with carShape

	# regularization parameter to improve numerical stability
	reg = 1e-7;
	##############################
	# cost function
	##############################
	# (min control inputs)+
	# (min input rate)
	# (min time)+
	# (regularization wrt HA* traj)
	##############################

	u0 = [0,0]
	#fix time objective
	if fixTime == 1
		@NLobjective(m, Min,sum(0.01*u[1,i]^2 + 0.5*u[2,i]^2 for i = 1:N) +
		 					 sum(0.1*((u[1,i+1]-u[1,i])/Ts)^2 + 0.1*((u[2,i+1]-u[2,i])/Ts)^2 for i = 1:N-1)+
		 					 	(0.1*((u[1,1]-u0[1])/(Ts))^2 + 0.1*((u[2,1]-u0[2])/(Ts))^2) +
		 					 sum(0.001*(x[1,i]-xWS[i,1])^2 + 0.001*(x[2,i]-xWS[i,2])^2 + 0.0001*(x[3,i]-xWS[i,3])^2 for i=1:N+1))
	else
	#varo time objective
		@NLobjective(m, Min,sum(0.01*u[1,i]^2 + 0.1*u[2,i]^2 for i = 1:N) +
		 					 sum(0.1*((u[1,i+1]-u[1,i])/(timeScale[i]*Ts))^2 + 0.1*((u[2,i+1]-u[2,i])/(timeScale[i]*Ts))^2 for i = 1:N-1) +
		 					    (0.1*((u[1,1]-u0[1])   /(timeScale[1]*Ts))^2 + 0.1*((u[2,1]-u0[2])   /(timeScale[1]*Ts))^2) +
		 					 sum(0.5*timeScale[i] + 1*timeScale[i]^2 for i = 1:N+1)+
		 					 sum(0.001*(x[1,i]-xWS[i,1])^2 + 0.001*(x[2,i]-xWS[i,2])^2 + 0.0001*(x[3,i]-xWS[i,3])^2 for i=1:N+1))
	end

	##############################
	# bounds on states, inputs,
	# and dual multipliers.
	##############################
	#input constraints
	@constraint(m, [i=1:N], -0.6 <= u[1,i] <= 0.6)
	@constraint(m, [i=1:N], -0.4 <= u[2,i] <= 0.4)

	#state constraints
	@constraint(m, [i=1:N+1], XYbounds[1] <= x[1,i] <= XYbounds[2])
	@constraint(m, [i=1:N+1], XYbounds[3] <= x[2,i] <= XYbounds[4])

	@constraint(m, [i=1:N+1], -1  <= x[4,i] <= 2)

	# bounds on time scaling
	if fixTime == 0
		@constraint(m, 0.8 .<= timeScale .<= 1.2)
	end

	# positivity constraints on dual multipliers
	@constraint(m, l .>= 0)
	@constraint(m, n .>= 0)

	##############################
	# start and finish point
	##############################

	#starting point
	@constraint(m, x[1,1] == x0[1])
	@constraint(m, x[2,1] == x0[2])
	@constraint(m, x[3,1] == x0[3])
	@constraint(m, x[4,1] == x0[4])

	#end point
	@constraint(m, x[1,N+1] == xF[1])
	@constraint(m, x[2,N+1] == xF[2])
	@constraint(m, x[3,N+1] == xF[3])
	@constraint(m, x[4,N+1] == xF[4])

	##############################
	# dynamics of the car
	##############################
	# - unicycle dynamic with euler forward
	# - sampling time scaling, is identical over the horizon

	for i in 1:N

 		if fixTime == 1  # sampling time is fixed
			@NLconstraint(m, x[1,i+1] == x[1,i] + Ts*(x[4,i] + Ts/2*u[2,i])*cos((x[3,i] + Ts/2*x[4,i]*tan(u[1,i])/L)))
		    @NLconstraint(m, x[2,i+1] == x[2,i] + Ts*(x[4,i] + Ts/2*u[2,i])*sin((x[3,i] + Ts/2*x[4,i]*tan(u[1,i])/L)))
		    @NLconstraint(m, x[3,i+1] == x[3,i] + Ts*(x[4,i] + Ts/2*u[2,i])*tan(u[1,i])/L)
		    @NLconstraint(m, x[4,i+1] == x[4,i] + Ts*u[2,i])
	    else			# sampling time is variable
		    @NLconstraint(m, x[1,i+1] == x[1,i] + timeScale[i]*Ts*(x[4,i] + timeScale[i]*Ts/2*u[2,i])*cos((x[3,i] + timeScale[i]*Ts/2*x[4,i]*tan(u[1,i])/L)))
		    @NLconstraint(m, x[2,i+1] == x[2,i] + timeScale[i]*Ts*(x[4,i] + timeScale[i]*Ts/2*u[2,i])*sin((x[3,i] + timeScale[i]*Ts/2*x[4,i]*tan(u[1,i])/L)))
		    @NLconstraint(m, x[3,i+1] == x[3,i] + timeScale[i]*Ts*(x[4,i] + timeScale[i]*Ts/2*u[2,i])*tan(u[1,i])/L)
		    @NLconstraint(m, x[4,i+1] == x[4,i] + timeScale[i]*Ts*u[2,i])
	    end

	    if fixTime == 0
	    	@constraint(m, timeScale[i] == timeScale[i+1])
    	end

	end

	u0 = [0,0]
	if fixTime == 1
		for i in 1:N
			if i==1
				@constraint(m,-0.6<=(u0[1]-u[1,i])/Ts <= 0.6)
			else
				@constraint(m,-0.6<=(u[1,i-1]-u[1,i])/Ts <= 0.6)
			end
		end
	else
		for i in 1:N
			if i==1
				@NLconstraint(m,-0.6<=(u0[1]-u[1,i])/(timeScale[i]*Ts) <= 0.6)
			else
				@NLconstraint(m,-0.6<=(u[1,i-1]-u[1,i])/(timeScale[i]*Ts) <= 0.6)
			end
		end
	end

	##############################
	# obstacle avoidance constraints
	##############################

	# width and length of ego set
	W_ev = ego[2]+ego[4]
	L_ev = ego[1]+ego[3]

	g = [L_ev/2,W_ev/2,L_ev/2,W_ev/2]

	# ofset from X-Y to the center of the ego set
	offset = (ego[1]+ego[3])/2 - ego[3]

	for i in 1:N+1 	# iterate over time steps
		for j = 1 : nOb 	# iterate over obstacles
			Aj = A[sum(vOb[1:j-1])+1 : sum(vOb[1:j]) ,:]	# extract obstacle matrix associated with j-th obstacle
			lj = l[sum(vOb[1:j-1])+1 : sum(vOb[1:j]) ,:]	# extract lambda dual variables associate j-th obstacle
			nj = n[(j-1)*4+1:j*4 ,:] 						# extract mu dual variables associated with j-th obstacle
			bj = b[sum(vOb[1:j-1])+1 : sum(vOb[1:j])]		# extract obstacle matrix associated with j-th obstacle

			# norm(A'*lambda) = 1
			@NLconstraint(m, (sum(Aj[k,1]*lj[k,i] for k = 1 : vOb[j]))^2 + (sum(Aj[k,2]*lj[k,i] for k = 1 : vOb[j]))^2  == 1   )

			# G'*mu + R'*A*lambda = 0
			@NLconstraint(m, (nj[1,i] - nj[3,i]) +  cos(x[3,i])*sum(Aj[k,1]*lj[k,i] for k = 1:vOb[j]) + sin(x[3,i])*sum(Aj[k,2]lj[k,i] for k = 1:vOb[j]) == 0  )
			@NLconstraint(m, (nj[2,i] - nj[4,i]) -  sin(x[3,i])*sum(Aj[k,1]*lj[k,i] for k = 1:vOb[j]) + cos(x[3,i])*sum(Aj[k,2]lj[k,i] for k = 1:vOb[j]) == 0  )

			# -g'*mu + (A*t - b)*lambda > 0
			@NLconstraint(m, (-sum(g[k]*nj[k,i] for k = 1:4) + (x[1,i]+cos(x[3,i])*offset)*sum(Aj[k,1]*lj[k,i] for k = 1:vOb[j])
							+ (x[2,i]+sin(x[3,i])*offset)*sum(Aj[k,2]*lj[k,i] for k=1:vOb[j]) - sum(bj[k]*lj[k,i] for k=1:vOb[j]))  >= dmin  )
		end
	end

	##############################
	# set initial guesses
	##############################

	if fixTime == 0
		setvalue(timeScale,1*ones(N+1,1))
	end
	setvalue(x,xWS')
	setvalue(u,uWS[1:N,:]')

	lWS,nWS = DualMultWS(N,nOb,vOb, A, b,xWS[:,1],xWS[:,2],xWS[:,3])

	setvalue(l,lWS')
	setvalue(n,nWS')

	##############################
	# solve problem
	##############################
	# ipopt has sometimes problems in the restoration phase,
	# it turns out that restarting ipopt with the previous solution
	# as an initial guess works well to achieve a high success rate.
	#
	# if restoration failure is reported by IPOPT, solution should be checked manually as it can still be feasible
	##############################

	# at most three attempts considered
	time1 = 0
	time2 = 0

	exitflag = 0

	tic()
	status = solve(m; suppress_warnings=true)
	time1 = toq();

	if status == :Optimal
	    exitflag = 1
	elseif status ==:Error || status ==:UserLimit

		Feasible = 0
		if Feasible == 0
		    tic()
		    status = solve(m; suppress_warnings=true)
		    time2 = toq();

		    if status == :Optimal
		        exitflag = 1
		    elseif status ==:Error || status ==:UserLimit
		    	xp = getvalue(x)
				up = getvalue(u)
				if fixTime == 1
					timeScalep = ones(1,N+1)
				else
					timeScalep = getvalue(timeScale)
				end
				lp = getvalue(l)
				np = getvalue(n)
				Feasible = 0
				Feasible = ParkingConstraints(x0,xF,N,Ts,L,ego,XYbounds,nOb,vOb, A, b,xp,up,lp,np,timeScalep,fixTime,1)

		        if Feasible == 1
		            exitflag = 1
		        else
		            exitflag = 0
		        end
	        end
	    else
	        exitflag = 1
	    end
	else
	    exitflag = 0
	end

	##############################
	# return values
	##############################

	# computation times is the sum of all trials
	time = time1+time2
	# print("  elapsed time: ")
	# print(time)
	# println(" seconds")

	xp = getvalue(x)
	up = getvalue(u)
	if fixTime == 1
		timeScalep = ones(1,N+1)
	else
		timeScalep = getvalue(timeScale)
	end

	lp = getvalue(l)
	np = getvalue(n)

	return xp, up, timeScalep, exitflag, time, lp, np

end

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
function foo()
    a = 10
    b = a
    @assert(b == 10)
    println(join([b], " "))
end

function main()
    foo()
end

main()


text:
using OptiMo, Bazinga
using Printf

include("nonconvex_qp_tiny.jl")

prob = PPNCVXQP()

solver = Bazinga.PANOC( verbose=true, freq=1 )
@time out = solver( prob )
print( out )

@printf "\n\n"
solver = Bazinga.ZEROFPR( verbose=true, freq=1 )
@time out = solver( prob )
print( out )


text:
module LinearInterpolations

@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end LinearInterpolations

export interpolate, Interpolate

using ArgCheck
import Adapt
import StaticArrays

struct TinyVector{T} <: AbstractVector{T}
    elements::NTuple{2,T}
end
TinyVector(x, y) = TinyVector(promote(x,y))
Base.size(o::TinyVector) = (2,)
Base.@propagate_inbounds function getindex(o::TinyVector, i::Integer)
    index::Int = Int(i)
    return o[index]
end
@inline function Base.getindex(o::TinyVector, i::Int)
    @boundscheck checkbounds(o, i)
    @inbounds ifelse(i === 1, o.elements[1], o.elements[2])
end

struct WeightsArray{T,N,Factors} <: AbstractArray{T,N}
    factors::Factors
end
@inline function WeightsArray(factors::NTuple{N, TinyVector}) where {N}
    T = promote_type(map(eltype, factors)...)
    Factors = typeof(factors)
    return WeightsArray{T,N,Factors}(factors)
end
function Base.size(o::WeightsArray{T,N}) where {T,N}
    return ntuple(_->2, Val(N))
end
function Base.getindex(o::WeightsArray{T,N}, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(o, I...)
    xs = map(o.factors, I) do v, i
        @inbounds v[i]
    end
    prod(xs)::T
end

struct NeighborsArray{T,N,Objs, Inds} <: AbstractArray{T,N}
    objs::Objs
    inds::Inds
end
function Base.size(o::NeighborsArray{T,N}) where {T,N}
    ntuple(_->2, Val(N))
end
function Base.getindex(o::NeighborsArray{T,N}, I::Vararg{Int, N}) where {T,N}
    @boundscheck checkbounds(o, I...)
    inds_big = map(o.inds, I) do inds, i
        @inbounds inds[i]
    end
    @inbounds o.objs[inds_big...]
end
function NeighborsArray(objs::AbstractArray{T,N}, inds::NTuple{N,TinyVector}) where {T,N}
    Objs = typeof(objs)
    Inds = typeof(inds)
    NeighborsArray{T,N,Objs, Inds}(objs,inds)
end

const EXTRAPOLATE_SYMBOLS = [:replicate, :reflect, :error, :fuzzy]

function project1d(extrapolate::Symbol, xs, x)
    if extrapolate === :replicate
        project1d(Replicate(), xs, x)
    elseif extrapolate === :reflect
        project1d(Reflect(), xs, x)
    elseif extrapolate === :error
        project1d(Error(), xs, x)
    elseif extrapolate === :fuzzy
        project1d(Fuzzy(), xs, x)
    else
        @argcheck extrapolate in EXTRAPOLATE_SYMBOLS
        error("""Unreachable
              extrapolate = $extrapolate
              EXTRAPOLATE_SYMBOLS = $EXTRAPOLATE_SYMBOLS
        """)
    end
end

"""
    Fuzzy(;atol, rtol)

Throw an error, for data far away from the grid, but project
it onto the grid, if it is approximate accoriding to the fields
of `Fuzzy`.
"""
struct Fuzzy{K <: NamedTuple}
    kw::K
end
function Fuzzy(;kw...)
    Fuzzy((;kw...))
end

function project1d(fuzzy::Fuzzy, xs, x)
    x_inside = clamp(x, first(xs), last(xs))
    if isapprox(x, x_inside; fuzzy.kw...)
        x_inside
    else
        msg = """
        x=$x is not approximately between first(xs)=$(first(xs)) and last(xs)=$(last(xs))
        Further information:
        x = $x
        x_inside = $x_inside
        extrapolate = $fuzzy
        """
        throw(ArgumentError(msg))
    end
end

"""
    Replicate()

Extrapolate data by replacing out of grid points with the closest in grid point.
"""
struct Replicate end

function project1d(::Replicate, xs, x)
    return clamp(x, first(xs), last(xs))
end

"""
    Reflect()

Extrapolate data by reflecting out of grid points into the grid.
"""
struct Reflect end

function project1d(::Reflect, xs, x)
    x_first = first(xs)
    x_last = last(xs)
    Δx = x_last - x_first
    x1 = if x <= x_first
        n = floor((x_last-x)/(2Δx))
        x + n*2Δx
    elseif x >= x_last
        n = floor((x-x_first)/(2Δx))
        x - n*2Δx
    else
        x
    end
    x2 = if x1 > x_last
        2x_last - x1
    elseif x1 < x_first
        2x_first - x1
    else
        x1
    end
    # It should now hold that:
    # x_first <= x <= x_last
    # but floating point issues might prevent that?
    x3 = clamp(x2, x_first, x_last)
    return x3
end

"""
    Error()

Throw an error when trying to interpolate outside of the grid.
"""
struct Error end

function project1d(::Error, xs, x)
    x_first = first(xs)
    x_last = last(xs)
    if (x_first <= x <= x_last)
        x
    else
        msg = """
        x=$x is not between first(xs)=$(first(xs)) and last(xs)=$(last(xs))
        You can suppress this error by passing the `extrapolate` argument.
        """
        throw(ArgumentError(msg))
    end
end

"""
    AssumeInside()

Assume without checking, that a point is inside the grid when interpolating.
Shoud the point lie outside of the grid, behaviour is undefined.
"""
struct AssumeInside end
project1d(::AssumeInside, xs, x) = x

"""
    Constant(value)

When evaluating at a point outside the grid, return value.
"""
struct Constant{C}
    value::C
end
function apply_extrapolate(o::Constant, pt)
    o.value
end

"""
    WithPoint(f)

When evaluating at a point `pt` outside the grid, return `f(pt)`.
"""
struct WithPoint{F}
    f::F
end
function apply_extrapolate(o::WithPoint, pt)
    o.f(pt)
end

function _neighbors_and_weights1d(xs, x, extrapolate)
    x = project1d(extrapolate, xs, x)
    # searchsortedfirst: index of first value in xs greater than or equal to x
    # since we called clamp, we are inbounds
    ixu = searchsortedfirst(xs, x)
    ixu = ifelse(ixu == firstindex(xs), ixu + 1, ixu)
    ixl = ixu - 1
    xu = @inbounds xs[ixu]
    xl = @inbounds xs[ixl]
    xu_eq_xl = xu == xl
    w_tot = ifelse(xu_eq_xl, one(xu - xl), xu - xl)
    l = ifelse(xu_eq_xl, one(xu - x), xu - x)
    u = ifelse(xu_eq_xl, zero(x - xl), x - xl)
    wl = l / w_tot
    wu = u / w_tot
    wts = TinyVector(wl, wu)
    nbs = TinyVector(ixl, ixu)
    return (nbs, wts)
end

function _neighbors_and_weights(axes::NTuple{N,Any}, pt::NTuple{N,Any}, extrapolate) where {N}
    nbs_wts = let extrapolate=extrapolate
        map(axes, pt) do xs, x
            _neighbors_and_weights1d(xs, x, extrapolate)
        end
    end
    wts = WeightsArray(map(last, nbs_wts))
    nbs = map(first, nbs_wts)
    return nbs, wts
end

"""
    combine(weights::AbstractArray, objects::AbstractArray)::Object

Take a weighted collection of objects and combine them into a single object.
The default behaviour is `sum(weights .* objects)`

This `combine` can be overloaded to allow interpolation of objects that do not implement `*` or `+`.
"""
@inline function combine(wts, objs)
    _combine(wts, objs)
end
function _combine(wts::AbstractVector, objs::AbstractVector)
    w1,w2 = wts
    o1,o2 = objs
    return muladd(w2,o2,w1*o1)
end
function _combine(wts::AbstractMatrix, objs::AbstractMatrix)
    w1,w2,w3,w4 = wts
    o1,o2,o3,o4 = objs
    (w1*o1) + (w2*o2) + (w3*o3) + (w4*o4)
end
function _combine(wts::AbstractArray{Wt,3}, objs::AbstractArray{Obj,3}) where {Wt, Obj}
    w1,w2,w3,w4,w5,w6,w7,w8, = wts
    o1,o2,o3,o4,o5,o6,o7,o8 = objs
    (w1*o1) + (w2*o2) + (w3*o3) + (w4*o4) + (w5*o5) + (w6*o6) + (w7*o7) + (w8*o8)
end

function _combine(wts, objs)
    #mapreduce(*, +, wts, objs) # this allocates urgs
    w1, state_wts = iterate(wts)
    o1, state_objs = iterate(objs)
    ret = w1*o1
    while true
        next_wts = iterate(wts, state_wts)
        next_wts === nothing && break
        wi, state_wts = next_wts
        oi, state_objs = iterate(objs, state_objs)
        ret = muladd(wi,oi,ret)
    end
    return ret
end

"""
    interpolate(axes, objs, pt; extrapolate = :error, [,combine])
    interpolate(xs::AbstractVector, ys::AbstractVector, pt; kw...)

Create an `Interpolate` and evaluate it at `pt`.
See [`Interpolate`](@ref) for details.
"""
function interpolate(axes, objs, pt; extrapolate = :error, combine::C = combine) where {C}
    itp = Interpolate(combine, axes, objs, extrapolate)
    return itp(pt)
end

function interpolate(xs::AbstractVector, ys::AbstractVector, pt; kw...)
    axes = (xs,)
    return interpolate(axes, ys, pt; kw...)
end

function unsafe_Interpolate end
struct Interpolate{C,A,V,O}
    combine::C
    axes::A
    values::V
    extrapolate::O
    function Interpolate(combine, axes::Tuple, values, extrapolate)
        @argcheck size(values) == map(length, axes)
        @argcheck all(size(values) .>= 2)
        @argcheck Base.axes(values) == map(eachindex, axes)
        if extrapolate isa Symbol
            @argcheck extrapolate in EXTRAPOLATE_SYMBOLS
        end
        # some sanity checks, all of them O(ndims(values))
        # we assume issorted(xs), which wouldbe O(length(values)) to check
        @argcheck all((!isempty).(axes))
        @argcheck all(isfinite.(first.(axes)))
        @argcheck all(isfinite.(last.(axes)))
        @argcheck all(first.(axes) .<= last.(axes))
        @argcheck all(ndims.(axes) .== 1)
        @argcheck all(eachindex.(axes) .== Base.axes(values))
        ret = LinearInterpolations.unsafe_Interpolate(combine, axes, values, extrapolate)
        return ret
    end
    function LinearInterpolations.unsafe_Interpolate(combine, axes::Tuple, values, extrapolate)
        C = typeof(combine)
        A = typeof(axes)
        V = typeof(values)
        O = typeof(extrapolate)
        return new{C,A,V,O}(combine, axes, values, extrapolate)
    end
end
function Adapt.adapt_structure(to, itp::Interpolate)
    unsafe_Interpolate(
        Adapt.adapt_structure(to,itp.combine),
        Adapt.adapt_structure(to,itp.axes),
        Adapt.adapt_structure(to,itp.values),
        Adapt.adapt_structure(to,itp.extrapolate),
    )
end

function Interpolate(
    xs::AbstractVector,
    ys::AbstractVector;
    combine = combine,
    extrapolate = :error,
)
    return Interpolate(combine, (xs,), ys, extrapolate)
end

"""

    itp = Interpolate(axes, values; combine = LinearInterpolations.combine, extrapolate = :error)

Create an `Interpolate` from the given arguments. The interpolate `itp` can be evaluated on points. E.g. `itp([1,2])`.

# Arguments

* axes: The coordinates of the grid points.
* values: An AbstractArray which whose size is consistent with the length each axis.
* extrapolate: How the interpolate behaves on points outside the grid limits. Possible values are:
  - [`Replicate`](@ref), `:replicate`
  - [`Reflect`](@ref), `:reflect`
  - [`Error`](@ref), `:error`
  - [`Fuzzy`](@ref), `:fuzzy`
  - [`WithPoint`](@ref)
  - [`Constant`](@ref)
* combine: A custom function that combines weights and values into the final result.
See [`LinearInterpolations.combine`](@ref)
"""
function Interpolate(axes, values; combine = combine, extrapolate = :error)
    return Interpolate(combine, axes, values, extrapolate)
end

axestype(o::Interpolate) = axestype(typeof(o))
axestype(::Type{<:Interpolate{C,A}}) where {C,A} = A
_length(::Type{<:NTuple{N,Any}}) where {N} = N
Base.ndims(L::Type{<:Interpolate}) = _length(axestype(L))
Base.ndims(o::Interpolate) = ndims(typeof(o))

NDims(o::Interpolate) = NDims(typeof(o))
NDims(::Type{<:Interpolate{C, <: NTuple{N, Any}}}) where {C,N} = Val(N)

function _make_NTuple(itr, ::Val{0})
    return ()
end
function _make_NTuple(itr, ::Val{1})
    x1, = itr
    return (x1,)
end
function _make_NTuple(itr, ::Val{2})
    x1,x2 = itr
    promote(x1,x2)
end
function _make_NTuple(itr, ::Val{3})
    x1,x2,x3 = itr
    promote(x1,x2,x3)
end
function _make_NTuple(itr, ::Val{4})
    x1,x2,x3,x4 = itr
    promote(x1,x2,x3,x4)
end
function _make_NTuple(itr, ::Val{N}) where {N}
    ret = promote(NTuple{N,Any}(itr)...)
    ret::NTuple{N, first(ret)}
end

tupelize(itp, pt) = _make_NTuple(pt, NDims(itp))

function check_dims(itp::Itp, pt) where {Itp <: Interpolate}
    _check_dims(itp, pt, NDims(Itp))
end
_check_dims(itp, pt::NTuple{N,Any}, ::Val{N}) where {N} = nothing
_check_dims(itp, pt::Number, ::Val{1}) = nothing
_check_dims(itp, pt::StaticArrays.StaticVector{N}, ::Val{N}) where {N} = nothing
function _check_dims(itp, pt, ndims_itp)
    if ndims(itp) != length(pt)
        msg = """
        Dimensions of interpolate and point do not match. Got:
        length(pt) = $(length(pt))
        ndims(itp) = $(ndims(itp))
        """
        throw(ArgumentError(msg))
    end
end

function (itp::Interpolate)(pt)
    check_dims(itp, pt)
    pt2 = tupelize(itp, pt)
    return itp(pt2)
end

function (itp::Interpolate)(pt::Tuple)
    dispatch_extrapol(itp, pt, itp.extrapolate)
end

function eval_interpolate(itp, pt::Tuple)
    dispatch_extrapol(itp, pt, itp.extrapolate)
end

function dispatch_extrapol(itp, pt::Tuple, extrapolate)
    nbs, wts = _neighbors_and_weights(itp.axes, pt, extrapolate)
    objs = NeighborsArray(itp.values, nbs)
    itp.combine(wts, objs)
end

function isinside(pt::Tuple, axes::Tuple)
    all(
        map(pt, axes) do x, r
            first(r) <= x <= last(r)
        end
    )
end

function dispatch_extrapol(itp, pt::Tuple, extrapolate::Union{Number, AbstractArray})
    return dispatch_extrapol(itp,pt,Constant(extrapolate))
end
function dispatch_extrapol(itp, pt::Tuple, extrapolate::Union{WithPoint,
                                                                        Constant,
                                                                        })
    if isinside(pt, itp.axes)
        return dispatch_extrapol(itp, pt, AssumeInside())
    else
        T = typeof(dispatch_extrapol(itp, pt, Replicate()))
        convert(T, apply_extrapolate(extrapolate, pt))
    end
end

end


text:
module TestDescentDirections

using LinearAlgebra
using Test
using OptimizationAlgorithms
using OptimizationAlgorithms: fixedpoint!, StoppingCriterion, update!, optimize!
using ForwardDiff: derivative, gradient, hessian

# TODO: separate test problems from algorithms
# TODO: tests for GaussNewton, NaturalGradient
@testset "gradient descent" begin
    n = 2
    A = randn(n, n)
    A = A'A + I
    μ = randn(n)
    f = x -> (x-μ)'A*(x-μ)
    x = randn(n)
    G = OptimizationAlgorithms.ScaledGradient(f, x)
    ε = 1e-6
    x, t = fixedpoint!(G, x, StoppingCriterion(x, dx = 1e-2ε))
    # x, t = optimize!(G, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε
end

@testset "Newton's method" begin
    # one step convergence for quadratics
    f = x -> x^2
    x = randn()
    N = OptimizationAlgorithms.Newton(f, x)
    x += N(x)
    @test derivative(f, x) ≈ 0

    # non-quadratic problem
    f = x -> x^2 + sin(x)
    x = 1.
    N = OptimizationAlgorithms.Newton(f, x)
    for i in 1:8
        x += N(x)
    end
    @test derivative(f, x) ≈ 0

    # one-step convergence for higher-dimensional quadratic
    n = 3
    A = randn(n, n) + I
    A = A'A
    μ = randn(n)
    f = x -> (x-μ)'A*(x-μ)
    ε = 1e-12

    x = randn(n)
    N = OptimizationAlgorithms.Newton(f, x)
    x += N(x)
    @test norm(gradient(f, x)) < ε

    # saddle-free Newton
    x = randn(n)
    SFN = OptimizationAlgorithms.SaddleFreeNewton(f, x)
    x += SFN(x)
    @test norm(gradient(f, x)) < ε

    # higher-dimensional non-quadratic
    f = x -> 1/2*(x-μ)'A*(x-μ) - sum(x->log(abs2(x)), x)

    x = 2randn(n)
    N = OptimizationAlgorithms.Newton(f, x)
    ε = 1e-6
    fixedpoint!(N, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε

    x = 2randn(n)
    SFN = OptimizationAlgorithms.SaddleFreeNewton(f, x)
    fixedpoint!(SFN, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε

    # test non-convex iteration
    f = x -> sum(cos, x)
    x = 1e-6randn(n)
    SFN = OptimizationAlgorithms.SaddleFreeNewton(f, x)
    fixedpoint!(SFN, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε
    @test isposdef(hessian(f, x)) # not a saddle point or local maximum
end

@testset "(L)BFGS" begin
    # convergence for quadratic
    n = 3
    A = randn(n, n)
    A = A'A
    μ = randn(n)
    f = x -> (x-μ)'A*(x-μ)

    x = randn(n)
    B = OptimizationAlgorithms.BFGS(f, x)
    ε = 1e-6
    x, t = fixedpoint!(B, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε

    # BFGS with ill-conditioned problem
    x = randn(n)
    α = [1e3, 1., 1e-3]
    g = x -> f(α.*x)
    B = OptimizationAlgorithms.BFGS(g, x)
    x, t = fixedpoint!(B, x, StoppingCriterion(x, dx = 1e-3ε, maxiter = 256))
    # test BFGS where inverse Hessian is initialized to Diagonal(1.0./α.^2)
    xs = randn(n)
    BS = OptimizationAlgorithms.BFGS(g, xs, Matrix(Diagonal(1.0./α.^2)))
    xs, ts = fixedpoint!(BS, xs, StoppingCriterion(x, dx = 1e-3ε, maxiter = 256))
    @test norm(gradient(g, x)) < ε
    @test norm(gradient(g, xs)) < ε

    # @test ts < t # scaling outperforms non-scaling for ill-conditioned problem
    # # WARNING: This test could fail, but most times doesn't

    # OptimizationAlgorithms of 1D non-convex objective (has to forgo BFGS update)
    x = [2.] # fill(2., 2)
    h = x -> 1-exp(-sum(abs2, x))
    B = OptimizationAlgorithms.BFGS(h, x, check = false)
    x, t = fixedpoint!(B, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(h, x)) < ε

    x = [2.] # fill(2., 2)
    B = OptimizationAlgorithms.LBFGS(h, x, 3, check = false)
    x, t = fixedpoint!(B, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(h, x)) < ε

    ###################### limited-memory BFGS ##################################
    m = 3
    x = randn(n)
    B = OptimizationAlgorithms.LBFGS(f, x, m)
    ε = 1e-6
    x, t = fixedpoint!(B, x, StoppingCriterion(x, dx = 1e-2ε))
    # println(t)
    @test norm(gradient(f, x)) < ε

    # higher-dimensional non-quadratic
    f = x -> 1/2*(x-μ)'A*(x-μ) - sum(x->log(abs2(x)), x)
    x = randn(n)
    B = OptimizationAlgorithms.BFGS(f, x)
    ε = 1e-6
    x, t = fixedpoint!(B, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε

    x = randn(n)
    B = OptimizationAlgorithms.LBFGS(f, x, m)
    ε = 1e-6
    x, t = fixedpoint!(B, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε
end

@testset "CustomDirection" begin
    # convergence for quadratic
    n = 3
    A = randn(n, n)
    A = 1/2*A'A + I
    μ = randn(n)
    f = x -> (x-μ)'A*(x-μ)
    function valdir(x::AbstractVector) # shows how evaluation and gradient computation can be pooled
        g = A*(x-μ)
        v = (x-μ)'g
        g .*= -g'g / (g'A*g) # optimal stepsize for quadratic problems
        v, g
    end
    x = 100randn(n) + μ # normally distributed around μ
    B = OptimizationAlgorithms.CustomDirection(f, valdir, x)
    ε = 1e-6
    fixedpoint!(B, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε

    x = 1e2(randn(n) .+ 1)
    @test norm(B(x)) > 1 # normal direction is large

    # testing direction modifiers
    U = OptimizationAlgorithms.UnitDirection(B)
    @test U(x) isa AbstractVector
    @test norm(U(x)) ≈ 1
    @test dot(U(x), B(x)) ≈ norm(B(x)) # co-linear with original direction

    TD1 = OptimizationAlgorithms.TrustedDirection(B, 1., Inf) # testing norm threshold
    TD2 = OptimizationAlgorithms.TrustedDirection(B, Inf, 1.) # testing elementwise threshold
    @test TD1(x) isa AbstractVector
    @test norm(TD1(x)) ≈ 1
    @test norm(TD2(x)) > 1 && maximum(abs, TD2(x)) ≈ 1

    @. x = μ + 1e-6*randn()
    @test norm(TD1(x)) ≈ norm(B(x)) # no change if gradient is sufficiently small
    @test norm(TD2(x)) ≈ norm(B(x)) # no change if gradient is sufficiently small

end

end # TestDescentDirections


text:
module Model
import Base: ==
using StructTypes, Dates

export Paper

mutable struct Paper
    id::Int64 # service-managed
    titile::String
    abst::String
    authors::Vector{String}
    publishdata::Date
    timespicked::Int64 # service-managed
end

==(x::Paper, y::Paper) = x.id == y.id
Paper() = Paper(0, "", "", String[], Date("1900-01-01", "y-m-d"), 0)
Paper(title, abst, authors, publishdate) = Paper(0,
                                                 title,
                                                 abst,
                                                 authors,
                                                 Date(publishdate, "y-m-d"),
                                                 0)

StructTypes.StructType(::Type{Paper}) = StructTypes.Mutable()

end


text:
using Pkg: TOML
using SHA

# Not resource serving, but other interesting dynamic things that the PkgServer can do.
const artifact_toml_re = Regex("""^(/artifact/$hash_re)/([^/]+)\$""")

function serve_artifact_toml(http::HTTP.Stream,
                             artifact_path::AbstractString,
                             artifact_name::AbstractString)
    artifact_tree_hash = basename(artifact_path)
    artifact_sha256_hash = bytes2hex(open(io -> SHA.sha256(io), artifact_path))
    # build a TOML dict to serialize and send out:
    toml = Dict(
        artifact_name => Dict(
            "git-tree-sha1" => artifact_tree_hash,
            "download" => [Dict(
                "sha256" => artifact_sha256_hash,
                "url" => joinpath(get_pkgserver_url(), "artifact", artifact_tree_hash),
            )]
        )
    )
    io = IOBuffer()
    TOML.print(io, toml)
    toml_string = String(take!(io))

    HTTP.setheader(http, "Content-Length" => string(length(toml_string)))
    HTTP.setheader(http, "Content-Type" => "application/toml")
    HTTP.setheader(http, "Content-Disposition" => "attachment; filename=Artifacts.toml")
    startwrite(http)
    return write(http, toml_string)
end

max_stars_repo_name: StefanKarpinski/PkgServer

max_issues_repo_name: StefanKarpinski/PkgServer

max_forks_repo_name: StefanKarpinski/PkgServer


text:
using DataFrames
using HDFS
using Blocks

function find_nrows(datafile)
    b = Block(datafile)
    b |> as_recordio |> as_bufferedio |> x->as_dataframe(x, separator='\t');

    pmapreduce(nrow, +, b)
end

function filter_smiley(df, s)
    (df[:(x4 .== ":)"), :])[:x3]
end

function find_smiley_counts(datafile)
    b = Block(datafile)
    b |> as_recordio |> as_bufferedio
    dt = dreadtable(b, separator='\t')
    by(dt, [:x4], x->sum(x[:x3]), sum)
end

function run_tests()
    datafile = HdfsURL("hdfs://localhost:9000/user/tan/streaminp/smileys_for_julia.tsv")
    nrows = find_nrows(datafile)
    smiley_counts = find_smiley_counts(datafile)

    println("num rows: $nrows")
    println(smiley_counts)
end


text:
import RobotDynamics: jacobian!

""""
Specifies whether the constraint is an equality or inequality constraint.
Valid subtypes are `Equality`, and `Inequality`.
The sense of a constraint can be queried using `sense(::AbstractConstraint)`
"""
abstract type ConstraintSense end
"""
Inequality constraints of the form ``h(x) \\leq 0``.
Type singleton, so it is created with `Inequality()`
"""
struct Equality <: ConstraintSense end
"""
Equality constraints of the form ``g(x) = 0`.
Type singleton, so it is created with `Equality()`.
"""
struct Inequality <: ConstraintSense end

"""
    AbstractConstraint

Abstract vector-valued constraint for a trajectory optimization problem.
May be either inequality or equality (specified by `sense(::AbstractConstraint)::ConstraintSense`),
and be function of single or adjacent knotpoints.

Interface:
Any constraint type must implement the following interface:
```julia
n = state_dim(::MyCon)
m = control_dim(::MyCon)
p = Base.length(::MyCon)
sense(::MyCon)::ConstraintSense
c = evaluate(::MyCon, args...)
jacobian!(∇c, ::MyCon, args...)
```

All constraints are categorized into the following type tree:
```text
                        AbstractConstraint
                        ↙                ↘
           StageConstraint               CoupledConstraint
            ↙        ↘                       ↙           ↘
StageConstraint ControlConstraint CoupledStateConstraint CoupledControlConstraint
```

The state and control dimensions (where applicable) can be queried using
`state_dim(::AbstractConstraint)` and `control_dim(::AbstractConstraint)`.
The dimensions of a constraint can be verified using [`check_dims`](@ref).
The width of the constraint Jacobian is given by [`get_inds`](@ref) or [`widths`](@ref).

The number of constraint values associated with the constraint (length of the constraint vector)
is given with `length(::AbstractConstraint)`.

# Evaluation methods
Refer to the doc strings for the following methods for more information on the required
signatures.
* [`evaluate`](@ref)
* [`jacobian!`](@ref)
* [`∇jacobian!`](@ref)
"""
abstract type AbstractConstraint end

"Only a function of states and controls at a single knotpoint"
abstract type StageConstraint <: AbstractConstraint end
"Only a function of states at a single knotpoint"
abstract type StateConstraint <: StageConstraint end
"Only a function of controls at a single knotpoint"
abstract type ControlConstraint <: StageConstraint end
"Only a function of states and controls at two adjacent knotpoints"
abstract type CoupledConstraint <: AbstractConstraint end
"Only a function of states at adjacent knotpoints"
abstract type CoupledStateConstraint <: CoupledConstraint end
"Only a function of controls at adjacent knotpoints"
abstract type CoupledControlConstraint <: CoupledConstraint end

const StateConstraints =
    Union{StageConstraint,StateConstraint,CoupledConstraint,CoupledStateConstraint}
const ControlConstraints =
    Union{StageConstraint,ControlConstraint,CoupledConstraint,CoupledControlConstraint}

"Get constraint sense (Inequality vs Equality)"
sense(::C) where {C<:AbstractConstraint} = throw(NotImplemented(:sense, Symbol(C)))

"Dimension of the state vector"
RobotDynamics.state_dim(::C) where {C<:StateConstraint} =
    throw(NotImplemented(:state_dim, Symbol(C)))

"Dimension of the control vector"
RobotDynamics.control_dim(::C) where {C<:ControlConstraint} =
    throw(NotImplemented(:control_dim, Symbol(C)))

"Return the constraint value"
evaluate(::C) where {C<:AbstractConstraint} = throw(NotImplemented(:evaluate, Symbol(C)))

"Length of constraint vector"
Base.length(::C) where {C<:AbstractConstraint} = throw(NotImplemented(:length, Symbol(C)))

# widths(con::StageConstraint, n=state_dim(con), m=control_dim(con)) = (n+m,)
# widths(con::StateConstraint, n=state_dim(con), m=0) = (n,)
# widths(con::ControlConstraint, n=0, m=control_dim(con)) = (m,)
# widths(con::CoupledConstraint, n=state_dim(con), m=control_dim(con)) = (n+m, n+m)
# widths(con::CoupledStateConstraint, n=state_dim(con), m=0) = (n,n)
# widths(con::CoupledControlConstraint, n=0, m=control_dim(con)) = (m,m)

"Upper bound of the constraint, as a vector, which is 0 for all constraints
(except bound constraints)"
@inline upper_bound(con::AbstractConstraint) =
    upper_bound(sense(con)) * @SVector ones(length(con))
@inline upper_bound(::Inequality) = 0.0
@inline upper_bound(::Equality) = 0.0

"Upper bound of the constraint, as a vector, which is 0 equality and -Inf for inequality
(except bound constraints)"
@inline lower_bound(con::AbstractConstraint) =
    lower_bound(sense(con)) * @SVector ones(length(con))
@inline lower_bound(::Inequality) = -Inf
@inline lower_bound(::Equality) = 0.0

"""
    primal_bounds!(zL, zU, con::AbstractConstraint)

Set the lower `zL` and upper `zU` bounds on the primal variables imposed by the constraint
`con`. Return whether or not the vectors `zL` or `zU` could be modified by `con`
(i.e. if the constraint `con` is a bound constraint).
"""
primal_bounds!(zL, zU, con::AbstractConstraint) = false

"Is the constraint a bound constraint or not"
@inline is_bound(con::AbstractConstraint) = false

"Check whether the constraint is consistent with the specified state and control dimensions"
@inline check_dims(con::StateConstraint, n, m) = state_dim(con) == n
@inline check_dims(con::ControlConstraint, n, m) = control_dim(con) == m
@inline check_dims(con::AbstractConstraint, n, m) =
    state_dim(con) == n && control_dim(con) == m

get_dims(con::Union{StateConstraint,CoupledStateConstraint}, nm::Int) =
    state_dim(con), nm - state_dim(con)
get_dims(con::Union{ControlConstraint,CoupledControlConstraint}, nm::Int) =
    nm - control_dim(con), control_dim(con)
get_dims(con::AbstractConstraint, nm::Int) = state_dim(con), control_dim(con)

con_label(::AbstractConstraint, i::Int) = "index $i"

"""
    get_inds(con::AbstractConstraint)

Get the indices of the joint state-control vector that are used to calculate the constraint.
If the constraint depends on more than one time step, the indices start from the beginning
of the first one.
"""
get_inds(con::StateConstraint, n, m) = (1:n,)
get_inds(con::ControlConstraint, n, m) = (n .+ (1:m),)
get_inds(con::StageConstraint, n, m) = (1:n+m,)
get_inds(con::CoupledConstraint, n, m) = (1:n+m, n+m+1:2n+2m)

"""
    widths(::AbstractConstraint)
    widths(::AbstractConstraint, n, m)

Return a tuple of the widths of the Jacobians for a constraint. If `n` and `m` are not passed
in, they are assumed to be consistent with those returned by `state_dim` and `control_dim`.
"""
@inline widths(con::AbstractConstraint, n, m) = length.(get_inds(con, n, m))
@inline widths(con::StageConstraint) = (state_dim(con) + control_dim(con),)
@inline widths(con::StateConstraint) = (state_dim(con),)
@inline widths(con::ControlConstraint) = (control_dim(con),)
@inline widths(con::CoupledConstraint) =
    (state_dim(con) + control_dim(con), state_dim(con) + control_dim(con))
@inline widths(con::CoupledStateConstraint) = (state_dim(con), state_dim(con))
@inline widths(con::CoupledControlConstraint) = (control_dim(con), control_dim(con))

############################################################################################
# 								EVALUATION METHODS 										   #
############################################################################################
"""
    evaluate!(vals, con::AbstractConstraint, Z, [inds])

Evaluate constraints for entire trajectory. This is the most general method used to evaluate
constraints along the trajectory `Z`, and should be the one used in other functions.
The `inds` argument determines at which knot points the constraint is evaluated.

If `con` is a `StageConstraint`, this will call `evaluate(con, z)` by default, or
`evaluate(con, z1, z2)` if `con` is a `CoupledConstraint`.
"""
function evaluate!(
    vals::Vector{<:AbstractVector},
    con::StageConstraint,
    Z::AbstractTrajectory,
    inds = 1:length(Z),
)
    for (i, k) in enumerate(inds)
        vals[i] .= evaluate(con, Z[k])
    end
end

function evaluate!(
    vals::Vector{<:AbstractVector},
    con::CoupledConstraint,
    Z::AbstractTrajectory,
    inds = 1:length(Z)-1,
)
    for (i, k) in enumerate(inds)
        vals[i] .= evaluate(con, Z[k], Z[k+1])
    end
end

"""
    jacobian!(∇c, con::AbstractConstraint, Z, [inds])

Evaluate constraints for entire trajectory. This is the most general method used to evaluate
constraints along the trajectory `Z`, and should be the one used in other functions.
The `inds` argument determines at which knot points the constraint is evaluated.

The values are stored in `∇c`, which should be a matrix of matrices. If `con` is a
`StageConstraint`, `size(∇c,2) = 1`, and `size(∇c,2) = 2` if `con` is a `CoupledConstraint`.

If `con` is a `StageConstraint`, this will call `jacobian!(∇c, con, z)` by default, or
`jacobian!(∇c, con, z1, z2, i)` if `con` is a `CoupledConstraint`.
"""
function jacobian!(
    ∇c::VecOrMat{<:AbstractMatrix},
    con::StageConstraint,
    Z::AbstractTrajectory,
    inds = 1:length(Z),
)
    for (i, k) in enumerate(inds)
        jacobian!(∇c[i], con, Z[k])
    end
end

function jacobian!(
    ∇c::VecOrMat{<:AbstractMatrix},
    con::CoupledConstraint,
    Z::AbstractTrajectory,
    inds = 1:size(∇c, 1),
)
    for (i, k) in enumerate(inds)
        jacobian!(∇c[i, 1], con, Z[k], Z[k+1], 1)
        jacobian!(∇c[i, 2], con, Z[k], Z[k+1], 2)
    end
end

"""
    ∇jacobian!(G, con::AbstractConstraint, Z, λ, inds, is_const, init)
    ∇jacobian!(G, con::AbstractConstraint, Z::AbstractKnotPoint, λ::AbstractVector)

Evaluate the second-order expansion of the constraint `con` along the trajectory `Z`
after multiplying by the lagrange multiplier `λ`.

The method for each constraint should calculate the Jacobian of the vector-Jacobian product,
    and therefore should be of size n × n if the input dimension is n.

Importantly, this method should ADD and not overwrite the contents of `G`, since this term
is dependent upon all the constraints acting at that time step.
"""
function ∇jacobian!(
    G::VecOrMat{<:AbstractMatrix},
    con::StageConstraint,
    Z::AbstractTrajectory,
    λ::Vector{<:AbstractVector},
    inds = 1:length(Z),
    is_const = ones(Bool, length(inds)),
    init::Bool = false,
)
    for (i, k) in enumerate(inds)
        if init || !is_const[i]
            is_const[i] = ∇jacobian!(G[i], con, Z[k], λ[i])
        end
    end
end

function ∇jacobian!(
    G::VecOrMat{<:AbstractMatrix},
    con::CoupledConstraint,
    Z::AbstractTrajectory,
    λ::Vector{<:AbstractVector},
    inds = 1:length(Z),
    is_const = ones(Bool, length(inds)),
    init::Bool = false,
)
    for (i, k) in enumerate(inds)
        if init || !is_const[i]
            is_const[i] = ∇jacobian!(G[i, 1], con, Z[k], Z[k+1], λ[i], 1)
            is_const[i] = ∇jacobian!(G[i, 2], con, Z[k], Z[k+1], λ[i], 2)
        end
    end
end

# Default methods for converting KnotPoints to states and controls for StageConstraints
"""
    evaluate(con::AbstractConstraint, z)       # stage constraint
    evaluate(con::AbstractConstraint, z1, z2)  # coupled constraint

Evaluate the constraint `con` at knot point `z`. By default, this method will attempt to call

    evaluate(con, x)

if `con` is a `StateConstraint`,

    evaluate(con, u)

if `con` is a `ControlConstraint`, or

    evaluate(con, x, u)

if `con` is a `StageConstraint`. If `con` is a `CoupledConstraint` the constraint should
define

    evaluate(con, z1, z2)

"""
@inline evaluate(con::StateConstraint, z::AbstractKnotPoint) = evaluate(con, state(z))
@inline evaluate(con::ControlConstraint, z::AbstractKnotPoint) = evaluate(con, control(z))
@inline evaluate(con::StageConstraint, z::AbstractKnotPoint) =
    evaluate(con, state(z), control(z))

"""
    jacobian!(∇c, con::AbstractConstraint, z, i=1)       # stage constraint
    jacobian!(∇c, con::AbstractConstraint, z1, z2, i=1)  # coupled constraint

Evaluate the constraint `con` at knot point `z`. By default, this method will attempt to call

    jacobian!(∇c, con, x)

if `con` is a `StateConstraint`,

    jacobian!(∇c, con, u)

if `con` is a `ControlConstraint`, or

    jacobian!(∇c, con, x, u)

if `con` is a `StageConstraint`. If `con` is a `CoupledConstraint` the constraint should
define

    jacobian!(∇c, con, z, i)

where `i` determines which Jacobian should be evaluated. E.g. if `i = 1`, the Jacobian
with respect to the first knot point's stage and controls is calculated.

# Automatic Differentiation
If `con` is a `StateConstraint` or `ControlConstraint` then this method is automatically
defined using ForwardDiff.
"""
jacobian!(∇c, con::StateConstraint, z::AbstractKnotPoint, i = 1) =
    jacobian!(∇c, con, state(z))
jacobian!(∇c, con::ControlConstraint, z::AbstractKnotPoint, i = 1) =
    jacobian!(∇c, con, control(z))
jacobian!(∇c, con::StageConstraint, z::AbstractKnotPoint, i = 1) =
    jacobian!(∇c, con, state(z), control(z))

# ForwardDiff jacobians that are of only state or control
function jacobian!(∇c, con::StageConstraint, x::StaticVector)
    eval_c(x) = evaluate(con, x)
    ∇c .= ForwardDiff.jacobian(eval_c, x)
    return false
end

@inline ∇jacobian!(G, con::StateConstraint, z::AbstractKnotPoint, λ, i = 1) =
    ∇jacobian!(G, con, state(z), λ)
@inline ∇jacobian!(G, con::ControlConstraint, z::AbstractKnotPoint, λ, i = 1) =
    ∇jacobian!(G, con, control(z), λ)
@inline ∇jacobian!(G, con::StageConstraint, z::AbstractKnotPoint, λ, i = 1) =
    ∇jacobian!(G, con, state(z), control(z), λ)

function ∇jacobian!(G, con::StageConstraint, x::StaticVector, λ)
    eval_c(x) = evaluate(con, x)'λ
    G_ = ForwardDiff.hessian(eval_c, x)
    G .+= G_
    return false
end

function ∇jacobian!(
    G,
    con::StageConstraint,
    x::StaticVector{n},
    u::StaticVector{m},
    λ,
) where {n,m}
    ix = SVector{n}(1:n)
    iu = SVector{m}(n .+ (1:m))
    eval_c(z) = evaluate(con, z[ix], z[iu])'λ
    G .+= ForwardDiff.hessian(eval_c, [x; u])
    return false
end

function gen_jacobian(con::AbstractConstraint, i = 1)
    ws = widths(con)
    p = length(con)
    C1 = SizedMatrix{p,ws[i]}(zeros(p, ws[i]))
end

function gen_views(∇c::AbstractMatrix, con::StateConstraint, n = state_dim(con), m = 0)
    view(∇c, :, 1:n), view(∇c, :, n:n-1)
end

function gen_views(∇c::AbstractMatrix, con::ControlConstraint, n = 0, m = control_dim(con))
    view(∇c, :, 1:0), view(∇c, :, 1:m)
end

function gen_views(
    ∇c::AbstractMatrix,
    con::AbstractConstraint,
    n = state_dim(con),
    m = control_dim(con),
)
    if size(∇c, 2) < n + m
        view(∇c, :, 1:n), view(∇c, :, n:n-1)
    else
        view(∇c, :, 1:n), view(∇c, :, n .+ (1:m))
    end
end


text:
function worldcupstages()
    games = ["12", "13", "14", "23", "24", "34"]
    results = "000000"

    function nextresult()
        if (results == "222222")
            return false
        end
        results = lpad(string(parse(Int, results, base=3) + 1, base=3), 6, '0')
        true
    end

    points = zeros(Int, 4, 10)
    while true
        records = zeros(Int, 4)
        for i in 1:length(games)
            if results[i] == '2'
                records[games[i][1] - '0'] += 3
            elseif results[i] == '1'
                records[games[i][1] - '0'] += 1
                records[games[i][2] - '0'] += 1
            elseif results[i] == '0'
                records[games[i][2] - '0'] += 3
            end
        end
        sort!(records)
        for i in 1:4
            points[i, records[i] + 1] += 1
        end
        if !nextresult()
            break
        end
    end

    for (i, place) in enumerate(["First", "Second", "Third", "Fourth"])
        println("$place place: $(points[5 - i, :])")
    end
end

worldcupstages()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using TenX,CSV
using Test

dn = joinpath(@__DIR__,"data")
gzd = joinpath(dn,"gz")
tsd = joinpath(dn,"tsv")
checkfn = joinpath(dn,"check.csv")

@testset "TenX.jl" begin
	gdf = TenX.loadFolder(gzd)

	check = CSV.read(checkfn,DataFrame)
    @test Matrix(check) == Matrix(gdf[:,2:end])

	tdf = TenX.loadFolder(tsd)
	@test Matrix(check) == Matrix(tdf[:,2:end])
end


text:
include(joinpath("..", "src", "SimpleMolecule.jl"))
include(joinpath("Common", "Common.jl"))
include(joinpath("TestAtom", "TestAtom.jl"))
include(joinpath("TestMolecule", "TestMolecule.jl"))

TestAtom.run_tests()
TestMolecule.run_tests()


text:
function string_time(str::String, T::Array{Y,1}) where Y<:Integer
  if isdigit(str[end])
    s = str*"~"
    io = IOBuffer(s)
  else
    io = IOBuffer(str)
  end
  seekstart(io)
  t = stream_time(io, T)
  close(io)
  return t
end
string_time(str::String) = string_time(str, BUF.date_buf)


text:
# Julia wrapper for header: /usr/include/scip/cons_components.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrComponents(scip)
    ccall((:SCIPincludeConshdlrComponents, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 120da254-1c28-11ec-1938-a55c8cde5727
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ c99f9fc5-0ab1-444a-80f6-b5241edf2f8f
begin
	using QuantumStateBase
	using Plots
	using LinearAlgebra
end

# ╔═╡ db69bedc-2803-4268-9f56-34e4a88a324e
state = SqueezedThermalState(Complex{Float64}, ξ(0.5, π/4), 0.5, )

# ╔═╡ cd19d00d-bf2b-479d-8c88-493311d97303
d = gaussian_state_sampler(Float64, state, 4096)

# ╔═╡ b9ea40f4-1ed1-4862-8858-be07f2c04633
scatter(d[1, :], d[2, :])

# ╔═╡ Cell order:
# ╟─120da254-1c28-11ec-1938-a55c8cde5727
# ╠═c99f9fc5-0ab1-444a-80f6-b5241edf2f8f
# ╠═db69bedc-2803-4268-9f56-34e4a88a324e
# ╠═cd19d00d-bf2b-479d-8c88-493311d97303
# ╠═b9ea40f4-1ed1-4862-8858-be07f2c04633