text:
# # Generating Data
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/gendata.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/gendata.ipynb)

# RHEOS has several convenience functions for generating arbitrarily complex loading patterns. These may be particularly useful for investigating the responses of viscoelastic models that users may be unfamiliar with.

using RHEOS
using PyPlot

# #### Step 1: Generate Timeline
# The first step requires the generation of a [`RheoTimeData`](@ref) struct with only the time data. This is achieved with the [`timeline`](@ref) function. If no parameter is specified, it generates a time data from `t_start = 0` to `t_end = 10`, with a default time step of `(t_start-t_end)/250 = 0.04`.

datat = timeline(t_start = 0, t_end = 20.0, step = 0.02)
rheotimedatatype(datat)

#md # !!! note "Note"
#md #     The same functionality is provided to generate a frequency spectrum by using [`frequencyspec`](@ref) (see an example of its usage in the [Fitting and Predicting - Frequency data](@ref) section).

#nb # Note that the same functionality is provided to generate a frequency spectrum by using frequencyspec (see an example of its usage in the 'Fitting and Predicting - Frequency data' section).

# #### Step 2: Generate Stress or Strain Data

# The user can decide either to generate strain data ([`strainfunction`](@ref)) or stress data ([`stressfunction`](@ref)). Both functions require a [`RheoTimeData`](@ref) with time data defined and a function that takes the time as an input to generate the desired set of data. Both function return a new [`RheoTimeData`](@ref) struct with time data and either stress or strain data. Below is an example.

## calculates strain data by applying a function of time
dsin = strainfunction(datat, t -> sin(2*t))

fig, ax = subplots(1, 1, figsize = (8, 2))
plot(dsin.t, dsin.ϵ)
rheotimedatatype(dsin)
#!nb fig #hide

# ## Pre-defined Functions

# RHEOS provides some functions that can be directly applied to a timeline. Below is the complete list with examples. Note that if the keyword arguments for the functions are not defined, default parameters are used (see the [API](@ref) section). In the examples below the strain data is generated, the same can be applied to the stress by replacing [`strainfunction`](@ref) with [`stressfunction`](@ref).

fig, ax = subplots(2, 3, figsize = (14, 7))

## Step generation
dstep = strainfunction(datat, hstep(offset = 5.0, amp = 5))
    ## Plotting
    ax[1, 1].plot(dstep.t, dstep.ϵ)
    ax[1, 1].set_ylim(bottom = -0.1, top = 5.4)
    ax[1, 1].set_title("Step")

## Ramp generation
dramp = strainfunction(datat, ramp(offset = 2.0, gradient = 0.8))
    ax[1, 2].plot(dramp.t, dramp.ϵ)
    ax[1, 2].set_ylim(bottom = -0.1)
    ax[1, 2].set_title("Ramp")

## Stairs generation
dstairs = strainfunction(datat, stairs(offset = 5.0, amp = 0.1, width = 1))
    ax[1, 3].plot(dstairs.t, dstairs.ϵ)
    ax[1, 3].set_ylim(bottom = -0.1)
    ax[1, 3].set_title("Stairs")

## Square generation
dsquare = strainfunction(datat, square(offset = 5.0, amp = 0.5, period = 4, width = 1))
    ax[2, 1].plot(dsquare.t, dsquare.ϵ)
    ax[2, 1].set_ylim(bottom = -0.05, top = 0.6)
    ax[2, 1].set_title("Square")

## Sawtooth generation
dsawtooth = strainfunction(datat, sawtooth(offset = 5.0, amp = 2, period = 5))
    ax[2, 2].plot(dsawtooth.t, dsawtooth.ϵ)
    ax[2, 2].set_ylim(bottom = -0.1, top = 2.2)
    ax[2, 2].set_title("Sawtooth")

## Triangle generation
dtriangle = strainfunction(datat, triangle(offset = 4.0, amp = 1, period = 4))
    ax[2, 3].plot(dtriangle.t, dtriangle.ϵ)
    ax[2, 3].set_ylim(bottom = -0.1, top = 1.2)
    ax[2, 3].set_title("Tringle")

#!nb fig #hide

# ## More Complex Patterns

# More complicated patterns can be generated by making use of the possibility to add, subtract, multiply [RheoTimeData](@ref) structs. Note that addition and subtraction can be performed only if the data have the same sample rate.

fig, ax = subplots(1, 3,figsize = (14, 4))

## Ramp & hold
dhold = dramp - strainfunction(datat, ramp(offset = 5.0, gradient = 0.8))
    ax[2].plot(dhold.t, dhold.ϵ)
    ax[2].set_ylim(bottom = -0.1, top = 2.8)
    ax[2].set_title("Ramp & hold")

## Step with oscillatory loading
doscil = dstep + dsquare
    ax[1].plot(doscil.t, doscil.ϵ)
    ax[1].set_ylim(bottom = -0.5, top = 6)
ax[1].set_title("Oscillations")

## Ramp & hold & oscillation
dcomplex = dhold + dsquare
    ax[3].plot(dcomplex.t, dcomplex.ϵ)
    ax[3].set_ylim(bottom = -0.1, top = 3.2)
    ax[3].set_title("Ramp & hold & oscillations")

#!nb fig #hide


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e9146374-dd10-11ea-1b83-c978edacc525
begin
	using Plots
	using Random
	using Distributions
	using Printf

	pyplot();
	Random.seed!(20200813);
end

# ╔═╡ 7396bcf4-dd11-11ea-2577-55daf9dbbd68
begin
	N = 15
	noise_scale = 5
	x = collect(range(1, 10, length=N))
	ϵ = noise_scale*randn(N)
	y = x .+ x.^2 .+ ϵ
end

# ╔═╡ d32c85c0-dd13-11ea-0476-8bec5f639c3a
begin
	# loss function
	loss(w) = mean((w.*x - y).^2)
	gradient(w) = mean(2*(w.*x - y).*x)
end

# ╔═╡ 3fdce730-dd39-11ea-216b-efd542f9528e
@bind w html"<input type='range' min='-10' max='30' value='$w' step='0.1'>"

# ╔═╡ e0649a0e-dd11-11ea-1e21-19d9fb8fe64b
begin
	η = 0.01
	current_loss = loss(w)
	current_loss_str = @sprintf "%.1f" current_loss
	w_str = @sprintf "%.2f" w
	grad = gradient(w)

	p1 = scatter(x, y, label="")
	ŷ = w .* x
	plot!(x, ŷ, label="ŷ")

	p2 = plot(loss, -10:25,  xlims=(-10, 30), ylims = (0,10000), label="loss(w)")
	annotate!((w, loss(w)*1.1, Plots.text("w=$w_str, loss=$current_loss_str", 10) ))
	vline!([w], label="")
	plot!([w, w-η*grad], fill(current_loss, 2), label="gradient", color=:red, arrow=:right)

	plot(p1, p2, layout=(1, 2))

end

# ╔═╡ Cell order:
# ╠═e9146374-dd10-11ea-1b83-c978edacc525
# ╠═7396bcf4-dd11-11ea-2577-55daf9dbbd68
# ╠═d32c85c0-dd13-11ea-0476-8bec5f639c3a
# ╠═3fdce730-dd39-11ea-216b-efd542f9528e
# ╟─e0649a0e-dd11-11ea-1e21-19d9fb8fe64b

ext: jl

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
ENV["MPLBACKEND"]="agg" # no GUI
using PyPlot, PyCall
@info("SimpleHypergraphs is using Matplotlib $(PyPlot.version) with Python $(PyCall.pyversion)")

using Test, SimpleHypergraphs, StatsBase
using Random
using DataStructures
import LightGraphs

h1 = Hypergraph{Float64, Int, String}(5,4)
h1[1:3,1] .= 1.5
h1[3,4] = 2.5
h1[2,3] = 3.5
h1[4,3:4] .= 4.5
h1[5,4] = 5.5
h1[5,2] = 6.5

@testset "SimpleHypergraphs Hypergraph          " begin

    h = hg_load("data/test1.hgf"; T=Int)
    @test size(h) == (4, 4)
    @test nhv(h) == 4
    @test nhe(h) == 4
    m = Matrix(h)
    @test m == h
    @test h == [1       nothing 4       nothing
                2       3       nothing nothing
                nothing nothing 5       nothing
                nothing nothing 6       nothing]
    mktemp("data") do path, _
        println(path)
        hg_save(path, h)

        loaded_hg = replace(read(path, String), r"\n*$" => "")

        @test loaded_hg ==
            reduce(replace,
                ["\r\n"=>"\n",
                r"^\"\"\"(?s).*\"\"\"\n"=>"", #remove initial comments
                r"\n*$"=>""], #remove final \n*
                init=read("data/test1.hgf", String)) #no comments

        @test loaded_hg ==
            reduce(replace,
                ["\r\n"=>"\n",
                r"^\"\"\"(?s).*\"\"\"\n"=>"", #remove initial comments
                r"\n*$"=>""], #remove final \n*
                init=read("data/test_singlelinecomment.hgf", String)) #single line comment

        @test loaded_hg ==
            reduce(replace,
                ["\r\n"=>"\n",
                r"^\"\"\"(?s).*\"\"\"\n"=>"", #remove initial comments
                r"\n*$"=>""], #remove final \n*
                init=read("data/test_multiplelinescomment.hgf", String)) #multiple lines comment

        for v=1:nhv(h)
            set_vertex_meta!(h1, v, v)
        end

        for he=1:nhe(h)
            set_hyperedge_meta!(h1, string(he), he)
        end

        hg_save(path, h1; format=JSON_Format())
        loaded_hg = hg_load(path; format=JSON_Format(), T=Float64, V=Int, E=String)

        @test h1 == loaded_hg
        @test h1.v_meta == loaded_hg.v_meta
        @test h1.he_meta == loaded_hg.he_meta

        @test get_vertex_meta(h1, 1) == get_vertex_meta(loaded_hg, 1)
        @test get_hyperedge_meta(h1, 2) == get_hyperedge_meta(loaded_hg, 2)

    end

    @test_throws ArgumentError hg_load("data/test_malformedcomment.hgf"; T=Int)
    @test_throws ArgumentError hg_load("data/test_argumenterror.hgf"; T=Int)

    h2 = Hypergraph{Float64}(0,0)
    @test h2 == Hypergraph{Float64,Nothing}(0,0)
    @test h2 == Hypergraph{Float64,Nothing,Nothing}(0,0)
    @test h2 == Hypergraph{Float64,Nothing,Nothing,Dict{Int,Float64}}(0,0)

    h3 = Hypergraph(0,0)
    @test h3 == Hypergraph{Bool, Nothing, Nothing, Dict{Int, Bool}}(0,0)

    for i in 1:4 add_vertex!(h2) end
    add_hyperedge!(h2;vertices=Dict(1:3 .=> 1.5))
    add_hyperedge!(h2)
    add_vertex!(h2;hyperedges=Dict(2=>6.5))
    add_hyperedge!(h2;vertices=Dict(2 => 3.5, 4 => 4.5))
    add_hyperedge!(h2;vertices=Dict(3:5 .=> (2.5,4.5,5.5)))
    @test h1 == h2
    m = Matrix(h1)
    @test  m == Matrix(h2)
    @test h1 == Hypergraph(m)
    @test h1 == Hypergraph{Float64}(m)
    @test h1 == Hypergraph{Float64,Nothing}(m)
    @test h1 == Hypergraph{Float64,Nothing, Nothing}(m)
    @test h1 == Hypergraph{Float64,Nothing, Nothing,Dict{Int,Float64}}(m)
    @test all(Matrix(h1) .== Matrix(Hypergraph{Float64,Nothing, Nothing,SortedDict{Int,Float64}}(m)))
    @test getindex(h1,3,1) == 1.5

    h3 = Hypergraph{Float64,String,Nothing}(1,1)
    @test add_vertex!(h3;v_meta="test") == 2
    @test set_vertex_meta!(h3,"t",1) == ["t","test"]
    @test get_vertex_meta(h3,2) == "test"
    @test get_hyperedge_meta(h3,1) == nothing
    @test_throws BoundsError get_hyperedge_meta(h3,2)

    h4 = Hypergraph{Float64,Nothing,String}(1,1)
    @test add_hyperedge!(h4;he_meta="test") == 2
    @test set_hyperedge_meta!(h4,"t",1) == ["t","test"]
    @test get_hyperedge_meta(h4,2) == "test"
    @test get_vertex_meta(h4,1) == nothing
    @test_throws BoundsError get_vertex_meta(h4,2)

    h5 = Hypergraph{Float64,String,String,SortedDict{Int,Float64}}(1,1)
    @test typeof(h5.v2he[1]) <: SortedDict{Int,Float64}
    @test typeof(h5.he2v[1]) <: SortedDict{Int,Float64}
    @test add_vertex!(h5;v_meta="test") == 2
    @test set_vertex_meta!(h5,"t",1) == ["t","test"]
    @test get_vertex_meta(h5,2) == "test"
    @test get_hyperedge_meta(h5,1) == nothing
    @test add_hyperedge!(h5;he_meta="test") == 2
    @test set_hyperedge_meta!(h5,"t",1) == ["t","test"]
    @test get_hyperedge_meta(h5,2) == "test"
    @test_throws BoundsError get_vertex_meta(h5,3)
    @test_throws BoundsError get_hyperedge_meta(h5,3)
    h5 .= [1.0 2.0;3.0 4.0]
    @test h5[2,2] == 4

    h1_0 = deepcopy(h1)
    @test add_vertex!(h1_0) == 6
    h1_0[6,:] = h1_0[5,:]
    @test remove_vertex!(h1_0,5) == h1
    setindex!(h1_0, nothing, 1, 1)
    @test h1_0[1,1] == nothing
    @test_throws BoundsError setindex!(h1_0, nothing, 10, 9)

    h1_1 = Hypergraph([nothing nothing nothing nothing
                       1       1       nothing nothing
                       nothing nothing 1       nothing
                       nothing nothing 1       nothing])
    @test add_hyperedge!(h1_1) == 5
    @test size(remove_hyperedge!(h1_1, 5))[2] == 4
    @test add_vertex!(h1_1) == 5
    @test add_hyperedge!(h1_1) == 5
    hp = prune_hypergraph(h1_1)
    @test size(hp)[1] == 3 && size(h)[1] == 4
    @test size(hp)[2] == 3 && size(h)[1] == 4
    prune_hypergraph!(h1_1)
    @test size(h1_1)[1] == 3
    @test size(h1_1)[2] == 3
end;

@testset "SimpleHypergraphs BipartiteView       " begin
    h2 = deepcopy(h1)

    @test LightGraphs.nv(LightGraphs.zero(BipartiteView{Int})) == 0

    b = BipartiteView(h2)
    @test LightGraphs.edgetype(b) == LightGraphs.SimpleGraphs.SimpleEdge{Int}
    @test LightGraphs.has_vertex(b, 0) == false
    @test LightGraphs.has_vertex(b, 1) == true
    @test LightGraphs.has_edge(b, 1, 1) == false
    @test LightGraphs.nv(LightGraphs.zero(b)) == 0

    @test LightGraphs.is_directed(b) == false
    @test LightGraphs.is_directed(typeof(b)) == false
    @test LightGraphs.eltype(b) == Int

    @test sum(LightGraphs.adjacency_matrix(LightGraphs.SimpleGraph(b))) == 18

    @test sort(collect(LightGraphs.outneighbors(b,5))) == [7,9]
    @test sort(collect(LightGraphs.outneighbors(b,1))) == [6]
    @test sort(collect(LightGraphs.inneighbors(b,9))) == [3,4,5]

    @test Set(LightGraphs.vertices(b)) == Set(1:LightGraphs.nv(b))

    @test shortest_path(b,1,5) == [1,3,5]
    @test LightGraphs.is_weakly_connected(b) == true

    @test add_vertex!(h2) == 6
    @test add_hyperedge!(h2) == 5
    h2[5,5] = 1
    h2[6,5] = 1

    @test shortest_path(b,1,6) == [1,3,5,6]

    bipartite_graph = LightGraphs.SimpleGraph(b)

    @test LightGraphs.SimpleGraphs.fadj(bipartite_graph)==LightGraphs.SimpleGraphs.fadj(b)
    @test LightGraphs.nv(b) == 11
    @test LightGraphs.ne(b) == 11

    @test sort!(LightGraphs.SimpleGraphs.fadj(b,1)) == [7]
    @test sort!(LightGraphs.SimpleGraphs.fadj(b,2)) == [7,9]
end;

@testset "SimpleHypergraphs TwoSectionView      " begin

    ht = Hypergraph{Float64}(3,3)
    ht[1:2,1:2] .= 2.
    ht[:, 3] .= 2.

    add_vertex!(h1)
    add_hyperedge!(h1)
    h1[5,5] = 1
    h1[6,5] = 1

    @test LightGraphs.nv(LightGraphs.zero(TwoSectionView{Int})) == 0

    t = TwoSectionView(h1)
    @test LightGraphs.edgetype(t) == LightGraphs.SimpleGraphs.SimpleEdge{Int}
    @test LightGraphs.has_vertex(t, 0) == false
    @test LightGraphs.has_vertex(t, 1) == true
    @test LightGraphs.nv(LightGraphs.zero(t)) == 0

    @test LightGraphs.is_directed(t) == false
    @test LightGraphs.is_directed(typeof(t)) == false
    @test LightGraphs.eltype(t) == Int

    @test LightGraphs.nv(t) == 6
    @test LightGraphs.ne(t) == 8

    @test sort(LightGraphs.all_neighbors(t, 1)) == [2,3]
    @test sort(LightGraphs.outneighbors(t, 5)) == [3,4,6]
    @test sort(LightGraphs.inneighbors(t, 4)) == [2,3,5]
    @inferred LightGraphs.all_neighbors(t, 1)

    @test LightGraphs.has_edge(t, 1, 2) == true
    @test LightGraphs.has_edge(t, 1, 5) == false

    @test sum(LightGraphs.adjacency_matrix(LightGraphs.SimpleGraph(t))) == 16
    @test shortest_path(t,1,5) == [1,3,5]
    @test LightGraphs.is_weakly_connected(t) == true

    @test SimpleHypergraphs.add_vertex!(h1) == 7
    h1[7,5] = 1

    @test shortest_path(t,1,6) == [1,3,5,6]

    @test LightGraphs.ne(t) == 10
    @test LightGraphs.nv(t) == 7
    @test sort(LightGraphs.outneighbors(t, 5)) == [3,4,6,7]

    @test sum(LightGraphs.adjacency_matrix(LightGraphs.SimpleGraph(t))) == 20

    Random.seed!(0);
    g = LightGraphs.erdos_renyi(8, 0.3)
    h_from_g = Hypergraph(g)
    @test LightGraphs.adjacency_matrix(g) == LightGraphs.adjacency_matrix(TwoSectionView(h_from_g))
    @test minimum([sum((h_from_g .== true)[:,n]) for n in 1:6] .== 2)
    @test LightGraphs.modularity(g,[1,1,2,2,3,3,4,4]) ≈ modularity(h_from_g, Set.([[1,2],[3,4],[5,6],[7,8]]))
    @test LightGraphs.SimpleGraphs.fadj(g) == LightGraphs.SimpleGraphs.fadj(TwoSectionView(h_from_g))
end;

@testset "SimpleHypergraphs random-models       " begin

    Hᵣ = random_model(5,5)
    @test nhv(Hᵣ) == 5
    @test nhe(Hᵣ) == 5
    @test  all(length.(Hᵣ.v2he) .> 0)
    @test  all(length.(Hᵣ.v2he) .<= 5)

    Hᵣ2 = random_model(5,0)
    add_hyperedge!(Hᵣ2;vertices=Dict(2 => true, 4 => true))
    @test nhv(Hᵣ2) == 5
    @test nhe(Hᵣ2) == 1

    Hᵣ3  = random_model(5,1)
    @test nhe(Hᵣ3) == 1

    Hκ = random_kuniform_model(5, 5, 3)
    @test nhv(Hκ) == 5
    @test nhe(Hκ) == 5
    @test all(length.(Hκ.he2v) .== 3)

    Hδ = random_dregular_model(5, 5, 3)
    @test nhv(Hδ) == 5
    @test nhe(Hδ) == 5
    @test all(length.(Hδ.v2he) .== 3)

    H∂ = random_preferential_model(20, 0.5)
    @test nhv(H∂) == 20
end;

@testset "SimpleHypergraphs modularity          " begin
    Random.seed!(1234);
    hg = Hypergraph{Bool}(10, 12)
    for i in eachindex(hg)
        if rand() < 0.2
            hg[i] = true
        end
    end

    cfmr = CFModularityRandom(3,10000)

    @test findcommunities(hg,cfmr).bm ≈ 0.21505688117829677
    @test modularity(hg,  Set.([1:10])) == 0.0
    Random.seed!(1234);
    @test typeof(randompartition(hg, 2)) == Vector{Set{Int}}
    @test length(randompartition(hg, 2)) == 2
    @test sort(vcat(collect.(randompartition(hg, 2))...))==1:nhv(hg)

    hh = Hypergraph{Bool}(7,4)
    hh[1,1] = true
    hh[2,1:2] .= true
    hh[3,1:3] .= true
    hh[4,4] = true
    hh[5:6,3] .= true
    @test nhv(hh) == 7
    @test nhe(hh) == 4

    @test modularity(hh,Set.([[1,2,3],[4],[5],[6],[7]])) ≈ 223/972
    @test modularity(hh,Set.([[1,2,3],[4,5,6,7]])) ≈ 14/72
    @test modularity(hh,Set.([[1,2,3,5,6],[4,7]])) ≈ 16/81
    @test modularity(hh,Set.([[1,2,3,5,6],[4],[7]])) ≈ 16/81
    @test modularity(hh, Set.([1:nhv(hh)])) == 0.0
    ha = SimpleHypergraphs.HypergraphAggs(hh)
    @test ha.hes == [3, 2, 3, 1]
    @test ha.max_hes == 3
    @test ha.deg_vs == [1, 2, 3, 1, 1, 1, 0]
    @test ha.volV == 9
    @test modularity(hh, Set.([[1,2,3],[4],[5],[6],[7]]), ha) ≈ 223/972
    cfmr = CFModularityRandom(2,10000)
    @test cfmr.n==2
    @test cfmr.reps == 10000
    @test findcommunities(hh,cfmr).bm ≈ 16/81
    Random.seed!(1234);
    cnm = CFModularityCNMLike(1000)
    @test cnm.reps == 1000
    @test findcommunities(hh, CFModularityRandom(4,10000)).bm ≈ findcommunities(hh, cnm).bm
    Random.seed!(0);
    @test any(isapprox.(findcommunities(hh, cnm).bm, [16/81, 223/972]))
end;

@testset "SimpleHypergraphs randomized tests    " begin
    Random.seed!(0)
    N = 500
    res = Vector{Bool}(undef, N)
    for i in 1:N
        m1 = CFModularityCNMLike(100)
        m2 = CFModularityRandom(2,100)
        r = rand([repeat([nothing],6)..., true], 12, 8)
        hh = Hypergraph(r)
        bm1 = findcommunities(hh, m1).bm
        bm2 = findcommunities(hh, m2).bm
        res[i] = (bm1 > bm2)
    end
    @test sum(res) >= N*0.75
end

@testset "SimpleHypergraphs label propagation   " begin
    Random.seed!(1234);
    hg = Hypergraph{Bool}(10, 12)
    for i in eachindex(hg)
        if rand() < 0.2
            hg[i] = true
        end
    end

    cflp = CFLabelPropagationFinder(100, 1234)
    @test_throws AssertionError findcommunities(hg, cflp)

    h = Hypergraph(11, 2)
    h[1:5, 1] .= true
    h[5:11, 2] .= true

    comms = findcommunities(h, cflp)
    @test comms.np == [Set([7, 9, 10, 11, 8, 5, 6]), Set([4, 2, 3, 1])]
    @test comms.hep == Set[Set([2]), Set([1])]

    add_hyperedge!(h)
    comms = findcommunities(h, cflp)

    @test comms.helabels == [4, 7, -1]
end;

@testset "SimpleHypergraphs randomwalk          " begin
    h1 = Hypergraph{Float64}(5,4)
    h1[1:3,1] .= 1.5
    h1[3,4] = 2.5
    h1[2,3] = 3.5
    h1[4,3:4] .= 4.5
    h1[5,4] = 5.5
    h1[5,2] = 6.5
    # randomized, unseeded tests
    w1 = countmap([random_walk(h1, 1) for _ in 1:10^6])
    @test keys(w1) == Set([1,2,3])
    @test -(extrema(values(w1))...) > -10000
    w2 = countmap([random_walk(h1, 2) for _ in 1:10^6])
    @test keys(w2) == Set([1,2,3,4])
    @test abs((w2[2]-w2[4]) - w2[1]) < 10000
    @test abs(w2[1]-w2[3]) < 10000
    w5 = countmap([random_walk(h1, 5) for _ in 1:10^6])
    @test keys(w5) == Set([3,4,5])
    @test abs(w5[3]-w5[4]) < 10000
    @test abs(w5[5]-w5[4]-500000) < 10000
    @test_throws ArgumentError random_walk(h1, 0)
end

@testset "SimpleHypergraphs connected components" begin
    bip = LightGraphs.SimpleGraph(BipartiteView(h1))
    cc = LightGraphs.connected_components(bip)
    filter!.(x -> x <= nhv(h1), cc)
    filter!(!isempty, cc)

    cc2 = SimpleHypergraphs.get_connected_components(h1)
    @test sort!(sort!.(cc)) == sort!(sort!.(cc2))
    @test typeof(cc2) == Vector{Vector{Int}}
end

@testset "SimpleHypergraphs hypernetx bridge    " begin

    if (!SimpleHypergraphs.support_hypernetx())
        @warn "HyperNetX is not installed. Skipping hypernetx tests"
        return
    end

    h_hnx = SimpleHypergraphs._convert_to_hnx(h1)
    data = Dict{String, Array{Int, 1}}(
        "1" => [1, 2, 3],
        "2" => [5],
        "3" => [2, 4],
        "4" => [3, 4, 5],
        "5" => [5, 6, 7]
    )
    h2 = SimpleHypergraphs.hnx.Hypergraph(data)
    @test h_hnx == h2

    h_hnx =
        SimpleHypergraphs._convert_to_hnx(
            h1,
            node_labels = Dict{Int, String}(
                1=>"A", 2=>"B", 3=>"C", 4=>"D", 5=>"E", 6=>"F", 7=>"G"),
            edge_labels = Dict{Int, String}(
                1=>"HE1", 2=>"HE2", 3=>"HE3", 4=>"HE4", 5=>"HE5"
            ))
    data = Dict{String, Array{String, 1}}(
        "HE1" => ["A", "B", "C"],
        "HE2" => ["E"],
        "HE3" => ["B", "D"],
        "HE4" => ["C", "D", "E"],
        "HE5" => ["E", "F", "G"]
    )
    h2 = SimpleHypergraphs.hnx.Hypergraph(data)
    @test h_hnx == h2

    @test SimpleHypergraphs.get_next_div_id() == 1
    @test SimpleHypergraphs.get_next_div_id() == 2
end;

@testset "SimpleHypergraphs conductance         " begin
  h = Hypergraph{Float64, Int}(5,4)
  h[1:3,1] .= 1
  h[3,4] = 1
  h[2,3] = 1
  h[4,3:4] .= 1
  h[5,4] = 1
  h[5,2] = 1
  @test SimpleHypergraphs.conductance(h, Set(1:3)) == 5 / 5
  @test SimpleHypergraphs.conductance(h, Set(1)) == 3 / 1
  @test SimpleHypergraphs.conductance(h, Set([1, 4])) == 8 / 3
  @test SimpleHypergraphs.conductance(h, Set(2:5)) == 3 / 8
  h = Hypergraph{Int}(6, 7)
  h[1:3, 1] .= 1
  h[1:2, 2] .= 1
  h[1, 3] = 1
  h[3, 3] = 1
  h[3:4, 4] .= 1
  h[4:6, 5] .= 1
  h[4:5, 6] .= 1
  h[5:6, 7] .= 1
  @test SimpleHypergraphs.conductance(h, Set(1:3)) == 2 / 8
  @test SimpleHypergraphs.conductance(h, Set([1, 4])) == 14 / 6
  @test_throws ErrorException SimpleHypergraphs.conductance(h, Set{Int}())
  @test_throws ErrorException SimpleHypergraphs.conductance(h, Set(1:nhv(h)))
end;

@testset "SimpleHypergraphs dual                " begin
    m = [
          1         nothing   nothing   4
          1         2         3         nothing
          1         2         3         4
          nothing   2         3         nothing
          nothing   nothing   3         nothing
          nothing   nothing   nothing   4
    ]

    v_meta = Array{Union{Nothing, Char}, 1}(collect('a':'f'))
    he_meta = Array{Union{Nothing, Symbol}, 1}(Symbol.(collect('A':'D')))

    h = Hypergraph{Int, Char, Symbol}(m; v_meta=v_meta, he_meta=he_meta)
    h_dual = dual(h)

    @test nhv(h_dual) == nhe(h)
    @test nhe(h_dual) == nhv(h)

    @test h.v_meta == h_dual.he_meta
    @test h.he_meta == h_dual.v_meta

    m_dual = Matrix(h_dual)
    m_dual[m_dual .== nothing] .= 0
    m[m .== nothing] .= 0

    @test m == transpose(m_dual)

    @test_throws AssertionError dual(Hypergraph(0, 0))
end;

@testset "SimpleHypergraphs s-distance          " begin
    h = Hypergraph{Int}(6,4)

    h[1:3, 1] .= 1
    h[2:4, 2] .= 2
    h[2:5, 3] .= 3
    h[1, 4] = 4
    h[3, 4] = 4
    h[6, 4] = 4

    @test adjacency_matrix(Matrix(h); s=1) == [
                                                0  1  2  0  0  1
                                                1  0  3  2  1  0
                                                2  3  0  2  1  1
                                                0  2  2  0  1  0
                                                0  1  1  1  0  0
                                                1  0  1  0  0  0
                                               ]

    @test adjacency_matrix(Matrix(h); s=1, weighted=false) == [
                                                            0  1  1  0  0  1
                                                            1  0  1  1  1  0
                                                            1  1  0  1  1  1
                                                            0  1  1  0  1  0
                                                            0  1  1  1  0  0
                                                            1  0  1  0  0  0
                                                        ]

    @test edge_adjacency_matrix(h; s=1, weighted=false) == [
                                                            0  1  1  1
                                                            1  0  1  1
                                                            1  1  0  1
                                                            1  1  1  0
                                                        ]

    @test edge_adjacency_matrix(h; s=1) == [
                                            0  2  2  2
                                            2  0  3  1
                                            2  3  0  1
                                            2  1  1  0
                                        ]

    @test adjacency_matrix(Matrix(h); s=2) == [
                                                 0  0  2  0  0  0
                                                 0  0  3  2  0  0
                                                 2  3  0  2  0  0
                                                 0  2  2  0  0  0
                                                 0  0  0  0  0  0
                                                 0  0  0  0  0  0
                                             ]

    @test distance(h, SnodeDistanceDijkstra(1, 5, 1)) == 2
    @test distance(h, SnodeDistanceDijkstra(3, 6, 1)) == 1

    @test distance(h, SnodeDistanceDijkstra(1, 2, 2)) == 2
    @test distance(h, SnodeDistanceDijkstra(1, 6, 2)) == typemax(Int)

    @test distance(h, SedgeDistanceDijkstra(1, 3, 1)) == 1
    @test distance(h, SedgeDistanceDijkstra(2, 3, 3)) == 1
    @test distance(h, SedgeDistanceDijkstra(1, 3, 3)) == typemax(Int)
end;


text:
function func1(a::Int)
    a
end

function SomeModule.func2(a::Int)
    a
end

module MyModule
function func3(a::Int)
    a
end

function func4(a::Int)
    function func5(b::Int)
        b
    end
    func5(a)
end

end


text:
Base.@kwdef struct ServerSecurity
    require_token_for_open_links::Bool = true
end

# ServerSecurity(val::Bool) = ServerSecurity(val, val, val)

# THIS IS NOT USED YET
Base.@kwdef struct ServerConfiguration
    host::AbstractString = "127.0.0.1"
    port::Integer = 1234
    skip_main_menu::Bool = false
    show_file_system::Bool = true
end

# ServerSecurity(val::Bool) = ServerSecurity(val)


text:
using DiffEqDocs, DifferentialEquations
using Test

# Tests are just for docs generation


text:
import GreyDecision.Utility: makeminmax
import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Wpm: wpm

@testset "WPM" begin
    tol = 0.00001
    weights = [0.221, 0.159, 0.175, 0.127, 0.117, 0.201]
    fns = [maximum, minimum, minimum, maximum, minimum, maximum]
    decmat = [
            3.0 12.5 2 120 14 3
            5.0 15 3 110 38 4
            3.0 13 2 120 19 3
            4.0 14 2 100 31 4
            3.0 15 1.5 125 40 4
        ]
    nrow, ncol = size(decmat)
    # Converting white decision matrix to
    # grey decision matrix
    greydecmat = Array{GreyNumber{Float64}, 2}(undef, nrow, ncol)

    for i in 1:nrow
        for j in 1:ncol
            greydecmat[i, j] = GreyNumber(decmat[i, j])
        end
    end

    result = wpm(greydecmat, weights, fns)

    scores = result[:scores]
    @test isapprox(scores[1].a, 0.7975224331331, atol = tol)
    @test isapprox(scores[2].a, 0.7532541470585, atol = tol)
    @test isapprox(scores[3].a, 0.7647463553356, atol = tol)
    @test isapprox(scores[4].a, 0.7873956894791, atol = tol)
    @test isapprox(scores[5].a, 0.7674278741782, atol = tol)

    @test result[:bestIndex] == 1
    @test result[:ranking] == [2, 3, 5, 4, 1]
end


text:
module LE1NAFEMS_examples
using FinEtools
using FinEtools.MeshExportModule
using PyCall

function LE1NAFEMS()
    println("LE1NAFEMS, plane stress.")
    t0 = time()

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD= 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Radius= 1.0*phun("m")
    n=64; # number of elements per side
    tolerance=1.0/n/1000.;#Geometrical tolerance

    fens,fes = Q4block(1.0,pi/2, n, n*2)

    bdryfes = meshboundary(fes);
    icl = selectelem(fens, bdryfes, box=[1.0,1.0,0.0,pi/2],inflate=tolerance);
    for i=1:count(fens)
        t=fens.xyz[i,1]; a=fens.xyz[i,2];
        fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a) (t*2.75+(1-t)*1)*sin(a)];
    end

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    l1 =selectnode(fens; box=[0.0 Inf 0.0 0.0], inflate = tolerance)
    setebc!(u,l1,true,2,0.0)
    l1 =selectnode(fens; box=[0.0 0.0 0.0 Inf], inflate = tolerance)
    setebc!(u,l1,true,1,0.0)

    applyebc!(u)
    numberdofs!(u)

    el1femm =  FEMMBase(IntegData(subset(bdryfes,icl), GaussRule(1, 2)))
    function pfun(forceout::FFltVec, x::FFltMat, J::FFltMat, l::FInt)
        pt= [2.75/3.25*x[1] 3.25/2.75*x[2]]
        copy!(forceout, vec(p*pt/norm(pt)))
        return forceout
    end
    fi = ForceIntensity(FFlt, 2, pfun);
    F2= distribloads(el1femm, geom, u, fi, 2);

    MR = DeforModelRed2DStress
    material=MatDeforElastIso(MR, E, nu)
    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(2, 2)), material)

    K =stiffness(femm, geom, u)
    K=cholesky(K)
    U=  K\(F2)
    scattersysvec!(u,U[:])

    nl=selectnode(fens, box=[2.0,2.0,0.0,0.0],inflate=tolerance);
    thecorneru=zeros(FFlt,1,2)
    gathervalues_asmat!(u,thecorneru,nl);
    thecorneru=thecorneru/phun("mm")
    println("$(time()-t0) [s];  displacement =$(thecorneru) [MM] as compared to reference [-0.10215,0] [MM]")

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
    println("  Target stress: $(fld.values[nl][1]/phun("MEGA*PA")) compared to $(sigma_yD/phun("MEGA*PA"))")

    File =  "a.vtk"
    vtkexportmesh(File, fes.conn, geom.values,
    FinEtools.MeshExportModule.Q4;
    vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)
    true

end # LE1NAFEMS

function LE1NAFEMS_Q4_convergence()
    println("LE1NAFEMS, 2D version."        )
    t0 = time()

    E = 210e3*phun("MEGA*PA");
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Radius = 1.0*phun("m")
    Thickness = 0.1*phun("m")

    sigyderrs = Dict{Symbol, FFltVec}()
    nelems = []

    for extrapolation in [:extrapmean]
        sigyderrs[extrapolation] = FFltVec[]
        nelems = []
        for n in [16, 32, 64, 128, 256, 512]
            tolerance = 1.0/n/1000.; # Geometrical tolerance

            fens,fes = Q4block(1.0, pi/2, n, n)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2], inflate=tolerance);
            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2];
                fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a)];
            end

            geom = NodalField(fens.xyz)
            u = NodalField(zeros(size(fens.xyz,1), 2)) # displacement field

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)

            applyebc!(u)
            numberdofs!(u)

            el1femm =  FEMMBase(IntegData(subset(bdryfes,icl), GaussRule(1, 2), Thickness))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2]]
                forceout .=    vec(p*pt/norm(pt));
                return forceout
            end
            fi = ForceIntensity(FFlt, 2, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            # Note that the material object needs to be created for plane stress.
            MR = DeforModelRed2DStress

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(2, 2), Thickness), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(F2)
            scattersysvec!(u, U[:])

            nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0],inflate=tolerance);
            thecorneru = zeros(FFlt,1,2)
            gathervalues_asmat!(u, thecorneru, nl);
            thecorneru = thecorneru/phun("mm")
            println("$(time()-t0) [s];  displacement =$(thecorneru) [MM] as compared to reference [-0.10215,0] [MM]")

            fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
            println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) vs $(sigma_yD/phun("MPa")) [MPa]")

            println("$(n), $(fld.values[nl,1][1]/phun("MPa"))")
            push!(nelems, count(fes))
            push!(sigyderrs[extrapolation], abs(fld.values[nl,1][1]/sigma_yD - 1.0))
            #     File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            #                FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            #                scalars=[("sigmay", fld.values)])
            # @async run(`"paraview.exe" $File`)
        end
    end

    @pyimport matplotlib.pyplot as plt
    plt.style[:use]("seaborn-whitegrid")
    fig = plt.figure()
    ax = plt.axes()
    ax[:plot](vec(nelems), vec(sigyderrs[:extrapmean]), linestyle="none", marker=:o, label="Extrapolated stress")
    ax[:grid](linestyle="--", linewidth=0.5, color=".25", zorder=-10)
    ax[:set_xlabel]("Number of elements")
    ax[:set_ylabel]("Norm. approx. error")
    plt.legend()
    plt.show()
    # plt.xlim(0.0, 1.0)
    # plt.axis("equal")
    #
    #

end # LE1NAFEMS_Q4_convergence

function LE1NAFEMS_Q8_stress()
    println("LE1NAFEMS, plane stress.")
    t0 = time()

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD= 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Radius= 1.0*phun("m")
    n=20; # number of elements per side
    tolerance=1.0/n/1000.;#Geometrical tolerance

    fens,fes = Q8block(1.0,pi/2, n, n*2)

    bdryfes = meshboundary(fes);
    icl = selectelem(fens, bdryfes, box=[1.0,1.0,0.0,pi/2],inflate=tolerance);
    for i=1:count(fens)
        t=fens.xyz[i,1]; a=fens.xyz[i,2];
        fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a) (t*2.75+(1-t)*1)*sin(a)];
    end

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    l1 =selectnode(fens; box=[0.0 Inf 0.0 0.0], inflate = tolerance)
    setebc!(u,l1,true,2,0.0)
    l1 =selectnode(fens; box=[0.0 0.0 0.0 Inf], inflate = tolerance)
    setebc!(u,l1,true,1,0.0)

    applyebc!(u)
    numberdofs!(u)

    el1femm =  FEMMBase(IntegData(subset(bdryfes,icl), GaussRule(1, 3)))
    function pfun(forceout::FFltVec, x::FFltMat, J::FFltMat, l::FInt)
        pt= [2.75/3.25*x[1] 3.25/2.75*x[2]]
        copy!(forceout, vec(p*pt/norm(pt)))
        return forceout
    end
    fi = ForceIntensity(FFlt, 2, pfun);
    F2= distribloads(el1femm, geom, u, fi, 2);

    MR = DeforModelRed2DStress
    material=MatDeforElastIso(MR, E, nu)

    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(2, 3)), material)

    K =stiffness(femm, geom, u)
    K=cholesky(K)
    U=  K\(F2)
    scattersysvec!(u,U[:])

    nl=selectnode(fens, box=[2.0,2.0,0.0,0.0],inflate=tolerance);
    thecorneru=zeros(FFlt,1,2)
    gathervalues_asmat!(u,thecorneru,nl);
    thecorneru=thecorneru/phun("mm")
    println("$(time()-t0) [s];  displacement =$(thecorneru) [MM] as compared to reference [-0.10215,0] [MM]")

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
    println("  Target stress: $(fld.values[nl][1]/phun("MEGA*PA")) compared to $(sigma_yD/phun("MEGA*PA"))")

    File =  "a.vtk"
    vtkexportmesh(File, fes.conn, geom.values,
    FinEtools.MeshExportModule.Q8;
    vectors=[("u", u.values)], scalars=[("sigy", fld.values)])
    @async run(`"paraview.exe" $File`)
    true

end # LE1NAFEMS_Q8_stress

function allrun()
    println("#####################################################")
    println("# LE1NAFEMS ")
    LE1NAFEMS()
    println("#####################################################")
    println("# LE1NAFEMS_Q4_convergence ")
    LE1NAFEMS_Q4_convergence()
    println("#####################################################")
    println("# LE1NAFEMS_Q8_stress ")
    LE1NAFEMS_Q8_stress()
    return true
end # function allrun

end # module LE1NAFEMS_examples


text:
abstract type Node end

setid!(n::Node, id::Symbol) = setid!(n, valid_id(string(id)))
setid!(n::Node, id::String) = (n.id = id)
getid(n::Node) = n.id

setnext!(n::Node, next::Node) = error("'$(getid(n)) = $(typeof(n))(...) -> next' not implemented for $(typeof(n))")
setnext!(n::Node, nexts::Vector{<:Node}) = error("'$(getid(n)) = $(typeof(n))(...) -> [next_1, ..., next_n]' not implemented for $(typeof(n))")
setnext!(n::Node, name::Symbol, next::Node) = error("'$(getid(n)) = $(typeof(n))(...) -> [next_kw = next, ...]' not implemented for $(typeof(n))")

getmsg(n::Node) = error("Not implemented!")

function valid_id(id::AbstractString)
	valid_char(c) = islowercase(c) | (c == '_') | isdigit(c)
	islowercase(id[1]) || error("First character in Node ID needs to be lower case character.")
	all(valid_char, id) || error("Node ID can only contain lower case letters, digits and underscores.")
	return id
end

abstract type NonInteractiveNode <: Node end
getnext(n::NonInteractiveNode) = error("Not implemented!")

abstract type InteractiveNode <: Node end
getopt(n::InteractiveNode) = error("Not implemented!")
getnext(n::InteractiveNode, opt) = error("Not implemented!")

abstract type StateInteractionNode <: Node end
getnext!(n::StateInteractionNode, state) = error("Not implemented!")
getmsg(n::StateInteractionNode, state) = getmsg(n)

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

mutable struct Start <: NonInteractiveNode
	id::String
	next::Node

	Start() = new()
end
getnext(n::Start) = n.next
setnext!(n::Start, next::Node) = (n.next = next)
getmsg(n::Start) = n.id

mutable struct End <: NonInteractiveNode
	id::String
	text::String

	End(text="End of Action") = (obj = new(); obj.text = strip(text); obj)
end
getnext(n::End) = error("$(typeof(n)) have no children.")
getmsg(n::End) = n.text

mutable struct Dummy <: NonInteractiveNode
	id::String
	next::Node

	Dummy() = new()
end
getnext(n::Dummy) = n.next
setnext!(n::Dummy, next::Node) = (n.next = next)
getmsg(n::Dummy) = "Dummy: $(getid(n))"

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

mutable struct JumpToGraph <: NonInteractiveNode
	id::String
	next::Node
	jump_graph::String

	JumpToGraph(jump_graph) = (obj = new(); obj.jump_graph = jump_graph; obj)
end
getnext(n::JumpToGraph) = n.next
setnext!(n::JumpToGraph, next::Node) = (n.next = next)
getmsg(n::JumpToGraph) = "Jump to: $(n.jump_graph)"

mutable struct ReturnFromGraph <: NonInteractiveNode
	id::String

	ReturnFromGraph() = new()
end
getnext(n::ReturnFromGraph) = error("$(typeof(n)) have no children.")
getmsg(n::ReturnFromGraph) = "Continue from where\nyou jumped to this\n graph."

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

mutable struct PerformAction <: InteractiveNode
	id::String
	action::String
	next::Node

	PerformAction(action="Continue from where\nyou jumped to this\n graph.") = (obj = new(); obj.action = strip(action); obj)
end
getopt(n::PerformAction) = [CMD.Blank()]
getnext(n::PerformAction, opt::CMD.Blank) = n.next
setnext!(n::PerformAction, next::Node) = (n.next = next)
getmsg(n::PerformAction) = n.action

mutable struct BinaryCondition <: InteractiveNode
	id::String
	condition::String
	n_true::Node
	n_false::Node

	BinaryCondition(condition) = (obj = new(); obj.condition = strip(condition); obj)
end
getopt(n::BinaryCondition) = [CMD.True(), CMD.False()]
getnext(n::BinaryCondition, opt::CMD.True) = n.n_true
getnext(n::BinaryCondition, opt::CMD.False) = n.n_false
setnext!(n::BinaryCondition, name::Symbol,next::Node) = setfield!(n,name,next)
getmsg(n::BinaryCondition) = n.condition

mutable struct MultipleChoice <: InteractiveNode
	id::String
	conditions::String
	nexts::Vector{<:Node}

	MultipleChoice(conditions) = (obj = new(); obj.conditions = strip(conditions); obj)
end
getopt(n::MultipleChoice) = CMD.Option.(1:length(n.nexts))
getnext(n::MultipleChoice, opt::CMD.Option) = n.nexts[opt.opt]
setnext!(n::MultipleChoice, nexts::Vector{<:Node}) = (n.nexts = nexts)
getmsg(n::MultipleChoice) = n.conditions

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

struct QuellerGraph
	root_node::Node
	nodes::Vector{<:Node}
	source_file::String
end

getid(g::QuellerGraph) = getid(root(g))
root(g::QuellerGraph) = g.root_node

function load_graphs_from_file(file)
	starts, nodes = include(file)
	return [QuellerGraph(start, nodes, file) for start in starts]
end

function load_graphs(files...)
	graphs = mapreduce(load_graphs_from_file, vcat, files)

	unique, conflicts = unique_root_ids(graphs)
	!unique && error("Error loading graphs, not all root node ids are unique. Conflictings graphs:\n$(strvec2str(conflicts))")

	jumps_exists, conflicts = all_jump_points_exists(graphs)
	!jumps_exists && error("Warning, not all jumps exists. Non-existing jumps:\n$(strvec2str(conflicts))")

	return Dict(getid(g) => g for g in graphs)
end

function unique_root_ids(graphs)
	conflicts = not_unique(graphs, (g,j) -> (getid(g) == getid(j)) )
	return isempty(conflicts), map(g -> getid(g)*" : "*g.source_file, conflicts)
end

function get_jump_points(g)
	jump_nodes = filter(n -> isa(n, JumpToGraph), g.nodes)
	return tuple.(getfield.(jump_nodes, :jump_graph),[ g.source_file])
end

function all_jump_points_exists(graphs)
	jump_points = mapreduce(get_jump_points, vcat, graphs)
	graph_ids = getid.(graphs)
	jump_not_exists(j) = !(j[1] in graph_ids)
	conflicts = filter(jump_not_exists, jump_points)
	return isempty(conflicts), map(c -> c[1]*" : "*c[2], conflicts)
end

function get_graphs_not_jumped_to(graphs)
	jump_points = getindex.(mapreduce(get_jump_points, vcat, graphs),1)
	graph_not_jumped_to(g) = !(getid(g) in jump_points)
	unjumped = filter(graph_not_jumped_to, collect(graphs))
	return map(g -> getid(g)*" : "*g.source_file, unjumped)
end

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

macro graphs(node_block)
	return quote
		__nodes = Dict{Symbol, Node}()

		__edges_single = Vector{Tuple{Node,Symbol}}()
		__edges_multi = Vector{Tuple{Node,Vector{Symbol}}}()
		__edges_named = Vector{Tuple{Node,Symbol,Symbol}}()

		$node_block

		try
			for (n,next) in __edges_single
				setnext!(n, __nodes[next])
			end
			for (n,nexts) in __edges_multi
				setnext!(n, [__nodes[next] for next in nexts])
			end
			for (n,name,next) in __edges_named
				setnext!(n, name, __nodes[next])
			end
		catch e
			e isa KeyError && error("Can't link to node '$(string(e.key))', node not found.")
			rethrow()
		end

		nodes = collect(values(__nodes))
		starts = filter(n -> n isa Start, nodes)

		undef_nodes = filter(!(fieldsdefined), nodes)
		for n in undef_nodes
			error("Not all edges set on node '$(getid(n))' (or it have another un-initialized field).")
		end

		(starts, nodes)
	end
end

macro node(exp)
	err_msg = """
	Invalid node format. Valid formats are:
		@node id = NodeType(...) -> []
		@node id = NodeType(...) -> next
		@node id = NodeType(...) -> [next_1, ..., next_n]
		@node id = NodeType(...) -> [next_kw = next, ...]
	"""
	exp.head != :(=) && error(err_msg)
	exp.args[2].head != :(->) && error(err_msg)

	id = exp.args[1]
	node = exp.args[2].args[1]
	next = exp.args[2].args[2].args[2]

	isa_vecexpr(e) = e isa Expr && e.head in [:vect, :hcat, :vcat]

	if next isa Symbol
		edge_push = :(push!(__edges_single, (n, $(Meta.quot(next)))))

	elseif isa_vecexpr(next) && isempty(next.args)
		edge_push = :()

	elseif isa_vecexpr(next) && all(a -> (a isa Symbol), next.args)
		next = [s for s in next.args]
		edge_push = :(push!(__edges_multi, (n, $(Meta.quot(next)))))

	elseif isa_vecexpr(next) && all(a -> a isa Expr, next.args) && all(a -> a.head == :(=), next.args)
		edge_push = :()
		for ass in next.args
			name = ass.args[1]
			next = ass.args[2]
			edge_push = :($edge_push; push!(__edges_named, (n, $(Meta.quot(name)), $(Meta.quot(next)))) )
		end

	else
		error(err_msg)
	end

	return quote
		id = $(Meta.quot(id))
		id in keys(__nodes) && error("Node '$(string(id))' already exists.")

		n = $node
		setid!(n,id)
		__nodes[id] = n

		$edge_push
	end
end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
# Autogenerated wrapper script for XPA_jll for i686-linux-gnu
export libxpa, xpaaccess, xpaget, xpainfo, xpamb, xpans, xpaset

JLLWrappers.@generate_wrapper_header("XPA")
JLLWrappers.@declare_library_product(libxpa, "libxpa.so.1")
JLLWrappers.@declare_executable_product(xpaaccess)
JLLWrappers.@declare_executable_product(xpaget)
JLLWrappers.@declare_executable_product(xpainfo)
JLLWrappers.@declare_executable_product(xpamb)
JLLWrappers.@declare_executable_product(xpans)
JLLWrappers.@declare_executable_product(xpaset)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libxpa,
        "lib/libxpa.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

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

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

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

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

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("../src/fdfd_local_field.jl")

using DelimitedFiles

fname="examples/data/test_data__0.01.csv"
data_array = readdlm(fname, ',')
start_index = findnext('.', fname, 1) - 1
end_index = findnext('.', fname, start_index+2) - 1
δ = parse(Float64, fname[start_index:end_index])

real_gd, imag_gd = get_val_gradient_from_data(data_array, δ=δ)
writedlm(fname[1:end-4]*"_real_gd.csv", real_gd, ',')
writedlm(fname[1:end-4]*"_imag_gd.csv", imag_gd, ',')


text:
using JLD2, Requires, Pkg.TOML, DataDeps
using Flux: loadparams!

export list_pretrains, load_pretrained, @pretrained_str

const pretrained_configs = open(TOML.parse, joinpath(@__DIR__, "pretrains.toml"))

"""
    readckpt(path)

Load weights from a tensorflow checkpoint file into a Dict.
"""
readckpt(path) = error("readckpt require TensorFlow.jl installed. run `Pkg.add(\"TensorFlow\"); using TensorFlow`")

@init @require TensorFlow="1d978283-2c37-5f34-9a8e-e9c0ece82495" begin
    import .TensorFlow
    function readckpt(path)
        weights = Dict{String, Array}()
        TensorFlow.init()
        ckpt = TensorFlow.pywrap_tensorflow.x.NewCheckpointReader(path)
        shapes = ckpt.get_variable_to_shape_map()

        for (name, shape) ∈ shapes
            # Ignore training related data - e.g. learning rates
            # Also ignore scalars and other variables that aren't stored correctly in the checkpoint (shape == Any[])
            (occursin("training", name) || shape == Any[]) && continue
            weight = ckpt.get_tensor(name)
            if length(shape) == 2
                weight = collect(weight')
            end
            weights[name] = weight
        end

        weights
    end
end

"""
    ckpt2_to_jld2(ckptpath::String, ckptname::String, savepath::String)

Loads the pre-trained model weights from a TensorFlow checkpoint and saves to JLD2
"""
function ckpt_to_jld2(ckptpath::String, ckptname::String; savepath::String="./")
    weights = readckpt(joinpath(ckptpath, ckptname))
    jld2name = normpath(joinpath(savepath, ckptname[1:end-5]*".jld2"))
    @info "Saving the model weights to $jld2name"
    JLD2.@save jld2name weights
end

function register_configs(configs)
    for (model_name, config) in pairs(configs)
        model_desc = Transformers.Pretrain.description(config["description"], config["host"], config["link"])
        checksum = config["checksum"]
        url = config["url"]
        dep = DataDep(model_name, model_desc, url, checksum;
                      fetch_method=Transformers.Datasets.download_gdrive)
        DataDeps.register(dep)
    end
end

"""
    list_pretrains()

List all the available pre-trained models.
"""
function list_pretrains()
    println.(keys(pretrained_configs))
    return
end

"""
    load_pretrained(path)

Loads a pre-trained Music Transformer model.
"""
function load_pretrained(model_name::String)
    if model_name ∉ keys(pretrained_configs)
        error("""Invalid model.
               Please try list_pretrains() to check the available pre-trained models""")
    end

    model_config = pretrained_configs[model_name]
    loader = loading_method(Val(Symbol(model_config["model_type"])))

    model_path = @datadep_str("$model_name/$model_name.jld2")
    if !endswith(model_path, ".jld2")
        error("""Invalid file. A .jld2 file is required to load the model.
                 If this is a tensorflow checkpoint file, run ckpt_to_jld2 to convert it.""")
    end

    JLD2.@load model_path weights
    loader(weights, model_config)
end

loading_method(::Val{:unconditional}) = load_unconditional_musictransformer
loading_method(::Val{:melodyconditioned}) = load_melodyconditioned_musictransformer

# Macro from Transformers.jl
macro pretrained_str(name)
    :(load_pretrained($(esc(name))))
end

function load_transformer_encoder!(model, weights, num_layers::Int)
    layers = keys(weights)
    for i = 1:num_layers
        # Get all layers in a block
        block = filter(l->occursin("layer_$(i-1)/", l), layers)
        for k ∈ block
            if occursin("self_attention", k)
                if occursin("q/kernel", k)
                    loadparams!(model[i].mh.iqproj.W, [weights[k]])
                elseif occursin("k/kernel", k)
                    loadparams!(model[i].mh.ikproj.W, [weights[k]])
                elseif occursin("v/kernel", k)
                    loadparams!(model[i].mh.ivproj.W, [weights[k]])
                elseif occursin("output_transform/kernel", k)
                    loadparams!(model[i].mh.oproj.W, [weights[k]])
                elseif occursin("layer_norm_scale", k)
                    loadparams!(model[i].mhn.diag.α, [weights[k]])
                elseif occursin("layer_norm_bias", k)
                    loadparams!(model[i].mhn.diag.β, [weights[k]])
                else
                    @warn "Unknown variable: $k"
                end
            elseif occursin("ffn", k)
                if occursin("conv1/kernel", k)
                    loadparams!(model[i].pw.din.W, [weights[k]])
                elseif occursin("conv1/bias", k)
                    loadparams!(model[i].pw.din.b, [weights[k]])
                elseif occursin("conv2/kernel", k)
                    loadparams!(model[i].pw.dout.W, [weights[k]])
                elseif occursin("conv2/bias", k)
                    loadparams!(model[i].pw.dout.b, [weights[k]])
                elseif occursin("layer_norm_scale", k)
                    loadparams!(model[i].pwn.diag.α, [weights[k]])
                elseif occursin("layer_norm_bias", k)
                    loadparams!(model[i].pwn.diag.β, [weights[k]])
                else
                    @warn "Unknown variable: $k"
                end
            else
                @warn "Unknown variable: $k"
            end
        end
    end

    # Captures both encoder/layer_prepostprocess and decoder/layer_prepostprocess
    output_norm = filter(l->occursin("coder/layer_prepostprocess/layer_norm", l), layers)
    block_index = num_layers + 1
    for k ∈ output_norm
        if occursin("layer_norm_scale", k)
            loadparams!(model[block_index].diag.α, [weights[k]])
        elseif occursin("layer_norm_bias", k)
            loadparams!(model[block_index].diag.β, [weights[k]])
        end
    end
end

function load_transformer_decoder!(model, weights, num_layers::Int)
    # Except for the encoder-decoder attention layer in the decoder,
    # the rest of the layers are similar for both the encoder and decoder
    # So we load the rest of the layers using the load_transformer_encoder
    layers = filter(kv -> !occursin("encdec_attention", kv.first), weights)
    load_transformer_encoder!(model, layers, num_layers)

    # Load the encdec_attention layer
    layers = keys(weights)
    for i = 1:num_layers
        # Get all layers in a block
        block = filter(l->occursin("layer_$(i-1)/", l), layers)
        for k ∈ block
            if occursin("encdec_attention", k)
                if occursin("q/kernel", k)
                    loadparams!(model[i].imh.iqproj.W, [weights[k]])
                elseif occursin("k/kernel", k)
                    loadparams!(model[i].imh.ikproj.W, [weights[k]])
                elseif occursin("v/kernel", k)
                    loadparams!(model[i].imh.ivproj.W, [weights[k]])
                elseif occursin("output_transform/kernel", k)
                    loadparams!(model[i].imh.oproj.W, [weights[k]])
                elseif occursin("layer_norm_scale", k)
                    loadparams!(model[i].imhn.diag.α, [weights[k]])
                elseif occursin("layer_norm_bias", k)
                    loadparams!(model[i].imhn.diag.β, [weights[k]])
                else
                    @warn "Unknown variable: $k"
                end
            end
        end
    end
end

function load_unconditional_musictransformer(weights, config)
    num_layers = config["num_layers"]
    heads = config["heads"]
    depth = config["depth"]
    ffn_depth = config["ffn_depth"]

    mt = UnconditionalMusicTransformer(depth, heads, ffn_depth, num_layers)

    # Load encoder layers
    load_transformer_encoder!(mt.body, weights, num_layers)

    # Load embedding
    embedding = weights["transformer/symbol_modality_310_512/shared/weights_0"]
    loadparams!(mt.embedding[1], [embedding])

    # This model shares embedding and softmax weights
    # TODO: use Base.lastindex once it's defined for Stack
    # 16 transformer body blocks, 1 output layer norm, + 1 is the last embedding-softmax layer
    last_layer = num_layers + 1 + 1
    loadparams!(mt.body[last_layer].W, [embedding'])

    mt
end

function load_melodyconditioned_musictransformer(weights, config)
    num_layers = config["num_layers"]
    heads = config["heads"]
    depth = config["depth"]
    ffn_depth = config["ffn_depth"]

    mt = MelodyConditionedMusicTransformer(depth, heads, ffn_depth, num_layers, num_layers)
    encoder_weights = filter(kv -> occursin("encoder", kv.first), weights)
    load_transformer_encoder!(mt.encoder, encoder_weights, num_layers)
    decoder_weights = filter(kv -> occursin("decoder", kv.first), weights)
    load_transformer_decoder!(mt.decoder, decoder_weights, num_layers)

    # Load embeddings
    input_embedding = weights["transformer/symbol_modality_92_512/input_emb/weights_0"]
    loadparams!(mt.encoder_embedding[1], [input_embedding])

    target_space_embedding = weights["transformer/body/target_space_embedding/kernel"]
    loadparams!(mt.encoder_embedding[2], [target_space_embedding])

    target_embedding = weights["transformer/symbol_modality_310_512/target_emb/weights_0"]
    loadparams!(mt.decoder_embedding[1], [target_embedding])

    softmax_embedding = weights["transformer/symbol_modality_310_512/softmax/weights_0"]
    # $num_layers transformer body blocks, 1 postprocess layernorm, + 1 is the last embedding layer
    last_layer = num_layers + 1 + 1
    loadparams!(mt.decoder[last_layer], [softmax_embedding'])
    mt
end


text:
module pindex

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

using Distributions, Parameters, LinearAlgebra, Compat, Test

function FixedLaspeyres(p,q)
N=size(p,1)
T=size(p,2)
Pl=zeros(N+1,T)
s=zeros(N,T)
   for t in 1:T
        for n in 2:N+1
        s[n-1,t] = (p[n-1,t]*q[n-1,t])/p[:,t]'q[:,t]
        Pl[n,t]=Pl[n-1,t]+p[n-1,t]/p[n-1,1]*s[n-1,1]
        end
    end
return Pl[N+1,:]
end

function FixedPaasche(p,q)
N=size(p,1)
T=size(p,2)
Pp=zeros(N+1,T)
s=zeros(N,T)
   for t in 1:T
        for n in 2:N+1
        s[n-1,t] = p[n-1,t]*q[n-1,t]/p[:,t]'q[:,t]
        Pp[n,t]=Pp[n-1,t]+p[n-1,1]/p[n-1,t]*s[n-1,t]
        end
    end
return Pp[N+1,:].^(-1)
end

function FixedFisher(p,q)
     Pf=(FixedPaasche(p,q).^(1/2)).*(FixedLaspeyres(p,q).^(1/2))
return Pf
end

function FixedTornqvist(p,q)
N=size(p,1)
T=size(p,2)
Pt=zeros(N+1,T)
s=zeros(N,T)
   for t in 1:T
        for n in 2:N+1
        s[n-1,t] = p[n-1,t]*q[n-1,t]/p[:,t]'q[:,t]
        Pt[n,t]=(1/2)*(s[n-1,1]+s[n-1,t])*log(p[n-1,t]/p[n-1,1])
        end
    end
return exp.(Pt[N+1,:])
end

function ChainedLaspeyres(p,q)
N=size(p,1)
T=size(p,2)
Pl=zeros(N+1,T)
Plc=ones(N+1,T)
s=zeros(N,T)
   for t in 2:T
        for n in 2:N+1
        s[n-1,t-1] = p[n-1,t-1]*q[n-1,t-1]/p[:,t-1]'q[:,t-1]
        Pl[n,t]=Pl[n-1,t]+p[n-1,t]/p[n-1,t-1]*s[n-1,t-1]
        Plc[n,t]=Plc[n,t-1]*Pl[n,t]
        end
    end
return Plc[N+1,:]
end

function ChainedPaasche(p,q)
N=size(p,1)
T=size(p,2)
Pp=zeros(N+1,T)
Ppc=ones(N+1,T)
s=zeros(N,T)
   for t in 2:T
        for n in 2:N+1
        s[n-1,t] = (p[n-1,t]*q[n-1,t])/p[:,t]'q[:,t]
        Pp[n,t]=Pp[n-1,t]+p[n-1,t-1]/p[n-1,t]*s[n-1,t]
        Ppc[n,t]=Ppc[n,t-1]*(Pp[n,t])^(-1)
        end
    end
return Ppc[N+1,:]
end

function ChainedFisher(p,q)
     Pfc=(ChainedPaasche(p,q).^(1/2)).*(ChainedLaspeyres(p,q).^(1/2))
return Pfc
end

function ChainedTornqvist(p,q)
N=size(p,1)
T=size(p,2)
Pt=zeros(N+1,T)
Ptc=ones(N+1,T)
s=zeros(N,T)
   for t in 2:T
        for n in 2:N+1
        s[n-1,t-1] = p[n-1,t-1]*q[n-1,t-1]/p[:,t-1]'q[:,t-1]
        s[n-1,t] = (p[n-1,t]*q[n-1,t])/p[:,t]'q[:,t]
        Pt[n,t]=Pt[n-1,t]+(1/2)*(s[n-1,t-1]+s[n-1,t])*log(p[n-1,t]/p[n-1,t-1])
        Ptc[n,t]=Ptc[n,t-1]*exp(Pt[n,t])
        end
    end
return Ptc[N+1,:]
end

function FixLaspeyres(p,q)
N=size(p,1)
T=size(p,2)
Pl=zeros(T)
    for t in 1:T
    Pl[t]=p[:,t]'q[:,1]/p[:,1]'q[:,1]
    end
return Pl
end

function FixPaasche(p,q)
N=size(p,1)
T=size(p,2)
Pl=zeros(T)
    for t in 1:T
    Pl[t]=p[:,t]'q[:,t]/p[:,1]'q[:,t]
    end
return Pl
end

function FixFisher(p,q)
     Pf=(FixPaasche(p,q).^(1/2)).*(FixLaspeyres(p,q).^(1/2))
return Pf
end

function ChainLaspeyres(p,q)
N=size(p,1)
T=size(p,2)
Pl=ones(T)
Plc=ones(T)
    for t in 2:T
    Pl[t]=p[:,t]'q[:,t-1]/p[:,t-1]'q[:,t-1]
    Plc[t]=Plc[t-1]*Pl[t]
    end
return Plc
end

function ChainPaasche(p,q)
N=size(p,1)
T=size(p,2)
Pp=ones(T)
Ppc=ones(T)
    for t in 2:T
    Pp[t]=p[:,t]'q[:,t]/p[:,t-1]'q[:,t]
    Ppc[t]=Ppc[t-1]*Pp[t]
    end
return Ppc
end

function ChainFisher(p,q)
    Pfc=(ChainPaasche(p,q).^(1/2)).*(ChainLaspeyres(p,q).^(1/2))
return Pfc
end

export FixedLaspeyres,  FixedPaasche, FixedFisher, FixedTornqvist, ChainedLaspeyres, ChainedPaasche, ChainedFisher, ChainedTornqvist, FixLaspeyres, FixPaasche, FixFisher, ChainLaspeyres, ChainPaasche, ChainFisher

end # module


text:
using Languages
using Test
using JSON
using InteractiveUtils

@testset "Languages Tests" begin

lang = Languages.English()

@testset "wordlists $lang" begin
  @test all(articles(lang) .== ["a", "an", "the"])
  @test all(indefinite_articles(lang) .== ["a", "an"])
  @test all(definite_articles(lang) .== ["the"])
  @test length(prepositions(lang)) == 33
  @test length(pronouns(lang)) == 36
  @test length(stopwords(lang)) == 488
end

lang = Languages.Romanian()

@testset "wordlists $lang" begin
  @test all(articles(lang) .== ["a", "aceea", "aceia", "acel", "ai", "ale",
                                "aoalea", "cea", "cei", "cel", "cele",
                                "lui", "nişte", "o", "un", "unei", "unor",
                                "unui", "unul"])
  @test all(indefinite_articles(lang) .== ["lui", "nişte", "o", "un", "unei",
                                           "unor", "unui", "unul"])
  # Test that all indefinite articles are present in the articles
  @test all(map(article->in(article, articles(lang)),
                indefinite_articles(lang)))
  @test isempty(definite_articles(lang))
  @test length(prepositions(lang)) == 19
  @test length(pronouns(lang)) == 257
  @test length(stopwords(lang)) == 101
end

lang = Languages.French()

@testset "wordlists $lang" begin
  @test all(articles(lang) .== ["au", "aux", "de l'", "de la", "des",
                                "du", "la", "le", "les", "un", "une",
                                "à la"])
  @test all(indefinite_articles(lang) .== ["des", "un", "une"])
  # Test that all definite, indefinite articles are present in the articles
  @test all(definite_articles(lang) .== ["la", "le", "les"])
  @test all(map(article->in(article, articles(lang)),
                definite_articles(lang)))
  @test all(map(article->in(article, articles(lang)),
                indefinite_articles(lang)))
  @test length(prepositions(lang)) == 258
  @test length(pronouns(lang)) == 122
  @test length(stopwords(lang)) == 665
end

lang = Languages.Italian()

@testset "wordlists $lang" begin
  @test all(articles(lang) .== ["gli", "i", "il", "l'", "la", "le",
                                "lo", "un", "un'", "una", "uno"])
  @test all(indefinite_articles(lang) .== ["un", "un'", "una", "uno"])
  # Test that all definite, indefinite articles are present in the articles
  @test all(definite_articles(lang) .== ["gli", "i", "il", "l'", "la",
                                         "le", "lo"])
  @test all(map(article->in(article, articles(lang)),
                definite_articles(lang)))
  @test all(map(article->in(article, articles(lang)),
                indefinite_articles(lang)))
  @test length(prepositions(lang)) == 45
  @test length(pronouns(lang)) == 27
  @test length(stopwords(lang)) == 642
end

@testset "lang code" begin
    @test Languages.from_code("ben") == Languages.Bengali()
    @test Languages.from_code("Ben") == Languages.Bengali()
    @test Languages.from_code("BEn") == Languages.Bengali()

    @test Languages.from_code("abc") == nothing
end

@testset "lang detect" begin
    include("whatlang.jl")
end

@testset "list languages" begin
    @test list_languages() isa Array
    @test length(list_languages()) == length(subtypes(Language))
end

end


text:
using ODLCommonTools
#using Base.Test
#using FactCheck
using Base.Test
using ArrayViews
import ArrayViews.view

function not_isapprox(args...; kwargs...)
  return !isapprox(args...; kwargs...)
end

@testset "----- Testing Interface and Boundary -----" begin
  # write your own tests here
  @test ( 1 )== 1

  b1 = Boundary(1, 2)
  b2 = Boundary(2, 3)

  @test ( b1 < b2 )== true

  @test ( b2 > b1 )== true

  @test ( getElementL(b1) )== 1
  @test ( getElementL(b2) )== 2
  @test ( getFaceL(b1) )== 2

  b2 = Boundary(1, 3)

  @test ( b1 < b2 )== false
  @test ( b2 > b1 )== false

  b2 = Boundary(1,2)

  @test ( (b1 < b2) )== false
  @test ( (b2 > b1) )== false

  i1 = Interface(1, 2, 3, 4, 5)
  i2 = Interface(2, 3, 4, 5, 6)

  @test ( (i1 < i2) )== true
  @test ( (i2 > i1) )== true

  @test ( getElementL(i1) )== 1
  @test ( getFaceL(i1) )== 3

end

include("test_tools.jl")
include("test_copy.jl")
include("test_roview.jl")
include("test_io.jl")
include("test_lapack.jl")
include("test_null_array.jl")

#FactCheck.exitstatus()


text:
const RTR = ReverseDiff.TrackedReal
const RTV = ReverseDiff.TrackedVector
const RTM = ReverseDiff.TrackedMatrix
using ReverseDiff: record_mul
using ReverseDiff: SpecialInstruction

_eps(::Type{<:RTR{T}}) where {T} = _eps(T)

Base.:*(A::RTM{<:Any, D}, B::AbstractTriangular) where {D} = record_mul(A, B, D)
Base.:*(A::AbstractTriangular, B::RTV{<:Any, D}) where {D} = record_mul(A, B, D)
Base.:*(A::AbstractTriangular, B::RTM{<:Any, D}) where {D} = record_mul(A, B, D)
function Base.:*(A::Adjoint{T, <:AbstractTriangular{T}}, B::RTM{<:Any, D}) where {T, D}
    return record_mul(A, B, D)
end
function Base.:*(A::Adjoint{T, <:AbstractTriangular{T}}, B::RTV{<:Any, D}) where {T, D}
    return record_mul(A, B, D)
end

function replace_diag(::typeof(log), X::RTM{<:Any, D}) where {D}
    tp = ReverseDiff.tape(X)
    X_value = ReverseDiff.value(X)
    f(i, j) = i == j ? log(X_value[i, j]) : X_value[i, j]
    out_value = f.(1:size(X_value, 1), (1:size(X_value, 2))')
    function back(∇)
        g(i, j) = i == j ? ∇[i, j]/X_value[i, j] : ∇[i, j]
        return g.(1:size(X_value, 1), (1:size(X_value, 2))')
    end
    out = ReverseDiff.track(out_value, D, tp)
    ReverseDiff.record!(tp, ReverseDiff.SpecialInstruction, replace_diag, (log, X), out, (back,))
    return out
end

function replace_diag(::typeof(exp), X::RTM{<:Any, D}) where {D}
    tp = ReverseDiff.tape(X)
    X_value = ReverseDiff.value(X)
    f(i, j) = i == j ? exp(X_value[i, j]) : X_value[i, j]
    out_value = f.(1:size(X_value, 1), (1:size(X_value, 2))')
    function back(∇)
        g(i, j) = i == j ? ∇[i, j]*exp(X_value[i, j]) : ∇[i, j]
        return g.(1:size(X_value, 1), (1:size(X_value, 2))')
    end
    out = ReverseDiff.track(out_value, D, tp)
    ReverseDiff.record!(tp, ReverseDiff.SpecialInstruction, replace_diag, (exp, X), out, (back,))
    return out
end

@noinline function ReverseDiff.special_reverse_exec!(instruction::SpecialInstruction{typeof(replace_diag)})
    output = instruction.output
    input = instruction.input
    input_deriv = ReverseDiff.deriv(input[2])
    P = instruction.cache[1]
    input_deriv .+= P(ReverseDiff.deriv(output))
    ReverseDiff.unseed!(output)
    return nothing
end

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(replace_diag)})
    output, input = instruction.output, instruction.input
    out_value = replace_diag(ReverseDiff.value(input[1]), ReverseDiff.value(input[2]))
    ReverseDiff.value!(output, out_value)
    return nothing
end


text:
#####################
#Credal upper_marginal
#####################
#TODO This code seems to assume logspace flows as floating point numbers. if so, enforca that on type F

##################### for local LP ###########################################
using JuMP
using Clp
#using LinearAlgebra

function minimi(coeff :: Array{Float64,1}, l_bounds :: Array{Float64,1}, u_bounds :: Array{Float64,1})
# cosi poi gli entro in broadcasting c, exp.(prob_origin(n).log_thetas), exp.(prob_origin(n).log_thetas_u)
    my_model = Model(Clp.Optimizer)
    set_optimizer_attribute(my_model, "LogLevel", 1)
     set_optimizer_attribute(my_model, "Algorithm", 4)

     @variable(my_model, l_bounds[i] <= x[i = 1:length(l_bounds)]<= u_bounds[i])
     @objective(my_model,Min,coeff' *x) # perché anche * ?
     @constraint(my_model, normalization, ones(Float64, length(l_bounds))'x == 1)

     optimize!(my_model)
        # println("Optimal Solutions:")
        # println("x = ", JuMP.value.(x))
        #  println("Optimal Value:")
        #  println("min =", JuMP.objective_value(my_model))
        #  #getobjectivevalue
     println("funzia")
     return JuMP.objective_value(my_model)
end

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

function credal_marginal_upper_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    resize_flows(circuit, num_examples(data))
    cache = zeros(Float64, num_examples(data)) #TODO: fix type later
    credal_marginal_upper_pass_up_node(n::UpFlowΔNode, ::PlainXData) = ()

    function credal_marginal_upper_pass_up_node(n::UpFlowLiteral{O,F}, cache::Array{Float64}, data::PlainXData{E}) where {E <: eltype(F)} where {O,F}
        pass_up_node(n, data)
        # now override missing values by 1
        npr = pr(n)
        npr[feature_matrix(data)[:,variable(n)] .< zero(eltype(F))] .= 1
        npr .= log.( npr .+ 1e-300 )
        return nothing
    end

    function credal_marginal_upper_pass_up_node(n::UpFlow⋀Cached, cache::Array{Float64}, ::PlainXData)
        pr(n) .= 0
        for i=1:length(n.children)
            # pr(n) .+= pr(n.children[i])
            broadcast!(+, pr(n), pr(n), pr(n.children[i]))
        end
        return nothing
    end

    function credal_marginal_upper_pass_up_node(n::UpFlow⋁Cached, cache::Array{Float64}, ::PlainXData)
        n.pr .= 1e-300
        #pr(n) .= 1e-300

        # for i=1:length(n.children)
        #     cache .= 0
        #     #broadcast reduced memory allocation, though accessing prob_origin(n).log_thetas[i] still allocates lots of extra memory,
        #     #it is proabably due to derefrencing the pointer
        #     broadcast!(+, cache, pr(n.children[i]), prob_origin(n).log_thetas_u[i])
        #     broadcast!(exp, cache, cache)
        #     broadcast!(+, pr(n), pr(n), cache)
        # end

        # broadcast!(log, pr(n), pr(n));
        #u = exp(prob_origin(n).log_thetas_u)
        #u = prob_origin(n).log_thetas_u
        u = Array{Float64}(undef, length(n.children), num_examples(data))
        c = Array{Float64}(undef, length(n.children), num_examples(data))

        # for i=1:length(n.children)
        #     for j=1:num_examples(data)
        #         u[i,j] =  exp(prob_origin(n).log_thetas_u[i])
        #     end
        # end

        ### sotto versione broadcastata

        for j=1:num_examples(data)
            u[:,j] .=  exp.(prob_origin(n).log_thetas_u)
        end

        #println(u)
        # for i=1:length(n.children)
        #     for j=1:num_examples(data)

        #      c[i,j] = exp(n.children[i].pr[j]) ## TODO log is monotone + broadcast
        #     end
        # end

        ### sotto versione broadcastata

        for j=1:num_examples(data)
            c[:,j] .= exp.(pr.(n.children)[j]) ## TODO log is monotone + broadcast
        end

        # for i=1:length(n.children)
        #     println(n.children[i].pr)
        # end
        # println(c)

        optx = Array{Float64}(undef, length(n.children), num_examples(data))
        optx = copy(u)

        for j=1:num_examples(data)

            for i=1:length(n.children)
                if c[i,j] == sort(c[:,j])[1]
                    optx[i,j] = 0.0
                    somma = sum(optx[:,j])
                    optx[i,j] = 1.0-somma ## TODO compatta questo pezzo con find o findall
                    break
                end
            end
        end

        # println(c)
        # println(u)
        # println(optx[:,1])
        # println(sum(optx[:,1]))

        for j=1:num_examples(data)
            pr(n)[j] = log(optx[:,j]'c[:,j])
        end
        return nothing

    end

    ## Pass Up on every node in order
    for n in circuit
        credal_marginal_upper_pass_up_node(n, cache, data)
    end
    return nothing
end

##########################################################################
##Credal lower marginalization
##########################################################################
function credal_marginal_lower_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    resize_flows(circuit, num_examples(data))
    cache = zeros(Float64, num_examples(data)) #TODO: fix type later
    credal_marginal_lower_pass_up_node(n::UpFlowΔNode, ::PlainXData) = ()

    function credal_marginal_lower_pass_up_node(n::UpFlowLiteral{O,F}, cache::Array{Float64}, data::PlainXData{E}) where {E <: eltype(F)} where {O,F}
        pass_up_node(n, data)
        # now override missing values by 1
        npr = pr(n)
        npr[feature_matrix(data)[:,variable(n)] .< zero(eltype(F))] .= 1
        npr .= log.( npr .+ 1e-300 )
        return nothing
    end

    function credal_marginal_lower_pass_up_node(n::UpFlow⋀Cached, cache::Array{Float64}, ::PlainXData)
        pr(n) .= 0
        for i=1:length(n.children)
            # pr(n) .+= pr(n.children[i])
            broadcast!(+, pr(n), pr(n), pr(n.children[i]))
        end
        return nothing
    end

    function credal_marginal_lower_pass_up_node(n::UpFlow⋁Cached, cache::Array{Float64}, ::PlainXData)
        pr(n) .= 1e-300
        l = Array{Float64}(undef, length(n.children), num_examples(data))
        c = Array{Float64}(undef, length(n.children), num_examples(data))

        # for i=1:length(n.children)
        #     for j=1:num_examples(data)
        #         l[i,j] =  exp(prob_origin(n).log_thetas[i])
        #     end
        # end

        for j=1:num_examples(data)
            l[:,j] .=  exp.(prob_origin(n).log_thetas)
        end

        # for i=1:length(n.children)
        #     for j=1:num_examples(data)

        #      c[i,j] = exp(n.children[i].pr[j])
        #     end
        # end

        for j=1:num_examples(data)
            c[:,j] .= exp.(pr.(n.children)[j]) ## TODO log is monotone + broadcast
        end

        optx = Array{Float64}(undef, length(n.children), num_examples(data))
        optx = copy(l)

        for j=1:num_examples(data)

            for i=1:length(n.children)
                if c[i,j] == sort(c[:,j])[1]
                    optx[i,j] = 0.0
                    somma = sum(optx[:,j])
                    optx[i,j] = 1.0-somma ## TODO compatta questo pezzo con find o findall
                    break
                end
            end
        end

        #println(sum(optx[:,1])) this is not exactly 1... approx stuff? or is there something wrong in the code?

        for j=1:num_examples(data)
            pr(n)[j] = log(optx[:,j]'c[:,j])
        end
        return nothing

    end

    # function credal_marginal_lower_pass_up_node(n::UpFlow⋁Cached, cache::Array{Float64}, ::PlainXData)
    #     pr(n) .= 1e-300
    #     for i=1:length(n.children)
    #         cache .= 0
    #         # broadcast reduced memory allocation, though accessing prob_origin(n).log_thetas[i] still allocates lots of extra memory,
    #         # it is proabably due to derefrencing the pointer
    #         broadcast!(+, cache, pr(n.children[i]), prob_origin(n).log_thetas_u[i])
    #         broadcast!(exp, cache, cache)
    #         broadcast!(+, pr(n), pr(n), cache)
    #     end
    #     broadcast!(log, pr(n), pr(n));
    #     return nothing
    # end

    ## Pass Up on every node in order
    for n in circuit
        credal_marginal_lower_pass_up_node(n, cache, data)
    end
    return nothing
end

##### marginal_pass_down

function marginal_pass_down(circuit::DownFlowΔ{O,F}) where {O,F}
    resize_flows(circuit, flow_length(origin(circuit)))
    for n in circuit
        reset_downflow_in_progress(n)
    end
    for downflow in downflow_sinks(circuit[end])
        # initialize root flows to 1
        downflow.downflow .= one(eltype(F))
    end
    for n in Iterators.reverse(circuit)
        marginal_pass_down_node(n)
    end
end

marginal_pass_down_node(n::DownFlowΔNode) = () # do nothing
marginal_pass_down_node(n::DownFlowLeaf) = ()

function marginal_pass_down_node(n::DownFlow⋀Cached)
    # todo(pashak) might need some changes, not tested, also to convert to logexpsum later
     # downflow(n) = EF_n(e), the EF for edges or leaves are note stored
    for c in n.children
        for sink in downflow_sinks(c)
            if !sink.in_progress
                sink.downflow .= downflow(n)
                sink.in_progress = true
            else
                sink.downflow .+= downflow(n)
            end
        end
    end
end

function marginal_pass_down_node(n::DownFlow⋁Cached)
    # todo(pashak) might need some changes, not tested, also to convert to logexpsum later
    # downflow(n) = EF_n(e), the EF for edges or leaves are note stored
    for (ind, c) in enumerate(n.children)
        for sink in downflow_sinks(c)
            if !sink.in_progress
                sink.downflow .= downflow(n) .* exp.(prob_origin(n).log_thetas[ind] .+ pr(origin(c)) .- pr(origin(n)) )
                sink.in_progress = true
            else
                sink.downflow .+= downflow(n) .* exp.(prob_origin(n).log_thetas[ind] .+ pr(origin(c)) .- pr(origin(n)))
            end
        end
    end
end

#### marginal_pass_up_down

function marginal_pass_up_down(circuit::DownFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    @assert !(E isa Bool)
    # marginal_pass_up(origin(circuit), data)
    marginal_pass_down(circuit)
end

"""
Conditional inference Upper bound
"""
function conditional_upper_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    resize_flows(circuit, num_examples(data))

    conditional_upper_pass_up_node(n::UpFlowΔNode, ::PlainXData) = ()

    # pr() now returns a vector of triples (a,b,c), where:
    #  * a:Float64 is the minValue,
    #  * b:Float64 is the maxValue,
    #  * c:Int8 is the type of sub-circuit (to check if contains query variables)

    function conditional_upper_pass_up_node(n::UpFlowLiteral{O,F}, data::PlainXData{E}) where {E <: eltype(F)} where {O,F}
        npr = pr(n)
        nsize = length(npr)
        val = zeros(nsize)
        # Type 0 observed, 1 query, 2 marginalize
        type = zeros( nsize)

        #copy values from data input
        #values form data input
        #   *1,0 observed variable values;
        #   *2,3 query variable values, where  2:0 and 3:1
        #   *-1 unobserved variables (marginalize variables)
        if positive(n)
            val = feature_matrix(data)[:,variable(n)]
        else
            val =  ones(Float64,nsize) - feature_matrix(data)[:,variable(n)]
        end

        # override query variable values by 0 and 1
        val[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 0

        val[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1

        type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+2.0)] .= 1

        type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+3.0)] .= 1

        #  override marginalize variable  by 1
        val[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1
        type[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 2

        # Construct the tuple to  propagate 4 values : conditional message, marginal upper message, marginal lower message, type
        pr(n) .= [x for x in zip(val,val,val,type)]
        return nothing
    end

    function conditional_upper_pass_up_node(n::UpFlow⋀Cached, ::PlainXData)
        npr = pr(n)
        #Initialize values for tuples
        val_min = ones( length(npr))
        val_max = ones( length(npr))
        type = zeros( length(npr))

        for i=1:length(n.children)
            #Obtain values from child i -> Message
            val_min_i = getindex.(pr(n.children[i]),1) # vector a: minValue len(a)=number of queries
            val_max_i = getindex.(pr(n.children[i]),2) # vector b: maxValue
            type_i = getindex.(pr(n.children[i]),3)    # vector c: type of subcircuit

            #...
        end

        pr(n) .= [x for x in zip(val_min,val_max,val_max,type)]
        return nothing
    end

    function conditional_upper_pass_up_node(n::UpFlow⋁Cached,  ::PlainXData)
        npr = pr(n)
        val_min = zeros( length(npr))
        val_max = zeros( length(npr))
        type = zeros( length(npr))
        for i=1:length(n.children)
            val_min_i = getindex.(pr(n.children[i]),1)
            val_max_i = getindex.(pr(n.children[i]),2)
            type_i = getindex.(pr(n.children[i]),3)

            #.....

        end
        pr(n) .= [x for x in zip(val_min,val_max,val_max,type)]
        return nothing
    end

    for n in circuit
        conditional_upper_pass_up_node(n,data)
    end
    return nothing
end

"""
Conditional inference lower bound
"""
function conditional_lower_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}, mu::Array{Float64,1}) where {E <: eltype(F)} where {O,F}
    resize_flows(circuit, num_examples(data))

    conditional_lower_pass_up_node(n::UpFlowΔNode, ::PlainXData, mu::Array{Float64,1}) = ()

    # pr() now returns a vector of triples (a,b,c), where:
    #  * a:Float64 is the minValue,
    #  * b:Float64 is the maxValue,
    #  * c:Int8 is the type of sub-circuit (to check if contains query variables)

    function conditional_lower_pass_up_node(n::UpFlowLiteral{O,F}, data::PlainXData{E}, mu::Array{Float64,1}) where {E <: eltype(F)} where {O,F}

        println("LITERAL NODE")

        npr = pr(n)
        #mu = 1.0 #  se metto 0.6 mi da ERROR perché ???   # dopo mu::Array{Float64,1}
        nsize = length(npr)
        val = zeros(nsize)           # primo el tuple
        marginal_lower = zeros(nsize) # secondo el tuple
        marginal_upper = zeros(nsize) # terzo el tuple
        # Type 0 observed, 1 query, 2 marginalize
        type = zeros( nsize)

        ################ MARGINAL MESSAGES ################################################
        ### qui preparo il messaggio marginal_lower
        #####pass_up_node(n, data) questo non funge
        # now override missing values by 1
        #######npr_marginal_lower = pr(n)
        #######npr_marginal_lower[feature_matrix(data)[:,variable(n)] .< zero(eltype(F))] .= 1
       ######## npr_marginal_lower .= log.( npr_marginal_lower .+ 1e-300 )

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

        #copy values from data input
        #values form data input
        #   *1,0 observed variable values;
        #   *2,3 query variable values, where  2:0 and 3:1
        #   *-1 unobserved variables (marginalize variables)
        if positive(n) ### qui non ci sono nodi top ricorda che i top diventano or !!! quindi n è o pos o neg

            marginal_lower = feature_matrix(data)[:,variable(n)] #marginal_lower
            marginal_upper = feature_matrix(data)[:,variable(n)] #marginal_lower
            val = feature_matrix(data)[:,variable(n)]
        else

            marginal_lower =  ones(Float64,nsize) - feature_matrix(data)[:,variable(n)] ### TODO in verità se non è queried posso mettere un valore qualsiasi come per quelle da marginalizzare, vedi pseudocode
            marginal_upper =  ones(Float64,nsize) - feature_matrix(data)[:,variable(n)] ### TODO in verità se non è queried posso mettere un valore qualsiasi come per quelle da marginalizzare, vedi pseudocode
            val =  ones(Float64,nsize) - feature_matrix(data)[:,variable(n)] ### TODO in verità se non è queried posso mettere un valore qualsiasi come per quelle da marginalizzare, vedi pseudocode
        end ### sistemato il caso in cui la variabile sia una var osservat

            # override query variable values by 0 and 1, depending both on the queried status and the literal status
        if positive(n)
            #println("LITERAL POSITIVE NODE")

            val[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= (zeros(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)]

            val[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= (ones(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)]

            marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 0.0
            marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 0.0

            marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0
            marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 1.0

        else

            #println("LITERAL NEGATIVE NODE")

            val[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= (ones(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)]

            val[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= (zeros(Float64,nsize) .- mu)[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)]

            marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0
            marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 2.0)] .= 1.0

            marginal_lower[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 0.0
            marginal_upper[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize) .+ 3.0)] .= 0.0
        end

         # type of queried variables (independent both from the query status and the literal status, correct)

         type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+2.0)] .= 1

         type[feature_matrix(data)[:,variable(n)] .== (zeros(Float64,nsize).+3.0)] .= 1

         #  override marginalize variable  by 1, independently from queried status and literal status
         val[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1
         marginal_lower[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1
         marginal_upper[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 1

         # type of variables to marginalize
         type[feature_matrix(data)[:,variable(n)] .< zeros(Float64,nsize)] .= 2

         # TODO sistemare i marginal per la questione log

         #marginal_lower .= log.(marginal_lower .+ 1e-300 )
         #marginal_upper .= log.(marginal_upper .+ 1e-300 )

         # Construct the tuple to  propagate 4 values : conditional message, marginal lower message, marginal upper message, type
         pr(n) .= [x for x in zip(val,marginal_lower,marginal_upper,type)]
         println("tupla literal : ", pr(n))

         return nothing
    end

    function conditional_lower_pass_up_node(n::UpFlow⋀Cached, ::PlainXData, mu::Array{Float64,1})

        println("AND NODE")

    ############ MARGINAL STUFF ############
    # function credal_marginal_lower_pass_up_node(n::UpFlow⋀Cached, cache::Array{Float64}, ::PlainXData)
    #     pr(n) .= 0
    #     for i=1:length(n.children)
    #         # pr(n) .+= pr(n.children[i])
    #         broadcast!(+, pr(n), pr(n), pr(n.children[i]))
    #     end
    #     return nothing
    # end

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

        npr = pr(n)
        #Initialize values for tuples
        val_min = ones( length(npr))
        marginal_lower = zeros( length(npr))
        marginal_upper = ones( length(npr))
        type = zeros( length(npr))
        vect = ones( length(npr))

        for i=1:length(n.children) # che per noi è 2 perché siamo in un AND node

            #### Marginal messages: TODO work with their log
            broadcast!(*, marginal_lower, marginal_lower, getindex.(pr(n.children[i]),2))
            broadcast!(*, marginal_upper, marginal_upper, getindex.(pr(n.children[i]),3))

            ##### conditional message (val_min) is the product of pi_lower and marginal_lower/upper of the children, depending of the types of the latter
            # val_min for now has all entries set to 1
            # start by giving pi_lower values when the type of the children i (in which we are in the for cycle) is 1
            val_min[getindex.(pr(n.children[i]),4) .== ones( length(npr))] .= getindex.(pr(n.children[i]),1)[getindex.(pr(n.children[i]),4) .== ones( length(npr))]

            ##### type. as soon as I see a children with type 1, set type = 1. otherwise it stays 0. IN AND/OR NODES TYPE IS ONLY 1 OR 0 !
            type[getindex.(pr(n.children[i]),4) .== ones( length(npr))] .= 1

        end

        # multiply val_min by the correct marginal
        val_min[(getindex.(pr(n.children[1]),4) .!== ones( length(npr))) .& (getindex.(pr(n.children[2]),1) .> zeros( length(npr))) ].*= getindex.(pr(n.children[1]),2)[(getindex.(pr(n.children[1]),4) .!== ones( length(npr))) .& (getindex.(pr(n.children[2]),1) .> zeros( length(npr))) ]
        val_min[(getindex.(pr(n.children[1]),4) .!== ones( length(npr))) .& (getindex.(pr(n.children[2]),1) .< zeros( length(npr))) ].*= getindex.(pr(n.children[1]),3)[(getindex.(pr(n.children[1]),4) .!== ones( length(npr))) .& (getindex.(pr(n.children[2]),1) .< zeros( length(npr))) ]

        val_min[(getindex.(pr(n.children[2]),4) .!== ones( length(npr))) .& (getindex.(pr(n.children[1]),1) .> zeros( length(npr))) ].*= getindex.(pr(n.children[2]),2)[(getindex.(pr(n.children[2]),4) .!== ones( length(npr))) .& (getindex.(pr(n.children[1]),1) .> zeros( length(npr))) ]
        val_min[(getindex.(pr(n.children[2]),4) .!== ones( length(npr))) .& (getindex.(pr(n.children[1]),1) .< zeros( length(npr))) ].*= getindex.(pr(n.children[2]),3)[(getindex.(pr(n.children[2]),4) .!== ones( length(npr))) .& (getindex.(pr(n.children[1]),1) .< zeros( length(npr))) ]

        pr(n) .= [x for x in zip(val_min,marginal_lower,marginal_upper,type)]
        return nothing
    end

    function conditional_lower_pass_up_node(n::UpFlow⋁Cached,  ::PlainXData, mu::Array{Float64,1})
        println("OR NODE")

        npr = pr(n)
        val_min = zeros( length(npr))
        marginal_lower = zeros( length(npr))
        marginal_upper = zeros( length(npr))
        type = zeros( length(npr))

        if n.children[1] isa UpFlowLiteral

            coeff = Array{Float64}(undef, length(npr), length(n.children)) # length(n.children) dev'essere sempre 2, assert?
            println("val_min_prima = ", val_min)

            for i=1:2

                # println("positivity of children ", i, positive(n.children[i]))
                # println("tripletta i-esimo children, per ogni instance, con i= ", i, " : " ,   pr(n.children[i]))
                # println("lower bounds : ", exp.(prob_origin(n).log_thetas)) ### questi sono gli stessi per ogni istanza
                # println("upper bounds : ", exp.(prob_origin(n).log_thetas_u))

                coeff[:,i] .= getindex.(pr(n.children[i]),1)

            end

            for i=1:length(npr) ### TODO broadcast this
                val_min[i] = minimi(coeff[i,:], exp.(prob_origin(n).log_thetas), exp.(prob_origin(n).log_thetas_u))

            end

             println("val_min_dopo = ", val_min)
             println("IT'S A TOP !")
        else

             # for j=1:length(n.children)
            #     val_min_j = getindex.(pr(n.children[j]),1)
            #     val_max_j = getindex.(pr(n.children[j]),2)
            #     type_j = getindex.(pr(n.children[j]),3) # è un vettore in cui la k-esima entry mi dice il tipo del i-esimo child wrt la k-esima istanza
            #     left_index = zeros( length(npr))
            #     right_index = zeros( length(npr))

            #     # deduco il type di n dal type dei suoi figli. ridondante rifarlo per ogni children
            #     # visto che tutti i type_j sono uguali (sono i tipi degli elementi del nodo or) ma vabbè

            #     type[type_j .== ones(Int8,length(npr))] .= 1

            #     # dato il j-esimo children, mi memorizzo, per ogni istanza, se la queried variable è una left var o right var.
            #     # se la queried var è in n (ie, se la i-esima entry di type di n è uguale a 1), left_index e right_index devono sommare a 1 in quella entry.
            #     # altrimenti la queried var non è né left né right (ie, la i-esima entry di type di n è uguale a 0).
            #     # In tal caso ( e solo in tal caso) left_index e right_index devono sommare a 0 in quella entry. condizione da usare dopo per settare npr !!! :)
            #     #println("pr del primo grandchildren del j-esimo child")
            #     #println(typeof(pr(n.children[j].children[1])))
            #     #println(pr(n.children[j].children[1]))
            #     #println(pr(n.children[j].children[1])[:,3]) cosi non gli piace ! è un array unidimensionale (ie un vettore) di tuples (triplette qui)

            #     #println(getindex(pr(n.children[j].children[1])[1],3)) cosi ok

            #     left_index[getindex.(pr(n.children[j].children[1]),3) .== 1] .= 1
            #     right_index[getindex.(pr(n.children[j].children[2]),3) .== 1] .= 1  #### TODO assertion for somma minore di uno

            #     ######## OK fino a qua apposto. ora posso creare i coeff per ogni nodo or che non sia un top:

            #     ### comincio a mettere il \pi   QUESTI FUNZIONAVA CON C MATRICE
            #     c[:,j][left_index .== 1.0] .= getindex.(pr(n.children[j].children[1])[left_index .== 1.0],1) ### qua mi calcola i pi_lower
            #     c[:,j][right_index .== 1.0] .= getindex.(pr(n.children[j].children[2])[right_index .== 1.0],1)

            #      ### ora aggiungo (in realtà moltiplico) il \sigma
            #     # c[:,j][c[:,j].< zeros(Float64,length(npr)) e left_index .== 1.0] .*= exp.(log_marginal_upper(subcircuito che parte da n.children[j].children[2], sub-batch che concerne solo le variabili di sto sub-cricuito))
            #     # c[:,j][c[:,j].< zeros(Float64,length(npr)) e right_index .== 1.0] .*= exp.(log_marginal_upper(subcircuito che parte da n.children[j].children[1], sub-batch che concerne solo le variabili di sto sub-cricuito))
            #     # c[:,j][c[:,j].> zeros(Float64,length(npr)) e left_index .== 1.0] .*= exp.(log_marginal_lower(subcircuito che parte da n.children[j].children[2], sub-batch che concerne solo le variabili di sto sub-cricuito))
            #     # c[:,j][c[:,j].> zeros(Float64,length(npr)) e right_index .== 1.0] .*= exp.(log_marginal_lower(subcircuito che parte da n.children[j].children[1], sub-batch che concerne solo le variabili di sto sub-cricuito))

            #    ### vettore = zeros( length(npr))
            #     ###println(vettore)
            #     ###println(c[:,j])
            #     ###c[:,j][(c[:,j].< zeros(Float64,length(npr))) .&& (left_index .== 1.0)] .= vettore[(c[:,j].< zeros(Float64,length(npr))) .&& (left_index .== 1.0)]
            # end

         ##### qui risolvo, in parallelo per ogni instance (query) LP e storo gli ottimi nei vari vettori

         ## conditional

            println("robetta da passare come coeff di tipo : ", typeof(getindex.(pr.(n.children),1)))
            println("tuple dei children : ", pr.(n.children))
            println("tuple della 1a inst  : ", getindex.(pr.(n.children),1))
            println("primi el delle tuple della 1a inst: ", getindex.(getindex.(pr.(n.children),1),1))

            c_co = Array{Float64}(undef, length(npr), length(n.children))
            c_marg_lo = Array{Float64}(undef, length(npr), length(n.children))
            c_marg_up = Array{Float64}(undef, length(npr), length(n.children))

            for i=1:length(npr)
                c_co[i,:] = getindex.(getindex.(pr.(n.children),i),1) # vettore dei coefficenti del i-esimo problema COND  (i-esima instance)
                c_marg_lo[i,:] = getindex.(getindex.(pr.(n.children),i),2) # vettore dei coefficenti del i-esimo problema MARG_LOWER (i-esima instance)
                c_marg_up[i,:] = getindex.(getindex.(pr.(n.children),i),3) # vettore dei coefficenti del i-esimo problema MARG_UPPER  (i-esima instance)
            end

            println("coeff cond  : ", c_co)
            println("coeff marg_lo  : ", c_marg_lo)
            println("coeff marg_up  : ", c_marg_up)

            for i=1:length(npr)
                val_min[i] = minimi(c_co[i,:], exp.(prob_origin(n).log_thetas), exp.(prob_origin(n).log_thetas) )
                marginal_lower[i] = minimi(c_marg_lo[i,:], exp.(prob_origin(n).log_thetas), exp.(prob_origin(n).log_thetas) )
                marginal_upper[i] = minimi(c_marg_up[i,:], exp.(prob_origin(n).log_thetas), exp.(prob_origin(n).log_thetas) )
            end

        end

        pr(n) .= [x for x in zip(val_min,marginal_lower,marginal_upper,type)]
        return nothing
    end

    for n in circuit
        conditional_lower_pass_up_node(n,data,mu)
       #println( "tupla : ", pr(n)) #non va
    end
    return nothing
end


text:
abstract type Day7 <: AocDay end
day7 = Day7

using LightGraphs
using MetaGraphs
using GraphPlot
using StatsBase
using SimpleWeightedGraphs

function readinput(day::Type{Day7})
    rules = readlines("input/day7.txt")
    parsed_rules = parse_rule.(rules)
    parsed_rules
end

function parse_rule(rule)
    # extract container
    re = r"^(?<container>.*) bags contain"
    m = match(re, rule)

    # extract contents
    container = m["container"]
    contents = replace(rule, container * " bags contain " => "")
    contents = contents[1:end-1] # strip dot at the end
    contents = split(contents, ", ")
    contents = [replace(c, r" bags?" => "") for c in contents]
    labels = [string(container)]
    if contents != ["no other"]
        for bagtype in contents
            m = match(r"^(?<num>\d+) (?<bt>.*)$", bagtype)
            for i in 1:parse(Int, m["num"])
                push!(labels, string(m["bt"]))
            end
        end
    end

    labels
end

get_nodelabels_from_idv(idv) = [string(i) * " " * idv[2][i] for i in sort(collect(keys(idv[2])))]

function get_graph(r)
    vertices = first.(r);
    idv = [
        Dict(j => i for (i, j) in enumerate(vertices)),
        Dict(i => j for (i, j) in enumerate(vertices))
    ]
    bag_contents = Dict(i[1] => countmap(i[2:end]) for i in r)
    edges = Dict(i[1] => countmap(i[2:end]) for i in r)
    g = SimpleWeightedDiGraph(length(vertices))
    for (k, v) in bag_contents
        src_id = idv[1][k]
        for (kk, vv) in v
            dst_id = idv[1][kk]
            weight = vv
            for w in 1:weight
                add_edge!(g, src_id, dst_id, w)
            end
        end
    end
    (idv, g)
end

function solve(day::Type{Day7})
    r = readinput(day)
    idv, g = get_graph(r)
    nodelabels = get_nodelabels_from_idv(idv)
    starting_vertex = idv[1]["shiny gold"]
    ps = []
    function traverse_until_outer(g, vertex, path)
        parents = inneighbors(g, vertex)
        for p in parents
            push!(ps, p)
        end
        length(parents) == 0 && return path
        [traverse_until_outer(g, parent, [path; parent]) for parent in parents]
    end
    traverse_until_outer(g, starting_vertex, [])
    ps |> unique |> length
end

function solve_part2(day::Type{Day7})
    r = readinput(day)
    idv, g = get_graph(r)
    nodelabels = get_nodelabels_from_idv(idv)

    starting_vertex = idv[1]["shiny gold"]
    cs = []
    function traverse_in(g, vertex, path)
        @show path
        children = outneighbors(g, vertex)
        for c in children
            w = g.weights[c, vertex]
            #@show [c, w]
            push!(cs, [length(path), c, w])
        end
        length(children) == 0 && return path
        [traverse_in(g, child, [path; child]) for child in children]
    end
    traverse_in(g, starting_vertex, [starting_vertex])
    s = 0
    for (idx, (level, v, w)) in enumerate(cs)
        #@show idx, level, v, w
        s += w^level
    end
    return r, idv, g, s

end