text:
module SIWPD
export
    siwpd,
    makesiwpdtree

using
    Wavelets

using
    ..Utils,
    ..DWT

"""
    siwpd(x, wt[, L=maxtransformlevels(x), d=L])

Computes the Shift-Invariant Wavelet Packet Decomposition originally developed
by Cohen, Raz & Malah on the vector `x` using the discrete wavelet filter `wt`
for `L` levels with depth `d`.
"""
function siwpd(x::AbstractArray{T}, wt::OrthoFilter,
        L::Integer=maxtransformlevels(x), d::Integer=L) where T<:Number

    g, h = WT.makereverseqmfpair(wt, true)
    si = Vector{T}(undef, length(wt)-1)                 # tmp filter vector
    return siwpd(x, wt, h, g, si, L, d)
end

function siwpd(x::AbstractVector{Tx}, wt::OrthoFilter, h::Vector{T},
        g::Vector{T}, si::Vector{T}, L::Integer, d::Integer=L) where
        {Tx<:Number, T<:Number}

    n = length(x)
    @assert isdyadic(n)
    @assert 0 <= L <= maxtransformlevels(n)
    @assert 1 <= d <= L

    # decomposed array y
    y = Array{T,2}(undef, (n,2^(L+1)-1))
    y[:,1] = x

    for node in axes(y,2)
        level = floor(Int, log2(node)) # current node level (starts from 0)
        nodelen = nodelength(n, level)
        # decomposition depth of current node
        dₗ = (0 <= level <= L-d) ? d : L-level
        decompose!(y, wt, h, g, si, node, nodelen,
            level, dₗ, 0, zeroshiftdecomposition=false)
    end

    return y
end

function decompose!(y::AbstractArray{Ty,2}, wt::OrthoFilter,
        h::Vector{T}, g::Vector{T}, si::Vector{T},
        node::Integer, nodelen::Integer, level::Integer, d::Integer,
        shift::Integer; zeroshiftdecomposition::Bool=true) where
        {Ty<:Number, T<:Number}

    # base case
    @assert d > 0 || return nothing
    @assert 0 <= level <= maxtransformlevels(size(y,1))
    @assert d <= maxtransformlevels(size(y,2)+1) - level
    @assert 0 <= shift < size(y,1)
    filtlen = length(wt)

    pstart = shift*nodelen + 1                                  # start index
    pend = (shift+1)*nodelen                                    # end index
    parent = y[pstart:pend, node]

    cnodelen = nodelen ÷ 2
    cstart = shift*cnodelen + 1
    lchild = @view y[:, getchildindex(node,:left)]
    rchild = @view y[:, getchildindex(node,:right)]

    # non-shifted subtree
    Wavelets.Transforms.filtdown!(      # scaling coefficients
        g, si, lchild, cstart, cnodelen, parent, 1, 0, false
    )
    Wavelets.Transforms.filtdown!(      # detail coefficients
        h, si, rchild, cstart, cnodelen, parent, 1, -filtlen+1, true
    )
    if zeroshiftdecomposition   # pursue the same thing with d-1 if shift > 0
        decompose!(                     # decompose left child
            y, wt, h, g, si, getchildindex(node,:left), cnodelen, level+1, d-1, shift
        )
        decompose!(                     # decompose right child
            y, wt, h, g, si, getchildindex(node,:right), cnodelen, level+1, d-1, shift
        )
    end

    # shifted subtree
    parent = circshift(parent,1)
    cstart = (shift + 1<<level)*cnodelen + 1
    Wavelets.Transforms.filtdown!(      # scaling coefficients
        g, si, lchild, cstart, cnodelen, parent, 1, 0, false
    )
    Wavelets.Transforms.filtdown!(      # detail coefficients
        h, si, rchild, cstart, cnodelen, parent, 1, -filtlen+1, true
    )
    # pursue the same thing with d-1
    decompose!(                         # decompose shifted left child
        y, wt, h, g, si, getchildindex(node,:left), cnodelen, level+1, d-1, shift+1<<level
    )
    decompose!(                         # decomposed shifted right child
        y, wt, h, g, si, getchildindex(node,:right), cnodelen, level+1, d-1, shift+1<<level
    )

    return nothing
end

"""
    makesiwpdtree(n, L, d)

Returns the multi-level, multi-depth binary tree corresponding to the Shift-
Invariant Wavelet Packet Decomposition.
"""
function makesiwpdtree(n::Integer, L::Integer, d::Integer)
    @assert 0 <= L <= maxtransformlevels(n)
    @assert 1 <= d <= L

    tree = Vector{BitVector}(undef, 2^(L+1)-1)
    for i in eachindex(tree)
        level = floor(Int, log2(i))
        len = nodelength(n, level)
        nshift = n ÷ len        # number of possible shifts for subspace Ω(i,j)
        node = falses(nshift)
        k = ceil(Int, nshift / 1<<d)  # "skips" between shifts
        node[1:k:end] .= true         # mark nodes present in tree
        tree[i] = node
    end
    return tree
end

end # end module


text:
"""
Created on Mon Apr 13 2020

@author: Yoann Pradat

Plot divergences
"""

include("./src/VNMF.jl")
using .VNMF
using LaTeXStrings
using Plots

#### # Plot β-divergence for different values of β
#### #################################################################################################################

βys = [-1, 0, 1, 2]
βzs = [-0.5, 0.5, 1.5, 2.5]
titles = [L"(a) \beta < 0", L"(b) 0 \leq \beta < 1", L"(c) 1 \leq \beta < 2", L"(d) 2 \leq \beta"]

h = 0.05
p = plot(layout=(2,2), dpi=300, size=(1000, 600))

for sp=1:4
    βy = βys[sp]
    βz = βzs[sp]

    if βz < 1
	X = h:h:8
    else
	X = h:h:4
    end

    Y = [β_divergence(hcat([1.]), hcat([x]), βy) for x in X]
    Z = [β_divergence(hcat([1.]), hcat([x]), βz) for x in X]

    plot!(
	X,
	hcat(Y, Z),
	title = titles[sp],
	labels = ["\$\\beta = $(βy)\$" "\$\\beta = $(βz)\$" ],
	ls    = [:solid :dash],
	lw    = 1,
	color = "black",
	subplot = sp
    )

    plot!(
	xlabel = "",
	xlims  = (0,X[end]),
	xticks = 0:X[end]/4:X[end],
	ylims  = (0,1.5),
	yticks = 0:0.5:1.5,
	subplot = sp
    )
end

folder   = "demo/plot/divergence"
mkpath(folder)
filename = "beta_divergences.pdf"
filepath = joinpath(folder, filename)

savefig(p, filepath)
println("plot saved at $(filepath)")

ext: jl


text:
struct NuclearNormPenalty <: ArrayPenalty end
function value(p::NuclearNormPenalty, A::AbstractMatrix{<:Number})
    >(size(A)...) ? tr(sqrt(A'A)) : tr(sqrt(A * A'))
end
function prox!(p::NuclearNormPenalty, A::AbstractMatrix{<:Number}, λ::Number)
    svdecomp = svd!(A)
    soft_thresh!(svdecomp.S, λ)
    copyto!(A, Matrix(svdecomp))
end

struct GroupLassoPenalty <: ArrayPenalty end
value(p::GroupLassoPenalty, A::AbstractArray{<:Number}) = norm(A)
function prox!(p::GroupLassoPenalty, A::AbstractArray{T}, λ::Number) where {T <: Number}
    denom = norm(A)
    if denom <= λ
        fill!(A, zero(T))
    else
        scaling = λ / denom
        for i in eachindex(A)
            @inbounds A[i] = (1 - scaling) * A[i]
        end
    end
    A
end

"""
    MahalanobisPenalty(C)

Supports a Mahalanobis distance penalty (`xᵀCᵀCx` for a vector `x`).
"""
mutable struct MahalanobisPenalty{T <: Number, S <: AbstractArray{T,2}} <: ArrayPenalty
    C::S
    CtC::S
    CtC_Iλ::LU{T, Matrix{T}} # LU factorization of C'C + I/λ
    λ::T
end
function MahalanobisPenalty(C::AbstractMatrix{T}, λ::T = one(T)) where {T<:Number}
    MahalanobisPenalty(C, C'C, lu(C'C + I), λ)
end
function value(p::MahalanobisPenalty{T}, x::AbstractVector{T}) where {T <: Number}
    inv(T(2)) * T(sum(abs2, p.C * x))
end
function prox!(p::MahalanobisPenalty{T}, A::AbstractArray{T, 1}, λ::Number) where {T <: Number}
    if λ != p.λ
        p.λ = λ
        p.CtC_Iλ = lu(p.CtC + I / λ)
    end
    rmul!(A, one(T) / λ)
    ldiv!(p.CtC_Iλ, A) # overwrites result in A
end

struct ScaledArrayPenalty{T, P <: ArrayPenalty} <: ArrayPenalty
    penalty::P
    λ::T
end
scaled(p::ArrayPenalty, λ::Number) = (_scale_check(λ); ScaledArrayPenalty(p, λ))

value(p::ScaledArrayPenalty, θ::AbstractArray{<:Number}) = p.λ * value(p.penalty, θ)
prox!(p::ScaledArrayPenalty, θ::AbstractArray{<:Number}) = prox!(p.penalty, θ, p.λ)


text:
using Images, WNNMDenoise, TestImages
using Random
using Test

@testset "WNNMDenoise.jl" begin
    imgfile = joinpath(pkgdir(WNNMDenoise), "house.png")
    img = imresize(Float64.(load(imgfile)) * 255, ratio=0.5);
    clean_img = img;
    σₙ = 40;
    Random.seed!(0)
    noisy_img = clean_img .+ σₙ .* randn(Float32, size(clean_img))

    @test assess_psnr(noisy_img, img, 255) ≈ 16.08294 atol=1e-5

    out = copy(noisy_img)
    f = WNNM(σₙ)
    f(out, noisy_img)

    println("PSNR: ", assess_psnr(out, img, 255))
    @test assess_psnr(out, img, 255) >= 28.8
end


text:
#

using LibArchive
using Base.Test

using Compat

## Version
info("Testing libarchive $(LibArchive.version()::VersionNumber)")

## Error
info("Test error translation")
@test_throws EOFError LibArchive._la_error(LibArchive.Status.EOF)
@test_throws ArchiveRetry LibArchive._la_error(LibArchive.Status.RETRY)
@test_throws ArchiveFailed LibArchive._la_error(LibArchive.Status.FAILED)
@test_throws ArchiveFatal LibArchive._la_error(LibArchive.Status.FATAL)
@test_throws ErrorException LibArchive._la_error(Cint(-1024))

## Reader error
info("Test reader error handling")
let
    reader = LibArchive.Reader()
    @test reader.ptr != C_NULL
    LibArchive.free(reader)
    @test reader.ptr == C_NULL
    LibArchive.free(reader)
    @test_throws ErrorException LibArchive.support_filter_all(reader)
end

LibArchive.Reader() do reader
    LibArchive.set_exception(reader, EOFError())
    @test Libc.errno(reader) == LibArchive.Status.EOF
    @test LibArchive.error_string(reader) == "end of file"
    LibArchive.clear_error(reader)

    LibArchive.set_exception(reader, ArchiveRetry("retry"))
    @test Libc.errno(reader) == LibArchive.Status.RETRY
    @test LibArchive.error_string(reader) == "retry"
    LibArchive.clear_error(reader)

    LibArchive.set_exception(reader, ArchiveFailed("failed"))
    @test Libc.errno(reader) == LibArchive.Status.FAILED
    @test LibArchive.error_string(reader) == "failed"
    LibArchive.clear_error(reader)

    LibArchive.set_exception(reader, ArchiveFatal("fatal"))
    @test Libc.errno(reader) == LibArchive.Status.FATAL
    @test LibArchive.error_string(reader) == "fatal"
    LibArchive.Reader() do reader2
        LibArchive.copy_error(reader2, reader)
        @test Libc.errno(reader2) == LibArchive.Status.FATAL
        @test LibArchive.error_string(reader2) == "fatal"
    end
    LibArchive.clear_error(reader)

    err_ex = ErrorException("error")
    LibArchive.set_exception(reader, err_ex)
    @test Libc.errno(reader) == LibArchive.Status.FAILED
    @test LibArchive.error_string(reader) == string(err_ex)
    LibArchive.clear_error(reader)
end

LibArchive.Reader("/this_file_does_not_exist") do reader
    local ex
    try
        LibArchive.next_header(reader)
    catch ex
    end
    @test isa(ex, ArchiveFatal)
    @test !isempty(ex.msg)
end

# Writer error
info("Test writer error handling")
let
    writer = LibArchive.Writer()
    @test writer.ptr != C_NULL
    LibArchive.free(writer)
    @test writer.ptr == C_NULL
    LibArchive.free(writer)
    @test_throws ErrorException LibArchive.add_filter_bzip2(writer)
end

LibArchive.Writer("/this_dir_does_not_exist/file") do writer
    local ex
    try
        LibArchive.finish_entry(writer)
    catch ex
    end
    @test isa(ex, ArchiveFatal)
    @test !isempty(ex.msg)
end

info("Test do block result and exception pass through")

@test_throws BoundsError LibArchive.Writer(_->throw(BoundsError([], 10)))
@test_throws ArgumentError LibArchive.Reader(_->throw(ArgumentError("a")))
@test LibArchive.Writer(_->10) === 10
@test LibArchive.Reader(_->0.1) === 0.1

info("Test availability of filters and formats")
info("    Reader")
LibArchive.Reader() do reader
    LibArchive.support_filter_all(reader)
end

LibArchive.Reader() do reader
    LibArchive.support_filter_bzip2(reader)
    LibArchive.support_filter_compress(reader)
    if LibArchive.version() >= v"3.1.0"
        LibArchive.support_filter_grzip(reader)
        LibArchive.support_filter_lrzip(reader)
        LibArchive.support_filter_lzip(reader)
        LibArchive.support_filter_lzop(reader)
    end
    LibArchive.support_filter_lzma(reader)
    LibArchive.support_filter_rpm(reader)
    LibArchive.support_filter_uu(reader)
    LibArchive.support_filter_xz(reader)
end

LibArchive.Reader() do reader
    LibArchive.support_format_all(reader)
end

LibArchive.Reader() do reader
    LibArchive.support_format_7zip(reader)
    LibArchive.support_format_ar(reader)
    LibArchive.support_format_by_code(reader, LibArchive.Format._7ZIP)
    LibArchive.support_format_cab(reader)
    LibArchive.support_format_cpio(reader)
    LibArchive.support_format_empty(reader)
    LibArchive.support_format_gnutar(reader)
    LibArchive.support_format_iso9660(reader)
    LibArchive.support_format_lha(reader)
    LibArchive.support_format_mtree(reader)
    LibArchive.support_format_rar(reader)
    LibArchive.support_format_raw(reader)
    LibArchive.support_format_tar(reader)
    LibArchive.support_format_xar(reader)
    LibArchive.support_format_zip(reader)
end

if LibArchive.version() >= v"3.1.0"
    LibArchive.Reader() do reader
        LibArchive.set_format(reader, LibArchive.Format.TAR)
        LibArchive.append_filter(reader, LibArchive.FilterType.BZIP2)
    end
end

info("    Writer")
if LibArchive.version() >= v"3.1.0"
    LibArchive.Writer() do writer
        LibArchive.add_filter(writer, LibArchive.FilterType.COMPRESS)
        LibArchive.add_filter(writer, "bzip2")
    end
end

LibArchive.Writer() do writer
    if LibArchive.version() >= v"3.1.0"
        LibArchive.add_filter_b64encode(writer)
        LibArchive.add_filter_grzip(writer)
        LibArchive.add_filter_lrzip(writer)
        LibArchive.add_filter_lzop(writer)
        LibArchive.add_filter_uuencode(writer)
    end
    LibArchive.add_filter_bzip2(writer)
    LibArchive.add_filter_compress(writer)
    LibArchive.add_filter_gzip(writer)
    LibArchive.add_filter_none(writer)
    if !is_apple()
        LibArchive.add_filter_lzip(writer)
        LibArchive.add_filter_lzma(writer)
        LibArchive.add_filter_xz(writer)
    end
end

LibArchive.Writer() do writer
    LibArchive.set_format(writer, LibArchive.Format.CPIO)
    if LibArchive.version() >= v"3.1.0"
        LibArchive.set_format(writer, "gnutar")
    end
end

LibArchive.Writer() do writer
    if LibArchive.version() >= v"3.1.0"
        LibArchive.set_format_mtree_classic(writer)
        LibArchive.set_format_v7tar(writer)
    end
    LibArchive.set_format_7zip(writer)
    LibArchive.set_format_ar_bsd(writer)
    LibArchive.set_format_ar_svr4(writer)
    LibArchive.set_format_cpio(writer)
    LibArchive.set_format_cpio_newc(writer)
    LibArchive.set_format_gnutar(writer)
    LibArchive.set_format_iso9660(writer)
    LibArchive.set_format_mtree(writer)
    LibArchive.set_format_pax(writer)
    LibArchive.set_format_pax_restricted(writer)
    LibArchive.set_format_shar(writer)
    LibArchive.set_format_shar_dump(writer)
    LibArchive.set_format_ustar(writer)
    LibArchive.set_format_xar(writer)
    LibArchive.set_format_zip(writer)
end

# Copy entry
info("Test deepcopy of Entry")
let
    entry = LibArchive.Entry()
    @test !LibArchive.size_is_set(entry)
    LibArchive.set_pathname(entry, "a.txt")
    LibArchive.set_size(entry, 10)
    LibArchive.set_perm(entry, 0o644)
    LibArchive.set_filetype(entry, LibArchive.FileType.REG)

    @test LibArchive.pathname(entry) == "a.txt"
    @test LibArchive.size(entry) == 10
    @test LibArchive.size_is_set(entry)
    @test LibArchive.perm(entry) == 0o644
    @test LibArchive.filetype(entry) == LibArchive.FileType.REG

    entry_cp = deepcopy(entry)
    @test LibArchive.pathname(entry_cp) == "a.txt"
    @test LibArchive.size(entry_cp) == 10
    @test LibArchive.size_is_set(entry_cp)
    @test LibArchive.perm(entry_cp) == 0o644
    @test LibArchive.filetype(entry_cp) == LibArchive.FileType.REG

    LibArchive.clear(entry)
    LibArchive.free(entry)
    LibArchive.free(entry_cp)
end

# Entry properties
info("Test Entry properties")
let
    # Time stamps
    entry = LibArchive.Entry()
    t = floor(Int, time())
    ns = rand(1:(10^8))

    @test !LibArchive.atime_is_set(entry)
    LibArchive.set_atime(entry, t, ns)
    @test LibArchive.atime_is_set(entry)
    @test LibArchive.atime(entry) == t
    @test LibArchive.atime_nsec(entry) == ns
    LibArchive.unset_atime(entry)
    @test !LibArchive.atime_is_set(entry)

    @test !LibArchive.birthtime_is_set(entry)
    LibArchive.set_birthtime(entry, t, ns)
    @test LibArchive.birthtime_is_set(entry)
    @test LibArchive.birthtime(entry) == t
    @test LibArchive.birthtime_nsec(entry) == ns
    LibArchive.unset_birthtime(entry)
    @test !LibArchive.birthtime_is_set(entry)

    @test !LibArchive.ctime_is_set(entry)
    LibArchive.set_ctime(entry, t, ns)
    @test LibArchive.ctime_is_set(entry)
    @test LibArchive.ctime(entry) == t
    @test LibArchive.ctime_nsec(entry) == ns
    LibArchive.unset_ctime(entry)
    @test !LibArchive.ctime_is_set(entry)

    @test !LibArchive.mtime_is_set(entry)
    LibArchive.set_mtime(entry, t, ns)
    @test LibArchive.mtime_is_set(entry)
    @test LibArchive.mtime(entry) == t
    @test LibArchive.mtime_nsec(entry) == ns
    LibArchive.unset_mtime(entry)
    @test !LibArchive.mtime_is_set(entry)

    LibArchive.clear(entry)
    @test !LibArchive.atime_is_set(entry)
    @test !LibArchive.birthtime_is_set(entry)
    @test !LibArchive.ctime_is_set(entry)
    @test !LibArchive.mtime_is_set(entry)
    LibArchive.free(entry)
    @test_throws ErrorException LibArchive.atime_is_set(entry)
end

let
    # dev number
    entry = LibArchive.Entry()
    dev1 = UInt64(rand(UInt32))
    # There doesn't seem to be a portable way to convert between minor and
    # major dev_t and the full dev_t
    devmajor2 = UInt64(rand(UInt8))
    devminor2 = UInt64(rand(UInt8))

    @test !LibArchive.dev_is_set(entry)
    LibArchive.set_dev(entry, dev1)
    @test LibArchive.dev_is_set(entry)
    @test LibArchive.dev(entry) == dev1
    LibArchive.set_devmajor(entry, devmajor2)
    LibArchive.set_devminor(entry, devminor2)
    @test LibArchive.dev_is_set(entry)
    @test LibArchive.devmajor(entry) == devmajor2
    @test LibArchive.devminor(entry) == devminor2

    LibArchive.set_rdev(entry, dev1)
    @test LibArchive.rdev(entry) == dev1
    LibArchive.set_rdevmajor(entry, devmajor2)
    LibArchive.set_rdevminor(entry, devminor2)
    @test LibArchive.rdevmajor(entry) == devmajor2
    @test LibArchive.rdevminor(entry) == devminor2

    LibArchive.clear(entry)
    @test !LibArchive.dev_is_set(entry)
    LibArchive.free(entry)
end

let
    # file type
    entry = LibArchive.Entry()

    for ft in (LibArchive.FileType.MT, LibArchive.FileType.REG,
               LibArchive.FileType.LNK, LibArchive.FileType.SOCK,
               LibArchive.FileType.CHR, LibArchive.FileType.BLK,
               LibArchive.FileType.DIR, LibArchive.FileType.IFO)
        LibArchive.set_filetype(entry, ft)
        @test LibArchive.filetype(entry) == ft
    end

    LibArchive.free(entry)
end

is_unix() && let
    # fflags
    entry = LibArchive.Entry()
    @test_throws ArgumentError LibArchive.fflags_text(entry)

    LibArchive.set_fflags(entry, 1, 2)
    @test LibArchive.fflags(entry) == (1, 2)
    flags_txt = LibArchive.fflags_text(entry)
    @test !isempty(flags_txt)

    LibArchive.free(entry)

    entry2 = LibArchive.Entry()
    @test_throws ArgumentError LibArchive.fflags_text(entry2)
    LibArchive.set_fflags(entry2, flags_txt)
    @test LibArchive.fflags(entry2) == (1, 2)
    @test LibArchive.fflags_text(entry2) == flags_txt

    LibArchive.free(entry2)
end

let
    # ids / names
    entry = LibArchive.Entry()
    @test_throws ArgumentError LibArchive.gname(entry)
    @test_throws ArgumentError LibArchive.uname(entry)

    LibArchive.set_gid(entry, 2000)
    LibArchive.set_uid(entry, 2002)
    @test LibArchive.gid(entry) == 2000
    @test LibArchive.uid(entry) == 2002
    @test_throws ArgumentError LibArchive.gname(entry)
    @test_throws ArgumentError LibArchive.uname(entry)

    LibArchive.set_gname(entry, "group_name1")
    @test LibArchive.gname(entry) == "group_name1"
    if is_unix()
        LibArchive.set_gname(entry, "Group αβ")
        @test LibArchive.gname(entry) == "Group αβ"
    end

    LibArchive.set_uname(entry, "user_name1")
    @test LibArchive.uname(entry) == "user_name1"
    if is_unix()
        LibArchive.set_uname(entry, "User γδ")
        @test LibArchive.uname(entry) == "User γδ"
    end

    LibArchive.clear(entry)
    @test_throws ArgumentError LibArchive.gname(entry)
    @test_throws ArgumentError LibArchive.uname(entry)

    LibArchive.free(entry)
end

let
    # hardlink, pathname, sourcepath, symlink
    entry = LibArchive.Entry()

    @test_throws ArgumentError LibArchive.hardlink(entry)
    LibArchive.set_hardlink(entry, "hard_link1")
    @test LibArchive.hardlink(entry) == "hard_link1"
    if is_unix()
        LibArchive.set_hardlink(entry, "Hard Link α")
        @test LibArchive.hardlink(entry) == "Hard Link α"
    end
    LibArchive.clear(entry)
    @test_throws ArgumentError LibArchive.hardlink(entry)

    @test_throws ArgumentError LibArchive.pathname(entry)
    LibArchive.set_pathname(entry, "path_name2")
    @test LibArchive.pathname(entry) == "path_name2"
    if is_unix()
        LibArchive.set_pathname(entry, "Path Name β")
        @test LibArchive.pathname(entry) == "Path Name β"
    end
    LibArchive.clear(entry)
    @test_throws ArgumentError LibArchive.pathname(entry)

    @test_throws ArgumentError LibArchive.sourcepath(entry)
    LibArchive.set_sourcepath(entry, "source_path3")
    @test LibArchive.sourcepath(entry) == "source_path3"
    if is_unix()
        LibArchive.set_sourcepath(entry, "Source Path γ")
        @test LibArchive.sourcepath(entry) == "Source Path γ"
    end
    LibArchive.clear(entry)
    @test_throws ArgumentError LibArchive.sourcepath(entry)

    @test_throws ArgumentError LibArchive.symlink(entry)
    LibArchive.set_symlink(entry, "sym_link4")
    @test LibArchive.symlink(entry) == "sym_link4"
    if is_unix()
        LibArchive.set_symlink(entry, "Sym Link δ")
        @test LibArchive.symlink(entry) == "Sym Link δ"
    end
    LibArchive.clear(entry)
    @test_throws ArgumentError LibArchive.symlink(entry)

    LibArchive.free(entry)
end

let
    # ino and nlink
    entry = LibArchive.Entry()

    @test !LibArchive.ino_is_set(entry)
    LibArchive.set_ino(entry, 2345)
    @test LibArchive.ino(entry) == 2345
    LibArchive.set_nlink(entry, 10)
    @test LibArchive.nlink(entry) == 10

    LibArchive.clear(entry)
    @test !LibArchive.ino_is_set(entry)

    LibArchive.free(entry)
end

let
    # perm and mode
    entry = LibArchive.Entry()

    LibArchive.set_perm(entry, 0o644)
    @test LibArchive.perm(entry) == 0o644
    mode = LibArchive.mode(entry)
    strmode = LibArchive.strmode(entry)
    @test mode != 0
    @test !isempty(strmode)
    LibArchive.clear(entry)

    LibArchive.set_perm(entry, 0o600)
    @test LibArchive.perm(entry) == 0o600
    LibArchive.set_mode(entry, mode)
    @test LibArchive.perm(entry) == 0o644
    @test LibArchive.mode(entry) == mode
    @test LibArchive.strmode(entry) == strmode

    LibArchive.free(entry)
end

let
    # size
    entry = LibArchive.Entry()

    @test !LibArchive.size_is_set(entry)
    LibArchive.set_size(entry, 100)
    @test LibArchive.size_is_set(entry)
    @test LibArchive.size(entry) == 100
    LibArchive.unset_size(entry)
    @test !LibArchive.size_is_set(entry)

    LibArchive.free(entry)
end

# Create archive
info("Test creating and reading archive")
function create_archive(writer)
    entry = LibArchive.Entry(writer)
    LibArchive.set_pathname(entry, "test.txt")
    LibArchive.set_size(entry, 10)
    LibArchive.set_perm(entry, 0o644)
    LibArchive.set_filetype(entry, LibArchive.FileType.REG)
    LibArchive.write_header(writer, entry)
    write(writer, "012345678")
    write(writer, UInt8('9'))
    LibArchive.finish_entry(writer)

    entry = LibArchive.Entry(writer)
    LibArchive.set_pathname(entry, "test_a.txt")
    LibArchive.set_filetype(entry, LibArchive.FileType.LNK)
    LibArchive.set_symlink(entry, "test.txt")
    LibArchive.set_perm(entry, 0o755)
    LibArchive.write_header(writer, entry)
    LibArchive.finish_entry(writer)

    @test LibArchive.file_count(writer) == 2
end

function verify_archive(reader)
    entry = LibArchive.next_header(reader)
    @test LibArchive.pathname(entry) == "test.txt"
    @test LibArchive.size(entry) == 10
    @test LibArchive.perm(entry) == 0o644
    @test LibArchive.filetype(entry) == LibArchive.FileType.REG
    @test readstring(reader) == "0123456789"
    LibArchive.free(entry)

    entry = LibArchive.next_header(reader)
    @test LibArchive.pathname(entry) == "test_a.txt"
    @test LibArchive.filetype(entry) == LibArchive.FileType.LNK
    @test LibArchive.symlink(entry) == "test.txt"
    @test LibArchive.perm(entry) == 0o755
    LibArchive.free(entry)

    @test_throws EOFError LibArchive.next_header(reader)
    @test LibArchive.file_count(reader) == 2
end

info("    Filename")
mktempdir() do d
    cd(d) do
        LibArchive.Writer("./test.tar.bz2") do writer
            LibArchive.set_format_gnutar(writer)
            LibArchive.add_filter_bzip2(writer)
            LibArchive.set_bytes_per_block(writer, 4096)
            @test LibArchive.get_bytes_per_block(writer) == 4096
            LibArchive.set_bytes_in_last_block(writer, 1)
            @test LibArchive.get_bytes_in_last_block(writer) == 1
            create_archive(writer)
            @test LibArchive.format(writer) == LibArchive.Format.TAR_GNUTAR
            @test LibArchive.format_name(writer) == "GNU tar"
        end

        LibArchive.Reader("./test.tar.bz2") do reader
            LibArchive.support_filter_bzip2(reader)
            LibArchive.support_format_gnutar(reader)
            verify_archive(reader)
            @test LibArchive.filter_count(reader) > 0
            LibArchive.filter_bytes(reader, 0)
            @test (LibArchive.filter_code(reader, 0) ==
                   LibArchive.FilterType.BZIP2)
            @test LibArchive.filter_name(reader, 0) == "bzip2"
        end
    end
end

is_unix() && mktempdir() do d
    cd(d) do
        info("    FD")
        fd = ccall(:open, Cint, (Cstring, Cint, Cint),
                   "./test.tar.gz",
                   Compat.Filesystem.JL_O_WRONLY | Compat.Filesystem.JL_O_CREAT,
                   0o644)
        LibArchive.Writer(fd) do writer
            LibArchive.set_format_gnutar(writer)
            LibArchive.add_filter_gzip(writer)
            create_archive(writer)
        end
        ccall(:close, Cint, (Cint,), fd)

        fd = ccall(:open, Cint, (Cstring, Cint),
                   "./test.tar.gz", Compat.Filesystem.JL_O_RDONLY)
        LibArchive.Reader(fd) do reader
            LibArchive.support_filter_gzip(reader)
            LibArchive.support_format_gnutar(reader)
            verify_archive(reader)
        end
        ccall(:close, Cint, (Cint,), fd)
    end
end

info("    In memory")
let
    buffer = Vector{UInt8}(4096)
    used_size = LibArchive.Writer(buffer) do writer
        LibArchive.set_format_gnutar(writer)
        LibArchive.add_filter_bzip2(writer)
        create_archive(writer)
        # We need to flush the writer for `get_used` to return useful result
        close(writer)
        LibArchive.get_used(writer)
    end

    LibArchive.Reader(buffer, used_size) do reader
        LibArchive.support_filter_bzip2(reader)
        LibArchive.support_format_gnutar(reader)
        verify_archive(reader)
    end
end

info("    In memory (C pointer)")
let
    buffer = Libc.malloc(4096)
    used_size = LibArchive.Writer(buffer, 4096) do writer
        LibArchive.set_format_gnutar(writer)
        LibArchive.add_filter_bzip2(writer)
        create_archive(writer)
        # We need to flush the writer for `get_used` to return useful result
        close(writer)
        LibArchive.get_used(writer)
    end

    LibArchive.Reader(buffer, used_size) do reader
        LibArchive.support_filter_bzip2(reader)
        LibArchive.support_format_gnutar(reader)
        verify_archive(reader)
    end
    Libc.free(buffer)
end

info("    IO Stream")
let
    io = IOBuffer()
    LibArchive.Writer(io) do writer
        LibArchive.set_format_gnutar(writer)
        LibArchive.add_filter_none(writer)
        create_archive(writer)
    end

    seek(io, 0)
    LibArchive.Reader(io) do reader
        LibArchive.support_filter_none(reader)
        LibArchive.support_format_gnutar(reader)
        verify_archive(reader)
    end
end


text:
function line_search!(mechanism::Mechanism, α, rvio, bvio, opts)
    scale = 0
    system = mechanism.system

    rvio_cand, bvio_cand = Inf * ones(2)
    for n = Base.OneTo(opts.max_ls)
        for contact in mechanism.contacts
            candidate_step!(α, contact, get_entry(system, contact.id), scale)
        end
        for joint in mechanism.joints
            candidate_step!(α, joint, get_entry(system, joint.id), scale)
        end
        for body in mechanism.bodies
            ϕmax = 3.9 / mechanism.timestep^2
            candidate_step!(α, mechanism, body, get_entry(system, body.id), scale, ϕmax = ϕmax)
            if dot(body.state.ϕsol[2], body.state.ϕsol[2]) > 3.91 / mechanism.timestep^2
                error("Excessive angular velocity. Body-ID: $(string(body.name)) " * string(body.id) * ", ω: " * string(body.state.ϕsol[2]) * ".")
            end
        end

        rvio_cand = residual_violation(mechanism)
        bvio_cand = bilinear_violation(mechanism)

        if (rvio_cand > rvio) && (bvio_cand > bvio)
            scale += 1
        else
            return rvio_cand, bvio_cand
        end
    end

    return rvio_cand, bvio_cand
end

function cone_line_search!(mechanism::Mechanism;
    τort::T=0.95,
    τsoc::T=0.95,
    scaling::Bool=false) where T

    system = mechanism.system

    α = 1.0
    for contact in mechanism.contacts
        α = cone_line_search!(α, mechanism, contact, get_entry(system, contact.id), τort, τsoc; scaling = scaling)
    end
    for joint in mechanism.joints
        α = cone_line_search!(α, mechanism, joint, get_entry(system, joint.id), τort, τsoc; scaling = scaling)
    end

    return α
end

function cone_line_search!(α, mechanism, contact::ContactConstraint{T,N,Nc,Cs,N½},
        vector_entry::Entry, τort, τsoc;
        scaling::Bool=false) where {T,N,Nc,Cs<:NonlinearContact{T,N},N½}

    s = contact.impulses_dual[2]
    γ = contact.impulses[2]
    Δs = vector_entry.value[1:N½]
    Δγ = vector_entry.value[N½ .+ (1:N½)]
    αs_ort = positive_orthant_step_length(s[1:1], Δs[1:1], τ = τort)
    αγ_ort = positive_orthant_step_length(γ[1:1], Δγ[1:1], τ = τort)
    αs_soc = second_order_cone_step_length(s[2:4], Δs[2:4]; τ = τsoc)
    αγ_soc = second_order_cone_step_length(γ[2:4], Δγ[2:4]; τ = τsoc)

    return min(α, αs_soc, αγ_soc, αs_ort, αγ_ort)
end

function cone_line_search!(α, mechanism, contact::ContactConstraint{T,N,Nc,Cs,N½},
        vector_entry::Entry, τort, τsoc;
        scaling::Bool=false) where {T,N,Nc,Cs<:Union{ImpactContact{T,N},LinearContact{T,N}},N½}

    s = contact.impulses_dual[2]
    γ = contact.impulses[2]
    Δs = vector_entry.value[1:N½]
    Δγ = vector_entry.value[N½ .+ (1:N½)]

    αs_ort = positive_orthant_step_length(s, Δs, τ = τort)
    αγ_ort = positive_orthant_step_length(γ, Δγ, τ = τort)

    return min(α, αs_ort, αγ_ort)
end

function cone_line_search!(α, mechanism, joint::JointConstraint{T,N,Nc},
        vector_entry::Entry, τort, τsoc;
        scaling::Bool=false) where {T,N,Nc}

    for (i, element) in enumerate([joint.translational, joint.rotational])
        s, γ = split_impulses(element, joint.impulses[2][joint_impulse_index(joint,i)])
        Δs, Δγ = split_impulses(element,  vector_entry.value[joint_impulse_index(joint,i)])

        αs_ort = positive_orthant_step_length(s, Δs, τ = τort)
        αγ_ort = positive_orthant_step_length(γ, Δγ, τ = τort)
        α = min(α, αs_ort, αγ_ort)
    end

    return α
end

function positive_orthant_step_length(λ::AbstractVector{T}, Δ::AbstractVector{T};
    τ::T = 0.99) where T

    α = 1.0
    for i in eachindex(λ)
        if Δ[i] < 0 # safer
            α = min(α, - τ * λ[i] / Δ[i])
        end
    end

    return α
end

function second_order_cone_step_length(λ::AbstractVector{T}, Δ::AbstractVector{T};
        τ::T=0.99,
        ϵ::T=1e-14) where T

    # check Section 8.2 CVXOPT
    λ0 = λ[1]
    λ_λ = max(λ0^2 - λ[2:end]' * λ[2:end], 1e-25)
    if λ_λ < 0.0
        # @show λ_λ
        @warn "should always be positive"
    end
    λ_λ += ϵ
    λ_Δ = λ0 * Δ[1] - λ[2:end]' * Δ[2:end] + ϵ

    ρs = λ_Δ / λ_λ
    ρv = Δ[2:end] / sqrt(λ_λ)
    ρv -= (λ_Δ / sqrt(λ_λ) + Δ[1]) / (λ0 / sqrt(λ_λ) + 1) * λ[2:end] / λ_λ
    α = 1.0
    if norm(ρv) - ρs > 0.0
        α = min(α, τ / (norm(ρv) - ρs))
    end

    return α
end

function candidate_step!(α, mechanism::Mechanism, body::Body, vector_entry::Entry, scale; ϕmax = Inf)
    body.state.vsol[2] = body.state.vsol[1] + 1 / (2^scale) * α * vector_entry.value[SA[1; 2; 3]]
    body.state.ϕsol[2] = body.state.ϕsol[1] + 1 / (2^scale) * α * vector_entry.value[SA[4; 5; 6]]
    ϕ = body.state.ϕsol[2]
    ϕdot = dot(ϕ, ϕ)
    if ϕdot > ϕmax
        println("clipping ", scale, scn((ϕdot - ϕmax) / ϕmax), " ", scn(ϕdot), " ", scn(ϕmax), " ", body.name)
        body.state.ϕsol[2] *= ϕmax / ϕdot # this is overkill, but works better than sqrt(ϕmax/ϕdot)
    end
    return
end

function candidate_step!(α, joint::JointConstraint, vector_entry::Entry, scale)
    joint.impulses[2] = joint.impulses[1] + 1.0 / (2^scale) * α * vector_entry.value
    return
end

function candidate_step!(α::T, contact::ContactConstraint{T,N,Nc,Cs,N½}, vector_entry::Entry, scale) where {T,N,Nc,Cs,N½}
    contact.impulses_dual[2] = contact.impulses_dual[1] + 1 / (2^scale) * α * vector_entry.value[SVector{N½,Int64}(1:N½)]
    contact.impulses[2] = contact.impulses[1] + 1 / (2^scale) * α * vector_entry.value[SVector{N½,Int64}(N½+1:N)]
    return
end


text:
type FullyConnectedLayer <: Layer
    filters::Vector{Vol}
    biases::Vol
    inputActivation::Vol
    outputActivation::Vol

    function FullyConnectedLayer(numOutputs::Int, bias::Float64, dimensions)
        filters = [Vol(dimensions) for _ in 1:numOutputs]
        biases = Vol((1, 1, numOutputs), bias)
        new(filters, biases, Vol(dimensions, 0.0), Vol((1, 1, numOutputs), 0.0))
    end
end

function forward(self::FullyConnectedLayer, inputActivation::Vol)
    @assert size(self.inputActivation) == size(inputActivation)
    self.inputActivation = inputActivation
    for i in 1:length(self.filters)
        filter = self.filters[i]
        @assert size(self.inputActivation) == size(filter)
        self.outputActivation.w[1, 1, i] = sum(filter.w .* inputActivation.w) + self.biases.w[1, 1, i]
    end
    return self.outputActivation
end

function backward(self::FullyConnectedLayer)
    fill!(self.inputActivation.dw, 0.0)
    for i in 1:length(self.filters)
        filter = self.filters[i]
        chainGrad = self.outputActivation.dw[1, 1, i]

        self.inputActivation.dw += filter.w * chainGrad
        filter.dw += self.inputActivation.w * chainGrad
        self.biases.dw[1, 1, i] += chainGrad
        @printf("Filter %s, chainGrad: %s, inputActivation.w: %s\n", i, chainGrad, flatten(self.inputActivation.w))
    end
end

function gradientStates(self::FullyConnectedLayer)
    grads = [GradientState(filter, 1.0, 1.0) for filter in self.filters]
    append!(grads, [GradientState(self.biases, 0.0, 0.0)])
    return grads
end


text:
function average_temperature!(r::Robot)
    move_till_border!(r)
    move_by_field!(r::Robot)
end

function move_till_border!(r::Robot)
    while isborder(r, Nord) == false
        move!(r, Nord)
    end
    while isborder(r, West) == false
        move!(r, West)
    end
end

function move_by_field!(r::Robot)
    sum = 0
    k = 0
    while isborder(r, Nord) && isborder(r, Ost) == false
        while isborder(r, Sud) == false
            if ismarker(r) == true
                k = k + 1
                sum = sum + temperature(r)
            end
            move!(r, Sud)
        end
        if isborder(r, Ost) == false
            move!(r, Ost)
        end
        while isborder(r, Nord) == false
            move!(r, Nord)
            if ismarker(r) == true
                k = k + 1
                sum = sum + temperature(r)
            end
        end
        if isborder(r, Ost) == false
            if ismarker(r) == true
                k = k + 1
                sum = sum + temperature(r)
            end
            move!(r, Ost)
        end
    end
    if k > 0
        print(sum / k)
    else
        print("0")
    end
end


text:
import Gtk

function GTKdrawingwindow(name, w, h, closecb=nothing)
    win = Gtk.Window(name, w, h)
    c = Gtk.Canvas(win)
    if closecb !== nothing
        Gtk.on_signal_destroy(win, closecb, win)
    end
    c, win
end

_saved_gtk_renderer = nothing
function _saved_gtk_destroyed(::Ptr, widget::Gtk.GTKWidget)
    global _saved_gtk_renderer = nothing
    nothing
end
function gtk(self::PlotContainer, args...)
    global _saved_gtk_renderer, _saved_gtk_win
    opts = Winston.args2dict(args...)
    width = get(opts, :width, Winston.config_value("window","width"))
    height = get(opts, :height, Winston.config_value("window","height"))
    reuse_window = isinteractive() #&& Winston.config_value("window","reuse")
    device = _saved_gtk_renderer
    if device === nothing || !reuse_window
        device, win = GTKdrawingwindow("Julia", width, height, _saved_gtk_destroyed)
        _saved_gtk_renderer = device
    end
    display(device, self)
    self
end

function display(c::Gtk.Canvas, pc::PlotContainer)
    c.draw = function(_)
        ctx = getgc(c)
        set_source_rgb(ctx, 1, 1, 1)
        paint(ctx)
        Winston.page_compose(pc, Gtk.cairo_surface(c))
    end
    Gtk.draw(c)
end


text:
using LinearAlgebra, MAT, SparseArrays
using Statistics, HNSW, Distances, MatrixMarket
using Graphs, SimpleWeightedGraphs
using Arpack, Clustering

function undirectedEdges(idxs,dists) #O(nk)
    num_node = length(idxs) #each node
    k = length(idxs[1]) #nearest neighbors
    edges_dict = Dict{Tuple{Int,Int},Float64}()
    #edges = Vector{Tuple{Int,Int,Float64}}(undef,num_node*(k-1))
    idxs = convert(Array{Array{Int}},idxs) #idxs to integers
    for i = 1:num_node
        for j = 2:k
            if(dists[i][1]!=0.0)
                @warn "Node $i could not find itself: ef_construcion, M or ef in HNSW construction may be too low."
                a = min(i,idxs[i][j-1])
                b = max(i,idxs[i][j-1])
                edges_dict[(a,b)] = dists[i][j-1]
            else
                a = min(idxs[i][1],idxs[i][j])
                b = max(idxs[i][1],idxs[i][j])
                edges_dict[(a,b)] = dists[i][j]
            end
        end
    end
    sort(collect(edges_dict), by = x -> x[1])
end

function getLaplacian(G)
    A = adjacency_matrix(G)
    D = sparse(Diagonal(vec(sum(A;dims = 2))))
    L = D-A
    return L
end

function FiedlerAddEdges(G,fea,samplePerc,chosenPercEachIter,sig,embedTol)

    num_eigs_f = 5
    num_node,dim = size(fea)

    #=
    @info "laplacian_matrix"
    @time begin
    L = laplacian_matrix(G)
    end
    =#
    #@info "laplacian_matrix"
    #@time begin
    L = getLaplacian(G)
    #end

    L = L + (1/(sig^2))*sparse(LinearAlgebra.I,L.m,L.n)
    #mmwrite(raw"C:\Users\jacec\Desktop\code\code\JuliaL.mtx",L)
    #L = mmread(raw"C:\Users\jacec\Documents\MATLAB\GRASPEL_Code\code\JuliaL.mtx")
    vals, vecs = eigs(L; nev = num_eigs_f, which = :SM, tol = 1e-6)
    #@info vals
    #@info vecs[1,:]
    for i in 2:num_eigs_f
        vecs[:,i] = vecs[:,i]/sqrt(vals[i])
    end
    #vecs[:,2:num_eigs_f] = vecs[:,2:num_eigs_f].*[sqrt(i) for i in vals[2:num_eigs_f]]
    #@info vecs[1,:]

    #Fiedler Vector
    l2 = vecs[:,2]

    #Graph node sorting
    I = sortperm(l2)
    #@info I

    top = I[1:Int64(floor(num_node*samplePerc))]
    bot = I[num_node-Int64(floor(num_node*samplePerc)+1):num_node]
    chosenSet = zeros(Float64,Int64(ceil(num_node*chosenPercEachIter)),5)
    minDelta, location = findmin(chosenSet[:,3])

    #sample candidate edges

    for i in 1:length(chosenSet)*1000
        a_idx = rand(top);
        b_idx = rand(bot);
        #always set p to have the bigger node index
        p = max(a_idx, b_idx)
        q = min(a_idx, b_idx)
        #e_pq = zeros(num_node,1)       unsure what these were for
        #e_pq[p] = 1
        #e_pq[q] = -1
        distX = pairwise(SqEuclidean(),reshape(fea[p,:],:,1),reshape(fea[q,:],:,1))[1]
        weight = 1/distX
        distZ = pairwise(SqEuclidean(),reshape(vecs[p,:],:,1),reshape(vecs[q,:],:,1))[1]*dim
        delta_fiedler = distZ*weight
        #check if edge has been selected
        flagP = p in chosenSet[:,1]
        flagQ = q in chosenSet[:,2]
        if delta_fiedler > minDelta && !(flagP || flagQ) && delta_fiedler > embedTol
            chosenSet[location,1] = p
            chosenSet[location,2] = q
            chosenSet[location,3] = delta_fiedler
            chosenSet[location,4] = weight*1 #scale weight by up to 10
            gradient = (1-1/delta_fiedler)*distZ
            chosenSet[location,5] = gradient #gradient
            minDelta, location = findmin(chosenSet[:,3])
            #@info "P:$p Q:$q distX:$distX weight:$weight
            #distZ:$distZ delta:$delta_fiedler gradient:$gradient"
        end
    end

    #alternative cost formula
    if false
        lambda = diag(vals)
        cost = log(prod(lambda[2:num_eigs]))-trace(L*fea*transpose(fea))/dim
    else
        cost = maximum(chosenSet[:,5])
    end
    variationRatio = mean(chosenSet[:,3])

    if minimum(chosenSet[:,3]) == 0
        return G, variationRatio, cost
    end

    for k in 1:length(chosenSet[:,1])
        add_edge!(G,chosenSet[k,1],chosenSet[k,2],chosenSet[k,3])
    end

    return G, variationRatio, cost
end

#FiedlerConstruct
function FiedlerConstruct(fea, samplePerc, edge_iter, sig, maxIterations, k, tol, num_node, dim)
    varRatio = Array{Float64}(undef,maxIterations)
    costVec = Array{Float64}(undef,maxIterations)

    #Generate Initial KNN

    #Intialize HNSW struct
    feaVec = [vec(fea[i,:]) for i = 1:num_node] #input should be a vector of vectors
    @info "$k NN HNSW:..."
    @time begin
        #hnsw = HierarchicalNSW(feaVec; efConstruction=100, M=45, ef=90)
        hnsw = HierarchicalNSW(feaVec; efConstruction=100, M=45, ef=100)
        #Add all data points into the graph
        #Optionally pass a subset of the indices in data to partially construct the graph
        add_to_graph!(hnsw)

        # Find k+1 (approximate) nearest neighbors for each of the queries; first result is usually itself
        # results (# edges) are slightly different each time
        idxs, dists = knn_search(hnsw, feaVec, k+1)

        e = undirectedEdges(idxs,dists)
        #Normalize edge weights by w = 1/dist^2
        we = [Pair(e[i][1],1/(e[i][2]^2)) for i = 1:length(e)]
        #Store as sparse graph, ~14000-14005 nodes for USPS dataset
        G = SimpleWeightedGraph(num_node)
        for i in we
            add_edge!(G,i[1][1],i[1][2],i[2])
        end

    end
    @info "Processed Graph" G

    #MatrixMarket.mmwrite("JuliaL.mtx",L)

    @info "GRASPEL Iterations..."
    @time begin
        numIter = 0
        for i in 1:maxIterations
            numIter = i
            G, varR, cost = FiedlerAddEdges(G,fea,samplePerc,edge_iter,sig,tol)

            varRatio[i] = varR
            costVec[i] = cost

            if varR < tol
                break
            end

        end
    end
    return G, varRatio[1:numIter], costVec[1:numIter], numIter
end

max_stars_repo_name: janticev/GRASPEL

max_issues_repo_name: janticev/GRASPEL

max_forks_repo_name: janticev/GRASPEL


text:
using ExprOptimization, ExprRules, LinearAlgebra
using Test

let
    grammar = @grammar begin
        R = R + R
        R = |(1:3)
    end

    pcfg = ProbabilisticExprRules.ProbabilisticGrammar(grammar)
    @test ProbabilisticExprRules.probabilities(pcfg, :R) == normalize(ones(length(grammar[:R])), 1)

    dmap = mindepth_map(grammar)
    r = rand(RuleNode, pcfg, :R, dmap)

    iter = ExpressionIterator(grammar, 2, :R)
    pop = collect(iter)

    ProbabilisticExprRules.fit_mle!(pcfg, pop)
    @test all(isapprox.(ProbabilisticExprRules.probabilities(pcfg, :R), [0.3, 0.23333, 0.23333, 0.23333]; atol=0.001))

    ProbabilisticExprRules.uniform!(pcfg)
    @test ProbabilisticExprRules.probabilities(pcfg, :R) == normalize(ones(length(grammar[:R])), 1)
end


text:
using SparseArrays
using GraphNeuralNetworks
using BenchmarkTools
import Random: seed!
using LinearAlgebra

n = 1024
seed!(0)
A = sprand(n, n, 0.01)
b = rand(1, n)
B = rand(100, n)

g = GNNGraph(
    A,
    ndata=(; b=b, B=B),
    edata=(; A=reshape(A.nzval, 1, :)),
    graph_type=:coo  # changing to :sparse has little effect on performance
)

function spmv(g)
    propagate(
        (xi, xj, e) -> e .* xj ,  # same as e_mul_xj
        g, +; xj=g.ndata.b, e=g.edata.A
        )
end

function spmm1(g)
    propagate(
        (xi, xj, e) -> e .* xj ,  # same as e_mul_xj
        g, +; xj=g.ndata.B, e=g.edata.A
        )
end
function spmm2(g)
    propagate(
        e_mul_xj,
        g, +; xj=g.ndata.B, e=vec(g.edata.A)
        )
end

# @assert isequal(spmv(g),  b * A)  # true
# @btime spmv(g)  # ~5 ms
# @btime b * A  # ~32 us

@assert isequal(spmm1(g), B * A)  # true
@assert isequal(spmm2(g), B * A)  # true
@btime spmm1(g)  # ~9 ms
@btime spmm2(g)  # ~9 ms
@btime B * A  # ~400 us

function spmm_copyxj_fused(g)
    propagate(
        copy_xj,
        g, +; xj=g.ndata.B
        )
end

function spmm_copyxj_unfused(g)
    propagate(
        (xi, xj, e) -> xj,
        g, +; xj=g.ndata.B
        )
end

Adj = map(x -> x > 0 ? 1 : 0, A)
@assert spmm_copyxj_unfused(g) ≈ B * Adj
@assert spmm_copyxj_fused(g) ≈ B * Adj # bug fixed in #107

@btime spmm_copyxj_fused(g)  # 268.614 μs (22 allocations: 1.13 MiB)
@btime spmm_copyxj_unfused(g)  # 4.263 ms (52855 allocations: 12.23 MiB)
@btime B * Adj  # 196.135 μs (2 allocations: 800.05 KiB)

println()


text:
using Interact
using Plots

gr()

x = -π:0.05:π
y = sin.(x)

shift = slider(
0.01:0.1:10,
label="shift"
)

mp=@manipulate for s in shift
    @. y=sin(x - s)
    plot(x,y)
    ylims!(-1,1)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
export ScaledBeta, TrGeometric, TrBinomial, StringCategorical, likelihood, log_likelihood
export normalising_const
## Scaled beta distribution
"""
Scaled Beta distribution. Instantiate with ScaledBeta() e.g.
```
d = ScaledBeta(α, β, 0.0, 10.0)
```
would make a ScaledBeta type d which is bounded by 0.0 and 10.0. Note
- (α,β) are usual parameters;
- lb = lower bound;
- ub = upper bound.
"""
struct ScaledBeta<:ContinuousUnivariateDistribution
    α::Float64
    β::Float64
    lb::Float64
    ub::Float64
end

Distributions.rand(d::ScaledBeta) = d.lb + (d.ub-d.lb) * rand(Beta(d.α, d.β))
Distributions.mean(d::ScaledBeta) = d.lb + (d.ub-d.lb) * ( d.α / (d.α + d.β) )
Distributions.pdf(d::ScaledBeta, x::Real) = (d.lb ≤ x ≤ d.ub) ? pdf(Beta(d.α, d.β), (x - d.lb)/(d.ub-d.lb))/(d.ub-d.lb) : 0
Distributions.minimum(d::ScaledBeta) = d.lb
Distributions.maximum(d::ScaledBeta) = d.ub

## Trucated geometric
"""
`TrGeomtric(p::Float64, lb::Int, ub::Int)`

Truncated Geometric distribution. Note:
- lb = lower bound = samllest possible value (default=0)
- ub = upper bound = largest possible value.
Instantiation:

1. `TrGeometric(p, ub)` -  here `lb` defaults to 0
2. `TrGeometric(p, lb, ub)` here `lb` is specified

"""
struct TrGeometric<:DiscreteUnivariateDistribution
    p::Float64 # Prob param
    lb::Int # Lower bound (smallest value)
    ub::Int  # Upper bound
end

function normalising_const(d::TrGeometric)
    return (1 - (1-d.p)^(d.ub-d.lb+1))
end

TrGeometric(p::Float64, ub::Int) = TrGeometric(p, 0, ub)

function Distributions.pdf(d::TrGeometric, x::Integer)
    _insupport = insupport(d, x)
    return _insupport ? d.p*(1-d.p)^(x-d.lb) / normalising_const(d) : 0.0
end

function Distributions.logpdf(d::TrGeometric, x::Integer)
    _insupport = insupport(d, x)
    return _insupport ? log(d.p) + (x-d.lb) * log(1-d.p) -  log(normalising_const(d)) : -Inf
end

function Distributions.cdf(d::TrGeometric, x::Number)
    if d.lb ≤ x ≤ d.ub
        return ( 1 - (1 - d.p)^(trunc(Int,x)-d.lb+1) )/normalising_const(d)
    else
        return float(x > d.ub)
    end
end

# Here I make use of closed form of generalised inverse of the TrGeom dist.
# Not sure if this is the fastest implementation
function Distributions.quantile(d::TrGeometric, q::Real)
    return 0.0<q<1.0 ? ceil(Int, log(1-q*normalising_const(d))/log(1-d.p) - 1 )+d.lb : d.ub*(p≥1.0) + d.lb*(q≤0.0)
end

Distributions.minimum(d::TrGeometric) = d.lb
Distributions.maximum(d::TrGeometric) = d.ub

#Truncated Binomial wrappper - solves the problem of when we have n=0

struct TrBinomial <: DiscreteUnivariateDistribution
    n::Int
    p::Real
    lb::Int
    ub::Int
    dist::Distribution
    function TrBinomial(n::Int, p::Real, lb::Int, ub::Int)
        @assert (lb ≤ ub) "Lower bound must be ≤ upper bound."
        if lb < ub
            tmp_dist = Truncated(Binomial(n, p), lb, ub)
        else
            tmp_dist = DiscreteUniform(lb, ub)
        end
        new(n, p, lb, ub, tmp_dist)
    end
end

function Distributions.rand(d::TrBinomial)
    return rand(d.dist)
end
function Distributions.rand(d::TrBinomial, n::Int)
    return rand(d.dist, n)
end

struct StringCategorical <: DiscreteUnivariateDistribution
    x::Vector{String}
    p::Vector{T} where T<:Real
    latent_dist::Distribution
    probability_map::Dict{String, Float64}
    function StringCategorical(x::Vector{String}, p::Vector{T} where T<: Real)
        latent_dist = Categorical(p)
        probability_map = Dict(val => prob for (val, prob) in zip(x, p))
        new(x, p, latent_dist, probability_map)
    end
end

Distributions.pdf(d::StringCategorical, val::String) = d.probability_map[val]
function likelihood(d::StringCategorical, val::Vector{String})
    z = 1.0
    for s in val
        z *= pdf(d, s)
    end
    return z
end

Distributions.logpdf(d::StringCategorical, val::String) = log(d.probability_map[val])
function log_likelihood(d::StringCategorical, val::Vector{String})
    z = 0.0
    for s in val
        z += logpdf(d, s)
    end
    return z
end

function Distributions.rand(d::StringCategorical)
    return d.x[rand(d.latent_dist)]
end

function Distributions.rand(d::StringCategorical, n::Int)
    return d.x[rand(d.latent_dist, n)]
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie, ArchGDAL
GLMakie.activate!() # HIDE
let
    dataset = ArchGDAL.read("./data/fuji.tif")
    # just in case you want to know more about the file structure
    #ArchGDAL.read("fuji.tif") do ds
    #    println(ArchGDAL.gdalinfo(ds))
    #end
    dataset_read = ArchGDAL.read(dataset)
    topo = dataset_read[:,end:-1:1];

    fig, ax, pltobj = lines(topo[561,50:end-30], linewidth = 0.85,
        figure = (resolution = (800,1000),backgroundcolor = :black,),
        axis = (backgroundcolor = :black,))
    c = 560
    for i in 1:1:560
        lines!(ax, topo[i,50:end-30] .+ 15*c, linewidth = 0.85, overdraw = true,
            colorrange = (minimum(topo), maximum(topo)),
            color = topo[i,50:end-30], colormap = :turbo) # :linear_bmy_10_95_c78_n256
        c -= 1
    end
    text!(ax,"Fuji", position = (250, 0),color = :white, textsize = 24)
    text!(ax,"Visualization by",position=(0, -20),color = :white,textsize = 12)
    text!(ax,"@LazarusAlon", position = (0, -200),color = :white,textsize = 12)
    hidedecorations!(ax)
    hidespines!(ax)
    hidedecorations!(ax)
    hidespines!(ax)
    fig
    #save("FigTOPOFuji.png", fig, px_per_unit = 2) # HIDE
    save(joinpath(@__DIR__, "output", "topoLinesFuji.png"), fig, px_per_unit = 2.0) # HIDE
end

using Pkg # HIDE
Pkg.status(["GLMakie","ArchGDAL"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
#########################################################################
##############################  SNAP  ###################################
#########################################################################
function potential_energy(c::Configuration, p::SNAP)
    A = get_snap(c, p)
    p_e = dot(A[1, :]  ,  p.β )
    return p_e
end

function potential_energy(r::Vector{Configuration}, p::SNAP)
    n = length(r)
    pe = zeros(n)
    for i = 1:n
        A = get_snap(r[i], p)
        pe[i] = dot(A[1, :] , p.β)
    end
    return pe
end


text:
include("../examples/smiley_examples.jl")
using .SmileyExample22, .SmileyExample52, .SmileyExample54, .SmileyExample55

using BenchmarkTools
using ForwardDiff
using IntervalArithmetic
using IntervalRootFinding

import Random

const SUITE = BenchmarkGroup()

Random.seed!(0)  # Seed the RNG to get consistent results
tol = 1e-15

include("dietmar_ratz_functions.jl")

S = SUITE["Smiley"] = BenchmarkGroup()
for example in (SmileyExample22, SmileyExample52, SmileyExample54) #, SmileyExample55)
    s = S[example.title] = BenchmarkGroup()
    for method in (Newton, Krawczyk)
        s[string(method)] = @benchmarkable roots($(example.f), $(example.region), $method, $tol)
    end
end

S = SUITE["Rastigrin stationary points"] = BenchmarkGroup()

# Rastrigin function:
const A = 10

f(x, y) = 2A + x^2 - A*cos(2π*x) + y^2 - A*cos(2π*y)
f(X) = f(X...)

ForwardDiff.gradient(f, X::IntervalBox) = ForwardDiff.gradient(f, X.v)
∇f = X -> ForwardDiff.gradient(f, X)

L = 5.0
X = IntervalBox(-L..(L+1), 2)

for method in (Newton, Krawczyk)
    S[string(method)] = @benchmarkable roots($(∇(f)), $X, $method, 1e-5)
end

S = SUITE["Linear equations"] = BenchmarkGroup()

sizes = (2, 5, 10)

for n in sizes
    s = S["n = $n"] = BenchmarkGroup()
    A = Interval.(randn(n, n))
    b = Interval.(randn(n))

    s["Gauss seidel"] = @benchmarkable gauss_seidel_interval($A, $b)
    s["Gauss seidel contractor"] = @benchmarkable gauss_seidel_contractor($A, $b)
    s["Gauss elimination"] = @benchmarkable gauss_elimination_interval($A, $b)
end

S = SUITE["Dietmar-Ratz"] = BenchmarkGroup()
X = Interval(0.75, 1.75)

for (k, dr) in enumerate(dr_functions)
    s = S["Dietmar-Ratz $k"] = BenchmarkGroup()

    if k != 8  # dr8 is excluded as it has too many roots
        for method in (Newton, Krawczyk)
            s[string(method)] = @benchmarkable roots($dr, $X, $method, $tol)
        end
    end

    s["Automatic differentiation"] = @benchmarkable ForwardDiff.derivative($dr, $X)
    s["Slope expansion"] = @benchmarkable slope($dr, $X, $(mid(X)))
end


text:
using Pkg
println("Activating environment in $(pwd())...")
pkg"activate ."
println("Installing packages..."); flush(stdout);
pkg"instantiate"
pkg"precompile"

using WebIO
WebIO.install_jupyter_nbextension()

println("Done!")


text:
"""
```
plot_forecast_sequence(models, input_types, cond_types,
                       model_realized, input_type_realized, cond_type_realized,
                       class, forecast_var;
                       forecast_strings = fill("", length(models)),
                       forecast_string_realized = "",
                       title = string(forecast_var),
                       start_date::Date = quartertodate("2014-Q4"),
                       end_date::Date = quartertodate("2021-Q4"),
                       forecast_display_length = 4,
                       filepath::String = "")
```
Plot a sequence of forecasts (pulled from the tuple (model, input_type, cond_type, class, forecast_string)) on the same figure along with a realized series given by (model_realized, ...).
This is useful for seeing the evolution of forecasts with respect to the true realization of a series over time, e.g. is the DSGE model constantly under-forecasting realized inflation?

### Arguments
- `models::AbstractVector`: A vector of model objects with the correct vintage/saveroot settings from which to load the sequence of forecast MeansBands.
- `input_types::Vector{Symbol}`: A vector of the input types of each of the forecasts that were run.
- `cond_types::Vector{Symbol}`: A vector of the cond types of each of the forecasts that were run.
- `model_realized::AbstractDSGEModel`: The model object from which to load the historical realization/smoothed history of the series to compare to the forecasts.
- `input_type_realized::Symbol`: The input type from which the historical realization/smoothed history was generated.
- `cond_type_realized::Symbol`: The cond type from which the historical realization/smoothed history was generated.
- `class::Symbol`: The class, either :obs or :pseudo, that `forecast_var` belongs to.
- `forecast_var::Symbol`: The forecasted variable that is being plotted, e.g. :obs_gdp, or :NaturalRate.

### Keyword Arguments
- `forecast_strings::Vector{String}`: A vector of the forecast_strings of each of the forecasts that were run.
- `forecast_string_realized::String`: The forecast_string from which the historical realization/smoothed history was generated.
- `title::String`: The title.
- `start_date::Date`: The start date.
- `end_date::Date`: The end date.
- `forecast_display_length::String`: The number of periods of forecast that you want plotted.
- `filepath::String`: The filepath that you want to save the figure to.
"""
function plot_forecast_sequence(models::AbstractVector,
                                input_types::Vector{Symbol},
                                cond_types::Vector{Symbol},
                                model_realized::AbstractDSGEModel,
                                input_type_realized::Symbol,
                                cond_type_realized::Symbol,
                                class::Symbol,
                                forecast_var::Symbol;
                                forecast_strings::Vector{String} = fill("", length(models)),
                                forecast_string_realized::String = "",
                                title::String = string(forecast_var),
                                start_date::Date = quartertodate("2014-Q4"),
                                end_date::Date = quartertodate("2021-Q4"),
                                forecast_display_length::Int = 4,
                                plot_to_realized_end_date::Bool = false,
                                realized_linecolor::Symbol = :gray,
                                realized_linestyle::Symbol = :dash,
                                constant_to_plot::Float64 = Inf,
                                constant_linecolor::Symbol = :black,
                                ylims::Tuple = Tuple{}(),
                                filepath::String = "", kwargs...)
    n_models = length(models)
    @assert length(input_types) == length(cond_types) == length(forecast_strings) == n_models

    mb_realized = read_mb(model_realized, input_type_realized,
                          cond_type_realized, Symbol("hist", class),
                          forecast_string = forecast_string_realized)

    # Read in Forecast MBs
    mbs_forecast     = Vector{MeansBands}(undef, n_models)
    for i in 1:n_models
        mbs_forecast[i] = read_mb(models[i], input_types[i], cond_types[i],
                                  Symbol("forecast", class),
                                  forecast_string = forecast_strings[i])
    end

    plot_forecast_sequence(mbs_forecast, mb_realized, forecast_var;
                           title = title, start_date = start_date,
                           forecast_display_length = forecast_display_length,
                           plot_to_realized_end_date = plot_to_realized_end_date,
                           realized_linecolor = realized_linecolor, realized_linestyle = realized_linestyle,
                           constant_to_plot = constant_to_plot, constant_linecolor = constant_linecolor,
                           ylims = ylims, filepath = filepath)
end

function plot_forecast_sequence(mbs_forecast::Vector{MeansBands},
                                mb_realized::MeansBands,
                                forecast_var::Symbol;
                                title::String = string(forecast_var),
                                start_date::Date = quartertodate("2014-Q4"),
                                end_date::Date = quartertodate("2021-Q4"),
                                forecast_display_length::Int = 4,
                                plot_to_realized_end_date::Bool = false,
                                realized_linecolor::Symbol = :gray,
                                realized_linestyle::Symbol = :dash,
                                constant_to_plot::Float64 = Inf,
                                constant_linecolor::Symbol = :black,
                                ylims::Tuple = Tuple{}(),
                                filepath::String = "", kwargs...)
    n_forecasts = length(mbs_forecast)

    df_realized = mb_realized.means

    dfs_forecast = Vector{DataFrame}(undef, n_forecasts)
    for i in 1:n_forecasts
        dfs_forecast[i] = mbs_forecast[i].means
    end

    plot_forecast_sequence(dfs_forecast, df_realized, forecast_var;
                           title = title, start_date = start_date,
                           end_date = end_date, forecast_display_length = forecast_display_length,
                           plot_to_realized_end_date = plot_to_realized_end_date,
                           realized_linecolor = realized_linecolor, realized_linestyle = realized_linestyle,
                           constant_to_plot = constant_to_plot, constant_linecolor = constant_linecolor,
                           ylims = ylims, filepath = filepath)
end

function plot_forecast_sequence(dfs_forecast::Vector{DataFrame},
                                df_realized::DataFrame,
                                forecast_var::Symbol;
                                title::String = string(forecast_var),
                                start_date::Date = quartertodate("2014-Q4"),
                                end_date::Date = quartertodate("2021-Q4"),
                                forecast_display_length::Int = 4,
                                plot_to_realized_end_date::Bool = false,
                                realized_linecolor::Symbol = :gray,
                                realized_linestyle::Symbol = :dash,
                                constant_to_plot::Float64 = Inf,
                                constant_linecolor::Symbol = :black,
                                ylims::Tuple = Tuple{}(),
                                filepath::String = "")
    n_forecasts = length(dfs_forecast)

    # Construct realized series date indices
    realized_date_range_raw = Base.filter(x -> start_date <= x <= end_date, df_realized[:date])
    realized_date_range = map(quarter_date_to_number, realized_date_range_raw)
    date_inds  = findall(start_date .<= df_realized[:date] .<= end_date)

    # Plot realized
    p = plot(realized_date_range, df_realized[forecast_var][date_inds],
             label = "Realized", title = title, linecolor = realized_linecolor,
             linestyle = realized_linestyle, ylims = ylims)

    if !isinf(constant_to_plot)
        constant_start_ind   = realized_date_range[1]
        constant_end_ind     = realized_date_range[end]
    end

    # Plot forecasts
    for (i, df) in zip(1:n_forecasts, dfs_forecast)

        if plot_to_realized_end_date
            realized_end_date = df_realized[:date][end]

            # If `df[:date][1]` (the first forecast quarter)
            # is less than `forecast_display_length` away from the
            # realized_end_date, then plot the forecast up to
            # `forecast_display_length` number of periods.
            # Else, plot the forecast up to realized_end_date.
            # E.g. For forecasts that originate close to (or past) realized_end_date,
            # we still want to see a full forecast, but for forecasts that originate
            # far earlier than the current realized_end_date, we want to see their
            # forecast up to the realized_end_date.
            if subtract_quarters(realized_end_date, df[:date][1]) < forecast_display_length
                # Construct forecast date indices
                forecast_date_range_raw = Base.filter(x -> start_date <= x <= end_date,
                                                      dfs_forecast[i][:date])
                forecast_date_range = map(quarter_date_to_number, forecast_date_range_raw)

                # Truncate
                forecast_date_range = forecast_date_range[1:min(forecast_display_length, length(forecast_date_range))]
                date_inds = findall(start_date .<= dfs_forecast[i][:date] .<= end_date)
                date_inds = date_inds[1:min(forecast_display_length, length(forecast_date_range))]
            else
                # Construct forecast date indices
                forecast_date_range_raw = Base.filter(x -> start_date <= x <= realized_end_date,
                                                      dfs_forecast[i][:date])
                forecast_date_range = map(quarter_date_to_number, forecast_date_range_raw)

                date_inds = findall(start_date .<= dfs_forecast[i][:date] .<= realized_end_date)
            end
        else
            # Construct forecast date indices
            forecast_date_range_raw = Base.filter(x -> start_date <= x <= end_date,
                                                  dfs_forecast[i][:date])
            forecast_date_range = map(quarter_date_to_number, forecast_date_range_raw)

            # Truncate
            forecast_date_range = forecast_date_range[1:min(forecast_display_length, length(forecast_date_range))]
            date_inds = date_inds[1:min(forecast_display_length, length(date_inds))]
        end

        if i == 1
            label = "Forecast"
        else
            label = ""
        end

        plot!(p, forecast_date_range, dfs_forecast[i][forecast_var][date_inds],
              label = label, linecolor = :red, legend = :bottomright)

        if !isinf(constant_to_plot)
            constant_end_ind = max(constant_end_ind, forecast_date_range[end])
        end
    end

    # Plot the constant
    if !isinf(constant_to_plot)
        plot!(p, constant_start_ind:constant_end_ind,
              fill(constant_to_plot, length(constant_start_ind:constant_end_ind)),
              linecolor = constant_linecolor, label = "")
    end

    if !isempty(filepath)
        savefig(p, filepath)
        println("Saved $filepath")
    end

    return p
end


text:
A = Set([2,7,2,3])
B = Set(1:6)
omega = Set(1:10)

AunionB = union(A, B)
AintersectionB = intersect(A, B)
BdifferenceA = setdiff(B,A)
Bcomplement = setdiff(omega,B)
AsymDifferenceB = union(setdiff(A,B),setdiff(B,A))
println("A = $A, B = $B")
println("A union B = $AunionB")
println("A intersection B = $AintersectionB")
println("B diff A = $BdifferenceA")
println("B complement = $Bcomplement")
println("A symDifference B = $AsymDifferenceB")
println("The element '6' is an element of A: $(in(6,A))")
println("Symmetric difference and intersection are subsets of the union: ",
        issubset(AsymDifferenceB,AunionB),", ", issubset(AintersectionB,AunionB))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
mutable struct PerfectBipartiteMatchingAlgosSolver <: PerfectBipartiteMatchingSolver
  # All other instances only work with complete bipartite graphs, do the same here.
  n_arms::Int

  function PerfectBipartiteMatchingAlgosSolver()
    new(-1)
  end
end

has_lp_formulation(::PerfectBipartiteMatchingAlgosSolver) = false
approximation_ratio(::PerfectBipartiteMatchingAlgosSolver) = 1.0
approximation_term(::PerfectBipartiteMatchingAlgosSolver) = 0.0
approximation_ratio_budgeted(::PerfectBipartiteMatchingAlgosSolver) = 0.5
approximation_term_budgeted(::PerfectBipartiteMatchingAlgosSolver) = 0.0
supports_solve_budgeted_linear(::PerfectBipartiteMatchingAlgosSolver) = true
supports_solve_all_budgeted_linear(::PerfectBipartiteMatchingAlgosSolver) = true

function build!(solver::PerfectBipartiteMatchingAlgosSolver, n_arms::Int)
  solver.n_arms = n_arms
  nothing
end

function solve_linear(solver::PerfectBipartiteMatchingAlgosSolver, rewards::Dict{Tuple{Int, Int}, Float64})
  g = complete_bipartite_graph(solver.n_arms, solver.n_arms)
  r = Dict{Edge{Int}, Float64}(Edge(k[1], solver.n_arms + k[2]) => v for (k, v) in rewards)
  i = BipartiteMatchingInstance(g, r)

  s = matching_hungarian(i).solution

  return [(src(e), dst(e) - solver.n_arms) for e in s]
end

function solve_budgeted_linear(solver::PerfectBipartiteMatchingAlgosSolver,
                               rewards::Dict{Tuple{Int, Int}, Float64},
                               weights::Dict{Tuple{Int, Int}, Int},
                               budget::Int)
  g = complete_bipartite_graph(solver.n_arms, solver.n_arms)
  r = Dict{Edge{Int}, Float64}(Edge(k[1], solver.n_arms + k[2]) => v for (k, v) in rewards)
  w = Dict{Edge{Int}, Int}(Edge(k[1], solver.n_arms + k[2]) => v for (k, v) in weights)
  i = BudgetedBipartiteMatchingInstance(g, r, w, budget)

  s = matching_hungarian_budgeted_lagrangian_approx_half(i).solution

  return [(src(e), dst(e) - solver.n_arms) for e in s]
end

function solve_all_budgeted_linear(solver::PerfectBipartiteMatchingAlgosSolver,
                                   rewards::Dict{Tuple{Int, Int}, Float64},
                                   weights::Dict{Tuple{Int, Int}, Int},
                                   max_budget::Int)
  g = complete_bipartite_graph(solver.n_arms, solver.n_arms)
  r = Dict{Edge{Int}, Float64}(Edge(k[1], solver.n_arms + k[2]) => v for (k, v) in rewards)
  w = Dict{Edge{Int}, Int}(Edge(k[1], solver.n_arms + k[2]) => v for (k, v) in weights)
  i = BudgetedBipartiteMatchingInstance(g, r, w, max_budget)

  s = matching_dp_budgeted(i).solution

  sol = Dict{Int, Vector{Edge{Int}}}()
  vl = last(instance.matching.vertex_left)
  vr = last(instance.matching.vertex_right)
  for budget in 0:max_budget
    sol[budget] = s.solutions[vl, vr, 0, budget]
  end
  return sol
end


text:
# Search for things in modules.

"""
Return the set of all submodules of a given module `mod`.
"""
function submodules(mod::Module)
    out = Set([mod])
    for name in names(mod, true)
        if isdefined(mod, name)
            object = getfield(mod, name)
            validmodule(mod, object) && union!(out, submodules(object))
        end
    end
    out
end

validmodule(mod::Module, object::Module) = object ≠ mod && object ≠ Main
validmodule(::Module, other)             = false

"""
Which source files are known to be included in a module.

**Note:** files with only constants will not be detected.
"""
function includedfiles(mod::Module, candidates::Set)
    out = Set{String}()
    for name in names(mod, true)
        if isdefined(mod, name)
            object = getfield(mod, name)
            if isgeneric(object)
                for def in methods(object)
                    file = location(def)
                    if Utilities.samemodule(mod, def) && file ∈ candidates
                        push!(out, file)
                    end
                end
            elseif isa(object, Function) && isdefined(object, :code)
                file = location(object)
                if Utilities.samemodule(mod, object) && file ∈ candidates
                    push!(out, file)
                end
            end
        end
    end
    out
end

"""
Path to definition of a julia object, only methods are searched for.
"""
location(object::Method) = Utilities.expandpath(string(Utilities.lsdfield(object, :file)))
location(func::Function) = Utilities.expandpath(string(Utilities.lsdfield(func,   :file)))
location(other)          = ""

"""
Is the file the root for a module `mod`. Check for `Expr(:module, ...)`.

**Note:** Assumes that all submodules have unique names.
"""
function isrootfile(mod::Symbol, parsed::Expr)
    result = false
    for arg in parsed.args
        result = isexpr(arg, :module) && arg.args[2] ≡ mod ?
            true : isrootfile(mod, arg)
        result && break
    end
    result
end
isrootfile(mod::Symbol, other) = false

function matching(pred::Function, root::AbstractString)
    out = Set{String}()
    for object in readdir(root)
        path = joinpath(root, object)
        pred(path)  && push!(out, path)
        isdir(path) && union!(out, matching(pred, path))
    end
    out
end


text:
# ENV["PYCALL_JL_RUNTIME_PYTHON"] = Sys.which("python") or ''

#using PyCall
using ElectronicStructure: Atom, Geometry, MolecularSpec,
    InteractionOperator, MolecularData

# using ElectronicStructurePySCF: PySCF

using QuantumOps: FermiSum

"""
    qiskit_geometry_to_Geometry(geometry::Matrix)

Convert a geometry specification that originated in Python in the qiskit-nature format to
an `ElectronicStructure.Geometry` object. pyjulia will have translated the python input to
a `Matrix`.

# Example from python:
```python
In [1]: geometry = [['O', [0., 0., 0.]],
   ...:             ['H', [0.757, 0.586, 0.]],
   ...:             ['H', [-0.757, 0.586, 0.]]]

In [2]: Main.qiskit_geometry_to_Geometry(geometry)
Out[2]: <PyCall.jlwrap Geometry{Float64}(Atom{Float64}[Atom{Float64}(:O, (0.0, 0.0, 0.0)), Atom{Float64}(:H, (0.757, 0.586, 0.0)), Atom{Float64}(:H, (-0.757, 0.586, 0.0))])>
```
"""
function qiskit_geometry_to_Geometry(geometry::Matrix)
    (num_atoms, num_components_per_atom) = size(geometry)
    if num_components_per_atom != 2
        throw(ArgumentError("Expecting two elements to specify geometry of a single Atom. Got ", num_components_per_atom))
    end
    atoms = []
    for atom_in in eachrow(geometry)
        element_symbol = atom_in[1]
        coords = atom_in[2]
        atom = Atom(Symbol(element_symbol), (coords...,))
        push!(atoms, atom)
    end
    return Geometry(atoms...)
end

function fermionic_hamiltonian(geometry::Geometry, basis)
    # Construct the specification of the electronic structure problem.
    mol_spec = MolecularSpec(geometry=geometry, basis=basis)

    # Do calculations using `PySCF` and populate a `MolecularData` object with results.
    # `mol_pyscf` holds a constant, a rank-two tensor, and a rank-four tensor.
    mol_pyscf = MolecularData(PySCF, mol_spec);
    iop = InteractionOperator(mol_pyscf);
    return FermiSum(iop)
end

# Not useful. Just call jordan_wigner(fermi_op) directly
# function qubit_hamiltonian(fermi_op)
#     pauli_op = jordan_wigner(fermi_op)
#     return pauli_op
# end

# function qubit_hamiltonian(geometry::Geometry, basis)
#     fermi_op = fermionic_hamiltonian(geometry::Geometry, basis)
#     pauli_op = jordan_wigner(fermi_op)
#     return pauli_op
# end

max_stars_repo_name: jlapeyre/qiskit-alt

max_issues_repo_name: jlapeyre/qiskit-alt

max_forks_repo_name: Qiskit-Extensions/qiskit-alt


text:
using LinearAlgebra

"""
    ABCD = stuffABCD(a, g, b, c, form="CRFB")

Compute the ABCD matrix for the specified structure.

See realizeNTF for a list of supported structures.

mapABCD is the inverse function.
"""
function stuffABCD(a, g, b, c, form="CRFB")
    order = length(a)
    odd = rem(order, 2)
    even = 1 - odd
    ABCD = zeros(order+1, order+2)
    if length(b) == 1
        b = [b; zeros(order)]
    end

    if form == "CRFB"
        #C=(0 0...c_n)
        # This is done as part of the construction of A, below
        #B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD[:,order+1] = b
        #B2 = -(a_1 a_2... a_n)
        ABCD[1:order,order+2] = -a
        diagonal = 1:(order+2):order*(order+1)
        ABCD[diagonal] = ones(order)
        subdiag = diagonal[(1+even):2:order] .+ 1
        ABCD[subdiag]= c[1+even:2:order]
        supdiag = subdiag[(1+odd):length(subdiag)] .- 2
        ABCD[supdiag] = -g
        dly = (2+odd):2:order   # row numbers of delaying integrators
        ABCD[dly,:] = ABCD[dly,:] + diagm(c[dly.-1]) * ABCD[dly.-1,:]

    elseif form == "CRFF"
        #B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD[:,order+1] = b
        #B2 = -(c_1 0... 0)
        ABCD[1,order+2] = -c[1]
        diagonal = 1:(order+2):order*(order+1)  # num of elements = order
        ABCD[diagonal] = ones(order)
        subdiag = diagonal[1:2:order-1] .+ 1
        ABCD[subdiag] = c[2:2:order]
        if even == 1
            multg = 1:2:order   # rows to have g*(following row) subtracted.
            ABCD[multg,:] = ABCD[multg,:] - diagm(g) * ABCD[multg.+1,:]
        elseif order >= 3
            supdiag = diagonal[3:2:order] .- 1
            ABCD[supdiag] = -g
        end
        multc = 3:2:order   # rows to have c*(preceding row) added.
        ABCD[multc,:] = ABCD[multc,:] + diagm(c[multc]) * ABCD[multc.-1,:]
        ABCD[order+1,1:2:order] = a[1:2:order]
        for i = 2:2:order
            ABCD[order+1,:] = ABCD[order+1,:] + a[i] * ABCD[i,:]
        end

    elseif form == "CIFB"
        #C=(0 0...c_n)
        # This is done as part of the construction of A, below
        #B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD[:,order+1] = b
        #B2 = -(a_1 a_2... a_n)
        ABCD[1:order,order+2] = -a
        diagonal = 1:(order+2):order*(order+1)
        ABCD[diagonal] = ones(order)
        subdiag = diagonal[1:order] .+ 1
        ABCD[subdiag] = c
        supdiag = diagonal[(2+odd):2:order] .- 1
        ABCD[supdiag] = -g

    elseif form == "CIFF"
        #B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD[:,order+1] = b
        #B2 = -(c_1 0... 0)
        ABCD[1,order+2] = -c[1]
        diagonal = 1:(order+2):order*(order+1)
        ABCD[diagonal] = ones(order)
        subdiag = diagonal[1:order-1] .+ 1
        ABCD[subdiag]= c[2:end]
        #C = (a_1 a_2... a_n)
        ABCD[order+1,1:order] = a[1:order]
        supdiag = diagonal[(2+odd):2:order] .- 1
        ABCD[supdiag] = -g

    elseif form == "CRFBD"
        #C=(0 0...c_n)
        ABCD[order+1,order] = c[order]
        #B1 = (b_1 b_2... b_n), D=(b_n+1 0)
        ABCD[:,order+1] = b
        #B2 = -(a_1 a_2... a_n)
        ABCD[1:order,order+2] = -a
        diagonal = 1:(order+2):order*(order+1)
        ABCD[diagonal] = ones(order)
        dly = (1+odd):2:order   # row numbers of delaying integrators
        subdiag = diagonal[dly] .+ 1
        ABCD[subdiag]= c[dly]
        supdiag = subdiag[(1+odd):length(subdiag)] .- 2
        ABCD[dly,:] = ABCD[dly,:] - diagm(g) * ABCD[dly.+1,:]
        if order > 2
            coupl = 2+even:2:order
            ABCD[coupl,:] = ABCD[coupl,:] + diagm(c[coupl.-1]) * ABCD[coupl.-1,:]
        end

    elseif form == "CRFFD"
        diagonal = 1:(order+2):order*(order+1)
        subdiag = diagonal[1:order-1] .+ 1
        supdiag = diagonal[2:order] .- 1
        #B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD[:,order+1] = b
        #B2 = -(c_1 0... 0)
        ABCD[1,order+2] = -c[1]
        ABCD[diagonal] = ones(order)
        multc = 2:2:order   # rows to have c*(preceding row) added.
        if order > 2
            ABCD[subdiag[2:2:end]] = c[3:2:end]
        end
        if even == 1
            ABCD[supdiag[1:2:end]] = -g
        else
            # subtract g*(following row) from the multc rows
            ABCD[multc,:] = ABCD[multc,:] - diagm(g) * ABCD[multc.+1,:]
        end
        ABCD[multc,:] = ABCD[multc,:] + diagm(c[multc]) * ABCD[multc.-1,:]
        # C
        ABCD[order+1,2:2:order] = a[2:2:order]
        for i = 1:2:order
            ABCD[order+1,:] = ABCD[order+1,:] + a[i] * ABCD[i,:]
        end
        # The above gives y(n+1) need to add a delay to get y(n).
        # Do this by augmenting the states. Note: this means that
        # the apparent order of the NTF is one higher than it acually is.
        A, B, C, D = partitionABCD(ABCD, 2)
        A = [A zeros(order); C 0]
        B = [B; D]
        C = [zeros(1,order) 1]
        D = [0 0]
        ABCD = [A B; C D]

    elseif form == "PFF"
        #B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD[:,order+1] = b
        odd_1 = odd # !! Bold assumption !!
        odd_2 = 0   # !! Bold assumption !!
        gc = g .* c[1+odd_1:2:end]
        theta = acos.(1 .- gc/2)
        if odd_1 == 1
            theta0 = 0
        else
            theta0 = theta[1]
        end
        order_2 = 2 * length(findall(x -> abs(x-theta0) > 0.5, theta))
        order_1 = order - order_2
        #B2 = -(c_1 0...0 c_n 0...0)
        ABCD[1,order+2] = -c[1]
        ABCD[order_1+1,order+2] = -c[order_1+1]
        diagonal = 1:(order+2):order*(order+1)	# # of elements = order
        ABCD[diagonal] = ones(order)
        i = [1:2:order_1-1; order-order_2+1:2:order]
        subdiag = diagonal[i] .+ 1
        ABCD[subdiag] = c[i.+1]
        if odd_1 == 1
            if order_1 >= 3
                supdiag = diagonal[3:2:order_1] .- 1
                ABCD[supdiag] = -g[1:(order_1-1)÷2]
            end
        else
            multg = 1:2:order_1 # rows to have g*(following row) subtracted.
            ABCD[multg,:] = ABCD[multg,:] - diagm(g[1:order_1÷2]) * ABCD[multg.+1,:]
        end
        if odd_2 == 1
            if order_2 >= 3
                supdiag = diagonal[order_1+2:2:order] .- 1
                ABCD[supdiag] = -g[1:(order_1-1)÷2]
            end
        else
            multg = order_1+1:2:order # rows to have g*(following row) subtracted.
            gg = g[(order_1-odd_1)÷2+1:end]
            ABCD[multg,:] = ABCD[multg,:] - diagm(gg) * ABCD[multg.+1,:]
        end
        # Rows to have c*(preceding row) added.
        multc = [3:2:order_1; order_1+3:2:order]
        ABCD[multc,:] = ABCD[multc,:] + diagm(c[multc]) * ABCD[multc.-1,:]
        # C portion of ABCD
        i = [1:2:order_1; order_1+1:2:order]
        ABCD[order+1,i] = a[i]
        for i = [2:2:order_1; order_1+2:2:order]
            ABCD[order+1,:] = ABCD[order+1,:] + a[i] * ABCD[i,:]
        end

    elseif form == "Stratos"
        # code copied from if form == "CIFF":
        #B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD[:,order+1] = b
        #B2 = -(c_1 0... 0)
        ABCD[1,order+2] = -c[1]
        diagonal = 1:(order+2):order*(order+1)
        ABCD[diagonal] = ones(order)
        subdiag = diagonal[1:order-1] .+ 1
        ABCD[subdiag] = c[2:end]
        # code based on if form == "CRFF":
        multg = 1+odd:2:order-1	# rows to have g*(following row) subtracted.
        ABCD[multg,:] = ABCD[multg,:] - diagm(g) * ABCD[multg.+1,:]
        # code copied from if form== 'CIFF':
        #C = (a_1 a_2... a_n)
        ABCD[order+1,1:order] = a[1:order]

    elseif form == "DSFB"
        A0 = diagm(ones(order))
        A1 = zeros(order, order)
        indices = 2:order+1:(order-1)*order # subdiagonal
        A1[indices] = c[1:order-1]
        indices = order+1 : 2*(order+1) : order*order-1 #supdiagonal
        A1[indices] = -g
        C0 = diagm(ones(order))
        C1 = diagm(ones(order))
        D0 = zeros(order, 2)
        if odd == 1
            C1[order,order] = 0
        else
            A1[order-1:order,order-1:order] = [0 -g[end]/2; c[end-1] 0]
            C0[order,order] = 0
            D0[order-1,2] = g[end] / 2
        end
        M = inv(I - A1*C1)
        A = M * (A0 + A1*C0)
        B = M * ([b[1:order] -a] + A1*D0)
        C = [zeros(1,order-1) c[order]]
        D = [b[order+1] 0]
        ABCD = [A B; C D]

    else
        println("Unsupported form $form")
    end

    return ABCD
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
using Revise
using Flipping
##
Directory_path = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Stimulations/Rbp4_gen_control"
Directory_path = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Stimulations/WT_new_protocols"
pokes, streaks, DataIndex  = Flipping.create_exp_dataframes(Directory_path)
d = Flipping.find_behavior(Directory_path)
names(pokes)
###
filepath = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Stimulations/WT_new_protocols/WJ16_200105a.csv"
Flipping.session_info(filepath)
findmax(streaks[:,:Day])
curr_data = FileIO.load(filepath) |> DataFrame
names(curr_data)
DataFrames.rename!(curr_data, Symbol("") => :Poke)
##
pokes, streaks, DataIndex = Flipping.create_exp_dataframes(d[1:2,:])
###
###


text:
@testset "DiagonalNormal API: Reals" begin
    dn = TMVDiagonalNormal([0, 0], [1, 1])
    dn = TMVDiagonalNormal([0., 0.], [1., 1.])
    @test length(dn) == 2
    @test size(dn) == (2,)
    @test eltype(dn) == Float64
    @test mean(dn) == [0, 0]
    @test var(dn) == fill(exp(2), 2)
    @test cov(dn) == diagm(0 => fill(exp(2), 2))
    @test cor(dn) == diagm(0 => ones(2))
    mv = convert(Distributions.MvNormal, dn)
    @test mean(mv) == zeros(2)
    @test cov(mv) ≈ Diagonal(fill(exp(2), 2))
    @test size(rand(dn, 3)) == (2, 3)
    rng = MersenneTwister(23)
    rng2 = MersenneTwister(23)
    @test rand(rng, dn, 3) ≈  rand(rng2, mv, 3)

    rng = MersenneTwister(23)
    @test size(rand(rng, dn)) == (2,)
    # Gaussian evaluated at peak
    max_gaussian = prod(1 ./ sqrt.(2π * var(dn)))
    # gaussian evaluated at [1, 1]
    value = max_gaussian * prod(exp.(- [1, 1]/(2 * var(dn)[1])))
    expected_vec = [max_gaussian, value]
    @test pdf(dn, [0, 0]) ≈ max_gaussian
    @test pdf(dn, [[0 1]; [0 1]]) ≈ expected_vec
    @test logpdf(dn, [0, 0]) == log(max_gaussian)
    @test logpdf(dn, [[0 1]; [0 1]]) == log.(expected_vec)
    @test entropy(dn) == entropy(mv)
    @test loglikelihood(dn, [[0 1]; [0 1]]) == sum(log.(expected_vec))

    @test_throws AssertionError TMVDiagonalNormal(zeros(2,2), ones(2,2))
    @test typeof(TMVDiagonalNormal(ones(2)[:, :], ones(2)[:, :])) == typeof(TMVDiagonalNormal(ones(2), ones(2)))

    # Verify that TMVDiagonalNormal is treated as a scalar during broadcast
    dn = TMVDiagonalNormal([0, 0], [1, 1])
    @test size(dn .== dn) == ()
end

@testset "DiagonalNormal API: TrackedArray" begin
    dn = TMVDiagonalNormal(TrackedArray([0., 0]), TrackedArray([1., 1]))
    @test length(dn) == 2
    @test size(dn) == (2,)
    @test eltype(dn) == Float64
    @test mean(dn) == TrackedArray([0, 0])
    @test isa(mean(dn), TrackedArray)
    @test var(dn) == exp(2) * TrackedArray([1, 1])
    @test isa(var(dn), TrackedArray)
    @test cov(dn) == TrackedArray(diagm(0 => fill(exp(2), 2)))
    @test isa(cov(dn), TrackedArray)
    @test cor(dn) == Diagonal(ones(2))
    rng = MersenneTwister(23)
    @test size(rand(rng, dn)) == (2,)
    @test logpdf(dn, [0, 0]) == logpdf(MvNormal(zeros(2), exp.(ones(2))), zeros(2))
    @test_throws MethodError (logpdf(dn, [[0 1]; [0 1]]))

    @test pdf(dn, [0, 0]) == pdf(MvNormal(zeros(2), exp.(ones(2))), zeros(2))
    @test_throws MethodError pdf(dn, [[0 1]; [0 1]])
    @test_throws MethodError convert(Distributions.MvNormal, dn)
    @test_throws ErrorException rand(rng, dn, 1)
    @test_throws ErrorException rand(dn, 3)

    @test_throws MethodError entropy(dn)
    @test_nowarn loglikelihood(dn, [[0 1]; [0 1]])

    dn = TMVDiagonalNormal(TrackedArray([0, 0]), TrackedArray([1, 1]))
    rng = MersenneTwister(23)
    @test size(rand(rng, dn)) == (2,)

    # Check comparison with previous implementations:
    rng1 = MersenneTwister(1)
    rng2 = MersenneTwister(1)
    @test TrackedDistributions.sample(rng1, dn)==rand(rng2, dn)
    x = rand(2)
    @test TrackedDistributions.log_pdf(dn, x) == logpdf(dn, x)

    @test isa(data(dn), TrackedDistributions.TMVDiagonalNormal{Array{Float64,1}})

    # Verify that TMVDiagonalNormal is treated as a scalar during broadcast
    dn = TMVDiagonalNormal(TrackedArray([0., 0]), TrackedArray([1., 1]))
    @test size(dn .== dn) == ()
    end

@testset "DiagonalNormal API: Univariate" begin
    # There are no univariates! Convert to multivariate. Only works for reals.
    @test typeof(TMVDiagonalNormal([0], [1]))==typeof(TMVDiagonalNormal(0, 1))
end

@testset "DiagonalNormal API: Non-trivial cases" begin
    # Test log_pdf for D=1 case; x=0, μ=0, logσ=1
    DN = TMVDiagonalNormal([0.0], [1.])
    logP = logpdf(DN, [0.])
    @test logP ≈ -0.5*(2 + log(2*pi))

    # Test log_pdf D=3 case; x=0, μ=0, logσ=1
    DN = TMVDiagonalNormal(zeros(3), ones(3))
    X = zeros(1,3)
    logP = logpdf(DN, X[1,:])
    @test logP ≈ -0.5*(6 + 3*log(2*pi))

    # Test if case N=3,D=2 where we pass in x = column vec
    rng = MersenneTwister(9000)
    N = 3
    D = 2
    x = rand(rng, Float64, (N, D))
    DN = TMVDiagonalNormal(zeros(D), ones(D))
    logP = logpdf(DN, x[3,:])
    @test typeof(logP)==Float64
end

@testset "Getters" begin

    μ = [0.0]
    σ = [1.6]
    DN = TMVDiagonalNormal(μ, log.(σ))
    @test logσ(DN) == log.(σ)
    @test typeof(logσ(DN)) <: Vector{Float64}
    DN = TMVDiagonalNormal(TrackedArray(μ), TrackedArray(log.(σ)))
    @test logσ(DN) == TrackedArray(log.(σ))
    @test typeof(logσ(DN)) <: TrackedArray{T, N, Array{T,N}} where {T, N}

end


text:
#=
secondaryStructure:
- Julia version:
- Author: Dan
- Date: 2021-06-01
=#

mutable struct SecondaryStructure

end


text:
# Abstracts
export value
import Base: promote_rule, convert, show, eltype

"""
Physics Constant are composite types have method value
All Physics Constant can be promoted to corresponding
subtype of `Number`
"""
abstract type AbstractPhysicsConst{T} <: Number end

value(x::AbstractPhysicsConst) = x.value
eltype(x::AbstractPhysicsConst{T}) where T = T

promote_rule(::Type{F}, ::Type{T}) where {F, T<:AbstractPhysicsConst} = promote_type(eltype(F), eltype(T))
convert(::Type{F}, x::T) where {F<:Number, T<:AbstractPhysicsConst} = convert(F, value(x))

show(io::IO, x::AbstractPhysicsConst) = show(io, value(x))


text:
module QJuliaSolvers

  import QJuliaEnums
  import QJuliaBlas

  function identity_op(out::AbstractArray, inp::AbstractArray)
#    out .=@. inp
     QJuliaBlas.cpy(out, inp)
  end

  mutable struct QJuliaSolverParam_qj

    # Which linear solver to use
    inv_type::QJuliaEnums.QJuliaInverterType_qj
    # Which inner solver to use
    inv_type_precondition::QJuliaEnums.QJuliaInverterType_qj

    # Solver tolerance in the L2 residual norm
    tol::Float64
    # Solver tolerance in the L2 residual norm (used to restart InitCG)
    tol_restart::Float64

    # Whether to compute the true residual post solve
    compute_true_res::Bool
    # Actual L2 residual norm achieved in solver
    true_res::Float64
    # Maximum number of iterations in the linear solver
    maxiter::Int
    # Reliable update tolerance
    reliable_delta::Float64
    # Reliable update tolerance used in post multi-shift solver refinement
    reliable_delta_refinement::Float64
    # Whether to keep the partial solution accumuator in sloppy precision
    use_sloppy_partial_accumulator::Bool

    # Whether to use a pipelined solver with less global sums, gives pipeline length
    pipeline::Int
    # Krylov subspace dim
    nKrylov::Int
    # Preserve the source or not in the linear solver (deprecated)
    preserve_source::Bool
    # Whether to use initial guess
    use_init_guess::Bool

    # The precision used by the QJULIA solver
    dtype::DataType
    # The precision used by the QJULIA sloppy operator
    dtype_sloppy::DataType
    # The precision of the sloppy gauge field for the refinement step in multishift
    dtype_refinement_sloppy::DataType
    # The precision used by the QJULIA preconditioner
    dtype_precondition::DataType

    # Relaxation parameter used in MR,GCR etc. (default = 1.0)
    omega::Float64
    # Do we want global reductions?
    global_reduction::Bool

    # Number of steps
    Nsteps::Int

    # Shift
    shift::Float64

    # Reliable updates parameter
    delta::Float64

    #defualt constructor
    QJuliaSolverParam_qj() = new(QJuliaEnums.QJULIA_CG_INVERTER,
                                 QJuliaEnums.QJULIA_INVALID_INVERTER,
				 0.0,
				 0.0,
				 true,
				 0.0,
				 256,
				 1e-1,
				 1e-1,
				 false,
				 0,
				 0,
				 true,
                                 false,
				 Float64, Float64, Float64, Float64,
				 1.0,
				 true,
				 0,
				 0.0,
				 0.01)

  end #QJuliaSolverParam_qj

  include("../main/solvers/QJuliaMR.jl")
  include("../main/solvers/QJuliaLanMR.jl")
  include("../main/solvers/QJuliaGMRes.jl")
  include("../main/solvers/QJuliaPLCG.jl")
  include("../main/solvers/QJuliaP1CG.jl")
  include("../main/solvers/QJuliaPLCG2.jl")
  include("../main/solvers/QJuliaPCG.jl")
  include("../main/solvers/QJuliaPipePCG.jl")
  include("../main/solvers/QJuliaPipeFCG.jl")
  include("../main/solvers/QJuliaPipeNFCG.jl")
  include("../main/solvers/QJuliaPipeCG.jl")
  include("../main/solvers/QJuliaCGPCG.jl")

  function solve(out::AbstractArray, inp::AbstractArray, m::Any, mSloppy::Any, param::QJuliaSolverParam_qj, K::Function = identity_op)

    if((K != identity_op) && (param.inv_type_precondition == QJuliaEnums.QJULIA_INVALID_INVERTER))
      error("Preconditioner is not specified..")
    end

    if param.inv_type == QJuliaEnums.QJULIA_MR_INVERTER
      QJuliaMR.solver(out, inp, m,mSloppy, param)
    elseif param.inv_type == QJuliaEnums.QJULIA_GMRESDR_INVERTER
      QJuliaGMRes.solver(out, inp, m,mSloppy, param)
      #QJuliaPLCG.solver(out, inp, m,mSloppy, param)
      #QJuliaP1CG.solver(out, inp, m,mSloppy, param)
    elseif param.inv_type == QJuliaEnums.QJULIA_LANMR_INVERTER
      QJuliaLanMR.solver(out, inp, m, mSloppy, param)
    elseif param.inv_type == QJuliaEnums.QJULIA_PCG_INVERTER
      QJuliaPCG.solver(out, inp, m, mSloppy, param, K)
    elseif param.inv_type == QJuliaEnums.QJULIA_FCG_INVERTER
      if param.nKrylov == 0
        QJuliaPipeFCG.solver(out, inp, m, mSloppy, param, K)
      else
        QJuliaPipeNFCG.solver(out, inp, m, mSloppy, param, K)
      end
    elseif param.inv_type == QJuliaEnums.QJULIA_PIPEPCG_INVERTER
      #QJuliaCGPCG.solver(out, inp, m,mSloppy, param, K)
      if K != identity_op
        QJuliaPipePCG.solver(out, inp, m,mSloppy, param, K)
      else
        QJuliaPipeCG.solver(out, inp, m,mSloppy, param, K)
      end
    else
      error("Solver ", param.inv_type," is not available.")
    end

  end

end #QJuliaSolvers

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
function finalpositions(
    pointstofvecs::Dict{Point2D,Array{Float64,1}},
    scaler::Scaler,
    deltat::Float64,
    fd::Function,
    geps::Float64,
    deps::Float64,
    h0::Float64
)::Tuple{Any, Float64}
    newp, dpoints = buildnewpoints(pointstofvecs, scaler, deltat, geps, fd)
    finalp = buildfinalpoints(newp, scaler, fd, deps)
    return finalp, moveindex(dpoints, deltat, h0)
end

function buildnewpoints(
    pointstofvecs::Dict{Point2D,Array{Float64,1}},
    scaler::Scaler,
    deltat::Float64,
    geps::Float64,
    fd::Function
)::Tuple{Array{Array{Float64, 1},1},  Array{Array{Float64, 1},1}}
    newp = Array{Array{Float64, 1},1}()
    dpoints = Array{Array{Float64, 1},1}()
    for (point, force) in pointstofvecs
        p = unscaledpoint2d(point, scaler)
        np = [getx(p), gety(p)] + deltat * force
        push!(newp, [np[1], np[2]])
        if fd(np) < -geps
            push!(dpoints, force)
        end
    end
    return newp, dpoints
end

function buildfinalpoints(
    newp::Array{Array{Float64, 1},1},
    scaler::Scaler,
    fd::Function,
    deps::Float64
)
    finalp = Array{Array{Float64, 1}, 1}()
    for p in newp
        x = p[1]
        y = p[2]
        d = fd([x, y])
        if d > 0
            dgradx = (fd([x + deps, y]) - d) / deps
            dgrady = (fd([x, y + deps]) - d) / deps
            res = [x, y] - [d * dgradx, d * dgrady]
            push!(finalp, res)
        else
            push!(finalp, [x, y])
        end
    end
    finalp = map(p -> scaledpoint([p[1], p[2]], scaler), unique(finalp))
    finalp = transpose(reshape(vcat(finalp...), 2, length(finalp)))
end

function moveindex(
    dpoints::Array{Array{Float64,1},1},
    deltat::Float64,
    h0::Float64
)::Float64
    d = map(row -> sum(deltat * row .^ 2), dpoints)
    push!(d, 0)
    return maximum(sqrt.(d) / h0)
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "CSFML"
version = v"2.5.0"

# Collection of sources required to build CImGui
sources = [
    GitSource("https://github.com/SFML/CSFML.git",
              "61f17e3c1d109b65ef7e3e3ea1d06961da130afc")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/CSFML
mkdir build && cd build
CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release"
if [[ "${target}" == *mingw* ]]; then
    CMAKE_FLAGS="${CMAKE_FLAGS} -DCSFML_LINK_SFML_STATICALLY=false"
fi
cmake .. ${CMAKE_FLAGS}
make -j${nproc}
make install
install_license ../license.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("i686", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libcsfml-graphics", "csfml-graphics"], :libcsfml_graphics),
    LibraryProduct(["libcsfml-system", "csfml-system"], :libcsfml_system),
    LibraryProduct(["libcsfml-network", "csfml-network"], :libcsfml_network),
    LibraryProduct(["libcsfml-window", "csfml-window"], :libcsfml_window),
    LibraryProduct(["libcsfml-audio", "csfml-audio"], :libcsfml_audio),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("SFML_jll", v"2.5.1")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil