text:
# Tests for the graph module
using Arbiter.ArbiterGraphs

facts("push!") do
    graph = ArbiterGraph()

    @fact nodes(graph) --> ImmutableNodeSet()
    @fact roots(graph) --> ImmutableNodeSet()
    @fact :foo in graph --> false

    push!(graph, :foo)

    @fact nodes(graph) --> ImmutableNodeSet((:foo,))
    @fact roots(graph) --> ImmutableNodeSet((:foo,))
    @fact :foo in graph --> true
    @fact children(graph, :foo) --> ImmutableNodeSet()
    @fact parents(graph, :foo) --> ImmutableNodeSet()
    @fact ancestor_of(graph, :foo, :foo) --> false

    push!(graph, :bar, (:foo, :baz))

    @fact nodes(graph) --> ImmutableNodeSet((:foo, :bar, :baz))
    @fact roots(graph) --> ImmutableNodeSet((:foo,))

    @fact :foo in graph --> true
    @fact :bar in graph --> true
    @fact :baz in graph --> true

    @fact children(graph, :foo) --> ImmutableNodeSet((:bar,))
    @fact children(graph, :bar) --> ImmutableNodeSet()
    @fact children(graph, :baz) --> ImmutableNodeSet((:bar,))

    @fact parents(graph, :foo) --> ImmutableNodeSet()
    @fact parents(graph, :bar) --> ImmutableNodeSet((:foo, :baz))
    @fact parents(graph, :baz) --> ImmutableNodeSet()

    @fact ancestor_of(graph, :foo, :foo) --> false
    @fact ancestor_of(graph, :foo, :bar) --> false
    @fact ancestor_of(graph, :foo, :baz) --> false

    @fact ancestor_of(graph, :bar, :foo) --> true
    @fact ancestor_of(graph, :bar, :bar) --> false
    @fact ancestor_of(graph, :bar, :baz) --> true

    @fact ancestor_of(graph, :baz, :foo) --> false
    @fact ancestor_of(graph, :baz, :bar) --> false
    @fact ancestor_of(graph, :baz, :baz) --> false

    @fact_throws ErrorException push!(graph, :baz, (:bar,))
    @fact_throws ErrorException push!(graph, :ouroboros, (:ouroboros,))
    @fact_throws ErrorException push!(graph, :foo)

    @fact nodes(graph) --> ImmutableNodeSet((:foo, :bar, :baz))
    @fact roots(graph) --> ImmutableNodeSet((:foo,))
end

facts("pop! orphan") do
    graph = ArbiterGraph()

    push!(graph, :node)
    push!(graph, :bar, (:foo,))
    push!(graph, :baz, (:bar,))
    push!(graph, :beta, (:alpha,))
    push!(graph, :bravo, (:alpha,))

    @fact nodes(graph) --> ImmutableNodeSet((:node, :foo, :bar, :baz, :alpha, :beta, :bravo))
    @fact roots(graph) --> ImmutableNodeSet((:node,))

    # node with no children/parents
    @fact pop!(graph, :node; strategy=Orphan) --> NodeSet((:node,))
    @fact nodes(graph) --> ImmutableNodeSet((:foo, :bar, :baz, :alpha, :beta, :bravo))

    # node with child, unique stub parent
    @fact pop!(graph, :bar; strategy=Orphan) --> NodeSet((:bar, :foo))
    @fact nodes(graph) --> ImmutableNodeSet((:baz, :alpha, :beta, :bravo))
    @fact roots(graph) --> ImmutableNodeSet((:baz,))

    # node with non-unique stub parent
    @fact pop!(graph, :bravo; strategy=Orphan) --> NodeSet((:bravo,))
    @fact nodes(graph) --> ImmutableNodeSet((:baz, :alpha, :beta))
    @fact roots(graph) --> ImmutableNodeSet((:baz,))

    # stub
    @fact pop!(graph, :alpha; strategy=Orphan) --> NodeSet((:alpha,))
    @fact nodes(graph) --> ImmutableNodeSet((:baz, :beta))
    @fact roots(graph) --> ImmutableNodeSet((:baz, :beta))

    @fact_throws KeyError pop!(graph, :fake; strategy=Orphan)
end

facts("pop! promote") do
    # promote is the default
    graph = ArbiterGraph()

    push!(graph, :aye)
    push!(graph, :insect)
    push!(graph, :bee, (:aye, :insect))
    push!(graph, :cee, (:bee,))
    push!(graph, :child, (:stub, :stub2))
    push!(graph, :grandchild, (:child,))

    @fact nodes(graph) --> (
        ImmutableNodeSet((:aye, :insect, :bee, :cee, :child, :stub, :stub2, :grandchild))
    )
    @fact roots(graph) --> ImmutableNodeSet((:aye, :insect))

    # two new parents
    @fact pop!(graph, :bee) --> NodeSet((:bee,))
    @fact nodes(graph) --> (
        ImmutableNodeSet((:aye, :insect, :cee, :child, :stub, :stub2, :grandchild))
    )
    @fact roots(graph) --> ImmutableNodeSet((:aye, :insect))
    @fact children(graph, :aye) --> ImmutableNodeSet((:cee,))
    @fact children(graph, :insect) --> ImmutableNodeSet((:cee,))
    @fact parents(graph, :cee) --> ImmutableNodeSet((:aye, :insect))

    # now with stubs
    @fact pop!(graph, :child) --> NodeSet((:child,))
    @fact nodes(graph) --> ImmutableNodeSet((:aye, :insect, :cee, :stub, :stub2, :grandchild))
    @fact roots(graph) --> ImmutableNodeSet((:aye, :insect))
    @fact children(graph, :stub) --> ImmutableNodeSet((:grandchild,))
    @fact children(graph, :stub2) --> ImmutableNodeSet((:grandchild,))
    @fact parents(graph, :grandchild) --> ImmutableNodeSet((:stub, :stub2))

    # delete a stub
    @fact pop!(graph, :stub2) --> NodeSet((:stub2,))
    @fact nodes(graph) --> ImmutableNodeSet((:aye, :insect, :cee, :stub, :grandchild))
    @fact roots(graph) --> ImmutableNodeSet((:aye, :insect))
    @fact children(graph, :stub) --> ImmutableNodeSet((:grandchild,))
    @fact parents(graph, :grandchild) --> ImmutableNodeSet((:stub,))
end

facts("pop! remove") do
    # remove a node (and its children) from a graph
    graph = ArbiterGraph()

    push!(graph, :node)
    push!(graph, :bar, (:foo,))
    push!(graph, :baz, (:bar,))
    push!(graph, :beta, (:alpha,))
    push!(graph, :bravo, (:alpha,))

    @fact nodes(graph) --> ImmutableNodeSet((:node, :foo, :bar, :baz, :alpha, :beta, :bravo))
    @fact roots(graph) --> ImmutableNodeSet((:node,))

    # node with no children/parents
    @fact pop!(graph, :node; strategy=Remove) --> NodeSet((:node,))
    @fact nodes(graph) --> ImmutableNodeSet((:foo, :bar, :baz, :alpha, :beta, :bravo))
    @fact roots(graph) --> ImmutableNodeSet()

    # node with child, unique stub parent
    @fact pop!(graph, :bar; strategy=Remove) --> NodeSet((:bar, :foo, :baz))
    @fact nodes(graph) --> ImmutableNodeSet((:alpha, :beta, :bravo))
    @fact roots(graph) --> ImmutableNodeSet()

    # node with non-unique stub parent
    @fact pop!(graph, :bravo; strategy=Remove) --> NodeSet((:bravo,))
    @fact nodes(graph) --> ImmutableNodeSet((:alpha, :beta))
    @fact roots(graph) --> ImmutableNodeSet()

    # stub
    @fact pop!(graph, :alpha; strategy=Remove) --> NodeSet((:alpha, :beta))
    @fact nodes(graph) --> ImmutableNodeSet()

    @fact_throws KeyError pop!(graph, :fake; strategy=Remove)
end

facts("prune!") do
    # prune a graph
    graph = ArbiterGraph()

    push!(graph, :node)
    push!(graph, :bar, (:foo,))
    push!(graph, :baz, (:bar,))
    push!(graph, :beta, (:alpha,))
    push!(graph, :bravo, (:alpha,))

    @fact nodes(graph) --> ImmutableNodeSet((:node, :foo, :bar, :baz, :alpha, :beta, :bravo))
    @fact roots(graph) --> ImmutableNodeSet((:node,))

    @fact prune!(graph) --> NodeSet((:bar, :baz, :beta, :bravo))
    @fact nodes(graph) --> ImmutableNodeSet((:node,))
    @fact roots(graph) --> ImmutableNodeSet((:node,))

    @fact prune!(graph) --> NodeSet()
    @fact nodes(graph) --> ImmutableNodeSet((:node,))
    @fact roots(graph) --> ImmutableNodeSet((:node,))
end

facts("equality") do
    graph = ArbiterGraph()

    @fact graph == 1 --> false
    @fact graph != 0 --> true

    other = ArbiterGraph()

    @fact graph --> other
    @fact graph != other --> false

    push!(graph, :foo)

    @fact graph --> not(other)
    @fact graph != other --> true

    push!(graph, :bar, (:foo,))
    push!(other, :bar, (:foo,))

    # still shouldn't match graph['foo'] is a stub
    @fact graph --> not(other)
    @fact graph != other --> true

    push!(other, :foo)

    @fact graph --> other
    @fact graph != other --> false
end


text:
"""
    MaxPool(k::NTuple; stride = k, padding = 0)
This is a layer for max pooling with kernel size `k`.

Expects as input an array with `ndims(x) == N+2`, i.e. channel and batch dimensions, after the `N` feature dimensions, where `N = length(out)`.

The default stride is the same as kernel size `k`.

# Example
```
julia> N = NetWork(Conv((2, 2), 5=>2, relu), MaxPool(2, 2))
Layer1 : Convolution(k:(2, 2), IO:5=>2, σ:relu)
Layer2 : MaxPool(k:(2, 2), stride:(2, 2) padding:(0, 0, 0, 0))

julia> x = rand(Float64, 10, 10, 5, 5) |> size
(10, 10, 5, 5)

julia> N(x) |> size
(4, 4, 5, 2)
```
"""
mutable struct MaxPool{N, M} <: NParam
    k::NTuple{N, Int}
    stride::NTuple{N, Int}
    padding::NTuple{M, Int}
end

function MaxPool(k::NTuple{N, Int}; stride = k, padding = 0) where N
    stride = expand(Val(N), stride)
    padding = expand_padding(padding, k, 1, stride)
    return MaxPool(k, stride, padding)
end

function Base.show(io::IO, MP::MaxPool)
    print("MaxPool(k:"*string(MP.k)*", stride:"*string(MP.stride)*", padding:"*string(MP.padding)*")")
end

function (mp::MaxPool)(x)
    pooldims = PoolDims(x, mp.k; padding=mp.padding, stride=mp.stride)
    return maxpool(x, pooldims)
end

ext: jl


text:
function determine_grid_size(p::Array{S,2}) where {S <: Integer}

  q = similar(p)

  for i in eachindex(p)
    if p[i] > 0
      q[i] = 2
    else
      q[i] = 1
    end
  end

  s = 0
  for i in axes(q,1)
    t = 1
    for j in axes(q,2)
      t *= q[i,j]
    end
    s += t
  end

  return (s, size(q,2))

end


text:
function lex_m!(g::AbstractGraph)
    origin_edges = edges(g)
    labels = zeros(Int, nv(g))
    order = zeros(Int, nv(g))
    unumbered = collect(vertices(g))
    for i in nv(g):-1:1
        old_labels = copy(labels)
        v = argmax(labels)
        order[i] = v
        filter!(x->x != v, unumbered)
        labels[v] = -1
        for u in unumbered
            lower_labels = filter(x -> old_labels[x] < old_labels[u], unumbered)
            sg, vmap = induced_subgraph(g, vcat([v, u], lower_labels))
            if has_path(sg, findall(==(v), vmap)[1], findall(==(u), vmap)[1])
                labels[u] += i
                add_edge!(g, v, u)
            end
        end
    end
    data = Dict("order" => order,
                "added_edges" => setdiff(edges(g), origin_edges))
    return data
end

function lex_m(g::AbstractGraph; kwargs...)
    cg = SimpleGraph(g)
    data = LEX_M!(cg)
    return cg, data
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using LinearAlgebra

# const slow_mrrk_methods =
#     ((LSRK54CarpenterKennedy, 4), (LSRK144NiegemannDiehlBusch, 4))
# const fast_mrrk_methods = (
#     (LSRK54CarpenterKennedy, 4),
#     (LSRK144NiegemannDiehlBusch, 4),
#     (SSPRK33ShuOsher, 3),
#     (SSPRK34SpiteriRuuth, 3),
# )
const explicit_methods = (
    (LSRK54CarpenterKennedy(), 4),
    (LSRK144NiegemannDiehlBusch(), 4),
    (LSRKEulerMethod(), 1),
    (SSPRK22Heuns(), 2),
    (SSPRK22Ralstons(), 2),
    (SSPRK33ShuOsher(), 3),
    (SSPRK34SpiteriRuuth(), 3),
)

 const imex_methods = (
#     (ARK1ForwardBackwardEuler, 1),
#     (ARK2ImplicitExplicitMidpoint, 2),
     (ARK2GiraldoKellyConstantinescu, 2),
#     (ARK437L2SA1KennedyCarpenter, 4),
#     (ARK548L2SA2KennedyCarpenter, 5),
 )

# const mis_methods =
#     ((MIS2, 2), (MIS3C, 2), (MIS4, 3), (MIS4a, 3), (TVDMISA, 2), (TVDMISB, 2))

# const mrigark_erk_methods = ((MRIGARKERK33aSandu, 3), (MRIGARKERK45aSandu, 4))

# const mrigark_irk_methods = (
#     (MRIGARKESDIRK24LSA, 2),
#     (MRIGARKESDIRK23LSA, 2),
#     (MRIGARKIRK21aSandu, 2),
#     (MRIGARKESDIRK34aSandu, 3),
#     (MRIGARKESDIRK46aSandu, 4),
# )

# const fast_mrigark_methods =
#     ((LSRK54CarpenterKennedy, 4), (LSRK144NiegemannDiehlBusch, 4))

struct DirectSolver end

DirectSolver(args...) = DirectSolver()

function (::DirectSolver)(x,A,b,matrix_updated; kwargs...)
    n = length(x)
    M = mapslices(y -> mul!(similar(y), A, y), Matrix{eltype(x)}(I,n,n), dims=1)
    x .= M \ b
end


text:
@show findfirst(["no", "?", "yes", "maybe", "yes"], "yes")
@show indexin(["yes"], ["no", "?", "yes", "maybe", "yes"])
@show findin(["no", "?", "yes", "maybe", "yes"], ["yes"])
@show find(["no", "?", "yes", "maybe", "yes"] .== "yes")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    n_man(sem_fit::SemFit)
    n_man(model::AbstractSemSingle)

Return the number of manifest variables.
"""
function n_man end

n_man(sem_fit::SemFit) = n_man(sem_fit.model)

n_man(model::AbstractSemSingle) = n_man(model.observed)


text:
include(joinpath(pwd(), "examples/contact_implicit/hopper_vertical.jl"))

using TrajectoryOptimization
using Altro
using RobotDynamics
using StaticArrays
using LinearAlgebra
using ForwardDiff
const TO = TrajectoryOptimization

n, m = size(model)

# Horizon
T = 101

# Time step
tf = 1.0
h = tf / (T - 1)

# Objective
include("nonlinear_objective.jl")
z_h = 0.0
q1 = [0.0, 0.5 + z_h, 0.0, 0.25]
x1 = [q1; q1; 0.0]
qT = [0.0, 0.5 + z_h, 0.0, 0.5]
xT = [qT; qT; 0.0]
X0 = linear_interpolation(x1, xT, T)

include("initial_torque.jl")
u1 = initial_torque(model, q1, h)
U0 = [copy(u1) + 0.0 * rand(m) for t = 1:T-1]
u_ref = u1 #[u1[1:model.nu]; zeros(m - model.nu)]

function obj_stage_1(x, u)
    J = 0.0
    J += (u - u_ref)' * Diagonal([1.0e-1 * ones(model.nu); 1.0e-1 * ones(m - model.nu)]) * (u - u_ref)
    # J += 10.0 * (kinematics(model, view(x, model.nq .+ (1:model.nq)))[2] - 0.5)^2
    # J += 1.0 * sum((view(x, model.nq .+ (1:model.nq)) - view(x, 1:model.nq)).^2) / h
    return J
end

function obj_stage_t(x, u)
    J = 0.0
    J += (u - u_ref)' * Diagonal([1.0e-1 * ones(model.nu); 1.0e-1 * ones(m - model.nu)]) * (u - u_ref)
    J += 100.0 * (kinematics(model, view(x, model.nq .+ (1:model.nq)))[2] - 0.0)^2
    # J += 1.0 * sum((view(x, model.nq .+ (1:model.nq)) - view(x, 1:model.nq)).^2) / h
    # J += 10.0 * (view(x, model.nq .+ (1:model.nq))[2] - 1.0)^2
    return J
end

function obj_T(x)
    (x - xT)' * (x - xT) * 0.0
end

nl_stage_1 = NonlinearCostFunction{n, m}(obj_stage_1)
nl_stage_t = NonlinearCostFunction{n, m}(obj_stage_t)
nl_T = NonlinearCostFunction{n, m}(obj_T)
obj = Objective([t != 1 ? nl_stage_t : nl_stage_1 for t = 1:T-1], nl_terminal)

# Constraints
include(joinpath(pwd(), "contact_constraints.jl"))
cons = ConstraintList(n, m, T)

add_constraint!(cons, GoalConstraint(xT, (model.nq .+ (1:model.nq))), T)
# add_constraint!(cons, GoalConstraint(xT, (1:2*model.nq)), T)

uu = Inf * ones(m)
ul = zeros(m)
ul[model.idx_u] .= -Inf
add_constraint!(cons, BoundConstraint(n, m,
    x_min = [model.qL; model.qL; 0.0], x_max = [model.qU; model.qU; Inf],
    u_min = ul, u_max = uu), 1:T-1)

con_sd = SD(n, model.nc, model)
con_ic = IC(n, model.nc, model)
con_fc = FC(m, model.nc, model)
con_ns = NS(n, model.nc, model, h)
add_constraint!(cons, con_sd, 1:T)
add_constraint!(cons, con_ic, 2:T)
add_constraint!(cons, con_fc, 1:T-1)
add_constraint!(cons, con_ns, 2:T)

# x0 = @SVector ones(n)
# u0 = @SVector ones(m)
# TO.evaluate(con_sd, x0)
# TO.evaluate(con_ic, x0)
# TO.evaluate(con_fc, u0)
# TO.evaluate(con_fc, x0)

# Create and solve problem
prob = Problem(model, obj, xT, tf,
    U0 = U0,
    x0 = x1, constraints = cons, integration = PassThrough)
rollout!(prob)

solver = ALTROSolver(prob, projected_newton = false,
    constraint_tolerance = 1.0e-2, verbose = 2,
    penalty_initial = 10.0, penalty_scaling = 10.0)
cost(solver)           # initial cost
@time solve!(solver)         # solve with ALTRO
max_violation(solver)  # max constraint violation
cost(solver)           # final cost
iterations(solver)     # total number of iterations

# Get the state and control trajectories
X = states(solver)
U = controls(solver)

# using Plots
# plot(hcat(X...)')
# plot(hcat(U...)', linetype=:steppost)

# vis = Visualizer()
# render(vis)
visualize!(vis, model, state_to_configuration(X), Δt = h)


text:
using jInv.Mesh
using jInv.LinearSolvers
using Test
using KrylovMethods
using LinearAlgebra

println("===  Example 2D DivSigGrad ====");

domain = [0.0, 1.0, 0.0, 1.0];
n      = [5,5];
Mr     = getRegularMesh(domain,n)
G      = getNodalGradientMatrix(Mr);
m      = sparse(Diagonal(exp.(randn(size(G,1)))));
Ar     = G'*m*G;
Ar     = Ar + 1e-1*norm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));
N      = size(Ar,2);
b      = Ar*rand(N);

sPCG   = getIterativeSolver(KrylovMethods.cg,out=1,sym=1);

x,  = solveLinearSystem(Ar,b,sPCG);
@test norm(Ar*x-b)/norm(b) < sPCG.tol
IterMethod = (A,b;M=M,tol=1e-1,maxIter=10,out=-1) ->
                  bicgstb(A,b,M1=M,tol=tol,maxIter=maxIter,out=out)
sBiCG = getIterativeSolver(IterMethod,out=1);
x, = solveLinearSystem(Ar,b,sBiCG);
@test norm(Ar*x-b)/norm(b) < sBiCG.tol

println("===  Test for nonsymmetric matrices ====");
n = 100
m = 4
sBiCG.PC = :jac
sBiCG.doClear=true
sBiCG.out=0
A = sprandn(n,n,5.0/n) + 10*sparse(1.0I,n,n)
B = randn(n,m)
X, = solveLinearSystem(A,B,sBiCG)
for k=1:m
	@test norm(A*X[:,k] - B[:,k])/norm(B[:,k]) < sBiCG.tol
end

sBiCG.isTranspose=true
Xt, = solveLinearSystem(sparse(A'),B,sBiCG)
for k=1:m
	@test norm(A*X[:,k] - B[:,k])/norm(B[:,k]) < sBiCG.tol
	@test norm(Xt[:,k] - X[:,k])/norm(X[:,k]) < 1e-13
end


text:
export Arguments
"""
    type Arguments
        positional::Tuple
        keyword::Dict{Symbol, Any}
    end

Will store positional and keyword arguments for later use.

Create with [`collect_arguments`](@ref). You can also [`merge`](@ref) two
`Arguments`, [`push`](@ref) or [`unshift`](@ref) in new
arguments, and run with [`run`](@ref).
"""
type Arguments
  positional::Tuple
  keyword::Dict{Symbol, Any}
end

export LazyCall
"""
    type LazyCall{T <: Function}
        arguments::Arguments
        function_call::T
    end

Will store a function along with its arguments for later use. Create
with [`collect_call`](@ref) and run with [`run`](@ref)
"""
type LazyCall{T <: Function}
    arguments::Arguments
    function_call::T
end

"""
    merge(a::Arguments, b::Arguments)

Merge two [`Arguments`](@ref) types.

Positional arguments are added at the end, and new keyword arguments
are added to old keyword arguments, or, if the keys match, overwrite
them.

# Examples
```julia
merge_test = @chain begin
    collect_arguments(1, a = 2, b = 3)
    merge(_, collect_arguments(4, a = 5, c = 6) )
end

@test merge_test == collect_arguments(1, 4, a = 5, b = 3, c = 6)
```
"""
function Base.merge(a::Arguments, b::Arguments)
    positional = (a.positional..., b.positional...)
    keyword = merge(a.keyword, b.keyword)
    Arguments(positional, keyword)
end

"""
    merge(lazy_call::LazyCall, arguments::Arguments)

`merge` `arguments` into the `arguments` of `lazy_call`.

# Examples
```julia
merge_test = @chain begin
    collect_arguments([1, 2])
    unshift(_, vcat)
    LazyCall(_, map)
    merge(_, collect_arguments([3, 4]) )
    run
end

@test merge_test == [[1, 3], [2, 4]]
```
"""
Base.merge(lazy_call::LazyCall, arguments::Arguments) = @chain begin
    lazy_call.arguments
    merge(_, arguments)
    LazyCall(_, lazy_call.function_call)
end

export push
"""
    push(arguments::Arguments, positional...; keyword...)

Add positional and keyword arguments to `arguments`.

Positional arguments are added at the end, and new keyword arguments
are added to old keyword arguments, or, if the keys match, overwrite
them.

# Examples
```julia
push_test = @chain begin
    collect_arguments(1, a = 2, b = 3)
    push(_, 4, a = 5, c = 6)
end

@test push_test == collect_arguments(1, 4, a = 5, b = 3, c = 6)
```
"""
push(a::Arguments, positional...; keyword...) =
  @chain begin
      keyword
      Dict
      Arguments(positional, _)
      merge(a, _)
  end

"""
    push(lazy_call::LazyCall, positional...; keyword...)

`push` to `lazy_call.arguments`.

# Examples
```julia
push_test = @chain begin
    collect_arguments([1, 2])
    unshift(_, vcat)
    LazyCall(_, map)
    push(_, [3, 4])
    run
end

@test push_test == [[1, 3], [2, 4]]
```
"""
push(lazy_call::LazyCall, positional...; keyword...) = @chain begin
    lazy_call.arguments
    push(_, positional...; keyword...)
    LazyCall(_, lazy_call.function_call)
end

export unshift
"""
    unshift(arguments::Arguments, positional...)

Add positional arguments to `arguments`.

New arguments are added at the start.

# Examples
```julia
unshift_test = @chain begin
    collect_arguments(2, a = 3)
    unshift(_, 1)
end

@test unshift_test == collect_arguments(1, 2, a = 3)
```
"""
unshift(a::Arguments, positional...) = @chain begin
    positional
    Arguments(_, Dict() )
    merge(_, a)
end

"""
    unshift(lazy_call::LazyCall, positional...)

`unshift` to `lazy_call.arguments`.

# Examples
```julia
unshift_test = @chain begin
    collect_arguments([1, 2], [3, 4])
    LazyCall(_, map)
    unshift(_, vcat)
    run
end

@test unshift_test == [[1, 3], [2, 4]]
```
"""
unshift(lazy_call::LazyCall, positional...) = @chain begin
    lazy_call.arguments
    unshift(_, positional...)
    LazyCall(_, lazy_call.function_call)
end

export collect_arguments
"""
    collect_arguments(positional...; keyword...)

Easy way to build an [`Arguments`](@ref) type.

# Examples
```julia
a = collect_arguments(1, 2, a = 3, b = 4)
@test a.positional == (1, 2)
@test a.keyword == Dict{Symbol, Any}(:a => 3, :b => 4)
```
"""
collect_arguments(positional...; keyword...) = @chain begin
    keyword
    Dict
    Arguments(positional, _)
end

export collect_call
"""
    collect_call(f::Function, positional...; keyword...)

Easy way to build a [`LazyCall`](@ref) type.

# Examples
```julia
l = collect_call(vcat, [1, 2], [3, 4])
@test l.function_call == vcat
@test l.arguments == collect_arguments([1, 2], [3, 4])
```
"""
collect_call(f, positional...; keyword...) = @chain begin
    keyword
    Dict
    Arguments(positional, _)
    LazyCall(_, f)
end

import Base.==

==(a::Arguments, b::Arguments) =
    (a.positional == b.positional) && (a.keyword == b.keyword)

==(a::LazyCall, b::LazyCall) =
    (a.function_call == b.function_call) && (a.arguments == b.arguments)

export run
"""
     run(a::Arguments)

Call `run` on the [`Arguments`](@ref) in `a`

# Examples
```julia
run_test = @chain begin
    collect_arguments([1, 2], [3, 4])
    unshift(_, vcat)
    collect_arguments(_, map)
    run
end

@test run_test == map(vcat, [1, 2], [3, 4])
```
"""
Base.run(a::Arguments) = run(a, run)

"""
     run(l::LazyCall)

Call `l.function_call` on the [`Arguments`](@ref) in `l.arguments`

# Examples
```julia
run_test = @chain begin
    collect_arguments([1, 2], [3, 4])
    unshift(_, vcat)
    LazyCall(_, map)
    run
end

@test run_test == map(vcat, [1, 2], [3, 4])
```
"""
Base.run(l::LazyCall) = run(l.arguments, l.function_call)

"""
     run(a::Arguments, f::Function)

Call `f` on the [`Arguments`](@ref) in `a`

# Examples
```julia
run_test = @chain begin
    collect_arguments([1, 2], [3, 4])
    unshift(_, vcat)
    run(_, map)
end

@test run_test == map(vcat, [1, 2], [3, 4])
```
"""
Base.run(a::Arguments, f::Function) = f(a.positional...; a.keyword...)

"""
     run(l::LazyCall, f::Function)

Insert `l.function_call` as the first positional argument in
`l.arguments`, the standard position for functional programming,
then `run` `f` on the result.

# Examples
```julia
run_test = @chain begin
    collect_arguments([1, 2], [3,4])
    LazyCall(_, vcat)
    run(_, map)
end

@test run_test == map(vcat, [1, 2], [3, 4])
```
"""
Base.run(l::LazyCall, f::Function) = @chain begin
    l.arguments
    unshift(_, l.function_call)
    run(_, f)
end

export lazy_call
"""
    @lazy_call(e)

Will break apart a function call into a [`LazyCall`](@ref) object.

# Examples
```julia
test_function(arguments...; keyword_arguments...) =
    (arguments, keyword_arguments)

@test ( @lazy_call test_function(1, 2, a = 3) ) ==
    collect_call(test_function, 1, 2, a = 3)
```
"""
lazy_call(e) =
    MacroTools.@match e begin
        a_(b__) => Expr(:call, :collect_call, a, b...)
        a_ => a
    end

@nonstandard lazy_call
export @lazy_call


text:
export Material
export matparam, kottke_avg_param
import Base:string

# Note that the output of tensorize is always a matrix: for K = 1, it is a 1×1 matrix.
tensorize(x::Number) = tensorize(@SVector(fill(complex(x), 3)))  # SComplexF{3}(x,x,x): vector
tensorize(v::AbsVecNumber) = diagm(Val(0)=>SComplexF{3}(v))  # complex(v) is allocating, so use SComplexF{K}(v)
tensorize(m::AbsMatNumber) = S²ComplexF3(m)  # complex(m) is allocating, so use S²ComplexF{K,K*K}(m)

struct Material
    name::String
    param::Tuple{S²ComplexF3,S²ComplexF3}  # (material parameter interacting with E, and H)
end

Material(name::String, param::Tuple2{MatParam}) = Material(name, tensorize.(param))
Material(name::String; ε::MatParam=1, μ::MatParam=1) = Material(name, (ε,μ))  # ε, μ: keyword arguments
# Consider implementing Material2 and Material3 for 2D and 3D solvers.

string(m::Material) = m.name
matparam(m::Material, ft::FieldType) = m.param[Int(ft)]

kottke_avg_param(param1::AbsMatNumber, param2::AbsMatNumber, n12::AbsVecReal, rvol1::Real) =
    (K = length(n12); kottke_avg_param(S²ComplexF{K,K^2}(param1), S²ComplexF{K,K^2}(param2), SFloat{K}(n12), rvol1))

# Implement the averaging scheme between two local material parameter tensors developed in
# the paper by Kottke, Farjadpour, Johnson entitled "Perturbation theory for anisotropic
# dielectric interfaces, and application to subpixel smoothing of discretized numerical
# methods", Physical Review E 77 (2008): 036611.
function kottke_avg_param(param1::S²ComplexF{K}, param2::S²ComplexF{K}, n12::SFloat{K}, rvol1::Real) where {K}
    Scomp = @SMatrix rand(K,K-1)  # directions complementary to n12; works even for K = 1
    Stemp = [n12 Scomp]  # SMatrix{K,K}
    S = qr(Stemp).Q  # nonallocating; 1st column is normalized n12

    τ1 = τ_trans(transpose(S) * param1 * S)  # express param1 in S coordinates, and apply τ transform
    τ2 = τ_trans(transpose(S) * param2 * S)  # express param2 in S coordinates, and apply τ transform

    τavg = τ1 .* rvol1 + τ2 .* (1-rvol1)  # volume-weighted average

    return S * τ⁻¹_trans(τavg) * transpose(S)  # apply τ⁻¹ and transform back to global coordinates
end

# Kottke's averaging scheme reduced for the field dimension orthogonal to the shape
# dimension (such as Hz in the TE mode).  In that case, the field is always parallel to the
# shape boundaries, so the averaging scheme reduces to simple arithmetic averaging.  For
# example for the case with Shape{2} and Kf = 1 (like Hz in the TE mode), the material
# parameter is an 1×1 tensor (= scalar) and therefore the situation is the same as averaging
# isotropic ε between E-field voxels, which is arithmetic averaging.  The only other case is
# the case with Shape{1} and Kf = 1  (like Ex in the slapes in the z-direction).  We can
# easily prove that in these cases again the averaging reduces to simple arithmetic
# averaging between material parameter tensors.  See Agenda > MaxwellBase > Feb/22/2020.
kottke_avg_param(param1::AbsMatNumber, param2::AbsMatNumber, rvol1::Real) =
    param1 .* rvol1 + param2 .* (1-rvol1)

# Equation (4) of the paper by Kottke et al.  Need to verify for K = 1 and 2.
function τ_trans(ε::S²ComplexF3)
    ε₁₁, ε₂₁, ε₃₁, ε₁₂, ε₂₂, ε₃₂, ε₁₃, ε₂₃, ε₃₃ = ε

    return S²ComplexF3(
        -1/ε₁₁, ε₂₁/ε₁₁, ε₃₁/ε₁₁,
        ε₁₂/ε₁₁, ε₂₂ - ε₂₁*ε₁₂/ε₁₁, ε₃₂ - ε₃₁*ε₁₂/ε₁₁,
        ε₁₃/ε₁₁, ε₂₃ - ε₂₁*ε₁₃/ε₁₁, ε₃₃ - ε₃₁*ε₁₃/ε₁₁
    )
    # t = zeros(3,3)
    # t[2:3, 2:3] = s[2:3, 2:3]
    # s11 = s[1,1]
    # s[1,1] = -1
    #
    # t = t - (s(:,1) * s(1,:)) ./ s11
    # return t
end

function τ_trans(ε::S²ComplexF2)
    ε₁₁, ε₂₁, ε₁₂, ε₂₂ = ε

    return S²ComplexF2(
        -1/ε₁₁, ε₂₁/ε₁₁,
        ε₁₂/ε₁₁, ε₂₂ - ε₂₁*ε₁₂/ε₁₁
    )
end

τ_trans(ε::S²ComplexF1) = -1 ./ ε  # S²ComplexF1

# function τ_trans(ε::S²ComplexF{K}) where {K}
#     τ_temp1 = @SMatrix [ε[i,j] for i = 2:K, j = 2:K]  # cannot use type parameter K here: error is generated
#     τ_temp2 = [@SMatrix(zeros(ComplexF,1,K-1)); τ_temp1]
#     τ_temp = [@SVector(zeros(ComplexF,K)) τ_temp2]
#
#     σcol_temp = @SVector [ε[i,1] for i = 2:K]
#     σcol = [SVec{1,ComplexF}(-1); σcol_temp]
#
#     σrow_temp = @SVector [ε[1,j] for j = 2:K]
#     σrow = [SVec{1,ComplexF}(-1); σrow_temp]
#
#     τ = τ_temp - (σcol * transpose(σrow)) ./ ε[1,1]
#
#     return τ
# end

# Equation (23) of the paper by Kottke et al.
function τ⁻¹_trans(τ::S²ComplexF3)
    τ₁₁, τ₂₁, τ₃₁, τ₁₂, τ₂₂, τ₃₂, τ₁₃, τ₂₃, τ₃₃ = τ

    return S²ComplexF3(
        -1/τ₁₁, -τ₂₁/τ₁₁, -τ₃₁/τ₁₁,
        -τ₁₂/τ₁₁, τ₂₂ - τ₂₁*τ₁₂/τ₁₁, τ₃₂ - τ₃₁*τ₁₂/τ₁₁,
        -τ₁₃/τ₁₁, τ₂₃ - τ₂₁*τ₁₃/τ₁₁, τ₃₃ - τ₃₁*τ₁₃/τ₁₁
    )

    # s = zeros(3,3)
    # s[2:3, 2:3] = t[2:3, 2:3]
    # t11 = t[1,1]
    # t[1,1] = 1
    #
    # s = s - (t(:,1) * t(1,:)) ./ t11
    # return s
end

function τ⁻¹_trans(τ::S²ComplexF2)
    τ₁₁, τ₂₁, τ₁₂, τ₂₂ = τ

    return S²ComplexF2(
        -1/τ₁₁, -τ₂₁/τ₁₁,
        -τ₁₂/τ₁₁, τ₂₂ - τ₂₁*τ₁₂/τ₁₁
    )
end

τ⁻¹_trans(τ::S²ComplexF1) = -1 ./ τ  # S²ComplexF1


text:
module UGM
  # exports
  export EdgeStruct

  export configPotential

  export sampleDiscrete
  export sampleExact

  export decodeExact

  export inferExact

  # includes

  include("misc/edgeStruct.jl")
  include("misc/configPotential.jl")

  include("sample/sampleDiscrete.jl")
  include("sample/sampleExact.jl")

  include("decode/decodeExact.jl")

  include("infer/inferExact.jl")

end


text:
include("1okp.jl")

"""

    Le dernier bug (j'éspère !)

    Quand on créer un sous problème avec moins de variable (B&B), si on veut étudier cette solution par rapport aux autres
    Il faut absolument lui rajouter toutes les variables qu'on lui a enlevé, Faut vraiment être con pour avoir oublier ca
    (Oups !)

"""

"""
    return the λ scalarization of the given problem prob

    @prob : MOKP
    @λ : Vector{Float64}
"""
function weightedScalarRelax(prob::Problem, λ::Vector{Float64})
    @assert length(λ) == prob.nbObj "Le vecteur λ ne convient pas"

    obj = Vector{Float64}(undef, prob.nbVar) # new merged objective
    # GOAL : compute the coef of each variable in the new objective by merging all the objectives
    for iterVar = 1:prob.nbVar
        obj[iterVar] = sum([λ[iter] * prob.objs[iter].profits[iterVar] for iter = 1:prob.nbObj])
    end
    # construction of the problem
    return Problem(
        1,
        prob.nbVar,
        [Obj(obj)],
        prob.constraint
    )
end

"""
    return y the point associated with the solution x and the problem prob
"""
function evaluate(prob::Problem, x::Vector{Bool})
    y = zeros(Float64, prob.nbObj)
    # GOAL : computing the image of x by prob
    for iterObj = 1:prob.nbObj
        for iter = 1:prob.nbVar
            if x[iter]
                y[iterObj] += prob.objs[iterObj].profits[iter]
            end
        end
    end
    # the resulting point
    return Solution(x, y)
end

"""
    whichFathomed(upperBound::DualSet, lowerBound::Vector{Solution}, S::Vector{Solution}, consecutivePoint::Vector{Tuple{Solution, Solution}})

    return the way the subproblem is fathomed : none, infeasibility, optimality, dominance
"""
function whichFathomed(upperBound::DualSet, lowerBound::Vector{Solution}, S::Vector{Solution}, consecutivePoint::Vector{Tuple{Solution, Solution}})
    noNadirUnderUB = true # becomes false if we found a nadir point under the upperBound

    if length(lowerBound) == 0 # no solutions supported
        return infeasibility
    elseif length(lowerBound) == 1 # only one feasible solution
        return optimality
    else
        iterConsecPoint = 0
        while noNadirUnderUB && iterConsecPoint < length(consecutivePoint) # going through the consecutive points
            iterConsecPoint += 1
            solR, solL = consecutivePoint[iterConsecPoint] # getting the two consecutive points
            nadirPoint = [min(solR.y[ind], solL.y[ind]) for ind = 1:length(solR.y)] # constructing the nadir point
            # projection of the nadir point on the constraints of the upper bound
            Ax = upperBound.A * nadirPoint
            # verifying that the nadir point is under each constraint of the upperBound
            indexConstr = 0
            while noNadirUnderUB && indexConstr < length(upperBound.b)
                indexConstr += 1
                noNadirUnderUB = noNadirUnderUB && Ax[indexConstr] > upperBound.b[indexConstr]
            end
        end
        # conclusion on the type of pruning
        if noNadirUnderUB
            return dominance
        else
            return none
        end
    end
end

"""
    updateBounds!(S::Vector{Solution}, consecutiveSet::Vector{Tuple{Solution, Solution}}, lowerBoundSub::Vector{Solution})

    update upper et lower bound sets according to the new feasible subproblem

    @S : Solution
    @consecutiveSet : Set of tuples of consecutive points
    @lowerBoundSub : Lower bound set for the subproblem
"""

""" LUCAS J'AI DES QUESTIONS SUR LA METHODE GENERALE ???"""
function updateBounds!(S::Vector{Solution}, consecutiveSet::Vector{Tuple{Solution, Solution}}, lowerBoundSub::Vector{Solution})

    # variables used to memorize where to delete items in consecutiveSet (points that are no longer consecutive because we added some points between them)
    indSuppr       = Vector{Int}(undef, length(consecutiveSet))
    indEndIndSuppr = 0

    for subSol in lowerBoundSub
        for iter = 1:length(consecutiveSet)
            solR, solL = consecutiveSet[iter]
            if subSol.y[1] > min(solR.y[1], solL.y[1]) && subSol.y[2] > min(solR.y[2], solL.y[2]) # subSol dominates the nadir point of solR and solL
                push!(S, subSol) # subSol is a solution
                push!(consecutiveSet, (solR, subSol)) # we update consecutiveSet
                push!(consecutiveSet, (subSol, solL))

                indEndIndSuppr += 1
                indSuppr[indEndIndSuppr] = iter
            end
        end
        # now deleting the previous consecutive points that are no longer consecutive
        for iter = 1:indEndIndSuppr
            deleteat!(consecutiveSet, indSuppr[iter]-iter+1)
        end
        # reinitializing the variables
        indSuppr = Vector{Int}(undef, length(consecutiveSet)) # maybe there is a better way of doing the initialization... idk
        indEndIndSuppr = 0
    end
end

"""
    newAssignments(A::Vector{Int},i::Int)

    return a tuple of the two vectors : assignments for the two subproblems.
    For the first (copyA), the variable i has been assigned to zero, for the second (A), i has been assigned to one
"""
function newAssignments(A::Vector{Int},i::Int)
    """ VRAIE COPIE ????"""
    copyA = A[1:end]
    copyA[i] = 0
    A[i] = 1
    return copyA, A
end

"""

    En cours de Travaux

"""

"""
    computeBoundDicho(prob::Problem, assignment::Vector{Int}, indEndAssignment::Int, assignmentWeight::Float64, assignmentProfit::Vector{Float64}, ϵ::Float64; verbose = false)

    computes bound dicho

    @prob : current problem
    @assignment : assignment of variables (non ordered) (length = nbVars)
    @indEndAssignment : index of the the last variable (non ordered) assigned in the problem
    @assignmentWeight : sum of the weights of the variables assigned to one
    @assignmentProfit : sum of the profits of the variables assigned to one
    @ϵ : small value to compute the vector lambda and compute the weighted scalarization

    @XSEm - Vector of solutions : memorize all the supported efficient and extreme solutions
    @consecutiveSet : memorize tuples of points that are consecutive
    @S : Queue of tuples of consecutive points to study
"""
function computeBoundDicho(prob::Problem, assignment::Vector{Int}, indEndAssignment::Int, assignmentWeight::Float64, assignmentProfit::Vector{Float64}, ϵ::Float64; verbose = false)

    XSEm = Vector{Solution}() # memorize all the supported efficient and extreme solutions
    consecutiveSet = Vector{Tuple{Solution, Solution}}() # memorize tuples of points that are consecutive
    S = Vector{Tuple{Solution, Solution}}() # queue of solutions to study

    # computing the two first points, sol1 is the best for the only the objective number 1, and sol2 for the objective number 2.
    sol1 = solve1OKPMain(weightedScalarRelax(prob, [1., ϵ]), assignment, indEndAssignment, assignmentWeight, sum([1., ϵ] .* assignmentProfit), verbose = verbose)
    sol2 = solve1OKPMain(weightedScalarRelax(prob, [ϵ, 1.]), assignment, indEndAssignment, assignmentWeight, sum([ϵ, 1.] .* assignmentProfit), verbose = verbose)

    # evaluating the two sols and constructing the associated solutions
    sol1 = evaluate(prob, sol1.x)
    sol2 = evaluate(prob, sol2.x)

    # computing the maximum of each objective
    max1 = sol1.y[1] # maximum of objective 1
    max2 = sol2.y[2] # maximum of objective 2

    verbose && println("On crée les solution lexico : $(sol1.y) et $(sol2.y)")

    if sol1.y == sol2.y # solutions are adenticals
        push!(XSEm, sol1)
    else # solutions are different
        push!(S, (sol1, sol2)) # we'll need to study (sol1,sol2)
        push!(XSEm, sol1, sol2) # sol1 and sol2 are supported efficient and extreme solutions

        # goal : compute the dichotomy
        while length(S) != 0 # while we have solutions to study
            solR, solL = pop!(S) # get the two consecutive points from the queue of tuples to study

            verbose && println("On étudie la paire de solution : $(solR.y) et $(solL.y)")
            # computing the direction of the search
            λ = [solL.y[2]-solR.y[2], solR.y[1]-solL.y[1]]
            # computing the resulting solution
            solE = solve1OKPMain(weightedScalarRelax(prob, λ), assignment, indEndAssignment, assignmentWeight, sum(λ .* assignmentProfit), verbose = false)
            solE = evaluate(prob, solE.x)

            """ SEE BELOW ??? """
            if sum(λ .* solE.y) > sum(λ .* solR.y) # solE is better than solR according to λ
                push!(XSEm, solE) # solE is solution we want
                push!(S, (solR, solE), (solE, solL)) # now we need to study (solR, solE) and (solE, solL)
                verbose && println("On a trouvé la solution : $(solE.y)")
            else # solE is equal to solR according to λ, so we didn't find new solutions
                push!(consecutiveSet, (solR, solL)) # we know we won't find solutions between solR, solL, so we memorize the tuple
                verbose && println("On a rien trouvé dans cette direction")
            end
        end
    end

    # GOAL : compute dual bound set
    # we represent the dual bound set not with points because there is an infinity of them but with constraints
    # the dual bound : A.x = B
    A = Array{Float64, 2}(undef, length(consecutiveSet) + 2, 2)
    B = Vector{Float64}(undef, length(consecutiveSet)+2)
    # construction of the maxtrix A and the vector B
    for iter = 1:length(consecutiveSet)
        a, b = consecutiveSet[iter]
        a, b = a.y, b.y
        A[iter, 1] = b[2] - a[2]
        A[iter, 2] = a[1] - b[1]

        B[iter] = (b[2] - a[2]) * a[1] - (b[1] - a[1]) * a[2]
    end

    # finishing the computation by adding the max1, max2 values
    A[end-1, 1] = 1
    A[end-1, 2] = 0
    B[end-1] = max1

    A[end, 1] = 0
    A[end, 2] = 1
    B[end] = max2

    return XSEm, consecutiveSet, DualSet(A, B)
end

"""
    branchAndBound(prob::Problem, assignment::Vector{Int}, assignmentWeight::Float64, assignmentProfit::Vector{Float64}, S::Vector{Solution}, consecutiveSet::Vector{Tuple{Solution, Solution}}, indEndAssignment::Int = 0, ϵ::Float64 =0.01; verbose = false)

"""
function branchAndBound(prob::Problem, assignment::Vector{Int}, assignmentWeight::Float64, assignmentProfit::Vector{Float64}, S::Vector{Solution}, consecutiveSet::Vector{Tuple{Solution, Solution}}, indEndAssignment::Int = 0, ϵ::Float64 =0.01; verbose = false)

    # computing the dichotomy for : LB, UB and consecutivePointsSet
    lowerBoundSub, consecutivePointSub, upperBoundSub = computeBoundDicho(prob, assignment, indEndAssignment, assignmentWeight, assignmentProfit, ϵ, verbose = verbose)

    # computing the type of pruning for the subproblem
    fathomed::Fathomed = whichFathomed(upperBoundSub, lowerBoundSub, S, consecutiveSet)

    # stop search because we found an optimal solution
    if fathomed == optimality
        updateBounds!(S, consecutiveSet, lowerBoundSub)
    end

    # keep branching
    if fathomed == none && indEndAssignment < prob.nbVar
        A0, A1 = newAssignments(assignment, indEndAssign]ment+1)

        verbose && println("Branch and Bound sur la variable $(indEndAssignment+1), on la fixe à 0")

        branchAndBound(prob, A0, assignmentWeight, assignmentProfit, S, consecutiveSet, indEndAssignment+1, ϵ, verbose = verbose) # exploring the first subproblem

        verbose && println("Branch and Bound sur la variable $(i+1), on la fixe à 1")

        variableProfit = broadcast(obj->obj.profits[indEndAssignment+1], prob.objs)

        branchAndBound(prob, A1, assignmentWeight + prob.constraint.weights[indEndAssignment+1], assignmentProfit + variableProfit, S, consecutiveSet, indEndAssignment+1, ϵ, verbose = verbose) # exploring the second subproblem
    end
end

"""
    reorderVariable(prob::Problem, reorderVect::Vector{Int})

    returns a new problem ordered with reorderVect
"""
function reorderVariable(prob::Problem, reorderVect::Vector{Int})
    @assert prob.nbVar==length(reorderVect) "Wrong argument for the function reorderVariable"

    return Problem(
        prob.nbObj,
        prob.nbVar,
        [Obj(prob.objs[iter].profits[reorderVect]) for iter=1:prob.nbObj],
        Const(prob.constraint.maxWeight, prob.constraint.weights[reorderVect])
    )
end

"""

    Attention à l'ordre des variables !

"""

"""
    main_BranchandBound(prob::Problem, orderName = "random", ϵ::Float64 = 0.01 ; verbose = false)

    returns the set of non dominated points coputed with the multi objective branch and bound

    @S : XSEm at first, then completed with the branch and bound to become XE
    @consecutivePoints = vector of tiples of consecutive points (the points are the points in S)
"""
function main_BranchandBound(prob::Problem, orderName = "random", ϵ::Float64 = 0.01 ; verbose = false)

    permVect, revPermVect = permOrder(prob, orderName) # order variables accordingly
    auxProb = reorderVariable(prob, permVect) # creating the new ordered problem

    # computing a first dichotomy to obtain a lower bound (the XSEm set, here S)
    S, consecutivePoints = computeBoundDicho(auxProb, Vector{Int}(), 0, 0., zeros(Float64, prob.nbObj), ϵ, verbose = verbose)[1:2]

    # initializing the assignment vector
    assignment = Vector{Int}(undef, prob.nbVar)
    for iter=1:prob.nbVar
        assignment[iter] = -1
    end

    # computing the branch and bound
    branchAndBound(auxProb, assignment, 0., zeros(Float64, prob.nbObj),  S, consecutivePoints, 0, ϵ, verbose = verbose)

    # reorder variables inside items of S
    S = broadcast(sol->Solution(sol.x[revPermVect], sol.y), S)

    return S
end


text:
using Plots

# Plot of linked scalar sequences
x = 0:0.001:5
s1 = sin.(x)'
s2 = sin.(log.(range(1;length=5001, stop=exp(5))))'
l, m = linkseq(s1', s2', 1000, 50)

p = plot(x, [s1[:] s2[:]])
plot!(p, l, m; legend=false)
savefig(p, "docs/assets/2d_linked_seqs.svg")

# Plot of linked vector sequences
x = 0:0.001:5
s1 = [0.3; -0.7] * sin.(x)'
s2 = [0.7; 0.3] * sin.(log.(range(1;length=5001, stop=exp(5))))'
l, m = linkseq(s1, s2, 1000, 50)

p = plot(x, [s1[1,:] s2[1,:]], [s1[2,:] s2[2,:]])
plot!(p, l, m[1,:,:], m[2,:,:]; legend=false)
savefig(p, "docs/assets/3d_linked_seqs.svg")

# Plot different linked sequences
s1 = [0, 1, 2, 4, 5, 2, 3, 1, 0, 0]
s2 = [2, 4, 5, 6, 6, 5, 5, 4, 3, 1, 0, 0]
l, m = linkseq(s1, s2)
p = plot(s1)
plot!(p, s2)
plot!(p, l.+1, m; color=:black, legend=false)
savefig(p, "docs/assets/different_seqs.svg")

max_stars_repo_name: gabriel-ss/DynamicTimeWarp

max_issues_repo_name: gabriel-ss/DynamicTimeWarp

max_forks_repo_name: gabriel-ss/DynamicTimeWarp


text:
function neworder_phylo(args...; n::Int64=Int[], e1::Array{Int64, 1}=Int[], e2::Array{Int64, 1}=Int[],
	N::Int64=Int[], neworder::Array{Int64, 1}=Int[], order::Int64=Int[])
#   n: nb of tips
#   m: nb of nodes
#   N: nb of edges
#   degrmax is the largest value that a node degree can be

#* L is a 1-d array storing, for each node, it indices of the rows of
#  the edge matrix where the node is ancestor (i.e., present in the 1st
#  column). It is used in the same way than a matrix (which is actually
#  a vector) is used in R as a 2-d structure.
#* pos gives the position for each 'row' of L, that is the number of elements
#  which have already been stored for that 'row'.
    m = N - n + 1
    degrmax = n - m + 1

    L = zeros(Int, m * degrmax)
    pos = zeros(Int, m)

#* we now go down along the edge matrix
    for i = 1:N
    	k = e1[i] - n
    	j = pos[k]
    	pos[k] = pos[k] + 1
    	L[k + m * j] = i - 1
    end

#* L is now ready: we can start the recursive calls.
#* We start with the root 'n + 1':

    if order == 1
    	iii = 1
    	 foo_reorder(node=n+1,n=n,m=m,e1=e1,e2=e2,neworder=neworder,L=L,pos=pos)
    elseif order == 2
    	iii = N
    	bar_reorder(node=n+1,n=n,m=m,e1=e1,e2=e2,neworder=neworder,iii=iii,L=L,pos=pos)
    end
end


text:
#
# Power generation expansion problem
#
using SDDP, GLPKMathProgInterface, JuMP, SDDiP, Ipopt

struct GenTech
    capacity::Float64   # maximum MW this generator technology can be built for
    b_cost::Float64     # cost to build
    g_cost::Float64     # cost to use
    nunits::Int         # maximum number of units we are happy to build
    init::Vector{Int}   # initial state
end
struct GEPData
    T::Int                                  # planning duration
    demand::Array{Float64,2}                # power demand in each stage
    penalty::Float64                        # penalty per MW not met
    hours::Float64                          # hours in each stage
    rho::Float64                            # discounting rate
    gentypes::Dict{Symbol, GenTech}         # all technologies available
end

function makedata()
    T = 5
    S = 8
    inv_c = 1e4 * ones(T)
    gen_c = 4 * ones(T)
    # Some random demand
    d = [ 5   5  5   5   5   5   5   5;
          4   3  1   3   0   9   8  17;
          0   9  4   2  19  19  13   7;
         25  11  4  14   4   6  15  12;
          6   7  5   3   8   4  17  13]
    p = 5e5
    hours = 1.0 # scaled for nice numbers...
    gentypes = Dict(:TechA => GenTech(1.0, 1e4, 4.0, 20.0, zeros(20)))
    GEPData(T, d, p, hours, 0.99, gentypes)
end

data = makedata()
gentypes = collect(keys(data.gentypes))

getgen(s::Symbol) = data.gentypes[s]
maxgen(s::Symbol) = getgen(s).capacity
nunits(s::Symbol) = getgen(s).nunits
build(s::Symbol)  = getgen(s).b_cost
use(s::Symbol)    = getgen(s).g_cost
init(s::Symbol)   = getgen(s).init

m=SDDPModel(stages=data.T, objective_bound=0.0, sense=:Min, solver=GLPKSolverMIP()) do sp, stage
    @binarystate(sp, invested[i = gentypes, j = 1:nunits(i)], invested0 == init(i)[j], Bin)
    @variables(sp, begin
        generation[gentypes] >= 0
        penalty >= 0
        demand
    end)

    @constraints(sp, begin
        # Can't un-invest
        investment[i = gentypes, j = 1:nunits(i)], invested[i, j] >= invested0[i, j]
        # Generation capacity
        [i = gentypes], sum(invested[i, j] * maxgen(i) for j = 1:nunits(i)) >= generation[i]
        # Meet demand or pay a penalty
        penalty >= demand - sum(generation)
        # Order the units to break symmetry, units are identical
        [i = gentypes, j = 1:nunits(i)-1], invested[i, j] <= invested[i, j+1]
    end)

    # Demand is uncertain
    @rhsnoise(sp, D=data.demand[stage,:], demand == D)

    # Handy calculation of the investment cost in this stage
    @expression(sp, investment_cost[i = gentypes], build(i) * sum(invested[i, j] - invested0[i, j] for j = 1:nunits(i)))

    @stageobjective(sp,
        sum(investment_cost[i] for i = gentypes) * data.rho ^ (stage - 1) +
        sum(use(i) * generation[i] for i = gentypes) * data.hours * data.rho ^ (stage - 1) +
        penalty * data.penalty * data.hours)

    # Solve with the level method as the Lagrangian solver
    setSDDiPsolver!(sp, method = KelleyMethod(),
                        pattern = Pattern(benders=1, lagrangian=5, strengthened_benders=1),
                        LPsolver = GLPKSolverLP()
                        )

end

srand(11111)
status = solve(m, max_iterations=60)
@assert isapprox(getbound(m), 460_533.0, atol=1e3)


text:
#Reference: https://jump.dev/JuMP.jl/0.17/refexpr.html
using CSV, DelimitedFiles, DataFrames
using JuMP, CPLEX

#Preparing optimization model

const N_GEN = 145
const N_HRS = 24
const INITIAL_DAY = 1

dfGenerator = CSV.read(".//inputs//data_generators.csv", DataFrame)
dfDemand = CSV.read(".//inputs//demand_reserves.csv", DataFrame)

for day in 1:2
  local dfDayDemand = dfDemand[(day-1)*N_HRS+1:(N_HRS*day)+1, [:Hour, :Demand, :Reserves]];

    model = Model(CPLEX.Optimizer)
    set_optimizer_attribute(model, "CPX_PARAM_EPINT", 1e-5)

    @variable(model, genOnOff[1:N_GEN, 1:N_HRS], Bin)
    @variable(model, genOut[1:N_GEN, 1:N_HRS] >= 0)
    @variable(model, genStartUp[1:N_GEN, 1:N_HRS], Bin)
    @variable(model, genShutDown[1:N_GEN, 1:N_HRS], Bin)

    #Setting the objective function
    @objective(model, Min, sum(sum(dfGenerator.FixedCost[g]*genOnOff[g,t] + dfGenerator.VariableCost[g]*genOut[g,t] + dfGenerator.StartUpCost[g]*genStartUp[g,t] + dfGenerator.ShutdownCost[g]*genShutDown[g, t] for g in 1:N_GEN) for t in 1:N_HRS) )

    #Initial constraints
    #Generators should have an initial and defined value t=0
    @constraint(model, constInit[g=1:N_GEN], genOnOff[g,1] == dfGenerator.Uini[g])

    #Initial values of Generator shutdown and startup values should be defined
    @constraint(model, constShutDn[g=1:N_GEN], genShutDown[g,1] == 0)
    @constraint(model, constShutUp[g=1:N_GEN], genStartUp[g,1] == 0)

    #Demand constraint
    #for t in 1:N_HRS
    @constraint(model, constDemand[t=1:N_HRS], sum(genOut[g, t] for g=1:N_GEN) == dfDayDemand.Demand[t])
    #end

    #Power bounds
    @constraint(model, constMaxPow[t=1:N_HRS, g=1:N_GEN],  genOut[g,t] <= dfGenerator.MaxPowerOut[g]*genOnOff[g,t] )
    @constraint(model, constMinPow[t=1:N_HRS, g=1:N_GEN],  genOut[g,t] >= dfGenerator.MinPowerOut[g]*genOnOff[g,t] )

    #Logical Conditions
    @constraint(model, constLog1[t=2:N_HRS, g=1:N_GEN], (genStartUp[g,t] - genShutDown[g,t]) == (genOnOff[g,t] - genOnOff[g,t-1]))
    @constraint(model, constLog2[t=1:N_HRS, g=1:N_GEN], genStartUp[g,t] + genShutDown[g, t] <= 1)

    #Ramp limits
    @constraint(model, constRampUp[t=2:N_HRS, g=1:N_GEN], genOut[g,t] - genOut[g,t-1] <= dfGenerator.RampUpLimit[g]*genOnOff[g, t-1] + dfGenerator.RampStartUpLimit[g]*genStartUp[g,t])
    @constraint(model, constRampDown[t=2:N_HRS, g=1:N_GEN], (genOut[g,t-1] - genOut[g,t]) <= (dfGenerator.RampDownLimit[g]*genOnOff[g,t] + dfGenerator.RampShutDownLimit[g]*genShutDown[g,t]))

    #Reserves
    @constraint(model, constReserves[t=1:N_HRS], sum(dfGenerator.MaxPowerOut[g]*genOnOff[g,t] for g=1:N_GEN) >= (dfDayDemand.Demand[t] + dfDayDemand.Reserves[t]))

    mystatus = optimize!(model)

    #Printing general results
    println("------------------------------------")
    println("------- OBJECTIVE VALUE -------")
    println("Objective value: ", JuMP.objective_value(model))

  println("------- DATA TYPE -------")
    print(typeof(JuMP.value.(genOut)))
    genOutTrans = transpose(JuMP.value.(genOut))

    open(".//outputs//resultsTrans.csv", "w") do io
                  writedlm(io, genOutTrans, ',')
            end;

end

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
using HPAT

@acc hpat function cumsum_test(file_name, file_name2)
    df1 = DataSource(DataTable{:id=Int64,:x=Float64,:y=Float64}, HDF5, file_name)
    df2 = DataSource(DataTable{:id=Int64,:x1=Float64,:y1=Float64}, HDF5, file_name2)
    df3 = join(df1, df2, :id==:id, :id)
    return sum(df3[:x])
end

cumsum_test(ENV["HOME"]*"/tmp/udf.hdf5", ENV["HOME"]*"/tmp/udf2.hdf5")

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
# main.jl - Animated Particles
# Modified version of Luxor
push!(LOAD_PATH, "C:\\Users\\mcbaron\\Documents\\Julia")
using Luxor
# Include helper functions
include("deps.jl")

# Constants
num_particles = 50
bounds = [1360 768]
vidlength = 5 # sec
framerate = 32
numFrames = vidlength * framerate

# Instantiate num_particles random particles
liveParts = genScreen(num_particles, bounds)

screenFrames = Movie(bounds[1], bounds[2], "Particles", 0:numFrames)

# animate the whole thing
animate(screenFrames, [
    Scene(screenFrames, backdrop, 0:numFrames)
    Scene(screenFrames, drawScreen, 0:numFrames, lineartween, liveParts)
], creategif=true)

# As of 20190102 this works, but the particles seem to settle as they evolve.

max_stars_repo_name: mcbaron/procedural-art

max_issues_repo_name: mcbaron/procedural-art

max_forks_repo_name: mcbaron/procedural-art


text:
function EvaluateTransferFunction(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    ConvolvedDensity::ConvolvedDensity,
    κTransferFunction::κTransferFunction,
    PowerSpectrum::PowerSpectrum)
    for iidx in 1:length(ConvolvedDensity.ZBinArray)-1
        FFTLog = FFTLog(XArray = BackgroundQuantities.rZArray, FXArray =
        κTransferFunction.LensingSourceFunction.SourceFunctionArray[iidx, :] .*
        PowerSpectrum.GrowthFactor)
        Kl, κTransferFunction.TransferFunctionArray[iidx, :, :] =
        CosmoCentral.EvaluateFFTLog(FFTLog, CosmologicalGrid.ℓBinCenters)
        κTransferFunction.TransferFunctionArray[iidx, :, :] =
        κTransferFunction.TransferFunctionArray[iidx, :, :] ./
        Kl ./ Kl .* CosmologicalGrid.ℓBinCenters .*
        (CosmologicalGrid.ℓBinCenters .+ 1)
        CosmologicalGrid.KBeyondLimberArray = Kl
    end
end


text:
using LazilyInitializedFields
const LI = LazilyInitializedFields
using Documenter

using Test
@lazy mutable struct Foo{T}
    a::T
    @lazy b::Int
    @lazy c::Union{Float64, Nothing}
    @lazy d::Union{Int, Nothing}
    e::Float64
end
f = Foo{Int}(1, uninit, 2.0, uninit, 3.0)

@lazy struct Mut{T}
    a::T
    @lazy b::Int
end
m = Mut(1, uninit)

@testset "LazilyInitializedFields" begin
    @test f.a == 1
    @test_throws UninitializedFieldException f.b
    @test f.c == 2.0
    @test_throws UninitializedFieldException f.d
    @test f.e == 3.0
    @test f.e == 3.0

    # @test map(i->@isinit(f.$i), (:a, :b, :c, :d, :e)) == (true, false, true, false, true)
    rt = Core.Compiler.return_type
    @test rt((f -> f.a), Tuple{Foo{Int}}) == Int
    @test rt((f -> f.b), Tuple{Foo}) == Int
    @test rt((f -> f.c), Tuple{Foo}) == Union{Float64, Nothing}
    @test rt((f -> f.d), Tuple{Foo}) == Union{Int, Nothing}
    @test rt((f -> f.e), Tuple{Foo}) == Float64

    @test LI.islazyfield.(Foo, (:a, :b, :c, :d, :e)) == (false, true, true, true, false)

    @test_throws NonLazyFieldException @init! f.a = 2
    @test_throws NonLazyFieldException @isinit f.a
    @test_throws NonLazyFieldException @uninit! f.a
    @init! f.b = 2
    @test_throws AlreadyInitializedException @init! f.b = 2.0

    @test_throws AlreadyInitializedException @init! f.c = 2.0
    @test @isinit(f.c)
    f.c = nothing
    @test f.c == nothing
    @uninit! f.c
    @test !@isinit(f.c)
    @test_throws UninitializedFieldException f.d

    @test_throws LoadError @macroexpand @lazy a::Int

    @test_throws ErrorException m.a = 2
    @test_throws ErrorException m.b = 2
end

DocMeta.setdocmeta!(LazilyInitializedFields, :DocTestSetup, :(using LazilyInitializedFields))
doctest(LazilyInitializedFields; manual=false)


text:
export AbstractRefShape, RefSimplex, RefHypercube, RefLine, RefTriangle, RefTetrahedron,
       RefSquare, RefCube

abstract type AbstractRefShape end
struct RefSimplex{K}   <: AbstractRefShape end
struct RefHypercube{K} <: AbstractRefShape end

const RefTriangle    = RefSimplex{2}
const RefTetrahedron = RefSimplex{3}

const RefLine        = RefHypercube{1}
const RefSquare      = RefHypercube{2}
const RefCube        = RefHypercube{3}

include("gauss_quadrature.jl")
include("angular_quadrature.jl")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
@testset "Galilean trajectory unit tests..." begin
    #test box movement
    pos=[.5,.5]
    box=[0 1.
         0 1.]
    allow_d=[.1,.1]
    box_d=[1.,1.]

    allowed_pos, allowed_d=box_move(pos,allow_d,box)
    @test allowed_pos==pos+allow_d
    @test allowed_d == allow_d

    boxed_pos, boxed_d=box_move(pos, box_d, box)
    @test boxed_pos==[1.,1.]
    @test boxed_d==[.5,.5]

    box_reflect_v=box_reflect([.5,.99999999],box,[.1,.1],0.)
    @test box_reflect_v==[.1,-.1]

    #test boundary norm finding
    norm=1000
    v=[.5,.5]
    bn=boundary_norm(v, norm)
    @test bn == normalize(norm./v)

    #test reflection given norm
    rv=reflect(v,bn,0.)
    @test isapprox(rv, -v)

    #test reflection perturbaton
    pv=r_perturb(v,.01)
    @test v!=pv
end


text:
# Import the low-level interface module.
using HTSLib: htslib

"""
    iterateall(filepath::AbstractString)

Iterate over all records stored in `filepath` using the low-level interface.
"""
function iterateall(filepath::AbstractString)
    # Open the file.
    file = htslib.hts_open(filepath, "r")
    @assert file != C_NULL

    # Read the header.
    header = htslib.sam_hdr_read(file)
    @assert header != C_NULL

    # Allocate a placeholder for records.
    record = htslib.bam_init1()
    @assert record != C_NULL

    # Iterate over all records.
    n = 0
    res::Cint = 0
    while res ≥ 0
        res = htslib.sam_read1(file, header, record)
        @assert res ≥ -1
        res ≥ 0 && (n += 1)
    end
    println("# records: $(n)")

    # Destroy the allocated objects.
    htslib.bam_destroy1(record)
    htslib.sam_hdr_destroy(header)

    # Close the file.
    res = htslib.hts_close(file)
    @assert res == 0
end

# Script's entry point.
if PROGRAM_FILE == relpath(@__FILE__)
    @assert length(ARGS) == 1
    filepath = ARGS[1]
    iterateall(filepath)
end


text:
function JuMP.parse_constraint_call(_error::Function, vectorized::Bool,
                                            ::Val{:⊂}, lhs, rhs)
    _error("Unrecognized symbol ⊂ you mean ⊆ ?")
end
function JuMP.parse_constraint_call(_error::Function, vectorized::Bool,
                                            ::Val{:⊆}, lhs, rhs)
    parse_code = :()
    if vectorized
        build_call = :(JuMP.build_constraint.($_error, $(esc(lhs)), $(esc(:(SetProg.PowerSet.($rhs))))))
    else
        build_call = :(JuMP.build_constraint($_error, $(esc(lhs)), $(esc(:(SetProg.PowerSet($rhs))))))
    end
    return parse_code, build_call
end
function JuMP.parse_constraint_call(_error::Function, vectorized::Bool,
                                            ::Val{:⊃}, lhs, rhs)
    _error("Unrecognized symbol ⊃, did you mean ⊇ ?")
end
function JuMP.parse_constraint_call(_error::Function, vectorized::Bool,
                                            ::Val{:⊇}, lhs, rhs)
    parse_one_operator_constraint(_error, vectorized, Val(:⊆), rhs, lhs)
end


text:
# ---
# title: 1015. Smallest Integer Divisible by K
# id: problem1015
# author: Indigo
# date: 2021-01-30
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/smallest-integer-divisible-by-k/description/>
# hidden: true
# ---
#
# Given a positive integer `K`, you need to find the **length** of the
# **smallest** positive integer `N` such that `N` is divisible by `K`, and `N`
# only contains the digit `1`.
#
# Return _the **length** of _`N`. If there is no such `N`, return -1.
#
# **Note:** `N` may not fit in a 64-bit signed integer.
#
#
#
# **Example 1:**
#
#
#
#     Input: K = 1
#     Output: 1
#     Explanation: The smallest answer is N = 1, which has length 1.
#
#
# **Example 2:**
#
#
#
#     Input: K = 2
#     Output: -1
#     Explanation: There is no such positive integer N divisible by 2.
#
#
# **Example 3:**
#
#
#
#     Input: K = 3
#     Output: 3
#     Explanation: The smallest answer is N = 111, which has length 3.
#
#
#
#
# **Constraints:**
#
#   * `1 <= K <= 105`
#
#
## @lc code=start
using LeetCode

function smallest_repunit_div_by_k(K::Int)
    (K % 2 == 0 || K % 5 == 0) && (return -1)
    r, len = 0, 1
    while true
        r = (10 * r % K + 1) % K
        r == 0 && return len
        len += 1
    end
end
## @lc code=end


text:
"""
    FEMMDeforLinearAbstractBaseModule

Base module for operations on interiors of domains to construct system matrices and
system vectors for linear deformation models.
"""
module FEMMDeforLinearBaseModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: gradN!, nodesperelem, manifdim
import FinEtools.IntegDataModule: IntegData, integrationdata, Jacobianvolume
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.AssemblyModule: SysvecAssemblerBase, SysmatAssemblerBase, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
import FinEtools.FEMMBaseModule: FEMMAbstractBase, inspectintegpoints
import FinEtools.CSysModule: CSys, updatecsmat!
import FinEtools.DeforModelRedModule: nstressstrain, nthermstrain, Blmat!
import FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, add_btv!, locjac!, add_nnt_ut_only!
import FinEtools.MatDeforModule: rotstressvec
import FinEtools.SurfaceNormalModule: SurfaceNormal, updatenormal!
import LinearAlgebra: Transpose, mul!
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
import LinearAlgebra: norm, dot

abstract type FEMMDeforLinearAbstract <: FEMMAbstractBase end

function buffers(self::FEMMDeforLinearAbstract, geom::NodalField, u::NodalField)
    fes = self.integdata.fes
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    nstrs = nstressstrain(self.mr);  # number of stresses
    elmatdim = ndn*nne;             # dimension of the element matrix
    # Prepare buffers
    elmat = fill(zero(FFlt), elmatdim, elmatdim);      # element matrix -- buffer
    dofnums = zeros(FInt, elmatdim); # degree of freedom array -- buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    gradN = fill(zero(FFlt), nne, mdim); # intermediate result -- buffer
    D = fill(zero(FFlt), nstrs, nstrs); # material stiffness matrix -- buffer
    B = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    DB = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    elvecfix = fill(zero(FFlt), elmatdim); # vector of prescribed displ. -- buffer
    elvec = fill(zero(FFlt), elmatdim); # element vector -- buffer
    return dofnums, loc, J, csmatTJ, gradN, D, B, DB, elmat, elvec, elvecfix
end

"""
    mass(self::FEMMDeforLinearAbstract,  assembler::A,
      geom::NodalField{FFlt},
      u::NodalField{T}) where {A<:SysmatAssemblerBase, T<:Number}

Compute the consistent mass matrix

This is a general routine for the abstract linear-deformation  FEMM.
"""
function mass(self::FEMMDeforLinearAbstract,  assembler::A,  geom::NodalField{FFlt}, u::NodalField{T}) where {A<:SysmatAssemblerBase, T<:Number}
    fes = self.integdata.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdata);
    dofnums, loc, J, csmatTJ, gradN, D, B, DB, elmat = buffers(self, geom, u)  # Prepare buffers
    rho::FFlt = self.material.mass_density; # mass density
    NexpTNexp = FFltMat[];# basis f. matrix -- buffer
    ndn = ndofs(u)
    Indn = [i==j ? one(FFlt) : zero(FFlt) for i=1:ndn, j=1:ndn] # "identity"
    for j = 1:npts # This quantity is the same for all quadrature points
        Nexp = fill(zero(FFlt), ndn, size(elmat,1))
        for l1 = 1:nodesperelem(fes)
            Nexp[1:ndn, (l1-1)*ndn+1:(l1)*ndn] = Indn * Ns[j][l1];
        end
        push!(NexpTNexp, Nexp'*Nexp);
    end
    startassembly!(assembler,  size(elmat,1),  size(elmat,2),  count(fes), u.nfreedofs,  u.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        fill!(elmat,  0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdata, J, loc, fes.conn[i], Ns[j]);
            thefactor::FFlt =(rho*Jac*w[j]);
            elmat .+= NexpTNexp[j]*thefactor
        end # Loop over quadrature points
        gatherdofnums!(u,  dofnums,  fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler,  elmat,  dofnums,  dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function mass(self::FEMMDeforLinearAbstract,  geom::NodalField{FFlt},  u::NodalField{T}) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return mass(self, assembler, geom, u);
end

"""
    stiffness(self::FEMMDeforLinearAbstract, assembler::A,
          geom::NodalField{FFlt},
          u::NodalField{T}) where {A<:SysmatAssemblerBase, T<:Number}

Compute and assemble  stiffness matrix.
"""
function stiffness(self::FEMMDeforLinearAbstract, assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:SysmatAssemblerBase, T<:Number}
    fes = self.integdata.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdata);
    dofnums, loc, J, csmatTJ, gradN, D, B, DB, elmat, elvec, elvecfix = buffers(self, geom, u)
    self.material.tangentmoduli!(self.material, D, 0.0, 0.0, loc, 0)
    startassembly!(assembler, size(elmat, 1), size(elmat, 2), count(fes),
    u.nfreedofs, u.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        fill!(elmat,  0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdata, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], csmatTJ);
            Blmat!(self.mr, B, Ns[j], gradN, loc, self.mcsys.csmat);
            add_btdb_ut_only!(elmat, B, Jac*w[j], D, DB)
        end # Loop over quadrature points
        complete_lt!(elmat)
        gatherdofnums!(u, dofnums, fes.conn[i]); # retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums); # assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function stiffness(self::FEMMDeforLinearAbstract, geom::NodalField{FFlt},  u::NodalField{T}) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return stiffness(self, assembler, geom, u);
end

"""
    nzebcloadsstiffness(self::FEMMDeforLinearAbstract,  assembler::A,
      geom::NodalField{FFlt},
      u::NodalField{T}) where {A<:SysvecAssemblerBase, T<:Number}

Compute load vector for nonzero EBC for fixed displacement.
"""
function nzebcloadsstiffness(self::FEMMDeforLinearAbstract,  assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:SysvecAssemblerBase, T<:Number}
    fes = self.integdata.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdata);
    dofnums, loc, J, csmatTJ, gradN, D, B, DB, elmat, elvec, elvecfix = buffers(self, geom, u)
    self.material.tangentmoduli!(self.material, D, 0.0, 0.0, loc, 0)
    startassembly!(assembler,  u.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gatherfixedvalues_asvec!(u, elvecfix, fes.conn[i]);# retrieve element displacement vector
        if norm(elvecfix, Inf) != 0.0   # Is the load nonzero?
            fill!(elmat,  0.0); # Initialize element matrix
            for j = 1:npts # Loop over quadrature points
                locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
                Jac = Jacobianvolume(self.integdata, J, loc, fes.conn[i], Ns[j]);
                updatecsmat!(self.mcsys, loc, J, fes.label[i]);
                At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
                gradN!(fes, gradN, gradNparams[j], csmatTJ);
                Blmat!(self.mr, B, Ns[j], gradN, loc, self.mcsys.csmat);
                add_btdb_ut_only!(elmat, B, Jac*w[j], D, DB)
            end # Loop over quadrature points
            complete_lt!(elmat)
            gatherdofnums!(u, dofnums, fes.conn[i]); # retrieve degrees of freedom
            A_mul_B!(elvec, elmat, elvecfix)
            assemble!(assembler,  -elvec,  dofnums); # assemble element load vector
        end
    end # Loop over elements
    return makevector!(assembler);
end

function nzebcloadsstiffness(self::FEMMDeforLinearAbstract, geom::NodalField{FFlt}, u::NodalField{T}) where {T<:Number}
    assembler = SysvecAssembler()
    return  nzebcloadsstiffness(self, assembler, geom, u);
end

"""
    thermalstrainloads(self::FEMMDeforLinearAbstract, assembler::A,
        geom::NodalField{FFlt}, u::NodalField{T},
        dT::NodalField{FFlt}) where {A<:SysvecAssemblerBase, T<:Number}

Compute the thermal-strain load vector.
"""
function  thermalstrainloads(self::FEMMDeforLinearAbstract, assembler::A, geom::NodalField{FFlt}, u::NodalField{T}, dT::NodalField{FFlt}) where {A<:SysvecAssemblerBase, T<:Number}
    fes = self.integdata.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdata);
    dofnums, loc, J, csmatTJ, gradN, D, B, DB, elmat, elvec, elvecfix = buffers(self, geom, u)
    t= 0.0
    dt = 0.0
    DeltaT = fill(zero(FFlt), nodesperelem(fes))
    strain = fill(zero(FFlt), nstressstrain(self.mr)); # total strain -- buffer
    thstrain = fill(zero(FFlt), nthermstrain(self.mr)); # thermal strain -- buffer
    thstress = fill(zero(FFlt), nstressstrain(self.mr)); # thermal stress -- buffer
    startassembly!(assembler,  u.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asvec!(dT, DeltaT, fes.conn[i]);# retrieve element temperatures
        if norm(DeltaT, Inf) != 0     # Is the thermal increment nonzero?
            fill!(elvec,  0.0); # Initialize element matrix
            for j = 1:npts # Loop over quadrature points
                locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
                Jac = Jacobianvolume(self.integdata, J, loc, fes.conn[i], Ns[j]);
                updatecsmat!(self.mcsys, loc, J, fes.label[i]);
                At_mul_B!(csmatTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
                gradN!(fes, gradN, gradNparams[j], csmatTJ);#Do: gradN = gradNparams[j]/csmatTJ;
                Blmat!(self.mr, B, Ns[j], gradN, loc, self.mcsys.csmat);# strains in mater cs, displ in global cs
                self.material.thermalstrain!(self.material, thstrain, dot(vec(Ns[j]), DeltaT))
                thstress = self.material.update!(self.material, thstress, thstress, strain, thstrain, t, dt, loc, fes.label[i], :nothing)
                add_btv!(elvec, B, thstress, (-1)*(Jac*w[j]))
            end
            gatherdofnums!(u, dofnums, fes.conn[i]); # retrieve degrees of freedom
            assemble!(assembler,  elvec,  dofnums); # assemble element load vector
        end
    end # Loop over elements
    return makevector!(assembler);
end

function thermalstrainloads(self::FEMMDeforLinearAbstract, geom::NodalField{FFlt}, u::NodalField{T}, dT::NodalField{FFlt}) where {T<:Number}
    assembler = SysvecAssembler()
    return  thermalstrainloads(self, assembler, geom, u, dT);
end

"""
    inspectintegpoints(self::FEMMDeforLinearAbstract,
      geom::NodalField{FFlt},  u::NodalField{T},
      dT::NodalField{FFlt},
      felist::FIntVec,
      inspector::F,  idat, quantity=:Cauchy;
      context...) where {T<:Number, F<:Function}

Inspect integration point quantities.

`geom` - reference geometry field
`u` - displacement field
`dT` - temperature difference field
`felist` - indexes of the finite elements that are to be inspected:
     The fes to be included are: `fes[felist]`.
`context`    - structure: see the update!() method of the material.
`inspector` - functionwith the signature
        idat = inspector(idat, j, conn, x, out, loc);
   where
    `idat` - a structure or an array that the inspector may
           use to maintain some state,  for instance minimum or maximum of
           stress, `j` is the element number, `conn` is the element connectivity,
           `out` is the output of the update!() method,  `loc` is the location
           of the integration point in the *reference* configuration.
### Return
The updated inspector data is returned.
"""
function inspectintegpoints(self::FEMM, geom::NodalField{FFlt},  u::NodalField{T}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {FEMM<:FEMMDeforLinearAbstract, T<:Number, F<:Function}
    fes = self.integdata.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdata);
    dofnums, loc, J, csmatTJ, gradN, D, B, DB, elmat, elvec, elvecfix = buffers(self, geom, u)
    # Sort out  the output requirements
    outputcsys = self.mcsys; # default: report the stresses in the material coord system
    for apair in pairs(context)
        sy, val = apair
        if sy == :outputcsys
            outputcsys = val
        end
    end
    t= 0.0
    dt = 0.0
    dTe = fill(zero(FFlt), nodesperelem(fes)) # nodal temperatures -- buffer
    ue = fill(zero(FFlt), size(elmat, 1)); # array of node displacements -- buffer
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    xe = fill(zero(FFlt), nne, sdim); # array of node coordinates -- buffer
    qpdT = 0.0; # node temperature increment
    qpstrain = fill(zero(FFlt), nstressstrain(self.mr), 1); # total strain -- buffer
    qpthstrain = fill(zero(FFlt), nthermstrain(self.mr)); # thermal strain -- buffer
    qpstress = fill(zero(FFlt), nstressstrain(self.mr)); # stress -- buffer
    out1 = fill(zero(FFlt), nstressstrain(self.mr)); # stress -- buffer
    out =  fill(zero(FFlt), nstressstrain(self.mr));# output -- buffer
    # Loop over  all the elements and all the quadrature points within them
    for ilist = 1:length(felist) # Loop over elements
        i = felist[ilist];
        gathervalues_asmat!(geom, xe, fes.conn[i]);# retrieve element coords
        gathervalues_asvec!(u, ue, fes.conn[i]);# retrieve element displacements
        gathervalues_asvec!(dT, dTe, fes.conn[i]);# retrieve element temp. increments
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdata, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            At_mul_B!(csmatTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], csmatTJ);
            Blmat!(self.mr, B, Ns[j], gradN, loc, self.mcsys.csmat);
            updatecsmat!(outputcsys, loc, J, fes.label[i]);
            # Quadrature point quantities
            A_mul_B!(qpstrain, B, ue); # strain in material coordinates
            qpdT = dot(vec(dTe), vec(Ns[j]));# Quadrature point temperature increment
            self.material.thermalstrain!(self.material, qpthstrain, qpdT)
            # Material updates the state and returns the output
            out = self.material.update!(self.material, qpstress, out, vec(qpstrain), qpthstrain, t, dt, loc, fes.label[i], quantity)
            if (quantity == :Cauchy)   # Transform stress tensor,  if that is "out"
                (length(out1) >= length(out)) || (out1 = zeros(length(out)))
                rotstressvec(self.mr, out1, out, transpose(self.mcsys.csmat))# To global coord sys
                rotstressvec(self.mr, out, out1, outputcsys.csmat)# To output coord sys
            end
            # Call the inspector
            idat = inspector(idat, i, fes.conn[i], xe, out, loc);
        end # Loop over quadrature points
    end # Loop over elements
    return idat; # return the updated inspector data
end

function inspectintegpoints(self::FEMM, geom::NodalField{FFlt},  u::NodalField{T}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {FEMM<:FEMMDeforLinearAbstract, T<:Number, F<:Function}
    dT = NodalField(fill(zero(FFlt), nnodes(geom), 1)) # zero difference in temperature
    return inspectintegpoints(self, geom, u, dT, felist, inspector, idat, quantity; context...);
end

end


text:
import Pkg; Pkg.add("Documenter")
using Documenter,Juqst

makedocs(sitename="Juqst Documentation",pages = [
        "index.md",
        "chp.md",
        "open-systems.md",
        "channels.md",
        "probability-algs.md"
    ])


text:
using DashBootstrapComponents

list_group = dbc_listgroup([
    dbc_listgroupitem([
        html_div(
            [
                html_h5("This item has a heading", className = "mb-1"),
                html_small("Yay!", className = "text-success"),
            ],
            className = "d-flex w-100 justify-content-between",
        ),
        html_p("And some text underneath", className = "mb-1"),
        html_small("Plus some small print.", className = "text-muted"),
    ]),
    dbc_listgroupitem([
        html_div(
            [
                html_h5("This item also has a heading", className = "mb-1"),
                html_small("Ok!", className = "text-warning"),
            ],
            className = "d-flex w-100 justify-content-between",
        ),
        html_p("And some more text underneath too", className = "mb-1"),
        html_small("Plus even more small print.", className = "text-muted"),
    ]),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
#########################################################
#=
2-d example of the quadratic estimate.
Corresponds to Figure 4 in the paper.
Run this script with the command `include("make_figure_4.jl")` entered into the Julia REPL in this directory.
This script saves figure4.pdf to the current directory.
To include figure4.pdf in the paper it must be copied to the directory paper/
=#
#########################################################

# --- load modules
using NonstationaryPhase
using PyPlot

# ----  set the seed
#seedstart = rand(UInt64)
seedstart = 0xa65c9cf99fcbf681
srand(seedstart)

# --- save figures to disk or not
savefigures = true

# --- set grid geometry
const d_const = 2
const period_const = 2π
const nside_const  = nextprod([2,3,5,7], 400)

# --- parameters of Z(x)
ν      = 1.5
tild_ν = 1.7
ρ      = 0.015*period_const
tild_ρ = 0.014*period_const
t_0    = 1.5

# --- noise parameters
σpixl    = 0.0
beamFWHM = 0.0

# --- prior parameters
νΦ = 5.0
ρΦ = 0.15*period_const
σΦ = 0.5*(period_const/2π)^2

#########################################################
#=
Generate the first set of simulations
=#
#########################################################

# --- define an instance of PhaseParms which carries all the parameters
parms = PhaseParms(d_const, period_const, nside_const)
ηkfun, Ckfun, C1kfun, C2kfun = gen_tangentMatern(ν=ν, tild_ν=tild_ν, ρ=ρ, tild_ρ=tild_ρ, t_0=t_0)
NonstationaryPhase.CNNkfun{dm,T1,T2}(p::PhaseParms{dm,T1,T2})   = CNNkfun(p.k, p.deltx; σpixl=σpixl, beamFWHM=beamFWHM)
Cϕϕkfun{dm,T1,T2}(p::PhaseParms{dm,T1,T2})  = maternk(p.k; ν=νΦ, ρ=ρΦ, σ=σΦ)
parms.ηk[:]       =  ηkfun(parms.k)
parms.Ck[:]       =  Ckfun(parms.k)
parms.C1k[:]      =  C1kfun(parms.k)
parms.C2k[:]      =  C2kfun(parms.k)
parms.CNNk[:]     =  CNNkfun(parms)
parms.Cϕϕk[:]     = Cϕϕkfun(parms)
parms.CZZmk[:]    =  CZZmkfun(parms.Cϕϕk, parms.Ck, parms.ηk, parms.ξk, parms.x, parms.k, parms.deltx, parms.deltk)
parms.CZZmobsk[:] = parms.CZZmk + parms.CNNk
parms.CZZmobsk[parms.r .< 1]              = Inf
parms.CZZmobsk[parms.r .>= 0.9*parms.nyq] = Inf
parms.ξk[:] = Array{Complex{Float64},2}[ im.*parms.k[2], -im.*parms.k[1] ]  # new

# Now a High resolution version of PhaseParms for simulation
parmsHR = PhaseParms(d_const, period_const, nside_const)
parmsHR.ηk[:]   =  ηkfun(parmsHR.k)
parmsHR.Ck[:]   =  Ckfun(parmsHR.k)
parmsHR.CNNk[:] =  CNNkfun(parmsHR)
parmsHR.Cϕϕk[:] =  Cϕϕkfun(parmsHR)
parms.ξk[:] = Array{Complex{Float64},2}[ im.*parmsHR.k[2], -im.*parmsHR.k[1] ]  # new

# ---- generate the  simulation and other plotting quantities
zkobs, zk, zx, zx_noϕ, ϕk, ϕx = simNPhaseGRF(parms, parmsHR)
Dθx          = NonstationaryPhase.ϕx_2_Dθx(ϕx, parms)
@time Cℓvar        = Cℓvarfun(parms)
@time approxCℓbias = approxCℓbiasfun(parms)
@time estϕk  = estϕkfun(zkobs, parms)
estϕx  = ifftdr(estϕk, parms.deltk)

# #########################################################
# #=
# Now plot
# =#
# #########################################################

# ------ figure ---- show the quadratic estimate
function makeplots()
    global fig1 = figure(figsize = (10,10))

    # ----- figure: plot Z(x)
    ax1 = subplot2grid((2,2), (1,1))
    imshow(zx,
        interpolation = "nearest",
        origin="lower",
        extent=[minimum(parms.x[1]), maximum(parms.x[1]),minimum(parms.x[2]), maximum(parms.x[2])],
    )
    title("Nonstationary phase simulation")

    # ----- figure: plot Cℓvarfun, Cℓbias and the signal
    ax2 = subplot2grid((2,2), (1,0))
    δ0 = 1 / parms.deltk ^ d_const
    kbins, MSE_radpower = radial_power(abs2(parms.r.*(estϕk-ϕk))./δ0, 1, parms)
    semilogy(parms.r[1:25,1], abs2(parms.r[1:25,1]).*parms.Cϕϕk[1:25,1], "k:", label = L"\ell^2 C^{\phi\phi}_\ell", linewidth = 3.0)
    semilogy(parms.r[1:25,1], abs2(parms.r[1:25,1]).*approxCℓbias[1:25,1], "b--", label = L"approx $\ell^2 C^{bias\, \hat\phi}_\ell$", linewidth = 2.0)
    semilogy(parms.r[1:25,1], abs2(parms.r[1:25,1]).*Cℓvar[1:25,1], "g-", label = L"\ell^2C^{var\, \hat\phi}_\ell", linewidth = 2.0)
    semilogy(kbins[1:25],MSE_radpower[1:25], "go", label = L"radial $\ell^2 |\phi_\ell - \hat\phi_\ell|^2\delta_0$", linewidth = 2.0)
    #ylabel("spectral density")
    #xlabel(L"frequency $\ell$")
    title("Variance, bias and empirical MSE")
    #axis("tight")
    legend(loc="best",fontsize = 10)

    # --- figure: plot est Φ(x)
    ax4 = subplot2grid((2,2), (0,0))
    imshow(estϕx,
        interpolation = "nearest",
        origin="lower",
        extent=[minimum(parms.x[1]), maximum(parms.x[1]),minimum(parms.x[2]), maximum(parms.x[2])],
        )
    title(L"Quadratic estimate $\hat\phi(x)$")

    # --- figure: plot true Φ(x)
    ax3 = subplot2grid((2,2), (0,1))
    imshow(ϕx,
        interpolation = "nearest",
        origin="lower",
        extent=[minimum(parms.x[1]), maximum(parms.x[1]),minimum(parms.x[2]), maximum(parms.x[2])],
        vmin = minimum(estϕx), vmax = maximum(estϕx)
        )
    title(L"Simulation truth $\phi(x)$")

    fig1[:subplots_adjust](wspace=0.1)

end

# --- now evaluate the above `script'
makeplots()

# --- save the figures to disk
if savefigures
    fig1[:savefig]("figure4.pdf", dpi=300, bbox_inches="tight", transparent=true)
    plt[:close](fig1)
end


text:
using SimulationLogger
using DifferentialEquations
using Transducers
using Plots
using Test

function test()
    @Loggable function dynamics!(dx, x, p, t; u)
        @onlylog state, input = x, u  # __LOGGER_DICT__[:state] = x, __LOGGER_DICT__[:input] = u
        dx .= u
    end
    @Loggable function custom_control(x)
        @log a = 1
        -a*x
    end
    @Loggable function feedback_dynamics!(dx, x, p, t)
        @onlylog time = t  # __LOGGER_DICT__[:time] = t
        @log x, t  # __LOGGER_DICT__[:x] = x
        u = @nested_log custom_control(x)  # __LOGGER_DICT__[:a] = 1
        @log u  # __LOGGER_DICT__[:u] = -a*x
        @nested_log :linear x
        @nested_log :linear dynamics!(dx, x, p, t; u=u)
    end
    t0, tf = 0.0, 0.1
    Δt = 0.01
    saved_values = SavedValues(Float64, Dict)
    cb = CallbackSet()
    if hasmethod(feedback_dynamics!, Tuple{Any, Any, Any, Any, __LOG_INDICATOR__})
        # to avoid undefined error when not adding @Loggable
        log_func(x, t, integrator::DiffEqBase.DEIntegrator; kwargs...) = feedback_dynamics!(zero.(x), copy(x), integrator.p, t, __LOG_INDICATOR__(); kwargs...)
        cb = SavingCallback(log_func, saved_values; saveat=t0:Δt:tf)
    end
    # # sim
    x0 = [1, 2, 3]
    tspan = (t0, tf)
    prob = ODEProblem(
                      feedback_dynamics!, x0, tspan;
                      callback=cb,
                     )
    _ = solve(prob)
    ts = saved_values.saveval |> Map(datum -> datum[:t]) |> collect
    xs = saved_values.saveval |> Map(datum -> datum[:x]) |> collect
    us = saved_values.saveval |> Map(datum -> datum[:u]) |> collect
    times = saved_values.saveval |> Map(datum -> datum[:time]) |> collect
    states = saved_values.saveval |> Map(datum -> datum[:linear][:state]) |> collect
    inputs = saved_values.saveval |> Map(datum -> datum[:linear][:input]) |> collect
    as = saved_values.saveval |> Map(datum -> datum[:a]) |> collect
    @test ts == saved_values.t
    @test ts == times
    @test xs == states
    @test us == inputs
    @test as == ones(length(ts))
    p_x = plot(ts, hcat(xs...)')
    p_u = plot(ts, hcat(us...)')
    dir_log = "figures"
    mkpath(dir_log)
    savefig(p_x, joinpath(dir_log, "state.png"))
    savefig(p_u, joinpath(dir_log, "input.png"))
end


text:
using SnpArrays
using CSV, DataFrames
#
# separate by population
#
populations = ["african", "asian", "bangladeshi", "british", "caribbean", "chinese",
    "indian", "irish", "pakistani", "white_asian", "white_black", "white"]
xdata = SnpData("ukb_gen_merged") # 407434×591513 (all sample populations with all snps passing QC)
p = size(xdata.snparray, 2)
sampleIDs = parse.(Int, xdata.person_info.iid)
for pop in populations
    ids = CSV.read("samples_$pop.txt", DataFrame)[:, 1]
    sample_idx = sort!(Vector{Int}(indexin(ids, sampleIDs)))
    SnpArrays.filter("ukb_gen_merged", sample_idx, 1:p; des = "ukb.$pop")
end

#
# Subset chromosome 10 for testing
#
chr = 10
chr_snp_idx = findall(x -> x == "10", xdata.snp_info.chromosome)
mkdir("chr$chr")
cd("chr$chr")
for pop in populations
    sample_idx = 1:size(SnpArray("../ukb.$pop.bed"), 1)
    SnpArrays.filter("../ukb.$pop", sample_idx, chr_snp_idx; des = "ukb.chr$chr.$pop")
end


text:
using Gaugefields
using Wilsonloop
using LinearAlgebra

function make_cloverloops(μ,ν;Dim=4)
    loops = Wilsonline{Dim}[]
    loop_righttop = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)])
    loop_lefttop = Wilsonline([(ν,1),(μ,-1),(ν,-1),(μ,1)])
    loop_rightbottom = Wilsonline([(ν,-1),(μ,1),(ν,1),(μ,-1)])
    loop_leftbottom= Wilsonline([(μ,-1),(ν,-1),(μ,1),(ν,1)])
    push!(loops,loop_righttop)
    push!(loops,loop_lefttop)
    push!(loops,loop_rightbottom)
    push!(loops,loop_leftbottom)
    return loops
end

function cloverloops(Dim)
    loops_μν= Matrix{Vector{Wilsonline{Dim}}}(undef,Dim,Dim)
    for μ=1:Dim
        for ν=1:Dim
            loops_μν[μ,ν] = Wilsonline{Dim}[]
            if ν == μ
                continue
            end
            loops_μν[μ,ν] = make_cloverloops(μ,ν,Dim=Dim)
        end
    end
    return  loops_μν
end

function make_energy_density!(Wmat,U::Vector{<: AbstractGaugefields{NC,Dim}},temps) where {NC,Dim}
    W_operator = cloverloops(Dim)
    calc_wilson_loop!(Wmat,W_operator,U,temps)
    return
end

function calc_wilson_loop!(W,W_operator,U::Vector{<: AbstractGaugefields{NC,Dim}},temps) where {NC,Dim}
    for μ=1:Dim
        for ν=1:Dim
            if μ == ν
                continue
            end
            evaluate_gaugelinks!(W[μ,ν],W_operator[μ,ν],U,temps)
            W[μ,ν] = Traceless_antihermitian(W[μ,ν])
        end
    end
    return
end

function  make_energy_density_core(Wmat::Matrix{<: AbstractGaugefields{NC,Dim}}) where {NC,Dim}
    @assert Dim == 4
    W = 0.0 + 0.0im
    for μ=1:Dim # all directions
        for ν=1:Dim
            if μ == ν
                continue
            end
            W += -tr(Wmat[μ,ν],Wmat[μ,ν])/2
        end
    end
    return W
end

function calculate_energy_density(U::Array{T,1}, Wmat,temps) where T <: AbstractGaugefields
    # Making a ( Ls × Lt) Wilson loop operator for potential calculations
    WL = 0.0+0.0im
    NV = U[1].NV
    NC = U[1].NC
    make_energy_density!(Wmat,U,temps) # make wilon loop operator and evaluate as a field, not traced.
    WL =  make_energy_density_core(Wmat) # tracing over color and average over spacetime and x,y,z.
    NDir = 4.0*3.0/2 # choice of 2 axis from 4.
    return real(WL)/(NV*4^2)
end

function test()
    NX = 8
    NY = 8
    NZ = 8
    NT = 8
    Nwing = 0
    NC = 3

    Dim = 4

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot")
    filename="./conf_08080808.ildg"

    ildg = ILDG(filename)
    i = 1
    L = [NX,NY,NZ,NT]
    load_gaugefield!(U,i,ildg,L,NC)

    temp1 = similar(U[1])
    temp2 = similar(U[1])
    temp3 = similar(U[1])

    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    β = 5.7
    W_temp = Matrix{typeof(U[1])}(undef,Dim,Dim)
    for μ=1:Dim
        for ν=1:Dim
            W_temp[μ,ν] = similar(U[1])
        end
    end

    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println(" plaq_t = $plaq_t")

    dt = 0.01

    g = Gradientflow(U,eps = dt)

    for itrj=1:10
        flow!(U,g)
        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        e = calculate_energy_density(U, W_temp,[temp1,temp2])
        println("$itrj $dt $plaq_t $e # itrj dt plaq energy")
    end

end
test()


text:
using Unitful
function unitful_testfunction(Vi)
    if Vi ≤ 0.0u"V"
        return 0.0u"V"
    elseif Vi ≥ 1.0u"V"
        return 1.0u"V"
    else
        return Vi
    end
end

register_primitive(unitful_testfunction) # must be outside testset

@testset "Unitful" begin
    @info "Testing Unitful"

    PT = Particles
    for PT in (Particles, StaticParticles)
        p1 = PT(100, Uniform(-0.5,1.5)) * 1u"V"
        p2 = PT(100, Uniform(-0.5,1.5)) * u"V"
        @test_nowarn println(p1)
        @test_nowarn display(p1)

        @test_nowarn println(0p1)
        @test_nowarn display(0p1)

        @test typeof(p1) == typeof(p2)

        p3 = unitful_testfunction(p1)
        @test pextrema(p3) == (0.0u"V", 1.0u"V")

        @test (1 ± 0.5)u"m" * (1 ± 0)u"kg" ≈ (1 ± 0.5)u"kg*m"
        @test (1 ± 0.5)u"m" * 1u"kg" ≈ (1 ± 0.5)u"kg*m"

        @test (1 ± 0.5)u"m" / (1 ± 0)u"kg" ≈ (1 ± 0.5)u"m/kg"
        @test (1 ± 0.5)u"m" / 1u"kg" ≈ (1 ± 0.5)u"m/kg"

        @test (1 ± 0.5)u"m" + (1 ± 0)u"m" ≈ (2 ± 0.5)u"m"
        @test (1 ± 0.5)u"m" + 1u"m" ≈ (2 ± 0.5)u"m"

        @test 1u"m" * (1 ± 0.5)u"kg" ≈ (1 ± 0.5)u"kg*m"
        @test 1u"m" / (1 ± 0.5)u"kg" ≈ (1 ± 0.5)u"m/kg"
        @test 1u"m" + (1 ± 0.5)u"m" ≈ (2 ± 0.5)u"m"

        typeof(promote(1u"V", (1.0 ± 0.1)u"V")) <: Tuple{Particles{<:Quantity}, Particles{<:Quantity}}

        ρ = (2.7 ± 0.2)u"g/cm^3"
        mass = (250 ± 10)u"g"
        width = (30.5 ± 0.2)u"cm"
        l = (14.24 ± 0.2)u"m"
        thickness = u"μm"(mass/(ρ*width*l))
        @test thickness ≈ (21.3 ± 1.8)u"μm"

    end

end


text:
using Distributions

export mutate!,
    MutationOperator,
    setup!

export BitFlipMutation,
    AdditiveMutation,
    LiftedMutation,
    PointwiseMutation

abstract type MutationOperator{G} end

setup!(operator::MutationOperator{G}, model::AbstractPopulationModel{G}) where {G} = operator

"""
    mutate!(genome, operator, generation) -> genome
    mutate!(individual, operator, generation) -> individual

Mutate `genome` or `individual.genome` in place.  You only need to define the needed genome form.
"""
function mutate!(individual::Individual{G}, operator::MutationOperator{G},
                 generation::Integer) where {G}
    mutate!(individual.genome, operator, generation)
    individual
end

struct NoMutation <: MutationOperator{Any} end

mutate!(genome::Any, operator::NoMutation) = genome

struct LiftedMutation{G, M, I} <: MutationOperator{G}
    inner::M

    LiftedMutation{G}(operator::M) where {G, I, M<:MutationOperator{I}} = new{G, M, I}(operator)
    LiftedMutation{G, M}(args...) where {G, I, M<:MutationOperator{I}} = new{G, M, I}(M(args...))
end

struct BitFlipMutation <: MutationOperator{AbstractVector{Bool}}
    rate::Rate
end

function mutate!(genome::AbstractVector{Bool}, operator::BitFlipMutation, generation::Integer)
    for i in eachindex(genome)
        (rand() < operator.rate(generation)) && (genome[i] = !genome[i])
    end

    genome
end

struct PointwiseMutation{G<:AbstractVector} <: MutationOperator{G}
    rate::Rate
    tweak::Distribution{Univariate}

    PointwiseMutation{G}(rate, tweak::Distribution{Univariate, Discrete}) where {G<:AbstractVector{<:Integer}} =
        new{G}(rate, tweak)
    PointwiseMutation{G}(rate, tweak::Distribution{Univariate, Continuous}) where {G<:AbstractVector{<:AbstractFloat}} =
        new{G}(rate, tweak)
end

function mutate!(genome::G, operator::PointwiseMutation{G},
                 generation::Integer) where {G<:AbstractVector}
    for i in eachindex(genome)
        (rand() < operator.rate(generation)) && (genome[i] = rand(operator.tweak))
    end

    genome
end

struct AdditiveMutation{G<:AbstractVector, D<:Distribution{Univariate}} <: MutationOperator{G}
    rate::Rate
    tweak::D
    min::Real
    max::Real

    function AdditiveMutation{G}(rate, tweak::D, min, max) where {G, D<:Distribution{Univariate}}
        @assert (mean(tweak) == 0) "`tweak` should have zero mean!"
        new{G, D}(rate, tweak, min, max)
    end

    AdditiveMutation{G, D}(rate, r, min, max) where {G, D<:Uniform} =
        new{G, D}(rate, D(-r, r), min, max)
    AdditiveMutation{G, D}(rate, σ, min, max) where {G, D<:Normal} =
        new{G, D}(rate, D(0, σ), min, max)
    AdditiveMutation{G, D}(rate, r, min, max) where {G, D<:DiscreteUniform} =
        new{G, D}(rate, D(-r, r), min, max)
end

function mutate!(genome::G, operator::AdditiveMutation{G},
                 generation::Integer) where {G<:AbstractVector}
    for i in eachindex(genome)
        if rand() ≤ operator.rate(generation)
            genome[i] = rand(Truncated(Shifted(operator.tweak, genome[i]),
                                       operator.min, operator.max))
        end
    end

    genome
end


text:
using CuArrays

@testset "CUDA" begin
    @info "Testing CUDA"
    for f ∈ readdir("./cuda/")
        include("./cuda/$f")
    end
end


text:
using Distributed

# ----------------------------------------
# CONFIGURATION
# ----------------------------------------

module CommandLine

    # ----------------------------------------
    # IMPORTED NAMESPACES
    # ----------------------------------------

    import ArgParse

    # ----------------------------------------
    # CONFIGURATION
    # ----------------------------------------

    export arguments

    # ----------------------------------------
    # EXPORTED FUNCTIONS
    # ----------------------------------------

    function arguments(arguments)
        parser = ArgParse.ArgParseSettings()

        @ArgParse.add_arg_table parser begin
            "--workers"
            arg_type = Int
            default  = 0
            help     = "number of worker processes"
        end

        return ArgParse.parse_args(arguments, parser)
    end
end

# ----------------------------------------
# DISTRIBUTED SET UP
# ----------------------------------------

arguments = CommandLine.arguments(ARGS)
worker_count = arguments["workers"]

@info "WORKERS ... $(worker_count)"

if worker_count > 0
    addprocs(worker_count)
end

# ----------------------------------------
# IMPORTS
# ----------------------------------------

@everywhere using SharedArrays

# ----------------------------------------
# FUNCTIONS
# ----------------------------------------

@everywhere function transfer(from::SharedArray{Int64}, to::SharedArray{Int64}, channel)
    function work(i, j)
        return i * j
    end

    for i ∈ localindices(from)
        snooze = rand()
        to[i] = work(from[i], myid())
        @info "IN $(myid()) ... $(i) ... $(snooze)."
        sleep(snooze)
        put!(channel, true)
    end

    return 0
end

# ----------------------------------------
# MAIN PROCESSING
# ----------------------------------------

const jobs_count = 16

S = SharedArray{Int64, 1}(jobs_count, )
T = SharedArray{Int64, 1}(jobs_count, )

channel = RemoteChannel(() -> Channel{Bool}(jobs_count))

for i ∈ 1:length(S)
    S[i] = i
end

@sync begin
    for worker in procs(S)
        @async begin
            remotecall_wait(transfer, worker, S, T, channel)
        end
    end

    for i ∈ 1:jobs_count
        v = take!(channel)
        @info "HERE ... $(i)/$(jobs_count)"
    end

    @info "DONE!!!"
end

@show S
@show T

max_stars_repo_name: jlegare/Parallel

max_issues_repo_name: jlegare/Parallel

max_forks_repo_name: jlegare/Parallel


text:
include("config.jl")
println("evalon=$evalon")
include("gtest1.jl")
abc()

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
__precompile__(false)
module DiffDRO

using Conda
using PyCall

export DiffBasicWeightConvexSet,
    DiffBertsimas,
    DiffBenTal,
    DiffDelage,
    utility,
    expected_return,
    std_adjusted_expected_return,
    var_adjusted_expected_return

Conda.add("pip")

run(PyCall.python_cmd(`-m pip install numpy`))
run(PyCall.python_cmd(`-m pip install scs==2.1.3`))
run(PyCall.python_cmd(`-m pip install cvxpy`))
run(PyCall.python_cmd(`-m pip install torch`))
run(PyCall.python_cmd(`-m pip install cvxpylayers`))

py"""
import numpy as np
import cvxpy as cp
import torch
from cvxpylayers.torch import CvxpyLayer
"""

include("DiffBasicWeightConvexSet.jl")
include("DiffBertsimas.jl")
include("DiffBenTal.jl")
include("DiffDelage.jl")
include("DiffObj.jl")

include("test_utils.jl")

end

max_stars_repo_name: molet/DiffDRO

max_issues_repo_name: molet/DiffDRO

max_forks_repo_name: molet/DiffDRO


text:
#########################################################################
#Forecast using random walk model, which is a special case of AR1 model
#with AR parameter set to be 1
#########################################################################
function rwm_oos_forc(J, h, h_lag, Y, Y_lag)
    ŷ = zeros(Float64, size(J - (h_lag + h):length(Y)))
    for (n, vec_î) in enumerate(J - (h_lag + h):length(Y))
        LHS = Y[1:(vec_î - h)]
        sample_size = length(LHS)
        LHS_lag = Y_lag[1:sample_size]
        ε̂ = mean(LHS - LHS_lag)
        ŷ[n] = Y_lag[vec_î] + ε̂
    end
    RMSE = sqrt(mean((Y[J - (h_lag + h):length(Y)] - ŷ) .^ 2))
    return ŷ, RMSE
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning