text:
# Really dumb variable name replacement macro.  Doesn't know about scopes or
# the gory details of any `Expr`s
function replace_var!(ex, varname, replacement)
    if ex isa Expr
        map!(e->replace_var!(e, varname, replacement), ex.args, ex.args)
        return ex
    elseif ex === varname
        return replacement
    else
        return ex
    end
end

function expand_underscores_!(ex)
    if ex === :_
        return true
    end
    if ex isa Expr
        if (ex.head == :call && ex.args[1] == :|>)
            map!(expand_underscores!, ex.args, ex.args)
            return false
        else
            return any(expand_underscores_!, ex.args)
        end
    else
        return false
    end
end

function expand_underscores!(ex)
    if expand_underscores_!(ex)
        return replace_var!(Expr(:->, :_, ex), :_, gensym("_"))
    else
        ex
    end
end

"""
    @_ ex

Expand an expression containing underscore placeholders into a lambda.  The
lambda head is placed at the outermost scope, except where interrupted by `|>`
pipe operators.

# Examples

`data |> reduce(+,_)`  expands to  `data |> x->reduce(+,x)`

`data |> foo(bar(_))`  expands to  `data |> x->foo(bar(x))`

"""
macro _(ex)
    esc(expand_underscores!(ex))
end


text:
##############################################################################################
module DeSolvers
##############################################################################################
# imports and exports
using Random, Distributions
export rk4_step!, tay2_step!, em_step!

##############################################################################################
##############################################################################################
# Type union declarations for multiple dispatch and type aliases

# vectors and ensemble members of sample
VecA = Union{Vector{Float64}, SubArray{Float64, 1}}

# dictionary for model parameters
ParamDict = Union{Dict{String, Array{Float64}}, Dict{String, Vector{Float64}}}

# dictionary containing key and index pairs to subset the state vector and
# merge with dx_params
ParamSample = Dict{String, Vector{UnitRange{Int64}}}

##############################################################################################
##############################################################################################
# four-stage Runge-Kutta scheme

function rk4_step!(x::VecA, t::Float64, kwargs::Dict{String,Any})
    """Step of integration rule for 4 stage Runge-Kutta as discussed in Grudzien et al. 2020

    The rule has strong convergence order 1.0 for generic SDEs and order 4.0 for ODEs
    Arguments are given as
    x            -- array or sub-array of a single state possibly including parameter values
    t            -- time point
    kwargs       -- should include dx_dt, the paramters for the dx_dt and optional arguments
    dx_dt        -- time derivative function with arguments x and dx_params
    dx_params    -- parameters necessary to resolve dx_dt, not including
                    parameters to be estimated in the extended state vector
    h            -- numerical discretization step size
    diffusion    -- tunes the standard deviation of the Wiener process,
                    equal to sqrt(h) * diffusion
    diff_mat     -- structure matrix for the diffusion coefficients,
                    replaces the default uniform scaling
    state_dim    -- keyword for parameter estimation, dimension of the
                    dynamic state < dimension of full extended state
    param_sample -- ParamSample dictionary for merging extended state with dx_params
    ξ            -- random array size state_dim, can be defined in kwargs
                    to provide a particular realization for method validation
    """

    # unpack the integration scheme arguments and the parameters of the derivative
    h = kwargs["h"]::Float64
    diffusion = kwargs["diffusion"]::Float64
    dx_dt = kwargs["dx_dt"]

    if haskey(kwargs, "dx_params")
        # get parameters for resolving dx_dt
        dx_params = kwargs["dx_params"]::ParamDict
    end

    # infer the (possibly) extended state dimension
    sys_dim = length(x)

    # check if extended state vector
    if haskey(kwargs, "param_sample")
        param_sample = kwargs["param_sample"]::ParamSample
        state_dim = kwargs["state_dim"]::Int64
        v = @view x[begin: state_dim]
        param_est = true
    else
        # the state dim equals the system dim
        state_dim = sys_dim
        v = @view x[begin: state_dim]
        param_est = false
    end

    # check if SDE formulation
    if diffusion != 0.0
	    if haskey(kwargs, "ξ")
            # pre-computed perturbation is provided for controlled run
	        ξ = kwargs["ξ"]::Array{Float64,2}
	    else
            # generate perturbation for brownian motion if not neccesary to reproduce
            ξ = rand(Normal(), state_dim)
        end
        if haskey(kwargs, "diff_mat")
            # diffusion is a scalar intensity which is applied to the
            # structure matrix for the diffusion coefficients
            diff_mat = kwargs["diff_mat"]::Array{Float64}
            diffusion = diffusion * diff_mat
        end
        # rescale the standard normal to variance h for Wiener process
        W = ξ * sqrt(h)
    end

    # load parameter values from the extended state into the derivative
    if param_est
        if haskey(kwargs, "dx_params")
            # extract the parameter sample and append to other derivative parameters
            for key in keys(param_sample)
                dx_params = merge(dx_params, Dict(key => x[param_sample[key][1]]))
            end
        else
            # set the parameter sample as the only derivative parameters
            dx_params = Dict{String, Array{Float64}}
            for key in keys(param_sample)
                dx_params = merge(dx_params, Dict(key => x[param_sample[key][1]]))
            end
        end
    end

    # pre-allocate storage for the Runge-Kutta scheme
    κ = Array{Float64}(undef, state_dim, 4)

    # terms of the RK scheme recursively evolve the dynamic state components alone
    if diffusion != 0.0
        # SDE formulation
        κ[:, 1] = dx_dt(v, t, dx_params) * h + diffusion * W
        κ[:, 2] = dx_dt(v + 0.5 * κ[:, 1], t + 0.5 * h, dx_params) * h + diffusion * W
        κ[:, 3] = dx_dt(v + 0.5 * κ[:, 2], t + 0.5 * h, dx_params) * h + diffusion * W
        κ[:, 4] = dx_dt(v + κ[:, 3], t + h, dx_params) * h + diffusion * W
    else
        # deterministic formulation
        κ[:, 1] = dx_dt(v, t, dx_params) * h
        κ[:, 2] = dx_dt(v + 0.5 * κ[:, 1], t + 0.5 * h, dx_params) * h
        κ[:, 3] = dx_dt(v + 0.5 * κ[:, 2], t + 0.5 * h, dx_params) * h
        κ[:, 4] = dx_dt(v + κ[:, 3], t + h, dx_params) * h
    end

    # compute the update to the dynamic variables
    x[begin: state_dim] = v + (1.0 / 6.0) * (κ[:, 1] + 2.0*κ[:, 2] + 2.0*κ[:, 3] + κ[:, 4])
    x
end

##############################################################################################
# deterministic 2nd order Taylor Method

function tay2_step!(x::Vector{Float64}, t::Float64, kwargs::Dict{String,Any})
    """Second order Taylor method for step size h and state vector x.

    Arguments are given as
    x          -- array of a single state possibly including parameter values
    kwargs     -- should include dx_dt, the paramters for the dx_dt and optional arguments
    dx_dt      -- time derivative function with arguments x and dx_params
    dx_params  -- parameters necessary to resolve dx_dt, not including
                  parameters to be estimated in the extended state vector
    h          -- numerical discretization step size
    """

    # unpack dx_params
    h = kwargs["h"]::Float64
    dx_params = kwargs["dx_params"]::ParamDict
    dx_dt = kwargs["dx_dt"]
    jacobian = kwargs["jacobian"]

    # calculate the evolution of x one step forward via the second order Taylor expansion

    # first derivative
    dx = dx_dt(x, t, dx_params)

    # second order taylor expansion
    x .= x + dx * h + 0.5 * jacobian(x, t, dx_params) * dx * h^2.0
end

##############################################################################################
# Euler-Murayama step

function em_step!(x::Vector{Float64}, t::Float64, kwargs::Dict{String,Any})
    """This will propagate the state x one step forward by Euler-Murayama

    Step size is h, the Wiener process is assumed to have a scalar diffusion coefficient"""

    # unpack the arguments for the integration step
    h = kwargs["h"]::Float64
    dx_params = kwargs["dx_params"]::ParamDict
    diffusion = kwargs["diffusion"]::Float64
    dx_dt = kwargs["dx_dt"]
    state_dim = length(x)

    # check if SDE or deterministic formulation
    if diffusion != 0.0
	    if haskey(kwargs, "ξ")
            # pre-computed perturbation is provided for controlled run
	        ξ = kwargs["ξ"]::Array{Float64,2}
	    else
            # generate perturbation for brownian motion if not neccesary to reproduce
            ξ = rand(Normal(), state_dim)
        end

    else
        # if deterministic Euler, load dummy ξ of zeros
        ξ = zeros(state_dim)
    end

    # rescale the standard normal to variance h for Wiener process
    W = ξ * sqrt(h)

    # step forward by interval h
    x .= x +  h * dx_dt(x, t, dx_params) + diffusion * W
end

##############################################################################################
# end module

end

ext: jl

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
#==
All differentials need to define + and *.
That happens here.

We just use @eval to define all the combinations for AbstractDifferential
subtypes, as we know the full set that might be encountered.
Thus we can avoid any ambiguities.

Notice:
    The precedence goes:
    `Wirtinger, Zero, DoesNotExist, One, AbstractThunk, Composite, Any`
    Thus each of the @eval loops creating definitions of + and *
    defines the combination this type with all types of  lower precidence.
    This means each eval loops is 1 item smaller than the previous.
==#

function Base.:*(a::Wirtinger, b::Wirtinger)
    error("""
          Cannot multiply two Wirtinger objects; this error likely means a
          `WirtingerRule` was inappropriately defined somewhere. Multiplication
          of two Wirtinger objects is not defined because chain rule application
          often expands into a non-commutative operation in the Wirtinger
          calculus. To put it another way: simply given two Wirtinger objects
          and no other information, we can't know "locally" which components to
          conjugate in order to implement the chain rule. We could pick a
          convention; for example, we could define `a::Wirtinger * b::Wirtinger`
          such that we assume the chain rule application is of the form `f_a ∘ f_b`
          instead of `f_b ∘ f_a`. However, picking such a convention is likely to
          lead to silently incorrect derivatives due to commutativity assumptions
          in downstream generic code that deals with the reals. Thus, ChainRulesCore
          makes this operation an error instead.
          """)
end

function Base.:+(a::Wirtinger, b::Wirtinger)
    return Wirtinger(+(a.primal, b.primal), a.conjugate + b.conjugate)
end

for T in (:Zero, :DoesNotExist, :One, :AbstractThunk, :Any)
    @eval Base.:+(a::Wirtinger, b::$T) = a + Wirtinger(b, Zero())
    @eval Base.:+(a::$T, b::Wirtinger) = Wirtinger(a, Zero()) + b

    @eval Base.:*(a::Wirtinger, b::$T) = Wirtinger(a.primal * b, a.conjugate * b)
    @eval Base.:*(a::$T, b::Wirtinger) = Wirtinger(a * b.primal, a * b.conjugate)
end

Base.:+(::Zero, b::Zero) = Zero()
Base.:*(::Zero, ::Zero) = Zero()
for T in (:DoesNotExist, :One, :AbstractThunk, :Any)
    @eval Base.:+(::Zero, b::$T) = b
    @eval Base.:+(a::$T, ::Zero) = a

    @eval Base.:*(::Zero, ::$T) = Zero()
    @eval Base.:*(::$T, ::Zero) = Zero()
end

Base.:+(::DoesNotExist, ::DoesNotExist) = DoesNotExist()
Base.:*(::DoesNotExist, ::DoesNotExist) = DoesNotExist()
for T in (:One, :AbstractThunk, :Any)
    @eval Base.:+(::DoesNotExist, b::$T) = b
    @eval Base.:+(a::$T, ::DoesNotExist) = a

    @eval Base.:*(::DoesNotExist, ::$T) = DoesNotExist()
    @eval Base.:*(::$T, ::DoesNotExist) = DoesNotExist()
end

Base.:+(a::One, b::One) = extern(a) + extern(b)
Base.:*(::One, ::One) = One()
for T in (:AbstractThunk, :Any)
    @eval Base.:+(a::One, b::$T) = extern(a) + b
    @eval Base.:+(a::$T, b::One) = a + extern(b)

    @eval Base.:*(::One, b::$T) = b
    @eval Base.:*(a::$T, ::One) = a
end

Base.:+(a::AbstractThunk, b::AbstractThunk) = unthunk(a) + unthunk(b)
Base.:*(a::AbstractThunk, b::AbstractThunk) = unthunk(a) * unthunk(b)
for T in (:Any,)
    @eval Base.:+(a::AbstractThunk, b::$T) = unthunk(a) + b
    @eval Base.:+(a::$T, b::AbstractThunk) = a + unthunk(b)

    @eval Base.:*(a::AbstractThunk, b::$T) = unthunk(a) * b
    @eval Base.:*(a::$T, b::AbstractThunk) = a * unthunk(b)
end

################## Composite ##############################################################

# We intentionally do not define, `Base.*(::Composite, ::Composite)` as that is not meaningful
# In general one doesn't have to represent multiplications of 2 differentials
# Only of a differential and a scaling factor (generally `Real`)
Base.:*(s::Any, comp::Composite) = map(x->s*x, comp)
Base.:*(comp::Composite, s::Any) = map(x->x*s, comp)

function Base.:+(a::Composite{P}, b::Composite{P}) where P
    data = elementwise_add(backing(a), backing(b))
    return Composite{P, typeof(data)}(data)
end
function Base.:+(a::P, d::Composite{P}) where P
    try
        return construct(P, elementwise_add(backing(a), backing(d)))
    catch err
        throw(PrimalAdditionFailedException(a, d, err))
    end
end
Base.:+(a::Composite{P}, b::P) where P = b + a


text:
"""

`bmm(A, B))` performs a batch matrix-matrix product of matrices stored in `A`
and `B`. `A` and `B` must be 3d and the last dimension represents the batch size.

If A is a (m,n,b) tensor, B is a (n,k,b) tensor, and the output is a (m,k,b) tensor.

"""
function bmm(A, B)
    m,n,bs = size(A)
    n,k,bs = size(B)
    return bmm!(A, B, similar(A, (m, k, bs)))
end

function bmm(A::KnetArray, B::KnetArray{T}) where {T}
    m,n,bs = size(A)
    n,k,bs = size(B)
    return bmm!('N','N',one(T),A,B,zero(T),similar(A, (m, k, bs)))
end

function bmm!(transA::AbstractChar, transB::AbstractChar, alpha::Number, A::KnetArray{T}, B::KnetArray{T}, beta::Number, C::KnetArray{T}) where {T}
    cublasop(c::Char)=(if c=='N'; 0; elseif c=='T'; 1; elseif c=='C'; 2; else error("Unknown cublas op $c"); end)
    if ndims(A) != 3 || ndims(B) != 3
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    ma,ka,bsa = size(A)
    kb,nb,bsb = size(B)
    if bsa != bsb
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    bs = bsa
    if transA == 'N'
        m=ma; k=ka;
    else
        m=ka; k=ma;
    end
    if transB == 'N'
        k == kb || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
        n=nb; k==kb;
    else
        k == nb || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
        n=kb; k==nb;
    end
    (m == size(C,1) && n == size(C,2) && bs == size(C,3)) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    lda,ldb,ldc=ma,ka,ma
    transa = cublasop(transA); transb = cublasop(transB)
    alpha = T[alpha]; beta = T[beta]
    strideA, strideB, strideC = m*k, k*n, m*n
    if T<:Float64
        @cublas(cublasDgemmStridedBatched, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Clonglong, Ptr{T}, Cint, Clonglong, Ptr{T}, Ptr{T}, Cint, Clonglong, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, strideA, B, ldb, strideB, beta, C, ldc, strideC, bs)
    elseif T<:Float32
        @cublas(cublasSgemmStridedBatched, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Clonglong, Ptr{T}, Cint, Clonglong, Ptr{T}, Ptr{T}, Cint, Clonglong, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, strideA, B, ldb, strideB, beta, C, ldc, strideC, bs)
    else
        error("CUBLAS does not support $T")
    end
    return C
end

function bmm!(A, B, C)
    if ndims(A) != 3 || ndims(B) != 3
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    ma,ka,bsa = size(A)
    kb,nb,bsb = size(B)

    (bsa == bsb && ka == kb && size(C,1) == ma && size(C,2) == nb && size(C,3) == bsa) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))

    for i=1:bsa
        C[:, :, i] = view(A, :, :, i) * view(B, :, :, i)
    end
    return C
end

@primitive bmm(x1,x2),dy,y bmm(dy,permutedims(x2, [2,1,3])) bmm(permutedims(x1,[2,1,3]), dy)
@zerograd bmm!(transA::AbstractChar, transB::AbstractChar, alpha::Number, A::KnetArray, B::KnetArray, beta::Number, C::KnetArray)
@zerograd bmm!(A, B, C)


text:
using Yao
using BenchmarkTools

using Yao.Intrinsics
using LuxurySparse
import Yao.Intrinsics: basis

include("swap.jl")

msk = bmask(2,5)
state = rand(ComplexF64, 1<<16, 1)

bench = BenchmarkGroup()
bg = bench["Basis"] = BenchmarkGroup()
bg["swap"] = @benchmarkable swapbits2(236, $msk)

bg = bench["SwapGate"] = BenchmarkGroup()
bg["matrix"] = @benchmarkable swapgate(ComplexF64, 16, 3, 7)
bg["apply-mat"] = @benchmarkable swapapply!($(state), 7, 3)
bg["apply-vec"] = @benchmarkable swapapply!($(vec(state)), 7, 3)

#bg = bench["SingleControlXYZ"] = BenchmarkGroup()
#@benchmark cyapply!(s, 7, 1, 3)
#@benchmark czapply!(s, 7, 1, 3)
#@benchmark cxapply!(s, 15, 13, 3)
#@benchmark cxapply!(s, (7, 2), (1, 0), 3)
#@benchmark cyapply!(s, (7, 2), (1, 0), 3)
#@benchmark czapply!(s, (7, 2), (1, 0), 3)

showall(run(bench, verbose=true))


text:
# 2x + 3y
@test f(2,1) == 7
@test f(2,3) == 13


text:
module DrivenJumpProblems

export DrivenJumpProblem

using DiffEqJump
using SciMLBase
using PathWeightSampling: Trajectory
using StaticArrays
using StochasticDiffEq

"""
    IdentityMap()

A simple helper type that always returns `i` when indexed as `IdentityMap()[i]`.
"""
struct IdentityMap end
@inline Base.getindex(x::IdentityMap, i::Integer) = i

mutable struct TrajectoryCallback{T,IndexMap}
    traj::T
    index::Int
    index_map::IndexMap
end

TrajectoryCallback(traj, index_map=IdentityMap()) = TrajectoryCallback(traj, 1, index_map)

function (tc::TrajectoryCallback)(integrator::DiffEqBase.DEIntegrator) # affect!
    traj = tc.traj
    tc.index = min(tc.index + 1, length(traj.t))
    cond_u = traj.u[tc.index]
    for i in eachindex(cond_u)
        if integrator.u isa StaticVector
            integrator.u = setindex(integrator.u, cond_u[i], tc.index_map[i])
        else
            integrator.u[tc.index_map[i]] = cond_u[i]
        end
    end
    # it is important to call this to properly update reaction rates
    DiffEqJump.reset_aggregated_jumps!(integrator; update_jump_params=false)
    nothing
end

function (tc::TrajectoryCallback)(u, t::Real, i::DiffEqBase.DEIntegrator)::Bool # condition
    tc_index = tc.index
    traj_t = tc.traj.t
    traj_len = length(traj_t)
    @inbounds tcb = traj_t[tc_index]
    while tc_index < traj_len && t > tcb
        tc_index += 1
        @inbounds tcb = traj_t[tc_index]
    end
    tc.index = tc_index
    t == tcb
end

"""
    DrivenJumpProblem(jump_problem, driving_trajectory; index_map = IdentityMap(), save_jumps=false)

Create a `DrivenJumpProblem` from a base `jump_problem` and a `driving_trajectory`.

The optional argument `index_map` specifies how the components from `driving_trajectory`
map onto the components of the `jump_problem`. The default `IdentityMap` maps the N components
of the `driving_trajectory` onto the first N components of the jump problem, leaving the
remaining components of the jump problem unaltered by the driving trajectory.
"""
struct DrivenJumpProblem{Prob,Cb,TG}
    prob::Prob
    callback::Cb
    make_trajectory::TG
end

function DrivenJumpProblem(jump_problem::JP, driving_trajectory::Trajectory; index_map=IdentityMap(), save_jumps=false) where {JP}
    tcb = TrajectoryCallback(driving_trajectory, index_map)
    callback = DiscreteCallback(tcb, tcb, save_positions=(false, save_jumps))
    function make_trajectory(prob::DrivenJumpProblem)
        driving_trajectory
    end
    DrivenJumpProblem(jump_problem, callback, make_trajectory)
end

function DrivenJumpProblem(jump_problem::JP, driving_problem::SDEProblem; index_map=IdentityMap(), save_jumps=false, sde_dt=0.01) where {JP}
    tcb = TrajectoryCallback(Trajectory(Vector{Float64}[], Float64[]), index_map)
    callback = DiscreteCallback(tcb, tcb, save_positions=(false, save_jumps))
    function make_trajectory(prob::DrivenJumpProblem)
        u0 = prob.prob.prob.u0[begin:length(driving_problem.u0)]
        sde_prob = remake(driving_problem, u0=u0, tspan=prob.prob.prob.tspan)
        sol = solve(sde_prob, SOSRI(), saveat=sde_dt, dt=sde_dt)
        Trajectory(sol.u[begin:end-1], sol.t[begin+1:end])
    end
    DrivenJumpProblem(jump_problem, callback, make_trajectory)
end

function SciMLBase.init(prob::DrivenJumpProblem; kwargs...)
    tspan = prob.prob.prob.tspan
    driving_trajectory = prob.make_trajectory(prob)

    prob.callback.condition.traj = driving_trajectory
    prob.callback.condition.index = 1

    # find tstops
    tstops = driving_trajectory.t
    from = searchsortedfirst(tstops, tspan[1])
    to = searchsortedlast(tstops, tspan[2])
    tstops_clipped = @view tstops[from:to]

    integrator = DiffEqJump.init(
        prob.prob,
        SSAStepper();
        callback=prob.callback,
        tstops=tstops_clipped,
        save_start=false,
        kwargs...
    )
    integrator
end

function SciMLBase.solve(prob::DrivenJumpProblem; kwargs...)
    integrator = init(prob; kwargs...)
    solve!(integrator)
    integrator.sol
end

# for remaking
Base.@pure remaker_of(prob::T) where {T<:DrivenJumpProblem} = DiffEqBase.parameterless_type(T)
function DiffEqBase.remake(thing::DrivenJumpProblem; kwargs...)
    T = remaker_of(thing)

    errmesg = """
    DrivenJumpProblems can currently only be remade with new u0, p, tspan or prob fields. To change other fields create a new JumpProblem.
    """
    !issubset(keys(kwargs), (:u0, :p, :tspan, :prob)) && error(errmesg)

    if :prob ∉ keys(kwargs)
        jprob = DiffEqBase.remake(thing.prob; kwargs...)
    else
        any(k -> k in keys(kwargs), (:u0, :p, :tspan)) && error("If remaking a DrivenJumpProblem you can not pass both prob and any of u0, p, or tspan.")
        jprob = kwargs[:prob]
    end

    T(jprob, thing.callback, thing.make_trajectory)
end

Base.summary(io::IO, prob::DrivenJumpProblem) = string(DiffEqBase.parameterless_type(prob), " with problem ", DiffEqBase.parameterless_type(prob.prob))
function Base.show(io::IO, mime::MIME"text/plain", A::DrivenJumpProblem)
    println(io, summary(A))
end

# =======================================
# THE CODE BELOW IS CURRENTLY UNUSED
# =======================================

# We define a new jump aggregator to deal with a jump process that is driven
# by an input signal. The external driving signal is supplied as a trajectory of
# jump times t_i and values u_i. The reaction rates are determined by mass
# action kinetics.

function time_to_next_jump(remaining_input_trajectory)
    t_jump = zero(Float64)
    r = randexp()

    while true
        u_s, t_s = current_values(remaining_input_trajectory)
        u_s_next, t_s_next = iterate(remaining_input_trajectory)
        sumrate = eval_sumrate(u_s, u_x)
        next_sumrate = eval_sumrate(u_s_next, u_x)

        Δt = t_s_next - t_s
        if Δt * (sumrate + next_sumrate) / 2 > r
            # the next jump does not occur before t_s_next
            r -= Δt * (sumrate + next_sumrate) / 2
            continue
        else
            # now we have to find the precise jump time between t_s and t_s_next
            # which is the solution to the equation r = ∫_[t, t+Δt] sumrate(t) dt
            # when solved for Δt
            sumrate_slope = (next_sumrate - sumrate) / Δt
            t_jump = t_s + solve_quadratic(0.5 * sumrate_slope, sumrate, -r)[2]
            break
        end
    end

    t_jump
end

# ================== UTILITY FUNCTIONS ==================

# computes ab-cd in a numerically stable way
@inline function diff_of_products(a, b, c, d)
    w = d * c
    e = fma(-d, c, w)
    f = fma(a, b, -w)
    f + e
end

# computes real roots of quadratic equation ax^2 + bx + c = 0, assuming discriminant is positive
@inline function solve_quadratic(a, b, c)
    discriminant = sqrt(diff_of_products(b, b, 4 * a, c))
    q = -0.5 * (b + copysign(discriminant, b))
    q / a, c / q
end

end


text:
using Pkg, Distributed, DelimitedFiles

const NP = ARGS[1]
const SOLVER = ARGS[2]
const VERBOSE = ARGS[3] == "true"
const QUICK = ARGS[4] == "true"
const GCOFF = ARGS[5] == "true"
const DECODE = ARGS[6] == "true"

addprocs(parse(Int,NP),exeflags="--project=.")
Pkg.instantiate()

if SOLVER == "master"
    Pkg.add(PackageSpec(name="MadNLP",rev="master"))
    Pkg.add(PackageSpec(name="MadNLPHSL",rev="master"))
    Pkg.build("MadNLPHSL")
elseif SOLVER == "current"
    Pkg.develop(path=joinpath(@__DIR__,".."))
    Pkg.develop(path=joinpath(@__DIR__,"..","lib","MadNLPHSL"))
    Pkg.build("MadNLPHSL")
elseif SOLVER == "ipopt"
elseif SOLVER == "knitro"
else
    error("Proper ARGS should be given")
end

# Set verbose option
if SOLVER == "ipopt"
    const PRINT_LEVEL = VERBOSE ? 5 : 0
elseif SOLVER == "knitro"
    const PRINT_LEVEL = VERBOSE ? 3 : 0
else
    using MadNLP
    const PRINT_LEVEL = VERBOSE ? MadNLP.INFO : MadNLP.ERROR
end

# Set quick option


text:
module RanDAGen

"""
Configuration Options

    This package provides the following configuration options for customizing the generated
    DAG's.

Local
- `fatness::Float64`: Defaults to 0.5
- `ccr::Int64`: Defaults to 0
- `data_range::UnitRange{Int64}`: Defaults to 2048:11264
- `alpha_range::UnitRange{Float64}`: Defaults to 0.0:0.2
- `regular::Float64`: Defaults to 0.9 | Must be between 0 and 1
- `jump::Int64`: Defaults to 0
- `density::Float64`: Defaults to 0.5
- `max_task::Int64`: Defaults to 100
- `output_stream::IOStream`: Defaults to stdout
"""
using Distributions

include("statistics.jl")
export Statistics, print_stats

greet() = print("Hello World!")

end # module


text:
using JuMP
using Convex
import MathOptInterface as MOI

"""
Solve the QAP using quadratic integer progamming.
"""
function qap_exact(A::AbstractArray{<:Real,2}, B::AbstractArray{<:Real,2}, obj::Symbol; seeds::Int=0)
	@assert obj in [:min, :max]
	n = size(A,1) - seeds

	if seeds != 0
		A = A[seeds+1:end, seeds+1:end]
		B = B[seeds+1:end, seeds+1:end]
	end

	optimizer = get_optimizer()
	model = Model(optimizer)
	set_silent(model)

	@variable(model, P[1:n,1:n], Bin)

	@constraint(model, sum(P, dims=1) .== 1)
	@constraint(model, sum(P, dims=2) .== 1)

	model_obj = (obj == :min) ? MOI.MIN_SENSE : MOI.MAX_SENSE
	@objective(model, model_obj, tr(A * P * B' * P'))

	optimize!(model)
	status = Int(termination_status(model))

	if status != 1
		error("No optimal solution found")
	end

	P = Int.(value.(P))

	matching = [findfirst(==(1), P[i,:]) for i in 1:n]
	matching = matching .+ seeds
	matching = vcat(1:seeds, matching)

	return P, matching
end

"""
Solve the QAP using quadratic integer progamming. Uses Convex.jl for modeling.
"""
function qap_exact_alt(A::AbstractArray{<:Real,2}, B::AbstractArray{<:Real,2}, obj::Symbol; seeds::Int=0)
	@assert obj in [:min, :max]
	n = size(A,1) - seeds

	if seeds != 0
		A = A[seeds+1:end, seeds+1:end]
		B = B[seeds+1:end, seeds+1:end]
	end

	P = Variable((n, n), :Bin)

	optim = (obj == :min) ? minimize : maximize
	problem = optim(tr(A * P * B' * P'))

	o = ones(Int, n)
	problem.constraints += P * o == o
	problem.constraints += P' * o == o

	optimizer = get_optimizer()
	solve!(problem, optimizer(OutputFlag=0))

	P̂ = Int.(P.value)

	matching = [findfirst(==(1), P̂[i,:]) for i in 1:n]
	matching = matching .+ seeds
	matching = vcat(1:seeds, matching)

	return P̂, matching
end


text:
using Optim, ControlSystems, ControlSystemIdentification

## Generate fake data
w = exp10.(LinRange(-2, 2, 500))
Gtest = tf([1.2], [1.1, 0.0009, 1.2])
data = FRD(w, Gtest)

p0 = (a = [1.0, 1.0, 1.0], b = [1.1]) # Initial guess

Gest = tfest(data, p0, freq_weight = 0)
@test norm(pole(Gtest) - pole(Gest)) < 0.01
@test dcgain(Gtest) ≈ dcgain(Gest)

Gest = tfest(data, tf(p0.b, p0.a), freq_weight = 0)
@test norm(pole(Gtest) - pole(Gest)) < 0.01
@test norm(pole(Gtest) - pole(minimum_phase(Gest))) < 0.01
@test dcgain(Gtest) ≈ dcgain(Gest)

if isinteractive()
    bodeplot(Gtest, w)
    bodeplot!(Gest, w)
    bodeplot!(minimum_phase(Gest), w, l=(:dash,))
end

## Test fitting with basis functions
a = exp10.(LinRange(-1, 1, 7))
poles = ωζ2complex.(a, 0.1)
poles = [poles; conj.(poles)]
# poles = a
# basis = kautz(poles, 1/(200))
basis = laguerre_oo(1, 50)

Gest,p = tfest(data, basis)
r_est = FRD(w, Gest)

@test mean(abs2, log.(abs.(r_est.r)) .- log.(abs.(data.r))) < 0.05

if isinteractive()
    Gr, Gram = baltrunc(Gest)
    Gmp = ControlSystemIdentification.minimum_phase(Gr)
    Gr = baltrunc(Gmp)[1]
    bodeplot(Gtest, w, show=false)
    bodeplot!(Gest, w, show=false)
    bodeplot!(Gr, w, show=false)
    bodeplot!(Gmp, w, show=true)
end

##
if isinteractive()
    t = 0:0.01:100
    u = randn(length(t))
    y,_ = lsim(Gtest, u, t)
    yest,_ = lsim(Gest, u, t)
    ymp,_ = lsim(Gmp, u, t)
    plot([y yest ymp])
end

##
# bands = w
# desired = data.r[1:2:end]
# fi = remez(100, bands, abs.(desired), Hz=200, filter_type=filter_type_hilbert)

## FIR design

# plotly(show=false)

# function DSP.freqz(filter::AbstractVector, w::AbstractVector, h::Real)
#     n = h/length(filter)
#     map(w) do w
#         sum(filter[i]*cis(-w*(i-1)*n) for i in eachindex(filter))
#     end
# end

# function firfit_generic_phase(data::FRD, M, Ts, λ=sqrt(eps()))
#     w = data.w .* Ts
#     maximum(w) ≤ π || throw(ArgumentError("Maximum frequency higher than nyquist"))
#     dvec(w, M) = [cis(-k*w) for k in 0:M-1]
#     Av = [dvec(w, M) for w in w]
#     A = reduce(hcat, Av)'
#     A[:, 1] .= 1
#     A = [A; conj.(A)]
#     Y = [data.r; conj.(data.r)]
#     b = A\Y
#     @show cond(A)
#     @show norm(imag(b))
#     b = real(b)
#     # [reverse(b[2:end]); b]
#     b
# end

# function firfit_linear_phase(data::FRD, M, Ts, λ=sqrt(eps()))
#     w = data.w .* Ts
#     maximum(w) ≤ π || throw(ArgumentError("Maximum frequency higher than nyquist"))
#     # dvec(w, M) = [cis(-k*w) for k in 0:M]
#     dvec(w, M) = [2cos(k*w) for k in 0:(M-1)÷2]
#     Av = [dvec(w, M) for w in w]
#     A = reduce(hcat, Av)'
#     A[:, 1] .= 1
#     A = [A; λ*I]
#     # A = [A; conj.(A)]
#     # Y = [data.r; conj.(data.r)]
#     Y = data.r
#     b = A\[abs.(Y); zeros((M-1)÷2+1)]
#     @show cond(A)
#     b = real(b)
#     [reverse(b[2:end]); b]
# end

# A = map(1:50) do l
#     wl = l*pi/51
#     (1+cos(4wl))/2
# end

# phi = map(1:50) do l
#     wl = l*pi/51
#     -50*wl*(1+0.2*sin(wl))
# end

# w = map(1:50) do l
#     wl = l*pi/51
# end

# data = FRD(w, A .* cis.(phi))
# ##
# w = 2pi .* exp10.(LinRange(-2, 2, 100))
# # w = range(0, pi, length=500)
# # w = 0:0.01:pi
# Gtest = tf([1.2], [1.1, 0.01, 1.2])
# data = FRD(w, Gtest)
# ##
# using PyCall
# sps = pyimport("scipy.signal")
# M = 101
# bs = sps.fir_filter_design.firls(M, data.w, abs.(data.r), fs=1000)

# ##

# # b = firfit_linear_phase(data, M, 0.001, 1)
# b = firfit_generic_phase(data, M, 0.001)
# # bodeplot(tf(b, [1], 1), w, xscale=:identity, yscale=:identity, ticks=:default)
# # bodeplot!(Gtest, w)

# plot([b bs])
# display(current())

# ##
# bodeplot(tf(b, [1], 1), w)
# bodeplot!(tf((bs), [1], 1), w, xscale=:log10, yscale=:log10, ticks=:default)
# display(current())

# ##
# bs = sps.fir_filter_design.firls(101, data.w, abs.(data.r), fs=2pi)
# Gfir = tf((bs), [1], 1)
# bodeplot(Gfir, w, ticks=:default)
# display(current())


text:
# Makes plot for paper showing that transit timing
# precision obeys Brouwer's law:

using PyPlot
using JLD

@load "../../test/test_ttv_cartesian2_4e4days.jld"

using PyPlot

clf()
planet = ["Planet b","Planet c"]
# Plot the difference in the TTVs:
for i=2:3
  diff1 = convert(Array{Float64,1},abs.(tt[i,2:count[i]].-tt_big[i,2:count[i]])./h);
  loglog(tt[i,2:count[i]] .-tt[i,1],diff1,label=planet[i-1]);
#  diff2 = abs.(tt2[i,2:count[i]].-tt2_big[i,2:count[i]])/elements[i,2];
#  loglog(tt[i,2:count[i]]-tt[i,1],diff2);
end
xlabel("Time since start [days]",fontsize=12)
ylabel(L"$abs(t\vert_\mathrm{dbl}-t\vert_\mathrm{big})/h$",fontsize=15)
#loglog([1.0,1024.0],2e-15*[1,2^15],":")
loglog([1.0,40000.0],0.5e-16*([1.0,40000.0]/h).^1.5,":",label=L"$2^{-52}N_S^{3/2}$")
legend()
tight_layout()
savefig("transit_time_errors_4e4.pdf",bbox_inches="tight")
read(stdin,Char)

clf()
nmed = 10
for i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtdq0[i,2:count[i],k,l])) > 0
#    diff1 = abs.(asinh.(dtdq0_big[i,2:count[i],k,l])-asinh.(dtdq0[i,2:count[i],k,l]));
#    diff1 = abs.(dtdq0_big[i,2:count[i],k,l]-dtdq0[i,2:count[i],k,l])/maximum(dtdq0[i,2:count[i],k,l]);
    diff1 = abs.(convert(Array{Float64,1},dtdq0_big[i,2:count[i],k,l]) .-dtdq0[i,2:count[i],k,l])
    diff2 = copy(diff1);
    diff3 = abs.(convert(Array{Float64,1},dtdelements_big[i,2:count[i],k,l]) .-dtdelements[i,2:count[i],k,l])
#    diff3 = abs.(asinh.(dtdelements_big[i,2:count[i],k,l])-asinh.(dtdelements[i,2:count[i],k,l]))
    diff4 = copy(diff3);
    ntt1 = size(diff1)[1];
    for it=nmed:ntt1-nmed
      # Divide by a median smoothed dtdq0:
      diff2[it] = diff1[it]/maximum(abs.(dtdq0[i,it+1-nmed+1:it+1+nmed,k,l]))
      maxit = maximum(abs.(dtdelements[i,it+1-nmed+1:it+1+nmed,k,l]))
      if maxit > 0
        diff4[it] = diff3[it]/maximum(abs.(dtdelements[i,it+1-nmed+1:it+1+nmed,k,l]))
      end
    end;
    diff2[1:nmed-1] .= diff2[nmed]; diff2[ntt1-nmed+1:ntt1] .= diff2[ntt1-nmed];
    diff4[1:nmed-1] .= diff4[nmed]; diff4[ntt1-nmed+1:ntt1] .= diff4[ntt1-nmed];
#    diff3 = abs.(asinh.(dtdq2_big[i,2:count[i],k,l])-asinh.(dtdq2[i,2:count[i],k,l]));
    telapsed = tt[i,2:count[i]] .-tt[i,1]
#    loglog(telapsed,diff2,".",alpha=0.2);
    # Only plot the points which are local maxima:
    indx = collect(1:nmed)
    for i=nmed+1:count[i]-nmed-1
      if diff2[i] == maximum(diff2[i-nmed:i+nmed])
        push!(indx,i)
      end
    end
    push!(indx,count[i]-1)
    loglog(telapsed[indx],diff2[indx]);
#    loglog(telapsed,diff4,linestyle="--");
#    loglog(tt[i,2:count[i]]-tt[i,1],diff3);
    println("xvs: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff2))); #read(STDIN,Char);
    println("els: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff4))); #read(STDIN,Char);
#    println(i," ",k," ",l," asinh error h/2: ",convert(Float64,maximum(diff3))); #read(STDIN,Char);
  end
end
loglog([1.0,40000.0],0.5e-16*([1.0,40000.0]/h).^1.5,linestyle=":")

#axis([1,1024,1e-19,1e-9])
axis([1,31600,1e-16,1e-7])
xlabel("Time since start [days]",fontsize=12)
#ylabel("median(asinh(dt/dq)_dbl - asinh(dt/dq)_big,20)")
ylabel(L"$\frac{maxabs(\frac{dt}{dq}\vert_\mathrm{dbl}-\frac{dt}{dq}\vert_\mathrm{big}[i-9:i+10])}{maxabs(\frac{dt}{dq}\vert_\mathrm{dbl}[i-9:i+10])}$",fontsize=15)
tight_layout()
savefig("Timing_derivative_errors_4e4days_maxabs20.pdf",bbox_inches="tight")

max_stars_repo_name: langfzac/nbg-papers

max_issues_repo_name: langfzac/nbg-papers

max_forks_repo_name: langfzac/nbg-papers


text:
module symbolic_repr
export filler, output, maxHarm, concatFiller,build_candidates
mutable struct filler
    sym::String
    voiced::Bool
    aspirated::Bool
    act::Float64
end

Base.show(io::IO, x::filler) = print(io, "$(x.sym) : $(x.act)")

mutable struct output
    root_initial::filler
    root_final::filler
    suffix::filler
    repr::String
    harmony::Float64
end

Base.show(io::IO, x::output) = print(io, "$(x.repr) : harmony: $(x.harmony) ($(x.root_initial.sym) : $(x.root_initial.act), $(x.root_final.sym) : $(x.root_final.act), $(x.suffix.sym) : $(x.suffix.act))")

"""Join representations of two or three fillers into a string"""
function concatFiller(x::filler, y::filler, z::filler=nothing)
    if z == nothing
        repr = x.sym * y.sym
    else
        repr = x.sym * y.sym * z.sym
    end
    return repr
end

"""Build all possible combinations

Each candidate will be a Vector with the following elements

initial, final, suffix, string_repr, placeholder for harmony

"""
function build_candidates(initial::Vector{filler}, medial::Vector{filler}, suffix::Vector{filler})
    candidates = []
    for c_i ∈ initial
        for c_f ∈ medial
            for c_s ∈ suffix
                cand = output(c_i, c_f, c_s, concatFiller(c_i, c_f, c_s),-Inf)
                push!(candidates, cand)
            end
        end
    end
    return candidates
end

"""Find the candidate with the highest harmony

Return a tuple winner, maxh where

__winner__ : the winner
__maxh__ : harmony of the winner
"""
function maxHarm(candidates::Vector{Any})
    max_h = -Inf
    winner = []
    for el in candidates
        if el.harmony > max_h
            max_h = el.harmony
            winner = el.repr
        end
    end
    return (winner, max_h)
end

"""Decrease gradient activations for wrong winners"""
function decrease_activations(x::output)
    (x.root_initial.act > 0.05 && x.root_initial.act < 1 ) ? x.root_initial.act -= 0.02 : nothing
    (x.root_final.act > 0.05 && x.root_final.act < 1) ? x.root_final.act -= 0.02 : nothing
    (x.suffix.act > 0.05 && x.suffix.act < 1) ? x.suffix.act -= 0.02 : nothing
end

"""Increase activations for the winner"""
function increase_activations(x::output)
    x.root_initial.act < 1 ? x.root_initial.act += 0.02 : nothing
    x.root_final.act < 1 ? x.root_final.act += 0.02 : nothing
    x.suffix.act < 1 ? x.suffix.act += 0.02 : nothing
end

end # module

max_stars_repo_name: amasotti/grascom

max_issues_repo_name: amasotti/grascom

max_forks_repo_name: amasotti/grascom


text:
function Base.convert(Q::Type{<:Quadrangle}, b::Box{2,T}) where {T}
    mincoords = coordinates(b.min)
    maxcoords = coordinates(b.max)
    coords = (
        b.min,
        b.min + Vec(first(maxcoords) - first(mincoords), zero(T)),
        b.max,
        b.min + Vec(zero(T), last(maxcoords) - last(mincoords)),
    )
    Q(coords...)
end

const P2 = Point{2,Int}


text:
#######################################
#
# Extract epochs
#
#######################################

@doc """
Extract epoch data from array of channels.

#### Input

* Array of raw data. Samples x Channels
* Dictionary of trigger information
* Vector of valid trigger numbers
* Number of first triggers to remove
* Number of end triggers to remove

#### Example

```julia
epochs = extract_epochs(data, triggers, [1,2], 0, 0)
```
""" ->
function extract_epochs{T <: Number}(data::Array{T, 2}, triggers::Dict, valid_triggers::Union{AbstractVector, Int},
    remove_first::Int, remove_last::Int, trigger_offset::Int = 252)

    debug("Extracting epochs for $(size(data)[end]) channels using triggers $(valid_triggers)")

    validate_triggers(triggers)

    newTidx = [Int(i) for i in triggers["Index"]]

    triggers = DataFrame(Code = triggers["Code"], Index = newTidx)
    # TODO Use convert function
    #=triggers = convert(DataFrame, triggers)=#

    # Change offset so numbers are manageable
    triggers[:Code] = triggers[:Code] - trigger_offset

    # Determine indices of triggers which are valid
    valid_triggers = any(triggers[:Code] .== valid_triggers', 2)
    valid_triggers = find(valid_triggers .== true)

    debug("Number of valid triggers: $(length(valid_triggers))")

    # Remove unwanted triggers
    triggers = triggers[valid_triggers, :]                # That aren't valid
    triggers = triggers[remove_first+1:end-remove_last,:] # Often the first trigger is rubbish

    lenEpochs = minimum(diff(triggers[:Index]))
    numChans  = size(data)[end]

    # Check we aren't looking past the end of the data
    start_indices = convert(Array{Int}, triggers[:Index])
    end_indices   = convert(Array{Int}, start_indices + lenEpochs - 1)
    while end_indices[end] > size(data, 1)
        pop!(start_indices)
        pop!(end_indices)
        debug("Removed end epoch as its not complete")
    end

    # Create variable for epochs
    numEpochs = length(start_indices)
    epochs = zeros(Float64, (Int(lenEpochs), Int(numEpochs), Int(numChans)))

    # User feedback
    debug("Creating epochs: $lenEpochs x $numEpochs x $numChans")

    for si = 1:length(start_indices)
        epochs[:, si, :] = data[start_indices[si]: end_indices[si], :]
    end

    Logging.info("Generated $numEpochs epochs of length $lenEpochs for $numChans channels")

    return epochs
end

#######################################
#
# Create average epochs
#
#######################################

function average_epochs(ep::Array)

    Logging.info("Averaging down epochs to 1 epoch of length $(size(ep,1)) from $(size(ep,2)) epochs on $(size(ep,3)) channels")

    squeeze(mean(ep, 2), 2)
end


text:
struct SortGroups <: Constraint

    UpGroup::Variable
    DownGroup::Variable

end

function SortGroups(model, UpGroup, DownGroup)
    constraint = SortGroups(UpGroup, DownGroup)

    push!(UpGroup.linkedConstraint, constraint)
    push!(DownGroup.linkedConstraint, constraint)

    push!(model.constraints, constraint)
end

function filtrage!(groups::SortGroups)

    changeVariable = Vector{Variable}(undef, 3)
    nbChange = 0

    changeUp = Change(groups.UpGroup)
    changeDown = Change(groups.DownGroup)

    #### filtering DownGroup

    if !groups.DownGroup.isFixed
        # v : minimum element of UpGroup.upperBound
        v = minimum(groups.UpGroup.upperBound)

        c = 0
        n_UB = copy(groups.DownGroup.upperBound) # new Upper Bound

        changed = false
        for e in 1:length(groups.DownGroup.upperBound)
            if groups.DownGroup.upperBound[e] <= v # élément plus petit que le plus petit élément possible d'un groupe au dessus => ne peut pas être placé
                deleteat!(n_UB, e-c)
                c += 1

                push!(changeDown.removed, groups.DownGroup.upperBound[e])
                changed = true
            end
        end

        if changed
            nbChange += 1
            changeVariable[nbChange] = groups.DownGroup
        end

        groups.DownGroup.upperBound = n_UB

        cardSup = min(groups.DownGroup.cardinalSup, length(groups.DownGroup.upperBound))
        changeDown.cardRemoved = cardSup - groups.DownGroup.cardinalSup

        groups.DownGroup.cardinalSup += changeDown.cardRemoved
    end

    ####################################################################################

    #### filtering UpGroup
    if length(groups.DownGroup.lowerBound) > 0 && length(groups.UpGroup.lowerBound) > 0
        v = minimum(groups.DownGroup.lowerBound)
        w = minimum(groups.UpGroup.lowerBound)

        if !groups.UpGroup.isFixed && groups.UpGroup.cardinalSup - length(groups.UpGroup.lowerBound) == 1 && v < w # il ne reste plus qu'un élément à placer

            c = 0
            n_UB = copy(groups.UpGroup.upperBound)
            changed = false
            for e in 1:length(groups.UpGroup.upperBound)
                if groups.UpGroup.upperBound[e] > v && !(groups.UpGroup.upperBound[e] in groups.UpGroup.lowerBound)
                    deleteat!(n_UB, e-c)
                    c += 1

                    push!(changeUp.removed, groups.UpGroup.upperBound[e])
                    changed = true
                end
            end

            if changed
                nbChange += 1
                changeVariable[nbChange] = groups.UpGroup
            end

            groups.UpGroup.upperBound = n_UB

            cardSup = min(groups.UpGroup.cardinalSup, length(groups.DownGroup.upperBound))
            changeUp.cardRemoved = cardSup - groups.UpGroup.cardinalSup

            groups.UpGroup.cardinalSup += changeUp.cardRemoved
        end
    end

    ####################################################################################

    stop = groups.UpGroup.cardinalInf > groups.UpGroup.cardinalSup
	stop |= groups.DownGroup.cardinalInf > groups.DownGroup.cardinalSup

    return changeVariable[1:nbChange], (changeUp, changeDown), stop
end

function quoi()
    v1 = Variable("v1", [8,10,11], [3,6,8,9,10,11,13], 4, 4, Set([]), false)
    v2 = Variable("v2", [7,14,16], [5,7,12,14,15,16], 4, 4, Set([]), false)
    println(v1.upperBound)
    println(v2.upperBound)
    model = ModelTest()
    sg = SortGroups(v1,v2)
    println(sg)
    filtrage!(sg)
    println(v1.upperBound)
    println(v2.upperBound)
end


text:
const ⊗ = kron

"""
Build Laplacian matrix using Kronecker products
"""
function get_Laplacian3d_kron( LF::LF3dGrid )
  IIx = speye(LF.Nx)
  IIy = speye(LF.Ny)
  IIz = speye(LF.Nz)
  return LF.LFx.D2jl ⊗ IIy ⊗ IIz +
         IIx ⊗ LF.LFy.D2jl ⊗ IIz +
         IIx ⊗ IIy ⊗ LF.LFz.D2jl
end

"""
A naive implementation to build Laplacian by manually, by searching for non-zero
indices and their values.
"""
function get_Laplacian3d_naive( LF::LF3dGrid )
  N = LF.Nx * LF.Ny * LF.Nz

  NNZ = calc_NNZ(LF)
  Tnz = [false, false, false]
  # Nabla2 matrix
  values = zeros(Float64,NNZ)
  column = zeros(Int64,NNZ)
  row    = zeros(Int64,NNZ)
  inz    = 0
  for ip2 = 1:N^3
    for ip1 = ip2:N^3
      Tnz[:] = false
      i1 = LF.lin2xyz[1,ip1]
      i2 = LF.lin2xyz[1,ip2]
      #
      j1 = LF.lin2xyz[2,ip1]
      j2 = LF.lin2xyz[2,ip2]
      #
      k1 = LF.lin2xyz[3,ip1]
      k2 = LF.lin2xyz[3,ip2]
      #
      nabla2 = 0.0
      #
      if j1 == j2 && k1 == k2
        nabla2 = nabla2 + LF.LFx.D2jl[i1,i2]
        Tnz[1] = true
      end
      #
      if i1 == i2 && k1 == k2
        nabla2 = nabla2 + LF.LFy.D2jl[j1,j2]
        Tnz[2] = true
      end
      #
      if i1 == i2 && j1 == j2
        nabla2 = nabla2 + LF.LFz.D2jl[k1,k2]
        Tnz[3] = true
      end
      #
      if any(Tnz) # update non-zero element
        inz = inz + 1
        # Notice that we do not exploit symmetry here
        values[inz] = nabla2
        column[inz] = ip1
        row[inz]    = ip2
        #@printf("%5d %5d: Non zero found, inz = %8d\n", ip1, ip2, )
        if ip1 != ip2 # not a diagonal element
          inz = inz + 1 # update NNZ again
          values[inz] = nabla2
          column[inz] = ip2
          row[inz]    = ip1
        end
      end # if Tnz
    end # for ip1
  end # for ip2
  return sparse(row, column, values)
end

function calc_NNZ( LF::LF3dGrid ; verbose=false )
  Npoints = LF.Nx * LF.Ny * LF.Nz
  #
  NNZ = 0
  Tnz = [false, false, false]
  # Nabla2 matrix
  for ip2 = 1:Npoints
    for ip1 = ip2:Npoints
      Tnz[:] = false
      i1 = LF.lin2xyz[1,ip1]
      i2 = LF.lin2xyz[1,ip2]
      #
      j1 = LF.lin2xyz[2,ip1]
      j2 = LF.lin2xyz[2,ip2]
      #
      k1 = LF.lin2xyz[3,ip1]
      k2 = LF.lin2xyz[3,ip2]
      #
      nabla2 = 0.0
      #
      if (j1 == j2 && k1 == k2) || (i1 == i2 && k1 == k2) || (i1 == i2 && j1 == j2)
        NNZ = NNZ + 1
        #@printf("%5d %5d: Non zero found, inz = %8d\n", ip1, ip2, NNZ)
        if ip1 != ip2 # not a diagonal element
          NNZ = NNZ + 1
          #@printf("%5d %5d: Non zero found, inz = %8d\n", ip1, ip2, NNZ)
        end
      end
    end # ip1
  end # ip2
  if verbose
    @printf("NNZ = %8d, nonzero percentage = %7.2f%%\n", NNZ, NNZ/(Npoints^2) * 100.0 )
  end
  return NNZ
end

function old_calc_NNZ( LF::LF3dGrid ; verbose=false )
  Npoints = LF.Nx * LF.Ny * LF.Nz
  #
  NNZ = 0
  Tnz = [false, false, false]
  # Nabla2 matrix
  for ip2 = 1:Npoints
    for ip1 = ip2:Npoints
      Tnz[:] = false
      i1 = LF.lin2xyz[1,ip1]
      i2 = LF.lin2xyz[1,ip2]
      #
      j1 = LF.lin2xyz[2,ip1]
      j2 = LF.lin2xyz[2,ip2]
      #
      k1 = LF.lin2xyz[3,ip1]
      k2 = LF.lin2xyz[3,ip2]
      #
      nabla2 = 0.0
      #
      if j1 == j2 && k1 == k2
        Tnz[1] = true
      end
      #
      if i1 == i2 && k1 == k2
        Tnz[2] = true
      end
      #
      if i1 == i2 && j1 == j2
        Tnz[3] = true
      end
      if any(Tnz)
        NNZ = NNZ + 1
        #@printf("%5d %5d: Non zero found, inz = %8d\n", ip1, ip2, NNZ)
        if ip1 != ip2 # not a diagonal element
          NNZ = NNZ + 1
          #@printf("%5d %5d: Non zero found, inz = %8d\n", ip1, ip2, NNZ)
        end
      end
    end # ip1
  end # ip2
  if verbose
    @printf("NNZ = %8d, nonzero percentage = %7.2f%%\n", NNZ, NNZ/(Npoints^2) * 100.0 )
  end
  return NNZ
end


text:
# Just the bits of Lazy.jl that we actually use.

macro errs(ex)
  :(try $(esc(ex))
    catch e
      showerror(stderr, e, catch_backtrace())
      println(stderr)
    end)
end

d(xs...) = Dict{Any, Any}(xs...)


text:
module RunTime

include("DebugInfo.jl")

const kRunTimeOut = Channel{AbstractString}(32)
const kRunTimeIn = Channel{AbstractString}(33)

mutable struct BreakPoints
  filepath::AbstractString
  lineno::Array{Int64,1}
end

mutable struct AstInfo
  asts::Array
  blocks::Array # save struct BlockInfo
end

mutable struct BlockInfo
  BlockType::Symbol
  startline::Int64
  endline::Int64
  raw_code::AbstractString
end

# FileLine["FileName"] = line number
FileLine = Dict()
# FileAst["filename"] = struct AstInfo
FileAst = Dict()
# FileBp["filename"] = [bp line array]
FileBp = Dict()
# run file stack for `include`
RunFileStack = []
RunBlockStatus = []

EntryFile = ""
errors = ""

struct NotImplementedError <: Exception end
struct BreakPointStop <: Exception end

# entry file for debugging
function setEntryFile(file)
  global RunFileStack
  global FileLine
  global EntryFile
  EntryFile = abspath(file)
  push!(RunFileStack, EntryFile)
  if !haskey(FileAst, EntryFile)
    readSourceToAST(EntryFile)
  end
end

# clear runblockstatus, this should be called
# after collect info for a block
function clearBlockStatus()
  global RunBlockStatus
  RunBlockStatus = []
end

function readSourceToAST(file)
  global FileAst
  global FileLine
  asts = []
  blocks = []
  line = 0
  s = ""
  # recording blocks such as a function
  block_start = 0
  block_end = 0
  open(file) do f
    for ln in eachline(f)
      line += 1
      s *= ln
      ex = parseInputLine(s)
      if (isa(ex, Expr) && ex.head === :incomplete)
        s *= "\n"
        if block_start == 0
          block_start = line
        end
        continue
      else
        if block_start != 0
          # add block line postion
          block_end = line
          s = replace(s, "\r" => "")
          blockinfo = BlockInfo(ex.head, block_start, block_end, s)
          push!(blocks, blockinfo)
          block_start = 0
        end
        ast = Meta.lower(Main, ex)
        push!(asts, ast)
        s = ""
      end
    end
  end
  FileAst[file] = AstInfo(asts, blocks)
  if !haskey(FileLine, file)
    FileLine[file] = 1
  end
end

# run whole program
# return a status contains ast, current_line
function run()
  global FileLine
  global FileAst
  global RunFileStack
  global EntryFile
  global kRunTimeOut
  # wait until 'launch'
  while take!(kRunTimeIn) != "launch" end
  # reset all file line
  RunFileStack = []
  push!(RunFileStack, EntryFile)
  cd(abspath(joinpath(EntryFile, "../")))
  FileLine = Dict(EntryFile => 1)
  # start running program
  while true
    cmd = take!(kRunTimeIn)
    notFinish = true
    if cmd == "continue"
      notFinish = continous()
    elseif cmd == "stepOver"
      notFinish = stepOver()
    end
    if !notFinish
      break
    end
  end
  return
end

# check whether is `include` call
# if it is & included file has breakpoints,
# we manually run it
function tryRunNewFile(ast, isStepOver = false)
  global RunFileStack
  local isIncludeCall = false
  local filename = ""
  try
    if ast.args[1].code[1].args[1] == (:include)
      filename = abspath(joinpath(RunFileStack[end], "../" * ast.args[1].code[1].args[2]))
      isIncludeCall = true
    end
  catch err
    # we may get bounderror
    return false
  end
  if isIncludeCall
    if !haskey(FileAst, filename)
      readSourceToAST(filename)
    end
    # skip include call (we run it manually)
    updateLine()
    push!(RunFileStack, filename)
    FileLine[filename] = 1
    continous(isStepOver)
  end
  return isIncludeCall
end

# go to next line
# if next line is an empty line,
# we ignore this line and go until it's a valid ast
function stepOver()
  global FileAst
  global FileLine
  global RunFileStack
  current_file = RunFileStack[end]
  asts = FileAst[current_file].asts
  ast_index = getAstIndex(current_file, FileLine[current_file])
  # run next line code
  if asts[ast_index] isa Nothing
    # if we find a blank line, we skip it
    updateLine()
    stepOver()
  end
  # true line
  try
    if !tryRunNewFile(asts[ast_index], true)
      Core.eval(Main, asts[ast_index])
      updateLine()
    else
      throw(BreakPointStop(""))
    end
  catch err
    if err isa BreakPointStop
      Break()
      return true
    else
      global errors
      errors = string(err)
      println(errors)
      return false
    end
  end
  if getAstIndex(current_file, FileLine[current_file]) > lastindex(asts)
    pop!(RunFileStack)
    if length(RunFileStack) == 0
      put!(kRunTimeOut, "finish")
      return false
    end
  end
  Break()
  return true
end

# go until we meet a bp
# if stopOnPopFile is true: we return (withou Break())
# when a file is finished executing
function continous(stopOnPopFile = false)
  global FileAst
  global FileLine
  global RunFileStack
  global kRunTimeOut
  if length(RunFileStack) == 0
    put!(kRunTimeOut, "finish")
    return false
  end
  current_file = RunFileStack[end]
  asts = FileAst[current_file].asts

  ast_index = getAstIndex(current_file, FileLine[current_file])
  while ast_index <= lastindex(asts)
    ast = asts[ast_index]
    try
      if !tryRunNewFile(ast)
        Core.eval(Main, ast)
        updateLine()
      else
        return true
      end
    catch err
      if err isa BreakPointStop
        Break()
        return true
      else
        global errors
        errors = "runtime errors: $(current_file): $(FileLine[current_file]): $(err)"
        println(errors)
        put!(kRunTimeOut, "error")
        return
      end
    end
    ast_index = getAstIndex(current_file, FileLine[current_file])
  end
  # exit normally
  pop!(RunFileStack)
  if !stopOnPopFile
    continous()
  end
  # since only stepOver will make stopOnPopFile true
  # we don't need put!(kRunTimeOut) here
end

function stepIn()
  throw(NotImplementedError("step has not been implemented"))
end

# update info from this point
function Break()
  # debug info
  global RunFileStack
  global FileLine
  global kRunTimeIn
  global kRunTimeOut
  current_line = FileLine[RunFileStack[end]]
  # println("hit breakpoint: $(RunFileStack[end]): $(current_line)")
  # collect variable info
  DebugInfo.collectVarInfo()
  # collect stack info
  DebugInfo.collectStackInfo()
  put!(kRunTimeOut, "collected")
  # wait until get "go on" info
  while take!(kRunTimeIn) != "go on" end
end

function inBlockBreak(file, lineno)
  global RunFileStack
  global FileLine
  global kRunTimeIn
  global kRunTimeOut
  # println("hit breakpoint in block: $(file): $(lineno)")
  # collect variable info
  DebugInfo.collectVarInfo()
  # collect stack info
  DebugInfo.collectStackInfo()
  # set block running info
  push!(RunBlockStatus, Dict(file => lineno))
  put!(kRunTimeOut, "collected")
  # wait until get "go on" info
  while take!(kRunTimeIn) != "go on" end
end

# set breakpoint
# if breakpoint is on blank line,
# it will not be set
function setBreakPoints(filepath, lineno)
  global FileBp
  global FileAst
  result = []
  filepath = abspath(filepath)
  readSourceToAST(filepath)
  FileBp[filepath] = []
  asts = FileAst[filepath].asts
  blocks = FileAst[filepath].blocks
  for line in lineno
    ast_idx = getAstIndex(filepath, line)
    # outside empty line
    if asts[ast_idx] isa Nothing
      push!(result, false)
    else
      ast, res = trySetBpInBlock(filepath, line)
      if res == 0 || res == 1
        push!(FileBp[filepath], line)
        push!(result, true)
        # set modified ast to filepath
        if res == 1
          FileAst[filepath].asts[ast_idx] = ast
        end
      elseif res == 2
        push!(result, false)
      end
    end
  end
  return result
end

# check and set breakpoint inside blocks
# return: 0: line not in block
#         1: line in block and valid
#         2: line in block but not valid
function trySetBpInBlock(filepath, line)
  global FileAst
  blocks = FileAst[filepath].blocks
  # for windows dir
  input_path = replace(filepath, "\\" => "\\\\")
  result = 0
  block_idx = 1
  for block in blocks
    if line <= block.startline
      # means it is not in a block
      # setting bp on the first line of block
      # will only be stopped once
      break
    end
    if block.startline < line <= block.endline
      # means it is in a block
      modified_code = ""
      ofs = 0
      for code_line in split(block.raw_code, "\n")
        if line == block.startline + ofs
          # means we detect the position to set bp
          try
            if Meta.parse(code_line) isa Nothing
              # means bp are set at an empty line, which is not allowed
              result = 2
            else
              result = 1
              code_line = "Judy.EventHandler.RunTime.inBlockBreak(\"$(input_path)\", $(line));" * code_line
            end
          catch err
            # some errors may cause when try to parse like `else`
            result = 2
          end
        end
        modified_code = modified_code * code_line * " \n"
        # save modified raw code since we can have multiple bps in same block
        FileAst[filepath].blocks[block_idx].raw_code = modified_code
        ofs += 1
      end
      ast = Meta.lower(Main, parseInputLine(modified_code))
      return ast, result
    end
    block_idx += 1
  end
  ast = Meta.lower(Main, parseInputLine(""))
  return ast, result
end

# update line for run/next/continous call
function updateLine()
  global RunFileStack
  global FileLine
  global FileAst
  current_file = RunFileStack[end]
  blocks = FileAst[current_file].blocks
  ofs = 1
  for blockinfo in blocks
    if FileLine[current_file] == blockinfo.startline
      ofs += blockinfo.endline - blockinfo.startline + 1
      break
    elseif blockinfo.startline > FileLine[current_file]
      break
    end
  end
  ofs = ofs > 1 ? ofs - 1 : ofs
  FileLine[current_file] += ofs
  # continue update until we meet a not empty line
  idx = getAstIndex(current_file, FileLine[current_file])
  if idx <= lastindex(FileAst[current_file].asts)
    if FileAst[current_file].asts[idx] isa Nothing
      updateLine()
    end
  end
  if checkBreakPoint()
    throw(BreakPointStop())
  end
end

# check whether current line meets a bp
function checkBreakPoint()
  global RunFileStack
  global FileLine
  global FileBp
  current_file = RunFileStack[end]
  current_line = FileLine[RunFileStack[end]]
  if !haskey(FileBp, current_file)
    return false
  end
  if current_line in FileBp[current_file]
    return true
  else
    return false
  end
end

# get AstIndex from current file and line number
function getAstIndex(filepath, lineno)
  global FileAst
  blocks = FileAst[filepath].blocks
  base = lineno
  ofs = 0
  for block in blocks
    if block.startline <= lineno <= block.endline
      base = block.startline
      break
    elseif block.startline > lineno
      break
    end
    ofs += block.endline - block.startline
  end
  return base - ofs
end

function parseInputLine(s::String; filename::String="none", depwarn=true)
  # For now, assume all parser warnings are depwarns
  ex = if depwarn
      ccall(:jl_parse_input_line, Any, (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Csize_t),
            s, sizeof(s), filename, sizeof(filename))
  else
      with_logger(NullLogger()) do
          ccall(:jl_parse_input_line, Any, (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Csize_t),
                s, sizeof(s), filename, sizeof(filename))
      end
  end
  return ex
end

end # module RunTime

max_stars_repo_name: judy-vscode/judy

max_issues_repo_name: judy-vscode/judy

max_forks_repo_name: judy-vscode/judy


text:
function download_dep(orig, dest, hash)
    dest_file = joinpath("data", dest)
    if isfile(dest_file)
        dest_hash = open(dest_file, "r") do f
            bytes2hex(sha256(f))
        end
        if dest_hash == hash
            return nothing
        end
    end
    mkpath("data")
    download(orig, dest_file)
    return nothing
end

@testset "reproject-core" begin
    download_dep("https://astropy.stsci.edu/data/galactic_center/gc_2mass_k.fits", "gc_2mass_k.fits",
                 "763ef344df3ac8fa80ff46f00ca1ec59946ca3f99502562d6fcfb73320b1cec3")
    download_dep("https://astropy.stsci.edu/data/galactic_center/gc_msx_e.fits", "gc_msx_e.fits",
                 "3687fb3763911825f981e74b6a9b82c0e618f7e592b1e0cb17e2c63164e28cd6")

    imgin = FITS(joinpath("data", "gc_msx_e.fits"))    # project this
    imgout = FITS(joinpath("data", "gc_2mass_k.fits")) # into this coordinate

    hdu1 = astropy.io.fits.open(joinpath("data", "gc_2mass_k.fits"))[1]
    hdu2 = astropy.io.fits.open(joinpath("data", "gc_msx_e.fits"))[1]

    @test isapprox(reproject(imgin, imgout, order = 0)[1]', rp.reproject_interp(hdu2, hdu1.header, order = 0)[1], nans = true, rtol = 1e-7)
    @test isapprox(reproject(imgout, imgin, order = 0)[1]', rp.reproject_interp(hdu1, hdu2.header, order = 0)[1], nans = true, rtol = 1e-6)
    @test isapprox(reproject(imgin, imgout, order = 1)[1]', rp.reproject_interp(hdu2, hdu1.header, order = 1)[1], nans = true, rtol = 1e-7)
    @test isapprox(reproject(imgin, imgout, order = 2)[1]', rp.reproject_interp(hdu2, hdu1.header, order = 2)[1], nans = true, rtol = 6e-2)
    @test isapprox(reproject(imgin[1], imgout[1], shape_out = (1000,1000))[1]',
            rp.reproject_interp(hdu2, astropy.wcs.WCS(hdu1.header), shape_out = (1000,1000))[1], nans = true, rtol = 1e-7)

    wcs = WCSTransform(2; ctype = ["RA---AIR", "DEC--AIR"], radesys = "UNK")
    @test_throws ArgumentError reproject(imgin, wcs, shape_out = (100,100))

    fname = tempname() * ".fits"
    f = FITS(fname, "w")
    inhdr = FITSHeader(["CTYPE1", "CTYPE2", "RADESYS", "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
                        "HISTORY"],
                       ["RA---TAN", "DEC--TAN", "UNK", 1.0, 1, true, "string value", nothing, nothing],
                       ["",
                        "",
                        "",
                        "floating point keyword",
                        "",
                        "boolean keyword",
                        "string value",
                        "this is a comment",
                        "this is a history"])

    indata = reshape(Float32[1:100;], 5, 20)
    write(f, indata; header=inhdr)
    @test_throws ArgumentError reproject(f[1], imgin, shape_out = (100,100))
    close(f)
end


text:
module CUSPARSE

# using
# exports
# imports

# find the library
const libcusparse = find_library(["libcusparse"])
isempty(libcusparse) && error("CUSPARSE library cannot be found")

# error handling
const cusparseErrorString = ["Not initialized", "Alloc failed", "Invalid value", "Arch mismatch", "Mapping error", "Execution failed", "Internal error", "Matrix type not supported", "Zero pivot"]

function cusparseCheck(status)
    if status == CUSPARSE_STATUS_SUCCESS
        return nothing
    end
    # Because try/finally may disguise the source of the problem,
    # let's show a backtrace here
    warn("CUSPARSE error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    println()
    throw(cusparseErrorString[status])
end

# These are semi-automatically generated using Clang with wrap_cusparse.jl:
include("types.jl")
include("libcusparse.jl")

# Setup the default cusparse handle
cusparseHandlePtr = cusparseHandle_t[0]
cusparseCreate(cusparseHandlePtr)
cusparseHandle = cusparseHandlePtr[1]
# destroy cusparse handle at julia exit
atexit(()->cusparseDestroy(cusparseHandle))

# Set the default pointer mode to host
cusparseSetPointerMode(cusparseHandle,CUSPARSE_POINTER_MODE_HOST)

# Setup the default MatDescr with index base = 1
cusparseMatDescrPtr = cusparseMatDescr_t[0]
cusparseCreateMatDescr(cusparseMatDescrPtr)
cusparseMatDescrDefault = cusparseMatDescrPtr[1]
cusparseSetMatType(cusparseMatDescrDefault, CUSPARSE_MATRIX_TYPE_GENERAL)
cusparseSetMatIndexBase(cusparseMatDescrDefault, CUSPARSE_INDEX_BASE_ONE)

end # module


text:
"Initialize LIKWIDs topology module."
function init_topology()
    ret = LibLikwid.topology_init()
    if ret == 0
        _cputopo[] = unsafe_load(LibLikwid.get_cpuTopology())
        _cpuinfo[] = unsafe_load(LibLikwid.get_cpuInfo())
        _build_jl_cputopo()
        _build_jl_cpuinfo()
        topo_initialized[] = true
        return true
    end
    return false
end

function _build_jl_cputopo()
    ct = _cputopo[]
    ncachelvls = Int(ct.numCacheLevels)
    nhwthreads = Int(ct.numHWThreads)
    threadpools = unsafe_wrap(Array, ct.threadPool, nhwthreads)
    threads = Vector{HWThread}(undef, nhwthreads)
    for (i, tp) in enumerate(threadpools)
        threads[i] = HWThread(
            tp.threadId,
            tp.coreId,
            tp.packageId,
            tp.apicId,
            tp.dieId,
            tp.inCpuSet,
        )
    end
    cachelvls = unsafe_wrap(Array, ct.cacheLevels, ncachelvls)
    caches = Vector{CacheLevel}(undef, ncachelvls)
    for (i, clvl) in enumerate(cachelvls)
        caches[i] = CacheLevel(
            clvl.level,
            clvl.type == LibLikwid.DATACACHE ? :data :
                          clvl.type == LibLikwid.INSTRUCTIONCACHE ? :instruction :
                          clvl.type == LibLikwid.UNIFIEDCACHE ? :unified :
                          clvl.type == LibLikwid.ITLB ? :itlb :
                          clvl.type == LibLikwid.DTLB ? :dtlb : :nocache,
            clvl.associativity,
            clvl.sets,
            clvl.lineSize,
            clvl.size,
            clvl.threads,
            clvl.inclusive,
        )
    end

    cputopo[] = CpuTopology(
        nhwthreads,
        ct.activeHWThreads,
        ct.numSockets,
        ct.numDies,
        ct.numCoresPerSocket,
        ct.numThreadsPerCore,
        ncachelvls,
        threads,
        caches,
    )
    return nothing
end

"Close and finalize LIKWIDs topology module."
function finalize_topology()
    LibLikwid.topology_finalize()
    topo_initialized[] = false
    _cputopo[] = nothing
    _cpuinfo[] = nothing
    cputopo[] = nothing
    cpuinfo[] = nothing
    return nothing
end

function _build_jl_cpuinfo()
    ci = _cpuinfo[]
    cpuinfo[] = CpuInfo(
        ci.family,
        ci.model,
        ci.stepping,
        ci.vendor,
        ci.part,
        ci.clock,
        ci.turbo,
        unsafe_string(ci.osname),
        unsafe_string(ci.name),
        unsafe_string(ci.short_name),
        unsafe_string(ci.features),
        ci.isIntel,
        join(Char(c) for c in ci.architecture if !iszero(c)),
        ci.supportUncore,
        ci.supportClientmem,
        ci.featureFlags,
        ci.perf_version,
        ci.perf_num_ctr,
        ci.perf_width_ctr,
        ci.perf_num_fixed_ctr,
    )
    return nothing
end

"""
    get_cpu_topology() -> CpuTopology
Get the CPU topology of the machine.

Automatically initializes the topology and NUMA modules,
i.e. calls [`LIKWID.init_topology`](@ref) and [`LIKWID.init_numa`](@ref).
"""
function get_cpu_topology()
    if !topo_initialized[]
        init_topology() || error("Couldn't init topology.")
    end
    if !numa_initialized[]
        init_numa() || error("Couldn't init numa.")
    end
    return cputopo[]
end

"""
    get_cpu_info() -> CpuInfo
Get detailed information about the CPU.

Automatically initializes the topology and NUMA modules,
i.e. calls [`LIKWID.init_topology`](@ref) and [`LIKWID.init_numa`](@ref).
"""
function get_cpu_info()
    if !topo_initialized[]
        init_topology() || error("Couldn't init topology.")
    end
    if !numa_initialized[]
        init_numa() || error("Couldn't init numa.")
    end
    return cpuinfo[]
end

"""
    print_supported_cpus(; cprint=true)
Print a list of all supported CPUs.

If `cprint=false`, LIKWID.jl
will first capture the stdout and then `print` the list.
"""
function print_supported_cpus(; cprint=true)
    if cprint
        LibLikwid.print_supportedCPUs()
    else
        buf = IOBuffer()
        capture_stdout!(LibLikwid.print_supportedCPUs, buf)
        s = String(take!(buf))
        print(s)
    end
    return nothing
end

"""
Graphical visualization of the CPU topology. Extracts the corresponding output of `likwid-topology -g`.
"""
function print_cpu_topology()
    out = _execute(`likwid-topology -g`)
    if out[:stderr] != ""
        @warn out[:stderr]
    end
    stdout = out[:stdout]
    idcs = findfirst("Graphical Topology", stdout)
    print(stdout[first(idcs):end])
    return nothing
end


text:
# include("diffUtils.jl")

"""
	2D IMPLEMENTATION: f is search area, g is interrogation area (translated function)
"""

#  SQERR =         ( IA - SA )²         / ( sqrt( sum(IA²) )*sqrt( sum(SA²) ) )
#        = ( IA²  +   SA²  -  2*IA*SA ) / ( sqrt( sum(IA²) )*sqrt( sum(SA²) ) )
#           -----    -----     -------           --------         ---------
#           const   intArr  FFTcross(IA,SA)        const            intArr

function crossCorrelation!( ::NSQECC, cmat::A{T,2}, shifted::A{T,2}, shifts,
                                    padf::A{C{T},2}, padg::A{C{T},2}, plan, iplan,
                                    intArr2::A{T,2}, sumG2, sizeF, sizeG
						  ) where {T<:AbstractFloat}

	# Computing IA*SA
    crossCorrelation!( FFT(), cmat, padf, padg, plan, iplan );
    Base.circshift!( shifted, cmat, shifts[1:2] )

	# Computing sqrt(sum(IA²))
	denG = sqrt( sumG2 )

	# Computing sum(SA²) and sqrt(sum(SA²)).
	fh, fw = sizeF
	gh, gw = sizeG .- 1

	for c in 0:size(cmat,2)-1; 		c1 = max(1, fw - c ); c2 = min( fw, c1 + gw )
		for r in 0:size(cmat,1)-1;  r1 = max(1, fh - r ); r2 = min( fh, r1 + gh )

			# sometimes the precision of float operations produces small negative values.
			# We take absolute value to compensate for this, which does not really change
			# the position of the maximum peak, but it removes sqrt( negValue ) errors.
			sumF2 = abs( integralArea( intArr2, (r1-1,c1-1), (r2,c2) ) )
			num   = sumF2 + sumG2 - 2*shifted[r+1,c+1]

			cmat[r+1,c+1] = 1 / ( 1 + num/( denG * sqrt(sumF2) ) )
	end	end
end

function PIV_2D( ::NSQECC, img1::A{<:Real,2}, img2::A{<:Real,2},
                         IA::III, SM::III, overlap::III, mpass::I, width::I,
                         peak::S, sigNoise::S, filtFun::Function, threshold::F ;
                         corrType=Float32, vfType=Float32 )

	# Calculating size of the vector field
    step   = IA .- overlap;
	imsize = ( size(img1,1), size(img1,2), size(img1,3) )
    VFsize = length.( StepRange.( IA, step, imsize ) );

    U  = zeros( vfType, VFsize[1:2] );
    V  = zeros( vfType, VFsize[1:2] );
    SN = zeros( vfType, VFsize[1:2] );

	ignoreSN = ( sigNoise == "" );
	SM = ( SM[1], SM[2], 0 );

    for mp in mpass:-1:1

        last_mp = ( mp == 1 );

        # Scaling IA, SM and step (overlap) to multi-pass iteration.
        IA_mp    =  IA  .* mp;
        SM_mp    =  SM  .* mp;
        step_mp  = step .* mp;
		IAranges = StepRange.( IA_mp, step_mp, imsize )

        # Initialize cross-correlation variables
        csize   = 2 .* ( IA_mp .+ SM_mp  );
        shifts  = div.( csize, 2 ) .+ SM_mp .- 1;
        cmatrix = zeros( corrType, csize[1:2] );
        shifted = zeros( corrType, csize[1:2] );
        pads    = zeros( Complex{corrType}, csize[1:2] );
        padi    = zeros( Complex{corrType}, csize[1:2] );
        plan    =  FFTW.plan_fft!( pads );
        iplan   = FFTW.plan_ifft!( pads );

		# Initializing integral array for search area
        ssize   = IA_mp .+ 2 .* SM_mp
        intArr2 = zeros( corrType, ssize[1:2] .+ 1 );

		# Nested loop establishing IV and UVW coordinates
		vfx = [ 1-mp, 0 ]
		for x2 in IAranges[2];     x1 = x2-IA_mp[2]+1; vfx .+= mp; vfy = [ 1-mp, 0 ];
			for y2 in IAranges[1]; y1 = y2-IA_mp[1]+1; vfy .+= mp;

				# 0-. Filtering undesirable IAs, ex background
				interr = view( img1, y1:y2, x1:x2 )

				if filtFun( interr ) < threshold
					SN[ vfy[1]:vfy[2], vfx[1]:vfx[2] ] .= NaN32
					continue;
				end

			    # 1.1-. Previously computed displacements shift the search area during multipass
			    u = round(Int32, U[ vfy[1], vfx[1] ]);
			    v = round(Int32, V[ vfy[1], vfx[1] ]);

				scoords, offset = SAcoords( (y1,x1,1), (y2,x2,1), (u,v,0), SM_mp, imsize );

				# 1.2-. Copying IV and SV into the padded arrays for FFT
				setTo0!( pads )
				setTo0!( padi )
			    putWithinPadded!( pads, img2, 0.0, offset, scoords );
			    putWithinPadded!( padi, img1, 0.0, (0,0,0), (y1,y2,x1,x2,1,1) );

				# 1.4-. Constructing the integral array of SV
				setTo0!( intArr2 )
				integralArraySQ!( pads, intArr2 )

				# 1.5 -. Computing sum(IA²)
				sumG2 = 0.0
				@inbounds for t2 in x1:x2
					@simd for t1 in y1:y2
						sumG2 += img1[t1,t2]*img1[t1,t2]
				end	end

				# 2-. Cross-correlation
				crossCorrelation!( NSQECC(), cmatrix, shifted, shifts, pads, padi, plan, iplan,
                                           intArr2, sumG2, ssize, IA_mp )

				# 3.1-. Calculation of displacement
		        ( r, c ) = approxTranslation( cmatrix, peak, Val(last_mp) )

				# 3.2-. Computing signal to noise ration, if needed
				if ( !ignoreSN && last_mp )
		        	SN[ vfy[1]:vfy[2], vfx[1]:vfx[2] ] = SNRatio( cmatrix, sigNoise, width )
				end

		        # 4-. Updating U, V matrices
		        U[ vfy[1]:vfy[2], vfx[1]:vfx[2] ] .= u - r;
		        V[ vfy[1]:vfy[2], vfx[1]:vfx[2] ] .= v - c;
		end	end
	end # mulitpass

    return U, V, SN
end

"""
	3D IMPLEMENTATION: f is search volume, g is interrogation volume (translated function)
"""

function crossCorrelation!( ::NSQECC, cmat::A{T,3}, shifted::A{T,3}, shifts,
                                    padF::A{C{T},3}, padG::A{C{T},3}, plan, iplan,
                                    intArr2::A{T,3}, sumG2, sizeF, sizeG
						  ) where {T<:AbstractFloat}

	# Computing IV*SV
    crossCorrelation!( FFT(), cmat, padF, padG, plan, iplan );
    Base.circshift!( shifted, cmat, shifts )

	# Computing sqrt(sum(IV²))
	denG = sqrt( sumG2 )

	# Computing sum(SV²) and sqrt(sum(SV²)).
	fh, fw, fd = sizeF
	gh, gw, gd = sizeG .- 1

    for z in 0:size(cmat,3)-1;         z1 = max( 1, fd - z ); z2 = min( fd, z1 + gd )
		for c in 0:size(cmat,2)-1;     c1 = max( 1, fw - c ); c2 = min( fw, c1 + gw )
			for r in 0:size(cmat,1)-1; r1 = max( 1, fh - r ); r2 = min( fh, r1 + gh )

				sumF2 = abs( integralArea( intArr2, (r1-1,c1-1,z1-1), (r2,c2,z2) ) )
				num   = sumF2 + sumG2 - 2*shifted[r+1,c+1,z+1]

				cmat[r+1,c+1,z+1] = 1/( 1 + num/( denG * sqrt(sumF2) ) )
	end	end	end
end

function PIV_3D( ::NSQECC, vol1::A{<:Real,3}, vol2::A{<:Real,3},
                         IA::III, SM::III, overlap::III, mpass::I, width::I,
						 peak::S, sigNoise::S, filtFun::Function, threshold::F ;
                         corrType=Float32, vfType=Float32 )

	# Calculating size of the vector field
    step   = IA .- overlap;
    VFsize = length.( StepRange.( IA, step, size(vol1) ) );

    U  = zeros( vfType, VFsize );
    V  = zeros( vfType, VFsize );
    W  = zeros( vfType, VFsize );
    SN = zeros( vfType, VFsize );
	df = zeros( vfType, VFsize );

	ignoreSN = ( sigNoise == "" );

    for mp in mpass:-1:1

        last_mp = ( mp == 1 );

        # Scaling IA, SM and step (overlap) to multi-pass iteration.
        IA_mp   =  IA  .* mp;
        SM_mp   =  SM  .* mp;
        step_mp = step .* mp;
		IAranges = StepRange.( IA_mp, step_mp, size(vol1) )

        # Initialize cross-correlation variables
        csize   = 2 .* ( IA_mp .+ SM_mp  );
        cmatrix = zeros( corrType, csize );
        shifted = zeros( corrType, csize );
        shifts  = div.(csize,2) .+ SM_mp .- 1;
        pads    = zeros( Complex{corrType}, csize );
        padi    = zeros( Complex{corrType}, csize );
        plan    =  FFTW.plan_fft!( pads );
        iplan   = FFTW.plan_ifft!( pads );

		# Initializing integral array for search area
        ssize   = IA_mp .+ 2 .* SM_mp
        intArr2 = zeros( corrType, ssize .+ 1 );

		# Nested loop establishing IV and UVW coordinates
		vfz = [ 1-mp, 0 ]
		for z2 in IAranges[3];         z1 = z2-IA_mp[3]+1; vfz .+= mp; vfx = [ 1-mp, 0 ];
			for x2 in IAranges[2];     x1 = x2-IA_mp[2]+1; vfx .+= mp; vfy = [ 1-mp, 0 ];
				for y2 in IAranges[1]; y1 = y2-IA_mp[1]+1; vfy .+= mp;

					interr = view( vol1, y1:y2, x1:x2, z1:z2 )

					if filtFun( interr ) < threshold # Filtering of undesirable IAs, ex background
						SN[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= -1.0
						continue;
					end

				    # 1-. Previously computed displacements shift the search area during multipass
				    u = round(Int32, U[ vfy[1], vfx[1], vfz[1] ]);
				    v = round(Int32, V[ vfy[1], vfx[1], vfz[1] ]);
				    w = round(Int32, W[ vfy[1], vfx[1], vfz[1] ]);

				    # 1.2-. Search volume coordinates after shifting.
					scoords, so = SAcoords( (y1,x1,z1), (y2,x2,z2), (u,v,w), SM_mp, size(vol1) );

					#println( (y1,y2,x1,x2,z1,z2), scoords, so, vfy, vfx, vfz, (u,v,w) )

					# 1.3-. Copying IV and SV into the padded arrays for FFT
					setTo0!( pads )
					setTo0!( padi )
				    putWithinPadded!( pads, vol2, 0.0, so, scoords );
				    putWithinPadded!( padi, vol1, 0.0, (0,0,0), (y1,y2,x1,x2,z1,z2) );

					# 1.4-. Constructing the integral array of SV
					setTo0!( intArr2 )
					integralArraySQ!( pads, intArr2 )

					# 1.5 -. Power of IA and SA
					sumG2 = 0.0
					@inbounds for t3 in z1:z2, t2 in x1:x2
						@simd for t1 in y1:y2
							sumG2 += vol1[t1,t2,t3]*vol1[t1,t2,t3]
					end	end

					# 2-. Cross-correlation
					crossCorrelation!( NSQECC(), cmatrix, shifted, shifts, pads, padi, plan, iplan,
                                               intArr2, sumG2, ssize, IA_mp )

				    # 3-. Extracting the displacement
				    ( r, c, z ) = approxTranslation( cmatrix, peak, Val(last_mp) )

					# 3.2 -. Computing signal to noise ration, if needed
		 			if ( !ignoreSN && last_mp )
						snres = SNRatio( cmatrix, sigNoise, width )
				    	SN[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= snres
					end

				    # 4-. Updating U, V, W matrices

				    U[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= u - r;
				    V[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= v - c;
				    W[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= w - z;
				   df[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= maxval( cmatrix );
		end	end	end
	end # multipass

    return U, V, W, SN, df
end

function PIV_3D( ::NSQECC, vol1::A{<:Real,3}, vol2::A{<:Real,3}, mask::A{Bool,1},
                         IA::III, SM::III, overlap::III, mpass::I, width::I,
						 peak::S, sigNoise::S, filtFun::Function, threshold::F ;
                         corrType=Float32, vfType=Float32 )

	# Calculating size of the vector field
    step   = IA .- overlap;
    VFsize = length.( StepRange.( IA, step, size(vol1) ) );

    U  = zeros( vfType, VFsize );
    V  = zeros( vfType, VFsize );
    W  = zeros( vfType, VFsize );
    SN = zeros( vfType, VFsize );
	df = zeros( vfType, VFsize );

	ignoreSN = ( sigNoise == "" );

    for mp in mpass:-1:1

        last_mp = ( mp == 1 );

        # Scaling IA, SM and step (overlap) to multi-pass iteration.
        IA_mp   =  IA  .* mp;
        SM_mp   =  SM  .* mp;
        step_mp = step .* mp;
		IAranges = StepRange.( IA_mp, step_mp, size(vol1) )

        # Initialize cross-correlation variables
        csize   = 2 .* ( IA_mp .+ SM_mp  );
        cmatrix = zeros( corrType, csize );
        shifted = zeros( corrType, csize );
        shifts  = div.(csize,2) .+ SM_mp .- 1;
        pads    = zeros( Complex{corrType}, csize );
        padi    = zeros( Complex{corrType}, csize );
        plan    =  FFTW.plan_fft!( pads );
        iplan   = FFTW.plan_ifft!( pads );

		# Initializing integral array for search area
        ssize   = IA_mp .+ 2 .* SM_mp
        intArr2 = zeros( corrType, ssize .+ 1 );

		# Nested loop establishing IV and UVW coordinates
		cont = 1;
		vfz = [ 1-mp, 0 ]
		for z2 in IAranges[3];         z1 = z2-IA_mp[3]+1; vfz .+= mp; vfx = [ 1-mp, 0 ];
			for x2 in IAranges[2];     x1 = x2-IA_mp[2]+1; vfx .+= mp; vfy = [ 1-mp, 0 ];
				for y2 in IAranges[1]; y1 = y2-IA_mp[1]+1; vfy .+= mp;

					if !mask[cont] # Filtering of undesirable IAs, ex background
						SN[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= -1.0
						cont += 1;
						continue;
					end

				    # 1-. Previously computed displacements shift the search area during multipass
				    u = round(Int32, U[ vfy[1], vfx[1], vfz[1] ]);
				    v = round(Int32, V[ vfy[1], vfx[1], vfz[1] ]);
				    w = round(Int32, W[ vfy[1], vfx[1], vfz[1] ]);

				    # 1.2-. Search volume coordinates after shifting.
					scoords, so = SAcoords( (y1,x1,z1), (y2,x2,z2), (u,v,w), SM_mp, size(vol1) );

					#println( (y1,y2,x1,x2,z1,z2), scoords, so, vfy, vfx, vfz, (u,v,w) )

					# 1.3-. Copying IV and SV into the padded arrays for FFT
					setTo0!( pads )
					setTo0!( padi )
				    putWithinPadded!( pads, vol2, 0.0, so, scoords );
				    putWithinPadded!( padi, vol1, 0.0, (0,0,0), (y1,y2,x1,x2,z1,z2) );

					# 1.4-. Constructing the integral array of SV
					setTo0!( intArr2 )
					integralArraySQ!( pads, intArr2 )

					# 1.5 -. Power of IA and SA
					sumG2 = 0.0
					@inbounds for t3 in z1:z2, t2 in x1:x2
						@simd for t1 in y1:y2
							sumG2 += vol1[t1,t2,t3]*vol1[t1,t2,t3]
					end	end

					# 2-. Cross-correlation
					crossCorrelation!( NSQECC(), cmatrix, shifted, shifts, pads, padi, plan, iplan,
                                               intArr2, sumG2, ssize, IA_mp )

				    # 3-. Extracting the displacement
				    ( r, c, z ) = approxTranslation( cmatrix, peak, Val(last_mp) )

					# 3.2 -. Computing signal to noise ration, if needed
		 			if ( !ignoreSN && last_mp )
						snres = SNRatio( cmatrix, sigNoise, width )
				    	SN[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= snres
					end

				    # 4-. Updating U, V, W matrices

				    U[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= u - r;
				    V[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= v - c;
				    W[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= w - z;
				   df[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= maxval( cmatrix );

				   cont += 1;
		end	end	end
	end # mulitpass
    return U, V, W, SN, df
end

"""
	out of place cross-correlation implementations
"""

function crossCorrelation( ::NSQECC, g::A{<:Real,N}, f::A{<:Real,N}; corrType=Float32 ) where {N}

	sizef = size(f)
	sizeg = size(g)

	csize   = sizef .+ sizeg;
    shifts  = div.( csize, 2 ) .+ div.( size(f) .- size(g), 2 ) .- 1;

	corr    = zeros( corrType, csize );
	shifted = zeros( corrType, csize );
	padf    = zeros( Complex{corrType}, csize );
	padg    = zeros( Complex{corrType}, csize );

    putWithinPadded!( padf, f, 0.0, (0,0,0), (1,size(f,1),1,size(f,2),1,size(f,3)) );
    putWithinPadded!( padg, g, 0.0, (0,0,0), (1,size(g,1),1,size(g,2),1,size(g,3)) );

	plan  =  plan_fft!( padf )
	iplan = plan_ifft!( padf )

	sumG2 = sum( g .* g )

	intArr2 = zeros( corrType, sizef .+ 1 )
	integralArraySQ!( padf, intArr2  )

	crossCorrelation!( NSQECC(), corr, shifted, shifts, padf, padg, plan, iplan,
                               intArr2, sumG2, sizef, sizeg )

	return corr;
end

max_issues_repo_name: ErichZimmer/quickPIV

max_forks_repo_name: ErichZimmer/quickPIV


text:
# ==============================================================================
# Xavier Gandibleux - November 2021
#   Implemented in Julia 1.6

# ==============================================================================
# Using ϵ-constraint method with vOptGeneric, compute the set of non-dominated
# points of the following problem:
#
#   Max  sum{j=1,...,n} C(k,j) x(j)                k=1,...,2
#   st   sum{j=1,...,n} A(i,j) x(j) <= B(i)        i=1,...,m
#                       x(j) = 0 or 1

# ==============================================================================
using vOptGeneric, JuMP, GLPK

# ==============================================================================
# Setting a Bi01IP instance and solving with vOptGeneric
function vSolveBi01IP(solverSelected, C, A, B)

  m, n = size(A)

  # ---- setting the model
  # println("Building...")
  Bi01IP = vModel(solverSelected)
  set_silent(Bi01IP)
  @variable(Bi01IP, x[1:n], Bin)
  @addobjective(Bi01IP, Max, sum(C[1, j] * x[j] for j = 1:n))
  @addobjective(Bi01IP, Max, sum(C[2, j] * x[j] for j = 1:n))
  @constraint(Bi01IP, cte[i = 1:m], sum(A[i, j] * x[j] for j = 1:n) <= B[i])

  # ---- Invoking the solver (epsilon constraint method)
  #println("Solving...")
  vSolve(Bi01IP, method = :epsilon, step = 0.5, verbose = false)

  # ---- Querying the results
  #println("Querying...")
  Y_N = getY_N(Bi01IP)

  return Y_N, x
end

# ==============================================================================

# Example on how to use it :

# solverSelected = GLPK.Optimizer
# Y_N, X_E = setBi01IP(solverSelected, C, A, B)

# # ---- Displaying the results (X_E and Y_N)
# for n = 1:length(Y_N)
#    X = value.(X_E, n)
#    print(findall(elt -> elt ≈ 1, X))
#    println("| z = ",Y_N[n])
# end

# for n = 1:length(Y_N)
#    println(Y_N[n])
# end


text:
using MultivariateStats
using LinearAlgebra
using Test
import Statistics: mean, cov
import Random

Random.seed!(34568)

## LinearDiscriminant

@assert LinearDiscriminant <: Discriminant

w = [1., 2., 3., 4., 5.]
b = 2.5
x = [4., 5., 2., 3., 1.]

f = LinearDiscriminant(w, b)
@test length(f) == 5
@test evaluate(f, x) == 39.5
@test evaluate(f, -x) == -34.5

@test predict(f, x) == true
@test predict(f, -x) == false

X = rand(5, 8)
Y = evaluate(f, X)
@test size(Y) == (8,)
for i = 1:8
    @test Y[i] ≈ evaluate(f, X[:,i])
end
@test predict(f, X) == (Y .> 0)

## prepare data

t1 = deg2rad(30)
t2 = deg2rad(75)

R1 = [cos(t1) -sin(t1); sin(t1) cos(t1)]
R2 = [cos(t2) -sin(t2); sin(t2) cos(t2)]

n = 20
Xp = Diagonal([1.2, 3.6]) * randn(2, n) .+ [1.0, -3.0]
Xn = Diagonal([2.8, 1.8]) * randn(2, n) .+ [-5.0, 2.0]

up = vec(mean(Xp, dims=2))
un = vec(mean(Xn, dims=2))
Cp = cov(Xp, dims=2)
Cn = cov(Xn, dims=2)
C = 0.5 * (Cp + Cn)

w_gt = C \ (up - un)
w_gt .*= (2 / (dot(w_gt, up) - dot(w_gt, un)))
b_gt = 1.0 - dot(w_gt, up)

@test dot(w_gt, up) + b_gt ≈ 1.0
@test dot(w_gt, un) + b_gt ≈ -1.0

## LDA

f = ldacov(C, up, un)
@test f.w ≈ w_gt
@test f.b ≈ b_gt

f = ldacov(Cp, Cn, up, un)
@test f.w ≈ w_gt
@test f.b ≈ b_gt

f = fit(LinearDiscriminant, Xp, Xn)
@test f.w ≈ w_gt
@test f.b ≈ b_gt

# input of Float32 type
L = fit(LinearDiscriminant, [1 2 3 4f0], [5 6 7 8f0])
@test isa(L.b, Float32)


text:
@testset "SmithWilson" begin

    @testset "InstrumentQuotes" begin

        maturities = [1.3, 2.7]
        prices = [1.1, 0.8]
        zcq = Yields.ZeroCouponQuote.(prices, maturities)
        @test isa(first(zcq), Yields.ZeroCouponQuote)

        @test_throws DimensionMismatch Yields.ZeroCouponQuote.([1.3, 2.4, 0.9], maturities)

        rates = [0.4, -0.7]
        swq = Yields.SwapQuote.(rates, maturities, 3)
        @test first(swq).frequency == 3

        @test_throws DimensionMismatch Yields.SwapQuote.([1.3, 2.4, 0.9], maturities, 3)
        @test_throws DomainError Yields.SwapQuote.(rates, maturities, 0)
        @test_throws DomainError Yields.SwapQuote.(rates, maturities, -2)

        rates = [0.4, -0.7]
        bbq = Yields.BulletBondQuote.(rates, prices, maturities, 3)
        @test first(bbq).frequency == 3
        @test first(bbq).yield == first(rates)

        @test_throws DimensionMismatch Yields.BulletBondQuote.([1.3, 2.4, 0.9], prices, maturities, 3)
        @test_throws DimensionMismatch Yields.BulletBondQuote.(rates, prices, [4.3, 5.6, 4.4, 4.4], 3)
        @test first(Yields.BulletBondQuote.(rates, [5.7], maturities, 3)).price == 5.7
        @test_throws DomainError Yields.BulletBondQuote.(rates, prices, maturities, 0)
        @test_throws DomainError Yields.BulletBondQuote.(rates, prices, maturities, -4)

    end

    @testset "SmithWilson" begin

        ufr = 0.03
        α = 0.1
        u = [5.0, 7.0]
        qb = [2.3, -1.2]

        # Basic behaviour
        sw = Yields.SmithWilson(u, qb; ufr = ufr, α = α)
        @test sw.ufr == ufr
        @test sw.α == α
        @test sw.u == u
        @test sw.qb == qb
        @test_throws DomainError Yields.SmithWilson(u, [2.4, -3.4, 8.9], ufr = ufr, α = α)

        # Empty u and Qb should result in a flat yield curve
        # Use this to test methods expected from <:AbstractYieldCurve
        # Only discount and zero are explicitly implemented, so the others should follow automatically
        sw_flat = Yields.SmithWilson(Float64[], Float64[], ufr = ufr, α = α)
        @test discount(sw_flat, 10.0) == exp(-ufr * 10.0)
        @test accumulation(sw_flat, 10.0) ≈ exp(ufr * 10.0)
        @test rate(convert(Yields.Continuous(), zero(sw_flat, 8.0))) ≈ ufr
        @test discount.(sw_flat, [5.0, 10.0]) ≈ exp.(-ufr .* [5.0, 10.0])
        @test rate(convert(Yields.Continuous(), forward(sw_flat, 5.0, 8.0))) ≈ ufr

        # A trivial Qb vector (=0) should result in a flat yield curve
        ufr_curve = Yields.SmithWilson(u, [0.0, 0.0], ufr = ufr, α = α)
        @test discount(ufr_curve, 10.0) == exp(-ufr * 10.0)

        # A single payment at time 4, zero interest
        curve_with_zero_yield = Yields.SmithWilson([4.0], reshape([1.0], 1, 1), [1.0], ufr = ufr, α = α)
        @test discount(curve_with_zero_yield, 4.0) == 1.0

        # In the long end it's still just UFR
        @test rate(convert(Yields.Continuous(), forward(curve_with_zero_yield, 1000.0, 2000.0))) ≈ ufr

        # Three maturities have known discount factors
        times = [1.0, 2.5, 5.6]
        prices = [0.9, 0.7, 0.5]
        cfs = [1 0 0
            0 1 0
            0 0 1]

        curve_three = Yields.SmithWilson(times, cfs, prices, ufr = ufr, α = α)
        @test transpose(cfs) * discount.(curve_three, times) ≈ prices

        # Two cash flows with payments at three times
        prices = [1.0, 0.9]
        cfs = [0.1 0.1
            1.0 0.1
            0.0 1.0]
        curve_nondiag = Yields.SmithWilson(times, cfs, prices, ufr = ufr, α = α)
        @test transpose(cfs) * discount.(curve_nondiag, times) ≈ prices

        # Round-trip zero coupon quotes
        zcq_times = [1.2, 4.5, 5.6]
        zcq_prices = [1.0, 0.9, 1.2]
        zcq = Yields.ZeroCouponQuote.(zcq_prices, zcq_times)
        sw_zcq = Yields.SmithWilson(zcq, ufr = ufr, α = α)
        @testset "ZeroCouponQuotes round-trip" for idx = 1:length(zcq_times)
            @test discount(sw_zcq, zcq_times[idx]) ≈ zcq_prices[idx]
        end

        # uneven frequencies
        swq_maturities = [1.2, 2.5, 3.6]
        swq_interests = [-0.02, 0.3, 0.04]
        frequency = [2, 1, 2]
        swq = Yields.SwapQuote.(swq_interests, swq_maturities, frequency)
        swq_times = 0.5:0.5:3.5   # Maturities are rounded down to multiples of 1/frequency, [1.0, 2.5, 3.5]
        swq_payments = [
            -0.01 0.3 0.02
            0.99 0 0.02
            0 0.3 0.02
            0 0 0.02
            0 1.3 0.02
            0 0 0.02
            0 0 1.02
        ]

        # Round-trip swap quotes
        swq_maturities = [1.2, 2.5, 3.6]
        swq_interests = [-0.02, 0.3, 0.04]
        frequency = 2
        swq = Yields.SwapQuote.(swq_interests, swq_maturities, frequency)
        swq_times = 0.5:0.5:3.5   # Maturities are rounded down to multiples of 1/frequency, [1.0, 2.5, 3.5]
        swq_payments = [-0.01 0.15 0.02
            0.99 0.15 0.02
            0.0 0.15 0.02
            0.0 0.15 0.02
            0.0 1.15 0.02
            0.0 0.0 0.02
            0.0 0.0 1.02]
        sw_swq = Yields.SmithWilson(swq, ufr = ufr, α = α)
        @testset "SwapQuotes round-trip" for swapIdx = 1:length(swq_interests)
            @test sum(discount.(sw_swq, swq_times) .* swq_payments[:, swapIdx]) ≈ 1.0
        end

        # Round-trip bullet bond quotes (reuse data from swap quotes)
        bbq_prices = [1.3, 0.1, 4.5]
        bbq = Yields.BulletBondQuote.(swq_interests, bbq_prices, swq_maturities, frequency)
        sw_bbq = Yields.SmithWilson(bbq, ufr = ufr, α = α)
        @testset "BulletBondQuotes round-trip" for bondIdx = 1:length(swq_interests)
            @test sum(discount.(sw_bbq, swq_times) .* swq_payments[:, bondIdx]) ≈ bbq_prices[bondIdx]
        end

        @testset "SW ForwardStarting" begin
            fwd_time = 1.0
            fwd = Yields.ForwardStarting(sw_swq, fwd_time)

            @test discount(fwd, 3.7) ≈ discount(sw_swq, fwd_time, fwd_time + 3.7)
        end

        # EIOPA risk free rate (no VA), 31 August 2021.
        # https://www.eiopa.europa.eu/sites/default/files/risk_free_interest_rate/eiopa_rfr_20210831.zip
        eiopa_output_qb = [-0.59556534586390800
            -0.07442224713453920
            -0.34193181987682400
            1.54054875814153000
            -2.15552046042343000
            0.73559290752221900
            1.89365225129089000
            -2.75927773116240000
            2.24893737130629000
            -1.51625404117395000
            0.19284859623817400
            1.13410725406271000
            0.00153268224642171
            0.00147942301778158
            -1.85022125156483000
            0.00336230229850928
            0.00324546553910162
            0.00313268874430658
            0.00302383083427276
            1.36047951448615000]
        eiopa_output_u = 1:20
        eiopa_ufr = log(1.036)
        eiopa_α = 0.133394
        sw_eiopa_expected = Yields.SmithWilson(eiopa_output_u, eiopa_output_qb; ufr = eiopa_ufr, α = eiopa_α)

        eiopa_eurswap_maturities = [1:12; 15; 20]
        eiopa_eurswap_rates = [-0.00615, -0.00575, -0.00535, -0.00485, -0.00425, -0.00375, -0.003145,
            -0.00245, -0.00185, -0.00125, -0.000711, -0.00019, 0.00111, 0.00215]   # Reverse engineered from output curve. This is the full precision of market quotes.
        eiopa_eurswap_quotes = Yields.SwapQuote.(eiopa_eurswap_rates, eiopa_eurswap_maturities, 1)
        sw_eiopa_actual = Yields.SmithWilson(eiopa_eurswap_quotes, ufr = eiopa_ufr, α = eiopa_α)

        @testset "Match EIOPA calculation" begin
            @test sw_eiopa_expected.u ≈ sw_eiopa_actual.u
            @test sw_eiopa_expected.qb ≈ sw_eiopa_actual.qb
        end
    end
end


text:
# Move the swarm

function move!(rng::AbstractRNG, state::ParticleSwarmState{T}, w, c1, c2) where {T}
    X, V = state.X, state.V
    @. V = w*V + c1*rand(rng, T)*(state.pbest_X - X) + c2*rand(rng, T)*(state.gbest_X - X)
    X .+= V
    for (r, idx) in zip(state.ranges, state.indices)
        constrain!(r, view(X, :, idx))
    end
    return state
end

# Constrain particles' positions

function constrain!(r::NominalRange, X::AbstractArray{T}) where {T}
    @. X = min(one(T), max(eps(T), X))
    return X ./= sum(X, dims=2)
end

function constrain!(r::NumericRange, X::AbstractArray{T}) where {T}
    return @. X = min(T(r.upper), max(T(r.lower), X))
end

# Update pbest

function pbest!(state::ParticleSwarmState, measurements, prob_shift)
    X, pbest, pbest_X = state.X, state.pbest, state.pbest_X
    improved = measurements .<= pbest
    pbest[improved] .= measurements[improved]
    for (r, p, i) in zip(state.ranges, state.parameters, state.indices)
        @views _pbest!(pbest_X[improved, i], X[improved, i], r, p[improved], prob_shift)
    end
    return state
end

_pbest!(pbest_X, X, r::ParamRange, p, prob_shift) = pbest_X .= X

function _pbest!(pbest_X::AbstractArray{T}, X, r::NominalRange, p, prob_shift) where {T}
    pbest_X .= X .* (one(T) - T(prob_shift))
    sampled = map(pᵢ -> findfirst(==(pᵢ), r.values), p)
    for (i, j) in zip(axes(pbest_X, 1), sampled)
        pbest_X[i, j] += T(prob_shift)
    end
    return pbest_X
end

# Update gbest

function gbest!(state::ParticleSwarmState)
    pbest, pbest_X, gbest, gbest_X = state.pbest, state.pbest_X, state.gbest, state.gbest_X
    best, i = findmin(pbest)
    gbest .= best
    gbest_X .= pbest_X[i, :]'
    return state
end


text:
r = @theory begin
    foo(x,y) => 2*x%y
    foo(x,y) => sin(x)
    sin(x) => foo(x,x)
end
@testset "Basic Equalities 1" begin
    @test (@areequal r foo(b,c) foo(d,d)) == false
end

r = @theory begin
    a * 1 => foo
    a * 2 => bar
    1 * a => baz
    2 * a => mag
end
@testset "Matching Literals" begin
    @test (@areequal r a * 1 foo) == true
    @test (@areequal r a * 2 foo) == false
    @test (@areequal r a * 1 bar) == false
    @test (@areequal r a * 2 bar) == true

    @test (@areequal r 1 * a baz) == true
    @test (@areequal r 2 * a baz) == false
    @test (@areequal r 1 * a mag) == false
    @test (@areequal r 2 * a mag) == true
end

comm_monoid = @commutative_monoid (*) 1
@testset "Basic Equalities - Commutative Monoid" begin
    @test true == (@areequal comm_monoid a*(c*(1*d)) c*(1*(d*a)) )
    @test true == (@areequal comm_monoid x*y y*x )
    @test true == (@areequal comm_monoid (x*x)*(x * 1) x*(x*x) )
end

comm_group = @abelian_group (+) 0 inv
t = comm_monoid ∪ comm_group ∪ distrib(:(*), :(+))
@testset "Basic Equalities - Comm. Monoid, Abelian Group, Distributivity" begin
    @test true == (@areequal t (a * b) + (a * c) a*(b+c) )
    @test true == (@areequal t a*(c*(1*d)) c*(1*(d*a)) )
    @test true == (@areequal t a+(b*(c*d)) ((d*c)*b)+a )
    @test true == (@areequal t (x+y)*(a+b) ((a*(x+y)) + b*(x+y)) ((x*(a+b)) + y*(a+b)) )
    @test true == (@areequal t (((x*a + x*b) + y*a) + y*b) (x+y)*(a+b) )
    @test true == (@areequal t a+(b*(c*d)) ((d*c)*b)+a )
    @test true == (@areequal t a+inv(a) 0 (x*y)+inv(x*y) 1*0 )
end

@testset "Basic Equalities - False statements" begin
    @test false == (@areequal t (a * b) + (a * c) a*(b+a))
    @test false == (@areequal t (a * c) + (a * c) a*(b+c) )
    @test false == (@areequal t a*(c*c) c*(1*(d*a)) )
    @test false == (@areequal t c+(b*(c*d)) ((d*c)*b)+a )
    @test false == (@areequal t (x+y)*(a+c) ((a*(x+y)) + b*(x+y)) )
    @test false == (@areequal t ((x*(a+b)) + y*(a+b)) (x+y)*(a+c) )
    @test false == (@areequal t (((x*a + x*b) + y*a) + y*b) (x+y)*(a+x) )
    @test false == (@areequal t a+(b*(c*a)) ((d*c)*b)+a )
    @test false == (@areequal t a+inv(a) a )
    @test false == (@areequal t (x*y)+inv(x*y) 1 )
end

# expr = cleanast(:(1 * 1 * 1 * 1 * 1 * zoo * 1 * 1 * foo * 1))
#
# G = EGraph(expr)
#
# @time saturate!(G, comm_monoid)
#
# G.H |> display