text:
# change the stored lendiv to len b/c it's more intuitive to work around length than the
# inferred divisor
"""
    StaticLinRange{T,B,E,L,D}
"""
abstract type StaticLinRange{T,B,S,E,L,D} <: StaticOrdinalRange{T,B,S,E,L} end

@pure lendiv(::StaticLinRange{T,B,E,L,D}) where {T,B,E,L,D<:SReal} = D()
@pure lendiv(::Type{<:StaticLinRange{T,B,E,L,D}}) where {T,B,E,L,D<:SReal} = D()

# TODO: changed mind, length! should change compute and change lendiv
"""
    LinSRange

Static parametric type variant of StaticLinRange
"""
struct LinSRange{T,B,S,E,L,D} <: StaticLinRange{T,B,S,E,L,D}
    LinSRange{T}(start::B, step::S, stop::E, len::L, lendiv::D) where {T,B,S,E,L,D} = new{T,B,S,E,L,D}()
end

LinSRange{T}(start::B, stop::E, len::L) where {T,B,E,L} =
    LinSRange{T}(start, stop, len, max(len-SOne, SOne))
LinSRange{T}(start::B, stop::E, len::L, lendiv::D) where {T,B,E,L,D} =
    LinSRange{T}(start, (stop-start)/lendiv, stop, len, lendiv)

"""
    LinMRange

Mutable variant of StaticLinRange. Notable difference from `LinRange` is the absence of a
`lendiv` field. This allows the length to be changed dynamically with `length!(range, length)`.
"""
mutable struct LinMRange{T,B,E,L} <: StaticLinRange{T,B,Dynamic,E,L,Dynamic}
    start::B
    stop::E
    len::L
end

first(r::LinMRange{T,B,E}) where {T,B<:BaseReal,E} = getfield(r, :start)

step(r::LinMRange) = (last(r)-first(r))/lendiv(r)

last(r::LinMRange{T,B,E}) where {T,B,E<:BaseReal} = getfield(r, :stop)

length(r::LinMRange{T,B,E,<:BaseInteger}) where {T,B,E} = getfield(r, :len)

length!(r::LinMRange{T,B,E,L}, len::L) where {T,B,E,L<:BaseNumber} = setfield!(r, :len, len)

lendiv(r::LinMRange) = max(length(r)-SOne, SOne)

isstatic(::LinSRange) = true
isstatic(::Type{<:LinSRange}) = true

isstatic(::LinMRange) = false
isstatic(::Type{<:LinMRange}) = false

#LinMRange{T}(start::B, stop::E, len::L, lendiv::D) where {T,B,E,L,D} =
#    LinMRange{T,B,E,L}(start, stop, len, lendiv)

function LinMRange{T}(start::B, stop::E, len::L) where {T,B,E,L}
    len >= 0 || throw(ArgumentError("srange($start, stop=$stop, length=$len): negative length"))
    if len == 0
        start == stop || throw(ArgumentError("srange($start, stop=$stop, length=$len): endpoints differ"))
        return LinMRange{T,B,E,L}(start, stop, len)
    end
    return LinMRange{T,B,E,L}(start, stop, len)
end

function StaticLinRange(start, stop, len::Integer)
    StaticLinRange{eltype((stop-start)/len)}(start, stop, len)
end

for (S,B) in S2B
    @eval begin
        function StaticLinRange{T}(start::$S, stop::$S, len::SInteger) where T
            len >= 0 || throw(ArgumentError("srange($start, stop=$stop, length=$len): negative length"))
            LinSRange{T}(start, stop, len, max(len-SOne, SOne))
        end

        function StaticLinRange{T}(start::$S, stop::$S, len::SIntegerZeroType) where T
            LinSRange{T}(start, stop, len, SOne, SOne)
        end

        function StaticLinRange{T}(start::$S{X}, stop::$S{X}, len::SIntegerZeroType) where {T,X}
            throw(ArgumentError("range($start, stop=$stop, length=$len): endpoints differ"))
        end

        function StaticLinRange{T}(start::$B, stop::$B, len::BaseInteger) where {T}
            LinMRange{T}(start, stop, len)
        end

        function StaticLinRange{T}(start::$S, stop::$S, len::BaseInteger) where {T}
            LinMRange{T}(start, stop, len)
        end

        function StaticLinRange{T}(start::$B, stop::$S, len::BaseInteger) where {T}
            LinMRange{T}(start, stop, len)
        end

        function StaticLinRange{T}(start::$S, stop::$B, len::BaseInteger) where {T}
            LinMRange{T}(start, stop, len)
        end

        function StaticLinRange{T}(start::$B, stop::$S, len::SInteger) where {T}
            LinMRange{T}(start, stop, len)
        end

        function StaticLinRange{T}(start::$S, stop::$B, len::SInteger) where {T}
            LinMRange{T}(start, stop, len)
        end

    end
end

@inline function Base.iterate(r::LinSRange{T,B,S,E,L,D}) where {T,B,S,E,L,D}
    if values(L) < 1
        return nothing
    else
        return (unsafe_getindex(r, 1)::T, 1)::Tuple{T,Int}
    end
end

@inline function Base.iterate(r::LinSRange{T,B,S,E,L,D}, i::Int) where {T,B,S,E,L,D}
    if i < values(L)
        return (unsafe_getindex(r, i), i + 1)::Tuple{T,Int}
    else
        return nothing
    end
end

#= TODO
function Base.iterate(r::LinSRange, i::Int=1)
end
=#


text:
"""
This maps (v, C, log_α) ↦ u = log.(α) - log.(Kb) = -log.(Kb./α) = log.(α./Kb)
"""
function sinkhorn_step_log!(u::AbstractArray{T}, v, C, εinv, log_α, temp) where T
    for i in eachindex(u)
        for j in eachindex(v)          # sum over j
            temp[j] = C[i,j]
        end
        rmul!(temp, -εinv)
        temp .+= v .- log_α[i]
        u[i] = -logsumexp(temp)
    end
    return u
end

"""
wasserstein_distance_log(x, y, log_α, log_β, c, u₀, v₀, log_d₁₀, log_d₂₀, SP, caches)

Differentiable by x only!
"""
function sinkhorn_dvg_particles(    x::AbstractMatrix{T},
                                    y::AbstractMatrix{V},
                                    log_α::Vector{V},
                                    log_β::Vector{V},
                                    c::Base.Callable,
                                    u₀, v₀,
                                    log_d₁₀, log_d₂₀,
                                    SP, caches   ) where {T,V}

    VC = caches.VC
    εinv = 1/SP.ε

    C_xy = ForwardDiff.value.( LazyCost(x, y, c) )
    C_xx = ForwardDiff.value.( LazyCost(x, x, c) )
    C_yy = LazyCost(y, y, c)

    VC[:u,V] .= u₀
    VC[:v,V] .= v₀
    if SP.debias
        VC[:log_d₁,V] .= log_d₁₀
        VC[:log_d₂,V] .= log_d₂₀
    end

    for l in 1:SP.L
        sinkhorn_step_log!(VC[:u,V], VC[:v,V], C_xy, εinv, log_α, VC[:t1,V]) # u = log.(α./Kb)
        sinkhorn_step_log!(VC[:v,V], VC[:u,V], C_xy', εinv, log_β, VC[:t1,V]) # v = log.(β./Ka) assuming K = Kᵀ

        if SP.debias # && (l%2==0 || l==SP.L)
            VC[:t3,V] .= VC[:log_d₁,V]
            VC[:t2,V] .= VC[:log_d₁,V] .+ log_α
            sinkhorn_step_log!(VC[:log_d₁,V], VC[:t3,V], C_xx, εinv, VC[:t2,V], VC[:t1,V]) .*= 0.5

            VC[:t3,V] .= VC[:log_d₂,V]
            VC[:t2,V] .= VC[:log_d₂,V] .+ log_β
            sinkhorn_step_log!(VC[:log_d₂,V], VC[:t3,V], C_yy, εinv, VC[:t2,V], VC[:t1,V]) .*= 0.5
        end
    end

    C_xy = LazyCost(x, y, c)
    C_xx = LazyCost(x, x, c)

    S_ε = 0

    # calculate ε⟨α,u-d₁⟩ + ε⟨β,v-d₂⟩ - ε ∑ e^(u⊕v-C/ε) + 0.5ε ∑ e^(d₁⊕d₁-C/ε) + 0.5ε ∑ e^(d₂⊕d₂-C/ε)

    for i in eachindex(log_α)
        S_ε += SP.ε * exp(log_α[i]) * VC[:u,V][i]
        if SP.debias
            S_ε -= SP.ε * exp(log_α[i]) * VC[:log_d₁,V][i]
            for i_ in eachindex(log_α)
                S_ε += 0.5*SP.ε * exp( VC[:log_d₁,V][i] + VC[:log_d₁,V][i_] - C_xx[i,i_]*εinv )
            end
        end
        for j in eachindex(log_β)
            S_ε -= SP.ε * exp( VC[:u,V][i] + VC[:v,V][j] - C_xy[i,j] * εinv )
        end
    end

    for j in eachindex(log_β)
        S_ε += SP.ε * exp(log_β[j]) * VC[:v,V][j]
        if SP.debias
            S_ε -= SP.ε * exp(log_β[j]) * VC[:log_d₂,V][j]
            for j_ in eachindex(log_β)
                S_ε += 0.5*SP.ε * exp( VC[:log_d₂,V][j] + VC[:log_d₂,V][j_] - C_yy[j,j_]*εinv )
            end
        end
    end

    if SP.update_potentials
        u₀ .= VC[:u,V]
        v₀ .= VC[:v,V]
        if SP.debias
            log_d₁₀ .= VC[:log_d₁,V]
            log_d₂₀ .= VC[:log_d₂,V]
        end
    end

    return S_ε
end

ext: jl

max_stars_repo_name: JuliaRCM/WassersteinDictionaries

max_issues_repo_name: JuliaRCM/WassersteinDictionaries

max_forks_repo_name: JuliaRCM/WassersteinDictionaries


text:
# ---
# title: 806. Number of Lines To Write String
# id: problem806
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/number-of-lines-to-write-string/description/>
# hidden: true
# ---
#
# We are to write the letters of a given string `S`, from left to right into
# lines. Each line has maximum width 100 units, and if writing a letter would
# cause the width of the line to exceed 100 units, it is written on the next
# line. We are given an array `widths`, an array where widths[0] is the width of
# 'a', widths[1] is the width of 'b', ..., and widths[25] is the width of 'z'.
#
# Now answer two questions: how many lines have at least one character from `S`,
# and what is the width used by the last such line? Return your answer as an
# integer list of length 2.
#
#
#
#
#
#     **Example :**
#     Input:
#     widths = [10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10]
#     S = "abcdefghijklmnopqrstuvwxyz"
#     Output: [3, 60]
#     Explanation:
#     All letters have the same length of 10. To write all 26 letters,
#     we need two full lines and one line with 60 units.
#
#
#
#     **Example :**
#     Input:
#     widths = [4,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10]
#     S = "bbbcccdddaaa"
#     Output: [2, 4]
#     Explanation:
#     All letters except 'a' have the same length of 10, and
#     "bbbcccdddaa" will cover 9 * 10 + 2 * 4 = 98 units.
#     For the last 'a', it is written on the second line because
#     there is only 2 units left in the first line.
#     So the answer is 2 lines, plus 4 units in the second line.
#
#
#
#
# **Note:**
#
#   * The length of `S` will be in the range [1, 1000].
#   * `S` will only contain lowercase letters.
#   * `widths` is an array of length `26`.
#   * `widths[i]` will be in the range of `[2, 10]`.
#
#
## @lc code=start
using LeetCode

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


text:
# ==============================================================================
# ITERATIVE HESSENBERG
#   * Type of A, and type of b has to be the same. For example, if A is complex,
#   then b vector/matrix has to be complex too.
#   Storage Format:
#       H: The upper Hessenberg form is stored as vectors of vectors, where
#       each element in the vector `H`, are columns of the Hessenberg matrix.

mutable struct Arnoldi

    A::Function
    b
    H::Vector  # vector of vector
    Q::Vector  # vector of orthogonal vector
    r0
    x0
    maxK::Int64
    itr_count::Int64

    function Arnoldi(A::Function, b; x0=nothing, max_k=typemax(Int64))
        this = new()
        this.A = A
        this.b = b
        this.x0 = x0 === nothing ? b : x0
        this.r0 = b - A(this.x0)
        this.H = Vector{Vector{Number}}()
        this.Q = Vector{typeof(A(b))}()
        push!(this.Q, this.r0./sqrt(dot(this.r0, this.r0)))
        this.maxK = max_k
        this.itr_count = 0
        return this
    end

    function Arnoldi(
            A::Matrix,
            b::VecOrMat;
            x0=nothing,
            max_k=typemax(Int64)
        )
        this = Arnoldi((x) -> A*x, b; x0=x0, max_k=max_k)
        return this
    end

end

## Functional Operator Override.
function (this::Arnoldi)()

    A = this.A
    Q = this.Q
    u = A(Q[end])
    h = Vector()
    for q ∈ Q
        push!(h, dot(q, u))
        u .-= h[end].*q
    end
    if norm(u) ≈ 0
        error("There is no residual left after orthogonalization.")
    end
    if any(isinf, u) || any(isnan, u)
        error("Vector contains invalid numbers after orthogonalization")
    end

    push!(h, sqrt(dot(u, u)))
    push!(this.H, h)
    push!(Q, u./h[end])

    if length(Q) > this.maxK
        popfirst!(Q)
    end
    this.itr_count += 1
    return Q[end], this.H[end][end]
end

# Access elements and objects from this algorithms.

function GetHessenberMatrix(this::Arnoldi)
    if this.itr_count == 0
        error("Hessenberg Decomposition is never runned")
    end

    n = length(this.H)
    m = n + 1
    H = Matrix{typeof(this.H[1][1])}(undef, m, n)
    for IdxJ ∈ 1:n
        for IdxI ∈ (IdxJ + 1):-1:1
            OffsetFromDiag = IdxJ + 1 - IdxI
            if OffsetFromDiag < length(this.H[IdxJ])
                H[IdxI, IdxJ] = this.H[IdxJ][end - OffsetFromDiag]
            else
                H[IdxI, IdxJ] = 0
            end
        end
        for IdxI ∈ (IdxJ + 2):m  # zeros below everything.
            H[IdxI, IdxJ] = 0
        end
    end
    return H
end

function GetHMatrix(this::Arnoldi)
    return GetHessenberMatrix(this)
end

function GetOrthogonalMatrix(this::Arnoldi)
    Q = Matrix{typeof(this.Q[1][1])}(undef, length(this.Q[1]), length(this.Q))
    for IdxJ ∈ 1:size(Q, 2)
        Q[:, IdxJ] = this.Q[IdxJ]
    end
    return Q
end

function GetQMatrix(this::Arnoldi)
    return GetOrthogonalMatrix(this)
end


text:
@testset "MPI tests" begin
    testdir = @__DIR__
    run(`$(Base.julia_cmd()) $(joinpath(testdir, "parallel.jl"))`)
    mpiexec(cmd ->run(`$cmd -np 2 $(Base.julia_cmd()) $(joinpath(testdir, "parallel.jl"))`))
    mpiexec(cmd ->run(`$cmd -np 3 $(Base.julia_cmd()) $(joinpath(testdir, "parallel.jl"))`))
    run(`$(Base.julia_cmd()) $(joinpath(testdir, "../examples/farmer_mpi.jl"))`)
    mpiexec(cmd ->run(`$cmd -np 2 $(Base.julia_cmd()) $(joinpath(testdir, "../examples/farmer_mpi.jl"))`))
    mpiexec(cmd ->run(`$cmd -np 3 $(Base.julia_cmd()) $(joinpath(testdir, "../examples/farmer_mpi.jl"))`))
end


text:
#Test start to finish

using SimpleNeuralNet
using LinearAlgebra #for norm()
using Printf
using StatsBase
using Statistics

#Define the training parameters
lr=0.1;      #learning rate
epochs=10;   #epochs
bs=100;      #batch size

#Define some important functions
SimpleNeuralNet.nnAct(z) = sigmoid.(z); #activation function (assume the same everywhere)
SimpleNeuralNet.nnDAct(z) = dsigmoid.(z); #derivative of activation function (assume the same everywhere)
SimpleNeuralNet.nnCost(datay,y) = 0.5*norm(y-datay,2)^2;
SimpleNeuralNet.nnDCost(datay,y) = y - datay;

#Get the mnist dataset
using MLDatasets
train_x, train_y = MNIST.traindata();
datax = Float64.(MNIST.convert2features(train_x));
#datay = Float64.(train_y);
datay = zeros(10,length(train_y));
for i=0:9; datay[i+1,:] = (train_y .== i); end

#Build the neural net
nIn  = size(datax,1);
nOut = size(datay,1);
nNeurons = [nIn; 300; nOut]; #number of neurons by layer
nn = neuralNet(nNeurons);

#Initialize the weights and biases to random values
nn.w = rand(length(nn.w));
nn.b = rand(length(nn.b));
#Normalize
nn.w ./= sum(nn.w);

#Compute performance
outputs = nnComputeOutputs(nn,datax);
nnComputeCosts(outputs,datay;print=true);
bestGuess=zeros(Int64,size(datax,2)); for i=1:size(datax,2); bestGuess[i] = findmax(outputs[:,i])[2]-1; end
@printf("Error rate: %8.6f\n",mean( bestGuess .!= train_y ));

#Train
nnTrain(nn,datax,datay,epochs,lr;batchSize=bs, verbose=true);

#Compute performance
outputs = nnComputeOutputs(nn,datax);
nnComputeCosts(outputs,datay;print=true);
bestGuess=zeros(Int64,size(datax,2)); for i=1:size(datax,2); bestGuess[i] = findmax(outputs[:,i])[2]-1; end
@printf("Error rate: %8.6f\n",mean( bestGuess .!= train_y ));

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet


text:
# 0.18 deprecations. Remove before releasing 0.19
@deprecate path(t::Trie, str::AbstractString) partial_path(t::Trie, str::AbstractString)
@deprecate find_root find_root!
@deprecate top first
@deprecate reverse_iter Iterators.reverse
# Deprecations from #700
Base.@deprecate_binding DisjointSets DisjointSet
Base.@deprecate_binding IntDisjointSets IntDisjointSet
@deprecate DisjointSets(xs...) DisjointSet(xs)
# Enqueue and dequeue deprecations
@deprecate enqueue!(q::Queue, x) Base.push!(q, x)
@deprecate dequeue!(q::Queue) Base.popfirst!(q)


text:
struct NoAcceleration end

acceleration_style(::Type{<:NoAcceleration}) = NoAccelerationStyle()

initialize(::NoAcceleration, ::Any) = nothing


text:
using LSODA, Printf

function fex(t::T1, y::T2, ydot::T3, data::T4) where {T1, T2, T3, T4}
	x = unsafe_wrap(Array,y,neq)
	xdot = unsafe_wrap(Array,ydot,neq)

	xdot[1]=1.0E4 * x[2] * x[3] - .04E0 * x[1]
	xdot[3]=3.0E7 * x[2] * x[2]
	xdot[2]= -xdot[1] - xdot[3]

	return Int32(0)
end

fex_c = @cfunction(fex,Cint,(Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Nothing}))

const atol = Array{Float64}(undef,3)
const rtol = Array{Float64}(undef,3)
t = Array{Float64}(undef,1)
tout = Array{Float64}(undef,1)
const y = Array{Float64}(undef,3)

const neq = Int32(3)

y[1] = 1.0E0
  y[2] = 0.0E0
  y[3] = 0.0E0
  t[1] = 0.0E0
  tout[1] = 0.4E0

opt = lsoda_opt_t()
  opt.ixpr = 0
  opt.rtol = pointer(rtol)
  opt.atol = pointer(atol)
  opt.itask = 1

rtol[1] = 1.0E-4
	rtol[2] = 1.0E-4
  rtol[3] = 1.0E-4
	atol[1] = 1.0E-6
	atol[2] = 1.0E-10
	atol[3] = 1.0E-6

ctx = lsoda_context_t()
  ctx.function_ = fex_c
  ctx.neq = neq
  ctx.state = 1

lsoda_prepare(ctx,opt)

@time for i=1:12
  lsoda(ctx,y,t,tout[1])
  Printf.@printf("at t = %12.4e y= %14.6e %14.6e %14.6e\n",t[1],y[1], y[2], y[3])
  if (ctx.state <= 0)
			error("error istate = ", ctx.state)
	end
	tout[1] *= 10.0E0
end
println("Done!")
lsoda_free(ctx)


text:
"""
# Usage

    solve(prob::FDDEMatrixProblem, h, MatrixForm())

### Reference

https://github.com/mandresik/system-of-linear-fractional-differential-delayed-equations
"""
struct MatrixForm <: FractionalDiffEqAlgorithm end

function solve(prob::FDDEMatrixProblem, h, ::MatrixForm)
    @unpack α, τ, A, B, f, x0, tspan = prob
    t0 = tspan[1]; T = tspan[2]
    limit = 100
    var_num = length(A[:, 1])
    m = ceil(Int, α)

    t_tau = t0-τ
    temp = collect(t0-h:-h:t_tau)
    t = [temp[end:-1:1]; collect(t0:h:T)]

    index = length(temp)
    N = length(t)-index-1

    x = zeros(length(t), var_num)
    b = zeros(1, N)
    F = zeros(var_num, N)
    x0TP = x0(t0)

    @fastmath @inbounds @simd for i=1:index+1
        x[i, :] = x0(t[i])*[1; zeros(m-1, 1)]
    end

    function function_values(M, limit)
        cols::Int64 = size(M, 2)
        risk_index = zeros(1, var_num)
        col = 0
        M0 = copy(M)

        @fastmath @inbounds @simd for r=1:var_num
            @fastmath @inbounds @simd for s=1:cols
                if M0[r, s]==Inf || M0[r, s]==-Inf
                    risk_index = hcat(risk_index[:, 1:col], [r; s])
                    col += 1
                end
            end
        end
        Mt = zeros(var_num, cols*length(t))

        for j=1:length(t)
            Mt[:, (cols*j-cols+1):(cols*j)] = M
        end

        if risk_index !== zeros(var_num, 1)
            risk_start_index = judgeeqsum(t, -1)
            inv_num = judgesum(t, 1)-judgeeqsum(t, -1)
            risks_num = size(risk_index, 2)
            @fastmath @inbounds @simd for j=1:inv_num
                @fastmath @inbounds @simd for k=1:risks_num
                    id = risk_start_index*cols + risk_index[2, k] + cols*(j-1)
                    if Mt[risk_index[1, k], id] > limit
                        Mt[risk_index[1, k], id] = limit
                    elseif Mt[risk_index[1, k], id] < -limit
                        M[risk_index[1, k], id] = -limit
                    end
                end
            end
        end
        return Mt
    end

    At = function_values(A, limit)
    Bt = function_values(B, limit)
    ft = function_values(f, limit)

    @fastmath @inbounds @simd for n=1:N
        col_num = size(A, 2)
        ind = col_num*(n+index)
        b[n] = (n+1)^α-n^α
        F[:, n] = At[:, ind-col_num+1:ind]*x[n+index, :] + Bt[:, ind-col_num+1:ind]*x[n, :]+ft[:, n+index]
        TP = zeros(var_num, 1)
        @fastmath @inbounds @simd for k=0:m-1
            TP += ((t[n+index+1]-t0)^k).*x0TP[:, k+1]./factorial(k)
        end
        sum = zeros(var_num, 1)
        @fastmath @inbounds @simd for j=1:n
            sum += b[n-j+1].*F[:, j]
        end
        x[n+index+1, :] = TP + h^α.*sum
    end
    return x
end

function judgesum(t, thres)
    result = 0
    for i in t
        i < thres ? result+=i : continue
    end
    return result
end

function judgeeqsum(t, thres)
    result = 0
    for i in t
        i <= thres ? result+=i : continue
    end
    return result
end


text:
mutable struct OdeResult
    t::Stack{Float64}
    q::MatrixStack{Float64}

    function OdeResult(x::Vector, start_size::Int = 16)
        L = length(x)
        t = datastructs.Stack{Float64}(start_size);
        q = datastructs.MatrixStack{Float64}(L, start_size);

        return new(t, q)
    end
end

function n_dim(ode_res::OdeResult)
    return size(ode_res.q.data, 1)
end

function Base.length(ode_res::OdeResult)
    return Base.length(ode_res.t)
end

function init_ode_result(x)
    return OdeResult(x)
end

function Base.push!(ode_res::OdeResult, t, q)
    push!(ode_res.t, t)
    push!(ode_res.q, q)
end

function time(ode_res::OdeResult)
    v = view(ode_res.t.data, 1:ode_res.t.ptr)
    return v
end

function variables(ode_res::OdeResult)
    v = view(ode_res.q.data, 1:size(ode_res.q.data, 1), 1:ode_res.t.ptr)
    return v
end

function shrink_to_fit!(ode_res::OdeResult)
    datastructs.shrink_to_fit!(ode_res.t)
    datastructs.shrink_to_fit!(ode_res.q)
    return ode_res
end

struct VariableIndex
    name::Symbol
    index::Int
end

function plot(ode_res::OdeResult; max_var = 5)

    n_vars = n_dim(ode_res)
    n_max_var = min(max_var, n_vars)

    figh = PyPlot.figure()
    ax1 = PyPlot.subplot(n_max_var, 1, 1)
    PyPlot.plot( ode_res.t.data, ode_res.q.data[1, :] )
    PyPlot.grid()

    for ii = 2:n_max_var
        PyPlot.subplot(n_max_var, 1, ii, sharex = ax1)
        PyPlot.plot( ode_res.t.data, ode_res.q.data[ii, :] )
        PyPlot.grid()
    end

    return figh

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# using Revise,
using Test, ForwardDiff, Parameters, Setfield, LinearAlgebra
# using Plots
using BifurcationKit, Test
const BK = BifurcationKit

norminf = x -> norm(x, Inf)
####################################################################################################
function COm(u, p)
	@unpack q1,q2,q3,q4,q5,q6,k = p
	x, y, s = u
	z = 1-x-y-s
	[
		2q1 * z^2 - 2q5 * x^2 - q3 * x * y,
		q2 * z - q6 * y - q3 * x * y,
		q4 * (z - k * s)
	]
end
jet = BK.getJet(COm; matrixfree=false)

par_com = (q1 = 2.5, q2 = 1., q3 = 10., q4 = 0.0675, q5 = 1., q6 = 0.1, k = 0.4)

z0 = [0.001137, 0.891483, 0.062345]

opts_br = ContinuationPar(pMin = 0.5, pMax = 2.3, ds = 0.002, dsmax = 0.01, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 25, nev = 3, maxSteps = 20000)

br, = continuation(jet[1], jet[2], z0, par_com, (@lens _.q2), opts_br;
	recordFromSolution = (x, p) -> (x = x[1], y = x[2], s = x[3]),
	normC = norminf, bothside = true)
####################################################################################################
@set! opts_br.newtonOptions = NewtonPar(maxIter = 10, tol = 1e-12)
sn_codim2, = continuation(jet[1:2]..., br, 2, (@lens _.k), ContinuationPar(opts_br, pMax = 3.2, pMin = 0., detectBifurcation = 1, dsmin=1e-5, ds = -0.001, dsmax = 0.05, nInversion = 6) ;
	normC = norminf,
	detectCodim2Bifurcation = 2,
	updateMinAugEveryStep = 1,
	startWithEigen = true,
	d2F = jet[3], d3F = jet[4],
	recordFromSolution = (u,p; kw...) -> (x = BK.getVec(u)[1] ),
	bothside=true,
	bdlinsolver = MatrixBLS()
	)

@test sn_codim2.specialpoint[1].type == :bt
@test sn_codim2.specialpoint[2].type == :cusp
@test sn_codim2.specialpoint[3].type == :bt
@test sn_codim2.specialpoint[1].printsol.k 	≈ 0.971397 rtol = 1e-4
@test sn_codim2.specialpoint[1].printsol.q2 ≈ 1.417628 rtol = 1e-4
@test sn_codim2.specialpoint[3].printsol.k 	≈ 0.722339 rtol = 1e-4
@test sn_codim2.specialpoint[3].printsol.q2 ≈ 1.161199 rtol = 1e-4

# cusp normal form
cp = computeNormalForm(jet..., sn_codim2, 2; verbose = true)

@test isapprox(cp.nf.c, 0.362; rtol = 1e-2)

# case of BT point
bt = computeNormalForm(jet..., sn_codim2, 1; autodiff = false)
bt = computeNormalForm(jet..., sn_codim2, 1; autodiff = true)

@test isapprox(bt.nf.a |> abs, 0.083784; rtol = 1e-4)
@test isapprox(bt.nf.b |> abs, 2.1363; rtol = 1e-4)
@test isapprox(bt.nfsupp.K2, [-13.113, 51.15]; rtol = 1e-4)
@test isapprox(bt.nfsupp.d, -0.1778; rtol = 1e-3)
@test isapprox(bt.nfsupp.e, -7.1422; rtol = 1e-4)
@test isapprox(bt.nfsupp.a1, abs( -0.8618 ); rtol = 1e-4)
@test isapprox(bt.nfsupp.b1, abs( -7.1176 ); rtol = 1e-4)

hp_codim2, = continuation(jet[1:2]..., br, 1, (@lens _.k), ContinuationPar(opts_br, pMin = 0., pMax = 2.8, detectBifurcation = 1, ds = -0.0001, dsmax = 0.02, dsmin = 1e-4, nInversion = 12, maxSteps = 150, maxBisectionSteps = 35 ) ;
	normC = norminf,
	tangentAlgo = BorderedPred(),
	# linearAlgo = MatrixBLS(),
	detectCodim2Bifurcation = 2,
	updateMinAugEveryStep = 1,
	startWithEigen = true,
	d2F = jet[3], d3F = jet[4],
	bothside = true,
	bdlinsolver = MatrixBLS())

show(hp_codim2)

# @test hp_codim2.specialpoint[1].type == :bt
@test hp_codim2.specialpoint[2].type == :gh
@test hp_codim2.specialpoint[3].type == :gh
if length(hp_codim2.specialpoint) == 4
	@test hp_codim2.specialpoint[4].type == :bt
end

@test hp_codim2.specialpoint[2].printsol.k 	≈ 0.305879 rtol = 1e-3
@test hp_codim2.specialpoint[2].printsol.q2 ≈ 0.924255 rtol = 1e-3
@test hp_codim2.specialpoint[3].printsol.k 	≈ 0.232 rtol = 1e-2
@test hp_codim2.specialpoint[3].printsol.q2 ≈ 0.896099 rtol = 1e-2

gh = computeNormalForm(jet..., hp_codim2, 3)
@test isapprox(gh.nf.G21, -24.9125im; rtol = 1e-4)
@test isapprox(gh.nf.G32, -9322.74711 + 51539.23090im; rtol = 1e-4)
@test isapprox(gh.nf.l2, -776.89; rtol = 1e-4)
####################################################################################################
# branch switching at BT
# hp_from_bt, = continuation(jet..., sn_codim2, 3 , setproperties(sn_codim2.contparams, ds = 0.001, maxSteps = 200);
# 		verbosity = 3, plot  = true,
# 		δp = 1e-5,
# 		normC = norminf,
# 		updateMinAugEveryStep = 1,
# 		bothside = true,
# 		recordFromSolution = (u,p; kw...) -> (x = BK.getVec(u)[1] ),
# 		)
# bt = computeNormalForm(jet..., sn_codim2, 3; verbose = true)
# 	display(bt.x0)
# 	show(bt)
# HC = predictor(bt, Val(:HopfCurve), 0.)
#
# plot(sn_codim2, vars=(:q2, :k))
# 	_S = LinRange(0, 1e-7, 1000)
# 	plot!([HC.hopf(s)[1] for s in _S], [HC.hopf(s)[2] for s in _S], linewidth=3, label = "Hpred")
# 	plot!(hp_codim2, vars=(:q2, :k), color = :black)
#
# HC.ω(0.01)
#
# HC = predictor(bt, Val(:HopfCurve), 0.)
# 	HC.EigenVec(0.01)
#
#
# plot(sn_codim2, vars=(:q2, :x))
# 	plot!(hp_codim2, vars=(:q2, :x))
# 	plot!(hp_from_bt, vars=(:q2, :x), marker = :d)


text:
module Positions

using FinancialInstruments, FixedPointDecimals

export Currencies, Currency
export Countries, Country
export FinancialInstruments, FinancialInstrument, Cash
export Position

abstract type AbstractPosition end

struct Position{F<:FinancialInstrument,A<:Real} <: AbstractPosition
    amount::A
end
Position(p,a) = Position{typeof(p),typeof(a)}(a)
Position(c::Cash{C},a) where C = Position{Cash{C},FixedDecimal{Int,Currencies.unit(C())}}(FixedDecimal{Int,Currencies.unit(C())}(a))

Base.promote_rule(::Type{Position{F,A1}}, ::Type{Position{F,A2}}) where {F,A1,A2} = Position{F,promote_type(A1,A2)}
Base.convert(::Type{Position{F,A}}, x::Position{F,A}) where {F,A} = x
Base.convert(::Type{Position{F,A}}, x::Position{F,<:Real}) where {F,A} = Position{F,A}(convert(A,x.amount))

Base.:+(p1::Position{F,A},p2::Position{F,A}) where {F,A} = Position{F,A}(p1.amount+p2.amount)
Base.:+(p1::AbstractPosition,p2::AbstractPosition) = +(promote(p1,p2)...)

Base.:-(p1::Position{F,A},p2::Position{F,A}) where {F,A} = Position{F,A}(p1.amount-p2.amount)
Base.:-(p1::AbstractPosition,p2::AbstractPosition) = -(promote(p1,p2)...)

Base.:/(p1::Position{F,A},p2::Position{F,A}) where {F,A} = p1.amount/p2.amount
Base.:/(p1::AbstractPosition,p2::AbstractPosition) = /(promote(p1,p2)...)

# TODO: Should scalar multiplication and division respect the Position type or do normal promotion?

# Base.:/(p::Position{F,A},k::R) where {F,A,R<:Real} = Position{F,promote_type(A,R)}(p.amount/k)
Base.:/(p::Position{F,A},k::R) where {F,A,R<:Real} = Position{F,A}(p.amount/k)

# Base.:*(k::R,p::Position{F,A}) where {F,A,R<:Real} = Position{F,promote_type(A,R)}(p.amount*k)
Base.:*(k::R,p::Position{F,A}) where {F,A,R<:Real} = Position{F,A}(p.amount*k)
Base.:*(p::Position{F,A},k::R) where {F,A,R<:Real} = k*p

Base.show(io::IO,c::Position{Cash{C}}) where C = print(io,c.amount," ",C())
Base.show(io::IO,::MIME"text/plain",c::Position{Cash{C}}) where C = print(io,c.amount," ",C())

Base.zero(::Type{Position{F,A}}) where {F,A} = Position{F,A}(zero(A))
Base.one(::Type{Position{F,A}}) where {F,A} = Position{F,A}(one(A))

for (sym,ccy) in Currencies.list
    @eval Positions begin
        $sym = Position(Cash($ccy),1)
    end
end

end # module


text:
module IndexTest

using HPAT
using ParallelAccelerator
using CompilerTools
using HDF5
using MPI
using Base.Test
using HPAT.Partitioning
using HPAT.SEQ
using HPAT.TWO_D
using HPAT.ONE_D

#CompilerTools.OptFramework.set_debug_level(3)
#CompilerTools.LivenessAnalysis.set_debug_level(5)
#CompilerTools.LambdaHandling.set_debug_level(3)
#CompilerTools.TransitiveDependence.set_debug_level(3)
#HPAT.DistributedPass.set_debug_level(3)
#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.ParallelIR.set_debug_level(3)

@acc hpat function index_test(numCenter, file_name)
    points = DataSource(Matrix{Float64},HDF5,"/points", file_name)
    D,N = size(points) # number of features, instances
    centroids :: Array{Float64,2} = rand(D, numCenter)

#    dist::Array{Array{Float64,1},1} = [ Float64[sum(points[:,i+j]) for j in 1:numCenter] for i in 1:(N-numCenter)]
#    dist::Array{Array{Float64,1},1} = [ Float64[sqrt(sum((points[:,i+j]-centroids[:,j]))) for j in 1:numCenter] for i in 1:(N-numCenter)]
    dist::Array{Array{Float64,1},1} = [ Float64[sqrt(sum((points[:,i+j]-centroids[:,j]).^2)) for j in 1:numCenter] for i in 1:(N-numCenter)]
    labels::Array{Int,1} = [indmin(dist[i]) for i in 1:(N-numCenter)]
    return sum(labels)
end

function main()
    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        h5write("test1_1.hdf5", "/points", [1 2 3 4; 5 6 7 8])
    end
    a = index_test(2, "test1_1.hdf5")
    println(HPAT.get_saved_array_partitionings())

    # The partitioning for points should be SEQ.
    @test HPAT.get_saved_partitioning_for_symbol(:points)==SEQ
#    for p in values(HPAT.get_saved_array_partitionings())
#        @test p==SEQ
#    end

    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        rm("test1_1.hdf5")
    end

end

end

IndexTest.main()


text:
for (op,cfunc) in ((:factorial,:arb_fac_ui), (:doublefactorial,:arb_doublefac_ui))
  @eval begin
    function ($op)(x::ArbFloat{P}) where {P}
      signbit(x) && throw(ErrorException("Domain Error: argument is negative"))
      y = trunc(UInt, x)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, UInt, Int), z, y, P)
      z
    end
  end
end

function doublefactorial(xx::R) where {R <: Real}
   P = precision(ArbFloat)
   x = convert(ArbFloat{P},xx)
   doublefactorial(x)
end

for (op,cfunc) in ((:risingfactorial,:arb_rising),)
  @eval begin
    function ($op)(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
      signbit(x) && throw(ErrorException("Domain Error: argument is negative"))
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      z
    end
    function ($op)(xx::R, y::ArbFloat{P}, prec::Int=P) where {R <: Real,P}
      x = convert(ArbFloat{P},xx)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      z
    end
    function ($op)(x::ArbFloat{P}, yy::R, prec::Int=P) where {R <: Real,P}
      y = convert(ArbFloat{P},yy)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      z
    end
    #=
    function ($op){R1<:Real,R2<:Real}(xx::R1, yy::R2)
      P = precision(ArbFloat)
      x = convert(ArbFloat{P},xx)
      y = convert(ArbFloat{P},yy)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Void, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      z
    end
    =#
  end
end

for (op,cfunc) in ((:agm, :arb_agm), (:polylog, :arb_polylog))
  @eval begin
    function ($op)(x::ArbFloat{P}, y::ArbFloat{P}, prec::Int=P) where {P}
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      z
    end
    function ($op)(xx::R, y::ArbFloat{P}, prec::Int=P) where {R <: Real,P}
      x = convert(ArbFloat{P},xx)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      z
    end
    function ($op)(x::ArbFloat{P}, yy::R, prec::Int=P) where {R <: Real,P}
      y = convert(ArbFloat{P},yy)
      z = initializer(ArbFloat{P})
     # ccall(@libarb($cfunc), Void, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      z
    end
    #=
    function ($op){R1<:Real,R2<:Real}(xx::R1, yy::R2)
      P = precision(ArbFloat)
      x = convert(ArbFloat{P},xx)
      y = convert(ArbFloat{P},yy)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Void, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      z
    end
    =#
  end
end


text:
@block AnthonyWang [tutorials] begin

    @cell "Tutorial simple scatter" [tutorial, scatter, "2d"] begin
        x = rand(10)
        y = rand(10)
        colors = rand(10)
        scene = scatter(x, y, color = colors)
    end

    @cell "Tutorial markersize" [tutorial, scatter, markersize, "2d"] begin
        x = 1:10
        y = 1:10
        sizevec = x ./ 10
        scene = scatter(x, y, markersize = sizevec)
    end

    @cell "Tutorial simple line" [tutorial, line, "2d"] begin
        x = range(0, stop = 2pi, length = 40)
        f(x) = sin.(x)
        y = f(x)
        scene = lines(x, y, color = :blue)
    end

    @cell "Tutorial adding to a scene" [tutorial, line, scene, markers, "2d"] begin
        x = range(0, stop = 2pi, length = 80)
        f1(x) = sin.(x)
        f2(x) = exp.(-x) .* cos.(2pi*x)
        y1 = f1(x)
        y2 = f2(x)

        scene = lines(x, y1, color = :blue)
        scatter!(scene, x, y1, color = :red, markersize = 0.1)

        lines!(scene, x, y2, color = :black)
        scatter!(scene, x, y2, color = :green, marker = :utriangle, markersize = 0.1)
    end

    @cell "Tutorial removing from a scene" [tutorial, line, scene, markers, "2d"] begin
        x = range(0, stop = 2pi, length = 80)
        y1 = sin.(x)
        y2 = exp.(-x) .* cos.(2pi * x)

        scene = lines(x, y1, color = :blue)
        scatter!(scene, x, y1, color = :red, markersize = 0.1)

        lines!(scene, x, y2, color = :black)
        scatter!(scene, x, y2, color = :green, marker = :utriangle, markersize = 0.1)
        # initialize the stepper and give it an output destination
        st = Stepper(scene, @replace_with_a_path)
        step!(st)

        # remove last plot in scene
        delete!(scene, scene[end])
        step!(st)
        # remove second plot in scene
        delete!(scene, scene[2])
        step!(st)
        st
    end

    @cell "Tutorial adjusting scene limits" [tutorial, scene, limits, "2d"] begin
        x = range(0, stop = 10, length = 40)
        y = x
        #= specify the scene limits, note that the arguments for FRect are
            x_min, y_min, x_dist, y_dist,
            therefore, the maximum x and y limits are then x_min + x_dist and y_min + y_dist
        =#
        limits = FRect(-5, -10, 20, 30)

        scene = lines(x, y, color = :blue, limits = limits)

    end

    @cell "Tutorial basic theming" [tutorial, scene, limits, "2d"] begin
        x = range(0, stop = 2pi, length = 40)
        y = cos.(x)
        scene = lines(x, y, color = :blue)

        axis = scene[Axis] # get the axis object from the scene
        axis.grid.linecolor = ((:red, 0.5), (:blue, 0.5))
        axis.names.textcolor = ((:red, 1.0), (:blue, 1.0))
        axis.names.axisnames = ("x", "y = cos(x)")
        scene
    end

    @cell "Tutorial heatmap" [tutorial, heatmap, "2d"] begin
        data = rand(50, 50)
        scene = heatmap(data)
    end

    @cell "Tutorial linesegments" [tutorial, linesegments, "2d"] begin
        points = [
            Point2f0(0, 0) => Point2f0(5, 5);
            Point2f0(15, 15) => Point2f0(25, 25);
            Point2f0(0, 15) => Point2f0(35, 5);
        ]
        linesegments(points, color = :red, linewidth = 2)
    end

    @cell "Tutorial barplot" [tutorial, barplot, "2d"] begin
        data = sort(randn(100))
        barplot(data)
    end

    @cell "Tutorial title" [tutorial, title, "2d"] begin
        scene = lines(rand(10), axis=(names=(title="Random lines",),))
    end

    @cell "Tutorial title function" [tutorial, title, "2d"] begin
        scene = lines(rand(10))
        sc_t = title(scene, "Random lines")
        sc_t
    end

    @cell "Tutorial plot transformation" [tutorial, transformation, "2d"] begin
        data = rand(10)

        scene = Scene()
        st = Stepper(scene, @replace_with_a_path)
        lineplot = lines!(scene, data)[end]    # gets the last defined plot for the Scene
        scatplot = scatter!(scene, data)[end]  # same thing but the last defined plot is scatter
        step!(st)

        rotate!(lineplot, 0.025π) # only the lines are rotated, not the scatter
        step!(st)
        st
    end

    # @cell "Tutorial markersizes" [tutorial, markersize, scatter] begin
    #
    #     scene = scatter(rand(10); markersize = 10px)
    #
    # end

end


text:
#export report, ReportConstraints

#TODO add ability to save report every X iterations
#TODO read report file to warm-start solver
#TODO report file can be configured with cache, inner vars, points etc.

"""
    report_finish(p::DSProblem)

**Not Implemented**

If not silent, give a short printout of the result of the problem.

States the duration, final point, final cost, and reason for stopping.
A more detailed report can be shown with `report`.
"""
function report_finish(p::DSProblem)

end

"""
    report(p::DSProblem; save::Union{Bool,String}=false)

Prints a summary of the problem to the console, or saves it as a text file.

`save=true` saves the report in the current REPL directory, but will error
if the default file name `report` is taken. Otherwise give a path as an argument
to use that name.

The report contains the solver configuration, the initial and final incumbent points,
the initial and final cost, the duration, the number of iterations, and the current
status.
"""
function report(p::DSProblem; save::Union{Bool,String}=false,
               print_problem=true, print_solver=true, print_config=true)
    problem =
"""
---------------------------------
Problem Status
---------------------------------
Incumbent Point:    $(p.x)
Current Cost:       $(p.x_cost)

Initial Point:      $(p.x_initial)
Initial Cost:       $(p.x_initial_cost)

Variables:          $(p.N)
Variable Scaling:   $(p.meshscale)

"""
    solver =
"""
---------------------------------
Solver Status
---------------------------------
Status:             $(p.status)
Iterations:         $(p.iteration)

Duration:           To be added

"""

    config =
"""
---------------------------------
Solver Configuration
---------------------------------
Sense:              $(p.sense)

Mesh:               $(typeof(p.mesh))
Search Directions:  To be added
Poll Directions:    $(typeof(p.poll))

"""
    report_str = "Direct Search Optimisation Report\n=================================\n"

    if print_problem
        report_str *= problem
    end
    if print_solver
        report_str *= solver
    end
    if print_config
        report_str *= config
    end

    report_str *= "================================="

    println(report_str)
end

"""
    export_points(p::DSProblem)

Gives the trace of considered points during optimisation. Note that this is unavailable
if tracking of points has been disabled
"""
function export_points(p::DSProblem)
    for x in p.cache.order
        println("$(x)\t$(p.cache.costs[x])")
    end
end

ReportConstraints(p::DSProblem) = ReportConstraints(p.constraints)
function ReportConstraints(c::Constraints{T}) where T
    constraint_report=
"""
---------------------------------
Constraints
---------------------------------
Constraint Collections:     $(c.count)
Populated Collections:      $(length([col for col in c.collections if col.count > 0]))
Total Constraints:          $(sum([col.count for col in c.collections]))
"""
    print(constraint_report)
    for (index, collection) in enumerate(c.collections)
        println("\n$(CollectionIndex(index)) : $(typeof(collection))")
        ReportConstraintCollection(collection)
    end
end

function ReportConstraintCollection(c::ConstraintCollection{T, C}) where {T,C}
    collection_report =
"""
Constraints:       $(c.count)
Ignored:           $(c.ignore)
h_max:             $(c.h_max)
h_max update
function handle:   $(c.h_max_update)
result aggregate
function handle:   $(c.result_aggregate)
"""
    print(collection_report)
end


text:
#=using Pkg
import Pkg; Pkg.add("MathOptInterface")
import Pkg; Pkg.add("Statistics")
import Pkg; Pkg.add("QuantEcon")
import Pkg; Pkg.add("Interpolations")
import Pkg; Pkg.add("NLsolve")
import Pkg; Pkg.add("Optim")
import Pkg; Pkg.add("Random")
import Pkg; Pkg.add("IterTools")
import Pkg; Pkg.add("JuMP")
import Pkg; Pkg.add("Ipopt")
import Pkg; Pkg.add("BenchmarkTools")
import Pkg; Pkg.add("GaussianProcesses")
import Pkg; Pkg.add("TimerOutputs")
=#

using
    MathOptInterface,
    BenchmarkTools,
    LinearAlgebra,
    Statistics,
    QuantEcon,
    Interpolations,
    NLsolve,
    Optim,
    Random,
    IterTools,
    JuMP,
    Ipopt,
    GaussianProcesses,
    Random,
    TimerOutputs;
    #DateTime;

δ = 0.1
β = 0.86
ϕ = 0.6
ξ = 0.5
μ = 0.5
ψ = 1/ϕ
b = 1/(ϕ^ϕ * (1-ϕ)^(1-ϕ))

numSectors = 2
numLin = 10
gridMin = 2
gridMax = 30

u(x,y) = (log(x) + log(y))
wInit(x,y) = (log(x) + log(y))/(1 - β)
function gradwInit(g, x, y)
    g[1] = 1/x/(1-β)
    g[2] = 1/y/(1-β)
end
function hesswInit(h, x, y)
    h[1,1] = -1/x^2/(1-β)
    h[1,2] = 0
    h[2,1] = 0
    h[2,2] = -1/y^2/(1-β)
end
yTargGrd = Vector{Float64}[]
sTargGrd = Vector{Float64}[]
mTargGrd = Matrix{Float64}[]
xTargGrd = Matrix{Float64}[]
wTargGrd = Float64[]
grid = Vector{Float64}[]
for p in product(LinRange(gridMin,gridMax,numLin), LinRange(gridMin, gridMax,numLin))
    push!(grid, collect(p))
end
rng = MersenneTwister(1234)
#for timing
tOGrd = TimerOutput()
#wval = w.(grid)

#the following function generates target values (of the value function) for the ML process
##note that (unlike Sargent and Stachursky) the grid is written in terms of kapital
#this simplified the problem substantially as otherwise intermediate variables became dynamic
##also note that this accepts the smooth function
function TargGrd(w, grid, β ; return_policy = false)
    for n in 1:length(grid)
        @timeit tOGrd string("Pointwise-for-", length(grid),"-pts") begin
        k = grid[n]
        modTrial = Model(Ipopt.Optimizer)
        set_silent(modTrial)
        @variable(modTrial, c[1:numSectors] >= 0.0001)
        @variable(modTrial, y[1:numSectors] >= 0.0001)
        @variable(modTrial, kp[1:numSectors]>= 0.0001)
        @variable(modTrial, x[1:numSectors, 1:numSectors]>=0.0001)
        @variable(modTrial, xComb[1:numSectors]>=0.0001)
        @variable(modTrial, m[1:numSectors, 1:numSectors]>=0.0001)
        @variable(modTrial, mComb[1:numSectors]>=0.0001)
        @variable(modTrial, u)
        @variable(modTrial, w)
        for i in 1:numSectors
            #output
            @NLconstraint(modTrial, y[i] == b * k[i]^ϕ * mComb[i]^(1-ϕ))
            #future capital
            @constraint(modTrial, kp[i] == xComb[i] + (1-δ) * k[i])
            #supply equals demand (eq'm constraint)
            @constraint(modTrial, 0 == y[i] - c[i] - sum(m[i,:]) - sum(x[i,:])) #done
        end
        @NLconstraint(modTrial, xComb[1] == x[1,1]^ξ * x[2,1]^(1-ξ))
        @NLconstraint(modTrial, xComb[2] == x[1,2]^ξ * x[2,2]^(1-ξ))
        @NLconstraint(modTrial, mComb[1] == m[1,1]^μ * m[2,1]^(1-μ))
        @NLconstraint(modTrial, mComb[2] == m[1,2]^μ * m[2,2]^(1-μ))
        @NLconstraint(modTrial, u == log(c[1]) + log(c[2]))
        @NLconstraint(modTrial, w == (log(kp[1]) + log(kp[2]))/(1-β))
#        register(modTrial, :wInit, 2, wInit, gradwInit)
#        @NLconstraint(modTrial, w == wInit(kp[1],kp[2]))
#        register(modTrial, :u, 2, u; autodiff = true)
        @NLobjective(modTrial, Max, u + β * w)
        JuMP.optimize!(modTrial)
        push!(wTargGrd, JuMP.objective_value(modTrial))
        if return_policy
           push!(xTargGrd, value.(x))
           push!(sTargGrd, value.(c))
           push!(mTargGrd, value.(m))
           push!(yTargGrd, value.(y))
        end
                        end #of timeit
#        global intK[n] = value.(f)
    end
    #global prevK = intK
    return wTargGrd
end

function evalBellGrd()
    wTargGrd  = TargGrd(wInit, grid, β; return_policy= false)
    show(tOGrd)
    #reset containers and grid
    #tOCbd = TimerOutput()
    #yTargGrd = Vector{Float64}[]
    #sTargGrd = Vector{Float64}[]
    #mTargGrd = Matrix{Float64}[]
    #xTargGrd = Matrix{Float64}[]
    #wTargGrd = Float64[]
    #grid = Vector{Float64}[]
    #for p in product(LinRange(gridMin,gridMax,numLin), LinRange(gridMin, gridMax, numLin))
        #push!(grid, collect(p))
    #end
    println("\n ---------")
    return wTargGrd
end
#=
XTargGrd = Vector{Float64}[]
       for j = 1:numLin^numSectors
           push!(XTargGrd, sum(xTargGrd[j][:,i] for i in 1:numSectors))
       end
MTargGrd = Vector{Float64}[]
       for j = 1:numLin^numSectors
           push!(MTargGrd, sum(mTargGrd[j][:,i] for i in 1:numSectors))
       end
fTargGrd = sTargGrd + MTargGrd + XTargGrd

MPCf = Vector{Float64}[]
       for j = 1:numLin^numSectors
           push!(MPCf, sTargGrd[j]./fTargGrd[j])
       end
MPC = Vector{Float64}[]
       for j = 1:numLin^numSectors
           push!(MPC, sTargGrd[j]./yTargGrd[j])
       end
MPM = Vector{Float64}[]
       for j = 1:numLin^numSectors
           push!(MPM, MTargGrd[j]./yTargGrd[j])
       end
MPX = Vector{Float64}[]
       for j = 1:numLin^numSectors
           push!(MPX, XTargGrd[j]./yTargGrd[j])
       end
# a function for turning vectors of vectors into matrices
#function vvm(x)
#           dim1 = length(x)
#           dim2 = length(x[1])
#           matrix = zeros(Float64, dim1, dim2)
#           @inbounds @fastmath for i in 1:dim1, for j in 1:dim2
#                   matrix[i, j] = x[i][j]
#               end
#           end
#           return matrix
#end
#alternatively
gridM = reduce(hcat, grid)

mz = MeanZero()
ml = MeanLin([2., 3.])
mp = MeanPoly([2. 3.; 1. 2.])
kern = SE(log(2), log(2000.0))
gp = GP(gridM, wTarg+0.2*randn(numLin^numSectors), mz, kern)
#gp = GP(gridM, wTargGrd, MeanLin([2., 3.]), kern)
GaussianProcesses.optimize!(gp)
predict_f(gp, gridM)[1]
fstar(x) = predict_f(gp,hcat(x))[1]

plot!(contour(gp), heatmap(gp); fmt=:png)
#wTargp(x, y) = predict_f(gp,hcat([x, y]))[1]
=#

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
"A type for the ensemble average scattering coefficients.
Here they are discretised in terms of the depth x of the halfspace"
mutable struct DiscretePlaneWaveMode{T<:AbstractFloat}
    basis_order::Int # largest hankel order
    x::Vector{T} # spatial mesh
    amplitudes::Array{Complex{T}} # a matrix of the scattering amplitudes, size(A_mat) = (length(x), 2basis_order +1, number_of_species)
    # Enforce that the dimensions are correct
    function DiscretePlaneWaveMode{T}(basis_order::Int, x::Vector{T}, amplitudes::Array{Complex{T}}) where T <: AbstractFloat
        if (length(x), 2*basis_order+1) != size(amplitudes)[1:2]
            error("The amplitudes of DiscretePlaneWaveMode does not satisfy size(amplitudes)[1:2] == (length(X), 2*basis_order+1)")
        end
        new(basis_order,x,amplitudes)
    end
end

DiscretePlaneWaveMode(M::Int, x::AbstractVector{T}, as::AbstractArray{Complex{T}}) where T<:AbstractFloat = DiscretePlaneWaveMode{T}(M,collect(x),collect(as))

function DiscretePlaneWaveMode(x::AbstractVector{T}, A_mat::Array{Complex{T}}) where T<:Number
    DiscretePlaneWaveMode(Int((size(A_mat,2)-1)/2), collect(x), A_mat)
end

"Approximates the error in DiscretePlaneWaveMode.amplitudes due to the mesh DiscretePlaneWaveMode.x."
function discretewave_error(avg_w::DiscretePlaneWaveMode)
    ddf = circshift(avg_w.amplitudes, (2,0,0)) - 2*circshift(avg_w.amplitudes, (1,0,0)) + avg_w.amplitudes
    h = (avg_w.x[2] - avg_w.x[1])
    ddf = ddf[3:end,:,:]/(h^2)
    max_ddf = maximum(abs.(ddf))

    # trapezoidal error
    return max_ddf*h^3/12
end

"Calculates an DiscretePlaneWaveMode from one EffectiveWave"
function DiscretePlaneWaveMode(xs::AbstractVector{T}, wave_eff::EffectivePlaneWaveMode{T}, halfspace::Halfspace{T}) where T<:Number

    amps = wave_eff.eigenvectors
    ho = wave_eff.basis_order

    θ_eff = transmission_angle(wave_eff, halfspace)
    kcos_eff = wave_eff.wavenumber * dot(-conj(halfspace.normal), wave_eff.direction)

    S = size(amps,2)
    P = size(amps,3)
    if P > 1
        @warn "The plane-wave has more than one eigenvector per wavenumber. This case has not been implemented so using only the first. This occured for: ω = $(wave_eff.ω) and wavenumber = $(wave_eff.wavenumber)"
    end

    average_amps = [
        im^T(m)*exp(-im*m*θ_eff) * amps[m+ho+1,s,1] * exp(im*kcos_eff*x)
    for x in xs, m=-ho:ho, s=1:S]

    return DiscretePlaneWaveMode(ho,xs,average_amps)
end

"Calculates an DiscretePlaneWaveMode from a vector of EffectiveWave"
function DiscretePlaneWaveMode(xs::AbstractVector{T}, wave_effs::Vector{EffectivePlaneWaveMode{T,Dim}}, halfspace::Halfspace{T}) where {T<:Number,Dim}

    N = wave_effs[1].basis_order
    if !isempty(findall([w.basis_order != N for w in wave_effs]))
        @error "The effective waves $wave_effs have different orders, so they can not be combined into a discrete wave"
    end

    discrete_waves = [DiscretePlaneWaveMode(xs, wave, halfspace) for wave in wave_effs]
    amps = sum(discrete_waves[i].amplitudes[:,:,:] for i in eachindex(discrete_waves))

    return DiscretePlaneWaveMode(N, xs, amps)
end

"Numerically solve the integral equation governing ensemble average waves. Optionally can use wave_eff to approximate the wave away from the boundary."
function DiscretePlaneWaveMode(ω::T, source::PlaneSource{T,2,1,Acoustic{T,2}}, material::Material{2,Halfspace{T,2}};
        x::AbstractVector{T} = [zero(T)],
        tol::T = T(1e-4),
        wave_effs::Vector{EffectivePlaneWaveMode{T,2}} = EffectivePlaneWaveMode{T,2}[],
        max_size::Int = 1000,
        kws...
    ) where T<:Number

    k = real(ω/source.medium.c)
    specie = material.species[1]

    if x == [zero(T)]
        if isempty(wave_effs)
            wave_effs = WaveModes(ω, source, material; tol=tol, mesh_points=2, kws...)
        end
        # estimate a large coarse non-dimensional mesh based on the lowest attenuating effective wave
        a12 = T(2) * specie.exclusion_distance * outer_radius(specie)
        x = x_mesh(wave_effs[1]; tol = tol,  a12 = a12, max_size=max_size)
    end

    X = x.*k
    (MM_quad,b_mat) = discrete_wave_system(ω, X, source, material; tol = tol, kws...);

    M = Int( (size(b_mat,2) - 1)/2 )
    J = length(collect(X)) - 1

    len = (J + 1) * (2M + 1)
    MM_mat = reshape(MM_quad, (len, len));
    b = reshape(b_mat, (len));

    As = MM_mat\b
    As_mat = reshape(As, (J+1, 2M+1, 1))

    return DiscretePlaneWaveMode(M, collect(X)./k, As_mat)
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
# -*- coding: utf-8 -*-
#
# Calculate the percentile of the gradient image, which
# used as the 'kappa' parameter of the anisotropic diffusion.
#
# Credits:
# [1] KAZE: nldiffusion_functions.cpp / compute_k_percentile()
#
# Aaron LI
# 2015/07/20
#

include("scharr.jl");

function calc_k_percentile(img, percent=0.7, nbins=300)
    rows, cols = size(img)
    # derivatives of the image
    img_gx = scharr(img, 1, 0)
    img_gy = scharr(img, 0, 1)
    img_modg = sqrt(img_gx.^2 + img_gy.^2)
    # histogram
    hmax = maximum(img_modg)
    hist_e, hist_counts = hist(reshape(img_modg, length(img_modg)), nbins)
    hist_cum = cumsum(hist_counts)
    # find the percent of the histogram percentile
    npoints = sum(img_modg .> 0.0)
    nthreshold = npoints * percent
    k = sum(hist_cum .<= nthreshold)
    kperc = (k == length(hist_cum)) ? 0.03 : (hmax * k / nbins)
    return kperc
end

max_stars_repo_name: liweitianux/atoolbox

max_issues_repo_name: liweitianux/atoolbox

max_forks_repo_name: liweitianux/atoolbox


text:
using DiffEqProblemLibrary, Sundials, Base.Test

prob = prob_ode_linear
dt = 1//2^(4)
sol = solve(prob,ARKODE())
@test sol.errors[:l2] < 5e-2

f1 = (du,u,p,t) -> du .= u
f2 = (du,u,p,t) -> du .= u

prob = SplitODEProblem(f1,f2,rand(4,2),(0.0,1.0))
function (::typeof(prob.f))(::Type{Val{:analytic}},u0,p,t)
    exp(2t)*u0
end
sol = solve(prob,ARKODE())
@test sol.errors[:l2] < 1e-2
sol = solve(prob,ARKODE(),reltol=1e-12,abstol=1e-12)
@test sol.errors[:l2] < 1e-6


text:
using FastTransforms, Test

FastTransforms.ft_set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))

@testset "libfasttransforms" begin
    n = 64
    for T in (Float32, Float64)
        c = one(T) ./ (1:n)
        x = collect(-1 .+ 2*(0:n-1)/T(n))
        f = similar(x)
        @test FastTransforms.horner!(c, x, f) == f
        fd = T[sum(c[k]*x^(k-1) for k in 1:length(c)) for x in x]
        @test f ≈ fd
        @test FastTransforms.clenshaw!(c, x, f) == f
        fd = T[sum(c[k]*cos((k-1)*acos(x)) for k in 1:length(c)) for x in x]
        @test f ≈ fd
        A = T[(2k+one(T))/(k+one(T)) for k in 0:length(c)-1]
        B = T[zero(T) for k in 0:length(c)-1]
        C = T[k/(k+one(T)) for k in 0:length(c)]
        phi0 = ones(T, length(x))
        c = cheb2leg(c)
        @test FastTransforms.clenshaw!(c, A, B, C, x, phi0, f) == f
        @test f ≈ fd
    end

    α, β, γ, δ, λ, μ = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6
    function test_1d_plans(p1, p2, x)
        y = p1*x
        z = p2*y
        @test z ≈ x
        y = p1*x
        z = p1'y
        y = transpose(p1)*z
        z = transpose(p1)\y
        y = p1'\z
        z = p1\y
        @test z ≈ x
        y = p2*x
        z = p2'y
        y = transpose(p2)*z
        z = transpose(p2)\y
        y = p2'\z
        z = p2\y
        @test z ≈ x
        P = p1*I
        Q = p2*P
        @test Q ≈ I
        P = p1*I
        Q = p1'P
        P = transpose(p1)*Q
        Q = transpose(p1)\P
        P = p1'\Q
        Q = p1\P
        @test Q ≈ I
        P = p2*I
        Q = p2'P
        P = transpose(p2)*Q
        Q = transpose(p2)\P
        P = p2'\Q
        Q = p2\P
        @test Q ≈ I
    end

    for T in (Float32, Float64, Complex{Float32}, Complex{Float64}, BigFloat, Complex{BigFloat})
        x = T(1)./(1:n)
        Id = Matrix{T}(I, n, n)
        for (p1, p2) in ((plan_leg2cheb(Id), plan_cheb2leg(Id)),
                         (plan_ultra2ultra(Id, λ, μ), plan_ultra2ultra(Id, μ, λ)),
                         (plan_jac2jac(Id, α, β, γ, δ), plan_jac2jac(Id, γ, δ, α, β)),
                         (plan_lag2lag(Id, α, β), plan_lag2lag(Id, β, α)),
                         (plan_jac2ultra(Id, α, β, λ), plan_ultra2jac(Id, λ, α, β)),
                         (plan_jac2cheb(Id, α, β), plan_cheb2jac(Id, α, β)),
                         (plan_ultra2cheb(Id, λ), plan_cheb2ultra(Id, λ)))
            test_1d_plans(p1, p2, x)
        end
    end

    for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
        x = T(1)./(1:n)
        Id = Matrix{T}(I, n, n)
        p = plan_associatedjac2jac(Id, 1, α, β, γ, δ)
        V = p*I
        @test V ≈ p*Id
        y = p*x
        @test V\y ≈ x
    end

    @testset "Modified classical orthonormal polynomial transforms" begin
        (n, α, β) = (16, 0, 0)
        P1 = plan_modifiedjac2jac(Float64, n, α, β, [0.9428090415820636, -0.32659863237109055, -0.42163702135578396, 0.2138089935299396]) # u1(x) = (1-x)^2*(1+x)
        P2 = plan_modifiedjac2jac(Float64, n, α, β, [0.9428090415820636, -0.32659863237109055, -0.42163702135578396, 0.2138089935299396], [1.4142135623730951]) # u2(x) = (1-x)^2*(1+x)
        P3 = plan_modifiedjac2jac(Float64, n, α, β, [-0.9428090415820636, 0.32659863237109055, 0.42163702135578396, -0.2138089935299396], [-5.185449728701348, 0.0, 0.42163702135578374]) # u3(x) = -(1-x)^2*(1+x), v3(x) = -(2-x)*(2+x)
        P4 = plan_modifiedjac2jac(Float64, n, α+2, β+1, [1.1547005383792517], [4.387862045841156, 0.1319657758147716, -0.20865621238292037]) # v4(x) = (2-x)*(2+x)

        @test P1*I ≈ P2*I
        @test P1\I ≈ P2\I
        @test P3*I ≈ P2*(P4*I)
        @test P3\I ≈ P4\(P2\I)

        P5 = plan_modifiedlag2lag(Float64, n, α, [2.0, -4.0, 2.0]) # u5(x) = x^2
        P6 = plan_modifiedlag2lag(Float64, n, α, [2.0, -4.0, 2.0], [1.0]) # u6(x) = x^2
        P7 = plan_modifiedlag2lag(Float64, n, α, [2.0, -4.0, 2.0], [7.0, -7.0, 2.0]) # u7(x) = x^2, v7(x) = (1+x)*(2+x)
        P8 = plan_modifiedlag2lag(Float64, n, α+2, [sqrt(2.0)], [sqrt(1058.0), -sqrt(726.0), sqrt(48.0)]) # v8(x) = (1+x)*(2+x)

        @test P5*I ≈ P6*I
        @test P5\I ≈ P6\I
        @test P7*I ≈ P6*(P8*I)
        @test P7\I ≈ P8\(P6\I)

        P9 = plan_modifiedherm2herm(Float64, n, [2.995504568550877, 0.0, 3.7655850551068593, 0.0, 1.6305461589167827], [2.995504568550877, 0.0, 3.7655850551068593, 0.0, 1.6305461589167827]) # u9(x) = 1+x^2+x^4, v9(x) = 1+x^2+x^4

        @test P9*I ≈ P9\I
    end

    function test_nd_plans(p, ps, pa, A)
        B = copy(A)
        C = ps*(p*A)
        A = p\(pa*C)
        @test A ≈ B
        C = ps'*(p'A)
        A = p'\(pa'C)
        @test A ≈ B
        C = transpose(ps)*(transpose(p)*A)
        A = transpose(p)\(transpose(pa)*C)
        @test A ≈ B
    end

    A = sphones(Float64, n, 2n-1)
    p = plan_sph2fourier(A)
    ps = plan_sph_synthesis(A)
    pa = plan_sph_analysis(A)
    test_nd_plans(p, ps, pa, A)
    A = sphones(Float64, n, 2n-1) + im*sphones(Float64, n, 2n-1)
    p = plan_sph2fourier(A)
    ps = plan_sph_synthesis(A)
    pa = plan_sph_analysis(A)
    test_nd_plans(p, ps, pa, A)

    A = sphvones(Float64, n, 2n-1)
    p = plan_sphv2fourier(A)
    ps = plan_sphv_synthesis(A)
    pa = plan_sphv_analysis(A)
    test_nd_plans(p, ps, pa, A)
    A = sphvones(Float64, n, 2n-1) + im*sphvones(Float64, n, 2n-1)
    p = plan_sphv2fourier(A)
    ps = plan_sphv_synthesis(A)
    pa = plan_sphv_analysis(A)
    test_nd_plans(p, ps, pa, A)

    A = diskones(Float64, n, 4n-3)
    p = plan_disk2cxf(A, α, β)
    ps = plan_disk_synthesis(A)
    pa = plan_disk_analysis(A)
    test_nd_plans(p, ps, pa, A)
    A = diskones(Float64, n, 4n-3) + im*diskones(Float64, n, 4n-3)
    p = plan_disk2cxf(A, α, β)
    ps = plan_disk_synthesis(A)
    pa = plan_disk_analysis(A)
    test_nd_plans(p, ps, pa, A)

    A = rectdiskones(Float64, n, n)
    p = plan_rectdisk2cheb(A, β)
    ps = plan_rectdisk_synthesis(A)
    pa = plan_rectdisk_analysis(A)
    test_nd_plans(p, ps, pa, A)
    A = rectdiskones(Float64, n, n) + im*rectdiskones(Float64, n, n)
    p = plan_rectdisk2cheb(A, β)
    ps = plan_rectdisk_synthesis(A)
    pa = plan_rectdisk_analysis(A)
    test_nd_plans(p, ps, pa, A)

    A = triones(Float64, n, n)
    p = plan_tri2cheb(A, α, β, γ)
    ps = plan_tri_synthesis(A)
    pa = plan_tri_analysis(A)
    test_nd_plans(p, ps, pa, A)
    A = triones(Float64, n, n) + im*triones(Float64, n, n)
    p = plan_tri2cheb(A, α, β, γ)
    ps = plan_tri_synthesis(A)
    pa = plan_tri_analysis(A)
    test_nd_plans(p, ps, pa, A)

    α, β, γ, δ = -0.1, -0.2, -0.3, -0.4
    A = tetones(Float64, n, n, n)
    p = plan_tet2cheb(A, α, β, γ, δ)
    ps = plan_tet_synthesis(A)
    pa = plan_tet_analysis(A)
    test_nd_plans(p, ps, pa, A)
    A = tetones(Float64, n, n, n) + im*tetones(Float64, n, n, n)
    p = plan_tet2cheb(A, α, β, γ, δ)
    ps = plan_tet_synthesis(A)
    pa = plan_tet_analysis(A)
    test_nd_plans(p, ps, pa, A)

    A = spinsphones(Complex{Float64}, n, 2n-1, 2) + im*spinsphones(Complex{Float64}, n, 2n-1, 2)
    p = plan_spinsph2fourier(A, 2)
    ps = plan_spinsph_synthesis(A, 2)
    pa = plan_spinsph_analysis(A, 2)
    test_nd_plans(p, ps, pa, A)
end


text:
DP = SymbolicUtils.DynamicPolynomials
Bijections = SymbolicUtils.Bijections

# extracting underlying polynomial and coefficient type from Polyforms
underlyingpoly(x::Number)    = x
underlyingpoly(pf::PolyForm) = pf.p
coefftype(x::Number)    = typeof(x)
coefftype(pf::PolyForm) = DP.coefficienttype(underlyingpoly(pf))

#=
Converts an array of symbolic polynomials
into an array of DynamicPolynomials.Polynomials
=#
function symbol_to_poly(sympolys::AbstractArray)
    @assert !isempty(sympolys) "Empty input."

    # standardize input
    stdsympolys = map(unwrap, sympolys)
    sort!(stdsympolys, lt=(<ₑ))

    pvar2sym  = Bijections.Bijection{Any,Any}()
    sym2term  = Dict{Sym,Any}()
    polyforms = map(f -> PolyForm(f, pvar2sym, sym2term), stdsympolys)

    # Discover common coefficient type
    commontype = mapreduce(coefftype, promote_type, polyforms, init=Int)
    @assert commontype <: Union{Integer,Rational} "Only integer and rational coefficients are supported as input."

    # Convert all to DP.Polynomial, so that coefficients are of same type,
    # and constants are treated as polynomials
    # We also need this because Groebner does not support abstract types as input
    polynoms = Vector{DP.Polynomial{true,commontype}}(undef, length(sympolys))
    for (i, pf) in enumerate(polyforms)
        polynoms[i] = underlyingpoly(pf)
    end

    polynoms, pvar2sym, sym2term
end

#=
Converts an array of AbstractPolynomialLike`s into an array of
symbolic expressions mapping variables w.r.t pvar2sym
=#
function poly_to_symbol(polys, pvar2sym, sym2term, ::Type{T}) where {T}
    map(f -> PolyForm{T}(f, pvar2sym, sym2term), polys)
end

"""

groebner_basis(polynomials; reduced=true)

Computes a Groebner basis of the ideal generated by the given `polynomials`.

If `reduced` is set, the basis is reduced and *unique*.

"""
function groebner_basis(polynomials; reduced=true)
    polynoms, pvar2sym, sym2term = symbol_to_poly(polynomials)

    basis = groebner(polynoms, reduced=reduced)

    # polynomials is nonemtpy
    T = symtype(first(polynomials))
    poly_to_symbol(basis, pvar2sym, sym2term, T)
end


text:
# estimate with MCMC (Turing) for prior (Laplas and half Cauchy
# https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.Laplace
# half-Cauchy not implemented yet
# https://github.com/JuliaStats/Distributions.jl/issues/124
using Turing
using Distributions
using Random
using LinearAlgebra
using StatsPlots

# Generate moc data
n = 50
dim = 1
Random.seed!(99)
u = rand(Normal(0, 0.7), n) # true value σ² = 0.49
x = rand(Uniform(-sqrt(3.0), sqrt(3.0)), n, dim) # x generated from Uniform[-√3, √3]
X = [ones(n) x] # dependent variable matrix with ones
β = Vector(1:1:dim+1) # true coefficients
y = X*β + u # true data generating process

# set prior (Laplace for β, half-Cauchy for σ²)


text:
# THIS FILE IS DEACTIVATED

using Random
using Plots
using VisualRegressionTests

# run a visual regression test comparing the output to the saved reference png
function dotest(testname, func)
    Random.seed!(1)
    reffn = joinpath(refdir, "$testname.png")
    vtest = VisualTest(func, reffn)
    result = test_images(vtest)
    @test success(result)
end

# builds a testable function which saves a png to the location `fn`
# use: VisualTest(@plotfunc(plot(rand(10))), "/tmp/tmp.png")
macro plottest(testname, expr)
    esc(quote
        dotest(string($testname), fn -> begin
            $expr
            png(fn)
        end)
    end)
end

# don't let pyplot use a gui... it'll crash
# note: Agg will set gui -> :none in PyPlot
# - This unglyness needs to change
had_mlp = haskey(ENV, "MPLBACKEND")
_old_env = get(ENV, "MPLBACKEND", "")
ENV["MPLBACKEND"] = "Agg"
import PyPlot
info("Matplotlib version: $(PyPlot.matplotlib[:__version__])")
pyplot(size=(200,150), reuse=true)

refdir = joinpath(dirname(@__FILE__), "refimg")

@plottest "dynmv" begin
	history = ValueHistories.MVHistory(QHistory)
	for i=1:100
		x = 0.1i
		push!(history, :a, x, sin(x))
		push!(history, :wrongtype, x, "$(sin(x))")
		if i % 10 == 0
			push!(history, :b, x, cos(x))
		end
	end
	plot(history)
end

@plottest "dynmv_sub" begin
	history = ValueHistories.MVHistory()
	for i=1:100
		x = 0.1i
		push!(history, :a, x, sin(x))
		push!(history, :wrongtype, x, "$(sin(x))")
		if i % 10 == 0
			push!(history, :b, x, cos(x))
		end
	end
	plot(history, layout=2)
end

@plottest "queueuv" begin
	history = ValueHistories.QHistory(Int)
	for i = 1:100
		push!(history, i, 2i)
	end
	plot(history)
end

@plottest "vectoruv" begin
	history = ValueHistories.History(Int)
	for i = 1:100
		push!(history, i, 100-i)
	end
	plot(history)
end

@plottest "uv_vector" begin
	history1 = ValueHistories.History(Int)
	history2 = ValueHistories.QHistory(Int)
	for i = 1:100
		push!(history1, i, 2i)
		push!(history2, i, 100-i)
	end
	plot([history1, history2], layout = 2)
end

if had_mlp
    ENV["MPLBACKEND"] = _old_env
else
    delete!(ENV, "MPLBACKEND")
end


text:
function write_analog(t::NIDAQ.AOTask, data::Matrix{Float64}; timeout::Float64=-1.0, wait_for_space = true)
    num_samples_per_chan::Int32 = size(data, 1)
    data = reshape(data, length(data))
    num_samples_per_chan_written = Int32[0]
    buf_space_free = Ref{UInt32}(1)
    if wait_for_space #the function for checking write space returns 0 if the task has not been started yet.
        while true
            NIDAQ.catch_error(NIDAQ.GetWriteSpaceAvail(t.th, buf_space_free))
            #NOTE: with usb-6002 devices there seems to be more space available than was set by NIDAQ.DAQmxCfgOutputBuffer.  Strange but eems harmless enough.
            if buf_space_free[] >= num_samples_per_chan
                break
            end
            sleep(0.2) #buf size is currently 1s worth of samples, so this shouldn't need to be changed unless that changes
        end
    end
    NIDAQ.catch_error(NIDAQ.WriteAnalogF64(t.th,
        num_samples_per_chan,
        reinterpret(Bool32, UInt32(false)),
        timeout,
        reinterpret(Bool32,NIDAQ.Val_GroupByChannel),
        Ref(data, 1),
        Ref(num_samples_per_chan_written, 1),
        reinterpret(Ptr{Bool32},C_NULL)))
    if num_samples_per_chan_written[1] != num_samples_per_chan
        error("Wrote $num_samples_per_chan_written samples instead of the $num_samples_per_chan samples requested")
    end
end
write_analog(t::NIDAQ.AOTask, data::Vector{Float64}; timeout::Float64=-1.0) = write_analog(t, reshape(data,(length(data),1)), timeout=timeout)

#Probably no need for two names.  Just make this and write_analog two implementations of a _write function
function write_digital(t::NIDAQ.DOTask, data::Matrix{UInt8}; timeout::Float64=-1.0, wait_for_space = true)
    error("Not yet implemented")
end
write_digital(t::NIDAQ.DOTask, data::Vector{UInt8}; timeout::Float64=-1.0) = write_digital(t, reshape(data,(length(data),1)), timeout=timeout)

#returns a DAQ task handle
function prepare_ao(coms, bufsz::Int, trigger_terminal::String, clock_source)
    print("preparing DAQ for ao...\n")
    chns = map(daq_channel, coms)
    nsamps = length(first(coms))
    nchans = length(coms)
    usb_req_size = Ref{UInt32}(0)
    rig = rig_name(first(coms))
    if !in(split(DEVICE_PREFIX, "/")[1], NIDAQ.devices())
        error("Device $DEVICE_PREFIX not detected")
    end
    tsk = analog_output(DEVICE_PREFIX * chns[1])
    for i = 2:length(chns)
        analog_output(tsk, DEVICE_PREFIX * chns[i])
        NIDAQ.setproperty!(tsk, DEVICE_PREFIX * chns[i], "Max", 10.0) #"/ao0:1", "Max", 10.0)
    end
    if rig == "dummy-6002"
        #note: default Xfer size seems to be 8000 bytes (4000 samples)
        NIDAQ.catch_error(NIDAQ.GetAOUsbXferReqSize(tsk.th, Vector{UInt8}(DEVICE_PREFIX * chns[1]), usb_req_size))
        #See article here http://digital.ni.com/public.nsf/allkb/B7B47F50F9813DFD862575890054EF7C
        usb_req_size = Ref{UInt32}(usb_req_size[] * UInt32(2))
        for i = 1:length(chns)
            print("    Set Usb transfer request size for AO\n")
            NIDAQ.catch_error(NIDAQ.SetAOUsbXferReqSize(tsk.th, Vector{UInt8}(DEVICE_PREFIX * chns[i]), usb_req_size[]))
        end
    end
    #creating clk_b this way works around a bizarre bug that prevents precompilatoin when converting clock_source with other methods
    clk_b = UInt8[]
    if isempty(clock_source)
        clk_b = C_NULL
    elseif clock_source[1:2] == "ao"
        clk_b = convert(Ref{UInt8}, b"OnboardClock")
    else
        for ch in clock_source
            push!(clk_b, ch)
        end
	clk_b = convert(Ref{UInt8}, clk_b)
    end
    NIDAQ.catch_error(NIDAQ.CfgSampClkTiming(tsk.th,
            clk_b,
            Float64(ustrip(samprate(first(coms)))),
            NIDAQ.Val_Rising,
            NIDAQ.Val_FiniteSamps,
            UInt64(nsamps)))
    disallow_regen = Int32(10158) #default is allow regen (10097).  We want to disallow it (10158)
    NIDAQ.catch_error(NIDAQ.DAQmxSetWriteRegenMode(tsk.th, disallow_regen))
    print("AO buffer size (nsamps): $bufsz\n")
    NIDAQ.catch_error(NIDAQ.DAQmxCfgOutputBuffer(tsk.th, UInt32(bufsz)))
    if trigger_terminal != "disabled"
        print("    setting up digital start trigger for analog output...\n")
	#(drop trailing "/" from device name)
        props = NIDAQ.getproperties(String(split(DEVICE_PREFIX, "/")[1]))#; warning=false)
        terms = props["Terminals"][1]
        if !in("/" * DEVICE_PREFIX * trigger_terminal, terms)
            error("Terminal $trigger_terminal is does not exist.  Check valid terminals with NIDAQ.getproperties(dev)[\"Terminals\"]")
        end
        #The below three commands should be equivalent to calling CfgDigEdgeStartTrig
        #NIDAQ.catch_error(NIDAQ.SetStartTrigType(tsk.th, NIDAQ.Val_DigEdge))
        #NIDAQ.catch_error(NIDAQ.SetDigEdgeStartTrigSrc(tsk.th, "/" * DEVICE_PREFIX * trigger_terminal))
        #NIDAQ.catch_error(NIDAQ.SetDigEdgeStartTrigEdge(tsk.th, NIDAQ.Val_Rising))
        NIDAQ.catch_error(NIDAQ.CfgDigEdgeStartTrig(tsk.th, Ref(codeunits("/" * DEVICE_PREFIX * trigger_terminal),1), NIDAQ.Val_Rising))
    end
    return tsk
end

function prepare_do(coms, bufsz::Int, trigger_terminal::String; clock_source = clock)
    error("Not yet implemented")
end

function _output_analog_signals(coms::AbstractVector{T}, samps_per_write::Int, trigger_terminal::String, ready_chan::RemoteChannel, clock) where T<:ImagineSignal
    if length(finddigital(coms)) != 0
        error("Please provide only analog ImagineSignals")
    end
    if !all(map(isoutput, coms))
        error("Please provide only output ImagineSignals")
    end
    tsk = prepare_ao(coms, 2*samps_per_write, trigger_terminal, clock)
    nchans = length(coms)
    finished = false
    curstart = 1
    nsamps = length(first(coms))
    write_buffer = zeros(Float64, samps_per_write, nchans)
    write_incr = samps_per_write - 1
    wrote_once = false
    try
        while !finished
            curstop = curstart + write_incr
            if curstop >= nsamps
                finished = true
                curstop = min(nsamps, curstop)
            end
            nsamps_to_write = curstop - curstart + 1
            for i = 1:nchans
                write_buffer[1:nsamps_to_write, i] = ustrip.(get_samples(coms[i], curstart, curstop; sampmap = :volts))
            end
            #now write to daq buffer
            if wrote_once
                write_analog(tsk, write_buffer[1:nsamps_to_write,:], timeout=-1.0)
            else
                write_analog(tsk, write_buffer[1:nsamps_to_write,:], timeout=-1.0, wait_for_space = false)
            end
            print("So far $curstop of $nsamps samples have been written...\n")
            curstart = curstop + 1
            if !wrote_once
                start(tsk)
                put!(ready_chan, myid())
                wrote_once = true
            end
        end
    catch
        clear(tsk)
        rethrow()
    end
    print("Waiting for AO task to complete\n")
    NIDAQ.catch_error(NIDAQ.WaitUntilTaskDone(tsk.th, -1.0))
    stop(tsk)
    clear(tsk)
    print("Finished.\n")
    return ImagineSignal[]
end

function _output_digital_signals(coms::AbstractVector{T}, samps_per_write::Int, trigger_terminal::String, clock::AbstractString) where T>:ImagineSignal
    if length(findanalog(coms)) != 0
        error("Please provide only digital ImagineSignals")
    end
    if !all(map(isoutput, coms))
        error("Please provide only output ImagineSignals")
    end
    tsk = prepare_do(coms, 2*samps_per_write, trigger_terminal; clock_source = clock)
    nchans = length(coms)
    finished = false
    curstart = 1
    nsamps = length(first(coms))
    write_buffer = zeros(Float64, samps_per_write, nchans)
    write_incr = samps_per_write - 1
    wrote_once = false
    try
        while !finished
            curstop = curstart + write_incr
            if curstop >= nsamps
                finished = true
                curstop = min(nsamps, curstop)
            end
            nsamps_to_write = curstop - curstart + 1
            for i = 1:nchans
                write_buffer[1:nsamps_to_write, i] = ustrip.(get_samples(coms[i], curstart, curstop; sampmap = :volts))
            end
            #now write to daq buffer
            if wrote_once
                write_digital(tsk, write_buffer[1:nsamps_to_write,:], timeout=-1.0)
            else
                write_digital(tsk, write_buffer[1:nsamps_to_write,:], timeout=-1.0, wait_for_space = false)
            end
            curstart = curstop + 1
            if !wrote_once
                start(tsk)
                set_ready()
                wrote_once = true
            end
        end
    catch
        clear(tsk)
        rethrow()
    end
    NIDAQ.WaitUntilTaskDone(tsk.th, -1.0) #wait forever
    stop(tsk)
    clear(tsk)
    return myid()
end

0 #Without this zero, the (fetched) remotecall that includes this file will return the _ttl_pulse function and cause an error.  Instead this makes it return 0


text:
seats = readlines("day5/input.txt")

function get_seat_id(seat::String)
    value = 128
    row_value = 0
    for i in 1:7
        value ÷= 2
        if seat[i] == 'B'
            row_value += value
        end
    end
    value = 8
    col_value = 0
    for i in 8:10
        value ÷= 2
        if seat[i] == 'R'
            col_value += value
        end
    end
    row_value*8 + col_value
end

println(maximum(get_seat_id, seats))

seat_ids = Set(map(get_seat_id, seats))

for i in 100:900
    if !(i in seat_ids)
        println(i)
        break
    end
end

# i coded this while watching a christmas romcom
# and it was far more pleasant than yesterday :)


text:
# context management and global state

# to avoid CUDA-style implicit state, where operations can fail if they are accidentally
# executed in the wrong context, ownership should always be encoded in each object.
# the accessors below should only be used to determine initial ownership.

export driver, driver!, device, device!, context, context!, global_queue, synchronize

function driver()
    get!(task_local_storage(), :ZeDriver) do
        first(drivers())
    end
end

function driver!(drv::ZeDriver)
    task_local_storage(:ZeDriver, drv)
    delete!(task_local_storage(), :ZeDevice)
end

function device()
    get!(task_local_storage(), :ZeDevice) do
        first(devices(driver()))
    end
end

function device!(drv::ZeDevice)
    task_local_storage(:ZeDevice, drv)
end

const global_contexts = Dict{ZeDriver,ZeContext}()

function context()
    get!(task_local_storage(), :ZeContext) do
        get!(global_contexts, driver()) do
            ZeContext(driver())
        end
    end
end

function context!(ctx::ZeContext)
    task_local_storage(:ZeContext, ctx)
end

# the global queue can be used as a default queue to execute operations on,
# guaranteeing expected semantics when using a device on a Julia task.

function global_queue(ctx::ZeContext, dev::ZeDevice)
    # NOTE: dev purposefully does not default to context() or device() to stress that
    #       objects should track ownership, and not rely on implicit global state.
    get!(task_local_storage(), (:ZeCommandQueue, ctx, dev)) do
        ZeCommandQueue(ctx, dev)
    end
end

function oneL0.synchronize()
    oneL0.synchronize(global_queue(context(), device()))
end


text:
tests = {[1, 2, 3],
         primes(10),
         0:2:6,
         [-Inf, 0.0, Inf],
         [π, e, φ, catalan],
         [2015, 5],
         [-sqrt(50.0), 50.0^2],
         }

println("Testing islexfirst:")
for (a, b) in combinations(tests, 2)
    tres = islexfirst(a, b) ? " is " : " is not "
    tres *= "lexically prior to\n    "
    println("\n    ", a, tres, b)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Affine coupling layer from Dinh et al. (2017)
# Includes 1x1 convolution from in Putzky and Welling (2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export AffineCouplingLayerSLIM

"""
    CS = AffineCouplingLayerSLIM(nx, ny, n_in, n_hidden, batchsize, Ψ; logdet=false, permute=false, k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

 Create an invertible affine SLIM coupling layer.

 *Input*:

 - `nx, ny`: spatial dimensions of input

 - `n_in`, `n_hidden`: number of input and hidden channels

 - `Ψ`: link function

 - `loget`: bool to indicate whether to return the logdet (default is `false`)

 - `permute`: bool to indicate whether to apply a channel permutation (default is `false`)

 - `k1`, `k2`: kernel size of convolutions in residual block. `k1` is the kernel of the first and third
    operator, `k2` is the kernel size of the second operator

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`)

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`)

 *Output*:

 - `CS`: Invertible SLIM coupling layer

 *Usage:*

 - Forward mode: `Y, logdet = CS.forward(X, D, A)`    (if constructed with `logdet=true`)

 - Inverse mode: `X = CS.inverse(Y, D, A)`

 - Backward mode: `ΔX, X = CS.backward(ΔY, Y, D, A)`

 - where `A` is a linear forward modeling operator and `D` is the observed data.

 *Trainable parameters:*

 - None in `CL` itself

 - Trainable parameters in residual block `CL.RB` and 1x1 convolution layer `CL.C`

 See also: [`Conv1x1`](@ref), [`ResidualBlock`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct AffineCouplingLayerSLIM <: NeuralNetLayer
    C::Union{Conv1x1, Nothing}
    RB::Union{ResidualBlock, FluxBlock}
    AN::ActNorm
    Ψ::Function
    logdet::Bool
end

@Flux.functor AffineCouplingLayerSLIM

# Constructor from input dimensions
function AffineCouplingLayerSLIM(nx::Int64, ny::Int64, n_in::Int64, n_hidden::Int64, batchsize::Int64, Ψ::Function;
    k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, logdet::Bool=false, permute::Bool=false)

    # 1x1 Convolution and residual block for invertible layer
    permute == true ? (C = Conv1x1(n_in)) : (C = nothing)
    RB = ResidualBlock(nx, ny, Int(n_in/2), n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, fan=true)
    AN = ActNorm(1)

    return AffineCouplingLayerSLIM(C, RB, AN, Ψ, logdet)
end

# Forward pass: Input X, Output Y
function forward(X::AbstractArray{Float32, 4}, D, J, CS::AffineCouplingLayerSLIM)

    # Get dimensions
    nx, ny, n_s, batchsize = size(X)

    # Permute and split
    isnothing(CS.C) ? (X_ = copy(X)) : (X_ = CS.C.forward(X))
    X1_, X2_ = tensor_split(X_)

    # Gradient
    g = J'*(J*reshape(CS.Ψ(X1_[:,:,1:1,:]), :, batchsize) - D)
    g = reshape(g, nx, ny, 1, batchsize)
    gn = CS.AN.forward(g)
    gs = tensor_cat(gn, X1_[:,:,2:end,:])

    # Coupling layer
    Y1_ = copy(X1_)
    logS_T = CS.RB.forward(gs)
    logS, T = tensor_split(logS_T)
    S = Sigmoid(logS)
    Y2_ = S.*X2_ + T
    Y_ = tensor_cat(Y1_, Y2_)

    isnothing(CS.C) ? (Y = copy(Y_)) : (Y = CS.C.inverse(Y_))
    CS.logdet == true ? (return Y, slim_logdet_forward(S)) : (return Y)
end

# Inverse pass: Input Y, Output X
function inverse(Y::AbstractArray{Float32, 4}, D, J, CS::AffineCouplingLayerSLIM; save=false)

    # Get dimensions
    nx, ny, n_s, batchsize = size(Y)
    isnothing(CS.C) ?  (Y_ = copy(Y)) : (Y_ = CS.C.forward(Y))

    # Gradient
    Y1_, Y2_ = tensor_split(Y_)
    X1_ = copy(Y1_)
    g = J'*(J*reshape(CS.Ψ(X1_[:,:,1:1,:]), :, batchsize) - D)
    g = reshape(g, nx, ny, 1, batchsize)
    gn = CS.AN.forward(g)
    gs = tensor_cat(gn, X1_[:,:,2:end,:])

    # Coupling layer
    logS_T = CS.RB.forward(gs)
    logS, T = tensor_split(logS_T)
    S = Sigmoid(logS)
    X2_ = (Y2_ - T) ./ (S .+ eps(1f0)) # add epsilon to avoid
    X_ = tensor_cat(X1_, X2_)

    isnothing(CS.C) ? (X = copy(X_)) : (X = CS.C.inverse(X_))
    save == true ? (return X, X1_, X2_, gn, gs, S) : (return X)
end

# Backward pass: Input (ΔY, Y), Output (ΔX, X)
function backward(ΔY::AbstractArray{Float32, 4}, Y::AbstractArray{Float32, 4}, D, J, CS::AffineCouplingLayerSLIM; permute=false)

    # Recompute forward states
    X, X1_, X2_, gn, gs, S  = inverse(Y, D, J, CS; save=true)
    nx, ny, n_s, batchsize = size(Y)

    # Backpropagation
    isnothing(CS.C) ? (ΔY_ = copy(ΔY)) : (ΔY_ = CS.C.forward((ΔY, Y))[1])
    ΔY1_, ΔY2_ = tensor_split(ΔY_)
    ΔT = copy(ΔY2_)
    ΔS = ΔY2_ .* X2_
    CS.logdet == true && (ΔS -= slim_logdet_backward(S))
    ΔX2_ = ΔY2_ .* S
    Δgs = CS.RB.backward(cat(SigmoidGrad(ΔS, S), ΔT; dims=3), gs)
    Δgn = Δgs[:,:,1:1,:]
    Δg = CS.AN.backward(Δgn, gn)[1]
    Jg = J*reshape(Δg, :, batchsize)
    ΔD = -Jg
    ΔX1_ = tensor_cat(reshape(J'*Jg, nx, ny, 1, batchsize), Δgs[:,:,2:end,:]) + ΔY1_
    ΔX_ = tensor_cat(ΔX1_, ΔX2_)
    isnothing(CS.C) ? (ΔX = copy(ΔX_)) : (ΔX = CS.C.inverse((ΔX_, tensor_cat(X1_, X2_)))[1])

    return ΔX, ΔD, X
end

# Clear gradients
function clear_grad!(CS::AffineCouplingLayerSLIM)
    ~isnothing(CS.C) && clear_grad!(CS.C)
    clear_grad!(CS.RB)
    clear_grad!(CS.AN)
    CS.AN.s.data = nothing
    CS.AN.b.data = nothing
end

# Get parameters
function get_params(CS::AffineCouplingLayerSLIM)
    p = get_params(CS.RB)
    ~isnothing(CS.C) && (p = cat(p, get_params(CS.C); dims=1))
    return p
end

# Put parameters
function put_params!(CS::AffineCouplingLayerSLIM, Params::Array{Any,1})
    put_params!(CS.RB, Params[1:5])
    ~isnothing(CS.C) && put_params!(CS.C, Params[6:end])
end

# Logdet
slim_logdet_forward(S) = sum(log.(abs.(S))) / size(S, 4)
slim_logdet_backward(S) = 1f0./ S / size(S, 4)


text:
# Auxillary functions
const _uint_bit_length = sizeof(UInt) * 8
const _div_uint_size_shift = Int(log2(_uint_bit_length))
@inline _mul2(i::Integer) = i << 1
@inline _div2(i::Integer) = i >> 1
@inline _map_to_index(i::Integer) = _div2(i - 1)
@inline _map_to_factor(i::Integer) = _mul2(i) + 1
@inline _mod_uint_size(i::Integer) = i & (_uint_bit_length - 1)
@inline _div_uint_size(i::Integer) = i >> _div_uint_size_shift
@inline _get_chunk_index(i::Integer) = _div_uint_size(i + (_uint_bit_length - 1))
@inline _get_bit_index_mask(i::Integer) = UInt(1) << _mod_uint_size(i - 1)

function clear_factors_while!(arr::Vector{UInt}, factor_index::Integer, max_index::Integer)
  factor = _map_to_factor(factor_index)
  index = _div2(factor * factor)
  while index <= max_index
    @inbounds arr[_get_chunk_index(index)] |= _get_bit_index_mask(index)
    index += factor
  end
  return arr
end
function clear_factors_turbo!(arr::Vector{UInt}, factor_index::Integer, max_index::Integer)
  factor = _map_to_factor(factor_index)
  factor < _uint_bit_length && error("Factor must be greater than UInt bit length to avoid memory pendencies")
  @turbo for index in _div2(factor * factor):factor:max_index
    @inbounds arr[(index + 63) >> _div_uint_size_shift] |= 1 << ((index - 1) & 63)
  end
  return arr
end
function clear_factors_turbo_u4!(arr::Vector{UInt}, factor_index::Integer, max_index::Integer)
  factor = _map_to_factor(factor_index)
  factor < _uint_bit_length && error("Factor must be greater than UInt bit length to avoid memory pendencies")
  @turbo unroll=4 for index in _div2(factor * factor):factor:max_index
    @inbounds arr[(index + 63) >> _div_uint_size_shift] |= 1 << ((index - 1) & 63)
  end
  return arr
end

@testset "steprange" begin
  x0 = rand(UInt, cld(500_000, sizeof(UInt) * 8));
  x1 = copy(x0);
  x2 = copy(x0);
  clear_factors_while!(x0, 202, 500_000)
  clear_factors_turbo!(x1, 202, 500_000)
  clear_factors_turbo_u4!(x2, 202, 500_000)
  @test x0 == x1 == x2
end


text:
using TopOpt

problem = PointLoadCantilever(Val{:Linear}, (160, 40), (1.0, 1.0), 1.0, 0.3, 1.0)
solver = FEASolver(Displacement, Direct, problem)
obj = ComplianceFunction(problem, solver, filtering = true, rmin = 3.0)
constr = Constraint(VolumeFunction(problem, solver, 0.5))
optimizer = MMAOptimizer(obj, constr, MMA.MMA87(), ConjugateGradient(), maxiter=1000);


text:
struct SpectralConv{N,M,F,A,V,S}
    σ::F
    weight::A
    bias::V
    stride::NTuple{N,Int}
    pad::NTuple{M,Int}
    dilation::NTuple{N,Int}
    groups::Int
    spectralnormalizer::S
end

function SpectralConv(
    k::NTuple{N,Integer},
    ch::Pair{<:Integer,<:Integer},
    σ = identity;
    stride = 1,
    pad = 0,
    dilation = 1,
    groups = 1,
    init = glorot_normal,
    iterations = 1,
    norm_multiplier = 0.95,
    weight = Flux.convfilter(k, (ch[1] ÷ groups => ch[2])),
    bias = true,
    device = cpu,
) where {N}

    stride = Flux.expand(Val(N), stride)
    dilation = Flux.expand(Val(N), dilation)
    pad = Flux.calc_padding(SpectralConv, pad, size(weight)[1:N], dilation, stride)
    bias = create_bias(weight, bias, size(weight, N + 2))
    norm_multiplier = convert(Float32, norm_multiplier)

    kernel_shape = size(reshape(weight, ch[2], :))
    spectralnormalizer = SpectralNormalizer(
        kernel_shape,
        iterations = iterations,
        norm_multiplier = norm_multiplier,
        init = init,
        device = device,
    )

    return SpectralConv(σ, weight, bias, stride, pad, dilation, groups, spectralnormalizer)
end

@functor SpectralConv (weight, bias)

function (c::SpectralConv)(x)
    # update
    kernel_shape = size(c.weight)
    out_channels = kernel_shape[end]
    kernel = reshape(c.weight, out_channels, :)
    # Bound the kernel
    bounded_kernel = c.spectralnormalizer(kernel)
    bounded_kernel = reshape(bounded_kernel, kernel_shape)

    σ, b = c.σ, reshape(c.bias, ntuple(_ -> 1, length(c.stride))..., :, 1)
    output = Flux.NNlib.conv(
        x,
        bounded_kernel;
        stride = c.stride,
        pad = c.pad,
        dilation = c.dilation,
        groups = c.groups,
    )
    σ.(output .+ b)
end


text:
using ArgParse
using DataFrames
using Statistics

import Colors

function parse_preview_commandline(args::Union{Nothing, Array{String, 1}}=nothing)
    s = ArgParseSettings(prog="baysor preview")
    @add_arg_table! s begin
        "--config", "-c"
            help = "TOML file with config"
        "--x-column", "-x"
            help = "Name of x column. Overrides the config value."
        "--y-column", "-y"
            help = "Name of gene column. Overrides the config value."
        "--z-column", "-z"
            help = "Name of gene column. Overrides the config value."
        "--gene-column", "-g"
            help = "Name of gene column. Overrides the config value."
        "--min-molecules-per-cell", "-m"
            help = "Minimal number of molecules for a cell to be considered as real. It's an important parameter, as it's used to infer several other parameters. Overrides the config value."
            arg_type = Int
        "--min-pixels-per-cell"
            help = "Minimal number of pixels per cell. Used to estimate size of the dataset plot."
            arg_type = Int
            default = 15
        "--exclude-genes"
            help = "Comma-separated list of genes to ignore during segmentation"
        "--output", "-o"
            help = "Name of the output file or path to the output directory"
            default = "preview.html"
        "coordinates"
            help = "CSV file with coordinates of transcripts and gene type"
            required = true
        # TODO: add gene-composition-neighborhood
        # TODO: add verbosity level
    end

    return (args === nothing) ? parse_args(s) : parse_args(args, s)
end

function parse_preview_configs(args::Union{Nothing, Array{String, 1}}=nothing)
    r = parse_preview_commandline(args)
    if r === nothing
        return nothing
    end

    if r["config"] !== nothing
        extend_params_with_config!(r, parse_toml_config(r["config"]))
    else
        extend_params_with_config!(r, get_default_config())
    end

    for k in ["gene-column", "x-column", "y-column", "z-column"]
        r[k] = Symbol(r[k])
    end

    if isdir(r["output"]) || isdirpath(r["output"])
        r["output"] = joinpath(r["output"], "preview.html")
    end

    return r
end

function run_cli_preview(args::Union{Nothing, Array{String, 1}}=nothing)
    args = parse_preview_configs(args)
    (args !== nothing) || return 1

    # Set up logger

    log_file = open(append_suffix(args["output"], "preview_log.log"), "w")
    Base.CoreLogging.global_logger(DoubleLogger(log_file, stdout; force_flush=true))

    # Run preview

    @info get_baysor_run_str()
    @info "Loading data..."
    args["min-molecules-per-gene"] = 0
    df_spatial, gene_names = load_df(args)

    @info "Loaded $(size(df_spatial, 1)) transcripts"
    if size(df_spatial, 1) != size(unique(df_spatial), 1)
        @warn "$(size(df_spatial, 1) - size(unique(df_spatial), 1)) records are duplicates. You may need to filter them beforehand."
    end

    @info "Estimating noise level"
    confidence_nn_id = default_param_value(:confidence_nn_id, args["min-molecules-per-cell"])
    edge_lengths, confidences, d1, d2 = append_confidence!(df_spatial, nn_id=confidence_nn_id) # TODO: use segmentation mask if available here
    @info "Done"

    @info "Estimating local neighborhoods"

    if args["gene-composition-neigborhood"] === nothing
        args["gene-composition-neigborhood"] = default_param_value(:composition_neighborhood, args["min-molecules-per-cell"], n_genes=length(gene_names))
    end

    @info "Estimating local colors"
    gene_colors = gene_composition_colors(df_spatial, args["gene-composition-neigborhood"])

    ## Plot

    @info "Building transcript plots"
    gc_plot = plot_dataset_colors(df_spatial, gene_colors; min_molecules_per_cell=args["min-molecules-per-cell"],
        min_pixels_per_cell=args["min-pixels-per-cell"], title="Local expression similarity")

    conf_colors = map_to_colors(confidences, lims=(0.0, 1.0), palette=Colors.diverging_palette(10, 250, s=0.75, w=1.0));
    cc_plot = plot_dataset_colors(df_spatial, conf_colors[:colors]; min_molecules_per_cell=args["min-molecules-per-cell"],
        min_pixels_per_cell=args["min-pixels-per-cell"], title="Transcript confidence")

    @info "Building gene structure plot"
    vega_plots = Dict{String, VegaLite.VLSpec}()

    vega_plots["vg_gene_structure"] = plot_gene_structure(df_spatial, gene_names)
    vega_plots["vg_num_trans"] = plot_num_transcript_overview(df_spatial.gene, confidences, gene_names)
    vega_plots["vg_noise_dist"] = plot_noise_estimation_diagnostics(edge_lengths, confidences, d1, d2, confidence_nn_id=confidence_nn_id)

    @info "Plotting"

    open(args["output"], "w") do io
        print(io, """
            <!DOCTYPE html>

            $(vega_header("Report", ""))
            $(vega_style())

Content

Transcript plots

Noise level

Gene structure

            """)

        println(io, "
Transcript plots

")
        show(io, MIME("text/html"), gc_plot)
        println(io, "
")

        println(io, "
\n
Noise level

")
        show(io, MIME("text/html"), cc_plot)

        println(io, "
")
        println(io, "

")

        println(io, "
")
        println(io, "Minimal noise level=$(round(100 * mean(confidences .< 0.01), sigdigits=3))%. ",
            "Expected noise level=$(round(100 * mean(1 .- confidences), sigdigits=2))%.")

        println(io, "
")
        println(io, "
\n
Gene structure

")
        println(io, "

")

        println(io, "
")
        println(io, "

")

        println(io, "")

        println(io, vega_plot_html(vega_plots))
        println(io, "")
    end

    @info "All done!"

    close(log_file)

    return 0
end

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
for T in union(InteractiveUtils.subtypes(ItoIsometryMethod))
    @eval (::$T)(val1::Real,val2::Real) = val1 * val2
end
function (iim::Trapezoidal)(f1::AbstractVector{<:Real},f2::AbstractVector{<:Real})
    sum(f1 .* f2 .* iim) * iim.dt
end
function (iim::Trapezoidal{k})(f::Function,t0::Real,t1::Real) where k
    toFunctionMXelements(f.(t0:iim.dt:(t1+100eps(t1))),iim)
end
###############################################################################
############## Submatrix from the stochastic proportional terms ###############
###############################################################################
# TODO: increase the performance for constant matrices!
function (method::SemiDiscretization)(A::stCoeffMX{d,<:ProportionalMX}, rst::AbstractResult{d}) where d
    stSubMX.(Ref(A.nID),
        Ref([subMxRange(0, d)]),
        [[rst.itoisometrymethod(t->exp(rst.A_avgs[i]*(t1-t))*A(t)*exp(rst.A_avgs[i]*(t-rst.ts[i])),rst.ts[i],t1)]
          for (i,t1) in enumerate(rst.ts[2:end])])
end
# function (method::SemiDiscretization)(A::stCoeffMX{<:ProportionalMX{<:AbstractMatrix}}, rst::AbstractResult)
#     ... function body ...
# end
# function (method::SemiDiscretization)(A::stCoeffMX{<:ProportionalMX{<:Function}}, rst::AbstractResult)
#     ... function body ...
# end

###############################################################################
################# Submatrix from the stochastic delayed terms #################
###############################################################################
# TODO: increase the performance for constant matrices!
function (method::SemiDiscretization{<:NumericSD})(τ::Real, B::stCoeffMX{d,<:DelayMX}, rst::AbstractResult{d}) where {d,T}
    (τerr, ranges) = method(τ, rst)
    MXs = [[rst.itoisometrymethod((t -> (exp(rst.A_avgs[i] * (rst.ts[i + 1] - t)) * B(t) * lagr_el0(methodorder(method), k, τerr, rst.ts[i + 1], t))), rst.ts[i], rst.ts[i + 1]) for k in 0:methodorder(method)] for i in 1:rst.n_steps]
    stSubMX.(Ref(B.nID),Ref(ranges), MXs)
end
function (method::SemiDiscretization{<:NumericSD})(τs::Vector{<:Real}, B::stCoeffMX{d,<:DelayMX}, rst::AbstractResult{d}) where {d,T}
    (τerrs, rangess) = method(τs, rst)
    MXs = [[rst.itoisometrymethod((t -> (exp(rst.A_avgs[i] * (rst.ts[i + 1] - t)) * B(t) * lagr_el0(methodorder(method), k, τerrs[i], rst.ts[i + 1], t))), rst.ts[i], rst.ts[i + 1]) for k in 0:methodorder(method)] for i in 1:rst.n_steps]
    stSubMX.(Ref(B.nID),rangess, MXs)
end
###############################################################################
# To distinguish between τ as constant vs τ as function
function (method::SemiDiscretization{<:NumericSD})(B::stCoeffMX{d,<:DelayMX{d,<:Real,T}}, rst::AbstractResult{d}) where {d,T}
    method(B.cMX.τ.τ, B, rst)
end

function (method::SemiDiscretization{<:NumericSD})(B::stCoeffMX{d,<:DelayMX{d,<:Function,T}}, rst::AbstractResult{d}) where {d,T}
    τis = [quadgk(B.cMX.τ, rst.ts[i], rst.ts[i + 1])[1] / method.Δt for i in 1:rst.n_steps]
    method(τis, B, rst)
end

###############################################################################
################ Subvectors from the stochastic additive terms ################
###############################################################################
# TODO: increase the performance for constant matrices and vectors!
function (method::SemiDiscretization{<:NumericSD})(σV::stAdditive{d,Additive{d,T}}, rst::AbstractResult{d}) where {d,T}
    [stSubV(σV.nID, rst.itoisometrymethod((t -> (exp(rst.A_avgs[i] * (rst.ts[i + 1] - t)) * σV(t))), rst.ts[i], rst.ts[i + 1])) for i in 1:rst.n_steps]
end