text:
struct DatabaseLookup <: Expanders.ExpanderPipeline end

Selectors.order(::Type{DatabaseLookup}) = 0.5
Selectors.matcher(::Type{DatabaseLookup}, node, page, doc) = match_kw(node)

const regex_pattern = r"\@example_database\(([\"a-zA-Z_0-9.+ ]+)(?:, )?(\d+|plot|code|stepperplot)?\)"

match_kw(x::String) = occursin(regex_pattern, x)
match_kw(x::Paragraph) = any(match_kw, x.content)
match_kw(x::Any) = false

# ============================================= Simon's implementation
function look_up_source(database_key)
    entries = findall(x-> lowercase(x.title) == lowercase(database_key), database)
    # current implementation finds titles, but we can also search for tags too
    isempty(entries) && error("No entry found for database reference $database_key")
    length(entries) > 1 && error("Multiple entries found for database reference $database_key")
    sprint() do io
        idx = entries[1]
        print_code(
            io, database[entries[1]],
            scope_start = "",
            scope_end = "",
            indent = "",
            resolution = (entry)-> "resolution = (500, 500)",
            outputfile = (entry, ending)-> string(entry.unique_name, ending)
        )
    end
end

function Selectors.runner(::Type{DatabaseLookup}, x, page, doc)
    matched = nothing
    for elem in x.content
        if isa(elem, AbstractString)
            matched = match(regex_pattern, elem)
            matched != nothing && break
        end
    end
    matched == nothing && error("No match: $x")
    capture = matched[1]
    embed = matched[2]

    subidx = nothing
    # we also allow to reference a subsection annotated with @substep with an
    # inter: @example_database("Title", 2)
    if embed isa AbstractString && all(isnumeric, embed)
        subidx = parse(Int, embed)
        embed = nothing
    end
    # The sandboxed module -- either a new one or a cached one from this page.
    database_keys = lowercase.(filter(x-> !(x in ("", " ")), split(capture, '"')))
    # info("$database_keys with embed option $embed")
    map(database_keys) do database_key
        is_stepper = false

        # buffer for content
        content = []

        # bibliographic stuff
        idx = findall(x-> lowercase(x.title) == database_key, database)
        isempty(idx) && error("can't find example $database_key")
        entry = database[idx[1]]
        uname = string(entry.unique_name)
        lines = entry.file_range
        tags = entry.tags
        # print source code for the example
        if embed == nothing || isequal(embed, "code")
            source = look_up_source(database_key)
            if subidx != nothing
                source = split(source, "@substep", keepempty = false)[subidx]
            end
            src_code = Markdown.MD(Markdown.Code("julia", source))
            push!(content, src_code)
        end
        media_root = "/home/sd/ReferenceImages/gallery"
        # embed plot for the example
        if (embed == nothing) || isequal(embed, "plot")
            # print to buffer
            str = sprint() do io
                mpath = joinpath(media_root, uname, "media")
                files = readdir(mpath)
                if subidx !== nothing
                    idx = findfirst(files) do file
                        num = match(r"\d+", file)
                        num === nothing && return false
                        subidx == parse(Int, num.match)
                    end
                    if idx !== nothing
                        embed_media(io, [master_url(media_root, joinpath(mpath, files[idx]))])
                    end
                else
                    embed_media(io, master_url.(media_root, joinpath.(mpath, files)))
                end
            end
            # print code for embedding plot
            src_plot = Documenter.Documents.RawHTML(str)
            push!(content, src_plot)
        end
        # finally, map the content back to the page
        page.mapping[x] = content
    end
end

"""
    print_table(io::IO, dict::Dict)

Print a Markdown-formatted table with the entries from `dict` to specified `io`.
"""
function print_table(io::IO, dict::Dict)
    # get max length of the keys
    k = string.("`", collect(keys(dict)), "`")
    maxlen_k = max(length.(k)...)

    # get max length of the values
    v = string.(collect(values(dict)))
    maxlen_v = max(length.(v)...)

    j = sort(collect(dict), by = x -> x[1])

    # column labels
    labels = ["Symbol", "Description"]

    # print top header
    print(io, "|")
    print(io, "$(labels[1])")
    print(io, " "^(maxlen_k - length(labels[1])))
    print(io, "|")
    print(io, "$(labels[2])")
    print(io, " "^(maxlen_v - length(labels[2])))
    print(io, "|")
    print(io, "\n")

    # print second line (toprule)
    print(io, "|")
    print(io, "-"^maxlen_k)
    print(io, "|")
    print(io, "-"^maxlen_v)
    print(io, "|")
    print(io, "\n")

    for (idx, entry) in enumerate(j)
        print(io, "|")
        print(io, "`$(entry[1])`")
        print(io, " "^(maxlen_k - length(string(entry[1])) - 2))
        print(io, "|")
        print(io, "$(entry[2])")
        print(io, " "^(maxlen_v - length(entry[2])))
        print(io, "|")
        print(io, "\n")
    end
end

"""
    get_video_duration(path::AbstractString)

Returns the duration of the video in seconds (Float32).
Accepted file types: mp4, mkv, and gif.

Requires `ffprobe` (usually comes installed with `ffmpeg`).

Note that while this accepts gif, it will not work to get duration of the gif
(`ffprobe` doesn't support that), so it will just fallback to return 0.5 sec.
"""
function get_video_duration(path::AbstractString)
    !isfile(path) && error("input is not a file!")
    accepted_exts = ("mp4", "gif", "mkv")
    filename = basename(path)
    !(split(filename, ".")[2] in accepted_exts) && error("accepted file types are mp4 and gif!")
    try
        dur = read(`ffprobe -loglevel quiet -print_format compact=print_section=0:nokey=1:escape=csv -show_entries format=duration -i "$(path)"`, String)
        dur = parse(Float32, dur)
    catch e
        @warn("`get_video_duration` on $filename did not work, using fallback video duration of 0.5 seconds")
        dur = 0.5
    end
end


text:
export FusedLasso
export fusedlasso

mutable struct FusedLasso{T} <: AbstractLinearSolver
    # System matrix
    A::Array{T,2}
    # Regularizer
    reg::Vector{Regularization}
    # Solution
    x::Vector{T}
    # Measurement vector
    b::Vector{T}
    # 3D dimensions of solution
    shape::Array{Int64,1}
    # Row energies of system matrix
    rowEnergy
    # user params
    iterations::Int64
    nhood::Array{Int64,2}
    omega::Vector{T}
    gamma::T
    lambda::T
end

function FusedLasso(S::Matrix{T};
                    λ=[0.0002,0.0016],
                    shape::Array{Int64,1}=[size(S,2),1,1],
                    iterations = Int64(50),
                    directions = 3,
                    kappa = T(0.5),
                    gamma = T(10.0^-3),
                    kargs...) where T

    reg = Regularization(["TV","L1"], λ)

    # Determine direction vectors and weights for proximal mapping
  	if directions == 3
              nhood = Int64[1 0 0;0 1 0; 0 0 1]
              omega = T[1,1,1]

  	elseif directions == 13
  		 nhood = Int64[
  		    1 0 0;
  		    0 1 0;
  		    0 0 1;
  		    1 1 0;
  		    1 -1 0;
  		    1 0 1;
  		    1 0 -1;
  		    0 1 1;
  		    0 1 -1;
  		    1 1 1;
  		    1 1 -1;
  		    1 -1 -1;
  		    -1 1 -1
  		   ]
  		omega = T[
  		(2/sqrt(3) - 1)
  		(2/sqrt(3) - 1)
  		(2/sqrt(3) - 1)
  		((3 * sqrt(2) - 2 * sqrt(3)) / 6)
  		((3 * sqrt(2) - 2 * sqrt(3)) / 6)
  		((3 * sqrt(2) - 2 * sqrt(3)) / 6)
  		((3 * sqrt(2) - 2 * sqrt(3)) / 6)
  		((3 * sqrt(2) - 2 * sqrt(3)) / 6)
  		((3 * sqrt(2) - 2 * sqrt(3)) / 6)
  		((3 - 3 * sqrt(2) + sqrt(3)) / 6)
  		((3 - 3 * sqrt(2) + sqrt(3)) / 6)
  		((3 - 3 * sqrt(2) + sqrt(3)) / 6)
  		((3 - 3 * sqrt(2) + sqrt(3)) / 6)
  		]
    else
      @warn "No corresponding neighbourhood for direction number. 3 directions will be used!"
      nhood = Int64[1 0 0;0 1 0; 0 0 1]
      omega = T[1,1,1]
  	end

    # Call constructor interface (function below) and get linear solver object
    solver = FusedLasso(copy(S), reg,
               zeros(eltype(S),shape[1]*shape[2]*shape[3]),
               zeros(eltype(S), div(length(S),shape[1]*shape[2]*shape[3])),shape,
               zeros(eltype(S),size(S,1)),
               iterations,nhood,omega,gamma,kappa
            )

   # Store row energies and normalize rows of system matrix
   normalize!(solver.A, solver.b, solver.rowEnergy)

   return solver
end

#==============================================================================#
# Parameter Settings
#==============================================================================#

function normalize!(A::Array{T,2},b::Vector{T},energy::Vector{T}) where {T<:Real}

    for i=1:size(A,1)
      energy[i] = norm(vec(A[i,:]));
      for j=1:size(A,2)
        A[i,j] = A[i,j]/energy[i]
      end
      b[i] = b[i]/energy[i]
    end

end

function normalize!(b::Vector{T},energy::Vector{T}) where T
    for i=1:length(energy)
      b[i] = b[i] / energy[i]
    end
end

function setMeasurement!(solver::FusedLasso, measurement::Vector)
   m = copy(measurement)
   # Multiply measurement vector with row energies of system matrix
   normalize!(m,solver.rowEnergy)
   # Store result
   solver.b = m
end

#==============================================================================#
# Solve interface
#==============================================================================#

function solve(solver::FusedLasso, measurement::Vector)
    setMeasurement!(solver, measurement)
    solve(solver)
end

function solve(solver::FusedLasso{T}) where T
      shape = solver.shape
      solver.x[:] .= 0.0
      #linearProblem.x = zeros(eltype(linearProblem.x),shape[1]*shape[2]*shape[3])
      solver.x = vec(fusedlasso(
                reshape(solver.A ,size(solver.A,1),shape[1],shape[2],shape[3]),
                solver.b,
                reshape(solver.x,shape[1],shape[2],shape[3]);
                iterations = solver.iterations,
                nhood = solver.nhood,
                omega = solver.omega,
                lambda = solver.lambda,
                alpha = T(solver.reg[1].λ),
                beta = T(solver.reg[2].λ),
                gamma = solver.gamma)   )
end

#==============================================================================#
# Algorithm
#==============================================================================#

"""This function implements the base version of fused lasso reconstruction.

### Keyword/Optional Arguments
* 'iterations::Int64' Maximum number of Iterations.
* 'tol::Float32' Tolerance for the stopping criterion.
* 'nhood::Array{Int64,N,3}' Neighborhood relationships for tv problem.
* 'omega::Array{Float32,N}' Vector with weights for the tv results.
* 'lambda::Float32' Scale for the update step.
* 'alpha::Float32' Weight of the tv term.
* 'beta::Float32' Weight of the l1 term.
* 'gamma::Float32' Weight for the gradient descent step.
"""
function fusedlasso(S::Array{T,4},u::Vector{T},c::Array{T,3};
  iterations = Int64(50),
  tol = T(5.0*10.0^-6),
  nhood = Int64[1 0 0;0 1 0; 0 0 1],
  omega = T[1 1 1],
  lambda = T(1),
  alpha = T(0.000017),
  beta = T(0.001),
  gamma = T(10.0^-3),
  kargs...
  ) where T
  gradientFunc! = gradient_base!

  cSize = size(c)
  N = size(collect(omega),1)
  z = zeros(eltype(c),N+1,cSize[1],cSize[2],cSize[3])

  y = zeros(eltype(c),cSize[1],cSize[2],cSize[3])

  t = ones(T,N+1)/(N+1)
  yTemp = Array{eltype(y)}(undef, size(y))
  #cOld = copy(c)
  zTemp = Array{eltype(yTemp)}(undef, size(yTemp))

  @debug "Residuum pre reconstruction: " residuum(S,c,u)
  for i=1:iterations

      # Calculate gradient
      gradientFunc!(c,S,u,y,gamma)

      # Proximal mapping
      proxmap!(y,yTemp,z,c,N,nhood,alpha,beta,gamma,omega,lambda,t)

      # Pointwise maximum
      pointwise_max!(c,z,N+1,y,yTemp)

      # Update step
      update!(z,N+1,yTemp,c,lambda)

      # Averaging over single results
      c = reshape(c,cSize[1],cSize[2],cSize[3])
      weighted_sum!(t,z,c)
      # Evaluate stopping criterion
      #check_stop(cOld,c,tol) && break;
      #cOld = copy(c)

  end
  @debug "Residuum post reconstruction: " residuum(S,c,u)
  return c

end

"""
This function calculates the error gradient according to y = γA*(Ac-u).
"""
function gradient_base!(c::Array{T,3},A::Array{T,4},u::Array{T,1},y::Array{T,3},α::T) where {T}

  aSize = size(A)
  cSize = size(c)
  uSize = size(u)
  ySize = size(y)

  A = reshape(A,aSize[1],aSize[2]*aSize[3]*aSize[4])
  c = reshape(c,cSize[1]*cSize[2]*cSize[3])
  y = reshape(y,ySize[1]*ySize[2]*ySize[3])

  # Update yTemp with Ac-u
  yTemp = copy(u)
  BLAS.gemv!('N',one(T),A,c,-one(T),yTemp);

  # Update y with γA'yTemp
  BLAS.gemv!('T',α,A,yTemp,zero(T),y);

end

"""
This function performs the proximal mapping.
"""
function proxmap!(y::Array{T,3},yTemp::Array{T,3},z::Array{T,4},c::Array{T,3},
    N,nhood,alpha,beta,gamma,omega,lambda,t) where {T<:Real}
  cSize = size(c)
  ySize = size(y)
  zSize = size(z)

  updateTemp = Array{eltype(y)}(undef, ySize[1]*ySize[2]*ySize[3])

  c = reshape(c,cSize[1]*cSize[2]*cSize[3])
  y = reshape(y,ySize[1]*ySize[2]*ySize[3])
  yTemp = reshape(yTemp,ySize[1]*ySize[2]*ySize[3])

  for s=1:N

      yTemp = 2*c-reshape(z[s,:,:,:],ySize[1]*ySize[2]*ySize[3])-y;

      # Solve TV problems
      yTemp = reshape(yTemp,ySize[1],ySize[2],ySize[3])
      tv_denoise_3d_condat!(yTemp, nhood[s,:], gamma*omega[s]*alpha/t[s])

      # Soft thresholding
      softthresh!(yTemp,gamma * beta/(N*t[s]))

      # Update step
      c = reshape(c,cSize[1],cSize[2],cSize[3])
      update!(z,s,yTemp,c,lambda)
      c = reshape(c,cSize[1]*cSize[2]*cSize[3])

  end
end

"""
This function applies soft thresholding on given data y.
"""
function softthresh!(y::Array{T,3}, threshold) where {T<:Real}

  ySize = size(y)
  y = reshape(y,ySize[1]*ySize[2]*ySize[3])

  y[:] = [sign(i)*max(i-threshold,0) for i in y]
  return
end

"""
This function performs the update step.
"""
function update!(z::Array{T,4}, s::Int64, y::Array{T,3}, c::Array{T,3}, lambda) where {T<:Real}
  dataSize = size(c)
  flatSize = dataSize[1]*dataSize[2]*dataSize[3];
  y = reshape(y,flatSize)
  c = reshape(c,flatSize)
  yTemp = Array{T}(undef, flatSize)
  zTemp = view(reshape(z,Val(2)),s,:)

  broadcast!(-,yTemp,y,c)

  #yTemp[:] = y[:] - c[:]
  rmul!(yTemp,lambda)
  #zTemp[:] = zTemp[:] + yTemp[:]
  broadcast!(+,zTemp,zTemp,yTemp)
  #z[s,:,:,:] = reshape(yTemp,dataSize)[:,:,:]
  return
end

"""
This function applies the pointwise maximum step.
"""
function pointwise_max!(c::Array{T,3},z::Array{T,4},s::Int64,y::Array{T,3},out::Array{T,3}) where {T<:Real}
  ySize = size(y)
  cSize = size(c)

  y = reshape(y,ySize[1]*ySize[2]*ySize[3])
  c = reshape(c,cSize[1]*cSize[2]*cSize[3])
  zTemp = view(reshape(z,Val(2)),s,:)
  out = reshape(out,cSize[1]*cSize[2]*cSize[3])

  yTemp = 2*c-zTemp-y;

  for j=1:length(yTemp)
    out[j] = max(yTemp[j],0)
  end
  return
end

"""
This function calculates the weighted sum c depending on weights and the proximal
mapping results z.
"""
function weighted_sum!(weights::Array{T,1},z::Array{T,4},c::Array{T,3}) where {T<:Real}
  zSize = size(z[1,:,:,:])
  c = reshape(c,zSize[1]*zSize[2]*zSize[3])
  cTemp = zeros(eltype(c),zSize[1]*zSize[2]*zSize[3])

  for i=1:length(weights)
  zTemp = view(reshape(z,Val(2)),i,:)
  cTemp[:] = cTemp[:] + weights[i]*zTemp
  end

  for i = 1:length(c)
  c[i] = cTemp[i]
  end
  return
end

"""
This function checks if the stopping criterion for the main loop is reached.
"""
function check_stop(cOld::Array{T,3},c::Array{T,3},tol::T) where {T<:Real}
  cSize = size(c)

  c = reshape(c,cSize[1]*cSize[2]*cSize[3])
  cOld = reshape(cOld,cSize[1]*cSize[2]*cSize[3])

  if (norm(cOld-c)/norm(cOld) + 1.0*10.0^-3.0) <= tol
    return true
  else
    return false
  end
end

function residuum(A::Array{T,4},c::Array{T,3},u::Array{T,1}) where {T<:Real}
  aSize = size(A)
  cSize = size(c)
  uSize = size(u)

  A = reshape(A,aSize[1],aSize[2]*aSize[3]*aSize[4])
  c = reshape(c,cSize[1]*cSize[2]*cSize[3])

  return (norm(A*c-u)^2/norm(u))
end

ext: jl


text:
versioninfo()

using BandedMatrices, Base.Test

include("test_indexing.jl")
include("test_bandedlu.jl")
include("test_bandedqr.jl")

A,B=brand(10,12,2,3),brand(10,12,3,4)

@test_approx_eq full(sparse(A)) full(A)

@test_approx_eq full(A') full(A)'
@test_approx_eq full(A.') full(A).'
@test_approx_eq full((A+im*A)') (full(A)+im*full(A))'
@test_approx_eq full((A+im*A).') (full(A)+im*full(A)).'

@test_approx_eq full(A+B) (full(A)+full(B))
@test_approx_eq full(A-B) (full(A)-full(B))

@test_approx_eq full(A.*B) (full(A).*full(B))

C,D=brand(10,10,2,3),brand(12,12,3,4)

@test_approx_eq full(C*A) full(C)*full(A)
@test_approx_eq full(A*D) full(A)*full(D)

v=rand(12)
w=rand(10)

@test_approx_eq A*v full(A)*v
@test_approx_eq A'*w full(A)'*w

A=brand(Float64,5,3,2,2)
v=rand(Complex128,3)
@test_approx_eq A*v full(A)*v

A=brand(Complex128,5,3,2,2)
v=rand(Complex128,3)
@test_approx_eq A*v full(A)*v

v=rand(Float64,3)
@test_approx_eq A*v full(A)*v

 # Test for errors in collect
A=brand(12,10,2,3)
collect(BandedMatrices.eachbandedindex(A))
for (k,j) in BandedMatrices.eachbandedindex(A)
     A[k,j]
end

A=brand(1000,1000,200,300)
B=rand(1000,1000)
@test_approx_eq A*B full(A)*B

A=brand(1200,1000,200,300)
B=rand(1000,1000)
@test_approx_eq A*B full(A)*B

A=brand(1000,1200,200,300)
B=rand(1200,1200)
@test_approx_eq A*B full(A)*B

B=brand(10,10,0,4)
@test_approx_eq B*[collect(1.0:10) collect(1.0:10)] full(B)*[collect(1.0:10) collect(1.0:10)]

@which B*[collect(1.0:10) collect(1.0:10)]

A=brand(10000,10000,2,3)
B=brand(1000,1000,200,300)
v=rand(10000)
w=rand(1000)

show(brand(10,10,3,3))
show(brand(100,80,3,2))
println()

A*v
@time for k=1:100
    A*v
end

@time for k=1:100
    B*w
end

for n in (10,50), m in (12,50), Al in (0,1,2,30), Au in (0,1,2,30)
    A=brand(n,m,Al,Au)
    kr,jr=3:10,5:12
    @test_approx_eq full(A[kr,jr]) full(A)[kr,jr]
end

# test gbmm! subpieces step by step and column by column
for n in (1,5,50), ν in (1,5,50), m in (1,5,50),
                Al in (0,1,2,30), Au in (0,1,2,30),
                Bl in (0,1,2,30), Bu in (0,1,2,30)
    A=brand(n,ν,Al,Au)
    B=brand(ν,m,Bl,Bu)
    α,β,T=0.123,0.456,Float64
    C=brand(Float64,n,m,A.l+B.l,A.u+B.u)
    a=pointer(A.data)
    b=pointer(B.data)
    c=pointer(C.data)
    sta=max(1,stride(A.data,2))
    stb=max(1,stride(B.data,2))
    stc=max(1,stride(C.data,2))

    sz=sizeof(T)

    mr=1:min(m,1+B.u)
    exC=(β*full(C)+α*full(A)*full(B))
    for j=mr
        BandedMatrices.A11_Btop_Ctop_gbmv!(α,β,
                                       n,ν,m,j,
                                       sz,
                                       a,A.l,A.u,sta,
                                       b,B.l,B.u,stb,
                                       c,C.l,C.u,stc)
   end
    @test_approx_eq C[:,mr] exC[:,mr]

    mr=1+B.u:min(1+C.u,ν+B.u,m)
    exC=(β*full(C)+α*full(A)*full(B))
    for j=mr
        BandedMatrices.Atop_Bmid_Ctop_gbmv!(α,β,
                                       n,ν,m,j,
                                       sz,
                                       a,A.l,A.u,sta,
                                       b,B.l,B.u,stb,
                                       c,C.l,C.u,stc)
   end
   if !isempty(mr)
       @test_approx_eq C[:,mr] exC[:,mr]
   end

   mr=1+C.u:min(m,ν+B.u,n+C.u)
   exC=(β*full(C)+α*full(A)*full(B))
   for j=mr
       BandedMatrices.Amid_Bmid_Cmid_gbmv!(α,β,
                                      n,ν,m,j,
                                      sz,
                                      a,A.l,A.u,sta,
                                      b,B.l,B.u,stb,
                                      c,C.l,C.u,stc)
  end
  if !isempty(mr)
      @test_approx_eq C[:,mr] exC[:,mr]
  end

  mr=ν+B.u+1:min(m,n+C.u)
  exC=(β*full(C)+α*full(A)*full(B))
  for j=mr
      BandedMatrices.Anon_Bnon_C_gbmv!(α,β,
                                     n,ν,m,j,
                                     sz,
                                     a,A.l,A.u,sta,
                                     b,B.l,B.u,stb,
                                     c,C.l,C.u,stc)
 end
 if !isempty(mr)
     @test_approx_eq C[:,mr] exC[:,mr]
 end
end

# test gbmm!

for n in (1,5,50), ν in (1,5,50), m in (1,5,50), Al in (0,1,2,30), Au in (0,1,2,30), Bl in (0,1,2,30), Bu in (0,1,2,30)
    A=brand(n,ν,Al,Au)
    B=brand(ν,m,Bl,Bu)
    α,β,T=0.123,0.456,Float64
    C=brand(Float64,n,m,A.l+B.l,A.u+B.u)

    exC=α*full(A)*full(B)+β*full(C)
    BandedMatrices.gbmm!(α,A,B,β,C)

    @test_approx_eq full(exC) full(C)
end

for n in (1,5,50), ν in (1,5,50), m in (1,5,50), Al in (0,1,2,30), Au in (0,1,2,30), Bl in (0,1,2,30), Bu in (0,1,2,30)
    A=brand(n,ν,Al,Au)
    B=brand(ν,m,Bl,Bu)
    @test_approx_eq full(A*B) full(A)*full(B)
end

A=brand(10000,10000,2,3)
B=brand(1000,1000,200,300)

@time for k=1:10
    A*A
end

@time for k=1:10
    B*B
end

#gc_enable(false)

A*v
@time for k=1:100
    A*v
end
println("Time should be   0.004986 seconds (400 allocations: 7.639 MB)")
@time for k=1:100
    B*w
end
println("Time should be   0.017208 seconds (300 allocations: 792.188 KB)")

A=brand(10000,10000,2,3)
B=brand(1000,1000,200,300)

@time for k=1:10
    A*A
end

println("Time should be   0.008115 seconds (30 allocations: 8.393 MB)")

@time for k=1:10
    B*B
end

println("Time should be   0.644119 seconds (30 allocations: 76.371 MB)")

#gc_enable(true)

## Banded Matrix of Banded Matrix

A=BandedMatrix(BandedMatrix{Float64},1,2,0,1)
A[1,1]=beye(1,1,0,1)
A[1,2]=bzeros(1,2,0,1)
A[1,2][1,1]=-1/3
A[1,2][1,2]=1/3
B=BandedMatrix(BandedMatrix{Float64},2,1,1,1)
B[1,1]=0.2beye(1,1,0,1)
B[2,1]=bzeros(2,1,1,0)
B[2,1][1,1]=-2/30
B[2,1][2,1]=1/3

@test_approx_eq (A*B)[1,1] 1/3

## Bug

# BandedMatrices.eachbandedindex(BandedMatrix(Float64,1,2,2,1))|>collect

## BigFloat

B=bzeros(BigFloat,5,5,2,3)
    for (k,j)=BandedMatrices.eachbandedindex(B)
        B[k,j]=randn()
    end

x=BigFloat[1:size(B,1)...]

@test_approx_eq full(B)*x B*x
@test_approx_eq full(B*B) full(B)*full(B)


text:
@testset "Persistence" begin

    # Intervals
    p = Interval(1.0 => 2.0)
    q = Interval(0.0 => 3.0)
    r = Interval(0.0 => 4.0)

    @test q < p
    @test q < r

    @test ComputationalHomology.birthx(p) == -1.
    @test ComputationalHomology.deathx(p) == 3.
    @test diag(p)  == Interval(1.5 => 1.5)

    #######
    # from "Computational Topology - An Introduction" by Edelsbrunner & Harer, p. 184
    #######
    flt = Filtration(SimplicialComplex, Float64)
    splxs=[ Simplex(1) => 1.0,
            Simplex(2) => 2.0,
            Simplex(3) => 3.0,
            Simplex(1,2) => 4.0,
            Simplex(2,3) => 5.0,
            Simplex(3,1) => 6.0,
            Simplex(1,2,3) => 7.0]
    for s in splxs
        push!(flt, s...)
    end
    @test size(complex(flt)) == (3,3,1)

    dgm = Dict( 0 => diagram(2.0=>4.0, 3.0=>5.0, 1.0=>Inf), 1 => diagram(6.0=>7.0))

    @testset "Intervals for TwistReduction" for (d, itrs) in diagram(TwistReduction, flt)
        for itr in itrs
            @test itr ∈ dgm[d]
        end
    end

    @testset "Intervals for Standard Reduction" for (d, itrs) in diagram(StandardReduction, flt)
        for itr in itrs
            @test itr ∈ dgm[d]
        end
    end

    @testset "Intervals for Standard Reduction" for (d, itrs) in diagram(PersistentCocycleReduction{Float64}, flt)
        for itr in itrs
            @test itr ∈ dgm[d]
        end
    end

    #######
    # from "Topology for Computing" by Zomorodian, pp.138-145
    #######
    flt = Filtration(SimplicialComplex)
    splxs=[ Simplex('a') => 0.0,
            Simplex('b') => 0.0,
            Simplex('c') => 1.0,
            Simplex('d') => 1.0,
            Simplex('a', 'b') => 1.0,
            Simplex('b', 'c') => 1.0,
            Simplex('c', 'd') => 2.0,
            Simplex('a', 'd') => 2.0,
            Simplex('a', 'c') => 3.0,
            Simplex('a', 'b', 'c') => 4.0,
            Simplex('a', 'c', 'd') => 5.0]
    for s in splxs
        push!(flt, s...)
    end

    @testset "PH Intervals " for (d, itrs) in diagram(flt, length0=true)
        titrs = d == 0 ? diagram(0.0=>Inf, 0.0=>1.0, 1.0=>1.0, 1.0=>2.0) : diagram(3.0=>4.0, 2.0=>5.0)
        for itr in itrs
            @test itr ∈ titrs
        end
    end

    # PCohomology
    gens = Dict(
        0 => [
                Chain(0, Float64) + (Simplex('a'), 1.0) + (Simplex('b'), 1.0) + (Simplex('c'), 1.0) + (Simplex('d'), 1.0),
                Chain(0, Float64) + (Simplex('b'), 1.0),
                Chain(0, Float64) + (Simplex('d'), 1.0),
            ],
        1 => [
                Chain(1, Float64) + (Simplex('a', 'd'), 1.0),
                Chain(1, Float64) + (Simplex('a', 'c'), 1.0),
            ]
    )
    @testset "PCH Intervals " for ((d1, itrs1), (d2, itrs2)) in zip(diagram(flt), diagram(PersistentCocycleReduction{Float64}, flt))
        d2 == 2 && continue
        for (i1,i2) in zip(
                sort(itrs1, by=ComputationalHomology.birthx),
                sort(itrs2, by=ComputationalHomology.birthx),
            )
            @test i1 == i2
        end
        for (ch1,ch2) in zip(itrs2, gens[d2])
            for (c1, c2) in zip(ch1.g,ch2)
                @test c1 == c2
            end
        end
    end

    #######
    # Ex.1
    #######
    flt = Filtration(SimplicialComplex)
    splxs=[ Simplex(1) => 0.0,
            Simplex(2) => 0.0,
            Simplex(3) => 0.0,
            Simplex(4) => 0.0,
            Simplex(5) => 1.0,
            Simplex(1,2) => 0.0,
            Simplex(2,3) => 0.0,
            Simplex(3,4) => 0.0,
            Simplex(4,1) => 0.0,
            Simplex(3,5) => 2.0,
            Simplex(4,5) => 3.0,
            Simplex(3,4,5) => 7.0]
    for s in splxs
        push!(flt, s...)
    end

    itr = diagram(flt, length0=true)
    @test Interval(0.0=>0.0) ∈ itr[0]
    @test Interval(1.0=>2.0) ∈ itr[0]
    @test Interval(3.0=>7.0) ∈ itr[1]

    itr = diagram(StandardReduction, flt)
    @test Interval(1.0=>2.0) ∈ itr[0]
    @test Interval(3.0=>7.0) ∈ itr[1]
    @test Interval(0.0=>Inf) ∈ itr[1]

    #######
    # Ex.2
    #######
    flt = Filtration(SimplicialComplex)
    splxs=[ Simplex(1) => 1.0,
            Simplex(2) => 2.0,
            Simplex(3) => 3.0,
            Simplex(4) => 4.0,
            Simplex(5) => 5.0,
            Simplex(1,2) => 6.0,
            Simplex(2,3) => 7.0,
            Simplex(3,4) => 8.0,
            Simplex(4,1) => 9.0,
            Simplex(3,5) => 10.0,
            Simplex(4,5) => 11.0,
            Simplex(3,4,5) => 12.0]
    for s in splxs
        push!(flt, s...)
    end
    @test size(complex(flt)) == (5,6,1)
    @test length(flt) == 12

    # Betti numbers
    ∂ = boundary(flt)
    @test sparse(∂)[5,11] == 5
    @test length(∂[9]) > 0
    R = reduce(StandardReduction, ∂)
    @test betti(flt, R, 0) == 1
    @test betti(flt, R, 1) == 1
    @test betti(flt, R, 2) == 0
    @test_throws AssertionError betti(flt, R, 3)
    @test length(R[9]) == 0
    reduce!(StandardReduction, ∂)
    @test length(∂[9]) == 0

    # PH object
    ph = persistenthomology(flt)
    @test eltype(ph) == Tuple{Float64,NTuple}
    @test length(ph) == 13

    # ph = persistenthomology(flt)
    # dgm = diagram(ph)
    # ph = persistenthomology(PersistentCocycleReduction{Float64}, flt)
    # dgm = diagram(ph)
    # generators(ph)

    # PH iterator
    @testset "Method Comparison" for (g1, g2) in zip(
            betti(homology(Int, complex(flt))),
            last(betti(persistenthomology(TwistReduction, flt))) )
        @test g1 == g2
        @test g1 == g2
    end
end


text:
export timer, TimerObservable

import Base: ==
import Base: show

"""
    timer(delay::Int)
    timer(delay::Int, period::Int)

Creation operator for the `TimerObservable`. Its like `interval`(@ref), but you can specify when should the emissions start.
`timer` returns an Observable that emits an infinite sequence of ascending integers,
with a constant interval of time, period of your choosing between those emissions.
The first emission happens after the specified `delay`.
If `period` is not specified, the output Observable emits only one value, 0.
Otherwise, it emits an infinite sequence.

# Arguments
- `delay`: the initial delay time specified as an integer denoting milliseconds to wait before emitting the first value of 0.
- `period`: the minimum period of time between emissions of the subsequent numbers (in milliseconds).

# Examples
```
using Rocket

source = timer(0, 50)

sleep(0.075)
subscription = subscribe!(source, logger())
sleep(0.105)
unsubscribe!(subscription)

close(source)
;

# output

[LogActor] Data: 2
[LogActor] Data: 3

```

See also: [`interval`](@ref), [`TimerObservable`](@ref), [`subscribe!`](@ref), [`logger`](@ref)
"""
function timer(delay::Int)
    @assert delay >= 0 "'delay' argument should be positive"
    return TimerObservable(delay, 0)
end

function timer(delay::Int, period::Int)
    @assert delay >= 0 "'delay' argument should be positive"
    @assert period >= 0 "'period' argument should be positive"
    return TimerObservable(delay, period)
end

"""
    TimerObservable

An Observable that starts emitting after an `delay` and emits
ever increasing numbers after each `period` of time thereafter.

# Parameters
- `delay`: The initial delay time specified as an integer denoting milliseconds to wait before emitting the first value of 0`.
- `period`: The minimum period of time in milliseconds between emissions of the subsequent numbers.

See also: [`timer`](@ref), [`Subscribable`](@ref)
"""
struct TimerObservable <: Subscribable{Int}
    delay  :: Int
    period :: Int
end

getdelay_ms(observable::TimerObservable)  = observable.delay
getperiod_ms(observable::TimerObservable) = observable.period

getdelay_sec(observable::TimerObservable)  = getdelay_ms(observable) / MILLISECONDS_IN_SECOND
getperiod_sec(observable::TimerObservable) = getperiod_ms(observable) / MILLISECONDS_IN_SECOND

struct TimerSubscription <: Teardown
    timer  :: Timer
end

as_teardown(::Type{ <: TimerSubscription }) = UnsubscribableTeardownLogic()

mutable struct TimerActorProps
    counter :: Int
end

struct TimerActor{A} <: Actor{ Nothing }
    actor :: A
    once  :: Bool
    props :: TimerActorProps
end

TimerActor(actor::A, once::Bool) where A = TimerActor{A}(actor, once, TimerActorProps(0))

getcounter(actor::TimerActor) = actor.props.counter
itcounter!(actor::TimerActor) = actor.props.counter = actor.props.counter + 1

isonce(actor::TimerActor) = actor.once

function on_next!(actor::TimerActor, ::Nothing)
    next!(actor.actor, getcounter(actor))
    if isonce(actor)
        complete!(actor)
    else
        itcounter!(actor)
    end
end

on_error!(actor::TimerActor, err) = error!(actor.actor, err)
on_complete!(actor::TimerActor)   = complete!(actor.actor)

function on_subscribe!(observable::TimerObservable, actor)

    tactor   = TimerActor(actor, getperiod_ms(observable) === 0)

    callback = let tactor = tactor
        (timer) -> begin
            try
                next!(tactor, nothing)
            catch err
                error!(tactor, err)
                close(timer)
            end
        end
    end

    timer = Timer(callback, getdelay_sec(observable), interval = getperiod_sec(observable))

    return TimerSubscription(timer)
end

function on_unsubscribe!(subscription::TimerSubscription)
    close(subscription.timer)
    return nothing
end

Base.:(==)(t1::TimerObservable, t2::TimerObservable) = getdelay_ms(t1) === getdelay_ms(t1) && getperiod_ms(t1) === getperiod_ms(t2)

Base.show(io::IO, observable::TimerObservable)   = print(io, "TimerObservable($(getdelay_ms(observable)), $(getperiod_ms(observable)))")
Base.show(io::IO, ::TimerSubscription)           = print(io, "TimerSubscription()")


text:
export GLMModel, LogisticModel, LinearModel, QuantileModel
export grad_scratch!

##############################
## Used to specify GLM kind ##
## (link and family pair)   ##
##############################
abstract GLMModel
#Subtypes should define the loss, predict!, and grad_scratch! functions
# if the defaults are not appropriate.

function linpred!(pr::DenseVector{Float64},
                  coefs::Vector{Float64},
                  x::DSMat{Float64};
                  offset::Vector{Float64} = emptyvector(Float64))
        A_mul_B!(pr, x, coefs)
        if !isempty(offset)
            add!(pr, offset)
        end
        pr
end

#
# default predict! for GLMModel uses identity link
#
predict!(m::GLMModel,
        pr::DenseVector{Float64},
        coefs::Vector{Float64},
        x::DSMat{Float64};
        offset::Vector{Float64} = emptyvector(Float64)) = linpred!(pr, coefs, x; offset=offset)

# grad_scratch! calculates the gradient without allocating new memory
# for the residual vector.
# resid is scratch space and should have dims equal to size(y)
# Default is appropriate for GLMs with canonical link
function grad_scratch!(m::GLMModel,
                      gr::Vector{Float64},
                      coefs::Vector{Float64},
                      x::DSMat{Float64},
                      y::Vector{Float64},
                      resid::Vector{Float64};
                      offset::Vector{Float64} = emptyvector(Float64))
    predict!(m, resid, coefs, x; offset=offset)
    subtract!(resid, y) # resid = prediction - y
    alpha = 1.0/size(x, 1)
    At_mul_B!(alpha, x, resid, 0.0, gr)
    #BLAS.gemv!('T', alpha, x, resid, 0.0, gr) #grad \propto x'resid = -x'(y - pred)
    gr
end

grad!(m::GLMModel,
    gr::Vector{Float64},
    coefs::Vector{Float64},
    x::DSMat{Float64},
    y::Vector{Float64};
    offset::Vector{Float64} = emptyvector(Float64)) = grad_scratch!(m, gr, coefs, x, y, similar(y), offset=offset)

#########
## OLS ##
#########
type LinearModel <: GLMModel end

loss(m::LinearModel, pr::DenseVector{Float64}, y::Vector{Float64}) = meansqdiff(pr, y)/2

#########################
## Logistic regression ##
#########################
type LogisticModel <: GLMModel end

function predict!(m::LogisticModel,
                           pr::DenseVector{Float64},
                           coefs::Vector{Float64},
                           x::DSMat{Float64};
                           offset::Vector{Float64} = emptyvector(Float64))
    linpred!(pr, coefs, x; offset=offset)
    map1!(LogisticFun(), pr)
    pr
end

loss(m::LogisticModel, pr::DenseVector{Float64}, y::Vector{Float64}) = mean(NBLL(), pr, y)

#########################
## Quantile regression ##
#########################
type QuantileModel <: GLMModel
    tau :: Float64
end

QuantileModel() = QuantileModel(0.5)

function grad_scratch!(m::QuantileModel,
                               gr::Vector{Float64},
                               coefs::Vector{Float64},
                               x::DSMat{Float64},
                               y::Vector{Float64},
                               scratch::Vector{Float64};
                               offset::Vector{Float64} = emptyvector(Float64))
    predict!(m, scratch, coefs, x; offset=offset)
    subtract!(scratch, y)
    map1!(Qgrad(), scratch, m.tau)
    alpha = 2.0/size(x, 1)
    BLAS.gemv!('T', alpha, x, scratch, 0.0, gr)
    gr
end

loss(m::QuantileModel, pr::DenseVector{Float64}, y::Vector{Float64}) = meanabsdiff(pr, y)


text:
function miplibdownloadfromcsv()
    files = readdir("../csv/")
    csvFiles = Vector{String}()
    for file in files
        if length(file) > 4
            if lowercase(file[(end-3):end]) == ".csv"
                push!(csvFiles, file)
            end
        end
    end

    url = "https://miplib.zib.de/WebData/instances/"
    for file in csvFiles
        dirName = file[1:(end-4)]
        mkdir("../"*dirName)
        open("../csv/"*file) do instances
            readline(instances)
            while !eof(instances)
                line = readline(instances)
                instance = strip(split(line, ',')[1], [' ', '"', '\n'])
                try
                    download(url*instance*".mps.gz", "../"*dirName*"/"*instance*".mps.gz")
                catch
                    @warn "Download from "*url*instance*".mps.gz did not work"
                end
            end
        end
    end
end

max_stars_repo_name: remi-garcia/MIPLIB-Downloader

max_issues_repo_name: remi-garcia/MIPLIB-Downloader

max_forks_repo_name: remi-garcia/MIPLIB-Downloader


text:
using PlmDCA,Test

tests = ["testdca","testcorr","testmi"]

printstyled("Running tests:\n", bold=true, color=:light_blue)

using Random
Random.seed!(345679)

res = map(tests) do t
    @eval module $(Symbol("Test_",t))
    using Random
    Random.seed!(345679)
    include($t * ".jl")
    end
    return
end
# print method ambiguities

ambiguities=Test.detect_ambiguities(PlmDCA)
printstyled("Potentially stale exports: ",bold=true, color=:light_blue)
isempty(ambiguities) ? (printstyled("\t none\n",bold=true, color=:light_green)) : (display(ambiguities))

max_stars_repo_name: giovact/PlmDCA

max_issues_repo_name: giovact/PlmDCA

max_forks_repo_name: giovact/PlmDCA


text:
##############################################################################
##
## Generalized reduced rank tests using the singularvalue decomposition
## Kleibergen Paap, Journal of Econometrics (2006)
##
## See also the Stata command ranktest
## RANKTEST: Stata module to test the rank of a matrix using the Kleibergen-Paap rk statistic
## Authors: Frank Kleibergen, Mark E Schaffer
##
## More precisely, it corresponds to the Stata command:  ranktest  (X) (Z), wald full
##############################################################################

function ranktest!(X::Matrix{Float64},
                    Z::Matrix{Float64},
                    Pi::Matrix{Float64},
                    vcov_method::CovarianceEstimator,
                    df_small::Int,
                    df_absorb::Int)

    if (size(X, 2) == 0) | (size(Z, 2) == 0)
        return NaN
    end

    # Compute theta
    Fmatrix = cholesky!(Symmetric(Z' * Z), check = false).U
    Gmatrix = cholesky!(Symmetric(X' * X), check = false).U
    theta = Fmatrix * (Gmatrix' \ Pi')'

    # compute lambda
    svddecomposition = svd(theta, full = true)
    u = svddecomposition.U
    vt = svddecomposition.Vt

    k = size(X, 2)
    l = size(Z, 2)

    u_sub = u[k:l, k:l]
    a_qq = u[1:l, k:l] * (u_sub \ sqrt(u_sub * u_sub'))

    vt_sub = vt[k,k]
    b_qq = sqrt(vt_sub * vt_sub') * (vt_sub' \ vt[1:k, k]')

    kronv = kron(b_qq, a_qq')
    lambda = kronv * vec(theta)

    # compute vhat
    if vcov_method isa Vcov.SimpleCovariance
        vlab = cholesky!(Hermitian((kronv * kronv') ./ size(X, 1)), check = false)
    else
        K = kron(Gmatrix, Fmatrix)'
        vcovmodel = Vcov.VcovData(Z, K, X, size(Z, 1) - df_small - df_absorb)
        matrix_vcov2 = Vcov.S_hat(vcovmodel, vcov_method)
        vhat = K \ (K \ matrix_vcov2)'
        vlab = cholesky!(Hermitian(kronv * vhat * kronv'), check = false)
    end
    r_kp = lambda' * (vlab \ lambda)
    return r_kp[1]
end


text:
recurs_sym_type(ex::Any) =
  (ex==None || typeof(ex)==Symbol || length(ex.args)==1) ? eval(ex) : Expr(ex.head, ex.args[1], recurs_sym_type(ex.args[2]))
macro c(ret_type, func, arg_types, lib)
  local _arg_types = Expr(:tuple, [recurs_sym_type(a) for a in arg_types.args]...)
  local _ret_type = recurs_sym_type(ret_type)
  local _args_in = Any[ symbol(string('a',x)) for x in 1:length(_arg_types.args) ]
  local _lib = eval(lib)
  quote
    $(esc(func))($(_args_in...)) = ccall( ($(string(func)), $(Expr(:quote, _lib)) ), $_ret_type, $_arg_types, $(_args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  real_t = recurs_sym_type(real_t)
  quote
    typealias $fake_t $real_t
  end
end

# header: /usr/local/include/cvode/cvode_band.h
@ctypedef CVDlsDenseJacFn Ptr{:Void}
@ctypedef CVDlsBandJacFn Ptr{:Void}
@c Int32 CVDlsSetDenseJacFn (Ptr{:None},:CVDlsDenseJacFn) shlib
@c Int32 CVDlsSetBandJacFn (Ptr{:None},:CVDlsBandJacFn) shlib
@c Int32 CVDlsGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 CVDlsGetNumJacEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVDlsGetNumRhsEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVDlsGetLastFlag (Ptr{:None},Ptr{:Int64}) shlib
@c Ptr{:Uint8} CVDlsGetReturnFlagName (:Int64,) shlib
@c Int32 CVBand (Ptr{:None},:Int64,:Int64,:Int64) shlib

# header: /usr/local/include/cvode/cvode_bandpre.h
@c Int32 CVBandPrecInit (Ptr{:None},:Int64,:Int64,:Int64) shlib
@c Int32 CVBandPrecGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 CVBandPrecGetNumRhsEvals (Ptr{:None},Ptr{:Int64}) shlib

# header: /usr/local/include/cvode/cvode_bbdpre.h
@ctypedef CVLocalFn Ptr{:Void}
@ctypedef CVCommFn Ptr{:Void}
@c Int32 CVBBDPrecInit (Ptr{:None},:Int64,:Int64,:Int64,:Int64,:Int64,:realtype,:CVLocalFn,:CVCommFn) shlib
@c Int32 CVBBDPrecReInit (Ptr{:None},:Int64,:Int64,:realtype) shlib
@c Int32 CVBBDPrecGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 CVBBDPrecGetNumGfnEvals (Ptr{:None},Ptr{:Int64}) shlib

# header: /usr/local/include/cvode/cvode_dense.h
@c Int32 CVDense (Ptr{:None},:Int64) shlib

# header: /usr/local/include/cvode/cvode_diag.h
@c Int32 CVDiag (Ptr{:None},) shlib
@c Int32 CVDiagGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 CVDiagGetNumRhsEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVDiagGetLastFlag (Ptr{:None},Ptr{:Int64}) shlib
@c Ptr{:Uint8} CVDiagGetReturnFlagName (:Int64,) shlib

# header: /usr/local/include/cvode/cvode_direct.h

# header: /usr/local/include/cvode/cvode.h
@ctypedef ptrdiff_t Void
@ctypedef size_t Void
@ctypedef wchar_t Int32
@ctypedef __u_char Uint8
@ctypedef __u_short Uint16
@ctypedef __u_int Uint32
@ctypedef __u_long Uint64
@ctypedef __int8_t Uint8
@ctypedef __uint8_t Uint8
@ctypedef __int16_t Int16
@ctypedef __uint16_t Uint16
@ctypedef __int32_t Int32
@ctypedef __uint32_t Uint32
@ctypedef __int64_t Int64
@ctypedef __uint64_t Uint64
@ctypedef __quad_t Int64
@ctypedef __u_quad_t Uint64
@ctypedef __dev_t Uint64
@ctypedef __uid_t Uint32
@ctypedef __gid_t Uint32
@ctypedef __ino_t Uint64
@ctypedef __ino64_t Uint64
@ctypedef __mode_t Uint32
@ctypedef __nlink_t Uint64
@ctypedef __off_t Int64
@ctypedef __off64_t Int64
@ctypedef __pid_t Int32
@ctypedef __fsid_t Void
@ctypedef __clock_t Int64
@ctypedef __rlim_t Uint64
@ctypedef __rlim64_t Uint64
@ctypedef __id_t Uint32
@ctypedef __time_t Int64
@ctypedef __useconds_t Uint32
@ctypedef __suseconds_t Int64
@ctypedef __daddr_t Int32
@ctypedef __swblk_t Int64
@ctypedef __key_t Int32
@ctypedef __clockid_t Int32
@ctypedef __timer_t Ptr{:None}
@ctypedef __blksize_t Int64
@ctypedef __blkcnt_t Int64
@ctypedef __blkcnt64_t Int64
@ctypedef __fsblkcnt_t Uint64
@ctypedef __fsblkcnt64_t Uint64
@ctypedef __fsfilcnt_t Uint64
@ctypedef __fsfilcnt64_t Uint64
@ctypedef __ssize_t Int64
@ctypedef __loff_t __off64_t
@ctypedef __qaddr_t Ptr{:__quad_t}
@ctypedef __caddr_t Ptr{:Uint8}
@ctypedef __intptr_t Int64
@ctypedef __socklen_t Uint32
@ctypedef FILE Void
@ctypedef __FILE Void
@ctypedef __mbstate_t Void
@ctypedef _G_fpos_t Void
@ctypedef _G_fpos64_t Void
@ctypedef _G_int16_t Int16
@ctypedef _G_int32_t Int32
@ctypedef _G_uint16_t Uint16
@ctypedef _G_uint32_t Uint32
@ctypedef va_list __builtin_va_list
@ctypedef __gnuc_va_list __builtin_va_list
@ctypedef _IO_lock_t None
# enum __codecvt_result
@ctypedef __codecvt_result Uint32
const __codecvt_ok = 0
const __codecvt_partial = 1
const __codecvt_error = 2
const __codecvt_noconv = 3
# end
@ctypedef _IO_FILE Void
@ctypedef __io_read_fn Void
@ctypedef __io_write_fn Void
@ctypedef __io_seek_fn Void
@ctypedef __io_close_fn Void
@c Int32 __underflow (Ptr{:_IO_FILE},) shlib
@c Int32 __uflow (Ptr{:_IO_FILE},) shlib
@c Int32 __overflow (Ptr{:_IO_FILE},:Int32) shlib
@c Int32 _IO_getc (Ptr{:_IO_FILE},) shlib
@c Int32 _IO_putc (:Int32,Ptr{:_IO_FILE}) shlib
@c Int32 _IO_feof (Ptr{:_IO_FILE},) shlib
@c Int32 _IO_ferror (Ptr{:_IO_FILE},) shlib
@c Int32 _IO_peekc_locked (Ptr{:_IO_FILE},) shlib
@c None _IO_flockfile (Ptr{:_IO_FILE},) shlib
@c None _IO_funlockfile (Ptr{:_IO_FILE},) shlib
@c Int32 _IO_ftrylockfile (Ptr{:_IO_FILE},) shlib
@c Int32 _IO_vfscanf (Ptr{:_IO_FILE},Ptr{:Uint8},Ptr{:__va_list_tag},Ptr{:Int32}) shlib
@c Int32 _IO_vfprintf (Ptr{:_IO_FILE},Ptr{:Uint8},Ptr{:__va_list_tag}) shlib
@c __ssize_t _IO_padn (Ptr{:_IO_FILE},:Int32,:__ssize_t) shlib
@c size_t _IO_sgetn (Ptr{:_IO_FILE},Ptr{:None},:size_t) shlib
@c __off64_t _IO_seekoff (Ptr{:_IO_FILE},:__off64_t,:Int32,:Int32) shlib
@c __off64_t _IO_seekpos (Ptr{:_IO_FILE},:__off64_t,:Int32) shlib
@c None _IO_free_backup_area (Ptr{:_IO_FILE},) shlib
@ctypedef off_t __off_t
@ctypedef ssize_t __ssize_t
@ctypedef fpos_t _G_fpos_t
@c Int32 remove (Ptr{:Uint8},) shlib
@c Int32 rename (Ptr{:Uint8},Ptr{:Uint8}) shlib
@c Int32 renameat (:Int32,Ptr{:Uint8},:Int32,Ptr{:Uint8}) shlib
@c Ptr{:FILE} tmpfile () shlib
@c Ptr{:Uint8} tmpnam (Ptr{:Uint8},) shlib
@c Ptr{:Uint8} tmpnam_r (Ptr{:Uint8},) shlib
@c Ptr{:Uint8} tempnam (Ptr{:Uint8},Ptr{:Uint8}) shlib
@c Int32 fclose (Ptr{:FILE},) shlib
@c Int32 fflush (Ptr{:FILE},) shlib
@c Int32 fflush_unlocked (Ptr{:FILE},) shlib
@c Ptr{:FILE} fopen (Ptr{:Uint8},Ptr{:Uint8}) shlib
@c Ptr{:FILE} freopen (Ptr{:Uint8},Ptr{:Uint8},Ptr{:FILE}) shlib
@c Ptr{:FILE} fdopen (:Int32,Ptr{:Uint8}) shlib
@c Ptr{:FILE} fmemopen (Ptr{:None},:size_t,Ptr{:Uint8}) shlib
@c Ptr{:FILE} open_memstream (Ptr{Ptr{:Uint8}},Ptr{:size_t}) shlib
@c None setbuf (Ptr{:FILE},Ptr{:Uint8}) shlib
@c Int32 setvbuf (Ptr{:FILE},Ptr{:Uint8},:Int32,:size_t) shlib
@c None setbuffer (Ptr{:FILE},Ptr{:Uint8},:size_t) shlib
@c None setlinebuf (Ptr{:FILE},) shlib
@c Int32 fprintf (Ptr{:FILE},Ptr{:Uint8}) shlib
@c Int32 printf (Ptr{:Uint8},) shlib
@c Int32 sprintf (Ptr{:Uint8},Ptr{:Uint8}) shlib
@c Int32 vfprintf (Ptr{:FILE},Ptr{:Uint8},Ptr{:__va_list_tag}) shlib
@c Int32 vprintf (Ptr{:Uint8},Ptr{:__va_list_tag}) shlib
@c Int32 vsprintf (Ptr{:Uint8},Ptr{:Uint8},Ptr{:__va_list_tag}) shlib
@c Int32 snprintf (Ptr{:Uint8},:size_t,Ptr{:Uint8}) shlib
@c Int32 vsnprintf (Ptr{:Uint8},:size_t,Ptr{:Uint8},Ptr{:__va_list_tag}) shlib
@c Int32 vdprintf (:Int32,Ptr{:Uint8},Ptr{:__va_list_tag}) shlib
@c Int32 dprintf (:Int32,Ptr{:Uint8}) shlib
@c Int32 fscanf (Ptr{:FILE},Ptr{:Uint8}) shlib
@c Int32 scanf (Ptr{:Uint8},) shlib
@c Int32 sscanf (Ptr{:Uint8},Ptr{:Uint8}) shlib
@c Int32 vfscanf (Ptr{:FILE},Ptr{:Uint8},Ptr{:__va_list_tag}) shlib
@c Int32 vscanf (Ptr{:Uint8},Ptr{:__va_list_tag}) shlib
@c Int32 vsscanf (Ptr{:Uint8},Ptr{:Uint8},Ptr{:__va_list_tag}) shlib
@c Int32 fgetc (Ptr{:FILE},) shlib
@c Int32 getc (Ptr{:FILE},) shlib
@c Int32 getchar () shlib
@c Int32 getc_unlocked (Ptr{:FILE},) shlib
@c Int32 getchar_unlocked () shlib
@c Int32 fgetc_unlocked (Ptr{:FILE},) shlib
@c Int32 fputc (:Int32,Ptr{:FILE}) shlib
@c Int32 putc (:Int32,Ptr{:FILE}) shlib
@c Int32 putchar (:Int32,) shlib
@c Int32 fputc_unlocked (:Int32,Ptr{:FILE}) shlib
@c Int32 putc_unlocked (:Int32,Ptr{:FILE}) shlib
@c Int32 putchar_unlocked (:Int32,) shlib
@c Int32 getw (Ptr{:FILE},) shlib
@c Int32 putw (:Int32,Ptr{:FILE}) shlib
@c Ptr{:Uint8} fgets (Ptr{:Uint8},:Int32,Ptr{:FILE}) shlib
@c Ptr{:Uint8} gets (Ptr{:Uint8},) shlib
@c __ssize_t __getdelim (Ptr{Ptr{:Uint8}},Ptr{:size_t},:Int32,Ptr{:FILE}) shlib
@c __ssize_t getdelim (Ptr{Ptr{:Uint8}},Ptr{:size_t},:Int32,Ptr{:FILE}) shlib
@c __ssize_t getline (Ptr{Ptr{:Uint8}},Ptr{:size_t},Ptr{:FILE}) shlib
@c Int32 fputs (Ptr{:Uint8},Ptr{:FILE}) shlib
@c Int32 puts (Ptr{:Uint8},) shlib
@c Int32 ungetc (:Int32,Ptr{:FILE}) shlib
@c size_t fread (Ptr{:None},:size_t,:size_t,Ptr{:FILE}) shlib
@c size_t fwrite (Ptr{:None},:size_t,:size_t,Ptr{:FILE}) shlib
@c size_t fread_unlocked (Ptr{:None},:size_t,:size_t,Ptr{:FILE}) shlib
@c size_t fwrite_unlocked (Ptr{:None},:size_t,:size_t,Ptr{:FILE}) shlib
@c Int32 fseek (Ptr{:FILE},:Int64,:Int32) shlib
@c Int64 ftell (Ptr{:FILE},) shlib
@c None rewind (Ptr{:FILE},) shlib
@c Int32 fseeko (Ptr{:FILE},:__off_t,:Int32) shlib
@c __off_t ftello (Ptr{:FILE},) shlib
@c Int32 fgetpos (Ptr{:FILE},Ptr{:fpos_t}) shlib
@c Int32 fsetpos (Ptr{:FILE},Ptr{:fpos_t}) shlib
@c None clearerr (Ptr{:FILE},) shlib
@c Int32 feof (Ptr{:FILE},) shlib
@c Int32 ferror (Ptr{:FILE},) shlib
@c None clearerr_unlocked (Ptr{:FILE},) shlib
@c Int32 feof_unlocked (Ptr{:FILE},) shlib
@c Int32 ferror_unlocked (Ptr{:FILE},) shlib
@c None perror (Ptr{:Uint8},) shlib
@c Int32 fileno (Ptr{:FILE},) shlib
@c Int32 fileno_unlocked (Ptr{:FILE},) shlib
@c Ptr{:FILE} popen (Ptr{:Uint8},Ptr{:Uint8}) shlib
@c Int32 pclose (Ptr{:FILE},) shlib
@c Ptr{:Uint8} ctermid (Ptr{:Uint8},) shlib
@c None flockfile (Ptr{:FILE},) shlib
@c Int32 ftrylockfile (Ptr{:FILE},) shlib
@c None funlockfile (Ptr{:FILE},) shlib
@ctypedef CVRhsFn Ptr{:Void}
@ctypedef CVRootFn Ptr{:Void}
@ctypedef CVEwtFn Ptr{:Void}
@ctypedef CVErrHandlerFn Ptr{:Void}
@c Ptr{:None} CVodeCreate (:Int32,:Int32) shlib
@c Int32 CVodeSetErrHandlerFn (Ptr{:None},:CVErrHandlerFn,Ptr{:None}) shlib
@c Int32 CVodeSetErrFile (Ptr{:None},Ptr{:FILE}) shlib
@c Int32 CVodeSetUserData (Ptr{:None},Ptr{:None}) shlib
@c Int32 CVodeSetMaxOrd (Ptr{:None},:Int32) shlib
@c Int32 CVodeSetMaxNumSteps (Ptr{:None},:Int64) shlib
@c Int32 CVodeSetMaxHnilWarns (Ptr{:None},:Int32) shlib
@c Int32 CVodeSetStabLimDet (Ptr{:None},:Int32) shlib
@c Int32 CVodeSetInitStep (Ptr{:None},:realtype) shlib
@c Int32 CVodeSetMinStep (Ptr{:None},:realtype) shlib
@c Int32 CVodeSetMaxStep (Ptr{:None},:realtype) shlib
@c Int32 CVodeSetStopTime (Ptr{:None},:realtype) shlib
@c Int32 CVodeSetMaxErrTestFails (Ptr{:None},:Int32) shlib
@c Int32 CVodeSetMaxNonlinIters (Ptr{:None},:Int32) shlib
@c Int32 CVodeSetMaxConvFails (Ptr{:None},:Int32) shlib
@c Int32 CVodeSetNonlinConvCoef (Ptr{:None},:realtype) shlib
@c Int32 CVodeSetIterType (Ptr{:None},:Int32) shlib
@c Int32 CVodeSetRootDirection (Ptr{:None},Ptr{:Int32}) shlib
@c Int32 CVodeSetNoInactiveRootWarn (Ptr{:None},) shlib
@c Int32 CVodeInit (Ptr{:None},:CVRhsFn,:realtype,:N_Vector) shlib
@c Int32 CVodeReInit (Ptr{:None},:realtype,:N_Vector) shlib
@c Int32 CVodeSStolerances (Ptr{:None},:realtype,:realtype) shlib
@c Int32 CVodeSVtolerances (Ptr{:None},:realtype,:N_Vector) shlib
@c Int32 CVodeWFtolerances (Ptr{:None},:CVEwtFn) shlib
@c Int32 CVodeRootInit (Ptr{:None},:Int32,:CVRootFn) shlib
@c Int32 CVode (Ptr{:None},:realtype,:N_Vector,Ptr{:realtype},:Int32) shlib
@c Int32 CVodeGetDky (Ptr{:None},:realtype,:Int32,:N_Vector) shlib
@c Int32 CVodeGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 CVodeGetNumSteps (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVodeGetNumRhsEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVodeGetNumLinSolvSetups (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVodeGetNumErrTestFails (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVodeGetLastOrder (Ptr{:None},Ptr{:Int32}) shlib
@c Int32 CVodeGetCurrentOrder (Ptr{:None},Ptr{:Int32}) shlib
@c Int32 CVodeGetNumStabLimOrderReds (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVodeGetActualInitStep (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 CVodeGetLastStep (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 CVodeGetCurrentStep (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 CVodeGetCurrentTime (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 CVodeGetTolScaleFactor (Ptr{:None},Ptr{:realtype}) shlib
@c Int32 CVodeGetErrWeights (Ptr{:None},:N_Vector) shlib
@c Int32 CVodeGetEstLocalErrors (Ptr{:None},:N_Vector) shlib
@c Int32 CVodeGetNumGEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVodeGetRootInfo (Ptr{:None},Ptr{:Int32}) shlib
@c Int32 CVodeGetIntegratorStats (Ptr{:None},Ptr{:Int64},Ptr{:Int64},Ptr{:Int64},Ptr{:Int64},Ptr{:Int32},Ptr{:Int32},Ptr{:realtype},Ptr{:realtype},Ptr{:realtype},Ptr{:realtype}) shlib
@c Int32 CVodeGetNumNonlinSolvIters (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVodeGetNumNonlinSolvConvFails (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVodeGetNonlinSolvStats (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Ptr{:Uint8} CVodeGetReturnFlagName (:Int64,) shlib
@c None CVodeFree (Ptr{Ptr{:None}},) shlib

# header: /usr/local/include/cvode/cvode_impl.h
@ctypedef CVodeMem Ptr{:Void}
@c Int32 CVEwtSet (:N_Vector,:N_Vector,Ptr{:None}) shlib
@c None CVProcessError (:CVodeMem,:Int32,Ptr{:Uint8},Ptr{:Uint8},Ptr{:Uint8}) shlib
@c None CVErrHandler (:Int32,Ptr{:Uint8},Ptr{:Uint8},Ptr{:Uint8},Ptr{:None}) shlib

# header: /usr/local/include/cvode/cvode_spbcgs.h
@ctypedef CVSpilsPrecSetupFn Ptr{:Void}
@ctypedef CVSpilsPrecSolveFn Ptr{:Void}
@ctypedef CVSpilsJacTimesVecFn Ptr{:Void}
@c Int32 CVSpilsSetPrecType (Ptr{:None},:Int32) shlib
@c Int32 CVSpilsSetGSType (Ptr{:None},:Int32) shlib
@c Int32 CVSpilsSetMaxl (Ptr{:None},:Int32) shlib
@c Int32 CVSpilsSetEpsLin (Ptr{:None},:realtype) shlib
@c Int32 CVSpilsSetPreconditioner (Ptr{:None},:CVSpilsPrecSetupFn,:CVSpilsPrecSolveFn) shlib
@c Int32 CVSpilsSetJacTimesVecFn (Ptr{:None},:CVSpilsJacTimesVecFn) shlib
@c Int32 CVSpilsGetWorkSpace (Ptr{:None},Ptr{:Int64},Ptr{:Int64}) shlib
@c Int32 CVSpilsGetNumPrecEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVSpilsGetNumPrecSolves (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVSpilsGetNumLinIters (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVSpilsGetNumConvFails (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVSpilsGetNumJtimesEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVSpilsGetNumRhsEvals (Ptr{:None},Ptr{:Int64}) shlib
@c Int32 CVSpilsGetLastFlag (Ptr{:None},Ptr{:Int64}) shlib
@c Ptr{:Uint8} CVSpilsGetReturnFlagName (:Int64,) shlib
@c Int32 CVSpbcg (Ptr{:None},:Int32,:Int32) shlib

# header: /usr/local/include/cvode/cvode_spgmr.h
@c Int32 CVSpgmr (Ptr{:None},:Int32,:Int32) shlib

# header: /usr/local/include/cvode/cvode_spils.h

# header: /usr/local/include/cvode/cvode_sptfqmr.h
@c Int32 CVSptfqmr (Ptr{:None},:Int32,:Int32) shlib


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#
"""
Proper cone definitions, oracles, and utilities.
"""
module Cones

using DocStringExtensions
using LinearAlgebra
import LinearAlgebra.copytri!
import LinearAlgebra.HermOrSym
import LinearAlgebra.BlasReal
import PolynomialRoots
using SparseArrays
import Hypatia.RealOrComplex
import Hypatia.outer_prod!
import Hypatia.update_eigen!
import Hypatia.spectral_outer!
import Hypatia.posdef_fact_copy!
import Hypatia.inv_fact!

include("arrayutilities.jl")

"""
$(TYPEDEF)

A proper cone.
"""
abstract type Cone{T <: Real} end

"""
$(SIGNATURES)

The real vector dimension of the cone.
"""
dimension(cone::Cone)::Int = cone.dim

"""
$(SIGNATURES)

The barrier parameter ``\\nu`` of the cone.
"""
get_nu(cone::Cone)::Real = cone.nu

"""
$(SIGNATURES)

Set the array equal to the initial interior point for the cone.
"""
function set_initial_point!(arr::AbstractVector, cone::Cone) end

"""
$(SIGNATURES)

Returns true if and only if the currently-loaded primal point is strictly
feasible for the cone.
"""
is_feas(cone::Cone)::Bool = (cone.feas_updated ? cone.is_feas : update_feas(cone))

"""
$(SIGNATURES)

Returns false only if the currently-loaded dual point is outside the interior of
the cone's dual cone.
"""
is_dual_feas(cone::Cone)::Bool = true

"""
$(SIGNATURES)

The gradient of the cone's barrier function at the currently-loaded primal point.
"""
grad(cone::Cone) = (cone.grad_updated ? cone.grad : update_grad(cone))

"""
$(SIGNATURES)

The Hessian (symmetric positive definite) of the cone's barrier function at the
currently-loaded primal point.
"""
function hess(cone::Cone)
    cone.hess_updated && return cone.hess
    return update_hess(cone)
end

"""
$(SIGNATURES)

The inverse Hessian (symmetric positive definite) of the cone's barrier function
at the currently-loaded primal point.
"""
function inv_hess(cone::Cone)
    cone.inv_hess_updated && return cone.inv_hess
    return update_inv_hess(cone)
end

"""
$(SIGNATURES)

Compute the product of the Hessian of the cone's barrier function at the
currently-loaded primal point with a vector or array, in-place.
"""
function hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::Cone)
    cone.hess_updated || update_hess(cone)
    mul!(prod, cone.hess, arr)
    return prod
end

"""
$(SIGNATURES)

Compute the product of the inverse Hessian of the cone's barrier function at the
currently-loaded primal point with a vector or array, in-place.
"""
function inv_hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::Cone)
    update_hess_fact(cone)
    # TODO try equilibration, iterative refinement etc like posvx/sysvx
    ldiv!(prod, cone.hess_fact, arr)
    return prod
end

"""
$(SIGNATURES)

Returns true if and only if the oracle for the third-order directional
derivative oracle [`dder3`](@ref) can be computed.
"""
use_dder3(::Cone)::Bool = true

"""
$(SIGNATURES)

Compute the third-order directional derivative, in the direction `dir`, the
cone's barrier function at the currently-loaded primal point.
"""
function dder3(cone::Cone, dir::AbstractVector) end

# other oracles and helpers

use_dual_barrier(cone::Cone)::Bool = cone.use_dual_barrier

function setup_data!(cone::Cone{T}) where {T <: Real}
    reset_data(cone)
    dim = dimension(cone)
    cone.point = zeros(T, dim)
    cone.dual_point = zeros(T, dim)
    cone.grad = zeros(T, dim)
    if hasproperty(cone, :dder3)
        cone.dder3 = zeros(T, dim)
    end
    cone.vec1 = zeros(T, dim)
    cone.vec2 = zeros(T, dim)
    setup_extra_data!(cone)
    return cone
end

setup_extra_data!(cone::Cone) = nothing

load_point(
    cone::Cone{T},
    point::AbstractVector{T},
    scal::T,
    ) where {T <: Real} = (@. cone.point = scal * point)

load_point(
    cone::Cone,
    point::AbstractVector,
    ) = copyto!(cone.point, point)

load_dual_point(
    cone::Cone,
    point::AbstractVector,
    ) = copyto!(cone.dual_point, point)

function alloc_hess!(cone::Cone{T}) where {T <: Real}
    dim = dimension(cone)
    cone.hess = Symmetric(zeros(T, dim, dim), :U)
    return
end

function alloc_inv_hess!(cone::Cone{T}) where {T <: Real}
    dim = dimension(cone)
    cone.inv_hess = Symmetric(zeros(T, dim, dim), :U)
    return
end

reset_data(cone::Cone) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_fact_updated = false)

# decide whether to use sqrt oracles
function use_sqrt_hess_oracles(arr_dim::Int, cone::Cone)
    if !cone.hess_fact_updated
        (arr_dim < dimension(cone)) && return false # array is small
        update_hess_fact(cone) || return false
    end
    return (cone.hess_fact isa Cholesky)
end

# only use if use_sqrt_hess_oracles is true
function sqrt_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::Cone,
    )
    @assert cone.hess_fact_updated
    mul!(prod, cone.hess_fact.U, arr)
    return prod
end

# only use if use_sqrt_hess_oracles is true
function inv_sqrt_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::Cone,
    )
    @assert cone.hess_fact_updated
    # TODO try equilibration, iterative refinement etc like posvx/sysvx
    ldiv!(prod, cone.hess_fact.U', arr)
    return prod
end

update_hess_aux(cone::Cone) = nothing

function update_use_hess_prod_slow(cone::Cone{T}) where {T <: Real}
    cone.hess_updated || update_hess(cone)
    @assert cone.hess_updated
    rel_viol = abs(1 - dot(cone.point, cone.hess, cone.point) / get_nu(cone))
    # TODO tune
    cone.use_hess_prod_slow = (rel_viol > dimension(cone) * sqrt(eps(T)))
    # cone.use_hess_prod_slow && println("switching to slow hess prod")
    cone.use_hess_prod_slow_updated = true
    return
end

hess_prod_slow!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::Cone,
    ) = hess_prod!(prod, arr, cone)

function update_hess_fact(cone::Cone{T}) where {T <: Real}
    cone.hess_fact_updated && return true
    cone.hess_updated || update_hess(cone)
    if !isdefined(cone, :hess_fact_mat)
        cone.hess_fact_mat = zero(cone.hess)
    end

    # do not modify the hessian during recovery
    cone.hess_fact = posdef_fact_copy!(cone.hess_fact_mat, cone.hess, false)

    cone.hess_fact_updated = true
    return issuccess(cone.hess_fact)
end

function update_inv_hess(cone::Cone)
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    update_hess_fact(cone)
    inv_fact!(cone.inv_hess.data, cone.hess_fact)
    cone.inv_hess_updated = true
    return cone.inv_hess
end

# number of nonzeros in the Hessian and inverse
hess_nz_count(cone::Cone) = dimension(cone) ^ 2
hess_nz_count_tril(cone::Cone) = svec_length(dimension(cone))
inv_hess_nz_count(cone::Cone) = dimension(cone) ^ 2
inv_hess_nz_count_tril(cone::Cone) = svec_length(dimension(cone))
# row indices of nonzero elements in column j
hess_nz_idxs_col(cone::Cone, j::Int) = 1:dimension(cone)
hess_nz_idxs_col_tril(cone::Cone, j::Int) = j:dimension(cone)
inv_hess_nz_idxs_col(cone::Cone, j::Int) = 1:dimension(cone)
inv_hess_nz_idxs_col_tril(cone::Cone, j::Int) = j:dimension(cone)

# check numerics of some oracles used in proximity check TODO tune
function check_numerics(
    cone::Cone{T},
    gtol::T = sqrt(sqrt(eps(T))),
    Htol::T = 10sqrt(gtol),
    ) where {T <: Real}
    g = grad(cone)
    dim = length(g)
    nu = get_nu(cone)

    # grad check
    (abs(1 + dot(g, cone.point) / nu) > gtol * dim) && return false

    # inv hess check
    Hig = inv_hess_prod!(cone.vec1, g, cone)
    (abs(1 - dot(Hig, g) / nu) > Htol * dim) && return false

    return true
end

# compute central path proximity for a cone; if using max proximity, proximity
# is computed differently if cone is not primitive, eg nonnegative cone
function get_proximity(
    cone::Cone{T},
    rtmu::T,
    ::Bool, # use sum proximity
    negtol::T = sqrt(eps(T)),
    ) where {T <: Real}
    g = grad(cone)
    vec1 = cone.vec1
    vec2 = cone.vec2

    @. vec1 = cone.dual_point + rtmu * g
    inv_hess_prod!(vec2, vec1, cone)
    prox_sqr = dot(vec2, vec1)
    (prox_sqr < -negtol * length(g)) && return T(NaN) # should be positive
    return sqrt(abs(prox_sqr)) / rtmu
end

include("nonnegative.jl")
include("possemideftri.jl")
include("doublynonnegativetri.jl")
include("possemideftrisparse/possemideftrisparse.jl")
include("linmatrixineq.jl")
include("epinorminf.jl")
include("epinormeucl.jl")
include("epipersquare.jl")
include("epinormspectral.jl")
include("matrixepipersquare.jl")
include("generalizedpower.jl")
include("hypopowermean.jl")
include("hypogeomean.jl")
include("hyporootdettri.jl")
include("hypoperlog.jl")
include("hypoperlogdettri.jl")
include("epipersepspectral/epipersepspectral.jl")
include("epirelentropy.jl")
include("epitrrelentropytri.jl")
include("wsosinterpnonnegative.jl")
include("wsosinterppossemideftri.jl")
include("wsosinterpepinormone.jl")
include("wsosinterpepinormeucl.jl")

end


text:
@testset "test dataset sample termination" begin
    @testset "5-bus case file, number of samples" begin
		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path="data")
        @test size(results["inputs"]["pd"],1) == K
		@test size(results["inputs"]["pd"],2) == 3
		@test size(results["outputs"]["p_gen"],1) == K
		@test size(results["outputs"]["p_gen"],2) == 5
		@test size(results["duals"]["v_max"],1) == K
		@test size(results["duals"]["v_max"],2) == 5
    end

    @testset "30-bus case file, number of samples" begin
		net_file = "pglib_opf_case30_ieee.m"
		K = 10
		results = create_samples(net_file, K, net_path="data")
        @test size(results["inputs"]["pd"],1) == K
		@test size(results["inputs"]["pd"],2) == 21
		@test size(results["outputs"]["p_gen"],1) == K
		@test size(results["outputs"]["p_gen"],2) == 6
		@test size(results["duals"]["v_max"],1) == K
		@test size(results["duals"]["v_max"],2) == 30
    end
end

@testset "test dataset sample variable specification" begin
    @testset "5-bus case file, default variables" begin
		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path="data")
        @test length(results["inputs"]) == length(OPFLearn.DEFAULT_INPUTS)
		@test length(results["outputs"]) == length(OPFLearn.DEFAULT_OUTPUTS)
		@test length(results["duals"]) == length(OPFLearn.DEFAULT_DUALS) + 3  # + k, uas, and variances
    end

    @testset "5-bus case file, specified variables" begin
		input_vars = ["pd", "qd"]
		output_vars = ["p_gen", "vm_gen"]
		dual_vars = ["v_min", "qg_min", "pto_max"]

		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path="data", input_vars=input_vars,
							     output_vars=output_vars, dual_vars=dual_vars)
        @test length(results["inputs"]) == length(input_vars)
		@test length(results["outputs"]) == length(output_vars)
		@test length(results["duals"]) == length(dual_vars) + 3  # Duals includes k, uas, and variances
		@test all([x in keys(results["inputs"]) for x in input_vars])
		@test all([x in keys(results["outputs"]) for x in output_vars])
		@test all([x in keys(results["duals"]) for x in dual_vars])
    end

    @testset "14-bus case file, specified variables" begin
		input_vars = ["pd"]
		output_vars = ["q_gen"]
		dual_vars = ["v_min", "qg_min", "pto_max"]

		net_file = "pglib_opf_case14_ieee.m"
		K = 10
		results = create_samples(net_file, K, net_path="data", input_vars=input_vars,
								 output_vars=output_vars, dual_vars=dual_vars)
        @test length(results["inputs"]) == length(input_vars)
		@test length(results["outputs"]) == length(output_vars)
		@test length(results["duals"]) == length(dual_vars) + 3  # Duals includes k, uas, and variances
		@test all([x in keys(results["inputs"]) for x in input_vars])
		@test all([x in keys(results["outputs"]) for x in output_vars])
		@test all([x in keys(results["duals"]) for x in dual_vars])
    end

	@testset "14-bus case file, empty duals" begin
		input_vars = ["pd", "qd"]
		output_vars = ["p_gen", "vm_gen"]
		dual_vars = []

		net_file = "pglib_opf_case14_ieee.m"
		K = 10
		results = create_samples(net_file, K, net_path="data", input_vars=input_vars,
							     output_vars=output_vars, dual_vars=dual_vars)
        @test length(results["inputs"]) == length(input_vars)
		@test length(results["outputs"]) == length(output_vars)
		@test length(results["duals"]) == length(dual_vars) + 3  # Duals includes k, uas, and variances
		@test all([x in keys(results["inputs"]) for x in input_vars])
		@test all([x in keys(results["outputs"]) for x in output_vars])
		@test all([x in keys(results["duals"]) for x in dual_vars])
    end
end

@testset "test additional dataset results" begin
    @testset "save polytope" begin
		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path="data", save_certs=true)
        @test size(results["inputs"]["pd"],1) == K
		@test size(results["inputs"]["pd"],2) == 3
		@test size(results["outputs"]["p_gen"],1) == K
		@test size(results["outputs"]["p_gen"],2) == 5
		@test size(results["duals"]["v_max"],1) == K
		@test size(results["duals"]["v_max"],2) == 5
		@test haskey(results, "polytope")
	end

    @testset "save stats" begin
		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path="data", stat_track=1)
        @test size(results["inputs"]["pd"],1) == K
		@test size(results["inputs"]["pd"],2) == 3
		@test size(results["outputs"]["p_gen"],1) == K
		@test size(results["outputs"]["p_gen"],2) == 5
		@test size(results["duals"]["v_max"],1) == K
		@test size(results["duals"]["v_max"],2) == 5
		@test haskey(results, "stats")
	end

    @testset "save infeasible" begin
		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path="data", save_infeasible=true)
        @test size(results["inputs"]["pd"],1) == K
		@test size(results["inputs"]["pd"],2) == 3
		@test size(results["outputs"]["p_gen"],1) == K
		@test size(results["outputs"]["p_gen"],2) == 5
		@test size(results["duals"]["v_max"],1) == K
		@test size(results["duals"]["v_max"],2) == 5
		@test haskey(results, "infeasible_inputs")
	end

    @testset "save found max loads" begin
		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path="data", save_max_load=true)
        @test size(results["inputs"]["pd"],1) == K
		@test size(results["inputs"]["pd"],2) == 3
		@test size(results["outputs"]["p_gen"],1) == K
		@test size(results["outputs"]["p_gen"],2) == 5
		@test size(results["duals"]["v_max"],1) == K
		@test size(results["duals"]["v_max"],2) == 5
		@test haskey(results, "load_constraints")
	end

	@testset "all save options" begin
		net_file = "pglib_opf_case5_pjm.m"
		K = 10
		results = create_samples(net_file, K, net_path="data",
							     save_certs=true, stat_track=1, save_infeasible=true)
		@test haskey(results, "inputs")
		@test haskey(results, "outputs")
		@test haskey(results, "duals")
		@test haskey(results, "infeasible_inputs")
		@test haskey(results, "polytope")
		@test haskey(results, "stats")
	end
end


text:
# An easy-to-construct AbstractEvent.
struct TestEvent <: AbstractEvent
    x
end
Discord.mock(::Type{TestEvent}) = TestEvent(rand())

# A module with event handlers.
module Handlers

export a, b

using Discord

a(::Client, ::AbstractEvent) = nothing
a(::Client, ::TypingStart) = nothing
b(::Client, ::WebhooksUpdate) = nothing
c(::Client, ::WebhooksUpdate) = nothing

end

@testset "Handlers" begin
    c = Client("token")
    f(c, e) = nothing
    g(c, e) = nothing
    badh(c, e::String) = nothing

    @testset "Basics" begin
        empty!(c.handlers)

        # Adding handlers without a tag means we can have duplicates.
        add_handler!(c, MessageCreate, f)
        add_handler!(c, MessageCreate, f)
        @test length(get(c.handlers, MessageCreate, Dict())) == 2

        # Deleting handlers without a tag clears all handlers for that type.
        delete_handler!(c, MessageCreate)
        @test !haskey(c.handlers, MessageCreate)

        # Using tags prevents duplicates.
        add_handler!(c, MessageCreate, f; tag=:f)
        add_handler!(c, MessageCreate, f; tag=:f)
        @test length(get(c.handlers, MessageCreate, Dict())) == 1
        @test haskey(c.handlers[MessageCreate], :f)

        # With tags, we can delete specific handlers.
        add_handler!(c, MessageCreate, g; tag=:g)
        @test length(get(c.handlers, MessageCreate, Dict())) == 2
        delete_handler!(c, MessageCreate, :g)
        @test length(get(c.handlers, MessageCreate, Dict())) == 1
        @test handler(c.handlers[MessageCreate][:f]) == f

        # We can also use do syntax.
        add_handler!(c, TypingStart) do c, e
            f(c, e)
        end
        @test length(get(c.handlers, TypingStart, Dict())) == 1
    end

    @testset "Adding from module" begin
        empty!(c.handlers)

        # We can add handlers from a module.
        add_handler!(c, Handlers)
        @test length(get(c.handlers, AbstractEvent, Dict())) == 1
        @test length(get(c.handlers, TypingStart, Dict())) == 1
        # Only exported functions are considered.
        @test length(get(c.handlers, WebhooksUpdate, Dict())) == 1
        @test handler(first(values(c.handlers[WebhooksUpdate]))) == Handlers.b

        # Adding a module handler with keywords propogates to all handlers.
        empty!(c.handlers)
        add_handler!(c, Handlers; tag=:h)
        @test all(hs -> all(p -> p.first === :h, hs), values(c.handlers))
    end

    @testset "Invalid handlers" begin
        # We can't add a handler without a valid method.
        @test_throws ArgumentError add_handler!(c, MessageCreate, badh)

        # We can't add a handler that's already expired.
        @test_throws ArgumentError add_handler!(c, Ready, f; count=0)
        @test_throws ArgumentError add_handler!(c, Ready, f; timeout=Day(-1))
    end

    @testset "Predicate functions" begin
        empty!(c.handlers)

        # By default, the predicate always returns true.
        add_handler!(c, TestEvent, f; tag=:f)
        @test predicate(c.handlers[TestEvent][:f]) == alwaystrue
        @test predicate(c.handlers[TestEvent][:f])(c, TestEvent(1))

        # But we can also specify our own.
        add_handler!(c, TestEvent, g; tag=:g, predicate=(c, e) -> e.x == 1)
        @test predicate(c.handlers[TestEvent][:g])(c, TestEvent(1))
        @test !predicate(c.handlers[TestEvent][:g])(c, TestEvent(2))
    end

    @testset "Fallback functions" begin
        # TODO
    end

    @testset "Expiries" begin
        empty!(c.handlers)

        # We can pass a number for a counting expiry.
        add_handler!(c, TestEvent, f; tag=:f, count=1)
        @test !isexpired(c.handlers[TestEvent][:f])
        dec!(c.handlers[TestEvent][:f])
        @test isexpired(c.handlers[TestEvent][:f])

        # Or we can use a timed expiry.
        add_handler!(c, TestEvent, g; tag=:g, timeout=Millisecond(100))
        @test !isexpired(c.handlers[TestEvent][:g])
        sleep(Millisecond(101))
        @test isexpired(c.handlers[TestEvent][:g])
    end

    @testset "Waiting" begin
        empty!(c.handlers)

        # Blocking handler with count should return when the count reaches 0.
        t = @async add_handler!(c, TestEvent, f; tag=:f, count=5, wait=true)
        sleep(Millisecond(1))
        h = c.handlers[TestEvent][:f]
        @test iscollecting(h)
        push!(h.results, 1)
        put!(h, results(h))
        h.remaining = 0
        sleep(Millisecond(150))
        @test istaskdone(t)
        @test fetch(t) == Any[1]

        # Blocking handler with timeout should return when the timeout elapses.
        t = @async add_handler!(c, TestEvent, f; tag=:g, timeout=Second(1), wait=true)
        sleep(Millisecond(1))
        h = c.handlers[TestEvent][:g]
        @test iscollecting(h)
        push!(h.results, 2)
        put!(h, results(h))
        sleep(1.1)
        @test istaskdone(t)
        @test fetch(t) == Any[2]

        # Blocking handler with both expiries should return when either is done.
        t = @async add_handler!(
            c, TestEvent, f;
            tag=:i, count=1, timeout=Day(1), wait=true,
        )
        sleep(Millisecond(1))
        h = c.handlers[TestEvent][:i]
        @test iscollecting(h)
        push!(h.results, 3)
        push!(h.results, 4)
        put!(h, results(h))
        h.remaining = -1
        sleep(Millisecond(150))
        @test istaskdone(t)
        @test fetch(t) == Any[3, 4]
    end

    @testset "Handler collection" begin
        delete_handler!(c, Ready)
        add_handler!(c, Ready, f)
        add_handler!(c, Ready, f)
        add_handler!(c, AbstractEvent, f)
        add_handler!(c, FallbackEvent, f)
        @test length(handlers(c, Ready)) == 2

        # No handlers means no handlers.
        empty!(c.handlers)
        @test isempty(allhandlers(c, MessageCreate))
        @test isempty(allhandlers(c, AbstractEvent))
        @test isempty(allhandlers(c, FallbackEvent))

        # Both the specific and catch-all handler should match.
        add_handler!(c, Ready, f)
        add_handler!(c, AbstractEvent, f)
        @test allhandlers(c, Ready) == [
            collect(c.handlers[AbstractEvent]);
            collect(c.handlers[Ready]);
        ]

        # The fallback handler should only match if there are non non-default handlers.
        add_handler!(c, FallbackEvent, f)
        @test allhandlers(c, Ready) == [
            collect(c.handlers[AbstractEvent]);
            collect(c.handlers[Ready]);
        ]
        delete_handler!(c, Ready)
        @test allhandlers(c, Ready) == collect(c.handlers[AbstractEvent])
        add_handler!(c, Ready, f)
        delete_handler!(c, AbstractEvent)
        @test allhandlers(c, Ready) == collect(c.handlers[Ready])
        delete_handler!(c, Ready)
        add_handler!(c, Ready, f; tag=DEFAULT_HANDLER_TAG)
        @test allhandlers(c, Ready) == [
            collect(c.handlers[Ready]);
            collect(c.handlers[FallbackEvent]);
        ]
        delete_handler!(c, Ready)
        @test allhandlers(c, Ready) == collect(c.handlers[FallbackEvent])
    end

    @testset "Default handler lookup" begin
        c = Client("token")
        @test hasdefault(c, MessageCreate)
        delete_handler!(c, MessageCreate, DEFAULT_HANDLER_TAG)
        @test !hasdefault(c, MessageCreate)
        add_handler!(c, MessageCreate, f; tag=DEFAULT_HANDLER_TAG)
        @test hasdefault(c, MessageCreate)
    end

    # TODO: Finish the tests.
end


text:
module Hanja

const 한자사전 = Dict

struct 한자타입
  코드::UInt16
  정체::String
  간체::String
  독음::String
  획수::Int
end

export 한자찾기
export 독음으로찾기
export 한글보기
export 정체보기
export 간체보기
export 한자옆에한글보기
export 한글옆에한자보기
export 한자옆에간체보기
export 간체옆에정체보기
export 간체옆에정체보기
include("finder.jl")

export 사전가져오기
include("dict.jl")

end # module


text:
using Folly

function matmul!(A, B, C)
    for i in 1:size(A, 1)
        for j in 1:size(A, 2)
            A[i, j] = zero(eltype(A))
            for k in 1:size(B, 2)
                A[i, j] += B[i, k] * C[k, j]
            end
        end
    end
end

loopnest = Loop(
    :i, [
        Loop(
            :j, [
                :(A[i, j] = zero(eltype(A))),
                Loop(:k, [
                    :(A[i, j] += B[i, k] * C[k, j])
                ])])
        ]
    )
indvars = [
    IndVar(:i, Outer(), :(i += 1), 1, :(size(A, 1))),
    IndVar(:j, Outer(), :(j += 1), 1, :(size(A, 2))),
    IndVar(:k, Inner(), :(k += 1), 1, :(size(B, 2))),
]

matmul2! = Folly.compile(Native(), loopnest, indvars, [:A, :B, :C])
matmul_ka! = Folly.compile(KATarget(CPU()), loopnest, indvars, [:A, :B, :C])

using Test

B = rand(64, 32)
C = rand(32, 64)
A = zeros(64, 64)

matmul!(A, B, C)
@test A ≈ B*C

A = zeros(64, 64)
matmul2!(A, B, C)
@test A ≈ B*C

A = zeros(64, 64)
event = matmul_ka!(A, B, C)
wait(event)
@test A ≈ B*C


text:
# Note that we tried and failed to automate the detection of "train"
# mode looking at the type of argument.  The argument is of type Value
# only during training and only if it is a value influenced by model
# weights.  However people typically apply dropout to the input which
# is not a Value.  So we are going back to no automation, make sure to
# supply p=0 during testing to stop dropout.

"""
    dropout(x, p)

Given an array `x` and probability `0<=p<=1`, just return `x` if
`p==0`, or return an array `y` in which each element is 0 with
probability `p` or `x[i]/(1-p)` with probability `1-p`.  Use
`seed::Number` to set the random number seed for reproducible
results. See [(Srivastava et al. 2014)](http://www.jmlr.org/papers/v15/srivastava14a.html)
for a reference.

"""
function dropout(x,p; seed=0)
    if 0 < p < 1
        if seed != 0; Knet.seed!(seed); end
        dropout!(p,x,similar(x))
    elseif p == 0
        x
    elseif p == 1
        zero(x)
    else
        error("Dropout probability not in [0:1]: $p")
    end
end

function dropback(dy,y,x,p)
    if 0 < p < 1
        dropback!(p,x,y,dy,similar(x))
    elseif p == 0
        dy
    elseif p == 1
        zero(x)
    else
        error("Dropout probability not in [0:1]: $p")
    end
end

# Turn dropout into an AutoGrad primitive
@primitive dropout(x,p;seed=0),dy,y dropback(value.((dy,y,x,p))...)

# GPU implementation
for S in (32,64)
    T = Symbol("Float$S")
    forw = Symbol("dropout_$S")
    back = Symbol("dropback_$S")
    @eval begin
        function dropout!(p::Number, x::KnetArray{$T}, y::KnetArray{$T})
            rand!(y)
            @knet8($forw,(Cint,$T,Ptr{$T},Ptr{$T}),length(y),$T(p),x,y)
            return y
        end
        function dropback!(p::Number, x::KnetArray{$T}, y::KnetArray{$T}, dy::KnetArray{$T}, dx::KnetArray{$T})
            @knet8($back,(Cint,$T,Ptr{$T},Ptr{$T},Ptr{$T},Ptr{$T}),length(dx),$T(p),x,y,dy,dx)
            return dx
        end
    end
end

# CPU implementation
function dropout!(p,x,y)
    rand!(y)
    p = convert(eltype(y),p)
    q = 1-p
    @inbounds for i=1:length(y)
        if y[i] > p
            y[i] = x[i] / q
        else
            y[i] = 0
        end
    end
    return y
end

function dropback!(p,x,y,dy,dx)
    p = convert(eltype(y),p)
    q = 1-p
    @inbounds for i=1:length(dx)
        if y[i] == 0
            dx[i] = 0
        else
            dx[i] = dy[i] / q
        end
    end
    return dx
end


text:
include(joinpath(@__DIR__, "../src/utils/data.jl"))
import .UtilsData: TestData

download(TestData; branch = "master")


text:
function process(srcpath::String, verbose::Bool = false)
	if verbose
		println("Extracting document chunks")
	end
	chunks = extractdocs(srcpath)

	n = length(chunks)

	summaries = Array(UTF8String, n)
	for i in 1:n
		if verbose
			println("Lexing document chunk $i")
		end
		lexemes = lex(chunks[i])

		if verbose
			println("Parsing document chunk $i")
		end
		doc = parse_plus_semantics(lexemes)

		if verbose
			println("Generating markdown for document $i")
		end
		summaries[i] = markdown(doc)
	end
	return summaries
end

function roxygenize(srcdir::String, destdir::String; verbose::Bool = false)
	srcpaths = readdir(srcdir)
	if !ispath(destdir)
		try
			mkdir(destdir)
		catch
			error("Failed to make directory $destdir")
		end
	end
	for srcpath in srcpaths
		if endswith(srcpath, ".jl")
			if verbose
				@printf("Processing file '%s'\n", srcpath)
			end
			destpath = joinpath(destdir, replace(srcpath, ".jl", ".md"))
			summary = process(joinpath(srcdir, srcpath), verbose)
			if length(summary) > 0
				io = open(destpath, "w")
				print(io, join(summary, "\n---\n\n"))
				close(io)
			end
		end
	end
end


text:
using Parameters
using Random
using Sobol
using SpecialFunctions

"""
Returns particles drawn from g for the bump-on-tail case.

input: nb. pf particles Nₚ, x-marginal function gₓ, parameters μ
output: Particles struct
"""
function draw_g_accept_reject(N::Int, fₓ::Function, params::NamedTuple) where {T}
    @unpack κ, ε, a, v₀, σ = params

    x = zeros(N); v = zeros(N); w = zeros(N)

    # Sobol sampling
    s = SobolSeq(2); skip(s, 2N)

    n = 1
    while n < N+1
        ###
        y = next!(s)
        x₀ = y[1]*2π/κ  # proposal, uniform

        # accept-reject in x
        if rand(1)[1] > ( (fₓ(x₀, params))/(1 + ε) )
            continue                            # reject
        end
        x[n] = x₀                               # accept
        w[n] = 2π / κ / N

        # inverse CDF sampling in v
        vₛ = 2y[2] - 1
        v[n] = √2 * erfinv(vₛ)  # draw from bulk
        if rand(1)[1] > 1-a
            v[n] = v[n] * σ + v₀  # draw from tail
        end
        n = n+1
    end

    return ParticleList(x, v, w)
end

"""
Returns particles drawn from g for the bump-on-tail case.

input: nb. pf particles Nₚ, x-marginal function gₓ, parameters μ
output: Particles struct
"""
function draw_g_importance_sampling(N::Int, fₓ::Function, params::NamedTuple) where {T}
    @unpack κ, ε, a, v₀, σ = params

    x = zeros(N); v = zeros(N); w = zeros(N)

    # Sobol sampling
    s = SobolSeq(2); skip(s, 2N)

    n = 1
    while n < N+1
        ###
        y = next!(s)
        x₀ = y[1] * 2π / κ  # proposal, uniform

        # importance sampling in x
        x[n] = x₀
        w[n] = fₓ(x₀, params) .* 2π ./ κ ./ N    # *L is equivalent to 1/(1/L) i.e. uniform proposal. f is NOT normalized!

        # inverse CDF sampling in v
        vₛ = 2y[2] - 1
        v[n] = √2 * erfinv(vₛ)  # draw from bulk
        if rand(1)[1] > 1-a
            v[n] = v[n]*σ + v₀  # draw from tail
        end
        n = n+1
    end

    return ParticleList(x, v, w)
end

"""
Returns re-weighted particles according to new parameters.

input: Particles struct P₀, proposal distribution function g, target distribution function f
output: Particles struct
"""
function weight_f(P₀, g, f)
    Px = zero(P₀.x); Pv = zero(P₀.v); Pw = zero(P₀.w)
    Px .= P₀.x; Pv .= P₀.v; Pw .= P₀.w
    P = ParticleList(Px, Pv, Pw)
    P.w .*= f(P₀.x, P₀.v) ./ g(P₀.x, P₀.v)
    return P
end


text:
# MCMC.jl
"""
My MCMC stuff

Currently the only method is gibbs. Eventually, I would
like to add a Metropolis Hastings method.
"""
module MCMC

import Distributions
export gibbs, metropolis
"""
Generalized Gibbs sampler.

# Arguments
* `init`: Initial values for gibbs sampler
* `update`: function which takes a state::State as the only input and returns the update
* `B`: number of MCMC samples desired
* `burn`: burn-in
* `printFreq`: shows progress bar every printFreq%. (default: 0 => don't print progress)

# Example:

```julia
using Distributions

# Generate data
const n = 1000
const (mu,sig2) = (5,2)
const y = rand(Normal(mu,sqrt(sig2)),n)

# precomputes
const ybar = mean(y)
const (sig2_a, sig2_b) = (2.0, 1.0)
rig(shape::Float64,rate::Float64) = 1 / rand(Gamma(shape,1/rate))

# Define State
immutable State
  mu::Float64
  sig2::Float64
end

# define an update function for State
function update(state::State)
  const newMu = rand(Normal(ybar,sqrt(state.sig2/n)))
  const newSig2 = rig(sig2_a+n/2, sig2_b + sum((y-newMu).^2)/2)
  State(newMu,newSig2)
end

# run gibbs
@time out = gibbs(State(0,1),update,10000,1000);

## post processing
const postMu = map(o -> o.mu, out)
const postSig2 = map(o -> o.sig2, out)

```
"""

function gibbs{T}(init::T, update, B::Int, burn::Int; printFreq::Int=0)
  const out = Array{T,1}( (B+burn) )
  out[1] = init
  for i in 2:(B+burn)
    out[i] = update(out[i-1])
    if printFreq > 0 && i % printFreq == 0
      print("\rProgress: ",i,"/",B+burn)
    end
  end
  out[ (burn+1):end ]
end

"""
Univariate Metropolis step (with Normal proposal)

    metropolis(curr, loglike_plus_logprior, candSig::Real;
               inbounds = x-> -Inf<x<Inf)

# Arguments
* `curr`: current value of parameter
* `loglike_plus_logprior`: log-likelihoodd plus log-prior as function of the parameter to be updated
* `candSig`: The sd of the normal proposal
* `inbounds`: A function which checks if the current parameter is in the support
"""
function metropolis(curr, loglike_plus_logprior, candSig::Real;
                    inbounds = x-> -Inf<x<Inf)
  # need to do autotune

  const cand = rand(Distributions.Normal(curr,candSig))

  if inbounds(cand)
    p = loglike_plus_logprior(cand) - loglike_plus_logprior(curr)
    p > log(rand()) ?  cand : curr
  else
    curr
  end

end

end # module


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.test.behaviour.util;
#
# import static org.junit.Assert.assertTrue;
# import static org.junit.Assert.fail;
#
# public class Util {
#
#     public static void assertThrows(Runnable function) {
#         try {
#             function.run();
#             fail();
#         } catch (RuntimeException e) {
#             assertTrue(true);
#         }
#     }
#
#     public static void assertThrowsWithMessage(Runnable function, String message) {
#         try {
#             function.run();
#             fail();
#         } catch (RuntimeException e) {
#             System.out.println(e.toString());
#             assert (e.toString().contains(message));
#         }
#     }
# }


text:
using MPI, CUDAapi, Requires

const cuarray_pkgid = Base.PkgId(Base.UUID("3a865a2d-5b23-5a0f-bc46-62713ec82fae"), "CuArrays")

"""
    CLIMA.array_type()

Returns the default array type to be used with CLIMA, either a `CuArray` if a GPU is
available, or an `Array` otherwise.

Use of GPUs can be explicitly disabled by setting the environment variable `CLIMA_GPU=false`.
"""
function array_type()
  if get(ENV, "CLIMA_GPU", "") != "false" && CUDAapi.has_cuda_gpu()
    CuArrays = Base.require(cuarray_pkgid)
    CuArrays.CuArray
  else
    Array
  end
end

"""
    CLIMA.init(ArrayType)

Initialize MPI and allocate GPUs among MPI ranks if using GPUs.
"""
function init(ArrayType=array_type())
  if !MPI.Initialized()
    MPI.Init()
  end
  _init_array(ArrayType)
  return nothing
end

function _init_array(::Type{Array})
  nothing
end

@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
  using .CuArrays, .CuArrays.CUDAdrv, .CuArrays.CUDAnative
  function _init_array(::Type{CuArray})
    comm = MPI.COMM_WORLD
    # allocate GPUs among MPI ranks
    local_comm = MPI.Comm_split_type(comm, MPI.MPI_COMM_TYPE_SHARED,  MPI.Comm_rank(comm))
    # we intentionally oversubscribe GPUs for testing: may want to disable this for production
    device!(MPI.Comm_rank(local_comm) % length(devices()))
    CuArrays.allowscalar(false)
  end
end

function gpu_allowscalar(val)
  if haskey(Base.loaded_modules, CLIMA.cuarray_pkgid)
    Base.loaded_modules[CLIMA.cuarray_pkgid].allowscalar(val)
  end
  return nothing
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
import Statistics

loss_simple(W, b, x) = sum(W * x .+ b)
loss_double_broadcast(W, b, x) = sum(sin.(W * x) .+ b)
loss_kw_mean(W, b, x) = Statistics.mean(W * x .+ b; dims=1)[1]

@testset "grad: basic" begin
    args = (rand(3, 4), rand(3), rand(4))
    @test gradcheck(loss_simple, args...)
    @test gradcheck(loss_double_broadcast, args...)

    val, g = grad(loss_kw_mean, args...)
    @test val == loss_kw_mean(args...)
    @test gradcheck(loss_kw_mean, args...)
end

@testset "grad: getindex" begin
    x = rand(3, 4, 5)
    x1 = zero(x); x1[1] = 1
    @test grad(x -> x[1], x)[2][1] == x1
    x2 = zero(x); x2[1, 2, 1] = 1
    @test grad(x -> x[1, 2, 1], x)[2][1] == x2
    x3 = zero(x); x3[:, 1, :] .= 1
    @test grad(x -> sum(x[:, 1, :]), x)[2][1] == x3
end

sum_bcast(x, y) = sum(x .+ y)

@testset "special bcast" begin
    for args in [
        (rand(3, 4), rand(3)),
        (rand(3, 4), rand(3, 1)),
        (rand(3, 4), rand(1, 4)),
        (rand(3), rand(3, 4)),
        (rand(3, 1), rand(3, 4)),
        (rand(1, 4), rand(3, 4)),
        # (rand(3, 4), rand()),
        # (rand(), rand(3, 4)),
    ]
        @test gradcheck(sum_bcast, args...)
        val, g = grad(sum_bcast, args...)
        for i=1:length(args)
            @test size(g[i]) == size(args[i])
        end
    end
end

sum_prod_bcast(x, y) = sum(x .* y)

@testset "special bcast 2" begin
    for args in [
        (rand(3, 4), rand(3)),
        (rand(3, 4), rand(3, 1)),
        (rand(3, 4), rand(1, 4)),
        (rand(3), rand(3, 4)),
        (rand(3, 1), rand(3, 4)),
        (rand(1, 4), rand(3, 4)),
        # (rand(3, 4), rand()),
        # (rand(), rand(3, 4)),
    ]
        @test gradcheck(sum_prod_bcast, args...)
        val, g = grad(sum_prod_bcast, args...)
        for i=1:length(args)
            @test size(g[i]) == size(args[i])
        end
    end
end

@testset "grad: literal_pow" begin
    @test gradcheck(x -> sum(x .^ 2), rand(3, 4))
end

@testset "grad: transpose" begin
    a = rand(5)
    b = rand(2, 3)
    @test grad(x -> sum(transpose(x)), a)[2][1] == ones(size(a))
    @test grad(x -> sum(adjoint(x)), a)[2][1] == ones(size(a))
    @test grad(x -> sum(transpose(x)), b)[2][1] == ones(size(b))
    @test grad(x -> sum(adjoint(x)), b)[2][1] == ones(size(b))
end

mutable struct Linear{T}
    W::AbstractArray{T,2}
    b::AbstractArray{T}
end

forward(m::Linear, X) = m.W * X .+ m.b

loss(m::Linear, X) = sum(forward(m, X))

mutable struct Point x; y end
constructor_loss(a) = (p = Point(a, a); p.x + p.y)

@testset "grad: structs" begin
    args = Any[Linear(rand(3,4), rand(3)), rand(4,5)]
    val, g = grad(loss, args...)

    @test val == loss(args...)

    play!(g.tape)
    val, g = grad(loss, args...)
    @test val == loss(args...)

    _, g = grad(constructor_loss, 4.0)
    @test g[1] == 2
end

mutable struct Line p1; p2 end

function add_points(x)
    p = Point(x, x)
    l = Line(p, p)
    return 2*l.p1.x + 5*l.p2.y
end

@testset "grad: structs/new" begin
    # find_field_source_var
    _, tape = trace(add_points, rand())
    src_op = find_field_source_var(tape, tape[15])
    @test src_op.id == 1
    @test src_op.val isa Real
    src_op = find_field_source_var(tape, tape[13])
    @test src_op.id == 3
    @test src_op.val isa Point

    @test grad(add_points, rand())[2][1] == 7
end

make_t(x) = tuple(x, x)

function make_t_loss(x)
    a, b = make_t(x)
    return a + b
end

@testset "grad: tuple unpack" begin
    _, g = grad(make_t_loss, 3.0)
    @test g[1] == 2.0
end

# HESS = randn(3,3)
# hessian_fun(x) = x'*(HESS*x)
# hessian_fun2(x) = 0.5*x'*(HESS*x)

# @testset "grad: adjoint" begin
#     H = randn(3,3)
#     x = rand(3)
#     val, g = grad(hessian_fun2, x)
#     @test val isa Real
# end

@testset "grad: compiled" begin
    args = Any[Linear(rand(3,4), rand(3)), rand(4,5)]

    val, g = grad(loss, args...)

    # compiled vs. non-compiled
    tape = g.tape
    args = Any[Linear(rand(3,4), rand(3)), rand(4,5)]

    val1 = play!(tape, args...; use_compiled=true)
    last_val1 = tape[end].val
    val2 = play!(tape, args...; use_compiled=false)
    last_val2 = tape[end].val
    @test val1 == val2
    @test last_val1 == last_val2

    # # (* -> mul!) for mixed array-scalar vars
    # x = rand(3)
    # val1, g1 = grad(hessian_fun, x)  # interpreted
    # val2, g2 = grad(hessian_fun, x)  # compiled
    # @test val1 == val2
    # @test g1[1] == g2[1]
end


text:
# # Chapter 4
using BenchmarkTools

##Globals

p = 2

function pow_array(x::Vector{Float64})
  s = 0.0
    for y in x
      s = s + y^p
    end
    return s
end
t=rand(100000);

@btime pow_array($t)

@code_warntype pow_array(t)

### const

const p2 = 2
function pow_array2(x::Vector{Float64})
    s = 0.0
    for y in x
        s = s + y^p2
    end
    return s
end

 @btime pow_array2($t)

 @code_warntype pow_array2(t)

 function pow_array3(x::Vector{Float64})
    return pow_array_inner(x, p)
end

function pow_array_inner(x, pow)
   s = 0.0
   for y in x
      s = s + y^pow
   end
   return s
end

@btime pow_array3($t)

## Inlining

function f(x)
    a=x*5
    b=a+3
end

g(x) = f(2*x)

@code_typed g(3)

@code_llvm g(3)

function f(x)
    a=x*5
    b=a+3
    c=a-4
    if c < 0
        throw(DomainError())
    elseif c < 2
        d=c^3
    else
        d=c^2
    end
end

@inline function f_in(x)
    a=x*5
    b=a+3
    c=a-4
    if c < 0
        throw(DomainError())
    elseif c < 2
        d=c^3
    else
        d=c^2
    end
end

g_in(x) = f_in(2*x)

@noinline function f_ni(x)
    a=x*5
    b=a+3
end

g_ni(x) = f_ni(2*x)

@code_typed g_ni(3)

#Constant Propagation

sqr(x) = x * x

sqr2() = sqr(2)

@code_typed sqr2()

#Macros

function poly_naive(x, a...)
  p=zero(x)
  for i = 1:length(a)
    p = p + a[i] *  x^(i-1)
  end
  return p
end

f_naive(x) = poly_naive(x, 1,2,3,4,5)

x=3.5

@btime f_naive(3.5)

function poly_horner(x, a...)
  b=zero(x)
    for i = length(a):-1:1
      b = a[i] + b * x
    end
    return b
end

f_horner(x) = poly_horner(x, 1,2,3,4,5)

@btime f_horner($x)

macro horner(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :(muladd(t, $ex, $(esc(p[i]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

f_horner_macro(x) = @horner(x, 1,2,3,4,5,6,7,8,9,10)

@btime f_horner_macro($x)

# Generated Functions

function prod_dim(x::Array{T, N}) where {T, N}
    s = 1
    for i = 1:N
        s = s * size(x, i)
    end
    return s
end

prod_dim(rand(10, 5, 5))

@generated function prod_dim_gen(x::Array{T, N}) where {T, N}
     ex = :(1)
     for i = 1:N
          ex = :(size(x, $i) * $ex)
     end
     return ex
end

prod_dim_gen(rand(10, 5, 5))

function prod_dim_gen_impl(x::Array{T, N}) where {T,N}
    ex = :(1)
    for i = 1:N
       ex = :(size(x, $i) * $ex)
    end
    return ex
end

#Keyword arguments
named_param(x; y=1, z=1) = x^y + x^z

pos_param(x,y,z) = x^y + x^z

@btime named_param(4; y=2, z=3)

@btime pos_param(4, 2, 3)