text:
# Based on Matlab source by Fredrik Fryklund and Erik Lehto

module PUX

mydir = Base.Filesystem.dirname(Base.source_path())
const PUX_MATLAB_SRC = "$(mydir)/../../external/pux"
const RBFQR_MATLAB_SRC = "$(mydir)/../../external/rbf-qr"

using MATLAB
using PyPlot
using CurveDiscretization

export pux_params, pux_setup!, pux_eval

function __init__()
    mat"clear()"
    mat"addpath(genpath($PUX_MATLAB_SRC))"
    mat"addpath(genpath($RBFQR_MATLAB_SRC))"
end

mutable struct PUXStorage
    numcenters::Int64
    PUXParams::Dict
    PUXData::Dict
    PUXDataInit::Bool
end

function pux_precompute(curves, dcurve, Ngrid, ep; P=0, R=0.0)
    Ngrid_int = Ngrid
    if P==0
        # Determine P from radius
        hgrid = 2*maximum(abs.(vec(dcurve.points))) / Ngrid
        P = ceil(R/hgrid)
    end
    # Setup volume grid
    Lgrid, Ngrid, X, Y = setup_volgrid(dcurve, Ngrid_int, P)
    xe = [vec(X) vec(Y)]
    xet = xe'
    interior = interior_points(dcurve, xet)
    # Setup PUX params, one setup per boundary
    if !isa(curves, Array)
        curves = [curves]
    end
    Ncurves = length(curves)
    PUXStor = Array{PUXStorage}(Ncurves)
    for i=1:Ncurves
        curve_i = curves[i]
        idx = findin(dcurve.curvenum, i)
        arclength = sum(dcurve.dS[idx])
        PUXStor[i] = pux_params(Lgrid/2, Ngrid, ep, P, arclength)
        # Get equispaced points
        _, z_equi = traparclen(curve_i, PUXStor[i].numcenters)
        z_grid = dcurve.points[1,idx] + 1im*dcurve.points[2,idx]
        # Precompute
        pux_setup!(xe, interior, z_grid, z_equi, PUXStor[i])
    end
    # Return everything that the solver needs
    return PUXStor, X, Y, Lgrid, Ngrid, interior
end

function pux_params(Lgrid, Ngrid, shape, support, arclength)
    M = Float64(Ngrid)
    LBox = Float64(Lgrid)
    ep = Float64(shape)
    P = Float64(support)
    arcL = Float64(arclength)
    # Percentage of parition radius when a zero-paritions will be removed
    R_ZeroPart_threshold = 0.8
    # radius of partition ratio (1 = coverage, >1 overlap)
    R_ratio = 2.5
    hgrid = 2.0*LBox/M

    R_part = P*hgrid
    R = R_part/R_ratio

    info("PUX PARAMS:")
    info("P = ", P)
    info("R = ", R)
    info("R_part = ", R_part)

    regularity = 0
    if regularity == 0
        regularity = min(floor(sqrt(P)-0.9), 5) # test (?)
    end
    # Number of Vogel points, i.e. RBF-centres.
    n_RBFCentres = round(min(0.8*pi/4*P^2, 4*P))
    n_RBFCentres = min(n_RBFCentres, 200)
    info("Number of Vogel nodes: ", n_RBFCentres)

    # Number of partitions
    n_Part = ceil(arcL/(2*R)) + 1

    # Pack into Dict
    PUXParams = Dict{String,Any}(
        "R_ZeroPart_threshold" => R_ZeroPart_threshold,
        "R_part" =>   R_part,
        "regularity" => regularity,
        "n_RBFCentres" =>         n_RBFCentres,
        "n_Part" =>                n_Part,
        "M" => M,
        "LBox" => LBox,
        "ep" => ep,
        "arcL" => arcL
    )
    # Return PUXStorgae
    stor = PUXStorage(PUXParams["n_Part"]*2, PUXParams, similar(PUXParams), false)
    return stor
end

function pux_setup!(xe::Array{Float64,2},
                    interior::BitArray{1},
                    z_grid::Vector{Complex{Float64}},
                    z_equi::Vector{Complex{Float64}},
                    stor::PUXStorage)
    params = stor.PUXParams
    mat"$PUXData = pux_setup($xe, $interior, $z_grid, $z_equi, $params)"
    stor.PUXData = PUXData
    stor.PUXDataInit = true
end

function pux_eval(f1, stor)
    # Dummy wrapper for vector PUX
    f2 = zeros(f1)
    u1, u2 = pux_eval_vec(f1, f2, stor)
    return u1
end

function pux_eval_vec(f1::Array{Float64},
                      f2::Array{Float64},
                      stor::PUXStorage);
                      limit_od::Float64=3.5
    pux_eval_vec(f1, f2, [stor])
end

function pux_eval_vec(f1::Array{Float64},
                      f2::Array{Float64},
                      stor_list::Array{PUXStorage};
                      limit_od::Float64=3.5)
    fe1 = zeros(size(f1))
    fe2 = zeros(size(f2))
    for curvenum = 1:length(stor_list)
        stor = stor_list[curvenum]
        @assert stor.PUXDataInit "Must run pux setup first!"
        data = stor.PUXData
        # Extraxt data from MATLAB struct
        # Everything is Float64 when we get it, convert explicitly
        A = data["A"]
        W = data["W"]
        If = map(Int64, data["If"])
        Jf = map(Int64, data["Jf"])
        n_Boxpnts = Int64(data["n_Boxpnts"])
        n_Part = Int64(data["n_Part"])
        n_InterpPart = Int64( data["n_InterpPart"] )
        n_extensionPnts = Int64(data["n_extensionPnts"])
        idx_xeOutOmegaInInterpPart = map(Int64, data["idx_xeOutOmegaInInterpPart"])
        idx_xeInOmega_boolean = map(Bool, data["idx_xeInOmega_boolean"])
        idx_InterpPartCentres = map(Int64, data["idx_InterpPartCentres"])
        idx_xeInInterpPart_stencil = map(Int64, data["idx_xeInInterpPart_stencil"])
        # Coarsen sets how to sample the evaluation grid to obtain the
        # interpolation grid. Set to 1 if they should be equal.
        coarsen = 1
        # Find values in collocation points
        f1_collocation = zeros(n_Boxpnts)
        f2_collocation = zeros(n_Boxpnts)
        f1_collocation[idx_xeInOmega_boolean] = f1[idx_xeInOmega_boolean]
        f2_collocation[idx_xeInOmega_boolean] = f2[idx_xeInOmega_boolean]
        # Create coarser grid for interpolation
        if (coarsen > 1)
            # NEEDS TO BE TESTED
            warn("Coarsening is untested in Julia implementation")
            M = Int64( stor.PUXParams["M"] )
            I, J = ndgrid(1:coarsen:M, 1:coarsen:M)
            I = I'
            J = J'
            i_crs = I+M*(J-1)
            i_crs = vec(i_crs)
        end
        # Vector for data used for building extension. The names S is
        # based on the input for MATLABS sparse:
        # S = sparse(i,j,s,m,n) uses vectors i, j, and s to generate an
        #     m-by-n sparse matrix such that S(i(k),j(k)) = s(k)
        # Sf stores information about the local extensions.
        Sf1 = zeros(n_extensionPnts)
        Sf2 = zeros(n_extensionPnts)
        # Counters used to fill in Sf.
        idx_endf = 0
        max_resnorm = 0.0
        ext_rel_magn = 0.0
        for i = 1:n_InterpPart
            ## Global indices
            # Indices of point in xe within R_Ip of parition center Ip(i).
            idx_xeInInterpPart_i = idx_InterpPartCentres[i] + idx_xeInInterpPart_stencil
            idx_xeInOmegaInInterpPart_i = idx_xeInInterpPart_i[idx_xeInOmega_boolean[idx_xeInInterpPart_i]]
            idx_xeInInterpPartOutOmega = idx_xeInInterpPart_i[.~idx_xeInOmega_boolean[idx_xeInInterpPart_i]]
            n_local_evals = length(idx_xeInInterpPartOutOmega)
            # Local indices
            idx_local_xeInOmega = find(idx_xeInOmega_boolean[idx_xeInInterpPart_i])
            idx_local_xeOutOmega = find(.~idx_xeInOmega_boolean[idx_xeInInterpPart_i])

            if (coarsen > 1)
                # Coarser grid. Samples points set by coarsen.
                ic, ij = intersect(idx_xeInInterpPart_i[idx_xeInOmega_boolean[idx_xeInInterpPart_i]],i_crs) # global
                idx_xeInOmegaInInterpPart_i = ic[idx_xeInOmega_boolean[ic]]  # global
                idx_local_xeInOmega = idx_local_xeInOmega[ij]
            end
            # Solve overdetermined system to obtain f in RBF centres.
            MAT = A[idx_local_xeInOmega,:]
            rhs = hcat(f1_collocation[idx_xeInOmegaInInterpPart_i], f2_collocation[idx_xeInOmegaInInterpPart_i])

            # Limit overdeterminedness
            rows, cols = size(MAT)
            if rows/cols > limit_od
                newrows = limit_od*cols
                idx = 1:floor(Integer, rows/newrows):rows
                f_RBFCentres = MAT[idx, :]\rhs[idx, :]
            else
                f_RBFCentres = MAT\rhs
            end
            # Compute residual
            resnorm = rel = norm(MAT*f_RBFCentres - rhs) / norm(rhs)
            max_resnorm = max(max_resnorm, resnorm)
            # Extrapolate
            Af = A[idx_local_xeOutOmega,:] * f_RBFCentres
            Sf1[idx_endf+1:idx_endf+n_local_evals] = Af[:,1]
            Sf2[idx_endf+1:idx_endf+n_local_evals] = Af[:,2]
            idx_endf = idx_endf + n_local_evals
            # Heuristic sanity check of extrapolation
            ext_rel_magn = max(ext_rel_magn, norm(Af, Inf) / norm(rhs, Inf))
        end
        #info("PUX: maximum residual: $max_resnorm")
        #info("PUX: extension rel. magnitude: $ext_rel_magn")
        if ext_rel_magn > 10
            #warn("PUX: Partitions too large?")
        end

        # cell_idx_xeOutOmegaInInterpPart is obtained by range search from
        #partition centers on uniform grid. Thus they are not the same amount as
        #idx_local_xeOutOmega_boolean. idx_local_xeOutOmega_boolean
        # comes from inpolygon.
        Sf1 = Sf1[1:idx_endf]
        Sf2 = Sf2[1:idx_endf]
        Imat1 = sparse(If, Jf, Sf1, n_Boxpnts, n_Part)
        Imat1 = Imat1[idx_xeOutOmegaInInterpPart,:]
        Imat2 = sparse(If, Jf, Sf2, n_Boxpnts, n_Part)
        Imat2 = Imat2[idx_xeOutOmegaInInterpPart,:]
        ## Combine local extensions into global extension
        fe1[idx_xeOutOmegaInInterpPart] = full(sum(Imat1.*W, 2))
        fe1[idx_xeInOmega_boolean] = f1_collocation[idx_xeInOmega_boolean]
        fe2[idx_xeOutOmegaInInterpPart] = full(sum(Imat2.*W, 2))
        fe2[idx_xeInOmega_boolean] = f2_collocation[idx_xeInOmega_boolean]
    end
    return fe1, fe2
end

function plot_partitions(stor::PUXStorage)
    function partplot(clist, Rlist, style)
        th = linspace(0, 2*pi)
        for i=1:size(clist,1)
            c = clist[i,:]
            plot(c[1] + Rlist[i]*cos.(th), c[2] + Rlist[i]*sin.(th), style)
        end
    end
    PUXData = stor.PUXData
    partplot(PUXData["Ip"], PUXData["R_Ip"]*ones(PUXData["n_InterpPart"]), "r")
    partplot(PUXData["Op"], PUXData["R_Op"], "g-")
end

function plot_partitions(arr::Array{PUXStorage})
    for s in arr
        plot_partitions(s)
    end
end

# Wrappers to original Matlab code
function pux_params_matlab(Lgrid, Ngrid, shape, support, arclength)
    # Convert to Float64 (MATLAB default type)
    pux_ep = Float64(shape)
    pux_P = Float64(support)
    Lgrid = Float64(Lgrid)
    Ngrid = Float64(Ngrid)
    arclength = Float64(arclength)
    mat"$PUXParams = pux_params($Ngrid, $Lgrid, $pux_ep, $pux_P, $arclength)"
    stor = PUXStorage(PUXParams["n_Part"]*2, PUXParams, similar(PUXParams), false)
    return stor
end

function pux_eval_matlab(fin::Array{Float64},
                         stor::PUXStorage)
    data = stor.PUXData
    mat"$fext = pux_eval($fin, $data)"
    return reshape(fext, size(fin))
end

function setup_volgrid(grid::DiscreteCurve, Ngrid, pux_P)
    # Setup minimal grid size, ***assuming curve is centered at origin***
    Dmax = maximum(abs.(vec(grid.points)))
    Lgrid = 2*Dmax
    hgrid = Lgrid/Ngrid
    padding = 2*pux_P+2 # This should always accomodate PUs (which is right, 2*P or 4*P ???))
    Ngrid += padding
    Ngrid += Ngrid%2 # Even FFT grid is better
    info("After padding Ngrid=$Ngrid")
    Lgrid += padding*hgrid
    Lgrid *= exp(0.001)  # fudge factor to reduce number of grid points exactly on bdry
    xgrid = linspace(-Lgrid/2, Lgrid/2, Ngrid+1)
    ygrid = linspace(-Lgrid/2, Lgrid/2, Ngrid+1)
    xgrid = xgrid[1:end-1]
    ygrid = ygrid[1:end-1]
    X, Y = ndgrid(xgrid, ygrid)
    Lgrid, Ngrid, X, Y
end

# Matlab compability
function ndgrid(x, y)
    nx = length(x)
    ny = length(y)
    vartype = typeof(x[1])
    @assert typeof(y[1])==vartype
    X = Array{vartype}(nx, ny)
    Y = Array{vartype}(nx, ny)
    for i=1:nx
        for j=1:ny
            X[i,j] = x[i]
            Y[i,j] = y[j]
        end
    end
    return X,Y
end

end


text:
using ImageFiltering.Models

@testset "solve_ROF_PD" begin
    # This testset is modified from its CPU version

    @testset "Numerical" begin
        # 2D Gray
        img = restrict(testimage("cameraman"))
        img_noisy = img .+ 0.05randn(MersenneTwister(0), size(img))
        img_smoothed = solve_ROF_PD(img_noisy, 0.05, 20)
        @test ndims(img_smoothed) == 2
        @test eltype(img_smoothed) <: Gray
        @test assess_psnr(img_smoothed, img) > 31.67
        @test assess_ssim(img_smoothed, img) > 0.90

        img_noisy_cu = CuArray(float32.(img_noisy))
        img_smoothed_cu = solve_ROF_PD(img_noisy_cu, 0.05, 20)
        @test img_smoothed_cu isa CuArray
        @test eltype(eltype(img_smoothed_cu)) == Float32
        @test Array(img_smoothed_cu) ≈ img_smoothed

        # 2D RGB
        img = restrict(testimage("lighthouse"))
        img_noisy = img .+ colorview(RGB, ntuple(i->0.05.*randn(MersenneTwister(i), size(img)), 3)...)
        img_smoothed = solve_ROF_PD(img_noisy, 0.03, 20)
        @test ndims(img_smoothed) == 2
        @test eltype(img_smoothed) <: RGB
        @test assess_psnr(img_smoothed, img) > 32.15
        @test assess_ssim(img_smoothed, img) > 0.90

        img_noisy_cu = CuArray(float32.(img_noisy))
        img_smoothed_cu = solve_ROF_PD(img_noisy_cu, 0.03, 20)
        @test img_smoothed_cu isa CuArray
        @test eltype(eltype(img_smoothed_cu)) == Float32
        @test Array(img_smoothed_cu) ≈ img_smoothed

        # 3D Gray
        img = Gray.(restrict(testimage("mri"), (1, 2)))
        img_noisy = img .+ 0.05randn(MersenneTwister(0), size(img))
        img_smoothed = solve_ROF_PD(img_noisy, 0.02, 20)
        @test ndims(img_smoothed) == 3
        @test eltype(img_smoothed) <: Gray
        @test assess_psnr(img_smoothed, img) > 31.78
        @test assess_ssim(img_smoothed, img) > 0.85

        img_noisy_cu = CuArray(float32.(img_noisy))
        img_smoothed_cu = solve_ROF_PD(img_noisy_cu, 0.02, 20)
        @test img_smoothed_cu isa CuArray
        @test eltype(eltype(img_smoothed_cu)) == Float32
        @test Array(img_smoothed_cu) ≈ img_smoothed

        # 3D RGB
        img = RGB.(restrict(testimage("mri"), (1, 2)))
        img_noisy = img .+ colorview(RGB, ntuple(i->0.05.*randn(MersenneTwister(i), size(img)), 3)...)
        img_smoothed = solve_ROF_PD(img_noisy, 0.02, 20)
        @test ndims(img_smoothed) == 3
        @test eltype(img_smoothed) <: RGB
        @test assess_psnr(img_smoothed, img) > 31.17
        @test assess_ssim(img_smoothed, img) > 0.79

        img_noisy_cu = CuArray(float32.(img_noisy))
        img_smoothed_cu = solve_ROF_PD(img_noisy_cu, 0.02, 20)
        @test img_smoothed_cu isa CuArray
        @test eltype(eltype(img_smoothed_cu)) == Float32
        @test Array(img_smoothed_cu) ≈ img_smoothed
    end
end

ext: jl


text:
## ----------  SingleId

"""
    SingleId

Id for one object or object vector.
Does not keep track of location in model.
Keeps track of `name` (Symbol identifier), `index`, and `description` (for display purposes).

`index` is used when there is a vector or matrix of objects of the same type
`index` is typically empty (scalar object) or scalar (vector objects)
"""
struct SingleId
    name :: Symbol
    index :: Vector{Int}
    description :: String
end

"""
	$(SIGNATURES)

Constructor for `SingleId` from `name` only.
"""
SingleId(name :: Symbol) = SingleId(name, Vector{Int}(), "");

SingleId(name :: Symbol, idx :: T1, descr = "") where T1 <: Integer =
    SingleId(name, [idx], descr);

"""
	$(SIGNATURES)

Constructor for a `SingleId` with index.

# Example
```
SingleId(:myId, [2, 1])
```
"""
SingleId(name :: Symbol, idxM :: Vector{I1}) where I1 <: Integer =
    SingleId(name, idxM, "");

SingleId(name :: Symbol, descr :: String) = SingleId(name, Vector{Int}(), descr);

## ----------  Access

Base.broadcastable(s :: SingleId) = Ref(s);

"""
	$(SIGNATURES)

Description of a `SingleId`.
"""
description(s :: SingleId) = s.description;

name(s :: SingleId) = s.name;
index(s :: SingleId) = s.index;

"""
	$(SIGNATURES)

Does `SingleId` have an index, indicating that it is part of a `Vector` of objects?
"""
has_index(this :: SingleId) = !Base.isempty(this.index)

function Base.:(==)(id1 :: SingleId, id2 :: SingleId)
    return (id1.name == id2.name)  &&  (id1.index == id2.index)
end

Base.isequal(id1 :: SingleId, id2 :: SingleId) =
    id1 == id2;

Base.:(==)(id1V :: Vector{SingleId},  id2V :: Vector{SingleId}) =
    all(isequal.(id1V, id2V));

Base.isequal(id1V :: Vector{SingleId}, id2V :: Vector{SingleId}) =
    all(isequal.(id1V, id2V));

# Must be implmented for storage in `Dict`.
Base.hash(id1 :: SingleId, h :: UInt) =
    hash(id1.name, hash(id1.index, hash(:SingleId, h)));

## ----------------  Show

# Make a string of the form "x[2, 1]"
# function show_string(s :: SingleId)
#     outStr = string(s.name);
#     if has_index(s)
#         outStr = outStr * "$(index(s))";
#     end
#     return outStr
# end

Base.show(io :: IO,  s :: SingleId) =
    print(io,  "SingleId:  $(make_string(s))");

"""
	$(SIGNATURES)

Make a string from a `SingleId`. Such as "x[2, 1]".
"""
function make_string(id :: SingleId)
    if !has_index(id)
        outStr = "$(name(id))"
    elseif length(id.index) == 1
        outStr = "$(name(id))$(index(id))"
    else
        outStr = "$(name(id))$(index(id))"
    end
    return outStr
end

"""
	$(SIGNATURES)

The inverse of [`make_string`](@ref).
"""
function make_single_id(s :: T1) where T1 <: AbstractString
    if occursin('[', s)
        # Pattern "id1[4, 3]"
        m = match(r"(.+)\[([0-9, ]+)+\]", s);
        idxV = parse.(Int, split(m[2], ","));
        sId = SingleId(Symbol(m[1]),  idxV);
    else
        sId = SingleId(Symbol(s));
    end
    return sId
end

# ----------------


text:
using IterTools

#functions to collapse multiple tree records into a single tree record with sum of expansions
Base.isless(lhs::Tree, rhs::Tree) = isless((lhs.diameter, lhs.height, lhs.species), (rhs.diameter, rhs.height, rhs.species))

function collapse(ta::Array{Tree,1})
    !issorted(ta) && sort(ta)
    #group by diameter,height,species-  needs to be sorted
    tgroups = groupby(s -> (s.diameter,s.height,s.species),sort(ta))
    # Extract the diameter,height,species from the group members
    trees = imap(s ->(s[1].diameter,s[1].height,s[1].species), tgroups)
    # Compute the sum of expansion
    expans = sum.(imap(x -> getfield.(x, :expansion), tgroups))
    #combine to new Array{Tree}
    return [Tree((t...,s)...) for (t,s) in zip(trees,expans)]
end

function collapse!(s::Stand)
    s.treelist = collapse(s.treelist)
end
#functions to split a single Tree instance into many depending on its expansion
function expand(t::Tree; min_exp = 1.0)
    ch = Channel() do c
        cnt = t.expansion
        while cnt > 0.0
            new_exp = cnt % min_exp != 0.0 ?  cnt % min_exp : min_exp
            put!(c, Tree(t.diameter, t.height, t.species, new_exp))
        cnt -= new_exp
        end
    end
    collect(ch)
end

function expand(ta::Array{Tree,1}; min_exp = 1.0)
    ch = Channel() do c
        for t in ta
        cnt = t.expansion
        while cnt > 0.0
            new_exp = cnt % min_exp != 0.0 ?  cnt % min_exp : min_exp
            put!(c, Tree(t.diameter, t.height, t.species, new_exp))
        cnt -= new_exp
        end
        end
    end
    collect(ch)
end

function expand!(s::Stand)
    s.treelist= expand(s.treelist)
end

function Base.:/(t1::Tree, denom::Number)
    expand(t1, min_exp =(t1.expansion/denom))
end


text:
# Autogenerated wrapper script for GnuPG_jll for i686-linux-gnu
export dirmngr, dirmngr_client, gpg, gpg_agent, gpg_connect_agent, gpgconf, gpgscm, gpgsm, gpgtar, gpgv, kbxutil

using GnuTLS_jll
using Libksba_jll
using Libgcrypt_jll
using Libgpg_error_jll
using nPth_jll
using Zlib_jll
using Libassuan_jll
using OpenLDAPClient_jll
using Bzip2_jll
using SQLite_jll
using libusb_jll
using Nettle_jll
JLLWrappers.@generate_wrapper_header("GnuPG")
JLLWrappers.@declare_executable_product(dirmngr)
JLLWrappers.@declare_executable_product(dirmngr_client)
JLLWrappers.@declare_executable_product(gpg)
JLLWrappers.@declare_executable_product(gpg_agent)
JLLWrappers.@declare_executable_product(gpg_connect_agent)
JLLWrappers.@declare_executable_product(gpgconf)
JLLWrappers.@declare_executable_product(gpgscm)
JLLWrappers.@declare_executable_product(gpgsm)
JLLWrappers.@declare_executable_product(gpgtar)
JLLWrappers.@declare_executable_product(gpgv)
JLLWrappers.@declare_executable_product(kbxutil)
function __init__()
    JLLWrappers.@generate_init_header(GnuTLS_jll, Libksba_jll, Libgcrypt_jll, Libgpg_error_jll, nPth_jll, Zlib_jll, Libassuan_jll, OpenLDAPClient_jll, Bzip2_jll, SQLite_jll, libusb_jll, Nettle_jll)
    JLLWrappers.@init_executable_product(
        dirmngr,
        "bin/dirmngr",
    )

    JLLWrappers.@init_executable_product(
        dirmngr_client,
        "bin/dirmngr-client",
    )

    JLLWrappers.@init_executable_product(
        gpg,
        "bin/gpg",
    )

    JLLWrappers.@init_executable_product(
        gpg_agent,
        "bin/gpg-agent",
    )

    JLLWrappers.@init_executable_product(
        gpg_connect_agent,
        "bin/gpg-connect-agent",
    )

    JLLWrappers.@init_executable_product(
        gpgconf,
        "bin/gpgconf",
    )

    JLLWrappers.@init_executable_product(
        gpgscm,
        "bin/gpgscm",
    )

    JLLWrappers.@init_executable_product(
        gpgsm,
        "bin/gpgsm",
    )

    JLLWrappers.@init_executable_product(
        gpgtar,
        "bin/gpgtar",
    )

    JLLWrappers.@init_executable_product(
        gpgv,
        "bin/gpgv",
    )

    JLLWrappers.@init_executable_product(
        kbxutil,
        "bin/kbxutil",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "memory allocator" begin

CuArrays.reclaim(1024)
CuArrays.reclaim()

@test CuArrays.@retry_reclaim(42, 42) == 42
@test CuArrays.@retry_reclaim(42, 41) == 41

@test_throws OutOfGPUMemoryError CuArray{Int}(undef, 10^20)

end


text:
using RegistryIndexer
using Test

@testset "RegistryIndexer.jl" begin
    # Write your tests here.
end


text:
# Gleb: the code taken from the ode elimination projects, tests are there

using Oscar

#---------Beginning of utils---------------------------------------------------

function eval_at_dict(poly::P, d::Dict{P, <: RingElem}) where P <: MPolyElem
    """
    Evaluates a polynomial on a dict var => val
    missing values are replaced with zeroes
    """
    point = [get(d, v, base_ring(parent(poly))(0)) for v in gens(parent(poly))]
    return evaluate(poly, point)
end

#------------------------------------------------------------------------------

function unpack_fraction(f::MPolyElem)
    return (f, one(parent(f)))
end

function unpack_fraction(f::Generic.Frac{<: MPolyElem})
    return (numerator(f), denominator(f))
end

#------------------------------------------------------------------------------

function str_to_var(s::String, ring::MPolyRing)
    ind = findfirst(v -> (string(v) == s), symbols(ring))
    if ind == nothing
        throw(Base.KeyError("Variable $s is not found in ring $ring"))
    end
    return gens(ring)[ind]
end

#------------------------------------------------------------------------------

function var_to_str(v::MPolyElem)
    ind = findfirst(vv -> vv == v, gens(parent(v)))
    return string(symbols(parent(v))[ind])
end

#------------------------------------------------------------------------------

function switch_ring(v::MPolyElem, ring::MPolyRing)
    """
    For a variable v, returns a variable in ring with the same name
    """
    ind = findfirst(vv -> vv == v, gens(parent(v)))
    return str_to_var(string(symbols(parent(v))[ind]), ring)
end

#---------End of utils---------------------------------------------------------
#------------------------------------------------------------------------------

function truncate_matrix(M::MatElem{<: Generic.AbsSeriesElem}, prec::Int)
    return map(e -> truncate(e, prec), M)
end

#------------------------------------------------------------------------------

function matrix_set_prec!(M::MatElem{<: Generic.AbsSeriesElem}, prec::Int)
    map(e -> set_prec!(e, prec), M)
end

#------------------------------------------------------------------------------

function ps_matrix_const_term(M::MatElem{<: Generic.AbsSeriesElem})
    return map(e -> coeff(e, 0), M)
end

#------------------------------------------------------------------------------

function _matrix_inv_newton_iteration(M::MatElem{T}, Minv::MatElem{T}) where T <: Generic.AbsSeriesElem{<: Generic.FieldElem}
    """
    Performs a single step of Newton iteration for inverting M with
    Minv being a partial result
    """
    return 2 * Minv - Minv * M * Minv
end

#--------------------------------------

function ps_matrix_inv(M::MatElem{<: Generic.AbsSeriesElem{<: Generic.FieldElem}}, prec::Int=-1)
    """
    Input:
        - M - a square matrix with entries in a univariate power series ring
          it is assumed that M(0) is invertible and all entries having the same precision
    Output: the inverse of M computed up to prec (the precision of the matrix if nothing)
    """
    const_term = ps_matrix_const_term(M)
    prec = (prec == -1) ? precision(M[1, 1]) : prec
    power_series_ring = base_ring(parent(M))
    result = map(a -> power_series_ring(a), Base.inv(const_term))
    cur_prec = 1
    while cur_prec < prec
        result = _matrix_inv_newton_iteration(M, result)
        cur_prec *= 2
    end
    return result
end

#------------------------------------------------------------------------------

function ps_diff(ps::Generic.AbsSeriesElem{<: Generic.RingElem})
    """
    Input: ps - (absolute capped) unvariate power series
    Output: the derivative of ps
    """
    result = zero(parent(ps))
    set_prec!(result, precision(ps))
    for exp in 1:(precision(ps) - 1)
        setcoeff!(result, exp - 1, coeff(ps, exp) * exp)
    end
    return result
end

#------------------------------------------------------------------------------

function ps_integrate(ps::Generic.AbsSeriesElem{<: Generic.FieldElem})
    """
    Input: ps - (absolute capped) unvariate power series
    Output: the integral of ps without constant term
    """
    result = zero(parent(ps))
    set_prec!(result, precision(ps) + 1)
    for exp in 0:(precision(ps) - 1)
        setcoeff!(result, exp + 1, coeff(ps, exp) // (exp + 1))
    end
    return result
end

#------------------------------------------------------------------------------

function _matrix_homlinear_de_newton_iteration(
        A::MatElem{<: Generic.AbsSeriesElem{T}},
        Y::MatElem{<: Generic.AbsSeriesElem{T}},
        Z::MatElem{<: Generic.AbsSeriesElem{T}},
        cur_prec::Int
    ) where T <: Generic.FieldElem
    Yprime = map(ps_diff, Y)
    Z = Z + truncate_matrix(Z * (one(parent(A)) - Y * Z), cur_prec)
    Y = Y - truncate_matrix(Y * map(ps_integrate, Z * (Yprime - truncate_matrix(A, 2 * cur_prec - 1) * Y)), 2 * cur_prec)
    return (Y, Z)
end

#--------------------------------------

function ps_matrix_homlinear_de(
        A::MatElem{<: Generic.AbsSeriesElem{T}},
        Y0::MatElem{<: T},
        prec::Int=-1
    ) where T <: Generic.FieldElem
    """
    Input:
        - A - a square matrix with entries in a univariate power series ring
        - Y0 - a square invertible matrix over the base field
    Output: matrix Y such that Y' = AY up to precision of A - 1 and Y(0) = Y0
    """
    prec = (prec == -1) ? precision(A[1, 1]) : prec
    ps_ring = base_ring(parent(A))
    cur_prec = 1
    Y = (one(parent(A)) + gen(ps_ring) * truncate_matrix(A, cur_prec)) * map(e -> truncate(ps_ring(e), cur_prec), Y0)
    Z = map(e -> truncate(ps_ring(e), cur_prec), Base.inv(Y0))
    while cur_prec < prec
        matrix_set_prec!(Y, 2 * cur_prec)
        matrix_set_prec!(Z, cur_prec)
        Y, Z = _matrix_homlinear_de_newton_iteration(A, Y, Z, cur_prec)
        cur_prec *= 2
    end
    return Y, Z
end

#------------------------------------------------------------------------------

function _variation_of_constants(
        A::MatElem{<: Generic.AbsSeriesElem{T}},
        B::MatElem{<: Generic.AbsSeriesElem{T}},
        Yh::MatElem{<: Generic.AbsSeriesElem{T}},
        Zh::MatElem{<: Generic.AbsSeriesElem{T}},
        Y0::MatElem{<: T},
        prec::Int
    ) where T <: Generic.FieldElem
    Zh += truncate_matrix(Zh * (one(parent(A)) - Yh * Zh), prec)
    Y_particular = truncate_matrix(Yh * map(ps_integrate, Zh * B), prec)
    return Y_particular + Yh * map(e -> parent(A[1, 1])(e), Y0)
end

#--------------------------------------

function ps_matrix_linear_de(
        A::MatElem{<: Generic.AbsSeriesElem{T}},
        B::MatElem{<: Generic.AbsSeriesElem{T}},
        Y0::MatElem{<: T},
        prec::Int=-1
    ) where T <: Generic.FieldElem
    """
    Input:
        - A, B - square matrices with entries in a univariate power series ring
        - Y0 - a matrix over the base field with the rows number the same as A
    Output: matrix Y such that Y' = AY + B up to precision of A - 1 and Y(0) = Y0
    """
    prec = (prec == -1) ? precision(A[1, 1]) : prec
    n = nrows(A)
    identity = one(MatrixSpace(base_ring(parent(Y0)), n, n))
    Yh, Zh = ps_matrix_homlinear_de(A, identity, prec)
    matrix_set_prec!(Zh, prec)
    return _variation_of_constants(A, B, Yh, Zh, Y0, prec)
end

#------------------------------------------------------------------------------

function ps_ode_solution(
        equations::Array{P, 1},
        ic::Dict{P, T},
        inputs::Dict{P, Array{T, 1}},
        prec::Int
    ) where {T <: Generic.FieldElem, P <: MPolyElem{T}}
    """
    Input
        - equations - a system of the form A(x, u, mu)x' - B(x, u, mu) = 0,
                      where A is a generically nonsingular square matrix
        - ic - initial conditions for x's (dictionary)
        - inputs - power series for inputs represented as arrays (dictionary)
        - prec - precision of the solution
        Assumption: A is nonzero at zero
    Output: power series solution of the system
    """
    n = length(equations)
    ring = parent(equations[1])
    S = MatrixSpace(ring, n, n)
    Sv = MatrixSpace(ring, n, 1)
    Svconst = MatrixSpace(base_ring(ring), n, 1)
    eqs = Sv(equations)

    x_vars = filter(v -> (var_to_str(v) * "_dot" in map(string, symbols(ring))), gens(ring))
    x_vars = [x for x in x_vars]
    x_dot_vars = [str_to_var(var_to_str(x) * "_dot", ring) for x in x_vars]

    Jac_dots = S([derivative(p, xd) for p in equations, xd in x_dot_vars])
    Jac_xs = S([derivative(p, x) for p in equations, x in x_vars])

    ps_ring, t = PowerSeriesRing(base_ring(ring), prec, "t"; model=:capped_absolute)
    solution = Dict()
    for (u, coeffs) in inputs
        solution[u] = sum([coeffs[i] * t^(i - 1) for i in 1:length(coeffs)])
    end
    for x in x_vars
        solution[x] = ps_ring(ic[x])
    end
    for xd in x_dot_vars
        solution[xd] = ps_ring(0)
        set_prec!(solution[xd], 1)
    end

    cur_prec = 1
    while cur_prec < prec
        for i in 1:length(x_vars)
            set_prec!(solution[x_vars[i]], 2 * cur_prec)
            set_prec!(solution[x_dot_vars[i]], 2 * cur_prec)
        end
        eval_point = [solution[v] for v in gens(ring)]
        map(ps -> set_prec!(ps, 2 * cur_prec), eval_point)
        eqs_eval = map(p -> evaluate(p, eval_point), eqs)
        J_eval = map(p -> evaluate(p, eval_point), Jac_xs)
        Jd_eval = map(p -> evaluate(p, eval_point), Jac_dots)
        Jd_inv = ps_matrix_inv(Jd_eval, cur_prec)

        X_err = ps_matrix_linear_de(-Jd_inv * J_eval, -Jd_inv * eqs_eval, zero(Svconst))
        for i in 1:length(x_vars)
            solution[x_vars[i]] = solution[x_vars[i]] + X_err[i]
            solution[x_dot_vars[i]] = ps_diff(solution[x_vars[i]])
        end
        cur_prec *= 2
    end

    return solution
end

#------------------------------------------------------------------------------

function ps_ode_solution(
        equations::Array{P, 1},
        ic::Dict{P, Int},
        inputs::Dict{P, Array{Int, 1}},
        prec::Int
    ) where P <: MPolyElem{<: Generic.FieldElem}
    bring = base_ring(parent(equations[1]))
    ps_ode_solution(
        equations,
        Dict(x => bring(v) for (x, v) in ic),
        Dict(u => map(v -> bring(v), vv) for (u, vv) in inputs),
        prec
    )
end

max_stars_repo_name: pogudingleb/ExperimentsBound

max_issues_repo_name: pogudingleb/ExperimentsBound

max_forks_repo_name: pogudingleb/ExperimentsBound


text:
module DummySepa

using SCIP

"""
A minimal no-op separator.
"""
mutable struct Sepa <: SCIP.AbstractSeparator
    called::Int64

    Sepa() = new(0)
end

function SCIP.exec_lp(sepa::Sepa)
    sepa.called += 1
    return SCIP.SCIP_DIDNOTRUN
end

end # module Dummy

module AddSingleCut

using SCIP

"""
A separator, that adds one predefined cut.
"""
mutable struct Sepa <: SCIP.AbstractSeparator
    scipd::SCIP.SCIPData
    varrefs::AbstractArray{SCIP.VarRef}
    coefs::AbstractArray{Float64}
    lhs::Float64
    rhs::Float64
end

function SCIP.exec_lp(sepa::Sepa)
    SCIP.add_cut_sepa(sepa.scipd.scip[], sepa.scipd.vars, sepa.scipd.sepas, sepa, sepa.varrefs, sepa.coefs,
                      sepa.lhs, sepa.rhs, removable=false)
    return SCIP.SCIP_SEPARATED
end

end # module SelectSol

"""
    sepa_set_scip_parameters(setter::Function)

Set all SCIP parameters, that are needed for the sepa-tests using `setter`.

# Arguments
- `setter(parameter::String, value::Int64)`
"""
function sepa_set_scip_parameters(setter::Function)
    setter("display/verblevel", 0)

    # All the following parameters make sure, that the separator is called and
    # the cut is added to the LP.

    setter("presolving/maxrounds", 0)
    setter("misc/usesymmetry", 0)

    setter("separating/minefficacy", 0)
    setter("separating/minefficacyroot", 0)
    setter("separating/minortho", 0)
    setter("separating/minorthoroot", 0)
    setter("separating/poolfreq", 1)

    # These parameters come from `SCIP> set sepa emph off`
    setter("separating/disjunctive/freq", -1)
    setter("separating/impliedbounds/freq", -1)
    setter("separating/gomory/freq", -1)
    setter("separating/strongcg/freq", -1)
    setter("separating/aggregation/freq", -1)
    setter("separating/clique/freq", -1)
    setter("separating/zerohalf/freq", -1)
    setter("separating/mcf/freq", -1)
    setter("separating/flowcover/freq", -1)
    setter("separating/cmir/freq", -1)
    setter("separating/rapidlearning/freq", -1)
    setter("constraints/cardinality/sepafreq", -1)
    setter("constraints/soc/sepafreq", -1)
    setter("constraints/SOS1/sepafreq", -1)
    setter("constraints/SOS2/sepafreq", -1)
    setter("constraints/varbound/sepafreq", -1)
    setter("constraints/knapsack/sepafreq", -1)
    setter("constraints/setppc/sepafreq", -1)
    setter("constraints/linking/sepafreq", -1)
    setter("constraints/or/sepafreq", -1)
    setter("constraints/and/sepafreq", -1)
    setter("constraints/xor/sepafreq", -1)
    setter("constraints/linear/sepafreq", -1)
    setter("constraints/orbisack/sepafreq", -1)
    setter("constraints/symresack/sepafreq", -1)
    setter("constraints/logicor/sepafreq", -1)
    setter("constraints/cumulative/sepafreq", -1)
    setter("constraints/abspower/sepafreq", -1)
    setter("constraints/bivariate/sepafreq", -1)
    setter("constraints/quadratic/sepafreq", -1)
    setter("constraints/nonlinear/sepafreq", -1)
    setter("constraints/indicator/sepafreq", -1)

    # These parameters come from `SCIP> set heuristics off`
    setter("heuristics/padm/freq", -1)
    setter("heuristics/ofins/freq", -1)
    setter("heuristics/trivialnegation/freq", -1)
    setter("heuristics/reoptsols/freq", -1)
    setter("heuristics/trivial/freq", -1)
    setter("heuristics/clique/freq", -1)
    setter("heuristics/locks/freq", -1)
    setter("heuristics/vbounds/freq", -1)
    setter("heuristics/shiftandpropagate/freq", -1)
    setter("heuristics/completesol/freq", -1)
    setter("heuristics/simplerounding/freq", -1)
    setter("heuristics/randrounding/freq", -1)
    setter("heuristics/zirounding/freq", -1)
    setter("heuristics/rounding/freq", -1)
    setter("heuristics/shifting/freq", -1)
    setter("heuristics/intshifting/freq", -1)
    setter("heuristics/oneopt/freq", -1)
    setter("heuristics/indicator/freq", -1)
    setter("heuristics/adaptivediving/freq", -1)
    setter("heuristics/farkasdiving/freq", -1)
    setter("heuristics/feaspump/freq", -1)
    setter("heuristics/conflictdiving/freq", -1)
    setter("heuristics/pscostdiving/freq", -1)
    setter("heuristics/fracdiving/freq", -1)
    setter("heuristics/nlpdiving/freq", -1)
    setter("heuristics/veclendiving/freq", -1)
    setter("heuristics/distributiondiving/freq", -1)
    setter("heuristics/objpscostdiving/freq", -1)
    setter("heuristics/rootsoldiving/freq", -1)
    setter("heuristics/linesearchdiving/freq", -1)
    setter("heuristics/guideddiving/freq", -1)
    setter("heuristics/rens/freq", -1)
    setter("heuristics/alns/freq", -1)
    setter("heuristics/rins/freq", -1)
    setter("heuristics/gins/freq", -1)
    setter("heuristics/lpface/freq", -1)
    setter("heuristics/crossover/freq", -1)
    setter("heuristics/undercover/freq", -1)
    setter("heuristics/subnlp/freq", -1)
    setter("heuristics/mpec/freq", -1)
    setter("heuristics/multistart/freq", -1)
    setter("heuristics/trysol/freq", -1)
end


text:
module HybridSystems

include("automata.jl")
include("switchings.jl")

using MathematicalSystems

export AbstractHybridSystem, HybridSystem

"""
    AbstractHybridSystem

Abstract supertype for a hybrid system.
"""
abstract type AbstractHybridSystem <: MathematicalSystems.AbstractSystem end

"""
    HybridSystem{A, S, R, W} <: AbstractHybridSystem

A hybrid system modelled as a hybrid automaton. The automaton `automaton` of
type `A` models the different discrete states and the allowed transitions with
corresponding labels.
The mode dynamic and domain are stored in a continuous dynamical system of type
`S` in the vector `modes`. They are indexed by the discrete states of the
automaton.
The reset maps and guards are given as discrete dynamical system or discrete map
of type `R` in the vector `resetmaps`. They are indexed by the labels of the
corresponding transition.
The switching of type `W` is given in the `switchings` vector,
indexed by the label of the transition.
Additional data can be stored in the `ext` field.

### Fields

- `automaton`  -- hybrid automaton of type `A`.
- `modes`      -- vector of modes of type `S` indexed by the discrete states,
                  both the domain and the dynamic are stored in this field.
                  See [`stateset`](@ref) to access the domain.
- `resetmaps`  -- vector of reset maps of type `R` indexed by the label of the
                  transition, the guard is stored as constraint of the map in
                  this field. See [`stateset`](@ref) to access the guard.
- `switchings` -- vector of switchings of type `W` indexed by the label of the
                  transition, see [`AbstractSwitching`](@ref).
- `ext`        -- dictionary that can be used by extensions.

## Examples

See [the Thermostat example](https://github.com/blegat/HybridSystems.jl/blob/master/examples/Thermostat.ipynb).
"""
struct HybridSystem{A, S, R, W} <: AbstractHybridSystem
    automaton::A
    modes::AbstractVector{S}
    resetmaps::AbstractVector{R}
    switchings::AbstractVector{W}
    # Can be used by extensions
    ext::Dict{Symbol, Any}
end
function HybridSystem(a, m, r, s)
    HybridSystem(a, m, r, s, Dict{Symbol, Any}())
end

if VERSION >= v"0.7-"
    Base.broadcastable(hs::HybridSystem) = Ref(hs)
end

# The hybrid system both acts like an automaton and a system

# Automaton
for f in [:state_property, :transition_property, :states, :nstates, :transitiontype, :transitions, :ntransitions, :source, :event, :target, :has_transition, :in_transitions, :out_transitions]
    @eval begin
        $f(hs::HybridSystem, args...) = $f(hs.automaton, args...)
    end
end
for f in [:state_property_type, :transition_property_type]
    @eval begin
        $f(::Type{<:HybridSystem{A}}, args...) where {A} = $f(A, args...)
    end
end

# System
export statedim, stateset, inputdim, inputset
"""
    statedim(hs::HybridSystem, u::Int)

Returns the dimension of the state space of the system at mode `u`.
"""
MathematicalSystems.statedim(hs::HybridSystem, u::Int) = statedim(hs.modes[u])

"""
    stateset(s::AbstractSystem, u::Int)

Returns the set of allowed states of the system at mode `u`.
"""
MathematicalSystems.stateset(hs::HybridSystem, u::Int) = stateset(hs.modes[u])

"""
    stateset(s::AbstractSystem, t)

Returns the guard for the transition `t`.
"""
MathematicalSystems.stateset(hs::HybridSystem, t) = stateset(hs.resetmaps[symbol(hs, t)])

"""
    inputdim(s::AbstractSystem, u::Int)

Returns the dimension of the input space of the system at mode `u`.
"""
MathematicalSystems.inputdim(hs::HybridSystem, u::Int) = inputdim(hs.modes[u])

"""
    inputset(s::AbstractSystem, u::Int)

Returns the set of allowed inputs of the system at mode `u`.
"""
MathematicalSystems.inputset(hs::HybridSystem, u::Int) = inputset(hs.modes[u])

"""
    inputset(s::AbstractSystem, t)

Returns the st of allowed inputs for the transition `t`.
"""
MathematicalSystems.inputset(hs::HybridSystem, t) = inputset(hs.resetmaps[symbol(hs, t)])

function Base.show(io::IO, hs::HybridSystem)
    print(io, "Hybrid System with automaton ")
    print(io, hs.automaton)
end

# Particular cases
include("switchedsystems.jl")

end # module


text:
"""
    SimpleNetHP

Hyperparameters for the simplenet architecture.

| Parameter                     | Description                                  |
|:------------------------------|:---------------------------------------------|
| `width :: Int`                | Number of neurons on each dense layer        |
| `depth_common :: Int`         | Number of dense layers in the trunk          |
| `depth_phead = 1`             | Number of hidden layers in the actions head  |
| `depth_vhead = 1`             | Number of hidden layers in the value  head   |
| `use_batch_norm = false`      | Use batch normalization between each layer   |
| `batch_norm_momentum = 0.6f0` | Momentum of batch norm statistics updates    |
"""
@kwdef struct SimpleNetHP
  width :: Int
  depth_common :: Int
  depth_phead :: Int = 1
  depth_vhead :: Int = 1
  use_batch_norm :: Bool = false
  batch_norm_momentum :: Float32 = 0.6f0
end

"""
    SimpleNet <: TwoHeadNetwork

A simple two-headed architecture with only dense layers.
"""
mutable struct SimpleNet <: TwoHeadNetwork
  gspec
  hyper
  common
  vhead
  phead
end

function SimpleNet(gspec::AbstractGameSpec, hyper::SimpleNetHP)
  bnmom = hyper.batch_norm_momentum
  function make_dense(indim, outdim)
    if hyper.use_batch_norm
      Chain(
        Dense(indim, outdim),
        BatchNorm(outdim, relu, momentum=bnmom))
    else
      Dense(indim, outdim, relu)
    end
  end
  indim = prod(GI.state_dim(gspec))
  outdim = GI.num_actions(gspec)
  hsize = hyper.width
  hlayers(depth) = [make_dense(hsize, hsize) for i in 1:depth]
  common = Chain(
    flatten,
    make_dense(indim, hsize),
    hlayers(hyper.depth_common)...)
  vhead = Chain(
    hlayers(hyper.depth_vhead)...,
    Dense(hsize, 1, tanh))
  phead = Chain(
    hlayers(hyper.depth_phead)...,
    Dense(hsize, outdim),
    softmax)
  SimpleNet(gspec, hyper, common, vhead, phead)
end

Network.HyperParams(::Type{SimpleNet}) = SimpleNetHP

function Base.copy(nn::SimpleNet)
  return SimpleNet(
    nn.gspec,
    nn.hyper,
    deepcopy(nn.common),
    deepcopy(nn.vhead),
    deepcopy(nn.phead)
  )
end


text:
function rule(line)
    outer, inner = split(line, r" bags contain ")
    list = Dict{String,Int}()
    for m in eachmatch(r"(\d+) (.*?) bag", inner)
        n, bag = m.captures
        list[bag] = parse(Int, n)
    end
    return outer => list
end

function part1(rules)
    queue = [k for (k, v) in rules if haskey(v, "shiny gold")]
    found = Set(queue)
    while !isempty(queue)
        bag = pop!(queue)
        new = (k for (k, v) in rules if haskey(v, bag))
        union!(found, new)
        append!(queue, new)
    end
    return length(found)
end

function part2(rules)
    queue = ["shiny gold" => 1]
    found = Pair{String,Int}[]
    while !isempty(queue)
        bag, n = pop!(queue)
        new = (k => n * v for (k, v) in rules[bag])
        append!(found, new)
        append!(queue, new)
    end
    return sum(last, found)
end

rules = Dict(rule(line) for line in eachline("./2020/data/day_7.txt"))

queue = ["shiny gold" => 1]
found = Pair{String,Int}[]
while !isempty(queue)
    lq = length(queue)
    lf = length(found)
    println("#Queue: $lq -- #QFound: $lf")
    bag, n = pop!(queue)
    new = (k => n * v for (k, v) in rules[bag])
    println("$rules[bag]")
    append!(found, new)
    append!(queue, new)
end
println(sum(last, found))

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
# phonopy_project.jl
# Demonstration / scratchpad for using JuliaPhonons module

push!(LOAD_PATH,"../src") # Temporary versions of modules in PWD
using JuliaPhonons

# The REPO comes with a basic MAPI cubic unit cell Phonopy calculation
P=read_POSCAR(open("../2015-11_CubicModeDecomposition/MAPbI/POSCAR"),expansion=[2,2,2])
eigenvectors,eigenmodes=read_meshyaml(open("../2015-11_CubicModeDecomposition/MAPbI/mesh.yaml"),P)

myf=open("mode_decomposition.dat","w")
gnuplot_header(P,f=myf)

mode=16

println("My eigenvector $mode, ")
show(eigenvectors[mode])
println("My eigenvalue for $mode, ")
show(eigenmodes[mode])
println()

# See: https://stackoverflow.com/a/11319919
# Essentially gen original + perturbed structure; best Quaternions minimises error on forwards / backwards transition
#

# Eeek! Hard coded.
println("C: ",eigenvectors[mode][1,:])
println("N: ",eigenvectors[mode][2,:])

end

for (count,(eigenvector,eigenmode)) in enumerate(zip(eigenvectors,eigenmodes))
    output_animated_xyz(P,count,eigenvector,eigenmode) # generates files anim_{count}.xyz

    decompose_eigenmode_atomtype(P,count,eigenvector,eigenmode)
    # Output again, to a file for later printing
    decompose_eigenmode_atomtype(P,count,eigenvector,eigenmode,f=myf)

    decompose_eigenmode_atom_contributions(P,count,eigenvector)
end

close(myf)

max_stars_repo_name: jarvist/Julia-Phonons

max_issues_repo_name: lucydot/Julia-Phonons

max_forks_repo_name: lucydot/Julia-Phonons


text:
"""A node which fundamentally represents a constant available over all time."""
struct Constant{T} <: NodeOp{T}
    value::T
end

Base.hash(x::Constant, h::UInt64) = hash(x.value, hash(:Constant, h))
function Base.:(==)(x::Constant{T}, y::Constant{T}) where {T}
    # Add short-circuit, in case Base.:(==)(::T, ::T) doesn't have one.
    x.value === y.value && return true
    return x.value == y.value
end

function Base.show(io::IO, op::Constant)
    return print(io, "$(typeof(op).name.name){$(value_type(op))}($(op.value))")
end

mutable struct ConstantState <: NodeEvaluationState
    ticked::Bool
    ConstantState() = new(false)
end

create_evaluation_state(::Tuple{}, ::Constant) = ConstantState()

_is_constant(::NodeOp) = false
_is_constant(::Constant) = true
_is_constant(node::Node) = _is_constant(node.op)

value(op::Constant) = op.value
value(node::Node) = value(node.op)

"""Identity if the argument is a node, otherwise wrap it in a constant node."""
_ensure_node(node::Node) = node
_ensure_node(value::Any) = obtain_node((), Constant(value))

function run_node!(
    op::Constant{T},
    state::ConstantState,
    time_start::DateTime,
    ::DateTime,  # time_end
) where {T}
    # The convention is that we emit a single knot at the start of the evaluation interval.
    return if state.ticked
        Block{T}()
    else
        state.ticked = true
        Block(unchecked, [time_start], [op.value])
    end
end

"""
    constant(value) -> Node

Explicitly wrap `value` into a `TimeDag` constant node, regardless of its type.

In many cases this isn't required, since many `TimeDag` functions which expect nodes will
automatically wrap non-node arguments into a constant node.

!!! warning
    If `value` is already a node, this will wrap it up in an additional node. This is very
    likely not what you want to do.
"""
constant(value::T) where {T} = obtain_node((), Constant(value))


text:
function showresults(sim::FastSim,head::Bool)
    net=sim.net
    nodes=nodecounts_f(net)
    links=linkcounts_f(net)
    push!(sim.results,[sim.t,nodes...,links...])
    ndig=sim.nodedigits
    ldig=sim.linkdigits
    tdig=sim.timedigits
    tdec=sim.timedec
    t=sim.t
    c=net.C-1
    l=net.L-2
    if head
        fastprint(sim,padright("Time",tdig))
        for i=1:c
            fastprint(sim,"    ",padright(net.nodealias[i],ndig))
        end
        for i=1:l
            fastprint(sim,"    ",padright(net.linkalias[i],ldig))
        end
        fastprint(sim,"\n")
    end
    fastprint(sim,showfloat(t,tdec,tdig))
    for i=1:c
        fastprint(sim,"    ",padleft(string(nodes[i]),ndig))
    end
    for i=1:l
        fastprint(sim,"    ",padleft(string(links[i]),ldig))
    end
    fastprint(sim,"\n")
    if length(sim.filename)!==0
        fl=open(sim.filename,"a")
        if sim.Nout===0
            print(fl,"Time")
            for i=1:c
                print(fl,",",net.nodealias[i])
            end
            for i=1:l
                print(fl,",",net.linkalias[i])
            end
            print(fl,"\n")
        end
        print(fl,t)
        for i=1:c
            print(fl,",",nodes[i])
        end
        for i=1:l
            print(fl,",",links[i])
        end
        print(fl,"\n")
        close(fl)
    end
end

"""
    simstep!(sim)

Simulate the next event in *sim*.

This function will always advance the sim by exactly one event.
If no event is possible it will advance time by one timeunit.
Output is generated at start time and directly after the event has occured.

See also [FastSim](#Fastnet.FastSim),[simstep](#Fastnet.simstep)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(2,1,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,2);

julia> function rates!(r,t)
         r[1]=countnodes(net,2)
         end
rates! (generic function with 1 method)

julia> function simpleproc!()
         node=randomnode(net,2)
         nodestate!(net,node,1)
         end
simpleproc! (generic function with 1 method)

julia> sim=FastSim(net,rates!,[simpleproc!]; output=false)
Simulation run, currently at time 0.0

julia> simstep!(sim)

julia> nodecounts(net)
2-element Vector{Int64}:
 2
 0
```
"""
function simstep!(sim::FastSim)
    rng=sim.net.rng
    rates=MVector{sim.Nproc}(zeros(sim.Nproc))
    sim.ratefunc(rates,sim.t)
    tot=sum(rates)
    if tot>0.0
         del=randexp(rng)/tot
    else
         del=1
    end
    initsim(sim,del,del)
    sim.repfunc(sim,true)
    sim.t+=del
    if tot>0.0
         r=rand(rng,EPS:tot)-rates[1]
         i=1
         while(r>0.0 && i<sim.Nproc)
            i+=1
            r-=rates[i]
        end
        sim.procfunc[i]()
    end
    sim.repfunc(sim,false)
end

"""
    runsim!(sim,dur,out=1.0)

Run the simulation *sim* for time *dur*, producing out at intervals *out*.

During the simulation the FastSim and its associated FastNet object will
be updated to reflect the current state of the network (though see notes
on the simulation time, [here](#Fastnet.FastSim)).

This function simulates the FastSim **at least** for a certain time. If there
are still events occuring in the simulation by the end of the simulation run
the simulation will stop directly after the first event that happens after *dur*.
So the simulation time will always be greater than *dur*. In general the difference and
the actual simulation time will be tiny, but in case events are extreley rare the
simulation may run significantly beyond *dur*. This behaviour is necessary to
avoid a watchdog-paradox artifact when repeatedly starting short runs.

Output is generated once at the start of the simulation and then at every multiple of out.
For example, if the simulation time is *t=3.12* at the start of the run, *dur=10* and *out=5* then
outputs will be generated at times 3.12, 5 and 10. As a result the network will
be left in a state that differs from the statistics in the last output.

See also [FastSim](#Fastnet.FastSim),[simstep!](#Fastnet.simstep!)

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(2,1,2,[])
Network of 0 nodes and 0 links

julia> n1=makenode!(net,1);

julia> n2=makenode!(net,2);

julia> function rates!(r,t)
         r[1]=countnodes(net,2)
         end
rates! (generic function with 1 method)

julia> function simpleproc!()
         node=randomnode(net,2)
         nodestate!(net,node,1)
         end
simpleproc! (generic function with 1 method)

julia> sim=FastSim(net,rates!,[simpleproc!]; output=false)
Simulation run, currently at time 0.0

julia> runsim!(sim,100)

julia> nodecounts(net)
2-element Vector{Int64}:
 2
 0
```
"""
function runsim!(sim::FastSim,dur,out=1.0)
    Tdur=Float64(dur)
    Tout=Float64(out)
    initsim(sim,Tdur,Tout)
    net=sim.net
    sim.repfunc(sim,true)
    Tend=sim.Tend
    t=sim.t
    Nproc=sim.Nproc
    rates=MVector{Nproc}(zeros(Nproc))
    nextout=(floor(t/Tout)+1.0)*Tout
    while(t<Tend)
        sim.ratefunc(rates,t)
        tot=sum(rates)
        if tot==0.0
            t=nextout
            sim.Nout+=1
            sim.t=nextout
            sim.repfunc(sim,false)
            nextout+=Tout
        else
            t+=randexp(net.rng)/tot
            while(nextout<t)
                sim.Nout+=1
                sim.t=nextout
                sim.repfunc(sim,false)
                nextout+=Tout
            end
            r=rand(net.rng,EPS:tot)-rates[1]
            i=1
            while(r>0 && i<Nproc)
                i+=1
                r-=rates[i]
            end
            sim.t=t
            sim.procfunc[i]()
        end
    end
    if sim.printresults==true
        display(quickplotnodes(sim))
    end
    nothing
end


text:
@registerN Hahn AbstractCDOP3 α β 𝐍
export Hahn

"""
    Hahn{α,β,𝐍}

References: [Koekoek and Swarttouw §1.5](https://arxiv.org/pdf/math/9602214.pdf)

!!! note
    In  [Koekoek and Swarttouw](https://arxiv.org/pdf/math/9602214.pdf) sections 1.3, 1.4, and 1.6  are other  Hahn-type polynomials,  not  implemented here.

"""
Hahn

abcde(::Type{<:Hahn{α,β,𝐍}})  where {α,β,𝐍} = NamedTuple{(:a,:b,:c,:d,:e)}((1, -(β+𝐍+1), 0,  α+β+2, -𝐍*(α+1)))

basis_symbol(::Type{<:Hahn{α,β,𝐍}}) where {α,β,𝐍} = "Q⁽ᵅᵝ⁾"
Polynomials.domain(::Type{<:Hahn{α, β, 𝐍}}) where {α, β,  𝐍} = Polynomials.Interval(0, 𝐍)
weight_function(::Type{<:Hahn{α,β,𝐍}}) where {α,β,𝐍} = x -> generalized_binomial(α+x,x) * generalized_binomial(𝐍+β-x,𝐍-x)

# use   recurrence relation  1.5.3 of Koekoek and Swarttouw
# -xQᵢ = AQᵢ₊₁ -  (A+C)Qᵢ + CQᵢ₋₁
#  or Qᵢ₊₁ = (1/A ⋅ x - (1 + C/A)) Qᵢ + C/A ⋅ Qᵢ₋₁
# Using kᵢ = A⁻¹ᵢ₋₁ ⋅ A⁻¹ᵢ₋₂ ⋯ A⁻¹₀; so kᵢ₊₁/kᵢ = A⁻¹ᵢ, kᵢ₊₁/kᵢ₋ᵢ = A⁻¹ᵢA⁻¹ᵢ₋₁
#
# function Aᵢ⁻¹(n,α,β,𝐍)
#     num = (2n+α+β + 1) * (2n+α+β+2)
#     den = (n+α + β + 1)*(n + α + 1) *(𝐍-n)
#     num/den
# end

# define kn, k1k_1 through defaults
function k1k0(P::Type{<:Hahn{α,β, 𝐍}}, n::Int) where {α,β,𝐍}

    num  = (2*n + α + β + 1)*(2*n + α + β + 2)
    den = (-n + 𝐍)*(n + α + 1)*(n + α + β + 1)

    -(one(eltype(P))*num)/den

end

# end

function classical_hypergeometric(::Type{<:Hahn{α,β,𝐍}}, n::Int, x)  where {α,β,𝐍}
    pFq((-n, -x, n+1+α+β),  (α+1, -𝐍),  1)
end


text:
using HCEstimator
using JuMP, Ipopt
using DataFrames, CSV

data = DataFrame(CSV.File("case33.csv"))

sub = DistSystem.Substation(
    12660,              # Nominal voltage (V)
    1,                  # Bus
    1.0,                # Vˢᴮ: Voltage (p.u.)
    4.0,                # Pˢᴮ: Active power capacity (MW)
    2.5,                # Qˢᴮ: Reactive power capacity (MVAr)
    [0.003, 12, 240]    # Costs (Not used)
    )

sys = DistSystem.factory_system(
    data,               # Dataframe of the distribution system
    0.93,               # V_low: Voltage lower bound (p.u.)
    1.05,               # V_upper: Voltage upper bound (p.u.)
    sub                 # DistSystem.Sustation
    )

der1 = DistSystem.DER(
    18,                  # Bus
    0.02,                # Sᴰᴱᴿ (MVA)
    1.00,                # Eᴰᴱᴿ (MWh). For DG's make Eᴰᴱᴿ a big number
    0.05,                # αᴰᴱᴿ
    0.05,                # βᴰᴱᴿ
    [0.0, 0.03],         # [Pᴰᴱᴿ_low, Pᴰᴱᴿ_upper]  (MW)
    [-0.03, 0.03],       # [Qᴰᴱᴿ_low, Qᴰᴱᴿ_upper]  (MVAr)
    [0.0, 1.0],          # Possible DER's Operation Scenarios ≠ μᴰᴱᴿ
    [0.0026, 10.26, 210] # Costs (Not used)
) #DG dispatchable

sys.dgs = [
    der1,
    DistSystem.DER(22, 0.02, 0.05, 0.0,  0.0,[0.0, 0.0], [-0.03, 0.03], [0.0, 1.0], [0.0026, 10.26, 210]), #DG non-active dispatchable
    DistSystem.DER(33, 0.02, 0.05, 0.05, 0.05,[-0.01, 0.01], [-0.01, 0.01], [-1.0, 0.0, 1.0], [0.0026, 10.26, 210]), #ESS
    DistSystem.DER(25, 0.02, 0.0, 0.00, 0.0,[0.0, 0.0], [0.0, 0.0], [0.0, 1.0], [0.0026, 10.26, 210]), #DG non-dispatchable
]

sys.m_load = [0.5, 0.8, 1.0] # μᴸ
sys.m_new_dg = [-1.0, 0.0, 1] # μᴴᶜ

sys.time_curr = 1.0 # Tᶜᵘʳ

model = build_model(Model(Ipopt.Optimizer), sys)
set_silent(model)

optimize!(model)
println("Hosting Capacity: ", round(objective_value(model), digits = 3), " MVA")

(Ω, bΩ, L, K, D, S) = Tools.Get.sets(sys)

## Voltage magnitude and power injection
dims_bus = Tuple(length(set) for set in (Ω, L, K, S))
V = zeros(Float64, dims_bus)
P = zeros(Float64, dims_bus)
Q = zeros(Float64, dims_bus)

## Power of DERs installed used by DisCO
dims_der = Tuple(length(set) for set in (D, L, K, S))
Pᴰᴱᴿ = zeros(Float64, dims_der)
Qᴰᴱᴿ = zeros(Float64, dims_der)

for b = Ω, l = L, k = K, s = S, d = D
    V[b, l, k, s] = Tools.Get.voltage_module(model, b, l, k, s)
    P[b, l, k, s] = Tools.Get.power_active(model, b, l, k, s)
    Q[b, l, k, s] = Tools.Get.power_reactive(model, b, l, k, s)
    Pᴰᴱᴿ[d, l, k, s] = Tools.Get.power_active_DER(model, d, l, k, s)
    Qᴰᴱᴿ[d, l, k, s] = Tools.Get.power_reactive_DER(model, d, l, k, s)
end


text:
"""
    NormalInverseGaussian(μ,α,β,δ)

The *Normal-inverse Gaussian distribution* with location `μ`, tail heaviness `α`, asymmetry parameter `β` and scale `δ` has probability density function

```math
f(x; \\mu, \\alpha, \\beta, \\delta) = \\frac{\\alpha\\delta K_1 \\left(\\alpha\\sqrt{\\delta^2 + (x - \\mu)^2}\\right)}{\\pi \\sqrt{\\delta^2 + (x - \\mu)^2}} \\; e^{\\delta \\gamma + \\beta (x - \\mu)}
```
where ``K_j`` denotes a modified Bessel function of the third kind.

External links

* [Normal-inverse Gaussian distribution on Wikipedia](http://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution)

"""
struct NormalInverseGaussian{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    α::T
    β::T
    δ::T
    γ::T
    function NormalInverseGaussian{T}(μ::T, α::T, β::T, δ::T) where T
        γ = sqrt(α^2 - β^2)

        new{T}(μ, α, β, δ, γ)
    end
end

NormalInverseGaussian(μ::T, α::T, β::T, δ::T) where {T<:Real} = NormalInverseGaussian{T}(μ, α, β, δ)
NormalInverseGaussian(μ::Real, α::Real, β::Real, δ::Real) = NormalInverseGaussian(promote(μ, α, β, δ)...)
function NormalInverseGaussian(μ::Integer, α::Integer, β::Integer, δ::Integer)
    return NormalInverseGaussian(float(μ), float(α), float(β), float(δ))
end

@distr_support NormalInverseGaussian -Inf Inf

#### Conversions
function convert(::Type{NormalInverseGaussian{T}}, μ::Real, α::Real, β::Real, δ::Real) where T<:Real
    NormalInverseGaussian(T(μ), T(α), T(β), T(δ))
end
function convert(::Type{NormalInverseGaussian{T}}, d::NormalInverseGaussian{S}) where {T <: Real, S <: Real}
    NormalInverseGaussian(T(d.μ), T(d.α), T(d.β), T(d.δ))
end

params(d::NormalInverseGaussian) = (d.μ, d.α, d.β, d.δ)
@inline partype(d::NormalInverseGaussian{T}) where {T<:Real} = T

mean(d::NormalInverseGaussian) = d.μ + d.δ * d.β / d.γ
var(d::NormalInverseGaussian) = d.δ * d.α^2 / d.γ^3
skewness(d::NormalInverseGaussian) = 3d.β / (d.α * sqrt(d.δ * d.γ))
kurtosis(d::NormalInverseGaussian) = 3 * (1 + 4*d.β^2/d.α^2) / (d.δ * d.γ)

function pdf(d::NormalInverseGaussian, x::Real)
    μ, α, β, δ = params(d)
    α * δ * besselk(1, α*sqrt(δ^2+(x - μ)^2)) / (π*sqrt(δ^2 + (x - μ)^2)) * exp(δ * d.γ + β*(x - μ))
end

function logpdf(d::NormalInverseGaussian, x::Real)
    μ, α, β, δ = params(d)
    log(α*δ) + log(besselk(1, α*sqrt(δ^2+(x-μ)^2))) - log(π*sqrt(δ^2+(x-μ)^2)) + δ*d.γ + β*(x-μ)
end

#### Sampling

# The Normal Inverse Gaussian distribution is a normal variance-mean
# mixture with an inverse Gaussian as mixing distribution.
#
# Ole E. Barndorff-Nielsen (1997)
# Normal Inverse Gaussian Distributions and Stochastic Volatility Modelling
# Scandinavian Journal of Statistics, Vol. 24, pp. 1--13
# DOI: http://dx.doi.org/10.1111/1467-9469.00045

function rand(rng::Random.AbstractRNG, d::NormalInverseGaussian)
    μ, α, β, δ = params(d)

    Z = InverseGaussian(δ/d.γ, δ^2)
    z = rand(rng, Z)
    X = Normal(μ + β*z, sqrt(z))
    return rand(rng, X)
end


text:
using YAXArrays, YAXArrayBase, Test, Dates

@testset "YAXArrays" begin
    data = collect(reshape(1:20,4,5))
    axlist = [RangeAxis("XVals",1.0:4.0), CategoricalAxis("YVals",[1,2,3,4,5])]
    props = Dict("att1"=>5, "att2"=>"Hallo")
    a = YAXArray(axlist, data, props)

    @testset "Array construction sanity checks" begin
        @test_throws ArgumentError YAXArray(axlist[1:1], data, props)
        @test_throws ArgumentError YAXArray([RangeAxis("XVals",1.0:4.0), CategoricalAxis("YVals",[1,2,3,4])],data,props)
    end

    @testset "Basic array Functions" begin
        @test size(a) == (4,5)
        @test size(a,1) == 4
        @test size(a,2) == 5
        @test eltype(a) == Int
        @test ndims(a) == 2
        pa = permutedims(a,(2,1))
        @test pa.axes == axlist[[2,1]]
        @test pa.data == transpose(data)
        @test a[2,3] == a.data[2,3]
        @test read(a) isa Array
        @test read(a) == collect(reshape(1:20,4,5))
    end

    @testset "YAXArray interface functions" begin

        a2 = readcubedata(a)
        @test a2.axes == a.axes
        @test a.data == a2.data
        @test a2.data isa Array

        @test caxes(a) == axlist

        @test YAXArrays.Cubes.cubechunks(a) == (4,5)

        cs = CartesianIndices.(collect(Iterators.product([1:2,3:5,6:6],[1:2,3:4])))
        @test YAXArrays.Cubes.common_size(cs) == (3,2)

        @test YAXArrays.Cubes.chunkoffset(a) == (0,0)

        renameaxis!(a2,"XVals"=>"Ax1")

        @test YAXArrayBase.dimname(a2,1) == :Ax1

        renameaxis!(a2,"Ax1"=>RangeAxis("XVals", 2:5))

        @test a2.axes[1].values == 2:5

        b = IOBuffer()
        show(b,a)
        # @test String(take!(b))=="""
        # YAXArray with the following dimensions
        # XVals               Axis with 4 Elements from 1.0 to 4.0
        # YVals               Axis with 5 elements: 1 2 3 4 5
        # Total size: 160.0 bytes
        # """
    end

    @testset "YAXArrayBase interface" begin
        using YAXArrayBase
        @test YAXArrayBase.dimvals(a,1) == 1.0:4.0
        @test YAXArrayBase.dimvals(a,2) == [1,2,3,4,5]
        @test YAXArrayBase.dimname(a,1) == :XVals
        @test YAXArrayBase.dimname(a,2) == :YVals
        @test YAXArrayBase.getattributes(a) == Dict("att1"=>5, "att2"=>"Hallo")
        @test YAXArrayBase.iscontdim(a,1) == true
        @test YAXArrayBase.iscontdim(a,2) == false
        @test YAXArrayBase.getdata(a) === data
        a2 = YAXArrayBase.yaxcreate(YAXArray, reshape(1:27,3,3,3),["d1","d2","d3"],
            [1:3,["One","Two","Three"], 0.1:0.1:0.3], Dict("att"=>5))
        @test a2 isa YAXArray
        @test a2.axes[1] == RangeAxis("d1",1:3)
        @test a2.axes[2] == CategoricalAxis("d2",["One","Two","Three"])
        @test a2.axes[3] == RangeAxis("d3",0.1:0.1:0.3)
        @test a2.data == reshape(1:27,3,3,3)
        @test a2.properties == Dict("att"=>5)
        @test YAXArrayBase.iscompressed(a) == false
    end

    @testset "Subsets" begin
        s = YAXArrays.Cubes.subsetcube(a,X=1.5..3.5)
        @test s.data == [2 6 10 14 18; 3 7 11 15 19]
        @test s.axes[1] == RangeAxis("XVals",2.0:3.0)
        @test s.axes[2] == CategoricalAxis("YVals",[1,2,3,4,5])
        ax = a.axes[1]
        @test YAXArrays.Cubes.interpretsubset(CartesianIndices((1:2,)), ax) == 1:2
        @test YAXArrays.Cubes.interpretsubset(CartesianIndex((2,)), ax) == 2
        @test YAXArrays.Cubes.interpretsubset(2.1, ax) == 2
        @test YAXArrays.Cubes.interpretsubset((3.5,1.5), ax) == 2:3
        @test YAXArrays.Cubes.interpretsubset(0.8..2.2, ax) == 1:2
        tax = RangeAxis("ADate",Date(2001):Day(1):Date(2003,2,28))
        @test YAXArrays.Cubes.interpretsubset((Date(2001,1,2),Date(2001,1,5)), tax) == 2:4
        @test YAXArrays.Cubes.interpretsubset(2001:2002, tax) == 1:730
        @test YAXArrays.Cubes.interpretsubset([1,3,5], a.axes[2]) == [1,3,5]

        s2 = a[X=0.5..3.5, Y=[1,5,4]]
        @test s2.data == [1 17 13; 2 18 14; 3 19 15]
        @test s2.axes[1] == RangeAxis("XVals",1.0:3.0)
        @test s2.axes[2] == CategoricalAxis("YVals",[1,5,4])
    end

end


text:
const MultiobjCoverProblemFitnessScheme =
    ParetoFitnessScheme{2,Float64,true,MultiobjCoverProblemScoreAggregator{MultiobjProblemSoftFold2d}}

"""
Parameters for the Borg-based optimization of [`MultiobjCoverProblem`](@ref).

See [`optimize`](@ref).
"""
struct MultiobjOptimizerParams <: AbstractOptimizerParams{MultiobjCoverProblem}
    pop_size::Int
    weight_digits::Union{Int, Nothing}
    max_steps::Int
    max_steps_without_progress::Int
    fitness_tolerance::Float64
    min_delta_fitness_tolerance::Float64
    fold_ratio_threshold::Float64
    #workers::Vector{Int}
    nworkers::Int
    trace_interval::Float64

    borg_params::BlackBoxOptim.ParamsDict

    function MultiobjOptimizerParams(;
        # default Borg/BBO Opt.Controller parameter overrides
        PopulationSize::Integer = 100,
        WeightDigits::Union{Integer, Nothing} = 2, # precision digits for set weights
        MaxSteps::Integer = 10_000_000,
        MaxStepsWithoutProgress::Integer = 50_000,
        FitnessTolerance::Real = 0.1,
        MinDeltaFitnessTolerance::Real = 1E-5,
        FoldRatioThreshold::Real = 1.0,
        NWorkers::Integer = max(1, Threads.nthreads() - 1), #Workers::AbstractVector{Int} = Vector{Int}(),
        TraceInterval::Real = 5.0,
        kwargs...
    )
        if #=isempty(Workers) &&=# NWorkers > max(1, Threads.nthreads()-1)
            @warn("Requested NWorkers=$NWorkers, while only $(Threads.nthreads()) worker thread(s) available, reducing")
            NWorkers = Threads.nthreads()-1
        end
        #if isempty(Workers) && NWorkers > 1
        #    Workers = workers()[1:NWorkers]
        #end
        new(PopulationSize, WeightDigits,
            MaxSteps, MaxStepsWithoutProgress,
            FitnessTolerance, MinDeltaFitnessTolerance, FoldRatioThreshold,
            NWorkers, TraceInterval,
            BlackBoxOptim.kwargs2dict(kwargs...))
    end
end

function borg_params(opt_params::MultiobjOptimizerParams,
                     problem::MultiobjCoverProblem)
    if haskey(opt_params.borg_params, :ϵ)
        eps = opt_params.borg_params[:ϵ]
    else
        #eps_scale = get(opt_params.borg_params, :ϵ_scale, [0.2, 0.1, 0.1])
        #eps = eps_scale .* [score_scales(problem)...]
        eps = 0.01
    end
    BlackBoxOptim.chain(BlackBoxOptim.BorgMOEA_DefaultParameters,
        ParamsDict(:PopulationSize=>opt_params.pop_size,
                   :NThreads=>opt_params.nworkers > 1 ? opt_params.nworkers : 0,
                   :ϵ=>eps),
        opt_params.borg_params)
end

bbo_ctrl_params(opt_params::MultiobjOptimizerParams) =
    BlackBoxOptim.chain(BlackBoxOptim.DefaultParameters,
        ParamsDict(:PopulationSize=>opt_params.pop_size,
                   :MaxSteps=>opt_params.max_steps,
                   :MaxStepsWithoutProgress=>opt_params.max_steps_without_progress,
                   :FitnessTolerance=>opt_params.fitness_tolerance,
                   :MinDeltaFitnessTolerance=>opt_params.min_delta_fitness_tolerance,
                   :TraceInterval=>opt_params.trace_interval))

genop_params(opt_params::MultiobjOptimizerParams) =
    BlackBoxOptim.ParamsDict()

"""
Wraps [`MultiobjCoverProblem`](@ref) as [`BlackBoxOptim.OptimizationProblem`]
"""
struct MultiobjCoverProblemBBOWrapper{P <: MultiobjCoverProblem, SS <: RectSearchSpace} <:
        BlackBoxOptim.OptimizationProblem{MultiobjCoverProblemFitnessScheme}
    orig::P
    fitness_scheme::MultiobjCoverProblemFitnessScheme

    search_space::SS

    function MultiobjCoverProblemBBOWrapper(
        orig::MultiobjCoverProblem,
        score_folding::DetailedScoreFolding,
        dimdigits::Union{Integer,Nothing}=2
    )
        fitness_scheme = MultiobjCoverProblemFitnessScheme(
            aggregator=MultiobjCoverProblemScoreAggregator(orig.params, score_folding))
        ss = RectSearchSpace(nvars(orig), (0.0, 1.0), dimdigits=dimdigits)
        new{typeof(orig), typeof(ss)}(orig, fitness_scheme, ss)
    end
end

Base.copy(problem::MultiobjCoverProblemBBOWrapper) =
    MultiobjCoverProblemBBOWrapper(problem.orig, dimdigits(search_space(problem), 1))
#=
BlackBoxOptim.show_fitness(io::IO, score::DetailedScore,
                           problem::MultiobjCoverProblemBBOWrapper{MultiobjProblemNoFolding}) =
    @printf(io, "(sets=%.3f set×set=%.3f uncovered=%.3f covered=%.3f)\n",
            score[1], score[2], score[3], score[4])

BlackBoxOptim.show_fitness(io::IO, score::NTuple{2,Float64},
                           problem::MultiobjCoverProblemBBOWrapper{MultiobjProblemDrop3}) =
    @printf(io, "(sets=%.3f set×set=%.3f)\n", score[1], score[2]^2)

BlackBoxOptim.show_fitness(io::IO, score::NTuple{3,Float64},
                           problem::MultiobjCoverProblemBBOWrapper{MultiobjProblemFold3d}) =
    @printf(io, "(sets+k(sets²+uncov)=%.3f sets²+k(sets+uncov)=%.3f uncov+k(sets+sets²)=%.3f)\n",
          score[1], score[2], score[3])
=#

BlackBoxOptim.show_fitness(io::IO, score::FoldedScore,
                           problem::MultiobjCoverProblemBBOWrapper) =
    @printf(io, "(sets+k·miscover=%.3f sets²=%.3f)\n", score[1], score[2])

BlackBoxOptim.show_fitness(io::IO, score::IndexedTupleFitness, problem::MultiobjCoverProblemBBOWrapper) =
    BlackBoxOptim.show_fitness(io, score.orig, problem)

fold_score(score::DetailedScore, p::MultiobjCoverProblemBBOWrapper) =
    p.fitness_scheme.aggregator.score_folding(score)

BlackBoxOptim.fitness(x::BlackBoxOptim.AbstractIndividual,
                      p::MultiobjCoverProblemBBOWrapper) =
    fold_score(score(x, p.orig), p)

generate_recombinators(problem::MultiobjCoverProblemBBOWrapper, params) =
    CrossoverOperator[BlackBoxOptim.DiffEvoRandBin1(chain(BlackBoxOptim.DE_DefaultOptions, params)),
                    SimplexCrossover{3}(1.05),
                    SimplexCrossover{2}(1.1),
                    SimulatedBinaryCrossover(0.05, 16.0),
                    SimulatedBinaryCrossover(0.05, 3.0),
                    SimulatedBinaryCrossover(0.1, 5.0),
                    SimulatedBinaryCrossover(0.2, 16.0),
                    UnimodalNormalDistributionCrossover{2}(chain(BlackBoxOptim.UNDX_DefaultOptions, params)),
                    UnimodalNormalDistributionCrossover{3}(chain(BlackBoxOptim.UNDX_DefaultOptions, params)),
                    ParentCentricCrossover{2}(chain(BlackBoxOptim.PCX_DefaultOptions, params)),
                    ParentCentricCrossover{3}(chain(BlackBoxOptim.PCX_DefaultOptions, params)),
]

"""
Switch from one activated var to another one that is connected to it
by the set overlap.
"""
struct VarSwitchMutation <: BlackBoxOptim.MutationOperator
    varXvar_probs::Vector{Tuple{Vector{Int}, Weights{Float64}}} # var -> var switching probabilities
end

function VarSwitchMutation(problem::MultiobjCoverProblem;
                           score_shape::Real = 0.1, vXv_cutoff::Real = -5.0)
    score_shape > 0 || throw(ArgumentError("score_shape must be positive"))
    varXvar_probs = Vector{Tuple{Vector{Int}, Weights{Float64}}}()
    for j in axes(problem.varXvar_scores, 2)
        vixs = Vector{Int}()
        w = Vector{Float64}()
        for i in axes(problem.varXvar_scores, 1)
            i == j && continue
            @inbounds vXv = problem.varXvar_scores[i, j]
            @assert vXv <= 0
            if vXv < vXv_cutoff
                push!(vixs, i)
                push!(w, 1 - exp(vXv * score_shape))
            end
        end
        push!(varXvar_probs, (vixs, weights(w)))
    end
    return VarSwitchMutation(varXvar_probs)
end

function BlackBoxOptim.apply!(m::VarSwitchMutation, v::AbstractVector{Float64}, index::Int)
    nnz = length(v) - sum(iszero, v)
    (nnz == 0) && return v
    # select a random vertex ix with non-zero weight to switch
    nzix = sample(1:nnz)
    ix = findfirst(!iszero, v)
    for _ in 1:(nzix-1)
        ix = findnext(!iszero, v, ix)
    end
    # select the neighbouring vertex to switch to
    jixs, jws = m.varXvar_probs[ix]
    isempty(jixs) && return v # nowhere to switch to
    jx = sample(jixs, jws)
    # do the switch
    v[jx] = max(v[jx], v[ix])
    v[ix] = 0.0
    return v
end

generate_modifier(problem::MultiobjCoverProblemBBOWrapper, params) =
    FixedGeneticOperatorsMixture(GeneticOperator[
                                 VarSwitchMutation(problem.orig),
                                 MutationClock(PolynomialMutation(search_space(problem), 30.0), 1/numdims(problem)),
                                 MutationClock(UniformMutation(search_space(problem)), 1/numdims(problem))],
                                 [0.5, 0.3, 0.2])

"""
The result of [`optimize`](@ref).
Contains the solutions on the Pareto front: weights of the annotation terms and corresponding cover scores.
"""
struct MultiobjCoverProblemResult
    var2set::Vector{Int}        # indices of sets in the original SetMosaic
    varscores::Vector{Float64}  # scores of the sets

    varweights::Matrix{Float64} # weights of the vars for each solution on Pareto front
    scores::Vector{DetailedScore} # detailed scores for each solution
    folded_scores::Vector{FoldedScore}
    agg_scores::Vector{Float64}
    best_ix::Int                # index of the best solution (agg_score minimum)
end

function MultiobjCoverProblemResult(problem::MultiobjCoverProblem,
                                    optres::BlackBoxOptim.OptimizationResults)
    # initialize solutions
    wmtx = similar(best_candidate(optres), nvars(problem),
                   length(pareto_frontier(optres)))
    orig_folded_scores = [archived_fitness(sol).orig for sol in pareto_frontier(optres)]
    front_perm = sortperm(orig_folded_scores) # sort lexicographically
    @inbounds for i in axes(wmtx, 2)
        wmtx[:, i] = BlackBoxOptim.params(pareto_frontier(optres)[front_perm[i]])
    end
    # remove small non-zero probabilities due to optimization method errors
    minw = problem.params.min_weight
    maxw = 1.0 - problem.params.min_weight
    @inbounds for i in eachindex(wmtx)
        if wmtx[i] <= minw
            wmtx[i] = 0.0
        elseif wmtx[i] >= maxw
            wmtx[i] = 1.0
        end
    end
    # calculate detailed scores for the corrected weights, then fold and aggregate
    scores = dropdims(mapslices(w -> score(w, problem), wmtx; dims=1), dims=1)
    folded_scores = fitness_scheme(optres).aggregator.score_folding.(scores)
    agg_scores = aggscore.(scores, Ref(problem.params))

    return MultiobjCoverProblemResult(problem.var2set, problem.var_scores,
                                      wmtx, scores, folded_scores, agg_scores,
                                      findmin(agg_scores)[2])
end

nvars(res::MultiobjCoverProblemResult) = length(res.var2set)
nsolutions(res::MultiobjCoverProblemResult) = size(res.varweights, 2)

Base.isempty(res::MultiobjCoverProblemResult) = isempty(res.varweights)

varweights(res::MultiobjCoverProblemResult, sol_ix::Integer) =
    view(res.varweights, :, sol_ix)

best_index(res::MultiobjCoverProblemResult, ::Nothing=nothing) = res.best_ix

function best_index(res::MultiobjCoverProblemResult, params::CoverParams)
    # FIXME replace with findmin(aggscore, res) when Julia 1.0 support is dropped
    best_ix = 1
    min_score = aggscore(res.scores[best_ix], params)
    @inbounds for i in 2:nsolutions(res)
        sol_score = aggscore(res.scores[i], params)
        if sol_score < min_score
            best_ix = i
            min_score = sol_score
        end
    end
    return best_ix
end

best_varweights(result::MultiobjCoverProblemResult) =
    @inbounds(varweights(result, best_index(result)))

best_aggscore(result::MultiobjCoverProblemResult) =
    @inbounds(result.agg_scores[best_index(result)])

function set2var(result::MultiobjCoverProblemResult, setix::Int)
    varix = searchsortedlast(result.var2set, setix)
    if varix > 0 && result.var2set[varix] == setix
        return varix
    else
        return 0
    end
end

"""
    optimize(problem::MultiobjCoverProblem,
             [opt_params::MultiobjOptimizerParams]) -> MultiobjCoverProblemResult

Optimize [`MultiobjCoverProblem`](@ref) and return the result.
Uses Borf multi-objective optimization method from *BlackBoxOptim.jl* package.
"""
function optimize(problem::MultiobjCoverProblem,
                  opt_params::MultiobjOptimizerParams = MultiobjOptimizerParams())
    if nvars(problem) == 0
        return MultiobjCoverProblemResult(
            Vector{Int}(), Vector{Float64}(), Matrix{Float64}(undef, 0, 0),
            Vector{DetailedScore}(), Vector{FoldedScore}(), Vector{Float64}(), 0)
    end

    score_folding = MultiobjProblemSoftFold2d(problem.params, opt_params.fold_ratio_threshold)
    bbowrapper = MultiobjCoverProblemBBOWrapper(problem, score_folding, opt_params.weight_digits)
    popmatrix = rand_individuals(search_space(bbowrapper), opt_params.pop_size, method=:latin_hypercube)
    # two extreme solutions and one totally neutral
    size(popmatrix, 2) > 0 && (popmatrix[:, 1] .= 0.0)
    size(popmatrix, 2) > 1 && (popmatrix[:, 2] .= 1.0)
    size(popmatrix, 2) > 2 && (popmatrix[:, 3] .= 0.5)
    size(popmatrix, 2) > 3 && (popmatrix[:, 4] .= sortperm(problem.var_scores, rev=true)./length(problem.var_scores))
    # solutions that select topN most significant sets
    score_qtls = quantile(problem.var_scores, range(0.0, stop=1.0, length=10))
    for i in 1:10
        (size(popmatrix, 2) < i + 4) && break
        popmatrix[:, i+4] .= ifelse.(problem.var_scores .<= score_qtls[i], 1.0, 0.0)
    end
    N = fieldcount(fitness_type(fitness_scheme_type(score_folding))) # FIXME numobjectives(bbowrapper) gives wrong result, looks like Julia problem
    population = FitPopulation(popmatrix, nafitness(IndexedTupleFitness{N,Float64}), ntransient=1)

    go_params = genop_params(opt_params)
    bboptimizer = BlackBoxOptim.BorgMOEA(bbowrapper, population,
            generate_recombinators(bbowrapper, go_params),
            generate_modifier(bbowrapper, go_params),
            RandomBound(search_space(bbowrapper)),
            borg_params(opt_params, problem))

    bboctrl = BlackBoxOptim.OptController(bboptimizer, bbowrapper,
                 bbo_ctrl_params(opt_params))
    bbores = bboptimize(bboctrl)
    BlackBoxOptim.isinterrupted(bbores) && throw(InterruptException())
    return MultiobjCoverProblemResult(problem, bbores)
end


text:
__precompile__()

module ReadXML
using LightXML

struct Fixed
    mkbound::Any
    mk::Any
    beg::Any
    count::Any
    property::Any
    bool::Any
end

struct Moving
    mkbound::Any
    mk::Any
    beg::Any
    count::Any
    property::Any
    refmotion::Any
    bool::Any
end

struct Floating
    mkbound::Any
    mk::Any
    beg::Any
    count::Any
    property::Any
    bool::Any
end

struct Fluid
    mkfluid::Any
    mk::Any
    beg::Any
    count::Any
    bool::Any
end

const BodyType =
    Dict("fixed" => Fixed, "moving" => Moving, "floating" => Floating, "fluid" => Fluid)

## Accesses the relevant part of the XML file array depending on BodyType ie. "k"
# is the array which "RelevantXmlSnip" produced. "BodySelect" is then used to
# grab the right structure.
function Access(k::Array{XMLElement,1}, Body::String)
    n = length(k)
    j = Array{BodyType[Body],1}(undef, n)
    for i = 1:n
        j[i] = BodySelect(k[i], Body)
    end
    return j
end

##To ease conversion... works for now.
function p(m)
    parse(Int64, m)
end

## Depending on choice of body seperate values are stored
function BodySelect(k, Body)
    mkbound = attribute(k, "mkbound")
    mkfluid = attribute(k, "mkfluid")
    mk = attribute(k, "mk")
    beg = attribute(k, "begin")
    count = attribute(k, "count")
    property = attribute(k, "property")
    refmotion = attribute(k, "refmotion")
    if Body == "fixed"
        Fixed(p(mkbound), p(mk), p(beg), p(count), property, false)
    elseif Body == "moving"
        Moving(p(mkbound), p(mk), p(beg), p(count), property, p(refmotion), false)
    elseif Body == "floating"
        Floating(p(mkbound), p(mk), p(beg), p(count), property, false)
    elseif Body == "fluid"
        Fluid(p(mkfluid), p(mk), p(beg), p(count), true)
    end
end

## This function returns the relevant xml snippet for each BodyType ie. "fluid",
## "floating" etc. Also returns the parsed xml file, so that it can be closed later.
function RelevantXmlSnip(
    filename::String,
    Bodies::Array{String,1} = ["fixed", "moving", "floating", "fluid"],
)
    nBody = length(Bodies)
    j = Array{Array{XMLElement,1}}(undef, nBody)
    xdoc = parse_file(filename)
    xroot = root(xdoc)
    ces = xroot["execution"]
    e1 = ces[1]
    t = find_element(e1, "particles")

    for i = 1:nBody
        j[i] = get_elements_by_tagname(t, Bodies[i])
    end

    return j, xdoc
end

## Generates final array based on former functions
function MkArray(
    filename::String,
    Bodies::Array{String,1} = ["fixed", "moving", "floating", "fluid"],
)
    y, xdoc = RelevantXmlSnip(filename)

    j = Array{Any,1}(undef, length(Bodies))

    for i = 1:length(j)
        j[i] = Access(y[i], Bodies[i])
    end

    free(xdoc)

    return j
end
end


text:
using SeisIO, SeisIO.RandSeis

# equivalence:
# FDSN wildcard         RegEx
#     ?           =>    .
#     *           =>    .*

function id_to_regex(cid::Array{UInt8,1})
  replace!(cid, 0x2e=>0x5c, 0x3f=>0x2e)
  L = length(cid)
  i = L
  while i > 0
    # replace: '*' => '.*'
    if cid[i] == 0x2a
      splice!(cid, i:i-1, 0x2e)
      i -= 1
    # now replace: '\' => '\.'
    elseif cid[i] == 0x5c
      if i == L
        push!(cid, 0x2e)
      else
        splice!(cid, i:i, [0x5c, 0x2e])
      end
    end
    i -= 1
  end
  return Regex(String(cid))
end

id_to_regex(id::AbstractString) = id_to_regex(copy(codeunits(id)))

function netsta_to_regex(id::AbstractString)
  cid = copy(codeunits(id))
  i = 0
  j = 0
  L = length(cid)
  while i < L
    i += 1
    if cid[i] == 0x2e
      j += 1
    end
    if j == 2
      deleteat!(cid, i:L)
      break
    end
  end
  return id_to_regex(cid)
end

function id_match(id::AbstractString, S::GphysData)
  j = findid(id, S.id)
  j > 0 && return [j]

  idr = id_to_regex(id)
  chans = Int64[]
  for (j, cid) in enumerate(S.id)
    if occursin(idr, cid)
      push!(chans, j)
    end
  end
  return chans
end


text:
using DataFrames
using Gadfly
using Distributions

function g(E)
  return (1./(sqrt(2*pi)*1.6e-9^3))*exp(-(E^2)/2.)
end

function s(g)
  prefac = (1./(sqrt(2*pi)*1.6e-9^3))
  return sqrt(2.*log(prefac/g))
end

Erand = DataFrame()

Erand[:evals] = collect(-5:0.001:5)
Erand[:gvals] = Array(Float64,size(Erand[:evals],1))
Erand[:svals] = Array(Float64,size(Erand[:evals],1))
Erand[:rands] = rand(Normal(0,1), size(Erand[:evals],1))
Erand[:ehist] = Array(Float64,size(Erand[:evals],1))
Erand[:label] = "E_rand"

prefac = (1./(sqrt(2*pi)*1.6e-9^3))

for i = 1:size(Erand[:evals],1)
  Erand[:gvals][i] = g(Erand[:evals][i])
  Erand[:svals][i] = prefac*exp(-(Erand[:rands][i]^2)/2)
  Erand[:ehist][i] = Erand[:evals]
end

p1 = plot(Erand,x="evals",y="gvals",Geom.line,Coord.Cartesian(xmin=-5.,xmax=5.))
p2 = plot(Erand,x="svals",color="label",Geom.histogram)
p3 = plot(Erand,x="rands",y="svals",Geom.line)
push!(p1,layer(Erand,x="rands",y="svals",color="label",Geom.point))

println("Ending script")

max_stars_repo_name: cnshingledecker/kmc-spin

max_issues_repo_name: cnshingledecker/kmc-spin

max_forks_repo_name: cnshingledecker/kmc-spin


text:
using Distributions, Random, Plots, LaTeXStrings; pyplot()
Random.seed!(1)

a, c, d = 2, 1, 5
sig = 0.02
next(x,y) = [a*x*(1-x) - x*y, -c*y + d*x*y]
equibPoint = [(1+c)/d ,(d*(a-1)-a*(1+c))/d]

initX = [0.8,0.05]
tEnd,tEndStoch = 100, 10

traj = [[] for _ in 1:tEnd]
trajStoch = [[] for _ in 1:tEndStoch]
traj[1], trajStoch[1] = initX, initX

for t in 2:tEnd
    traj[t] =  next(traj[t-1]...)
end

for t in 2:tEndStoch
    trajStoch[t] = next(trajStoch[t-1]...) + [rand(Normal(0,sig)),0.0]
end

scatter([traj[1][1]], [traj[1][2]],
	c=:black, ms=10, label="Initial state")
plot!(first.(traj), last.(traj),
	c=:blue, ls=:dash, m=(:dot, 5, Plots.stroke(0)),
	label="Deterministic trajectory")
plot!(first.(trajStoch),last.(trajStoch),
	c=:green, ls=:dash, m=(:dot, 5, Plots.stroke(0)),
	label="Stochastic trajectory")
scatter!([equibPoint[1]], [equibPoint[2]],
	c=:red, shape=:cross, ms=10, label="Equlibrium point",
	xlabel=L"X_1", ylabel=L"X_2")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# test.jl
# (c) 2016 E. Adrián Garro S. Costa Rica Institute of Technology.

include("tree.jl")

tree = Tree{Int64, Int64}()

seq = [16, 17, 11, 9, 5, 10, 12, 3, 19, 14, 13, 1, 4, 8, 15, 18, 7, 0, 6, 2]
for i in seq
	insert!(tree, i, i)
end

delete!(tree, 16)
delete!(tree, 14)
delete!(tree, 2)
insert!(tree, 23, 23)

for value in @task inorderwalk(tree)
	println(value)
end

prettyprint(tree)

println(isbst(tree))

println(isbalanced(tree))

max_stars_repo_name: adriangarro/llrb

max_issues_repo_name: adriangarro/llrb

max_forks_repo_name: adriangarro/llrb


text:
time_modeling(model::Model, srcGeometry, srcData, recGeometry, recData, perturbation, srcnum::UnitRange{Int64}, op::Char, mode::Int64) =
    time_modeling(model, srcGeometry, srcData, recGeometry, recData, perturbation, srcnum, op, mode, Options())

function time_modeling(model::Model, srcGeometry, srcData, recGeometry, recData, dm, srcnum::UnitRange{Int64}, op::Char, mode::Int64, options)
# time_modeling function for multiple sources. Depending on the operator and mode, this function distributes the sources
# and if applicable the input data amongst the available workers.

    results = judipmap(j -> time_modeling(model, subsample(srcGeometry,j), subsample(srcData, j),
                                          subsample(recGeometry,j), subsample(recData, j), dm,
                                          op, mode, subsample(options, j)), srcnum)

    if op=='F' || (op=='J' && mode==1)
        argout1 = vcat(results...)
    elseif op=='J' && mode==-1
        argout1 = sum(results)
    else
        error("operation no defined")
    end
    return argout1
end


text:
using Test
using SafeTestsets

@testset verbose=true "TensorNetworkCodes.jl" begin
    @safetestset "types.jl" begin include("types.jl") end
    @safetestset "pauli_functions.jl" begin include("pauli_functions.jl") end
    @safetestset "code_functions.jl" begin include("code_functions.jl") end
    @safetestset "code_graph_functions.jl" begin include("code_graph_functions.jl") end
    @safetestset "contraction_functions.jl" begin include("contraction_functions.jl") end
    @safetestset "plotting_functions.jl" begin include("plotting_functions.jl") end
    @safetestset "examples/simple.jl" begin include("examples/simple.jl") end
    @safetestset "examples/surface.jl" begin include("examples/surface.jl") end
    @safetestset "itensors_functions.jl" begin include("itensors_functions.jl") end
    @safetestset "TNDecode.jl" begin include("TNDecode.jl") end
    @safetestset "TNDistance.jl" begin include("TNDistance.jl") end
    @safetestset "QecsimAdaptors.jl" begin include("QecsimAdaptors.jl") end
end


text:
using BenchmarkTools, LinearAlgebra
import Zygote
import Zygote: @grad
using Test
import AutoGrad, YAAD

Zygote.@grad LinearAlgebra.tr(x) = LinearAlgebra.tr(x), Δ-> (Δ * Matrix(I, size(x)), )
Zygote.@grad Base.:(*)(lhs::Matrix, rhs::Matrix) = gemm(lhs, rhs), grad -> (grad * transpose(rhs), transpose(lhs) * grad)

function bench_tr_mul_yaad(x1, x2)
    z = tr(x1 * x2)
    YAAD.backward(z)
    x1.grad, x2.grad
end

function bench_tr_mul_autograd(x1, x2)
    z = AutoGrad.@diff tr(x1 * x2)
    AutoGrad.grad(z, x1), AutoGrad.grad(z, x2)
end

function bench_tr_mul_zygote(x1, x2)
    Zygote.gradient((x1, x2)->tr(x1 * x2), x1, x2)
end

xv, yv = rand(30, 30), rand(30, 30)

yaad_x = YAAD.Variable(xv)
yaad_y = YAAD.Variable(yv)

autograd_x = AutoGrad.Param(xv)
autograd_y = AutoGrad.Param(yv)

println("Check gradient:")

@test bench_tr_mul_yaad(yaad_x, yaad_y) == bench_tr_mul_autograd(autograd_x, autograd_y)

#  == bench_tr_mul_zygote(xv, yv)
println("----------------------------------------------------------------------")
println("YAAD:")
display(@benchmark bench_tr_mul_yaad(yaad_x, yaad_y))

println()
println("----------------------------------------------------------------------")
println("AutoGrad:")
display(@benchmark bench_tr_mul_autograd(autograd_x, autograd_y))

println()
println("----------------------------------------------------------------------")
println("Zygote:")
display(@benchmark bench_tr_mul_zygote(xv, yv))


text:
"""
    StraightLineHomotopy(start::AbstractPolySystem, target::AbstractPolySystem)

Constructs the homotopy `t * start + (1-t) * target`.
"""
struct StraightLineHomotopy{T<:Number, Start<:AbstractPolySystem{T}, Target<:AbstractPolySystem{T}} <: AbstractHomotopy{T}
    start::Start
    target::Target

    function StraightLineHomotopy{T, Start, Target}(start::Start, target::Target) where {T<:Number, Start<:AbstractPolySystem{T}, Target<:AbstractPolySystem{T}}
        if (homogenized(start) != homogenized(target))
            return error("start and target have to be both either homogenized or not")
        end
        if (nvariables(start) != nvariables(target))
            return error("The number of unknowns of start and target system doesn't match." *
                "The start system of the homotopy has $(nvariables(start)) unknowns " *
                "and the target system has $(nvariables(target)) unknowns.")
        end

        if (length(target) != length(start))
            return error("The number of equations of start and target system doesn't match." *
                "The start system of the homotopy has $(length(start)) equations " *
                "and the target system has $(length(target)) equations.")
        end

        new(start, target)
    end
end

function StraightLineHomotopy(start::Start,target::Target) where {T<:Number, Start<:AbstractPolySystem{T}, Target<:AbstractPolySystem{T}}
    StraightLineHomotopy{T, Start, Target}(start,target)
end
# promote if they don't have the same coefficients
function StraightLineHomotopy(start::AbstractPolySystem,target::AbstractPolySystem)
    T = promote_type(coefftype(start), coefftype(target))
    Start = promote_type(T, typeof(start))
    Target = promote_type(T, typeof(target))
    StraightLineHomotopy{T, Start, Target}(convert(Start, start), convert(Target, target))
end
StraightLineHomotopy(s, t) = StraightLineHomotopy(PolySystem(s), PolySystem(t))

function evaluate(H::StraightLineHomotopy, x, t)
    (1-t) .* evaluate(H.target, x) .+ t .* evaluate(H.start, x)
end
(H::StraightLineHomotopy)(x,t) = evaluate(H,x,t)

function substitute(H::StraightLineHomotopy, pair::Pair{Symbol,<:Number})
    StraightLineHomotopy(substitute(H.start, pair), substitute(H.target, pair))
end

function removepoly(H::StraightLineHomotopy, i::Int)
    StraightLineHomotopy(removepoly(H.start,i), removepoly(H.target,i))
end

function Base.show(io::IO, H::StraightLineHomotopy)
    print(io, typeof(H), ":\n")
    println(io, "* start:")
    println(io, H.start)
    println(io, "* target:")
    print(io, H.target)
end

function Base.promote_type(::Type{StraightLineHomotopy{T, Start, Target}}, ::Type{S}) where {S,T, Start, Target}
    StraightLineHomotopy{promote_type(T,S), promote_type(Start, S), promote_type(Target, S)}
end
function Base.promote_type(::Type{StraightLineHomotopy{T, Start1, Target1}}, ::StraightLineHomotopy{S, Start2, Target2}) where {S,T, Start1, Start2, Target1, Target2}
    StraightLineHomotopy{promote_type(T,S), promote_type(Start1, Start2), promote_type(Target1, Target2)}
end

function Base.convert(::Type{StraightLineHomotopy{T, Start, Target}}, H::StraightLineHomotopy) where {T, Start, Target}
    StraightLineHomotopy{T, Start, Target}(convert(Start, H.start), convert(Target, H.target))
end

"""
    differentiate(H::StraightLineHomotopy)

Differentiates `H` w.r.t. the place and returns an evaluation function
`(x,t) -> t * J_s(x) + (1 - t) * J_t(x)`
where `J_s` is the jacobian of the start system and
`J_t` is the jacobian of the target system.
"""
function differentiate(H::StraightLineHomotopy)
    J_s = differentiate(H.start)
    J_t = differentiate(H.target)

    (x, t) -> t .* J_s(x) .+ (1 - t) .* J_t(x)
end
"""
    dt(H::StraightLineHomotopy)

Differentiates `H` w.r.t. the the time and returns an evaluation function
`(x,t) -> T`
"""
dt(H::StraightLineHomotopy) = (x, ::Number) -> evaluate(H.start, x) .- evaluate(H.target, x)

"""
    homogenize(H::StraightLineHomotopy)

Make the start and target system of `H` homogenous.
"""
function homogenize(H::StraightLineHomotopy)
    StraightLineHomotopy(homogenize(H.start), homogenize(H.target))
end

"""
    homogenized(H::StraightLineHomotopy)

Checks whether `H` was homogenized.
"""
homogenized(H::StraightLineHomotopy) = homogenized(H.start)

ishomogenous(H::StraightLineHomotopy) = ishomogenous(H.start) && ishomogenous(H.target)

"""
    nvariables(H::StraightLineHomotopy)

The number variables of `H`.
"""
function nvariables(H::StraightLineHomotopy)
    # Constructor guarantees that start and target have the same number of variables
    nvariables(H.start)
end
function length(H::StraightLineHomotopy)
    # Constructor guarantees that start and target have the same length
    length(H.start)
end
"""
    degrees(H::StraightLineHomotopy)

The (total) degrees of the polynomials of the system `P`.
"""
degrees(H::StraightLineHomotopy) = max.(degrees(H.start), degrees(H.target))
"""
    startsystem(H::StraightLineHomotopy)

The startsystem of `H`.
"""
startsystem(H::StraightLineHomotopy) = H.start
"""
    targetsystem(H::StraightLineHomotopy)

The targetsystem of `H`.
"""
targetsystem(H::StraightLineHomotopy) = H.target

"""
    weylnorm(H, t)

Computes the weyl_norm of the homotopy to the given time `t`.

## Explanation
For ``H = (1-t)F+tG`` we have
```math
\begin{align*}
<H,H> &= <(1-t)F+tG,(1-t)F+tG> \\
      &= <(1-t)F,(1-t)F+tG> + <tG,(1-t)F+tG> \\
      &= <(1-t)F,(1-t)F> + <(1-t)F,tG> + <tG,(1-t)F> + <tG,tG> \\
      &= <(1-t)F,(1-t)F> + 2*Real(<(1-t)F,tG>) + <tG,tG> \\
      &= (1-t)^2<F,F> + (1-t)t*2*Real(<F,G>) + t^2<G,G> \\
\end{align*}
```
"""
function weylnorm(H::StraightLineHomotopy{T}, t::Float64) where {T<:Complex}
    F = H.target
    G = H.start

    sqrt((1-t)^2 * weyldot(F,F) + (1-t)*t*2*real(weyldot(F,G)) + t^2 * weyldot(G,G))
end


text:
### Algoritms and Recursion #####

# Guess-and-check algoritm

function guessandckeck(x)
    #initialize variable
    counter = 1
    # algoritm
    for guess in 0:abs(x)
        if guess^3 > abs(x)
            println("\n$x is not a perfect cube!")
            break
        end
        if guess^3 != abs(x)
            println(counter, "\tGuess = $guess\tGuess Cubed = $(guess^3)")
            counter += 1
        else
            if x < 0
                guess = -guess
            end
            println(counter, "\tGuess = $guess\tGuess cubed = $(guess^3)")
            println("\nThe cube root of $x is $guess.")
            break
        end
    end
end

guessandckeck(7)
guessandckeck(8)
guessandckeck(-8)

# guess and check only works for cubed value and not perfect for negative value. So lets try another algoritm

# Approximate solution

function approximation(x)

    # initialize
    guess = 0.0
    counter = 1
    increment = 0.01
    sensetivity = 0.1

    # algoritm
    while abs(guess^3 - x) >= sensetivity && abs(guess^3) <= abs(x)
        println(counter, "\tGuess = $guess\tGuess cubed = $(guess^3)")
        guess += increment
        counter += 1
    end

    println(counter, "\tGuess = $guess\tGuess cubed = $(guess^3)")
    guess = round(guess, digits = 2)
    x < 0 ? guess = -guess : guess = guess
    println("\nThe cube root of $x is approximately $guess.")
end

approximation(6)

# approximation algoritm works fine but not efficient for large numbers.  So,  Lets look at another algoritm

# Binary search algoritm

function binarysearch(x)

    # initialize

    low = 0
    high = x
    guess = (low + high) /2
    counter = 1
    sensetivity = 0.01

    #algoritm

    while abs(guess^3 -x) ≥ sensetivity
        println(counter, "\tGuess = $guess\tGuess cubed = $(guess^3)")

        if abs(guess^3) < abs(x)
            low = guess
        else
            high = guess
        end

        guess  = (low + high) /2
        counter += 1

    end
    println(counter, "\tGuess = $guess\tGuess cubed = $(guess^3)")
    guess = round(guess, digits = 3)
    println("\nThe cube root of $x is approximately $guess.")
end

binarysearch(6)

# performance comparation

using BenchmarkTools

@benchmark guessandckeck(1_000_000)
@benchmark approximation(1_000_000)
@benchmark binarysearch(1_000_000)

# Recursion

# recursion is a process that divides a problem into sub problem to solve larger one

# think on multipication, Z × 5  =  z + z + z + z + z

# multipication as iteravive solution

function multi_iterative(x,y)
    result = 0

    for i in 1:y
        result += x
        println("$i\t$x * $i\t= $result")
    end

    println("\n$x * $y = $result\n")
    return result
end

multi_iterative(5,10)

# multipication as recursion

function multi_recursive(x,y)
    if y == 1

        println("multi_recursive($x, $y)\t = $x\n")
        return x
    else

        println("multi_recursive($x, $y)\t = $x + multi_recursive($x, $(y - 1))")
        return x + multi_recursive(x, y - 1)
    end

end

multi_recursive(5,10)

# factorial as iterative method

function FactorialIterative(x)
    result = 1

    for i in 1:x
        result *= i

        println("$i factorial = $result")
    end

    println("\nThe Factorial of $x is $result\n")
    return result
end

FactorialIterative(5)

# factorial as recursion method

function FactorialRecursive(x)

    if x == 1

        println("FactorialRecursive($x) = $x\n")
        return 1
    else
        println("FactorialRecursive($x) = $x * FactorialRecursive($(x - 1))")
        return x * FactorialRecursive(x - 1)
    end
end


text:
using VLConstraintBasedModelGenerationUtilities

# setup path to sequence file -
path_to_gene_sequence_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/G-MZ373340.fasta"

# load -
result = build_gene_table(path_to_gene_sequence_file) |> check