text:
# Adios functions

export Adios
"""
    mutable struct Adios

Holds a C pointer `adios2_adios *`.

This value is finalized automatically. It can also be explicitly
finalized by calling `finalize(adios)`.
"""
mutable struct Adios
    ptr::Ptr{Cvoid}
    Adios(ptr) = finalizer(adios_finalize, new(ptr))
end
Adios() = Adios(C_NULL)

function Base.show(io::IO, ::MIME"text/plain", adios::Adios)
    return print(io, "Adios(0x$(string(UInt(adios.ptr); base=16)))")
end

export adios_init_mpi
"""
    adios = adios_init_mpi(comm::MPI.Comm)
    adios = adios_init_mpi(config_file::AbstractString, comm::MPI.Comm)
    adios::Union{Adios,Nothing}

Starting point for MPI apps. Creates an ADIOS handler. MPI collective
and it calls `MPI_Comm_dup`.
"""
function adios_init_mpi(comm::MPI.Comm)
    ptr = ccall((:adios2_init_mpi, libadios2_c_mpi), Ptr{Cvoid},
                (MPI.MPI_Comm,), comm)
    return ptr == C_NULL ? nothing : Adios(ptr)
end
function adios_init_mpi(config_file::AbstractString, comm::MPI.Comm)
    ptr = ccall((:adios2_init_config_mpi, libadios2_c_mpi), Ptr{Cvoid},
                (Cstring, MPI.MPI_Comm), config_file, comm)
    return ptr == C_NULL ? nothing : Adios(ptr)
end

export adios_init_serial
"""
    adios = adios_init_serial()
    adios = adios_init_serial(config_file::AbstractString)
    adios::Union{Adios,Nothing}

Initialize an Adios struct in a serial, non-MPI application. Doesn’t
require a runtime config file.

See also the [ADIOS2
documentation](https://adios2.readthedocs.io/en/latest/api_full/api_full.html#_CPPv418adios2_init_serialv).
"""
function adios_init_serial()
    ptr = ccall((:adios2_init_serial, libadios2_c), Ptr{Cvoid}, ())
    return ptr == C_NULL ? nothing : Adios(ptr)
end
function adios_init_serial(config_file::AbstractString)
    ptr = ccall((:adios2_init_config_serial, libadios2_c), Ptr{Cvoid},
                (Cstring,), config_file)
    return ptr == C_NULL ? nothing : Adios(ptr)
end

export declare_io
"""
    io = declare_io(adios::Adios, name::AbstractString)
    io::Union{AIO,Nothing}

Declare a new IO handler.

See also the [ADIOS2
documentation](https://adios2.readthedocs.io/en/latest/api_full/api_full.html#_CPPv417adios2_declare_ioP12adios2_adiosPKc).
"""
function declare_io(adios::Adios, name::AbstractString)
    ptr = ccall((:adios2_declare_io, libadios2_c), Ptr{Cvoid},
                (Ptr{Cvoid}, Cstring), adios.ptr, name)
    return ptr == C_NULL ? nothing : AIO(ptr, adios)
end

export adios_finalize
"""
    err = adios_finalize(adios::Adios)
    err::Error

Finalize the ADIOS context `adios`. It is usually not necessary to
call this function.

Instead of calling this function, one can also call the finalizer via
`finalize(adios)`. This finalizer is also called automatically when
the Adios object is garbage collected.

See also the [ADIOS2
documentation](https://adios2.readthedocs.io/en/latest/api_full/api_full.html#_CPPv415adios2_finalizeP12adios2_adios)
"""
function adios_finalize(adios::Adios)
    adios.ptr == C_NULL && return error_none
    err = ccall((:adios2_finalize, libadios2_c), Cint, (Ptr{Cvoid},), adios.ptr)
    adios.ptr = C_NULL
    return Error(err)
end


text:
using GLVisualize, GeometryTypes, GLAbstraction
using Colors, Reactive, FileIO
if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0, 1f0, 360))
end

description = """
Showing off the flexibility of the particle system by animating
all kind of atributes for an arbitrary mesh as the particle.
"""

cat    = GLNormalMesh(loadasset("cat.obj"))
sphere = GLNormalMesh(Sphere{Float32}(Vec3f0(0), 1f0), 12)

function scale_gen(v0, nv)
    l = length(v0)
    @inbounds for i=1:l
        v0[i] = Vec3f0(1,1,sin((nv*l)/i))/2
    end
    v0
end
function color_gen(v0, t)
    l = length(v0)
    @inbounds for x=1:l
        v0[x] = RGBA{N0f8}(x/l,(cos(t)+1)/2,(sin(x/l/3)+1)/2.,1.)
    end
    v0
end

t            = const_lift(x->x+0.1, timesignal)
ps           = sphere.vertices
scale_start  = Vec3f0[Vec3f0(1,1,rand()) for i=1:length(ps)]
scale        = foldp(scale_gen, scale_start, t)
colorstart   = color_gen(zeros(RGBA{N0f8}, length(ps)), value(t))
color        = foldp(color_gen, colorstart, t)
rotation     = sphere.normals
cats = visualize((cat, ps), scale = scale, color = color, rotation = rotation)

_view(cats, window)

if !isdefined(:runtests)
    renderloop(window)
end

ext: jl


text:
module SemidefiniteProgramming

import
    Base.print,
    Base.show,
    Base.getindex,
    Base.setindex!,
    Base.size,
    Base.start,
    Base.next,
    Base.done,
    Base.copy,
    Base.(*),
    Base.delete!,
    Base.(==),
    Base.isapprox

export
    SparseSDP,
    ismaximizationproblem,
    setmaximizationproblem!,
    setobj!,
    obj,
    setrhs!,
    rhs,
    setcon!,
    cons,
    ncons,
    blocksizes,

    writesdpasparse,
    readsdpasparse,

    SparseSDPSolution,
    dualobj,
    primalmatrix,
    dualvector,
    dualmatrix,

    SDPSolver,
    SDPA,
    SDPAQD,
    SDPAGMP,
    CSDP,
    solve

include("sparsesymmetricmatrix.jl")
include("sparsesymmetricblockmatrix.jl")
include("sparsesdp.jl")
include("normalization.jl")
include("sdpasparse.jl")
include("solvers.jl")
include("sdpsolution.jl")

end


text:
# Clayton, N. S., Yu, K. S. and Dickinson, A. (2003), Interacting Cache
# Memories: Evidence for Flexible Memory Use by Western Scrub-Jays (Aphelocoma
# californica) Journal of Experimental Psychology: Animal Behavior Processes,
# 1:14-22

function _summarize(exp::Experiment{:Clayton03_exp1}, results)
    res = similar(exp.data, 0)
    for g in groupby(@where(results, :action .== "cache", :set .== 7),
                     [:foodtype, :set])
        push!(res, [missing, g.foodtype[1], "cache", missing, missing, g.set[1],
                    mean(g.counts), sem(g.counts), missing, length(g.counts)])
    end
    res = [res; combine(groupby(@where(results, :action .== "inspect", :foodtype .== "other", :set .== 7),
                                [:group, :foodtype, :action, :RI, :trial, :set]),
             df -> DataFrame(μ = missing, sem = missing,
                             firstinspection = mean(df.firstinspection),
                             n = length(df.counts)))]
    res = [res; combine(groupby(cricketproportions(results), [:group, :RI]),
                    d -> DataFrame(foodtype = "cricketproportion", action = "inspect",
                                   set = 7, μ = mean(d.counts), sem = sem(d.counts),
                                   n = length(d.counts), trial = "pc",
                                   firstinspection = missing))]
    res
end
function statistical_tests(exp::Experiment{:Clayton03_exp1}, data)
    tests = Test[]
    datasetinspectClayton03_exp1 = @where(data, (|).(:set .== 4, :set .== 7),
                                         :action .== "inspect", :trial .== "pc",
                                         :RI .> 4)
    for g in groupby(datasetinspectClayton03_exp1, :RI)
        push!(tests, Test("first inspectes RI $(g.RI[1]) Clayton03 exp1",
                          fisherexacttest(g, cond1 = :group, cond2 = :foodtype,
                                          data = :firstinspection),
                          exp.tests))
    end
    cricketprop = cricketproportions(data)
    for g in groupby(cricketprop, :RI)
        id = "inspect RI $(g.RI[1]) Clayton03 exp1"
        push!(tests, Test(id,
                          anova(g, @formula(counts ~ group), exp.tests[id].locals),
                          exp.tests))
    end
#   How did they analyse this in the paper: RI 76 was measured twice which
#   yields unbalanced within factors.
    id = "cache Clayton03 exp1"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :action .== "cache", :set .== 7),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = [:foodtype, :RI],
                                     id = :id, data = :counts),
                      exp.tests))
    id = "inspect Clayton03 exp1"
    push!(tests, Test(id,
                      splitplotanova(cricketprop,
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = :RI, id = :id,
                                     data = :counts),
                      exp.tests))
#   TODO: "inspect degrade Clayton03 exp1" "inspect replenish Clayton03 exp1"
    TestCollection(tests)
end
function run!(::Experiment{:Clayton03_exp1}, models)
    results = DataFrame(group = Union{String,Missing}[],
                        foodtype = Union{String,Missing}[],
                        id = Int[], set = Int[], trial = Union{String,Missing}[],
                        RI = Union{Int,Missing}[],
                        action = String[], counts = Int[],
                        firstinspection = Union{Missing,Int}[])
    id = 0
    for group in ["degrade", "replenish"]
        for i in 1:8
            id += 1
            m = models[id]
            # training before Clayton 2003, experiment 1, January 2000
            wait!(m, 30u"d")
            for RI in (52u"hr", 76u"hr")
                Clayton01_trial!(m, results, group, id, 7, Peanut, Cricket, RI,
                                 degrade = false, pilfer = true, bothtrayspresent = true)
                #interleaved retraining
                Clayton01_trial!(m, results, group, id, 0, Peanut, Cricket, 4u"hr",
                                 degrade = false, pilfer = false)
            end
            wait!(m, 30u"d")
            #interleaved retraining
            Clayton01_trial!(m, results, group, id, 0, Peanut, Cricket, 100u"hr",
                             degrade = isdegrade(group, Cricket, 100u"hr"), pilfer = false)
            for RI in (76u"hr", 124u"hr")
                Clayton01_trial!(m, results, group, id, 7, Peanut, Cricket, RI,
                                 degrade = false, pilfer = true, bothtrayspresent = true)
                #interleaved retraining
                Clayton01_trial!(m, results, group, id, 0, Peanut, Cricket, 4u"hr",
                                 degrade = false, pilfer = false)
            end
        end
    end
    results
end

function _summarize(::Experiment{:Clayton03_exp2}, results)
    res = combine(groupby(@where(results, :action .== "inspect"),
                          [:group, :foodtype, :action, :condition, :RI]),
             d -> DataFrame(action = "inspect", μ = mean(d.counts),
                            sem = sem(d.counts),
                            firstinspection = mean(d.firstinspection),
                            n = length(d.counts)))
    res = [res; combine(groupby(@where(results, :action .== "cache"),
                                [:group, :foodtype, :action]),
                   d -> DataFrame(action = "cache", μ = mean(d.counts),
                                  sem = sem(d.counts),
                                  firstinspection = missing,
                                  condition = missing, RI = missing,
                                  n = length(d.counts)))]
end
function statistical_tests(exp::Experiment{:Clayton03_exp2}, data)
    tests = Test[]
    for g in groupby(@where(data, :action .== "inspect"), :group)
        group = g.group[1]
        id = "inspect $group"
        push!(tests, Test(id,
                          splitplotanova(g,
                                         exp.tests[id].locals,
                                         betweenfactors = :condition,
                                         withinfactors = [:foodtype, :RI],
                                         id = :id, data = :counts),
                         exp.tests))
        for subg in groupby(g, :RI)
            RI = subg.RI[1]
            id = "inspect $group RI $RI"
            push!(tests, Test(id,
                              splitplotanova(subg,
                                             exp.tests[id].locals,
                                             betweenfactors = :condition,
                                             withinfactors = :foodtype,
                                             id = :id, data = :counts),
                              exp.tests))
            if group == "degrade" && RI == 76
                for subsubg in groupby(subg, :condition)
                    id = "inspect $group RI $RI $(subsubg.condition[1])"
                    push!(tests, Test(id,
                                      anova(subsubg, @formula(counts ~ foodtype + id), exp.tests[id].locals),
                                      exp.tests))
                end
                for subsubg in groupby(subg, :foodtype)
                    id = "inspect $group RI $RI $(subsubg.foodtype[1])"
                    push!(tests, Test(id,
                                      anova(subsubg, @formula(counts ~ condition), exp.tests[id].locals),
                                      exp.tests))
                end
            end
        end
        if group == "replenish"
            for subg in groupby(g, :condition)
                id = "inspect $group $(subg.condition[1])"
                push!(tests, Test(id,
                                  splitplotanova(subg,
                                                 exp.tests[id].locals,
                                                 betweenfactors = [],
                                                 withinfactors = [:foodtype, :RI],
                                                 id = :id, data = :counts),
                                  exp.tests))
            end
            for subg in groupby(g, :foodtype)
                id = "inspect $group $(subg.foodtype[1])"
                push!(tests, Test(id,
                                  splitplotanova(subg,
                                                 exp.tests[id].locals,
                                                 betweenfactors = :condition,
                                                 withinfactors = [:RI],
                                                 id = :id, data = :counts),
                                  exp.tests))
            end
        end
        gsorted = sort(g, [:id, :RI, :foodtype])
        gsortedconsistent = @where(gsorted, :condition .== "consistent")
        gsortedreversed = @where(gsorted, :condition .== "reversed")
        push!(tests, Test("first inspectes $group",
                          chisqtest(firstinspect_dist(gsortedconsistent,
                                                      n = group == "degrade" ? 4 : 6),
                                    firstinspect_dist(gsortedreversed,
                                                      n = group == "degrade" ? 4 : 6)),
                          exp.tests))
    end
    for g in groupby(@where(data, :action .== "cache"), :group)
        id = "cache $(g.group[1])"
        push!(tests, Test(id,
                          splitplotanova(g,
                                         exp.tests[id].locals,
                                         betweenfactors = :condition,
                                         withinfactors = [:foodtype, :RI],
                                         id = :id, data = :counts),
                          exp.tests))
    end
    TestCollection(tests)
end
hack(i) = i
@inline function Clayton03_trial!(results, m, group, id, condition, RI)
    traypairs = [(Tray(hack(1)), Tray(hack(2))) for d in 1:3]
    # caching
    for (i, (tray1, tray2)) in enumerate(traypairs)
        remove!(m, Any)
        wait!(m, 4u"hr")
        cover!(tray2)
        add!(m, tray1); add!(m, tray2)
        add!(m, Peanut, 50)
        wait!(m, 15u"minute")
        remove!(m, Any)
        cover!(tray1); uncover!(tray2)
        add!(m, tray1); add!(m, tray2)
        add!(m, Cricket, 50)
        wait!(m, 15u"minute")
        remove!(m, Any)
        uncover!(tray1)
        if i == 3
            push!(results, [group, "peanut", id, condition, value(RI), "cache",
                            countcache(tray1), missing])
            push!(results, [group, "cricket", id, condition, value(RI), "cache",
                            countcache(tray2), missing])
        end
        add!(m, MaintenanceDiet)
        wait!(m, 25u"minute" + 15u"hr")
    end
    # recovery
    wait!(m, RI - 3u"d" - 4u"hr")
    for (i, (tray1, tray2)) in enumerate(traypairs)
        if i == 3
            pilfer!(tray1); pilfer!(tray2)
        elseif condition == "reversed"
            degrade!(tray2)
        end
        remove!(m, Any)
        wait!(m, 4u"hr")
        add!(m, tray1)
        add!(m, tray2)
        inspectionobserver = InspectionObserver()
        add!(m, inspectionobserver)
        wait!(m, 5u"minute")
        remove!(m, Any)
        if i == 3
            push!(results, [group, "peanut", id, condition, value(RI), "inspect",
                            countinspections(inspectionobserver, tray1),
                            firstinspection(inspectionobserver, tray1)])
            push!(results, [group, "cricket", id, condition, value(RI), "inspect",
                            countinspections(inspectionobserver, tray2),
                            firstinspection(inspectionobserver, tray2)])
        end
        add!(m, MaintenanceDiet)
        wait!(m, 55u"minute" + 15u"hr")
    end
end
function run!(::Experiment{:Clayton03_exp2}, models)
    results = DataFrame(group = Union{String,Missing}[],
                        foodtype = Union{String,Missing}[],
                        id = Int[], condition = Union{String,Missing}[],
                        RI = Union{Int,Missing}[],
                        action = String[], counts = Int[],
                        firstinspection = Union{Missing,Int}[])
    id = 0
    for group in ["degrade", "replenish"]
        for i in 1:8
            id += 1
            m = models[id]
            # retraining Feb 2000
            for wait in (4u"hr", 100u"hr")
                Clayton01_trial!(m, [], group, id, 0, Peanut, Cricket, wait,
                       degrade = isdegrade(group, Cricket, wait), pilfer = false)
            end
            condition = i > 4 ? "reversed" : "consistent"
            Clayton03_trial!(results, m, group, id, condition, 76u"hr")
            group == "replenish" && Clayton03_trial!(results, m, group, id, condition, 124u"hr")
            tmpres = [] # since Clayton01_trial! expects another form of results table
            Clayton01_trial!(m, tmpres, group, id, 1, Peanut, Cricket, 28u"hr",
                             degrade = false, pilfer = true, specificothername = true)
            for row in tmpres
                push!(results, [row[1:3]; condition; row[end-3:end]])
            end
        end
    end
    results
end


text:
initialize(self::ForcingBase, grid, state) = nothing

function initialize(self::ForcingBase{ForcingLES}, grid, state, LESDat::LESData)
    aux_gm = center_aux_grid_mean(state)
    nt = NC.Dataset(LESDat.les_filename, "r") do data
        imin = LESDat.imin
        imax = LESDat.imax

        zc_les = Array(get_nc_data(data, "zc"))

        dTdt_hadv = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dtdt_hadv", imin, imax))
        H_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "thetali_mean", imin, imax))
        dTdt_fluc = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dtdt_fluc", imin, imax))
        dqtdt_hadv = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dqtdt_hadv", imin, imax))
        qt_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "qt_mean", imin, imax))
        dqtdt_fluc = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dqtdt_fluc", imin, imax))
        subsidence = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "ls_subsidence", imin, imax))
        u_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "u_mean", imin, imax))
        v_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "v_mean", imin, imax))
        (; dTdt_hadv, H_nudge, dTdt_fluc, dqtdt_hadv, qt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge)
    end
    for k in real_center_indices(grid)
        aux_gm.dTdt_hadv[k] = nt.dTdt_hadv[k]
        aux_gm.H_nudge[k] = nt.H_nudge[k]
        aux_gm.dTdt_fluc[k] = nt.dTdt_fluc[k]
        aux_gm.dqtdt_hadv[k] = nt.dqtdt_hadv[k]
        aux_gm.qt_nudge[k] = nt.qt_nudge[k]
        aux_gm.dqtdt_fluc[k] = nt.dqtdt_fluc[k]
        aux_gm.subsidence[k] = nt.subsidence[k]
        aux_gm.u_nudge[k] = nt.u_nudge[k]
        aux_gm.v_nudge[k] = nt.v_nudge[k]
    end
end


text:
########################################################################################
# Caches
########################################################################################
abstract type ODEFiltersCache <: OrdinaryDiffEq.OrdinaryDiffEqCache end
mutable struct GaussianODEFilterCache{
    RType, ProjType, SolProjType, FP, uType, xType, AType, QType, matType, diffusionType, diffModelType,
    measType, llType,
} <: ODEFiltersCache
    # Constants
    d::Int                  # Dimension of the problem
    q::Int                  # Order of the prior
    A::AType
    Q::QType
    diffusionmodel::diffModelType
    R::RType
    Proj::ProjType
    SolProj::SolProjType
    Precond::FP
    # Mutable stuff
    u::uType
    u_pred::uType
    u_filt::uType
    tmp::uType
    x::xType
    x_pred::xType
    x_filt::xType
    x_tmp::xType
    x_tmp2::xType
    measurement::measType
    H::matType
    du::uType
    ddu::matType
    K::matType
    G::matType
    covmatcache::matType
    diffusion::diffusionType
    err_tmp::uType
    log_likelihood::llType
end

function OrdinaryDiffEq.alg_cache(
    alg::GaussianODEFilter, u, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck, IIP)
    initialize_derivatives=true

    if !(u isa AbstractVector)
        error("Problems which are not scalar- or vector-valued (e.g. u0 is a scalar or a matrix) are currently not supported")
    end

    if (alg isa EK1 || alg isa IEKS) && isnothing(f.jac)
        error("""EK1 requires the Jacobian. To automatically generate it with ModelingToolkit.jl use ODEFilters.remake_prob_with_jac(prob).""")
    end

    q = alg.order
    u0 = u
    t0 = t
    d = length(u)

    uType = typeof(u0)
    uElType = eltype(u0)
    matType = Matrix{uElType}

    # Projections
    Proj(deriv) = kron([i==(deriv+1) ? 1 : 0 for i in 1:q+1]', diagm(0 => ones(d)))
    SolProj = Proj(0)

    # Prior dynamics
    @assert alg.prior == :ibm "Only the ibm prior is implemented so far"
    Precond = preconditioner(d, q)
    A, Q = ibm(d, q, uElType)

    # Measurement model
    R = PSDMatrix(LowerTriangular(zeros(d, d)))
    # Initial states
    m0, P0 = initialize_derivatives ?
        initialize_with_derivatives(u0, f, p, t0, q) :
        initialize_without_derivatives(u0, f, p, t0, q)
    @assert iszero(P0)
    P0 = PSDMatrix(LowerTriangular(zero(P0)))
    x0 = Gaussian(m0, P0)

    # Pre-allocate a bunch of matrices
    h = Proj(0) * x0.μ
    H = copy(Proj(0))
    du = copy(u0)
    ddu = zeros(uElType, d, d)
    v, S = copy(h), copy(ddu)
    measurement = Gaussian(v, S)
    K = copy(H')
    G = copy(Matrix(P0))
    covmatcache = copy(G)

    diffusion_models = Dict(
        :dynamic => DynamicDiffusion(),
        :dynamicMV => MVDynamicDiffusion(),
        :fixed => FixedDiffusion(),
        :fixedMV => MVFixedDiffusion(),
        :fixedMAP => MAPFixedDiffusion(),
    )
    diffmodel = diffusion_models[alg.diffusionmodel]
    initdiff = initial_diffusion(diffmodel, d, q, uEltypeNoUnits)

    return GaussianODEFilterCache{
        typeof(R), typeof(Proj), typeof(SolProj), typeof(Precond),
        uType, typeof(x0), typeof(A), typeof(Q), matType, typeof(initdiff),
        typeof(diffmodel), typeof(measurement), uEltypeNoUnits,
    }(
        # Constants
        d, q, A, Q, diffmodel, R, Proj, SolProj, Precond,
        # Mutable stuff
        copy(u0), copy(u0), copy(u0), copy(u0),
        copy(x0), copy(x0), copy(x0), copy(x0), copy(x0),
        measurement,
        H, du, ddu, K, G, covmatcache, initdiff,
        copy(u0),
        zero(uEltypeNoUnits),
    )
end


text:
export adapt_rho_gamma_parallel

"""
Barzilai-Borwein scaling for Douglash-Rachford splitting on the dual problem related to standard ADMM.
Updates relaxation and augmented-Lagrangian penalty parameter.
To be used in parallel version of PARSDMM.
"""
function adapt_rho_gamma_parallel(
                                  gamma           ::Vector{TF},
                                  rho             ::Vector{TF},
                                  adjust_gamma    ::Bool,
                                  adjust_rho      ::Bool,
                                  y               ::Vector{Vector{TF}},
                                  y_old           ::Vector{Vector{TF}},
                                  s               ::Vector{Vector{TF}},
                                  s_0             ::Vector{Vector{TF}},
                                  l               ::Vector{Vector{TF}},
                                  l_hat_0         ::Vector{Vector{TF}},
                                  l_0             ::Vector{Vector{TF}},
                                  l_old           ::Vector{Vector{TF}},
                                  y_0             ::Vector{Vector{TF}},
                                  l_hat           ::Vector{Vector{TF}},
                                  d_l_hat         ::Vector{Vector{TF}},
                                  d_H_hat         ::Vector{Vector{TF}},
                                  d_l             ::Vector{Vector{TF}},
                                  d_G_hat         ::Vector{Vector{TF}}
                                  ) where {TF<:Real}

eps_correlation = TF(0.3) #hardcoded and suggested value by the paper based on numerical evidence

      if TF==Float64
          safeguard = 1e-10
      elseif TF==Float32
          safeguard = 1f-6
      end

      #standard loop-fusion
       @. l_hat[1]   = l_old[1] + rho[1]* ( -s[1] + y_old[1] )
      # @. d_l_hat[1] = l_hat[1] - l_hat_0[1]
      # @. d_H_hat[1] = s[1] - s_0[1]
      # @. d_l[1]     = l[1] - l_0[1]
      # @. d_G_hat[1] = -(y[1] - y_0[1])

      #wrap multi-theading inside a function
      a_is_b_min_c_MT!(d_l_hat[1],l_hat[1],l_hat_0[1])
      a_is_b_min_c_MT!(d_H_hat[1],s[1],s_0[1])
      a_is_b_min_c_MT!(d_l[1],l[1],l_0[1])
      a_is_b_min_c_MT!(d_G_hat[1],y_0[1],y[1])

      d_dHh_dlh   = dot(d_H_hat[1],d_l_hat[1])
      d_dGh_dl    = dot(d_G_hat[1],d_l[1])

      n_d_H_hat   = norm(d_H_hat[1])
      n_d_l_hat   = norm(d_l_hat[1])
      n_d_l       = norm(d_l[1])
      n_d_G_hat   = norm(d_G_hat[1])

      alpha_comp = false
      #print(alpha_comp)
      if (n_d_H_hat*n_d_l_hat) > safeguard && (n_d_H_hat.^2) > safeguard && d_dHh_dlh>safeguard
        alpha_correlation = d_dHh_dlh./( n_d_H_hat*n_d_l_hat  )
        if alpha_correlation > eps_correlation
          alpha_comp = true
          alpha_hat_MG  = d_dHh_dlh./(n_d_H_hat.^2);
          alpha_hat_SD  = (n_d_l_hat^2)./d_dHh_dlh;
          if (TF(2.0)*alpha_hat_MG) > alpha_hat_SD
            alpha_hat = alpha_hat_MG;
          else
            alpha_hat = alpha_hat_SD - alpha_hat_MG/TF(2.0);
          end
        end
      end

      beta_comp = false
      if (n_d_G_hat*n_d_l) > safeguard && (n_d_G_hat^2) > safeguard && d_dGh_dl > safeguard
        beta_correlation  = d_dGh_dl./( n_d_G_hat*n_d_l )
        if beta_correlation > eps_correlation
          beta_comp = true
          beta_hat_MG = d_dGh_dl./(n_d_G_hat^2);
          beta_hat_SD = (n_d_l^2)./d_dGh_dl;
          if (TF(2.0)*beta_hat_MG) > beta_hat_SD
            beta_hat = beta_hat_MG;
          else
            beta_hat = beta_hat_SD - beta_hat_MG/TF(2.0);
          end
        end
      end

      #update rho and or gamma
    if adjust_rho == true  && adjust_gamma == true
      if alpha_comp == true && beta_comp == true
        rho[1] = sqrt(alpha_hat*beta_hat);
        gamma[1]=TF(1.0)+(( TF(2.0)*sqrt(alpha_hat*beta_hat) )./( alpha_hat+beta_hat ));
      elseif alpha_comp == true && beta_comp == false
        rho[1] = alpha_hat;
        gamma[1]=TF(1.9);
      elseif  alpha_comp == false && beta_comp == true
        rho[1] = beta_hat;
        gamma[1]=TF(1.1);
      else
        #rho = rho; #do nothing
        gamma[1]=TF(1.5);
      end
    elseif adjust_rho == true  && adjust_gamma == false
        if alpha_comp == true && beta_comp == true
          rho[1] = sqrt(alpha_hat*beta_hat);
        elseif alpha_comp == true && beta_comp == false
          rho[1] = alpha_hat;
        elseif alpha_comp == false && beta_comp == true
          rho[1] = beta_hat;
        else
          #rho = rho; #do nothing
        end
      elseif adjust_rho == false  && adjust_gamma == true
        if alpha_comp == true && beta_comp == true
          gamma[1]=TF(1.0)+(( TF(2.0)*sqrt(alpha_hat*beta_hat) )./( alpha_hat+beta_hat ));
        elseif alpha_comp == true && beta_comp == false
          gamma[1]=TF(1.9);
        elseif  alpha_comp == false && beta_comp == true
          gamma[1]=TF(1.1);
        else
          gamma[1]=TF(1.5);
        end
      end #end compute new rho and gamma

# copy!(l_hat_0[1],l_hat[1])
# copy!(y_0[1],y[1])
# copy!(s_0[1],s[1])
# copy!(l_0[1],l[1])

#return rho,gamma,l_hat,d_l_hat,d_H_hat,d_l,d_G_hat
end #function adapt_rho_gamma


text:
@testset "misc" begin
    @test J.val(nothing) == 1
    @test J.val(2) == 1
    @test J.Option{Int} == Union{Nothing,Int}
end

@testset "isdef" begin

    b = board(name="b")
    b ++ @jsf foo(x) = 5x
    b ++ @jsf bar(x) = 7x
    b ++ @jsf foo(x) = 2x

    @jsf foo(x) = 2x
    @jsf bar(x) = 7x

    @test J.isdef(foo, b)
    @test J.isdef(bar, b)

    @jsf foo(x) = 5x

    @test !J.isdef(foo, b)
    @test !J.isdef(bar, nothing)
end


text:
const VERSION_DATASET_CHECKSUM = Dict(
    v"1.0.0" => DATASET_CHECKSUM_V1_0_0,
    v"1.1.0" => DATASET_CHECKSUM_V1_1_0,
)
const VERSION_DATASET_DESCRIPTION = Dict(
    v"1.0.0" => DATASET_DESCRIPTION_V1_0_0,
    v"1.1.0" => DATASET_DESCRIPTION_V1_1_0,
)
const VERSION_DESCRIPTION = Dict(
    v"1.0.0" => VERSION_DESCRIPTION_V1_0_0,
    v"1.1.0" => VERSION_DESCRIPTION_V1_1_0,
)
const VERSION_HASH = Dict(
    v"1.0.0" => "ff3b621f-6d31-40a8-b578-e87d7041f27a",
    v"1.1.0" => "d03ca12b-ac05-4fd9-8cd4-0be805c1fc99",
)


text:
module Keras

using Compat
using PyCall

using Compat: @__MODULE__
import PyCall: PyObject, @pydef

# We store our renamed aliases to aid with
# automatic docstring updating.
global const keras_func_aliases = Dict(
    "compile" => "compile!",
    "add" => "add!"
)

"""
Copied from LazyHelp in PyPlot.jl.

define a documentation object that lazily looks up
help from a PyObject via zero or more keys.
This saves us time when loading Keras, since we don't have
to load up all of the documentation strings right away.
"""
struct PyDoc
    obj::PyObject
    names::Tuple{Vararg{Symbol}}

    PyDoc(obj::PyObject, name::Symbol) = new(obj, (name,))
end

function Base.show(io::IO, ::MIME"text/plain", doc::PyDoc)
    obj = doc.obj

    for name in doc.names
        obj = obj[name]
    end

    if haskey(obj, "__doc__")
        print(io, convert(AbstractString, obj[:__doc__]))
    else
        print(io, "no Python docstring found for PyDoc($(doc.obj), $(doc.names))")
    end
end

Base.show(io::IO, doc::PyDoc) = show(io, "text/plain", doc)

function Base.Docs.catdoc(docs::PyDoc...)
    Base.Docs.Text() do io
        for doc in docs
            show(io, MIME"text/plain"(), doc)
        end
    end
end

# We need to handle our python dependencies carefully here
global const _keras = PyNULL()
global const _backend = PyNULL()
global const _layers = PyNULL()
global const _models = PyNULL()
global const _regularizers = PyNULL()
global const _optimizers = PyNULL()
global const _callbacks = PyNULL()
global const _constraints = PyNULL()
global const _initializers = PyNULL()

function __init__()
    copy!(_keras, pyimport("keras"))
    copy!(_backend, pyimport("keras.backend"))
    copy!(_layers, pyimport("keras.layers"))
    copy!(_models, pyimport("keras.models"))
    copy!(_regularizers, pyimport("keras.regularizers"))
    copy!(_optimizers, pyimport("keras.optimizers"))
    copy!(_callbacks, pyimport("keras.callbacks"))
    copy!(_constraints, pyimport("keras.constraints"))
    copy!(_initializers, pyimport("keras.initializers"))

    include(joinpath(@__DIR__, "utils.jl"))
end

include("tensors.jl")
include("callbacks.jl")
include("constraints.jl")
include("initializers.jl")
include("optimizers.jl")
include("regularizers.jl")
include("layers.jl")
include("models.jl")

end


text:
using Test

function quick_sort(arr)
    if length(arr) < 2 return arr end
    pivot = arr[1]
    less = [i for i in arr[2:end] if i <= pivot]
    greater = [i for i in arr[2:end] if i > pivot]
    return vcat(quick_sort(less), [pivot], quick_sort(greater))
end

@test quick_sort([3,5,2,1,4]) == [1,2,3,4,5]


text:
module FieldValues

using TensorValues

export FieldValue
export Point
export VectorValue
export TensorValue
export MultiValue
export inner
export outer
export meas
export normalvec
export trace
export tr
export symmetic_part
export mutable

import TensorValues: meas

"""
Type representing all possible field value types
"""
const FieldValue = Number

"""
Type representing a point of D dimensions with coordinates of type T
"""
const Point{D,T} = VectorValue{D,T}

function normalvec(v::MultiValue{Tuple{1,2}})
  n1 = v[1,2]
  n2 = -1*v[1,1]
  VectorValue(n1,n2)
end

function normalvec(v::MultiValue{Tuple{2,3}})
  n1 = v[1,2]*v[2,3] - v[1,3]*v[2,2]
  n2 = v[1,3]*v[2,1] - v[1,1]*v[2,3]
  n3 = v[1,1]*v[2,2] - v[1,2]*v[2,1]
  VectorValue(n1,n2,n3)
end

function meas(v::MultiValue{Tuple{1,2}})
  n = normalvec(v)
  sqrt(n*n)
end

function meas(v::MultiValue{Tuple{2,3}})
  n = normalvec(v)
  sqrt(n*n)
end

end # module


text:
using Base.Test
import Ogre

module test_util_resources
    immutable my_immutable_type
        field1
        field2
    end

    type my_mutable_type
        field1
        field2
    end
end

@testset "Utility function tests" begin
    res = test_util_resources

    @testset "Copy-modification of types" begin
        t1 = res.my_immutable_type("hello", 12)
        t2 = Ogre.cpmod(t1, field1="world")
        @test t2.field1 == "world"
        @test t2.field2 == 12

        t3 = res.my_mutable_type("foo", 13)
        t4 = Ogre.cpmod(t3, field2=14)
        @test t4.field1 == "foo"
        @test t4.field2 == 14
    end

    @testset "Mapping rows of a function" begin
        M = [1 2 3;
             4 5 6;
             7 8 9]
        @test vec(Ogre.maprows(sum, M)) == [6, 15, 24]
        @test vec(Ogre.maprows(mean, M)) == [2., 5., 8.]
    end

    @testset "Testing for NaN values" begin
        @test !Ogre.hasnan([1, 2, 3, 4])
        @test Ogre.hasnan([1, NaN, 3, NaN])
    end
end


text:
# reduced kernel for the s channel
function compute_s_reduced!(
    Λ    :: Float64,
    buff :: Matrix{Float64},
    v    :: Float64,
    dv   :: Float64,
    s    :: Float64,
    vsp  :: Float64,
    r    :: Reduced_lattice,
    m    :: Mesh,
    a    :: Action_su2,
    temp :: Array{Float64, 3}
    )    :: Nothing

    # get propagator
    p = -get_propagator(Λ, v + 0.5 * s, 0.5 * s - v, m, a)

    # get buffers for right vertex (left vertex is given by bare)
    bs2 = get_buffer_s(s, v, vsp, m)
    bt2 = get_buffer_t(-v - vsp, 0.5 * (s + v - vsp), 0.5 * (s - v + vsp), m)
    bu2 = get_buffer_u( v - vsp, 0.5 * (s + v + vsp), 0.5 * (s - v - vsp), m)

    # cache vertex values for all lattice sites in temporary buffer
    get_Γ_avx!(r, bs2, bt2, bu2, a, temp, 2)

    # compute contributions for all lattice sites
    @turbo unroll = 1 for i in eachindex(r.sites)
        # read cached values for site i
        v1s = a.Γ[1].bare[i]; v1d = a.Γ[2].bare[i]
        v2s =  temp[i, 1, 2]; v2d = temp[i, 2, 2]

        # compute contribution at site i
        Γs = -p * (-2.0 * v1s * v2s + v1s * v2d + v1d * v2s)
        Γd = -p * ( 3.0 * v1s * v2s + v1d * v2d)

        # parse result to output buffer
        buff[1, i] += dv * Γs
        buff[2, i] += dv * Γd
    end

    return nothing
end

# integration kernel for loop function
function compute_Σ_kernel(
    Λ  :: Float64,
    w  :: Float64,
    v  :: Float64,
    r  :: Reduced_lattice,
    m  :: Mesh,
    a1 :: Action_su2,
    a2 :: Action_su2
    )  :: Float64

    # get buffers for non-local vertex
    b1s = get_buffer_s(v + w, Inf, 0.5 * (v - w), m)
    b1t = get_buffer_empty()
    b1u = get_buffer_empty()

    # get buffers for local vertex
    b2s = get_buffer_s(v + w, Inf, 0.5 * (-v + w), m)
    b2t = get_buffer_empty()
    b2u = get_buffer_empty()

    # compute local contributions
    val = 3.0 * get_Γ_comp(1, 1, b2s, b2t, b2u, r, a2, apply_flags_su2, ch_t = false, ch_u = false) +
                get_Γ_comp(2, 1, b2s, b2t, b2u, r, a2, apply_flags_su2, ch_t = false, ch_u = false)

    # compute contributions for all lattice sites
    for j in eachindex(r.sites)
        val -= 2.0 * r.mult[j] * (2.0 * a2.S) * get_Γ_comp(2, j, b1s, b1t, b1u, r, a2, apply_flags_su2, ch_t = false, ch_u = false)
    end

    # multiply with full propagator
    val *= -get_G(Λ, v, m, a1) / (2.0 * pi)

    return val
end


text:
module CVImgCodecs

export imread, imwrite

using LibOpenCV
using CVCore
using Cxx

libopencv_imgcodecs = LibOpenCV.find_library_e("libopencv_imgcodecs")
try
    Libdl.dlopen(libopencv_imgcodecs, Libdl.RTLD_GLOBAL)
catch e
    warn("You might need to set DYLD_LIBRARY_PATH to load dependencies proeprty.")
    rethrow(e)
end

cxx"""
#include <opencv2/imgcodecs.hpp>
"""

for name in [
    # Imread flags
    :IMREAD_UNCHANGED,
    :IMREAD_GRAYSCALE,
    :IMREAD_COLOR,
    :IMREAD_ANYCOLOR,
    :IMREAD_LOAD_GDAL,
    :IMREAD_REDUCED_GRAYSCALE_2,
    :IMREAD_REDUCED_COLOR_2,
    :IMREAD_REDUCED_GRAYSCALE_4,
    :IMREAD_REDUCED_COLOR_4,
    :IMREAD_REDUCED_GRAYSCALE_8,
    :IMREAD_REDUCED_COLOR_8,

    # Imwrite flags
    :IMWRITE_JPEG_QUALITY,
    :IMWRITE_JPEG_PROGRESSIVE,
    :IMWRITE_JPEG_OPTIMIZE,
    :IMWRITE_JPEG_RST_INTERVAL,
    :IMWRITE_JPEG_LUMA_QUALITY,
    :IMWRITE_JPEG_CHROMA_QUALITY,
    :IMWRITE_PNG_COMPRESSION,
    :IMWRITE_PNG_STRATEGY,
    :IMWRITE_PNG_BILEVEL,
    :IMWRITE_PXM_BINARY,
    :IMWRITE_WEBP_QUALITY,

    # ImwritePNGFlags
    :IMWRITE_PNG_STRATEGY_DEFAULT,
    :IMWRITE_PNG_STRATEGY_FILTERED,
    :IMWRITE_PNG_STRATEGY_HUFFMAN_ONLY,
    :IMWRITE_PNG_STRATEGY_RLE,
    :IMWRITE_PNG_STRATEGY_FIXED,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), string("cv::", name, ";"))
    @eval global const $name = $ex
end

function imread(name::AbstractString, flag=IMREAD_UNCHANGED)
    Mat(@cxx cv::imread(pointer(name), flag))
end

imwrite(filename::AbstractString, img::AbstractCvMat) =
    icxx"cv::imwrite($(pointer(filename)), $(img.handle));"
imwrite(filename::AbstractString, img::AbstractCvMat, params) =
    icxx"cv::imwrite($(pointer(filename)), $(img.handle), $params);"

end # module


text:
# Create a model with const parameters ****************
#
export Model2Fit_Const_Parameters
export get_model, get_model_θ

@doc raw"""

Set some parameters to constant values

# Example

Let's assume that we have an initial model `model` with `θ` as parameter vector.
If we want θ1, θ3, θ5 to be constant, proceed as follows

```julia
indices = [1,3,5]
values = Float64[1,2,3]

model_with_const_params = Model2Fit_Const_Parameters(model, indices, values)
```

Do not forget to udpate parameter vector θ, this can be done as
follows:

```julia
deleteat!(θ, indices)
```

# Also see

"""
struct Model2Fit_Const_Parameters{MODEL <: Abstract_Model2Fit,
                                  INDICES <: AbstractVector{Int},
                                  VALUES <: AbstractVector}  <: Abstract_Model2Fit
    _model::MODEL
    _indices::INDICES
    _values::VALUES

    function Model2Fit_Const_Parameters(model::MODEL,
                                        indices::INDICES,
                                        values::VALUES) where {MODEL <: Abstract_Model2Fit,
                                                               INDICES <: AbstractVector{Int},
                                                               VALUES <: AbstractVector}
        @assert (length(indices) == 0) || ( (1 ≤ minimum(indices)) && (maximum(indices) ≤ parameter_size(model)) )
        @assert allunique(indices)
        @assert length(indices) == length(values)

        new{MODEL,INDICES,VALUES}(model, indices, values)

    end

end

# Convenience methods ================
#
@doc raw"""
```julia
get_model(model::Model2Fit_Const_Parameters) -> Absatrct_Model2Fit
```

Get back the wrapped model
"""
get_model(model::Model2Fit_Const_Parameters) = model._model

@doc raw"""
```julia
get_model_θ(mp::Model2Fit_Const_Parameters,θ::AbstractVector) -> θ::AbstractVector
```

Retrieve the parameter vector θ associated to the wrapped model [`get_model`](@ref).
"""
function get_model_θ(model::Model2Fit_Const_Parameters,θ::AbstractVector)
    @assert length(θ) == parameter_size(model)

    # Create the new θ by inserting constant values
    #
    insert_some_elements(θ, model._indices, model._values)
end

# Visit  ================
#
# note: no need to check
#   @assert 1 ≤ submodel_idx ≤ visit_submodel_size(model)
#   @assert length(X) == length(Y)
# as this is done higher in the visit() method
#
visit_submodel_size(model::Model2Fit_Const_Parameters) = 1
visit_get_submodel(model::Model2Fit_Const_Parameters,submodel_idx::Int) = model._model
visit_get_X(model::Model2Fit_Const_Parameters,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = X
visit_get_Y(model::Model2Fit_Const_Parameters,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = Y
visit_get_θ(model::Model2Fit_Const_Parameters,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = get_model_θ(model,θ)

# Interface ================
#
parameter_size(mp::Model2Fit_Const_Parameters) = parameter_size(mp._model) - length(mp._indices)

function accumulate_y!(model::Model2Fit_Const_Parameters,Y::AbstractVector,X::AbstractVector,θ::AbstractVector)
    @assert length(θ) == parameter_size(model)

    inner_model = get_model(model)
    θ_inner_model = get_model_θ(model,θ)

    accumulate_y!(inner_model,Y,X,θ_inner_model)
end


text:
module SpikeNetOpt

using PyCall
using OrderedCollections
using LinearAlgebra
using SpikeSynchrony
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
using Evolutionary
using DataStructures
using UnicodePlots
using Statistics
using JLD
using Plots
using Distributions
using LightGraphs
using Metaheuristics
#using SignalAnalysis
# using Distributed
include("current_search.jl")
include("utils.jl")
include("sdo_network.jl")

end


text:
module BEASTXMLConstructor

export BEASTXMLElement,
       MCMCOptions,
       save_xml,
       make_residual_xml,
       make_pfa_xml,
       make_joint_xml,
       set_mbd_precision,
       set_residual_precision,
       set_loadings,
       set_factor_precisions,
       DiffusionModel,
       ResidualVarianceModel,
       IntegratedFactorModel,
       JointProcessModel,
       DataModel,
       tip_dimension,
       data_dimension,
       set_options!,
       set_chain_length!,
       get_operators,
       set_file_logEvery,
       add_factor_proportion

using LightXML, LinearAlgebra, UnPack, StatsBase

include("BeastNames.jl")
bn = BeastNames

XMLOrNothing = Union{XMLElement, Nothing}

abstract type MyXMLElement end
abstract type SimpleXMLElement <: MyXMLElement end
abstract type AbstractDataXMLElement <: MyXMLElement end

import Base: copy

import LightXML: name
function name(x::MyXMLElement)
    error("Not yet implemented for object of type " * string(typeof(x)))
end

mutable struct BEASTXMLElement
    el::XMLOrNothing
    components::Vector{MyXMLElement}
end

mutable struct MCMCOptions
    chain_length::Int
    file_log_every::Int
    screen_log_every::Int
    likelihood_check_count::Int
end

function MCMCOptions(;
                     chain_length::Int = 10_000,
                     file_log_every::Int = max(1, div(chain_length, 1_000)),
                     screen_log_every::Int = max(1, div(chain_length, 100)),
                     likelihood_check_count::Int = -1)

    if likelihood_check_count == -1
        check_count = div(chain_length, 10)
        check_count = max(check_count, 100)
        check_count = min(check_count, 1000)
    else
        check_count = likelihood_check_count
    end

    return MCMCOptions(chain_length, file_log_every, screen_log_every,
                       likelihood_check_count)
end

import Base.show

function Base.show(io::IO, bx::BEASTXMLElement)
    disp = "$(typeof(bx)) with components:"
    for comp in bx.components
        disp *= "\n\t$(typeof(comp))"
    end
    println(io, disp)
end

function BEASTXMLElement()
    return BEASTXMLElement(nothing, MyXMLElement[])
end

abstract type ModelExtensionXMLElement <: MyXMLElement end
abstract type OperatorXMLElement <: MyXMLElement end

const element_dir = joinpath(@__DIR__, "XMLConstructor", "XMLElements")
const dir = joinpath(@__DIR__, "XMLConstructor")

include(joinpath(dir, "constructor_helpers.jl"))
include(joinpath(dir, "integration.jl"))

include(joinpath(element_dir, "NothingXMLElement.jl"))
include(joinpath(element_dir, "CompoundXMLElement.jl"))
include(joinpath(element_dir, "MatrixParameter.jl"))
include(joinpath(element_dir, "DataXMLElement.jl"))
include(joinpath(element_dir, "NewickXMLElement.jl"))
include(joinpath(element_dir, "TreeModelXMLElement.jl"))
include(joinpath(element_dir, "WishartPriorXMLElement.jl"))
include(joinpath(element_dir, "MBDXMLElement.jl"))
include(joinpath(element_dir, "RepeatedMeasuresXMLElement.jl"))
include(joinpath(element_dir, "MatrixShrinkage.jl"))
include(joinpath(element_dir, "IntegratedFactorsXMLElement.jl"))
include(joinpath(element_dir, "TraitLikelihoodXMLElement.jl"))
include(joinpath(element_dir, "LatentFactorModelXMLElement.jl"))
include(joinpath(element_dir, "LogPredictiveDensity.jl"))
include(joinpath(element_dir, "NormalGammaPrecisionOperatorXMLElement.jl"))
include(joinpath(element_dir, "PrecisionGibbsOperatorXMLElement.jl"))
include(joinpath(element_dir, "OperatorsXMLElement.jl"))
include(joinpath(element_dir, "LoggablesXMLElement.jl"))
include(joinpath(element_dir, "MCMCXMLElement.jl"))
include(joinpath(element_dir, "TimerXMLElement.jl"))
include(joinpath(element_dir, "TraitValidationXMLElement.jl"))
include(joinpath(element_dir, "OldLoadingsGibbsOperatorXMLElement.jl"))
include(joinpath(element_dir, "HMCOperatorXMLElement.jl"))
include(joinpath(element_dir, "LoadingsGibbsOperatorXMLElement.jl"))
include(joinpath(element_dir, "LatentFactorPrecisionOperatorXMLElement.jl"))
include(joinpath(element_dir, "FactorTreeGibbsOperatorXMLElement.jl"))
include(joinpath(element_dir, "TraitLoggerXMLElement.jl"))
include(joinpath(element_dir, "LKJPrecisionXMLElement.jl"))
include(joinpath(element_dir, "NormalMatrixNormLikelihood.jl"))
include(joinpath(element_dir, "ScaledOrthogonalMatrix.jl"))
include(joinpath(element_dir, "MultivariateGammaLikelihood.jl"))
include(joinpath(element_dir, "NormalDistributionLikelihood.jl"))
include(joinpath(element_dir, "MultiplicativeScalePrior.jl"))
include(joinpath(element_dir, "ForceOrderedLikelihood.jl"))
include(joinpath(element_dir, "IntegratedFactorsGibbsOperator.jl"))
include(joinpath(element_dir, "IndependentNormalDistributionModel.jl"))
include(joinpath(element_dir, "LoadingsScaleGibbsOperator.jl"))
include(joinpath(element_dir, "LatentLiability.jl"))
include(joinpath(element_dir, "FactorProportionStatistic.jl"))

include(joinpath(dir, "utils.jl"))
include(joinpath(dir, "constructors.jl"))
include(joinpath(dir, "getters.jl"))
include(joinpath(dir, "setters.jl"))

function make_xml(bx::BEASTXMLElement)
    xdoc = XMLDocument()
    bx.el = new_element(bn.BEAST)
    set_root(xdoc, bx.el)

    for el in bx.components
    add_el(bx, el)
    end

    return xdoc
end

function save_xml(path::String, bx::BEASTXMLElement;
                  change_filename::Bool = true)

    nm = basename(path)
    s = split(nm, '.')

    if length(s) == 1
        path = path * ".xml"
    elseif length(s) == 2
        @assert s[2] == "xml"
    end

    if change_filename
        filename = s[1]
        set_filename(bx, filename)
    end

    xdoc = make_xml(bx)
    save_file(xdoc, path)
    free(xdoc)
    return nothing
end

import LightXML: add_child
function add_child(xel::XMLElement, no::Nothing)
    # do nothing
end

function add_child(bx::BEASTXMLElement, el::MyXMLElement)
    push!(bx.components, el)
end

function add_child(bx::BEASTXMLElement, el::MyXMLElement, ind::Int)
    insert!(bx.components, ind, el)
end

import LightXML: find_element

function find_element(bx::BEASTXMLElement, type::Type)
    ind = find_element_ind(bx, type)
    if isnothing(ind)
        return nothing
    end
    return bx.components[find_element_ind(bx, type)]
end

function find_element_ind(bx::BEASTXMLElement, type::Type)
    if !(type <: MyXMLElement)
        error("Can only find objects of type MyXMLElement.")
    end

    found = findall(x -> typeof(x) <: type, bx.components)

    if length(found) == 0
        return nothing
    elseif length(found) > 1
        error("Multiple elements of type $type. " *
              "Please use function 'find_elements'.")
    end

    return found[1]
end

function find_elements(bx::BEASTXMLElement, type::Type)
    found = findall(x -> typeof(x) <: type, bx.components)
    return bx.components[found]
end

function add_parameter_id(pel::XMLElement, id::String)
    el = new_element(bn.PARAMETER)
    set_attribute(el, bn.ID, id)
    add_child(pel, el)
    return el
end

function add_loggable(bx::BEASTXMLElement, el::MyXMLElement;
                        already_made::Bool = false)

    loggables = LoggablesXMLElement([el], [already_made])
    add_loggables(bx, loggables)
end

function add_loggables(bx::BEASTXMLElement, loggables::LoggablesXMLElement)
    old_loggables = find_elements(bx, LoggablesXMLElement)
    if length(old_loggables) == 0
        mcmc = get_mcmc(bx)
        ind = findfirst(x -> x === mcmc, bx.components)
        mcmc_loggables = mcmc.loggables
        push!(mcmc_loggables, loggables)
        insert!(bx.components, ind, mcmc_loggables)
    elseif length(old_loggables) == 1
        push!(old_loggables[1], loggables)
    else
        error("Cannot join with more than one LoggablesXMLElement")
    end
    return nothing
end

function make_xml(n::Nothing)
    # do nothing
end

function add_el(bx::BEASTXMLElement, no::Nothing)
    #do nothing
end

function add_el(bx::BEASTXMLElement, el::TextXMLElement)
    add_child(bx.el, el.el)
end

function add_el(bx::BEASTXMLElement, el::MyXMLElement)
    make_xml(el)
    add_child(bx.el, el.el)
end

function add_el(bx::BEASTXMLElement, mbd_el::MBDXMLElement)
    make_xml(mbd_el)
    add_child(bx.el, mbd_el.el)
    add_el(bx, mbd_el.precision_prior)
end

function add_el(bx::BEASTXMLElement, n_el::NothingXMLElement)
    # do nothing
end

function add_el(bx::BEASTXMLElement, cp_el::CompoundXMLElement)
    for el in cp_el.els
        add_el(bx, el)
    end
end

function add_el(bx::BEASTXMLElement, rm_el::RepeatedMeasuresXMLElement)
    make_xml(rm_el)
    add_child(bx.el, rm_el.el)
    add_child(bx.el, rm_el.precision_prior.el)
end

function add_el(bx::BEASTXMLElement, if_el::IntegratedFactorsXMLElement)
    make_xml(if_el)
    add_el(bx, if_el.loadings)
    add_el(bx, if_el.loadings_prior)
    add_child(bx.el, if_el.el)
    add_child(bx.el, if_el.precision_prior_el)
end

function add_el(bx::BEASTXMLElement, lfm_el::LatentFactorModelXMLElement)
    make_xml(lfm_el)
    if !lfm_el.parameters_already_made
        add_el(bx, lfm_el.loadings)
        add_el(bx, lfm_el.loadings_prior)
    end
    add_child(bx.el, lfm_el.el)
    if !lfm_el.parameters_already_made
        add_child(bx.el, lfm_el.precision_prior_el)
    end
end

function add_el(bx::BEASTXMLElement, el::FactorLogPredictiveDensity)
    make_xml(el)
    add_child(bx.el, el.fac_el)
    add_child(bx.el, el.like_el)
    add_child(bx.el, el.cl_el)
end

function add_el(bx::BEASTXMLElement, els::Array{MyXMLElement})
    for el in els
        add_el(bx, el)
    end
end

function add_el(bx::BEASTXMLElement, loggables::LoggablesXMLElement)
    for i = 1:length(loggables.els)
        if !loggables.already_made[i]
            l = loggables.els[i]
            add_el(bx, l)
        end
    end
end

function add_el(bx::BEASTXMLElement, el::ScaledOrthogonalMatrix)
    make_xml(el)
    add_child(bx.el, el.scale.el)
    add_child(bx.el, el.U.el)
    add_child(bx.el, el.el)
end

function add_el(bx::BEASTXMLElement, el::MatrixShrinkageLikelihoods)
    make_xml(el)
    for mult in el.mult_els
        add_child(bx.el, mult)
    end
    for gp in el.gp_els
        add_child(bx.el, gp)
    end
    for ls in el.local_scale_els
        add_child(bx.el, ls)
    end
    for gs in el.global_scale_els
        add_child(bx.el, gs)
    end
    for bb in el.bb_els
        add_child(bx.el, bb)
    end
    add_child(bx.el, el.ms_el)
    for gp in el.global_prior_els
        add_child(bx.el, gp)
    end
end

function add_el(bx::BEASTXMLElement, el::MultiplicativeScalePrior)
    add_el(bx, el.mults)
    add_el(bx, el.precs)
    add_el(bx, el.mult_prior)
    add_el(bx, el.scale_prior)
end

end


text:
using DataFrames
using SemiparametricAUC

# Data analysis examples
fasd = DataFrames.readtable(joinpath(Pkg.dir("SemiparametricAUC"), "data/fasd.csv"))
# fasd = readtable("ds.csv")

# Define factor/categorical variable
function convert_to_factor(x)
    return(DataFrames.pool(x))
end

fasd[:group] = convert_to_factor(fasd[:group])
fasd[:x1]    = convert_to_factor(fasd[:x1])
fasd[:x2]    = convert_to_factor(fasd[:x2])
# fasd[:x3]  = convert_to_factor(fasd[:x3])

# Perform analysis
one_covariates_results = SemiparametricAUC.semiparametricAUC(model_formula = y ~ x1,
                                              treatment_group = :group, data = fasd)

two_covariates_results = SemiparametricAUC.semiparametricAUC(model_formula = y ~ x1 + x2,
                                            treatment_group = :group, data = fasd)


text:
# generally, Documenter.jl assumes markdown files ends with `.md`
const markdown_exts = [".md",]

"""
    struct MarkdownDemoCard <: AbstractDemoCard
    MarkdownDemoCard(path::String)

Constructs a markdown-format demo card from existing markdown file `path`.

# Fields

Besides `path`, this struct has some other fields:

* `path`: path to the source markdown file
* `cover`: path to the cover image
* `id`: cross-reference id
* `title`: one-line description of the demo card
* `description`: multi-line description of the demo card

# Configuration

You can pass additional information by adding a YAML front matter to the markdown file.
Supported items are:

* `cover`: relative path to the cover image. If not specified, it will use the first available image link, or all-white image if there's no image links.
* `description`: a multi-line description to this file, will be displayed when the demo card is hovered. By default it uses `title`.
* `id`: specify the `id` tag for cross-references. By default it's infered from the filename, e.g., `simple_demo` from `simple demo.md`.
* `title`: one-line description to this file, will be displayed under the cover image. By default, it's the name of the file (without extension).

An example of the front matter:

```text
---
title: passing extra information
cover: cover.png
id: non_ambiguious_id
description: this demo shows how you can pass extra demo information to DemoCards package.
---
```

See also: [`DemoSection`](@ref DemoCards.DemoSection), [`DemoPage`](@ref DemoCards.DemoPage)
"""
struct MarkdownDemoCard <: AbstractDemoCard
    path::String
    cover::Union{String, Nothing}
    id::String
    title::String
    description::String
end

function MarkdownDemoCard(path::String)::MarkdownDemoCard
    # first consturct an incomplete democard, and then load the config
    card = MarkdownDemoCard(path, "", "", "", "")

    cover = load_config(card, "cover")
    id    = load_config(card, "id")
    title = load_config(card, "title")
    card = MarkdownDemoCard(path, cover, id, title, "")

    # default description requires a title
    description = load_config(card, "description")
    MarkdownDemoCard(path, cover, id, title, description)
end

"""
    parse(card::MarkdownDemoCard)

parse out configuration of markdown files and return it as a `Dict`.

Possible configuration resources are:

* YAML front matter
* image links

!!! note

    Users of this function need to use `haskey` to check if keys are existed.
    They also need to validate the values.
"""
function parse(card::MarkdownDemoCard)
    contents = split(read(card.path, String), "---\n")

    if length(contents) == 1
        config = Dict()
        body = join(contents)
    else
        config = YAML.load(strip(contents[2]))
        body = join(contents[3:end])
    end

    if !haskey(config, "cover")
        # set the first valid image path as cover
        # TODO: only markdown syntax is supported now
        image_paths = map(eachmatch(regex_md_img, body)) do m
            m.captures[1]
        end
        filter!(x->isfile(dirname(card.path), x), image_paths)
        if !isempty(image_paths)
            config["cover"] = first(image_paths)
        end
    end

    if haskey(config, "cover")
        config["cover"] = replace(config["cover"],
                                  r"[/\\]" => Base.Filesystem.path_separator) # windows compatibility
    end

    return config
end


text:
# This file is part of Kpax3. License is MIT.

function test_selection()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  state = [Kpax3.AminoAcidState(data, [1; 1; 1; 1; 1; 1], priorR, priorC, settings);
           Kpax3.AminoAcidState(data, [1; 2; 3; 4; 5; 6], priorR, priorC, settings);
           Kpax3.AminoAcidState(data, [1; 1; 2; 2; 3; 3], priorR, priorC, settings);
           Kpax3.AminoAcidState(data, [1; 1; 1; 1; 2; 2], priorR, priorC, settings);
           Kpax3.AminoAcidState(data, [1; 2; 3; 1; 2; 3], priorR, priorC, settings)]

  popsize = 5

  logpp = Float64[state[i].logpp for i in 1:popsize]

  # (1, 2) => 2 wins
  # (1, 3) => 1 wins
  # (1, 4) => 1 wins
  # (1, 5) => 1 wins
  # (2, 3) => 2 wins
  # (2, 4) => 2 wins
  # (2, 5) => 2 wins
  # (3, 4) => 4 wins
  # (3, 5) => 3 wins
  # (4, 5) => 4 wins
  #
  # possible candidates:
  #
  # (1, 2, 3, 4) => (2) (4)
  # (1, 2, 3, 5) => (2) (3)
  # (1, 2, 4, 5) => (2) (4)
  # (1, 3, 2, 4) => (1) (2)
  # (1, 3, 2, 5) => (1) (2)
  # (1, 3, 4, 5) => (1) (4)
  # (1, 4, 2, 3) => (1) (2)
  # (1, 4, 2, 5) => (1) (2)
  # (1, 4, 3, 5) => (1) (3)
  # (1, 5, 2, 3) => (1) (2)
  # (1, 5, 2, 4) => (1) (2)
  # (1, 5, 3, 4) => (1) (4)
  # (2, 3, 4, 5) => (2) (4)
  # (2, 4, 3, 5) => (2) (3)
  # (2, 5, 3, 4) => (2) (4)
  #
  # each winning index is repeated 8 times
  # there are a total of 240 possible combinations (15 * 8 * 2)
  pr = [72 / 240; 96 / 240; 24 / 240; 48 / 240; 0.0]

  N = 1000000
  tmp = zeros(Float64, popsize)
  i1 = 0
  i2 = 0

  for i in 1:N
    (i1, i2) = Kpax3.selection(logpp)

    tmp[i1] += 1
    tmp[i2] += 1
  end

  tmp /= 2 * N

  @test isapprox(tmp, pr, rtol=0.01)

  nothing
end

test_selection()


text:
#=
Create training file
=#

function set_up_training(han::Tracker_Handles,get_mean=true)
    woi=get_woi_array(han)
    set_up_training(han.nn,han.wt.vid_name,han.end_frame,woi,han.wt.pad_pos,han.frame_list,get_mean)
end

function set_up_training(nn,vid_name,max_frames,woi,pad_pos,frame_list,get_mean=true)

    (w,h,fps)=get_vid_dims(vid_name)

    if get_mean
        (mean_img,std_img)=mean_std_video_gpu(vid_name,max_frames)
        nn.norm.min_ref = 0
        nn.norm.max_ref = 255
        nn.norm.mean_img = mean_img
        nn.norm.std_img = std_img

        #Rotate and Reshape to 256 256
        nn.norm.mean_img = reshape(imresize(nn.norm.mean_img[:,:,1]',(256,256)),(256,256,1))
    end

    (new_woi, new_frame_list) = check_whiskers(woi,frame_list,max_frames)

    WT_reorder_whisker(new_woi,pad_pos)

    nn.labels = make_heatmap_labels(new_woi,pad_pos)
    nn.imgs = get_labeled_frames(vid_name,new_frame_list);

    #Normalize
    nn.imgs = StackedHourglass.normalize_new_images(nn.imgs,nn.norm.mean_img);

    (nn.imgs,nn.labels)=augment_images(nn.imgs,nn.labels);
end

function create_new_weights(nn)
    hg=HG2(size(nn.labels,1),size(nn.labels,3),4);
    load_hourglass(nn.weight_path,hg)
    StackedHourglass.change_hourglass(hg,size(nn.labels,1),1,size(nn.labels,3))
    nn.features=features(hg)
    nn.hg = hg
    nn.use_existing_weights=true
    nothing
end

#=
Run Training methods
=#

function make_training_batch(img,ll,batch_size=8)
    dtrn=minibatch(img,ll,batch_size,xtype=KnetArray,ytype=KnetArray)
end

function run_training(hg,trn::Knet.Data,this_opt,p,epochs=100,ls=Array{Float64,1}())

    total_length=length(trn) * epochs
    minimizer = Knet.minimize(hg,ncycle(trn,epochs),this_opt)
    last_update = 0.0
    count=0
    Gtk.set_gtk_property!(p, :fraction, 0)
    sleep(0.0001)

        for x in takenth(minimizer,1)
            push!(ls,x)
            count+=1
            complete=round(count/total_length,digits=2)
            if complete > last_update
                Gtk.set_gtk_property!(p, :fraction, complete)
                last_update = complete
                reveal(p,true)
            end
            sleep(0.0001)
        end

    ls
end

function run_training_no_gui(hg,trn::Knet.Data,this_opt,epochs=100,ls=Array{Float64,1}())

    total_length=length(trn) * epochs
    minimizer = Knet.minimize(hg,ncycle(trn,epochs),this_opt)

    for x in takenth(minimizer,1)
        push!(ls,x)
        sleep(0.0001)
    end

    ls
end

#=
Error Checking before training
=#

#Make sure that dimensions of model weights match label dimensions
function check_hg_features(nn)

    #Output features
    if size(nn.labels,3) != features(nn.hg)
        StackedHourglass.change_hourglass_output(nn.hg,size(nn.labels,1),size(nn.labels,3))
        nn.features = features(nn.hg)
    end
end

function check_whiskers(woi,frame_list,max_frames)

    #Make sure all of the frames actually have a tracked whisker
    empty=trues(length(woi))
    for i=1:length(woi)
        if empty_whiskers(woi[i])
            empty[i] = false
            println("You have an untraced whisker at index: ", i)
        end
    end

    #Last frames don't always work
    if (frame_list[end] == max_frames)
        empty[end] = false
    end

    new_woi = woi[empty]
    new_frame_list = frame_list[empty]

    (new_woi, new_frame_list)
end

function empty_whiskers(w)

    (length(w.x)==0)|(length(w.y)==0)

end


text:
import Base: size, getindex, Matrix, *, \, /
import LinearAlgebra: logdet

"""
    WoodburyMat{T} <: AbstractMatrix{T}

A lazily-represented matrix to which the Woodbury-Sherman-Morrison formula can be applied
for efficient operations. Is an `AbstractMatrix` `A`, given by `A := Xᵀ Σ^{-1} X + σ^2 I`.
Note that this is a rather specific instantiation of the types of matrices considered by the
matrix inversion and determinant lemmas. Should probably be generalised at some point.
"""
struct WoodburyMat{T<:Real, TX<:AM{T}, TΣ<:LazyPDMat{T}} <: AM{T}
    X::TX
    Σ::TΣ
    σ::T
    Γ::TX
    Ω::TΣ
    function WoodburyMat(X::TX, Σ::TΣ, σ::T) where {T<:Real, TX<:AM{T}, TΣ<:LazyPDMat{T}}
        Γ = chol(Σ)' \ X ./ σ
        Ω = LazyPDMat(Γ * Γ' + I, 0)
        return new{T, TX, TΣ}(X, Σ, σ, Γ, Ω)
    end
end

# AbstractArray interface.
size(A::WoodburyMat) = (size(A.X, 2), size(A.X, 2))
function getindex(A::WoodburyMat, p::Int, q::Int)
    return Xt_invA_Y(view(A.X, :, p), A.Σ, view(A.X, :, q)) + (p == q) * A.σ^2
end

# Conversions.
LazyPDMat(A::WoodburyMat) = Xt_invA_X(A.Σ, A.X) + A.σ^2 * I
Matrix(A::WoodburyMat) = Matrix(LazyPDMat(A))

# Some unary operations.
logdet(A::WoodburyMat) = logdet(A.Ω) + 2 * size(A, 1) * log(A.σ)

# Efficient matrix multiplication.
*(A::WoodburyMat, B::VecOrMat) = A.X' * (A.Σ \ (A.X * B)) .+ (A.σ^2) .* B
*(A::Matrix, B::WoodburyMat) = ((A * B.X') / B.Σ) * B.X .+ (B.σ^2) .* A

# Efficient linear system solving.
\(A::WoodburyMat, B::VecOrMat) = (B .- A.Γ' * (A.Ω \ (A.Γ * B))) ./ A.σ^2
/(B::VecOrMat, A::WoodburyMat) = (B .- ((B * A.Γ') / A.Ω) * A.Γ) ./ A.σ^2

# Efficient inverse-quadratic-form-like computations.
Xt_invA_X(A::WoodburyMat, X::VecOrMat) = (X'X .- Xt_invA_X(A.Ω, A.Γ * X)) ./ (A.σ^2)


text:
#slowwwwwwwwwwwwwwwwwww
#Curent iteration Patches `CurPat`, Noise level of Patches `Sigma_arr`
function patEstimation(NL_mat::Array{Int32,2}, Self_arr::Array{Int32,1}, Sigma_arr::Array{Float32,3}, CurPat::Array{Float32,3}, Par::PAR)
    Depth = size(CurPat)[3]

    EPat   = zeros(Float32, size(CurPat))
    W      = zeros(Int32, size(CurPat))
    #mean
    Temp = zeros(Float32, size(CurPat))
    M_Temp = zeros(Float32, size(CurPat))

    for i = 1 : length(Self_arr)
        Temp    =   CurPat[:, NL_mat[:, i], :]
        M_Temp  =   Float32.(mean(Temp, dims = 2))
        Temp    =   Temp .- M_Temp
        E_Temp 	=   WNNM(Temp, Par.c, Sigma_arr[1, Self_arr[i], :], M_Temp)

        EPat[:,NL_mat[:, i],:]    += E_Temp
        W[:,NL_mat[:, i], :]      += ones(Par.patsize^2, size(NL_mat[:, i])[1], Depth)
    end

    (EPat, W)
end


text:
module Wilsonloops_module

    mutable struct Wilson_loop
        wilson_loop::Array{Tuple{Int8,Int8},1}
        origin::NTuple{4,Int8}
        Wilson_loop() = new([],(0,0,0,0))
        #Wilson_loop(loop::Array{Tuple{Int,Int},1},origin) = new(loop,origin)
        #Wilson_loop(loop,origin) = new([Tuple(loop)],origin)
        #Wilson_loop(loop::Array{Tuple{Int,Int},1}) = new(loop,(0,0,0,0))

        function Wilson_loop(loop::Array{Tuple{T,T},1}) where T <: Integer
            return new(make_links(loop),(0,0,0,0))
        end
        function  Wilson_loop(loop::Array{Tuple{T,T},1},origin) where T<: Integer
            return new(make_links(loop),origin)
        end

        Wilson_loop(loop) = new([Tuple(loop)],(0,0,0,0))
    end

    mutable struct Tensor_wilson_lines
        left::Wilson_loop
        #Array{Tuple{Int8,Int8},1}
        right::Wilson_loop
        #Array{Tuple{Int8,Int8},1}
        position::NTuple{4,Int8}
        origin::NTuple{4,Int8}
        dUμ::NTuple{2,Int8}
        rightstart::NTuple{4,Int8}
    end

    function get_leftstartposition(w::Tensor_wilson_lines)
        return w.origin
    end

    function get_rightstartposition(w::Tensor_wilson_lines)
        return w.rightstart
    end

    mutable struct Tensor_wilson_lines_set
        loopset::Array{Tensor_wilson_lines,1}
        Tensor_wilson_lines_set() = new([])
    end

    function Base.getindex(g::Tensor_wilson_lines_set,i)
        return g.loopset[i]
    end

    function Base.length(x::Tensor_wilson_lines_set)
        return length(x.loopset)
    end

    function Base.push!(g::Tensor_wilson_lines_set,w::Tensor_wilson_lines)
        push!(g.loopset,w)
        return
    end

    function Base.display(w::Tensor_wilson_lines)
        println("tensor differenciated by ",w.dUμ)
        println("relative origin = ",w.origin)
        println("position = ",w.position)
        println("startpoint for leftpart: ",w.origin)
        println("startpoint for rightpart: ",w.rightstart)
        print("[")
        for i=1:length(w.left)
            print(w.left[i],",")
        end
        print("] ")
        print(" ⊗ ")
        print(" [")
        for i=1:length(w.right)
            print(w.right[i],",")
        end
        println("]\t")
    end

    function Base.display(w::Wilson_loop)
        println("origin = ",w.origin)
        print("[")
        for i=1:length(w)
            print(w.wilson_loop[i],",")
        end
        println("]\t")
    end

    function set_origin!(w::Wilson_loop,origin)
        w.origin = origin
    end

    mutable struct Wilson_loop_set
        #loopset::Set{Wilson_loop}
        #loopset::Array{Wilson,1}
        loopset::Array{Wilson_loop,1}
        Wilson_loop_set() = new([])
    end

    mutable struct Tensor_derivative_set
        tensors::Array{Array{Dict{Tuple{Int8,Int8},Tensor_wilson_lines_set},1},1}
        function Tensor_derivative_set(staples::Array{Wilson_loop_set,1})
            tensor_derivative = Array{Array{Dict{Tuple{Int8,Int8},Tensor_wilson_lines_set},1},1}(undef,4)
            for mu=1:4
                #println("$mu -direction")
                #display(staples[i][mu])
                tensor_derivative[mu] = Array{Dict{Tuple{Int8,Int8},Tensor_wilson_lines_set},1}(undef,length(staples[mu].loopset))

                for (k,link) in enumerate(staples[mu].loopset)
                    tensor_derivative[mu][k] = make_tensor_derivative(link)
                    #display(tensor_derivative[i][mu][k])
                end
            end
            return new(tensor_derivative)
        end
    end

    function Base.getindex(g::Tensor_derivative_set,i)
        return g.tensors[i]
    end

    function Base.display(w::Wilson_loop_set)
        for i=1:length(w)
            println("$i-th loops")
            display(w[i])
        end
        println("\t")
    end

    function Base.length(x::Wilson_loop_set)
        return length(x.loopset)
    end

    function Base.length(x::Wilson_loop)
        return length(x.wilson_loop)
    end

    function Base.getindex(g::Wilson_loop_set,i)
        return g.loopset[i]
    end

    function Base.getindex(g::Wilson_loop,i)
        return g.wilson_loop[i]
    end

    function Base.push!(g::Wilson_loop,link)
        push!(g.wilson_loop,Tuple(link))
        return
    end

    function Base.push!(g::Wilson_loop_set,w::Wilson_loop)
        push!(g.loopset,w)
        return
    end

    function make_cloverloops(μ,ν)
        loops = Wilson_loop_set()
        loop_righttop = Wilson_loop([(μ,1),(ν,1),(μ,-1),(ν,-1)])
        loop_lefttop = Wilson_loop([(ν,1),(μ,-1),(ν,-1),(μ,1)])
        loop_rightbottom = Wilson_loop([(ν,-1),(μ,1),(ν,1),(μ,-1)])
        loop_leftbottom= Wilson_loop([(μ,-1),(ν,-1),(μ,1),(ν,1)])
        push!(loops,loop_righttop)
        push!(loops,loop_lefttop)
        push!(loops,loop_rightbottom)
        push!(loops,loop_leftbottom)
        return loops
    end

    function make_originalactions_fromloops(coupling_loops)
        actions = Array{Wilson_loop_set,1}(undef,length(coupling_loops))
        for (i,loopset) in enumerate(coupling_loops)
            actions[i] = make_originalaction(loopset)
        end
        return actions
    end

    function make_loopforactions(couplinglist,L)
        loops = Array{Wilson_loop_set,1}(undef,length(couplinglist))
        for (i,name) in enumerate(couplinglist)
            if  name == "plaq"
                loops[i] = make_plaq()
            elseif name == "rect"
                loops[i] = make_rect()
            elseif name == "polyt"
                loops[i] = make_polyakov(4,L[4])
            elseif name == "polyx"
                loops[i] = make_polyakov(3,L[3])
            elseif name == "polyy"
                loops[i] = make_polyakov(2,L[2])
            elseif name == "polyz"
                loops[i] = make_polyakov(1,L[1])
            else
                error("$name is not supported!")
            end
        end
        return loops
    end

    function make_originalaction(loopset::Array{T,1}) where T <: Any
        loop_wilson = Wilson_loop_set()
        for loops in loopset
            @assert eltype(loops) == Tuple{<: Int,<: Int} "The loop should be like [(1, 1), (2, 1), (1, -1), (2, -1)]"
            loop1 = Wilson_loop(loops)
            push!(loop_wilson,loop1)
        end
        return loop_wilson
    end

    function make_plaqloops()
        loopset = []
        for μ=1:4
            for ν=μ:4
                #for ν=1:4
                if ν == μ
                    continue
                end
                push!(loopset,[(μ,1),(ν,1),(μ,-1),(ν,-1)])
            end
        end
        return loopset
    end

    function make_rectloops(N)
        loopset = []
        for μ=1:4
            for ν=μ:4
                if ν == μ
                    continue
                end
                push!(loopset,[(μ,1),(ν,N),(μ,-1),(ν,-N)])
                push!(loopset,[(μ,N),(ν,1),(μ,-N),(ν,-1)])
            end
        end
        return loopset
    end

    function make_polyakovloops(μ,Nμ)
        loopset = []
        push!(loopset,[(μ,Nμ)])
        return loopset
    end

    function make_plaq()
        loopset = make_plaqloops()
        loops = make_originalaction(loopset)
        return loops
    end

    function make_plaq(origin)
        loops = Wilson_loop_set()
        for μ=1:4
            #for ν=1:4
            for ν=μ:4
                if ν == μ
                    continue
                end
                loop = make_links([(μ,1),(ν,1),(μ,-1),(ν,-1)])
                loop1 = Wilson_loop(loop,Tuple(origin))
                push!(loops,loop1)
            end
        end

        return loops
    end

    function make_plaq(μ,ν)
        loops = Wilson_loop_set()
        origin = zeros(Int8,4)

        loop = make_links([(μ,1),(ν,1),(μ,-1),(ν,-1)])
        loop1 = Wilson_loop(loop,Tuple(origin))
        push!(loops,loop1)

        return loops
    end

    function make_plaq(μ,ν,origin)
        loops = Wilson_loop_set()
        #origin = zeros(Int8,4)

        loop = make_links([(μ,1),(ν,1),(μ,-1),(ν,-1)])
        loop1 = Wilson_loop(loop,Tuple(origin))
        push!(loops,loop1)

        return loops
    end

    function make_rect()
        loops = Wilson_loop_set()
        origin = zeros(Int8,4)

        for μ=1:4
            for ν=μ:4
                if ν == μ
                    continue
                end
                #loop = make_links([(μ,1),(ν,2),(μ,-1),(ν,-2)])

                loop1 = Wilson_loop([(μ,1),(ν,2),(μ,-1),(ν,-2)])
                #loop1 = Wilson_loop(loop,Tuple(origin))
                push!(loops,loop1)

                loop1 = Wilson_loop([(μ,2),(ν,1),(μ,-2),(ν,-1)])
                push!(loops,loop1)
            end
        end

        return loops
    end

    function make_polyakov(μ,Nμ)
        loops = Wilson_loop_set()
        origin = zeros(Int8,4)
        loop = make_links([(μ,Nμ)])
        loop1 = Wilson_loop(loop,Tuple(origin))
        push!(loops,loop1)
        return loops
    end

    function make_polyakov_xyz(Nμ)
        loops = Wilson_loop_set()
        for μ=1:3
            origin = zeros(Int8,4)
            loop = make_links([(μ,Nμ)])
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loops,loop1)
        end

        return loops
    end

    function make_chair()
        #loopset = []
        loopset = Wilson_loop_set()
        set1 = (1,2,3)
        set2 = (1,2,4)
        set3 = (2,3,4)
        set4 = (1,3,4)

        for set in (set1,set2,set3,set4)
            mu,nu,rho = set
            origin = zeros(Int8,4)
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loopset,loop1)

            mu,rho,nu = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loopset,loop1)

            nu,rho,mu = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loopset,loop1)

            nu,mu,rho = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loopset,loop1)

            rho,mu,nu = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loopset,loop1)

            rho,nu,mu = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loopset,loop1)

        end
        return loopset
    end

    function make_plaq_staple(μ)
        loops = Wilson_loop_set()
        origin = zeros(Int8,4)
        origin[μ] = 1
        for ν=1:4
            if ν == μ
                continue
            end
            loop = make_links([(ν,1),(μ,-1),(ν,-1)])
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loops,loop1)

            loop = make_links([(ν,-1),(μ,-1),(ν,1)])
            loop2 = Wilson_loop(loop,Tuple(origin))
            push!(loops,loop2)
        end
        return loops
    end

    function Base.adjoint(x::Wilson_loop_set)
        xout = Wilson_loop_set()
        num = length(x)

        for i=1:num
            xi = x[i]
            #display(xi)
            numloops = length(xi)
            coordinates = calc_coordinate(xi)

            #println("endpoint")
            endpoint = collect(coordinates[end])
            lastlink = xi[numloops]
            #println(lastlink)
            endpoint[lastlink[1]] += lastlink[2]
            #println(endpoint)

            loop = Tuple{Int64,Int64}[]
            for k=numloops:-1:1
                link = xi[k]
                push!(loop,(link[1],-link[2]))
            end
            #println(loop)
            loopw = make_links(loop)

            loop2 = Wilson_loop(loopw,Tuple(endpoint))

            #println(xi)
            #println("inv")
            #display(loop2)
            push!(xout,loop2)
        end
        #exit()
        return xout
    end

    function make_plaq_staple_prime(μ)
        loops = Wilson_loop_set()
        origin = zeros(Int8,4)
        #origin[μ] = 1
        for ν=1:4
            if ν == μ
                continue
            end
            loop = make_links([(ν,1),(μ,1),(ν,-1)])
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loops,loop1)

            loop2 = Wilson_loop([(ν,-1),(μ,1),(ν,1)],Tuple(origin))
            push!(loops,loop2)
        end
        return loops
    end

    function make_plaq_staple_prime(μ,Dim)
        loops = Wilson_loop_set()
        origin = zeros(Int8,4)
        #origin[μ] = 1
        for ν=1:Dim
            if ν == μ
                continue
            end
            loop = make_links([(ν,1),(μ,1),(ν,-1)])
            loop1 = Wilson_loop(loop,Tuple(origin))
            push!(loops,loop1)

            loop2 = Wilson_loop([(ν,-1),(μ,1),(ν,1)],Tuple(origin))
            push!(loops,loop2)
        end
        return loops
    end

    function calc_coordinate(w::Wilson_loop)
        numloops = length(w)
        coordinates = Array{NTuple{4,Int8},1}(undef,numloops)

        for k=1:numloops-1
            if k == 1
                coordinates[k] = w.origin
            end
            loopk = w[k]
            if loopk[1] == 1
                coordinates[k+1] = (coordinates[k][1]+loopk[2],coordinates[k][2],coordinates[k][3],coordinates[k][4])
            elseif loopk[1] == 2
                coordinates[k+1] = (coordinates[k][1],coordinates[k][2]+loopk[2],coordinates[k][3],coordinates[k][4])
            elseif loopk[1] == 3
                coordinates[k+1] = (coordinates[k][1],coordinates[k][2],coordinates[k][3]+loopk[2],coordinates[k][4])
            elseif loopk[1] == 4
                coordinates[k+1] = (coordinates[k][1],coordinates[k][2],coordinates[k][3],coordinates[k][4]+loopk[2])
            end
        end
        return coordinates
    end

    function calc_shift(w::Wilson_loop)
        numloops = length(w)
        coordinates = calc_coordinate(w)
        shifts = Array{NTuple{4,Int8},1}(undef,numloops)
        for k=1:numloops
            loopk = w[k]
            if loopk[2] == 1
                shifts[k] = coordinates[k]
            elseif loopk[2] == -1
                if loopk[1] == 1
                    shifts[k] = (coordinates[k][1]+loopk[2],coordinates[k][2],coordinates[k][3],coordinates[k][4])
                elseif loopk[1] == 2
                    shifts[k] = (coordinates[k][1],coordinates[k][2]+loopk[2],coordinates[k][3],coordinates[k][4])
                elseif loopk[1] == 3
                    shifts[k] = (coordinates[k][1],coordinates[k][2],coordinates[k][3]+loopk[2],coordinates[k][4])
                elseif loopk[1] == 4
                    shifts[k] = (coordinates[k][1],coordinates[k][2],coordinates[k][3],coordinates[k][4]+loopk[2])
                end
            end
        end
        return shifts
    end

    const kindsof_loops = ["plaquette","rectangular","chair","polyakov_x","polyakov_y","polyakov_z","polyakov_t"]

    function make_loops(loops_list,L)
        loops = []
        for loopname in loops_list
            if loopname == "plaquette"
                loop = make_plaq()
            elseif loopname == "rectangular"
                loop = make_rect()
            elseif loopname == "chair"
                loop = make_chair()
            elseif loopname == "polyakov_x"
                loop = make_polyakov(1,L[1])
            elseif loopname == "polyakov_y"
                loop = make_polyakov(2,L[2])
            elseif loopname == "polyakov_z"
                loop = make_polyakov(3,L[3])
            elseif loopname == "polyakov_t"
                loop = make_polyakov(4,L[4])
            else
                error("loop name $(loopname) is not supported!")
            end
            push!(loops,loop)
        end
        return loops
    end

    function make_links(segments::Array{Tuple{T,T},1}) where T <: Integer
        links = Tuple{Int8,Int8}[]
        for segment in segments
            s = sign(segment[2])
            if segment[2] == 0
                push!(links,(segment[1],0))
            else
            #@assert segment[2] != 0
                for i=1:abs(segment[2])
                    push!(links,(segment[1],s*1))
                end
            end
        end
        return links
    end

    function make_staples(w::Wilson_loop_set)
        num = length(w)
        staplesset = Array{Wilson_loop_set,1}(undef,4)
        tempset = Array{Array{Any,1},1}(undef,4)
        #tempset = Array{Set,1}(undef,4)
        for mu=1:4
            staplesset[mu] = Wilson_loop_set()
            #tempset[mu] = Set()
            tempset[mu] = []
        end

        for i=1:num
            staples = make_staples(w[i])
            for mu=1:4
                origin = zeros(Int8,4)
                origin[mu] = 1

                for k=1:length(staples[mu])
                    push!(tempset[mu],staples[mu][k])
                    #Wilson_loop(staples[mu][k],Tuple(origin)))
                end

                #for k=1:length(staples[mu])
                #    push!(staplesset[mu],Wilson_loop(staples[mu][k],Tuple(origin))  )
                #end
            end
        end

        for mu=1:4
            #display(tempset[mu])
            #exit()
            origin = zeros(Int8,4)
            origin[mu] = 1
            for temp in tempset[mu]
                push!(staplesset[mu],Wilson_loop(temp,Tuple(origin)))
            end
        end

        return staplesset
    end

    function get_coordinate(w::Wilson_loop,position)
        coor = collect(w.origin)
        if position != 1
            for i=1:position-1
                link = w.wilson_loop[i]
                coor[link[1]] += link[2]
            end
        end

        return coor
    end

    #=
    mutable struct Tensor_wilson_lines
        left::Array{Tuple{Int8,Int8},1}
        right::Array{Tuple{Int8,Int8},1}
        position::NTuple{4,Int8}
    end
    =#

    function make_tensor_derivative(w::Wilson_loop)
        dVdU = Dict{Tuple{Int8,Int8},Tensor_wilson_lines_set}()
        #display(w)
        for μ=1:4
            for sign = [+1,-1]
                u = (μ,sign)
                #println("u = $u")
                dVdU[u] = make_tensor_Wilson_line(w,u)
            end
        end
        return dVdU
    end

    function make_tensor_Wilson_line(w::Wilson_loop,u::Tuple{T,T}) where T <: Int
        left = Tuple{Int8,Int8}[]
        right = Tuple{Int8,Int8}[]
        count,positions = count_terms(w,u)
        len = length(w)

        tw = Tensor_wilson_lines_set()
        #println("u = $u count = ",count," position = $positions")
        if count != 0
            #println("u = $u")
            for i=1:count
                left = Tuple{Int8,Int8}[]
                right = Tuple{Int8,Int8}[]

                position = positions[i]
                coor = get_coordinate(w,position)
                #println("position ", coor)
                #insertpoint = Tuple(collect(w.origin) .- coor)
                insertpoint = Tuple(- collect(coor))
                #println("position2 ", insertpoint )

                #println(positions[i])
                if position == 1
                    #println(u)
                    push!(left,(0,0))
                    #display(w)
                else
                    for k=1:position-1
                        push!(left,w.wilson_loop[k])
                    end
                end
                #println("left = ",left)

                if position == len
                    #if u == (1,1)
                    #    println(right)
                    #end
                    push!(right,(0,0))
                else
                    for k=position+1:len
                        push!(right,w.wilson_loop[k])
                    end
                end
                #println("right = ",right)
                relativeorigin = [0,0,0,0]
                if left[1] == (0,0)
                else
                    for i=length(left):-1:1
                        link = left[i]
                        relativeorigin[link[1]] -= link[2]
                    end
                end

                temp_right = [0,0,0,0]
                temp_right[u[1]] += u[2]

                left_wilson = Wilson_loop(make_links(left),Tuple(relativeorigin))
                right_wilson = Wilson_loop(make_links(right),Tuple(temp_right))

                #tensorlink = Tensor_wilson_lines(left,right,insertpoint,Tuple(relativeorigin),u,Tuple(temp_right))
                tensorlink = Tensor_wilson_lines(left_wilson,right_wilson,insertpoint,Tuple(relativeorigin),u,Tuple(temp_right))

                #=
                if position == len
                    if u == (1,1)
                        display(tensorlink)
                    end
                end
                =#

                #display(tensorlink)
                push!(tw,tensorlink)
            end

        end
        return tw
    end

    function count_terms(w::Wilson_loop,u)
        count = 0
        positions = []
        #origin = w.origin
        #x = collect(origin)
        for (i,link) in enumerate(w.wilson_loop)

            if link == u
                count += 1
                push!(positions,i)
            end
            #x[link[1]] += link[2]
            #println(x)
            #println("$link,$count,$positions")
        end
        return count,positions
    end

    function make_staple_staple(staples::Array{Wilson_loop,1})
        staplestaple = Array{Wilson_loop,2}(undef,4,4)
        for ν=1:4
            staple = make_staples(staples[ν])
            for μ=1:4
                staplestaple[ν,μ] = staple[μ]
            end
        end
        return staplestaple
    end

    function make_staples(w::Wilson_loop)
        staples = make_staples(w.wilson_loop)
        return staples[1],staples[2],staples[3],staples[4]
    end

    function make_staples(links)
        numlinks = length(links)
        #println("num. of links: $numlinks")
        staples = Array{Array{Array{Tuple{Int,Int},1}}}(undef,4)
        for mu=1:4
            linkmu = []
            for (i,link) in enumerate(links)
                if link[1] == mu
                    append!(linkmu,i)
                end
            end
            #println(mu,"\t",linkmu)

            staple_mu = []

            if length(linkmu) == 0
                staple = Tuple{Int,Int}[]
                staple_mu = Array{Tuple{Int,Int},1}[]
            end

            staple_mu = Array{Array{Tuple{Int,Int},1}}(undef,length(linkmu))
            count = 0
            for i in linkmu
                count += 1
                staple = Array{Tuple{Int,Int}}(undef,numlinks-1)
                if i != numlinks
                    staple[1:numlinks-i] = links[i+1:end]
                end

                if i != 1
                    staple[numlinks-i+1:end] = links[1:i-1]
                end

                if links[i][2] == -1
                    staple = staple[end:-1:1]
                    for j = 1:numlinks-1
                        staple[j] = (staple[j][1],-staple[j][2])
                    end
                end
                #println(staple)
                staple_mu[count] = staple[:]
            end
            staples[mu] = staple_mu
            #println(staple_mu)
        end

        #println(staples)
        return staples
    end

    function calc_loopset_μν_name(name)
        loops_μν= Array{Wilson_loop_set,2}(undef,4,4)
        if name == "plaq"
            numofloops = 1
            for μ=1:4
                for ν=1:4
                    if ν == μ
                        continue
                    end
                    loops = Wilson_loop_set()
                    loop = Wilson_loop([(μ,1),(ν,1),(μ,-1),(ν,-1)])
                    push!(loops,loop)
                    loops_μν[μ,ν] = loops
                end
            end
        elseif name == "clover"
            numofloops = 4
            for μ=1:4
                for ν=1:4
                    if ν == μ
                        continue
                    end
                    loops = Wilson_loop_set()

                    loop_righttop = Wilson_loop([(μ,1),(ν,1),(μ,-1),(ν,-1)])
                    loop_lefttop = Wilson_loop([(ν,1),(μ,-1),(ν,-1),(μ,1)])
                    loop_rightbottom = Wilson_loop([(ν,-1),(μ,1),(ν,1),(μ,-1)])
                    loop_leftbottom= Wilson_loop([(μ,-1),(ν,-1),(μ,1),(ν,1)])
                    push!(loops,loop_righttop)
                    push!(loops,loop_lefttop)
                    push!(loops,loop_rightbottom)
                    push!(loops,loop_leftbottom)

                    loops_μν[μ,ν] = loops
                end
            end
        elseif name == "rect"
            numofloops = 8
            for μ=1:4
                for ν=1:4
                    if ν == μ
                        continue
                    end
                    loops = Wilson_loop_set()

                    loop_righttop = Wilson_loop([(μ,2),(ν,1),(μ,-2),(ν,-1)])
                    loop_lefttop = Wilson_loop([(ν,1),(μ,-2),(ν,-1),(μ,2)])
                    loop_rightbottom = Wilson_loop([(ν,-1),(μ,2),(ν,1),(μ,-2)])
                    loop_leftbottom= Wilson_loop([(μ,-2),(ν,-1),(μ,2),(ν,1)])
                    push!(loops,loop_righttop)
                    push!(loops,loop_lefttop)
                    push!(loops,loop_rightbottom)
                    push!(loops,loop_leftbottom)

                    loop_righttop = Wilson_loop([(μ,1),(ν,2),(μ,-1),(ν,-2)])
                    loop_lefttop = Wilson_loop([(ν,2),(μ,-1),(ν,-2),(μ,1)])
                    loop_rightbottom = Wilson_loop([(ν,-2),(μ,1),(ν,2),(μ,-1)])
                    loop_leftbottom= Wilson_loop([(μ,-1),(ν,-2),(μ,1),(ν,2)])
                    push!(loops,loop_righttop)
                    push!(loops,loop_lefttop)
                    push!(loops,loop_rightbottom)
                    push!(loops,loop_leftbottom)

                    loops_μν[μ,ν] = loops
                end
            end
        else
            error("$name is not supported")
        end
        return loops_μν,numofloops
    end

end


text:
module DynamicOverlay

using Mixtape
import Mixtape: CompilationContext, transform, allow, show_after_inference,
                show_after_optimization, debug, @load_call_interface
using MacroTools
using BenchmarkTools

foo(x) = x^5
bar(x) = x^10
apply(f, x1, x2::Val{T}) where {T} = f(x1, T)

function f(x)
    g = x < 5 ? foo : bar
    return g(2)
end

struct MyMix <: CompilationContext end

allow(ctx::MyMix, m::Module) = m == DynamicOverlay
show_after_inference(ctx::MyMix) = false
show_after_optimization(ctx::MyMix) = false
debug(ctx::MyMix) = false

swap(e) = e
function swap(e::Expr)
    new = MacroTools.postwalk(e) do s
        isexpr(s, :call) || return s
        s.args[1] == Base.literal_pow || return s
        return Expr(:call, apply, Base.:(*), s.args[3:end]...)
    end
    return new
end

function transform(::MyMix, b)
    for (v, st) in b
        replace!(b, v, swap(st))
    end
    return b
end

# JIT compile an entry and call.
fn = Mixtape.jit(MyMix(), f, Tuple{Int64})
display(fn(3))
display(fn(6))
@btime fn(6)

# Mixtape cached call.
Mixtape.@load_call_interface()
display(call(MyMix(), f, 3))
display(call(MyMix(), f, 6))
@btime call(MyMix(), f, 6)

# Native.
f(5)
@btime f(5)

end # module


text:
using Blink

function blink_show(v::VegaVisualization)
    w = Window()
    path(url...) = joinpath(dirname(@__FILE__), "..", "assets", "bower_components", url...)
    d3 = path("d3","d3.min.js")
    topojson = path("topojson","topojson.js")
    cloudlayout = path("d3-cloud", "build", "d3.layout.cloud.js")
    vega = path("vega", "vega.js")
    for file in (d3, topojson, cloudlayout, vega)
        load!(w, file)
    end
    body!(w, body(v))
end

max_stars_repo_name: jpthor/vega


text:
#logmod.jl

using JuMP
using Ipopt

M = parse(Int,ARGS[1])
N = parse(Int,ARGS[2])

include(string("log",M,"_",N,"_dat.jl"))

m = Model()

simplify = true
@show simplify
if parse(Int, ARGS[3]) == 1
    m = Model(solver=IpoptSolver(), simplify_nonlinear_expressions=simplify)
    #use jump
elseif parse(Int, ARGS[3]) == 2
    include("interface_tape.jl")
    using TapeInterface
    m = Model(solver=TapeSolver(IpoptSolver()),simplify_nonlinear_expressions=simplify)
elseif parse(Int, ARGS[3]) == 3
    include("interface_tape_dict.jl")
    using TapeInterface
    m = Model(solver=TapeSolver(IpoptSolver()),simplify_nonlinear_expressions=simplify)
else
    @assert false
end

@variable(m, theta[1:N], start=-0.0001)
# @variable(m, theta2[1:N], start=-0.0001)

# @variable(m, sum1[1:M], start=1.0)
# for i = 1:M
#     @constraint(m, y[i] * sum{ theta[j] * x[i,j], j = 1:N } == sum1[i])
# end

# testing expression
# @NLobjective(m, Min, sum{ log(1+1/exp(sum{theta[j], j=1:N})) ,i=1:M})
#@NLobjective(m, Min, sum{ log(1+1/exp(sum{theta[j]+theta2[j], j=1:N})) ,i=1:M})
# @NLobjective(m, Min, sum{ log(1+1/exp(sum{theta[j] + theta2[j], j=1:N})) ,i=1:M})

#original logrithmic regression (no parameter)
#@NLobjective(m, Min, sum{ theta[i]^2, i = 1:N} + sum{ log(1+1/exp(sum{theta[j], j=1:N})), i =1:M})
#original logrithmic regression
@NLobjective(m, Min, lambda*(sum{ theta[i]^2, i = 1:N}) + sum{ log(1+1/exp(y[i] * sum{ theta[j] * x[i,j]  ,j = 1:N } )), i =1:M})

#remvoing parameter y
#@NLobjective(m, Min, lambda*(sum{ theta[i]^2, i = 1:N}) + sum{ log(1+1/exp(sum{ theta[j] * x[i,j]  ,j = 1:N } )), i =1:M})

#log x  -> x^2
#@NLobjective(m, Min, lambda*(sum{ theta[i]^2, i = 1:N}) + sum{ (1+1/exp(y[i] * sum{ theta[j] * x[i,j]  ,j = 1:N } ))^2, i =1:M})

#exp x - > x^2
#@NLobjective(m, Min, lambda*(sum{ theta[i]^2, i = 1:N}) + sum{ log(1+1/(y[i] * sum{ theta[j] * x[i,j]  ,j = 1:N } )^2), i =1:M})

# 1/exp x -> x^2
# @NLobjective(m, Min, lambda*(sum{ theta[i]^2, i = 1:N}) + sum{ log(1+(y[i] * sum{ theta[j] * x[i,j]  ,j = 1:N } )^2 ), i =1:M})

solve(m)


text:
module LloydsMirror
function lloydsmirror_singlereflection(c, f, r_src, r_hyd, z_src, z_hyd)
	λ = c/f
	k = 2π/λ

	s_dir = sqrt((r_src - r_hyd)^2 + (z_src - z_hyd)^2)
	s_rfl = sqrt((r_src - r_hyd)^2 + (z_src + z_hyd)^2)

	pressure(s) = exp(im*k*s)/s
	p_dir = pressure(s_dir)
	p_rfl = -pressure(s_rfl)

	p = p_dir + p_rfl
	TL = -20log10(abs(p))
end
end


text:
tests = [
    "tables",
    "randtn"]

print_with_color(:blue, "Running tests:\n")

srand(345678)

for t in tests
    test_fn = "$t.jl"
    print_with_color(:green, "* $test_fn\n")
    let
        include(test_fn)
    end
end

max_stars_repo_name: jdrugo/fasttn

max_issues_repo_name: jdrugo/fasttn

max_forks_repo_name: jdrugo/fasttn