text:
# CLO15
# Ideals, Varieties, and Algorithms
# Cox, Little and O'Shea, Fourth edition, 2015

# Taken from CLO15
# They have been adapted to the grlex ordering

#@testset "Equality between algebraic sets" begin
#    @polyvar x y
#    @test (@set x^2 == y && y + x^2 == 4) == (@set x^2 == y && x^2 == 2)
#    @test (@set y == x^2 && x*z == y^2) == (@set y == x^2 && x*z == x^4)
#end

@testset "S-polynomial" begin
    Mod.@polyvar x y z
    @test spolynomial(x^3*y^2 - x^2*y^3 + x, 3x^4*y + y^2) == -x^3*y^3 + x^2 - y^3/3
    @test spolynomial(y - x^2, z - x^3) == z - x*y
    @test spolynomial(x^3 - 2x*y, x^2*y - 2y^2 + x) == -x^2
end

@testset "Groebner basis" begin
    Mod.@polyvar x y z
    function testg(G, H)
        presort!(G)
        _isz(f) = isapproxzero(rem(f, H))
        @test all(_isz.(G - H))
    end
    function testf(F, H)
        testg(gröbnerbasis(F), H)
    end
    testf([4x^2 + 5x, 3x^3], [x])
    testf([y - x^2, z - x^3], [x*z - y^2, x*y - z, x^2 - y, y^3 - z^2])
    testf([x^3 - 2x*y, x^2*y - 2y^2 + x], [y^2 - 0.5x, x*y, x^2])
    #p = 9x^11 + 27x^14 + x
    # root x=>-0.83005, y=>-3*(-0.83005)^4)
    F = [x^3*y^2 - x^2*y^3 + x, 3x^4*y + y^2]
    H = [x*y^3 - y^4 - 3x^3,
         x^3*y^2 - x^2*y^3 + x,
         x^4*y + y^2/3,
         x^5 + x*y/3,
         y^6 + 3y^5 + 9x^4 + 9x^3*y + y^3/3 - x^2 - x*y - y^2 - 3x]
    function testb(pre, sel)
        testg(groebnerbasis(F, Buchberger(Base.rtoldefault(Float64), pre, sel)), H)
    end
    testb(identity, dummyselection)
    testb(identity, normalselection)
    testb(presort!, dummyselection)
    testb(presort!, normalselection)
end

@testset "Reduce" begin
    Mod.@polyvar x y
    V1 = @set x == 1 && y == x^2
    @test rem(x^2 + 3x*y + 2y, ideal(V1)) == 6
    p = x^2 + x
    V2 = FullSpace()
    @test rem(p, ideal(V2)) === p
    # Needs MP v0.1.1
    #V3 = @set x == y^2
    #@test rem(x^2 + 3x*y + 2y + y^4, ideal(V3)) == rem(2y^4 + 3y^3 + 2y, ideal(V3))
end


text:
module PNGFiles
# Started as a fork of https://github.com/FugroRoames/LibPNG.jl

using ImageCore
using IndirectArrays
using OffsetArrays
using libpng_jll

libpng_wrap_dir = joinpath(@__DIR__, "..", "gen", "libpng")
using CEnum
include(joinpath(libpng_wrap_dir, "ctypes.jl"))
include(joinpath(libpng_wrap_dir, "libpng_common.jl"))
include(joinpath(libpng_wrap_dir, "libpng_api.jl"))

const readcallback_c = Ref{Ptr{Cvoid}}(C_NULL)
const writecallback_c = Ref{Ptr{Cvoid}}(C_NULL)
const png_error_fn_c = Ref{Ptr{Cvoid}}(C_NULL)
const png_warn_fn_c = Ref{Ptr{Cvoid}}(C_NULL)

include("wraphelpers.jl")
include("utils.jl")
include("io.jl")

function __init__()
    readcallback_c[] = @cfunction(_readcallback, Cvoid, (png_structp, png_bytep, png_size_t));
    writecallback_c[] = @cfunction(_writecallback, Csize_t, (png_structp, png_bytep, png_size_t));
    png_error_fn_c[] = @cfunction(png_error_handler, Cvoid, (Ptr{Cvoid}, Cstring))
    png_warn_fn_c[] = @cfunction(png_warn_handler, Cvoid, (Ptr{Cvoid}, Cstring))
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(load, (String,))
    # save is harder to usefully precompile because of the diversity of array types
end
VERSION >= v"1.4.2" && _precompile_() # https://github.com/JuliaLang/julia/pull/35378

end # module

ext: jl


text:
@kwargs function reduce!(
    op::MonoidUnion,
    w::GBVector,
    A::GBMatOrTrans)
    A, desc, _ = _handletranspose(A, desc, nothing)
    op = getoperator(op, eltype(w))
    libgb.GrB_Matrix_reduce_Monoid(w, mask, accum, op, A, desc)
end

@kwargs function Base.reduce(op::MonoidUnion, A::GBMatOrTrans;
    dims = 2, typeout = nothing, init = nothing
)
    if typeout === nothing
        typeout = eltype(A)
    end
    if dims == 2
        w = GBVector{typeout}(size(A, 1))
        reduce!(op, w, A; desc, accum, mask)
        return w
    elseif dims == 1
        desc = desc + Descriptors.T0
        w = GBVector{typeout}(size(A, 2))
        reduce!(op, w, A; desc, accum, mask)
        return w
    elseif dims == (1,2)
        if init === nothing
            c = Ref{typeout}()
            typec = typeout
        else
            c = Ref(init)
            typec = typeof(init)
        end
        op = getoperator(op, typec)
        A, desc, _ = _handletranspose(A, desc, nothing)
        libgb.scalarmatreduce[typeout](c, accum, op, A, desc)
        return c[]
    end
end

@kwargs function Base.reduce(
    op::MonoidUnion, v::GBVector;
    typeout = nothing, init = nothing
)
    if typeout === nothing
        typeout = eltype(v)
    end
    if init === nothing
        c = Ref{typeout}()
        typec = typeout
    else
        c = Ref(init)
        typec = typeof(init)
    end
    op = getoperator(op, typec)
    libgb.scalarvecreduce[typeout](c, accum, op, v, desc)
    return c[]
end


text:
# ---
# title: 1536. Minimum Swaps to Arrange a Binary Grid
# id: problem1536
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/minimum-swaps-to-arrange-a-binary-grid/description/>
# hidden: true
# ---
#
# Given an `n x n` binary `grid`, in one step you can choose two **adjacent
# rows** of the grid and swap them.
#
# A grid is said to be **valid** if all the cells above the main diagonal are
# **zeros**.
#
# Return _the minimum number of steps_ needed to make the grid valid, or **-1**
# if the grid cannot be valid.
#
# The main diagonal of a grid is the diagonal that starts at cell `(1, 1)` and
# ends at cell `(n, n)`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/28/fw.jpg)
#
#
#
#     Input: grid = [[0,0,1],[1,1,0],[1,0,0]]
#     Output: 3
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/16/e2.jpg)
#
#
#
#     Input: grid = [[0,1,1,0],[0,1,1,0],[0,1,1,0],[0,1,1,0]]
#     Output: -1
#     Explanation: All rows are similar, swaps have no effect on the grid.
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/16/e3.jpg)
#
#
#
#     Input: grid = [[1,0,0],[1,1,0],[1,1,1]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `n == grid.length`
#   * `n == grid[i].length`
#   * `1 <= n <= 200`
#   * `grid[i][j]` is `0` or `1`
#
#
## @lc code=start
using LeetCode

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


text:
using SpecialFunctions, ForwardDiff, FiniteDifferences
include("shared.jl")

const BIG_FD    = central_fdm(10,1)
const BIG_FD_O2 = central_fdm(10,2)

besselkdv(v,x)  = ForwardDiff.derivative(_v->BesselK._besselk(_v,x,100,1e-12,false), v)
besselkdv2(v,x) = ForwardDiff.derivative(_v->besselkdv(_v,x), v)

dbesselkdv(v, x)  = BIG_FD(_v->besselk(_v,x), v)
dbesselkdv2(v, x) = BIG_FD_O2(_v->besselk(_v,x), v)

fastdbesselkdv2(v, x) = (besselk(v+2e-6, x) - 2*besselk(v+1e-6, x) + besselk(v, x))/1e-12

const BASELINE = [dbesselkdv2(z[1], z[2])  for z in Iterators.product(VGRID, XGRID)]
const OURSOL   = [besselkdv2(z[1], z[2]) for z in Iterators.product(VGRID, XGRID)]
const FASTFD   = [fastdbesselkdv2(z[1], z[2])  for z in Iterators.product(VGRID, XGRID)]
const TOLS     = atolfun.(zip(BASELINE, OURSOL))
const FDTOLS   = atolfun.(zip(BASELINE, FASTFD))
const DIFTOLS  = log10.(TOLS) .- log10.(FDTOLS)

let res = findmax(DIFTOLS)
  println("Worst case derivative difference: $(res[1])")
  println("Value of dfun: $(BASELINE[res[2]])")
end

gnuplot_save_matrix!("../plotdata/atols_deriv2_fd.csv",    FDTOLS, VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/atols_deriv2_ad.csv",      TOLS, VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/atols_deriv2_fdad.csv", DIFTOLS, VGRID, XGRID)


text:
"""
    LRhoSparseOpProblem <: AbstractProblem

Problem or finding the steady state of a ℒdagℒ matrix by computing
𝒞 = ∑|ρ(σ)|²|⟨⟨σ|ℒ |ρ⟩⟩|² only storing H and c_ops.

DO NOT USE WITH COMPLEX-WEIGHT NETWORKS, AS IT DOES NOT WORK
"""
struct LRhoSparseOpProblem{B, SM} <: LRhoSquaredProblem where {B<:Basis,
                                                 SM<:SparseMatrixCSC}
    HilbSpace::B            # 0
    HnH::SM
    HnH_t::SM
    L_ops::Vector{SM}       # 4
    L_ops_h::Vector{SM}     # 4
    L_ops_t::Vector{SM}     # 5
    ρss
end

"""
    LRhoSparseOpProblem([T=STD_REAL_PREC], args...)

Creates a problem for minimizing the cost function 𝒞 = ∑|ρ(σ)|²|⟨⟨σ|ℒ |ρ⟩⟩|².
Computes |⟨⟨σ|ℒ |ρ⟩⟩| by computing on the fly the commutator with the
Hamiltonian and with the collapse operators.

`args...` can either be a `GraphLindbladian`, or the Hamiltonian and a vector
of collapse operators.

`T=STD_REAL_PREC` by default is the numerical precision used. It should match that of
the network.
"""
LRhoSparseOpProblem(args...) = LRhoSparseOpProblem(STD_REAL_PREC, args...)
LRhoSparseOpProblem(T::Type{<:Number}, gl::GraphLindbladian) =
    LRhoSparseOpProblem(T, basis(gl), SparseOperator(hamiltonian(gl)), jump_operators(gl))
LRhoSparseOpProblem(T::Type{<:Number}, Ham::DataOperator, cops::Vector) =
    LRhoSparseOpProblem(T, Ham.basis_l, Ham, cops)
function LRhoSparseOpProblem(T::Type{<:Number}, Hilb::Basis, Ham::DataOperator, c_ops_q::Vector)
    # Fix complex numbers
    if real(T) == T
        T = Complex{T}
    end

    # Generate H_eff
    H_eff   = deepcopy(T.(Ham.data))
    ST = typeof(H_eff)

    c_ops       = Vector{ST}(undef, length(c_ops_q))
    c_ops_h     = Vector{ST}(undef, length(c_ops_q))
    c_ops_trans = Vector{ST}(undef, length(c_ops_q))
    for i=1:length(c_ops)
        c_ops[i]       = c_ops_q[i].data
        c_ops_h[i]     = c_ops[i]'
        c_ops_trans[i] = transpose(c_ops[i])
        H_eff         -= 0.5im * (c_ops[i]'*c_ops[i])
    end

    LRhoSparseOpProblem{typeof(Hilb), ST}(Hilb,                  # 0
                    H_eff,
                    transpose(H_eff),
                    c_ops,
                    c_ops_h,
                    c_ops_trans,
                    0.0)
end

basis(prob::LRhoSparseOpProblem) = prob.HilbSpace

function compute_Cloc!(LLO_i, ∇lnψ, prob::LRhoSparseOpProblem, net::MatrixNet, 𝝝,
                      lnψ=net(𝝝), 𝝝p=deepcopy(𝝝))
    HnH = prob.HnH
    HnH_t = prob.HnH_t
    c_ops = prob.L_ops
    c_ops_h = prob.L_ops_h
    c_ops_trans = prob.L_ops_t

    σ  = row(𝝝)
    σt = col(𝝝)
    set_index!(𝝝p, index(𝝝))
    𝝝p_row = row(𝝝p)
    𝝝p_col = col(𝝝p)

    for el=LLO_i
      el .= 0.0
    end

    i_σt = index(σt)
    i_σ  = index(σ)

    C_loc = zero(Complex{real(out_type(net))})

    # ⟨σ|Hρ|σt⟩ (using hermitianity of HdH)
    set!(𝝝p_col, toint(σt))
    for row_id = HnH_t.colptr[i_σ]:(HnH_t.colptr[i_σ+1]-1)
      i_σ_p = HnH_t.rowval[row_id]
      set_index!(𝝝p_row, i_σ_p)
      lnψ_i, ∇lnψ_i = logψ_and_∇logψ!(∇lnψ, net, 𝝝p)
      C_loc_i  =  -1.0im * HnH_t.nzval[row_id] * exp(lnψ_i - lnψ)

      for (LLOave, _∇lnψ)= zip(LLO_i, ∇lnψ_i.tuple_all_weights)
        LLOave .+= C_loc_i .* _∇lnψ
      end
      C_loc  += C_loc_i
     end

    # ⟨σ|ρHᴴ|σt⟩
    set!(𝝝p_row, toint(σ))
    for row_id = HnH.colptr[i_σt]:(HnH.colptr[i_σt+1]-1)
      i_σ_p = HnH.rowval[row_id]
      set_index!(𝝝p_col, i_σ_p)
      lnψ_i, ∇lnψ_i = logψ_and_∇logψ!(∇lnψ, net, 𝝝p)
      C_loc_i  =  1.0im * conj(HnH_t.nzval[row_id]) * exp(lnψ_i - lnψ)

      for (LLOave, _∇lnψ)= zip(LLO_i, ∇lnψ_i.tuple_all_weights)
        LLOave .+= C_loc_i .* _∇lnψ
      end
      C_loc  += C_loc_i
    end

    # L rho Ldag H #ok
    # -im ⟨σ|L ρ Lᴴ|σt⟩
    for i=1:length(c_ops_h)
      Ld  = c_ops_h[i]
      L   = c_ops_trans[i]

      for ext_row_id = L.colptr[i_σ]:(L.colptr[i_σ+1]-1)
        i_σ_p = L.rowval[ext_row_id]
        set_index!(𝝝p_row, i_σ_p)
        val_σ_p = L.nzval[ext_row_id]

        for int_row_id = Ld.colptr[i_σt]:(Ld.colptr[i_σt+1]-1)
          i_σt_p = Ld.rowval[int_row_id]
          set_index!(𝝝p_col, i_σt_p)

          lnψ_i, ∇lnψ_i = logψ_and_∇logψ!(∇lnψ, net, 𝝝p)
          C_loc_i  =  val_σ_p * Ld.nzval[int_row_id] *  exp(lnψ_i - lnψ)

          for (LLOave, _∇lnψ)= zip(LLO_i, ∇lnψ_i.tuple_all_weights)
            LLOave .+= C_loc_i .* _∇lnψ
          end
          C_loc  += C_loc_i
        end
      end
    end

    return C_loc
end

# pretty printing
Base.show(io::IO, p::LRhoSparseOpProblem) = print(io,
    "LRhoSparseOpProblem on space : $(basis(p)) computing the variance of Lrho using sparse H, c_ops")


text:
Base.MainInclude.include("ToPTformat/toRegName.jl")
rootdir = "D:/99-RAWData/xxxx-2019-0224~03.02/data"

Main.regFiles( rootdir  )

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
@recipe function f(
    sol::ODESolution,
    H::OpenQuantumBase.AbstractHamiltonian,
    lvl::Integer,
    t_axis=sol.t,
)
    y = instantaneous_population(sol, H, lvl=lvl, t_axis=t_axis)
    lab = ["\$E_$x\$" for x in (0:lvl - 1)']
    label --> lab
    (t_axis, y)
end

@recipe function f(
    sol::ODESolution,
    H::OpenQuantumBase.AbstractHamiltonian,
    lvl::AbstractArray{T,1},
    t_axis=sol.t,
) where {T <: Integer}
    y = instantaneous_population(sol, H, lvl=maximum(lvl) + 1, t_axis=t_axis)
    lab = ["\$E_$x\$" for x in (lvl)']
    label --> lab
    (t_axis, y[:, lvl .+ 1])
end

@recipe function f(
    sol::ODESolution,
    lvl::AbstractArray{T,1},
    s_axis=sol.t,
) where {T <: Integer}
    y = []
    if ndims(sol.u[1]) == 1
        for x in s_axis
            push!(y, abs2.(sol(x)[lvl .+ 1]))
        end
    elseif ndims(sol.u[1]) == 2
        for x in s_axis
            push!(y, real.(diag(sol(x))[lvl .+ 1]))
        end
    else
        throw(ArgumentError("Solution needs to be state vector or density matrix."))
    end

    y = hcat(y...)'
    (s_axis, y[:, lvl .+ 1])
end


text:
# Compares with Pal's computation:

if VERSION >= v"0.7"
  using DelimitedFiles
end
using PyPlot
include("../../../src/transit_poly_struct.jl")

data_pal = readdlm("ntiq_test.out")

# Now, run the same with transit_poly:

u = [0.2,0.3]
nb = 1001; b=zeros(nb)
r = 0.1  # radius ratio is 0.1
flux = zeros(nb); fgrad = zeros(4)
flux_big = zeros(nb); fgrad_big = zeros(BigFloat,4)
flux_grad = zeros(nb,4)
flux_grad_big = zeros(nb,4)
t1 = time_ns()
for i=0:1000
  b[i+1] = sqrt(((i-500.0)/500.0*(1.0+2.0*r))^2)
  flux[i+1]= transit_poly!(r,b[i+1],u,fgrad)
  flux_grad[i+1,:] = fgrad
end
println("elapsed time: " * string(1e-9 * (time_ns() - t1)) * " seconds")
for i=0:1000
  flux_big[i+1]= convert(Float64,transit_poly!(big(r),big(b[i+1]),big.(u),fgrad_big))
  flux_grad_big[i+1,:] = convert(Array{Float64,1},fgrad_big)
end

# Now, plot results:

fig,axes = subplots(3,1, figsize=(8,12))

ax = axes[1]
ax[:plot](b,flux,label="ALFM (2019)",linestyle="-",lw=2,color="C0")
ax[:plot](b,1 .- data_pal[:,1],label="Pal (2008)",linestyle="--",lw=1,color="C1")
ax[:set_ylabel]("Flux", fontsize=14)
ax[:legend](loc="upper left", fontsize=8)
ax[:set_title](L"$r = 0.1$", fontsize=14)

ax = axes[2]
ax[:plot](b,flux_grad[:,1],linestyle="-",lw=2,label="ALFM (2019)", color="C0")
ax[:plot](b,flux_grad[:,2],linestyle="-",lw=2, color="C0")
ax[:plot](b,50.0*flux_grad[:,3],linestyle="-",lw=2, color="C0")
ax[:plot](b,50.0*flux_grad[:,4],linestyle="-",lw=2, color="C0")
ax[:plot](b,-data_pal[:,2],ls="--",lw=1, label="Pal (2008)", color="C1")
ax[:plot](b,-data_pal[:,3],ls="--",lw=1, color="C1")
ax[:plot](b,-50.0*data_pal[:,4],ls="--",lw=1, color="C1")
ax[:plot](b,-50.0*data_pal[:,5],ls="--",lw=1, color="C1")
ax[:legend](loc="upper left",fontsize=8)
ax[:set_ylabel]("Derivatives", fontsize=14)

ax[:annotate](L"$\partial f / \partial r$", xy=(0.7, -0.185), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=11, fontsize=10)

ax[:annotate](L"$\partial f / \partial b$", xy=(0.1, 0.025), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=0)

ax[:annotate](L"$50\times\partial f / \partial u_1$", xy=(0.1, -0.18), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=4)

ax[:annotate](L"$50\times\partial f / \partial u_2$", xy=(0.1, -0.085), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=0)

ax = axes[3]
ax[:semilogy](b,abs.(1 .- data_pal[:,1]-flux_big),label=L"$f$")
ax[:plot](b,abs.(-data_pal[:,2]-flux_grad_big[:,1]),label=L"$\frac{\partial f}{\partial r}$")
ax[:plot](b,abs.(-data_pal[:,3]-flux_grad_big[:,2]),label=L"$\frac{\partial f}{\partial b}$")
ax[:plot](b,abs.(-data_pal[:,4]-flux_grad_big[:,3]),label=L"$\frac{\partial f}{\partial u_1}$")
ax[:plot](b,abs.(-data_pal[:,5]-flux_grad_big[:,4]),label=L"$\frac{\partial f}{\partial u_2}$")
ax[:set_xlabel](L"$b$", fontsize=18)
ax[:set_ylabel]("Difference", fontsize=14)
ax[:legend](loc="upper right",fontsize=10)

savefig("compare_pal.pdf", bbox_inches="tight")

max_stars_repo_name: dfm/limbdark

max_issues_repo_name: dfm/limbdark

max_forks_repo_name: dfm/limbdark


text:
using Plots,FFTW
using JLD2

for (ind, arg) in enumerate(ARGS)
    t = typeof(arg)
    msg = "$(ind) -> $arg :$t"
    println(msg)
end

function sqaure_windowing_signal(signal::Array{Float64})
    windowed_signal = signal
    window_correction_factor = 1
    return (windowed_signal,window_correction_factor)
end

function hanning_windowing_signal(signal::Array{Float64})
    windowed_signal::Array{Float64} = []
    window::Array{Float64} = []
    total_window_power = 0
    j = 1
    for w = signal
        window_val = 0.5*(1 - cos(2*pi*j/length(signal)))
        push!(window,window_val)
        windowed = window_val*w
        push!(windowed_signal, windowed)
        total_window_power += window_val
        j += 1
    end
    window_correction_factor = length(signal)/total_window_power
    return (windowed_signal,window_correction_factor)
end

function main()
    sampling_rate = 0.0
    window_size = 64*1024*1024;
    w_m_vals::Array{Float64} = []

    jldopen(ARGS[1], "r") do f
        sampling_rate = f["samplingrate"]
        w_m_vals = f["w_m"]
    end

    println(sampling_rate)
    println(length(w_m_vals))

    #(windowed_signal, window_correction_factor)  = sqaure_windowing_signal(w_m_vals[1:window_size])
    (windowed_signal, window_correction_factor) = hanning_windowing_signal(w_m_vals[1:window_size])
    println(window_correction_factor)
    powers = map(x->20*log10(abs(x)*window_correction_factor),fft(windowed_signal)|>fftshift)
    freqs = fftfreq( window_size, sampling_rate )|>fftshift

    plot_graph=plot(freqs, powers, xlims=(0,2000*1000), fmt =:png)
    savefig(plot_graph, "fft1.png")
    plot_graph=plot(freqs, powers, xlims=(944*1000,964*1000), fmt =:png)
    savefig(plot_graph, "fft2.png")
end

main()

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
#=
timeseries decompositions, i.e. separating some data into components
=#
export seasonal_decomposition
"""
    seasonal_decomposition(A::ClimArray, fs = [1, 2]) → seasonal, residual
Decompose `A` into a seasonal and residual components, where the seasonal contains the
periodic parts of `A`, with frequencies given in `fs`, and residual contains what's left.

`fs` is measured in 1/year. This function works even for non-equispaced time axis (e.g.
monthly averages) and uses LPVSpectral.jl and SignalDecomposition.jl.
"""
function seasonal_decomposition(A::AbDimArray, fs::Vector)
    @assert hasdim(A, Time)
    E = _numbertype(eltype(A))
    method = Sinusoidal(E.(fs ./ DAYS_IN_ORBIT))
    seasonal = DimensionalData.basetypeof(A)(copy(Array(A)), dims(A))
    residual = DimensionalData.basetypeof(A)(copy(Array(A)), dims(A))

    t = dims(A, Time).val
    truetime = realtime_days(t, E)
    for i in otheridxs(A, Time)
        y = Array(A[i...])
        sea, res = SignalDecomposition.decompose(truetime, y, method)
        seasonal[i...] .= sea
        residual[i...] .= res
    end
    return seasonal, residual
end

function seasonal_decomposition(A::AbDimArray{T, 1}, fs::Vector) where {T}
    @assert hasdim(A, Time)
    E = _numbertype(eltype(A))
    method = Sinusoidal(E.(fs ./ DAYS_IN_ORBIT))
    seasonal = DimensionalData.basetypeof(A)(copy(Array(A)), dims(A))
    residual = DimensionalData.basetypeof(A)(copy(Array(A)), dims(A))

    t = dims(A, Time).val
    truetime = realtime_days(t, E)
    y = Array(A)
    sea, res = SignalDecomposition.decompose(truetime, y, method)
    seasonal .= sea
    residual .= res
    return seasonal, residual
end


text:
using FactCheck
using LLRBTrees
#using LLRBVisualize

include("checkheight.jl")

function checktree{K,V}(tree::LLRBTree{K,V}, elements::Array{TreeNode{K,V},1})

    FactCheck.@fact tree --> check_height "not balanced"
    #order
    #no consecutive red nodes
    #left leaning

end
function randomscenarios(seed::Int=1, scenarios::Int=10, elNum=1000, var=10)
    #Fixed seed for reproducibility if seed=0 then random
    if(seed==0)
        seed=rand()
        srand(seed)
    else
        srand(seed)
    end

    FactCheck.facts("Test differents scenarios") do
        for i in 1:scenarios
            #Define a number of nodes that varies with scenario, and a random keyrange
            #that's at least as big as the number of nodes
            nodes = rand(elNum-var:elNum)
            mult = rand(1:var)
            keysmax = rand(mult*elNum:((mult+1)*elNum)+var)

            tree = LLRBTree{Int, ASCIIString}()
            keys = []
            elements = TreeNode{Int, ASCIIString}[]

            FactCheck.context("Add random values") do
                for i in 1:nodes
                    key = rand(mult*elNum:keysmax)
                    node = TreeNode(key, string(key) )
                    #println(node)
                    push!(tree, node)

                    indexes = Int[]
                    #Remove repeated notes
                    for j in 1:size(elements)[1]
                        nd = elements[j]
                        if nd.key ==  node.key.value
                            push!(indexes, j)
                        end
                    end
                    deleteat!(elements, indexes) #Remove repeated nodes
                    push!(elements, node)
                    push!(keys, node.key.value)

                    checktree(tree, elements)
                end
            end
            #println(keys)

            FactCheck.context("Search, min and max") do
                min = minimum(keys)
                max = maximum(keys)
                #FactCheck.@fact_throws KeyException tree[min] "getindex failed"
                FactCheck.@fact minimum(tree).key.value == min --> true "Minimum mismatch"
                FactCheck.@fact maximum(tree).key.value == max --> true "Maximum mismatch"

            end

            #vis=drawtree(tree)
            #img = SVG("test1.svg", 10cm, 10cm)
            #draw(img, vis)

            FactCheck.context("Remove keys not present") do
                lower = elNum + 1 + var
                upper = elNum + 1 + 2*var
                removesMax = rand(lower:upper)
                for i in lower:removesMax
                    #println("delete: ",i, " ", )
                    delete!(tree, i)
                    checktree(tree, elements)
                end
            end
            #println(tree)
            #vis=drawtree(tree)
            #img = SVG("test2.svg", 10cm, 10cm)
            #draw(img, vis)

            #Remove all keys in random order
            FactCheck.context("Remove all keys") do
                i=1 #for debugging
                while( size(elements)[1] > 0 )
                    checktree(tree, elements)

                    index = rand(1:size(elements)[1])

                    #Debugging
                    #println("removed ",elements[index])
                    #vis=drawtree(tree)
                    #img = SVG("test"*string(i)*".svg", 10cm, 10cm)
                    #draw(img, vis)
                    #i+=1

                    FactCheck.@fact isa(delete!( tree, elements[index].key.value ), LLRBTree{Int, ASCIIString}) --> true  "Did not return the tree with deleted node"
                    deleteat!(elements, index)
                    deleteat!(keys, index)
                end

                FactCheck.@fact isleaf(tree.root) -->  true "root was not leaf after deleting all keys"
            end

        end
    end
end


text:
# ---
# title: 1169. Invalid Transactions
# id: problem1169
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, String
# link: <https://leetcode.com/problems/invalid-transactions/description/>
# hidden: true
# ---
#
# A transaction is _possibly invalid_ if:
#
#   * the amount exceeds /$/1000, or;
#   * if it occurs within (and including) 60 minutes of another transaction with the same name in a different city.
#
# Each transaction string `transactions[i]` consists of comma separated values
# representing the name, time (in minutes), amount, and city of the transaction.
#
# Given a list of `transactions`, return a list of transactions that are
# possibly invalid.  You may return the answer in any order.
#
#
#
# **Example 1:**
#
#
#
#     Input: transactions = ["alice,20,800,mtv","alice,50,100,beijing"]
#     Output: ["alice,20,800,mtv","alice,50,100,beijing"]
#     Explanation: The first transaction is invalid because the second transaction occurs within a difference of 60 minutes, have the same name and is in a different city. Similarly the second one is invalid too.
#
# **Example 2:**
#
#
#
#     Input: transactions = ["alice,20,800,mtv","alice,50,1200,mtv"]
#     Output: ["alice,50,1200,mtv"]
#
#
# **Example 3:**
#
#
#
#     Input: transactions = ["alice,20,800,mtv","bob,50,1200,mtv"]
#     Output: ["bob,50,1200,mtv"]
#
#
#
#
# **Constraints:**
#
#   * `transactions.length <= 1000`
#   * Each `transactions[i]` takes the form `"{name},{time},{amount},{city}"`
#   * Each `{name}` and `{city}` consist of lowercase English letters, and have lengths between `1` and `10`.
#   * Each `{time}` consist of digits, and represent an integer between `0` and `1000`.
#   * Each `{amount}` consist of digits, and represent an integer between `0` and `2000`.
#
#
## @lc code=start
using LeetCode

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


text:
__precompile__()

module DispatchError

immutable Foo{A,B}
    a::A
    b::B
    c::Int
    # disable default outer constructor
    Foo(a, b, c) = new(a, b, c)
end

Foo(args...) = error("boo!")

Foo{A<:AbstractArray,B<:AbstractArray}(a::A, b::B, c::Int) = info("yay!")

end # module


text:
module Voronoi

include("Geometry.jl")
include("Diagram.jl")
include("EventQueue.jl")
include("BeachLine.jl")
include("Fortune.jl")
include("Intersect.jl")
include("Optimization.jl")

end

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
@testset "get_n_words!" begin

    line = "#=GF AC PF00571"
    @test get_n_words(line, 1) == String[line]
    @test get_n_words(line, 2) == String["#=GF", "AC PF00571"]
    @test get_n_words(line, 3) == String["#=GF", "AC", "PF00571"]
    @test get_n_words(line, 4) == String["#=GF", "AC", "PF00571"]

    @test get_n_words("\n",1) == String["\n"]
    @test get_n_words("#", 1) == String["#"]

    # ASCII
    str = "#=GR O31698/18-71 SS    CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH"
    @test get_n_words(str, 3) ==
        String["#=GR", "O31698/18-71", "SS    CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH"]

    # UTF-8
    str = "#=GF CC   (Römling U.  and Galperin M.Y. “Bacterial cellulose"
    @test get_n_words(str, 3) ==
        String["#=GF", "CC", "(Römling U.  and Galperin M.Y. “Bacterial cellulose"]

    str = "#=GF CC   not present in all SecA2–SecY2 systems. This family of Asp5 is"
    @test get_n_words(str, 3) ==
        String["#=GF", "CC", "not present in all SecA2–SecY2 systems. This family of Asp5 is"]
end

@testset "hascoordinates" begin

    @test hascoordinates("O83071/192-246")
    @test !hascoordinates("O83071")
end

@testset "select_element" begin

    @test select_element([1]) == 1
    @test_throws ErrorException select_element([])
end

@testset "Matrices to and from lists" begin

    @testset "matrix2list" begin

        mat = [ 1 2 3
                4 5 6
                7 8 9 ]

        @test matrix2list(mat) == [2, 3, 6]
        @test matrix2list(mat, diagonal=true) == [1, 2, 3, 5, 6, 9]
        @test matrix2list(mat, part="lower") == [4, 7, 8]
        @test matrix2list(mat, part="lower", diagonal=true) == [1, 4, 7, 5, 8, 9]
    end

    @testset "list2matrix" begin

        mat = [ 1 2 3
                2 5 6
                3 6 9 ]

        @test triu(list2matrix([2, 3, 6], 3), 1) == triu(mat, 1)
        @test list2matrix([1, 2, 3, 5, 6, 9], 3, diagonal=true) == mat
    end
end

@testset "General IO" begin

    @testset "lineiterator" begin

        # Julia 0.6: eachline return lines without line endings by default
        ppap = "pen\npineapple\napple\npen\n"
        @test collect(lineiterator(ppap)) == collect(eachline(IOBuffer(ppap)))

        @test collect(lineiterator("Hola")) == ["Hola"]
        @test collect(lineiterator("Hola\n")) == ["Hola"]
        @test collect(lineiterator("\n")) == [""]
        @test collect(lineiterator("Hola\nMundo")) == ["Hola", "Mundo"]
        @test collect(lineiterator("Hola\nMundo\n")) == ["Hola", "Mundo"]
        @test collect(lineiterator("Hola\nMundo\n\n")) == ["Hola", "Mundo", ""]
    end

    @testset "File checking" begin

        file_path = joinpath(DATA, "simple.fasta")

        @test isnotemptyfile(file_path)
        @test !isnotemptyfile(joinpath(DATA, "emptyfile"))

        @test check_file(file_path) == file_path
        @test_throws ErrorException check_file("nonexistentfile")
    end

    @testset "Download file" begin

        try
            @test ".tmp" == download_file("http://www.uniprot.org/uniprot/P69905.fasta",
                ".tmp",
		headers = Dict("User-Agent" => "Mozilla/5.0 (compatible; MSIE 7.01; Windows NT 5.0)"),
                redirect =true)
        finally
            if isfile(".tmp")
                rm(".tmp")
            end
        end
    end
end


text:
using PlotlyJS

x = 1:0.01:10
y = sin.(x)

p = plot(x, y);

savefig(p, "plotly-test.html")


text:
"""
    to_matrix(points)

Convert collection of d-dimensional points to d×n matrix
"""
function to_matrix(points)
    dim = length(points[1])
    T = eltype(points[1])
    n = length(points)
    result = zeros(T, (dim, n))

    for (i, p) in enumerate(points)
        length(p) == dim || throw(ArgumentError("points must have the same length"))
        result[:, i] .= p
    end

    return result
end

"""
    distances(points, metric=Euclidean(1e-12))

Return distance matrix calculated from `points` with `metric`.
"""
function distances(points, metric=Euclidean(1e-12))
    points_mat = to_matrix(points)
    dists = pairwise(metric, points_mat; dims=2)
    return dists
end

"""
    radius(dists)
    radius(points[, metric=Euclidean(1e-12)])

Calculate the radius of the space. This is used for default `thresholds`.
"""
function radius(dists::AbstractMatrix)
    return minimum(maximum(abs, dists[:, i]) for i in 1:size(dists, 1))
end
function radius(dists::SparseMatrixCSC)
    return maximum(dists)
end
function radius(points, metric=Euclidean(1e-12))
    radius = Inf
    for p in points
        p_max = 0.0
        for q in points
            p == q && continue
            p_max = max(p_max, metric(SVector(p), SVector(q)))
        end
        radius = min(p_max, radius)
    end
    return radius
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule spectra_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("spectra")
JLLWrappers.@generate_main_file("spectra", UUID("ef34ea38-0a9e-523f-b228-55b6d5c13eb1"))
end  # module spectra_jll


text:
module LineNumbers

export SourceFile, compute_line, LineBreaking
import Base: getindex, setindex!, length

# Offsets are 0 based
struct SourceFile
    data::Vector{UInt8}
    offsets::Vector{UInt64}
end
length(file::SourceFile) = length(file.offsets)

function SourceFile(data)
    offsets = UInt64[0]
    buf = IOBuffer(data)
    local line
    while !eof(buf)
        line = readuntil(buf,'\n')
        !eof(buf) && push!(offsets, position(buf))
    end
    if !isempty(offsets) && line[end] == '\n'
        push!(offsets, position(buf))
    end
    SourceFile(Vector{UInt8}(data),offsets)
end

function compute_line(file::SourceFile, offset)
    ind = searchsortedfirst(file.offsets, offset)
    ind <= length(file.offsets) && file.offsets[ind] == offset ? ind : ind - 1
end

function getindex(file::SourceFile, line::Int)
    if line == length(file.offsets)
        return file.data[(file.offsets[end]+1):end]
    else
        # - 1 to skip the '\n'
        return file.data[(file.offsets[line]+1):(file.offsets[line+1]-1)]
    end
end
getindex(file::SourceFile, arr::AbstractArray) = [file[x] for x in arr]

# LineBreaking

"""
Indexing adaptor to map from a flat byte offset to a `[line][offset]` pair.
Optionally, off may specify a byte offset relative to which the line number and
offset should be computed
"""
struct LineBreaking{T}
    off::UInt64
    file::SourceFile
    obj::T
end

function indtransform(lb::LineBreaking, x::Int)
    offline = compute_line(lb.file, lb.off)
    line = compute_line(lb.file, x)
    lineoffset = lb.file.offsets[line]
    off = x - lineoffset + 1
    if lineoffset < lb.off
        off -= lb.off - lineoffset
    end
    (line - offline + 1), off
end

function getindex(lb::LineBreaking, x::Int)
    l, o = indtransform(lb, x)
    lb.obj[l][o]
end

function setindex!(lb::LineBreaking, y, x::Int)
    l, o = indtransform(lb, x)
    lb.obj[l][o] = y
end
function setindex!(lb::LineBreaking, y, x::AbstractArray)
    for i in x
        lb[i] = y
    end
end

end # module


text:
DiffEqBase.has_reinit(i::AbstractNoiseProcess) = true
function DiffEqBase.reinit!(W::Union{NoiseProcess,NoiseApproximation},dt;
                            t0 = W.t[1],
                            erase_sol = true,
                            setup_next = false)

  if erase_sol
    resize!(W.t,1)
    resize!(W.W,1)
    if W.Z != nothing
      resize!(W.Z,1)
    end
  end

  W.curt = t0
  W.dt = dt
  if typeof(W) <: NoiseApproximation
    reinit!(W.source1)
    if W.source2 != nothing
      reinit!(W.source2)
    end
  end

  if isinplace(W)
    W.curW .= first(W.W)
    if W.Z != nothing
      W.curZ .= first(W.Z)
    end
  else
    W.curW = first(W.W)
    if W.Z != nothing
      W.curZ = first(W.Z)
    end
  end

  if typeof(W) <: NoiseProcess
    while !isempty(W.S₁)
      pop!(W.S₁) # Get a reset for this stack?
    end
    ResettableStacks.reset!(W.S₂)
  end
  setup_next && setup_next_step!(W)
end

function DiffEqBase.reinit!(W::AbstractNoiseProcess,dt;
                            t0 = W.t[1],
                            erase_sol = true,
                            setup_next = false)
  W.curt = t0
  W.dt = dt
  if typeof(W) <: NoiseGrid
    if isinplace(W)
      W.curW .= W.W[1]
      if W.Z != nothing
        W.curZ .= W.Z[1]
      end
    else
      W.curW = W.W[1]
      if W.Z != nothing
        W.curZ = W.Z[1]
      end
    end
    W.step_setup = true
  end
  setup_next && setup_next_step!(W)
end


text:
# Davis, Yin. "A Three-Operator Splitting Scheme and its Optimization
# Applications", Set-Valued and Variational Analysis, vol. 25, no. 4,
# pp. 829–858 (2017).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

struct DYS_iterable{R <: Real, C <: Union{R, Complex{R}}, T <: AbstractArray{C}, Tf, Tg, Th, TA}
    f::Tf
    g::Tg
    h::Th
    A::TA
    x0::T
    gamma::R
    lambda::R
end

mutable struct DYS_state{T, S}
    z::T
    xg::T
    y::S
    grad_h_y::S
    z_half::T
    xf::T
    res::T
end

function Base.iterate(iter::DYS_iterable{R, C, T}) where {R, C, T}
    z = copy(iter.x0)
    xg, = prox(iter.g, z, iter.gamma)
    y = iter.A * xg
    grad_h_y, = gradient(iter.h, y)
    z_half = 2 .* xg .- z .- iter.gamma .* (iter.A' * grad_h_y)

    # TODO: check the following
    xf = similar(z_half)
    prox!(xf, iter.f, z_half, iter.gamma)
    # xf, = prox(iter.f, z_half, iter.gamma)  # this fails the tests
    ############################################################################

    res = xf - xg
    z .+= iter.lambda .* res
    state = DYS_state{T, typeof(y)}(z, xg, y, grad_h_y, z_half, xf, res)
    return state, state
end

function Base.iterate(iter::DYS_iterable, state::DYS_state)
    prox!(state.xg, iter.g, state.z, iter.gamma)
    mul!(state.y, iter.A, state.xg)
    gradient!(state.grad_h_y, iter.h, state.y)
    mul!(state.z_half, iter.A', state.grad_h_y)
    state.z_half .*= -iter.gamma
    state.z_half .+= 2 .* state.xg .- state.z
    prox!(state.xf, iter.f, state.z_half, iter.gamma)
    state.res .= state.xf .- state.xg
    state.z .+= iter.lambda .* state.res
    return state, state
end

# Solver

struct DavisYin{R}
    gamma::Maybe{R}
    lambda::R
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int

    function DavisYin{R}(; gamma::Maybe{R}=nothing, lambda::R=R(1.0),
        maxit::Int=10000, tol::R=R(1e-8), verbose::Bool=false, freq::Int=100
    ) where R
        @assert gamma === nothing || gamma > 0
        @assert lambda > 0
        @assert maxit > 0
        @assert tol > 0
        @assert freq > 0
        new(gamma, lambda, maxit, tol, verbose, freq)
    end
end

function (solver::DavisYin{R})(x0::AbstractArray{C};
    f=Zero(), g=Zero(), h=Zero(), A=I, L::Maybe{R}=nothing
) where {R, C <: Union{R, Complex{R}}}

    stop(state::DYS_state) = norm(state.res, Inf) <= solver.tol
    disp((it, state)) = @printf("%5d | %.3e\n", it, norm(state.res, Inf))

    if solver.gamma === nothing
        if L !== nothing
            gamma = R(1)/L
        else
            error("You must specify either L or gamma")
        end
    else
        gamma = solver.gamma
    end

    iter = DYS_iterable(f, g, h, A, x0, gamma, solver.lambda)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose iter = tee(sample(iter, solver.freq), disp) end

    num_iters, state_final = loop(iter)

    return state_final.xf, state_final.xg, num_iters

end

# Outer constructors

"""
    DavisYin([gamma, lambda, maxit, tol, verbose, freq])

Instantiate the Davis-Yin splitting algorithm (see [1]) for solving
convex optimization problems of the form

    minimize f(x) + g(x) + h(A x),

where `h` is smooth and `A` is a linear mapping (for example, a matrix).
If `solver = DavisYin(args...)`, then the above problem is solved with

    solver(x0; [f, g, h, A])

Optional keyword arguments:

* `gamma::Real` (default: `nothing`), stepsize parameter.
* `labmda::Real` (default: `1.0`), relaxation parameter, see [1].
* `maxit::Integer` (default: `1000`), maximum number of iterations to perform.
* `tol::Real` (default: `1e-8`), absolute tolerance on the fixed-point residual.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `100`), frequency of verbosity.

If `gamma` is not specified at construction time, the following keyword
argument must be specified at solve time:

* `L::Real`, Lipschitz constant of the gradient of `h(A x)`.

References:

[1] Davis, Yin. "A Three-Operator Splitting Scheme and its Optimization
Applications", Set-Valued and Variational Analysis, vol. 25, no. 4,
pp. 829–858 (2017).
"""
DavisYin(::Type{R}; kwargs...) where R = DavisYin{R}(; kwargs...)
DavisYin(; kwargs...) = DavisYin(Float64; kwargs...)


text:
include("bonddimension.jl")

function optimizecontract(tensors::AbstractVector,labels::AbstractVector;maximalcost=-1,initialcost=0)
    length(tensors)==length(labels) || throw("number of tensors and number of labels don't match")
    numtensors=length(tensors)
    numtensors > 1 || throw(ArgumentError("a tensor network requires at least two tensors"))

    # determine type of spaces and of cost
    if all(x->isa(x,Tuple),tensors)
        T=mapreduce(x->promote_type(typeof(x)...),promote_type,tensors)
        if T<:Number || T<:Power
            P=(T,T...)
            T=typeof(T(0)+T(0))
            spaces=Array(P,numtensors)
            for n=1:numtensors
                spaces[n]=tensors[n]
            end
        else
            throw(ArgumentError("Unrecognized type of tensors"))
        end
    elseif all(x->isa(x,AbstractTensor),tensors)
        TT=mapreduce(typeof,promote_type,tensors)
        P=tensortype(TT)
        spaces=Array(P,numtensors)
        for n=1:numtensors
            spaces[n]=space(tensors[n])
        end
        T=Float64
    elseif all(x->isa(x,TensorSpace),tensors)
        P=mapreduce(typeof,promote_type,tensors)
        spaces=Array(P,numtensors)
        for n=1:numtensors
            spaces[n]=tensors[n]
        end
        T=Float64
    else
        throw(ArgumentError("Unrecognized type of tensors"))
    end

    # determine minimal and maximal cost
    if maximalcost==-1
        if T<:Poly
            maxcost=T(1,typemax(Int))
        else
            maxcost=typemax(T)
        end
    else
        maxcost=T(maximalcost)
    end
    if initialcost==0
        initcost=T(0)
    else
        initcost=T(initialcost)
    end

    return _optimizecontract(spaces,labels,initcost,maxcost)
end

# Auxiliary:
type SubGraph{T}
    tensorset::BitVector
    labelset::BitVector
    dualset::BitVector
    tree::Union{Int,Tuple}
    cost::T
end

function _optimizecontract{P,T}(spaces::Vector{P},labels::AbstractVector,initcost::T,maxcost::T)
    numtensors=length(labels)

    # Analyze contraction graph and detect components via adjacency matrix
    alllabels=unique(vcat(labels...))
    numlabels=length(alllabels)
    labelsets=Array(BitVector,numtensors) # map labels of tensor n to integers i and store in IntSet for faster set operations
    dualsets=Array(BitVector,numtensors) # store whether labelspaces[i] holds the space or its dual for tensor n

    tabletensor=zeros(Int,(numlabels,2))
    tableindex=zeros(Int,(numlabels,2))

    adjacencymatrix=falses(numtensors,numtensors)

    for n=1:numtensors
        length(labels[n])==length(unique(labels[n])) || throw(ArgumentError("duplicate labels for tensor $n: handle inner contraction first with tensortrace"))
        indn=findin(alllabels,labels[n])
        labelsets[n]=falses(numlabels)
        labelsets[n][indn]=true
        dualsets[n]=falses(numlabels)
        for i in indn
            if tabletensor[i,1]==0
                tabletensor[i,1]=n
                tableindex[i,1]=findfirst(labels[n],alllabels[i])
            elseif tabletensor[i,2]==0
                tabletensor[i,2]=n
                tableindex[i,2]=findfirst(labels[n],alllabels[i])
                n1=tabletensor[i,1]
                ind1=tableindex[i,1]
                ind2=tableindex[i,2]
                if isa(spaces[n1][ind1],Number) || isa(spaces[n1][ind1],Power)
                    spaces[n1][ind1]==spaces[n][ind2] || throw(ArgumentError("spaces don't match for label $(alllabels[i])"))
                else
                    spaces[n1][ind1]==dual(spaces[n][ind2]) || throw(ArgumentError("spaces don't match for label $(alllabels[i])"))
                end
                adjacencymatrix[n1,n]=true
                adjacencymatrix[n,n1]=true
                dualsets[n][i]=true
            else
                throw(ArgumentError("no label should appear more than two times"))
            end
        end
    end
    componentlist=connectedcomponents(adjacencymatrix)
    numcomponent=length(componentlist)

    # construct a list of the spaces of all labels
    labelspaces=P(ntuple(j->spaces[tabletensor[j,1]][tableindex[j]],numlabels))

    # generate output structures
    costlist=Array(T,numcomponent)
    treelist=Array(Tuple,numcomponent)
    labellist=Array(BitVector,numcomponent)

    # compute cost and optimal contraction order for every component
    for c=1:numcomponent
        # find optimal contraction for every component
        componentsize=length(componentlist[c])
        setcost=Array(Dict{BitVector,SubGraph{T}},componentsize)
        for k=1:componentsize
            setcost[k]=Dict{BitVector,SubGraph{T}}()
        end

        for n in componentlist[c]
            tensorset=falses(numtensors)
            tensorset[n]=true
            setcost[1][tensorset]=SubGraph{T}(tensorset,labelsets[n],dualsets[n],n,T(0))
        end
        currentcost=initcost
        while true
            nextcost=maxcost+1 # corresponds to uninitialized value
            for n=2:componentsize
                println("Component $c of $numcomponent: current cost = $currentcost, partition size $n out of $componentsize tensors")
                for k=1:div(n,2)
                    ksets=collect(values(setcost[k])) # subset of k tensors
                    if n-k==k
                        nmksets=ksets
                    else
                        nmksets=collect(values(setcost[n-k])) # subset of n-k tensors
                    end
                    for i=1:length(ksets)
                        for j=(n==2*k ? i+1 : 1):length(nmksets)
                            sc1=ksets[i]
                            sc2=nmksets[j]
                            if !anyand(sc1.tensorset,sc2.tensorset) # only select subsets with no common tensors
                                if anyand(sc1.labelset,sc2.labelset) # only select subsets with shared contraction lines
                                    clabels=sc1.labelset & sc2.labelset
                                    olabels1=sc1.labelset $ clabels
                                    olabels2=sc2.labelset $ clabels

                                    cost=sc1.cost+sc2.cost+contractcost(labelspaces,olabels1,olabels2,clabels,sc1.dualset,sc2.dualset)
                                    if cost<=currentcost
                                        tensorset=sc1.tensorset | sc2.tensorset
                                        if !haskey(setcost[n],tensorset)
                                            labelset=olabels1 | olabels2
                                            dualset=(sc1.dualset | sc2.dualset) & labelset
                                            setcost[n][tensorset]=SubGraph{T}(tensorset,labelset,dualset,(sc1.tree,sc2.tree),cost)
                                        else
                                            sc=setcost[n][tensorset]
                                            if sc.cost>cost
                                                sc.tree=(sc1.tree,sc2.tree)
                                                sc.cost=cost
                                            end
                                        end
                                    elseif cost<nextcost
                                        nextcost=cost
                                    end
                                end
                            end
                        end
                    end
                end
            end
            if nextcost>maxcost
                error("no solution with cost < maximalcost = $maxcost was found")
            end
            if !isempty(setcost[componentsize])
                break
            end
            currentcost=min(2*nextcost,maxcost)
        end
        component=falses(numtensors)
        component[componentlist[c]]=true
        sc=setcost[componentsize][component]
        costlist[c]=sc.cost
        treelist[c]=sc.tree
        labellist[c]=sc.labelset
    end

    # collect different components
    totaltree=treelist[1]
    totalcost=costlist[1]
    outputlabels=labellist[1]

    for c=2:numcomponent
        totaltree=(totaltree,treelist[c])
        totalcost=totalcost+costlist[c]+contractcost(labelspaces,outputlabels,labellist[c],IntSet(),IntSet(),IntSet())
        outputlabels=union(outputlabels,labellist[c])
    end

    return totaltree, totalcost
end

@inline function anyand(a::BitArray,b::BitArray)
    for n in 1:length(a)
        a[n] && b[n] && return true
    end
    return false
end
@inline function contractcost{D}(spaces::Tuple{Power{D,Int}, Vararg{Power{D,Int}}},olabels1::BitVector,olabels2::BitVector,clabels::BitVector,dualset1::BitVector,dualset2::BitVector)
    cost=Power{D,Int}(1,0)
    for n in find(olabels1 | olabels2 | clabels)
        cost*=spaces[n]
    end
    return cost
end
@inline function contractcost{T<:Number}(spaces::Tuple{T, Vararg{T}},olabels1::BitVector,olabels2::BitVector,clabels::BitVector,dualset1::BitVector,dualset2::BitVector)
    cost=T(0)
    for n in find(olabels1 | olabels2 | clabels)
        cost*=spaces[n]
    end
    return cost
end

function connectedcomponents(A::AbstractMatrix{Bool})
    # connectedcomponents(A::AbstractMatrix{Bool})
    #
    # For a given adjacency matrix of size n x n, connectedcomponents returns
    # a list componentlist that contains integer vectors, where every integer
    # vector groups the indices of the vertices of a connected component of the
    # graph encoded by A. The number of connected components is given by
    # length(componentlist).
    #
    # Used as auxiliary function to analyze contraction graph.

    n=size(A,1)
    assert(size(A,2)==n)

    componentlist=Array(Vector{Int},0)
    assignedlist=falses((n,))

    for i=1:n
        if !assignedlist[i]
            assignedlist[i]=true
            checklist=[i]
            currentcomponent=[i]
            while !isempty(checklist)
                j=pop!(checklist)
                for k=1:n
                    if A[j,k] && !assignedlist[k]
                        push!(currentcomponent,k)
                        push!(checklist,k)
                        assignedlist[k]=true;
                    end
                end
            end
            push!(componentlist,currentcomponent)
        end
    end
    return componentlist
end


text:
using Images: ImageShow
using Gtk: Window
using FileIO: action, open
using FileIO: action, save
using FileIO
using ImageIO
using Images
using ImageView

using Colors
using ImageMagick

using Gtk
using Gtk.ShortNames, GtkReactive

using LinearAlgebra
using FFTW
using ColorTypes
using Colors

include("bmpCompressor.jl")

# inFilePath = ARGS[1]
# d = ARGS[2]
# f = ARGS[3]
# outFilePath = inFilePath * "_" * d * "_" * f * ".bmp"
f_value = 0
d_value = 0
global file_global = ""

# Put your GUI code here

win = GtkWindow("BMP Compressor Tool", 300, 150)
vbox = GtkBox(:v)
push!(win, vbox)

f_text = GtkLabel("F value")
f_entry = GtkSpinButton(0,1000,1)
d_text = GtkLabel("D value")
d_entry = GtkSpinButton(0,1998,1)
filebox = GtkBox(:v)
file_text = GtkLabel("Choose a BMP file")

# file_entry = GtkFileChooserDialog(:action, open)
# set_gtk_property!(file_entry, :action, open)
choose_file = GtkButton("Pick a FIle")

compress_button = GtkButton("Compress")

push!(vbox, f_text)
push!(vbox, f_entry)
push!(vbox, d_text)
push!(vbox, d_entry)
push!(vbox, filebox)

push!(filebox,file_text)
push!(filebox, choose_file)

push!(vbox,compress_button)

showall(win)

function on_button_clicked(w)
    f_value = get_gtk_property(f_entry,:value,Int)
    d_value = get_gtk_property(d_entry,:value,Int)
    file = file_global
    #check f and d value
    # filename = get_gtk_property(file_entry, :filename, String)
    println("test")
    println("clicked  ", f_value, " and ", d_value)
    println(file)
    if (((2*f_value)-2) < d_value )
        println("Error!")
        errwin = GtkWindow("ERROR", 10, 5)
        box = GtkBox(:v)
        err = GtkLabel("D must be <= 2F-2")
        push!(box, err)
        push!(errwin, box)
        showall(errwin)

    else
        compressedImage = compress(file, f_value, d_value)
            # qui ho image out compressa
        imagestd = load(file)
        println("ok compressed")

    # VERSION1
        imshow(imagestd, name="Original")
        imshow(compressedImage, name = "Compressed")
    end

    # gui = imshow_gui((300, 300), (2, 1), name="Comparator")  # 2 columns, 1 row of images (each initially 300×300)
    # canvases = gui["canvas"]

    # im1 = Signal(imagestd[:,:,1])
    # im2 = Singal(compressedImage[:,:,1])

    # showall(gui["window"])

    # imshow(canvases[1, 1], im1)
    # imshow(canvases[1, 2], im2)

    # push!(im1, imagestd[:,:,1] )
    # push!(im2, compressedImage[:,:,1])

end

function on_choose_button_clicked(w)
    global file_global = open_dialog("Pick a BMP file", win, ("*.bmp",))
    println(file_global)

end

signal_connect(on_button_clicked, compress_button, "clicked")
signal_connect(on_choose_button_clicked, choose_file, "clicked")

if !isinteractive()
    c = Condition()
    signal_connect(win, :destroy) do widget
        notify(c)
    end
    @async Gtk.gtk_main()
    wait(c)
end

max_stars_repo_name: illoxian/bitmap-compressor-tool

max_issues_repo_name: illoxian/bitmap-compressor-tool

max_forks_repo_name: illoxian/bitmap-compressor-tool


text:
using LatBo.LB: collision, collision!

facts("Collision kernels") do
    context("Single relaxation time") do
        context("Linear vs τ⁻¹") do
            @fact 2collision(1, [1], [0]) --> collision(2, [1], [0])
            @fact 3collision(1, [0], [1]) --> collision(3, [0], [1])
        end
        context("Linear vs fᵢ (and feq)") do
            f = rand(Int64, 10)
            @fact 2collision(1, f, zeros(f)) --> collision(1, 2f, zeros(f))
            @fact 3collision(1, zeros(f), f) --> collision(1, zeros(f), 3f)
        end
        context("Anti-symmetric vs fᵢ and feq") do
            fᵢ, feq = rand(Int64, 10), rand(Int64, 10)
            @fact collision(1, fᵢ, feq) --> -collision(1, feq, fᵢ)
        end
        context("Lock-down the sign") do
            @fact collision(1, collect(1:10), zeros(Int64, 10)) .< 0 --> all
        end

        context("collision! <=--> collision") do
            τ⁻¹ = 3
            fᵢ = rand(Int64, (10, 5))
            fᵢ_copy = deepcopy(fᵢ)
            feq = rand(Int64, 10)
            for i in 1:size(fᵢ, 2)
                fᵢ_copy[:, i] += collision(τ⁻¹, fᵢ_copy[:, i], feq)
                collision!(fᵢ, τ⁻¹, feq, i)
                @fact fᵢ[:, i] --> fᵢ_copy[:, i] "$(fᵢ_copy[:, i])"
            end
        end
    end
end


text:
hascas(::Computed) = true  # maybe
maysync(::Computed) = true  # maybe
# TODO: check if `maysync(::Computed) == maysync(::Map) == true` is required/correct

function tryreact!(actr::Reactor{<:Computed}, a, rx::Reaction, offer::Union{Offer,Nothing})
    (; f) = actr.reagent
    r = f(a)::Union{Nothing,Failure,Reagent}
    if r === nothing
        return Block()
    elseif r isa Reagent
        return tryreact!(then(r, actr.continuation), a, rx, offer)
    else
        return r
    end
end

hascas(::Return) = false
maysync(::Return) = false  # See: maysync(::Map)

function tryreact!(actr::Reactor{<:Return}, _, rx::Reaction, offer::Union{Offer,Nothing})
    (; value) = actr.reagent
    value isa Failure && return value
    return tryreact!(actr.continuation, value, rx, offer)
end

struct ReturnIfBlocked{T} <: Reagent
    value::T
end

hascas(::ReturnIfBlocked) = false
maysync(::ReturnIfBlocked) = false  # See: maysync(::Map)

function tryreact!(
    actr::Reactor{<:ReturnIfBlocked},
    a,
    rx::Reaction,
    offer::Union{Offer,Nothing},
)
    (; value) = actr.reagent
    maysync(actr.continuation) || error("synchronizing continuation is required")
    ans = tryreact!(actr.continuation, a, rx, offer)
    @trace(
        label = :rib_tried,
        offerid = offerid(offer),
        taskid = objectid(current_task()),
        ans = ans,
        a,
        value,
    )
    if ans isa SomehowBlocked
        offer === nothing && return Block()  # require `offer` to try other branches
        if ans isa NeedNack
            rx = @set rx.postcommithooks = combine(rx.postcommithooks, ans.postcommithooks)
        end
        # NOTE: Using `Commit()` is OK only when there is no `Swap` in the
        # upstream reagents.  So, `ReturnIfBlocked` is only useful/usable with
        # the usage like `trysync!`.
        return tryreact!(Commit(), value, rx, offer)
    elseif ans isa Retry
        # Assuming the continuation has `Swap`, it must eventually return
        # `Block`.  So, retring the reaction should be fine.
        return ans
    else
        return ans
    end
end

hascas(::Map) = false
maysync(::Map) = false
# `Map`, `Return`, etc. may return `Block` for nudging reactions to try other
# branches (and also to indicate that `offer` is required). But these reagents
# themselves do not attempt to sync. It *seems* like retruning `false` in this
# case is more useful; see how `ReturnIfBlocked` use it for deadlock detection.

function tryreact!(actr::Reactor{<:Map}, a, rx::Reaction, offer::Union{Offer,Nothing})
    (; f) = actr.reagent
    b = f(a)
    b isa Failure && return b
    return tryreact!(actr.continuation, b, rx, offer)
end

Reagents.Map(::Type{T}) where {T} = Reagents.Map{Type{T}}(T)


text:
using Revise
using VCFTools
using MendelImpute
using GeneticVariation
using Random
using StatsBase
using CodecZlib
using ProgressMeter
using JLD2, FileIO, JLSO
using BenchmarkTools
using GroupSlices
using LinearAlgebra
using ThreadPools

function compress()
        max_d = 1000
	chr = 10
	minwidth = 0
	overlap = 0.0

        tgtfile = "target.chr$chr.typedOnly.masked.vcf.gz"
        reffile = "ref.chr$chr.excludeTarget.vcf.gz"
        outfile = "ref.chr$chr.excludeTarget.jlso"
        @time compress_haplotypes(reffile, tgtfile, outfile, max_d, minwidth, overlap)
end
compress()

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
using Test

@testset "BIC"  begin

    using StatsModelComparisons

    v1 = bic(-200.0, 10, 100)
    v2 = bic(-200.0, 2, 100)
    v3 = bic(-200.0, 2, 10)
    v4 = bic(-100.0, 2, 10)
    @test v2 < v1
    @test v3 < v2
    @test v4 < v3
end

@testset "AIC"  begin

    using StatsModelComparisons

    v1 = aic(-200.0, 10)
    v2 = aic(-200.0, 2)
    v3 = aic(-100.0, 2)
    @test v2 < v1
    @test v3 < v2
end

@testset "DIC"  begin

    using StatsModelComparisons
    x = randn(100)
    LLs = x .- mean(x)
    v1 = dic(LLs)
    v2 = dic(LLs*2)
    v3 = dic(LLs .+ 1)
    @test v1 < v2
    @test v3 < v1
end

@testset "WAIC and LOO" begin

    using StatsModelComparisons, StanSample
    using StatsFuns, RDatasets, Random, CSV
    import StanSample: read_csv_files

    cd(@__DIR__)
    Random.seed!(89905)
    # df = RDatasets.dataset("datasets", "cars")

    # stream = open("cars_model.stan", "r")
    # cars_stan = read(stream, String)
    # close(stream)

    # cars_stan_model = SampleModel("cars.model", cars_stan)
    # data = (N = size(df, 1), speed = df.Speed, dist = df.Dist)
    # rc = stan_sample(cars_stan_model; data)

    # stan_summary(cars_stan_model, true)
    # nt_cars = read_samples(cars_stan_model);
    # log_lik = nt_cars.log_lik'

    # files = readdir("test_chains")
    # files = map(f->"test_chains/"*f, files)
    # samples = read_csv_files(files, 2000, :namedtuple; include_internals=false)
    # df = vcat(dfs...)

    files = readdir("test_chains/")

    dfs = map(f->CSV.read("test_chains/"*f, DataFrame), files)
    dfs = map(x->x[1001:end,:], dfs)
    df = vcat(dfs...)

    col_idx = findall(x->occursin("log_lik", x), names(df))
    pw_log_liks = Array(df[:,col_idx])
    WAIC = waic(pw_log_liks)
    loo,_,_ = psisloo(pw_log_liks)

    @test WAIC.WAIC ≈ 421.0 rtol = .01
    @test WAIC.std_err ≈ 16.3 rtol = .01
    @test loo ≈ -210.5 rtol = .05
end


text:
#########################################################################
#
#   Julia Package for reverse mode automated differentiation (from source)
#
#########################################################################

module ReverseDiffSource

    import Base.show, Base.copy, Base.length

    using Compat  # for compatibility across julia versions

    # naming conventions
    const TEMP_NAME = "_tmp"   # prefix of new variables
    const DERIV_PREFIX = "d"   # prefix of gradient variables

    ## misc functions
    dprefix(v::Union{Symbol, AbstractString, Char}) = Symbol("$DERIV_PREFIX$v")

    isSymbol(ex)   = isa(ex, Symbol)
    isDot(ex)      = isa(ex, Expr) && ex.head == :.   && isa(ex.args[1], Symbol)
    isRef(ex)      = isa(ex, Expr) && ex.head == :ref && isa(ex.args[1], Symbol)

    ## temp var name generator
    let
        vcount = Dict()
        global newvar
        function newvar(radix::Union{AbstractString, Symbol}=TEMP_NAME)
            vcount[radix] = haskey(vcount, radix) ? vcount[radix]+1 : 1
            return Symbol("$(radix)$(vcount[radix])")
        end

        global resetvar
        function resetvar()
            vcount = Dict()
        end
    end

    ######  Includes  ######
    include("node.jl")
    include("bidict.jl")
    include("graph.jl")
    include("plot.jl")
    include("simplify.jl")
    include("tograph.jl")
    include("tocode.jl")
    include("zeronode.jl")
    include("reversegraph.jl")
    include("deriv_rule.jl")
    include("base_rules.jl")
    include("rdiff.jl")
    include("frdiff.jl")

    ######  Exports  ######
    export
        rdiff,
        @deriv_rule, deriv_rule

end # module ReverseDiffSource


text:
struct ImbensKalyanaraman <: BandwidthSelector end

_string(::ImbensKalyanaraman) = "Imbens Kalyanaraman bandwidth"

function kernel_constant(::ImbensKalyanaraman, kernel)
    kernel = EquivalentKernel(kernel)
    νs = [kernel_moment(kernel, Val(j)) for j = 0:3]
    πs = [squared_kernel_moment(kernel, Val(j)) for j = 0:3]
    C1 = 1 / 4 * abs2((abs2(νs[3]) - νs[2] * νs[4]) / (νs[3] * νs[1] - abs2(νs[2])))
    C2_num = abs2(νs[3]) * πs[1] - 2 * νs[2] * νs[3] * πs[2] + abs2(νs[2]) * πs[3]
    C2_denom = abs2(νs[3] * νs[1] - abs2(νs[2]))
    C2 = C2_num / C2_denom
    (C2 / (4 * C1))^(1 / 5)
end

function bandwidth(ik::ImbensKalyanaraman, kernel::SupportedKernels, ZsR::RDData)
    ZsR_untreated = ZsR[Untreated()]
    ZsR_treated = ZsR[Treated()]

    cutoff = ZsR.cutoff

    N = nobs(ZsR)

    N_untreated = nobs(ZsR_untreated)
    N_treated = nobs(ZsR_treated)

    # Step 1: Density and conditional variance at 0
    h₁ = 1.84 * Statistics.std(ZsR.Zs) * N^(-1 / 5)

    interval_h₁ = Interval{:closed,:closed}(cutoff-h₁, cutoff+h₁)

    ZsR_untreated_h₁ = ZsR_untreated[interval_h₁]
    ZsR_treated_h₁ = ZsR_treated[interval_h₁]

    N_h₁_untreated = nobs(ZsR_untreated_h₁)
    N_h₁_treated = nobs(ZsR_treated_h₁)

    Ȳ_h₁_untreated = mean(ZsR_untreated_h₁.Ys)
    sd_Y_h₁_untreated = Statistics.std(ZsR_untreated_h₁.Ys)
    Ȳ_h₁_treated = mean(ZsR_treated_h₁.Ys)
    sd_Y_h₁_treated = Statistics.std(ZsR_treated_h₁.Ys)

    f̂₀ = (N_h₁_untreated + N_h₁_treated) / 2 / N / h₁

    # Step 2: Estimation of second derivatives
    global_cubic_lm = lm(@formula(Ys ~ 1 + (ZsC >= 0) + ZsC + ZsC^2 + ZsC^3), ZsR)

    m̂₀_triple_prime = 6 * coef(global_cubic_lm)[5]

    h₂_untreated = 3.56 * (sd_Y_h₁_untreated^2 / f̂₀ / m̂₀_triple_prime^2)^(1 / 7) * N_untreated^(-1 / 7)
    h₂_treated =
        3.56 * (sd_Y_h₁_treated^2 / f̂₀ / m̂₀_triple_prime^2)^(1 / 7) * N_treated^(-1 / 7)

    ZsR_untreated_h₂ = ZsR_untreated[Interval{:closed,:closed}(cutoff-h₂_untreated, cutoff+h₂_untreated)]
    ZsR_treated_h₂ = ZsR_treated[Interval{:closed,:closed}(cutoff-h₂_treated, cutoff+h₂_treated)]

    N_h₂_untreated = nobs(ZsR_untreated_h₂)
    N_h₂_treated = nobs(ZsR_treated_h₂)

    quadratic_fit_untreated = lm(@formula(Ys ~ 1 + ZsC + ZsC^2), ZsR_untreated_h₂)
    m̂₀_double_prime_untreated = 2 * last(coef(quadratic_fit_untreated))

    quadratic_fit_treated = lm(@formula(Ys~ 1 + ZsC + ZsC^2), ZsR_treated_h₂)
    m̂₀_double_prime_treated  = 2 * last(coef(quadratic_fit_treated ))

    # Step 3: Calculation of regularization terms

    r̂_untreated = 2160 * sd_Y_h₁_untreated^2 / N_h₂_untreated/ h₂_untreated^4
    r̂_treated = 2160 * sd_Y_h₁_treated^2 / N_h₂_treated / h₂_treated^4

    Ckernel = kernel_constant(ik, kernel)# 3.4375

    regularized_squared_diff =
        abs2(m̂₀_double_prime_treated - m̂₀_double_prime_untreated) + r̂_untreated + r̂_treated
    ĥ_IK =
        Ckernel *
        ((sd_Y_h₁_untreated^2 + sd_Y_h₁_treated^2) / f̂₀ / regularized_squared_diff)^(1 / 5) *
        N^(-1 / 5)
    ĥ_IK
end


text:
insts = OrderedDict()
insts["minimal"] = [
    ((1, 2, 2, 2, true),),
    ((1, 2, 2, 2, false),),
    ]
insts["fast"] = [
    ((2, 2, 3, 2, true),),
    ((2, 2, 3, 2, false),),
    ((3, 3, 3, 3, true),),
    ((3, 3, 3, 3, false),),
    ((3, 3, 5, 4, true),),
    ((5, 2, 5, 3, true),),
    ((1, 30, 2, 30, true),),
    ((1, 30, 2, 30, false),),
    ((10, 1, 3, 1, true),),
    ((10, 1, 3, 1, false),),
    ]
return (PolyEnvelopeNative, insts)


text:
@testset "minimal_intersection_test" begin

    @test intersect(interval(1.0,3.0), interval(2.1,4.0)) === Interval(2.1,3.0)

    @test intersect(interval(1.0,3.0), interval(3.0,4.0)) === Interval(3.0,3.0)

    @test intersect(interval(1.0,3.0), emptyinterval()) === emptyinterval()

    @test intersect(entireinterval(), emptyinterval()) === emptyinterval()

    @test intersect(interval(1.0,3.0), entireinterval()) === Interval(1.0,3.0)

end

@testset "minimal_intersection_dec_test" begin

    @test intersect(DecoratedInterval(interval(1.0,3.0), com), DecoratedInterval(interval(2.1,4.0), com)) === DecoratedInterval(Interval(2.1,3.0), trv)

    @test intersect(DecoratedInterval(interval(1.0,3.0), dac), DecoratedInterval(interval(3.0,4.0), def)) === DecoratedInterval(Interval(3.0,3.0), trv)

    @test intersect(DecoratedInterval(interval(1.0,3.0), def), DecoratedInterval(emptyinterval(), trv)) === DecoratedInterval(emptyinterval(), trv)

    @test intersect(DecoratedInterval(entireinterval(), dac), DecoratedInterval(emptyinterval(), trv)) === DecoratedInterval(emptyinterval(), trv)

    @test intersect(DecoratedInterval(interval(1.0,3.0), dac), DecoratedInterval(entireinterval(), dac)) === DecoratedInterval(Interval(1.0,3.0), trv)

end

@testset "minimal_convex_hull_test" begin

    @test hull(interval(1.0,3.0), interval(2.1,4.0)) === Interval(1.0,4.0)

    @test hull(interval(1.0,1.0), interval(2.1,4.0)) === Interval(1.0,4.0)

    @test hull(interval(1.0,3.0), emptyinterval()) === Interval(1.0,3.0)

    @test hull(emptyinterval(), emptyinterval()) === emptyinterval()

    @test hull(interval(1.0,3.0), entireinterval()) === entireinterval()

end

@testset "minimal_convex_hull_dec_test" begin

    @test hull(DecoratedInterval(interval(1.0,3.0), trv), DecoratedInterval(interval(2.1,4.0), trv)) === DecoratedInterval(Interval(1.0,4.0), trv)

    @test hull(DecoratedInterval(interval(1.0,1.0), trv), DecoratedInterval(interval(2.1,4.0), trv)) === DecoratedInterval(Interval(1.0,4.0), trv)

    @test hull(DecoratedInterval(interval(1.0,3.0), trv), DecoratedInterval(emptyinterval(), trv)) === DecoratedInterval(Interval(1.0,3.0), trv)

    @test hull(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(emptyinterval(), trv)) === DecoratedInterval(emptyinterval(), trv)

    @test hull(DecoratedInterval(interval(1.0,3.0), trv), DecoratedInterval(entireinterval(), dac)) === DecoratedInterval(entireinterval(), trv)

end


text:
@testset "Test add_component Deprecation" begin
    data = IS.SystemData()

    component = IS.TestComponent("component1", 5)
    @test_deprecated IS.add_component!(data, component; deserialization_in_progress=true)
end


text:
# Include the known function
using Random
using Distributions
using IntervalSets
using SparseArrays

using CautiousEngine

# Get the absolute path of the experiment dir
EXPERIMENT_DIR = @__DIR__
experiment_type = "verification-example-single"
global_exp_dir = "$EXPERIMENT_DIR/$experiment_type/global"
specification_file = "foo"

#======================================================================================================
0. Define the Unknown Function, Process Noise, and the Full Dynamics Function
======================================================================================================#
system_mat = [0.4 0.1; 0. 0.5]
f = (x) -> system_mat*x

system = "mode1"
known_part = nothing

# Define the zero-mean process noise
σ_proc = 0.01
measurement_noise_dist = Normal(0., σ_proc)
process_noise_dist = nothing                                # Process noise is only used in synthesis, not verification

#======================================================================================================
1.  Perform GP regression and construct IMDP for each unknown mode.
======================================================================================================#
exp_dir = global_exp_dir
res_dirs = Array{String,1}()                        # Array that holds the single-mode result directories
X = Dict("x1" => [-2., 2.], "x2" => [-2., 2.])      # Defines the compact safe set
grid_sizes = Dict("x1" => 0.25, "x2" => 0.25)       # Specifies the discretization coarseness in each dimension
safety_dims = ["x1", "x2"]                          # Dimensions that have safety bounds
dependency_dims = Dict("x1" => [1, 1],              # Defines dependencies on other states
                       "x2" => [1, 1])
number_of_datapoints = 100                          # Number of data points per mode to generate
run_exps_flag = true                                # Set to true to rerun the single-mode constructions
m_opt = -1                                          # Unused
bound_type = "rkhs-tight"                           # GP bound type
epsilon = -1.                                       # Distance bound parameter, set to -1 for automatic determination

Lf_bound = 1.0                                  # Lipschitz constant bound for the mode over X
random_seed = 11                             	# Random seed for the random twister
eta = σ_proc                                    # η parameter set in a naiive manner

# Setup experiment parameter structures (defined in CautiousEngine.jl)
system_params = SystemParameters(system, f, known_part, measurement_noise_dist, process_noise_dist, Lf_bound, dependency_dims, 2, 2)
data_params = DataParameters(number_of_datapoints, m_opt, bound_type, σ_proc, epsilon, eta, safety_dims, -1., -1.)
experiment_params = ExperimentParameters(exp_dir, experiment_type, specification_file, X, grid_sizes, "foo", random_seed, system_params, data_params)

_, trans_mats = CautiousEngine.end_to_end_transition_bounds(experiment_params, reuse_gps_flag=true, reuse_regions_flag=true)
results_path = experiment_params.experiment_directory

#======================================================================================================
2.  Construct the IMDP Labels and call the appropriate verification solver wrapper.
======================================================================================================#
obs_extents = [Dict("x1" => -2.0.. -1.0, "x2" => -2.0.. -1.0),
                   Dict("x1" => -1.0..0.0, "x2" => 0.0..1.0),
                   Dict("x1" => 0.0..1.0, "x2" => -1.0..0.0),
                   Dict("x1" => 1.0..2.0, "x2" => 1.0..2.0)]
des_extents = [Dict("x1" => -1.0..0.0, "x2" => -1.0..0.0),
			   Dict("x1" => 0.0..1.0, "x2" => 0.0..1.0)]

labels_dict = Dict("!a∧b" => obs_extents, "a∧!b" => des_extents)

lbl = (state_extent; unsafe=false) -> general_label_fcn(state_extent, "!a∧!b", "!a∧b", labels_dict, unsafe=unsafe)

imdp = CautiousEngine.create_simple_imdp(trans_mats["minPr"], trans_mats["maxPr"])
CautiousEngine.create_imdp_labels(lbl, imdp, "$results_path/regions.bson")

k_vals = [1]
for k in k_vals
	verification_result_mat = CautiousEngine.BoundedUntil(imdp, "!a∧!b", "a∧!b", k, "$results_path/imdp.txt")
	CautiousEngine.save_legacy_mats(verification_result_mat, results_path, k)
end

annotation_dict = Dict()
annotation_dict["O"] = labels_dict["!a∧b"]
annotation_dict["D"] = labels_dict["a∧!b"]
extents_labels_colors = Dict()
extents_labels_colors["O"] = :white
extents_labels_colors["D"] = :black
CautiousEngine.plot_2d_verification_results(results_path, extents_dict=annotation_dict, extents_labels_colors=extents_labels_colors)


text:
#############
# SddMgr
#############

"Root of the SDD manager node hierarchy"
abstract type SddMgrNode <: VtreeNode end

const SddMgr = AbstractVector{<:SddMgrNode}

#############
# SddMgrNode
#############

"Root of the SDD circuit node hierarchy"
abstract type SddNode{V<:SddMgrNode} <: StructLogicalΔNode{V} end

"A SDD logical leaf node"
abstract type SddLeafNode{V} <: SddNode{V} end

"A SDD logical inner node"
abstract type SddInnerNode{V} <: SddNode{V} end

"A SDD logical literal leaf node, representing the positive or negative literal of its variable"
mutable struct SddLiteralNode{V} <: SddLeafNode{V}
    literal::Lit
    vtree::V
    data
    bit::Bool
    SddLiteralNode(l,v::V) where V = SddLiteralNode{V}(l,v)
    SddLiteralNode{V}(l,v::V) where V = new{V}(l,v,nothing,false)
end

"""
A structured logical constant leaf node, representing true or false.
These are the only structured nodes that don't have an associated vtree node (cf. SDD file format)
"""
abstract type SddConstantNode{V} <: SddInnerNode{V} end
mutable struct SddTrueNode{V} <: SddConstantNode{V}
    data
    bit::Bool
    SddTrueNode{V}() where V = new{V}(nothing,false)
end
mutable struct SddFalseNode{V} <: SddConstantNode{V}
    data
    bit::Bool
    SddFalseNode{V}() where V = new{V}(nothing,false)
end

"A structured logical conjunction node"
mutable struct Sdd⋀Node{V} <: SddInnerNode{V}
    prime::SddNode{<:V}
    sub::SddNode{<:V}
    vtree::V
    data
    bit::Bool
    Sdd⋀Node{V}(p,s,v::V) where V = new{V}(p,s,v,nothing,false)
end

"A structured logical disjunction node"
mutable struct Sdd⋁Node{V} <: SddInnerNode{V}
    children::Vector{<:Sdd⋀Node{<:V}}
    vtree::V
    data
    bit::Bool
    negation::Sdd⋁Node{V}
    Sdd⋁Node{V}(ch,v::V) where V = new{V}(ch, v, nothing, false) # leave negation uninitialized
    Sdd⋁Node{V}(ch,v::V,neg) where V = new{V}(ch, v, nothing, false, neg)
end

SddHasVtree = Union{Sdd⋁Node,Sdd⋀Node,SddLiteralNode}

"A structured logical circuit represented as a bottom-up linear order of nodes"
const Sdd{V} = AbstractVector{<:SddNode{<:V}} where {V <: SddMgrNode}

Base.eltype(::Type{Sdd}) = SddNode

#####################
# traits
#####################

@inline GateType(::Type{<:SddLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:SddConstantNode}) = ConstantGate()
@inline GateType(::Type{<:Sdd⋀Node}) = ⋀Gate()
@inline GateType(::Type{<:Sdd⋁Node}) = ⋁Gate()

#####################
# methods
#####################

@inline literal(n::SddLiteralNode)::Lit = n.literal
@inline constant(::SddTrueNode)::Bool = true
@inline constant(::SddFalseNode)::Bool = false
@inline children(n::Sdd⋀Node) = [n.prime,n.sub]
@inline children(n::Sdd⋁Node) = n.children

@inline vtree(n::SddHasVtree) = n.vtree

# alias some SDD terminology: primes and subs
@inline prime(n::Sdd⋀Node)::SddNode = n.prime
@inline sub(n::Sdd⋀Node)::SddNode = n.sub

Base.show(io::IO, c::SddTrueNode) = print(io, "⊤")
Base.show(io::IO, c::SddFalseNode) = print(io, "⊥")
Base.show(io::IO, c::SddLiteralNode) = print(io, literal(c))
Base.show(io::IO, c::Sdd⋀Node) = begin
    recshow(c::Union{SddConstantNode,SddLiteralNode}) = "$c"
    recshow(c::Sdd⋁Node) = "D$(hash(c))"
    print(io, "($(recshow(prime(c))),$(recshow(sub(c))))")
end
Base.show(io::IO, c::Sdd⋁Node) = begin
    elems = ["$e" for e in children(c)]
    print(io, "[$(join(elems,','))]")
end


text:
# Quadratic terms & constraints
#

function add_qpterms!(model::Model, qr::IVec, qc::IVec, qv::FVec)
    n = num_var(model)
    ((m = length(qr)) == length(qc) == length(qv)) || error("Inconsistent argument dimensions.")
    nqv = copy(qv)
    Q = sparse(qr, qc, nqv, n, n)
    if istriu(Q) || istril(Q) || issymmetric(Q)
        if VERSION >= v"0.7.0-DEV.3382"
            diag_matrix = spdiagm(0 => diag(Q))
        else
            diag_matrix = spdiagm(diag(Q))
        end
        Q = Q + Q' - diag_matrix # reconstruct full matrix like CPLEX wants
    else
        error("Matrix Q must be either symmetric or triangular")
    end
    qmatcnt = Vector{Cint}(undef, n)
    for k = 1:n
      qmatcnt[k] = Q.colptr[k+1] - Q.colptr[k]
    end
    stat = @cpx_ccall(copyquad, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cint},
                      Ptr{Cint},
                      Ptr{Cint},
                      Ptr{Cdouble}
                      ),
                      model.env.ptr, model.lp, convert(Array{Cint,1}, Q.colptr[1:end-1].-1), convert(Array{Cint,1},qmatcnt), convert(Array{Cint,1}, Q.rowval.-1), Q.nzval)
    if stat != 0
        throw(CplexError(model.env, stat))
        end
    model.has_qc = true
    nothing
end

function add_qpterms!(model::Model, qr::Vector, qc::Vector, qv::Vector)
    add_qpterms!(model, ivec(qr), ivec(qc), fvec(qv))
end

function add_qpterms!(model, H::SparseMatrixCSC{Float64}) # H must be symmetric
    n = num_var(model)
    (H.m == n && H.n == n) || error("H must be an n-by-n symmetric matrix.")

    nnz_h = nnz(H)
    qr = Vector{Cint}(undef, nnz_h)
    qc = Vector{Cint}(undef, nnz_h)
    qv = Vector{Float64}(undef, nnz_h)
    k = 0

    colptr::Vector{Int} = H.colptr
    nzval::Vector{Float64} = H.nzval

    for i = 1 : n
        qi::Cint = convert(Cint, i)
        for j = colptr[i]:(colptr[i+1]-1)
            qj = convert(Cint, H.rowval[j])

            if qi <= qj
                k += 1
                qr[k] = qi
                qc[k] = qj
                qv[k] = nzval[j]
            end
        end
    end

    add_qpterms!(model, qr[1:k], qc[1:k], qv[1:k])
end

function add_qpterms!(model, H::Matrix{Float64}) # H must be symmetric
    n = num_var(model)
    size(H) == (n, n) || error("H must be an n-by-n symmetric matrix.")

    nmax = div(n * (n + 1), 2)
    qr = Vector{Cint}(undef, nmax)
    qc = Vector{Cint}(undef, nmax)
    qv = Vector{Float64}(undef, nmax)
    k::Int = 0

    for i = 1 : n
        qi = convert(Cint, i)
        for j = i : n
            v = H[j, i]
            if v != 0.
                k += 1
                qr[k] = qi
                qc[k] = convert(Cint, j)
                qv[k] = v
            end
        end
    end

    add_qpterms!(model, qr[1:k], qc[1:k], qv[1:k])
end

function add_diag_qpterms!(model, H::Vector)  # H stores only the diagonal element
    n = num_var(model)
    n == length(H) || error("Incompatible dimensions.")
    q = [convert(Cint,1):convert(Cint,n)]
    add_qpterms!(model, q, q, fvec(h))
end

function add_diag_qpterms!(model, hv::Real)  # all diagonal elements are H
    n = num_var(model)
    q = [convert(Cint,1):convert(Cint,n)]
    add_qpterms!(model, q, q, fill(float64(hv), n))
end

# add_qconstr!

function add_qconstr!(model::Model, lind::IVec, lval::FVec, qr::IVec, qc::IVec, qv::FVec, rel::Cchar, rhs::Float64)
    qnnz = length(qr)
    qnnz == length(qc) == length(qv) || error("Inconsistent argument dimensions.")

    lnnz = length(lind)
    lnnz == length(lval) || error("Inconsistent argument dimensions.")

    if qnnz > 0 || lnnz > 0
        stat = @cpx_ccall(addqconstr, Cint, (
                          Ptr{Cvoid},   # env
                          Ptr{Cvoid},   # model
                          Cint,         # lnnz
                          Cint,         # qnnz
                          Float64,      # rhs
                          Cchar,        # sense
                          Ptr{Cint},    # lind
                          Ptr{Float64}, # lval
                          Ptr{Cint},    # qrow
                          Ptr{Cint},    # qcol
                          Ptr{Float64}, # qval
                          Ptr{UInt8}    # name
                          ),
                          model.env.ptr, model.lp, lnnz, qnnz, rhs, rel,
                          lind .- Cint(1), lval, qr .- Cint(1), qc .- Cint(1),
                          qv, C_NULL)
        if stat != 0
            throw(CplexError(model.env, stat))
        end
        model.has_qc = true
    end
    nothing
end

const sensemap = Dict('=' => 'E', '<' => 'L', '>' => 'G')
function add_qconstr!(model::Model, lind::Vector, lval::Vector, qr::Vector, qc::Vector, qv::Vector{Float64}, rel::GChars, rhs::Real)
    add_qconstr!(model, ivec(lind), fvec(lval), ivec(qr), ivec(qc), fvec(qv), cchar(sensemap[rel]), float(rhs))
end

function num_qconstr(model::Model)
    return @cpx_ccall(getnumqconstrs, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
                      model.env.ptr, model.lp)
end

function c_api_getqconstr(model::Model, row::Int)
    # In the first call, we ask CPLEX how many non-zero elements there are in
    # the affine (-linsurplus_p ) and quadratic (-quadsurplus_p) components.
    rhs_p = Ref{Cdouble}()
    sense_p = Ref{Cchar}()
    linsurplus_p = Ref{Cint}()
    quadsurplus_p = Ref{Cint}()
    stat = @cpx_ccall(
        getqconstr,
        Cint, (
            Ptr{Cvoid}, Ptr{Cvoid},
            Ptr{Cint}, Ptr{Cint}, Ptr{Float64}, Ptr{Cchar},
            Ptr{Cint}, Ptr{Float64}, Cint, Ptr{Cint},
            Ptr{Cint}, Ptr{Cint}, Ptr{Float64}, Cint, Ptr{Cint},
            Cint),
        model.env.ptr, model.lp,
        C_NULL, C_NULL, rhs_p, sense_p,
        C_NULL, C_NULL, 0, linsurplus_p,
        C_NULL, C_NULL, C_NULL, 0, quadsurplus_p,
        Cint(row-1))
    # In the second call, we initialize arrays to contain the number of non-zero
    # elements computed in the first part and then actually query the
    # coefficients.
    linspace = -linsurplus_p[]
    quadspace = -quadsurplus_p[]
    linind = fill(Cint(0), linspace)
    linval = fill(Cdouble(0.0), linspace)
    quadrow = fill(Cint(0), quadspace)
    quadcol = fill(Cint(0), quadspace)
    quadval = fill(Cdouble(0.0), quadspace)
    linnzcnt_p = Ref{Cint}()
    quadnzcnt_p = Ref{Cint}()
    stat = @cpx_ccall(
        getqconstr,
        Cint, (
            Ptr{Cvoid}, Ptr{Cvoid},
            Ptr{Cint}, Ptr{Cint}, Ptr{Float64}, Ptr{Cchar},
            Ptr{Cint}, Ptr{Float64}, Cint, Ptr{Cint},
            Ptr{Cint}, Ptr{Cint}, Ptr{Float64}, Cint, Ptr{Cint},
            Cint),
        model.env.ptr, model.lp,
        linnzcnt_p, quadnzcnt_p, rhs_p, sense_p,
        linind, linval, linspace, linsurplus_p,
        quadrow, quadcol, quadval, quadspace, quadsurplus_p,
        Cint(row-1))
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    if quadsurplus_p[] < 0 || linsurplus_p[] < 0
        error("Unable to query quadratic constraint, there were more " *
              "non-zero elements than expected.")
    end
    return linind, linval, quadrow, quadcol, quadval, sense_p[], rhs_p[]
end

function c_api_getquad(model::Model)
    num_variables = num_var(model)
    # In the first call, we ask CPLEX how many non-zero elements there are.
    surplus_p = Ref{Cint}()
    stat = @cpx_ccall(
        getquad,
        Cint, (
             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cint},
             Ptr{Cint}, Ptr{Cint}, Ptr{Float64},
             Cint, Ptr{Cint}, Cint, Cint),
             model.env.ptr, model.lp, C_NULL,
             C_NULL, C_NULL, C_NULL,
             0, surplus_p, 0, num_variables - 1)
    # In the second call, we initialize arrays to contain the number of
    # non-zero elements computed in the first part and then actually query the
    # coefficients.
    nzcnt_p = Ref{Cint}()
    qmatbeg = fill(Cint(0), num_variables)
    qmatind = fill(Cint(0), -surplus_p[])
    qmatval = fill(0.0, -surplus_p[])
    stat = @cpx_ccall(
        getquad,
        Cint, (
            Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cint},
            Ptr{Cint}, Ptr{Cint}, Ptr{Float64},
            Cint, Ptr{Cint}, Cint, Cint),
        model.env.ptr, model.lp, nzcnt_p,
        qmatbeg, qmatind, qmatval,
        -surplus_p[], surplus_p, 0, num_variables - 1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    if surplus_p[] < 0 || nzcnt_p[] != length(qmatind)
        error("Unable to query quadratic constraint, there were more " *
              "non-zero elements than expected.")
    end
    return qmatbeg, qmatind, qmatval
end

# int CPXdelqconstrs( CPXCENVptr env, CPXLPptr lp, int begin, int end )
function c_api_delqconstrs(model::Model,  first::Cint, last::Cint)
    stat = @cpx_ccall(
        delqconstrs,
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint),
        model.env.ptr, model.lp, first, last
    )
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return
end


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
import StatsBase: fit!, predict
using StatsBase
using BSON, FileIO
using Flux
using GenerativeModels
using DistributionsAD
using ValueHistories

s = ArgParseSettings()
@add_arg_table! s begin
	"modelname"
		default = "vae"
		arg_type = String
		help = "model name"
	"dataset"
		default = "iris"
		arg_type = String
		help = "dataset"
	"--seed"
		default = nothing
		help = "if specified, only results for a given seed will be recomputed"
	"--anomaly_class"
		default = nothing
		help = "if specified, only results for a given anomaly class will be recomputed"
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, modelname, seed, anomaly_class = parsed_args

masterpath = datadir("experiments/images_leave-one-in/$(modelname)/$(dataset)")
files = GenerativeAD.Evaluation.collect_files(masterpath)
mfiles = filter(f->occursin("model", f), files)
if seed != nothing
	filter!(x->occursin("/seed=$seed/", x), mfiles)
end
if anomaly_class != nothing
	filter!(x->occursin("/ac=$(anomaly_class)/", x), mfiles)
end

sample_score_batched(m,x,L,batchsize) =
	vcat(map(y-> Base.invokelatest(GenerativeAD.Models.reconstruction_score, m, y, L),
		Flux.Data.DataLoader(x, batchsize=batchsize))...)
sample_score_batched_gpu(m,x,L,batchsize) =
	vcat(map(y-> cpu(Base.invokelatest(GenerativeAD.Models.reconstruction_score, m,
		gpu(Array(y)), L)), Flux.Data.DataLoader(x, batchsize=batchsize))...)

function save_sample_score(f::String, data, seed::Int, ac=nothing)
	# get model
	savepath = dirname(f)
	mdata = load(f)
	model = mdata["model"]
	L = 100

	# setup entries to be saved
	save_entries = (
		modelname = modelname,
		fit_t = mdata["fit_t"],
		history = mdata["history"],
		dataset = dataset,
		npars = sum(map(p->length(p), Flux.params(model))),
		model = nothing,
		seed = seed
		)
	save_entries = (ac == nothing) ? save_entries : merge(save_entries, (anomaly_class=ac,))
	result = (x -> sample_score_batched_gpu(gpu(model), x, L, 256),
		merge(mdata["parameters"], (L = L, score = "reconstruction-sampled")))

	# if the file does not exist already, compute the scores
	savef = joinpath(savepath, savename(result[2], "bson", digits=5))
	if !isfile(savef)
		@info "computing sample score for $f"
		GenerativeAD.experiment(result..., data, savepath; save_entries...)
	end
end

for f in mfiles
	# get data
	savepath = dirname(f)
	local seed = parse(Int, replace(basename(savepath), "seed=" => ""))
	ac = occursin("ac=", savepath) ? parse(Int, replace(basename(dirname(savepath)), "ac=" => "")) : nothing
	data = GenerativeAD.load_data(dataset, seed=seed, anomaly_class_ind=ac,
		method="leave-one-in")

	# compute and save the score
	save_sample_score(f, data, seed, ac)
end
@info "DONE"


text:
module TestDeprecated

using CUDA
using FoldsCUDA
using Test
using Transducers

function test_foldx_cuda()
    result = Ref{Any}(nothing)
    xs = CUDA.ones(3)
    @test_deprecated result[] = foldx_cuda(TeeRF(min, max), xs)
    if result[] !== nothing
        @test result[] == (1, 1)
    end
end

function test_transduce_cuda()
    result = Ref{Any}(nothing)
    xs = CUDA.ones(3)
    @test_deprecated result[] =
        transduce_cuda(Map(identity), TeeRF(min, max), (0.0f0, 2.0f0), xs)
    if result[] !== nothing
        @test result[] == (0, 2)
    end
end

end  # module


text:
module BumpOnTail

using Parameters
using ..ParticleMethods: draw_g_accept_reject, draw_g_importance_sampling

# x-part of distribution function
function fₓ(x::T, params::NamedTuple) where {T}
    @unpack κ, ε = params
    @. ( 1 - ε * cos(κ*x) )
end

# v-part of distribution function
function fᵥ(v::T, params::NamedTuple) where {T}
    @unpack a, v₀, σ = params
    @. ( (1-a) * exp( - v^2 / 2 ) + a / σ * exp( - (v-v₀)^2 / (2σ^2) ) ) / √(2π)
end

# distribution function
function f(x::T, v::T, params) where {T}
    fₓ(x, params) .* fᵥ(v, params)
end

draw_accept_reject(Nₚ, params) = draw_g_accept_reject(Nₚ, fₓ, params)
draw_importance_sampling(Nₚ, params) = draw_g_importance_sampling(Nₚ, fₓ, params)

end


text:
@testset "Linear Spaces" begin
    @testset "LinearSubspace" begin
        A = LinearSubspace([1 0 3; 2 1 3], [5, -2])
        @test dim(A) == dim(intrinsic(A)) == dim(extrinsic(A)) == 1
        @test codim(A) == codim(intrinsic(A)) == codim(extrinsic(A)) == 2
        @test ambient_dim(A) == 3
        @test startswith(sprint(show, A), "1-dim. affine linear subspace")
        @test identity.(A) == A

        @test intrinsic(A) isa IntrinsicDescription
        @test startswith(sprint(show, intrinsic(A)), "IntrinsicDescription")
        @test extrinsic(A) isa ExtrinsicDescription
        @test startswith(sprint(show, extrinsic(A)), "ExtrinsicDescription")

        @test identity.(extrinsic(A)) == extrinsic(A)
        @test identity.(intrinsic(A)) == intrinsic(A)
        @test identity.(Intrinsic) == Intrinsic

        u = rand(1)
        x = A(u, Intrinsic)
        @test coord_change(A, Extrinsic, Intrinsic, x) ≈ u rtol = 1e-12
        @test coord_change(A, Intrinsic, Extrinsic, u) ≈ x rtol = 1e-12
        @test norm(A(x, Extrinsic)) ≈ 0 atol = 1e-14
        @test norm(ExtrinsicDescription(intrinsic(A))(x)) ≈ 0.0 atol = 1e-12

        B = rand_subspace(3; codim = 2)
        @test B != A
        copy!(B, A)
        @test extrinsic(B) == extrinsic(A)
        @test intrinsic(B) == intrinsic(A)
        @test B == A
        @test copy(A) == A
        @test copy(A) !== A

        C = rand_subspace(3, dim = 1, real = true)
        @test geodesic_distance(A, C) > 0
        γ = geodesic(A, C)
        @test size(γ(1)) == size(intrinsic(A).Y)

        A2 = translate(A, [1, 1], Extrinsic)
        A3 = LinearSubspace(extrinsic(A).A, extrinsic(A).b + [1, 1])
        @test A2.intrinsic.Y ≈ A3.intrinsic.Y

        # rand subspace through a point
        x = randn(ComplexF64, 5)
        L = rand_subspace(x; dim = 2)
        @test norm(L(x)) ≈ 0 atol = 1e-8
        L2 = rand_subspace(x; dim = 2, affine = false)
        @test is_linear(L2)
        @test norm(L2(x)) ≈ 0 atol = 1e-8
    end

    @testset "Intersect subspaces" begin
        L₁ = rand_subspace(7; codim = 2)
        L₂ = rand_subspace(7; codim = 3)
        L₃ = L₁ ∩ L₂
        @test codim(L₃) == 5
        @test dim(L₃) == 2
        @test ambient_dim(L₃) == 7
        E₃ = extrinsic(L₃)
        @test norm(L₁(E₃.A \ E₃.b, Extrinsic)) ≈ 0 atol = 1e-12
        @test norm(L₂(E₃.A \ E₃.b, Extrinsic)) ≈ 0 atol = 1e-12
        @test norm(L₃(E₃.A \ E₃.b, Extrinsic)) ≈ 0 atol = 1e-12
    end

    @testset "IntrinsicSubspaceHomotopy" begin
        @var x[1:4]
        f1 = rand_poly(x, 2)
        f2 = rand_poly(x, 2)
        F = System([f1, f2], x)

        A = rand_subspace(4; codim = 2)
        # Commpute witness set
        L = A(x, Extrinsic)
        res = track.(total_degree(F ∩ L)...)
        W = solution.(res)

        B = rand_subspace(4, codim = 2)

        graff_tracker = Tracker(
            IntrinsicSubspaceHomotopy(F, A, B),
            options = (automatic_differentiation = 3,),
        )
        graff_result = track.(graff_tracker, W)
        @test all(is_success, graff_result)

        C = rand_subspace(4, codim = 2)
        copy_A = copy(A)
        set_subspaces!(graff_tracker.homotopy, B, C)
        @test A == copy_A
        @test all(is_success, track.(graff_tracker, solution.(graff_result)))

        graff_path_tracker = EndgameTracker(IntrinsicSubspaceHomotopy(F, A, B))
        graff_path_result = track.(graff_path_tracker, W)
        @test all(is_success, graff_path_result)

        @test all(
            isapprox.(solution.(graff_path_result), solution.(graff_result), rtol = 1e-12),
        )
    end
end


text:
import Pkg

##Pkg.add("Optim")
Pkg.add("ForwardDiff")
Pkg.add("XMLDict")
Pkg.add("Plots")
Pkg.add("PyPlot")
Pkg.add("FFTW")
Pkg.add("GZip")
Pkg.add("SpecialFunctions")
Pkg.add("Suppressor")

##Pkg.add("Libz")
#Pkg.add("JLD2")
#Pkg.add("EzXML")
#Pkg.add("FileIO")


text:
# The Haar measure is a left- or right-invariant measure on a group manifold.
# That is, given p ∈ G and q ∈ H ⊆ G, where τ(H, p) is the set of all q τ-translated
# (left- or right-) by p, then μ(τ(H, p)) = μ(H). This is the τ-invariant Haar measure
# on G.

"""
    Haar{G,D}

The Haar measure on a group manifold.

The Haar measure on a group ``G`` is a measure that is invariant to all left- and/or right-
group translations on the manifold. That is, given the ``τ``-invariant Haar measure ``μ``
and some subset ``H ⊆ G``, then ``∫_H \\mathrm{d}μ(p) = ∫_H \\mathrm{d}μ(τ_q(p))`` for all
``q ∈ G``, where ``τ_q(p)`` is group translation of ``p`` by ``q`` in the ``τ`` direction
(left- or right-).

The convenient aliases [`LeftHaar`](@ref) and [`RightHaar`](@ref) are also provided.

# Constructors

    Haar(G::AbstractManifold, D::ActionDirection = LeftAction())

Construct the `D`-invariant Haar measure on group manifold `G`.
"""
struct Haar{G,D} <: PrimitiveMeasure
    group::G
    direction::D
end
Haar(G) = Haar(G, LeftAction())

"""
    LeftHaar{G<:AbstractManifold}

Alias for a left-[`Haar`](@ref) measure.
"""
const LeftHaar{G} = Haar{G,LeftAction}
LeftHaar(G) = Haar(G, LeftAction())

"""
RightHaar{G<:AbstractManifold}

Alias for a right-[`Haar`](@ref) measure.
"""
const RightHaar{M} = Haar{M,RightAction}
RightHaar(G) = Haar(G, RightAction())

function Base.show(io::IO, mime::MIME"text/plain", μ::LeftHaar)
    print(io, "LeftHaar(")
    show(io, mime, base_manifold(μ))
    return print(io, ")")
end
function Base.show(io::IO, mime::MIME"text/plain", μ::RightHaar)
    print(io, "RightHaar(")
    show(io, mime, base_manifold(μ))
    return print(io, ")")
end

Manifolds.base_manifold(μ::Haar) = μ.group
Manifolds.base_group(μ::Haar) = μ.group
Manifolds.direction(μ::Haar) = μ.direction

MeasureTheory.logdensity(::Haar, x) = zero(eltype(x))


text:
using DifferentialEquations, Plots, LaTeXStrings

##
function Hénon_Heiles(du,u,p,t)
    x  = u[1]
    y  = u[2]
    dx = u[3]
    dy = u[4]
    du[1] = dx
    du[2] = dy
    du[3] = -x - 2x*y
    du[4] = y^2 - y -x^2
end

##
initial1 = [1/6,0.0,-0.2,0.0]
initial2 = [1/6,0.0 + 0.01,-0.2+0.01,0.0]

##
e,x,y,ydot = initial1
xdot = sqrt(2*e - ydot^2 - y^2 + (2//3)*y^3)
ini = [x,y,xdot,ydot]
tspan = (0,150.)
prob = ODEProblem(Hénon_Heiles, ini, tspan)
sol = solve(prob, Vern9(), abs_tol=1e-16, rel_tol=1e-16,saveat = 0.1);
plot(sol, vars=(1,2), title = "The orbit of the Hénon-Heiles system", xaxis = "x", yaxis = "y", label = false,linecolor = :blue)
scatter!([x],[y],label = false,markercolor = :blue)

##
e2,x2,y2,ydot2 = initial2
xdot2 = sqrt(2*e2 - ydot2^2 - y2^2 + (2//3)*y2^3)
ini2 = [x2,y2,xdot2,ydot2]
prob2 = ODEProblem(Hénon_Heiles, ini2, tspan)
sol2 = solve(prob2, Vern9(), abs_tol=1e-16, rel_tol=1e-16,saveat = 0.1);
plot!(sol2, vars=(1,2), label = false,linecolor = :red)
scatter!([x2],[y2],label = false,markercolor = :red)

max_stars_repo_name: Ved-Mahajan/Henon-Heiles-System

max_issues_repo_name: Ved-Mahajan/Henon-Heiles-System

max_forks_repo_name: Ved-Mahajan/Henon-Heiles-System