text:
@testset "rbs.jl" begin
    @testset "expandRbs!()" begin
        ### Checking status ###
        # :NodeLimitStop

        search=SeaPearl.staticRBSearch{SeaPearl.InfeasibleNodeCriteria}(10, 10)
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.numberOfNodes = 1
        toCall = Stack{Function}()
        @test SeaPearl.expandRbs!(toCall, model, 1, search, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :NodeLimitStop
        @test isempty(toCall)

        # :SolutionLimitStop
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.numberOfSolutions = 0

        toCall = Stack{Function}()
        @test SeaPearl.expandRbs!(toCall, model, 1, search, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :SolutionLimitStop
        @test isempty(toCall)

        #:TimeLimitStop
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.searchingTime = 0

        toCall = Stack{Function}()
        SeaPearl.tic()
        @test SeaPearl.expandRbs!(toCall, model, 1, search, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic())  == :TimeLimitStop
        @test isempty(toCall)

        # :Infeasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(3, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.expandRbs!(toCall, model, 1, search, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :Infeasible
        @test isempty(toCall)

        # :Feasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 2, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.expandRbs!(toCall, model, 1, search, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :FoundSolution
        @test isempty(toCall)

        ### Checking stack ###
        search=SeaPearl.staticRBSearch{SeaPearl.InfeasibleNodeCriteria}(10, 10)

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.expandRbs!(toCall, model, 10, search, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :Feasible
        @test length(toCall) == 6

        @test pop!(toCall)(model) == :SavingState
        @test length(model.trailer.prior) == 1 # saveState!()

        @test pop!(toCall)(model) == :FoundSolution
        @test length(model.statistics.solutions) == 1 # Found a solution

        @test pop!(toCall)(model) == :BackTracking
        @test length(model.trailer.prior) == 0 # restoreState!()

        @test pop!(toCall)(model) == :SavingState
        @test length(model.trailer.prior) == 1 # saveState!()

        @test pop!(toCall)(model) == :FoundSolution
        @test length(model.statistics.solutions) == 2 # Found another solution

        @test pop!(toCall)(model) == :BackTracking
        @test length(model.trailer.prior) == 0 # restoreState!()
    end

    @testset "initroot(::staticRBSearch)" begin
        # :Feasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 2, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.initroot!(toCall, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10, 5), model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :FoundSolution
        @test length(toCall) == 4
        @test pop!(toCall).nodeLimitList[2] == 10
        @test pop!(toCall).nodeLimitList[3] == 10
        @test pop!(toCall).nodeLimitList[4] == 10
        @test pop!(toCall).nodeLimitList[5] == 10

    end

    @testset "initroot(::geometricRBSearch)" begin
        # :Feasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 2, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.initroot!(toCall, SeaPearl.geometricRBSearch{SeaPearl.VisitedNodeCriteria}(10, 5, 1.1), model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic())== :FoundSolution
        @test length(toCall) == 4
        @test pop!(toCall).nodeLimitList[2] == 11
        @test pop!(toCall).nodeLimitList[3] == 13
        @test pop!(toCall).nodeLimitList[4] == 14
        @test pop!(toCall).nodeLimitList[5] == 15

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 2, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.initroot!(toCall, SeaPearl.geometricRBSearch{SeaPearl.VisitedNodeCriteria}(10, 1, 1.1), model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic())== :FoundSolution
        @test isempty(toCall)
    end

    @testset "initroot(::lubyRBSearch)" begin   # The Luby sequence is a sequence of the following form: 1,1,2,1,1,2,4,1,1,2,1,1,2,4,8, . .
        # :Feasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 2, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        toCall = Stack{Function}()
        @test SeaPearl.initroot!(toCall, SeaPearl.lubyRBSearch{SeaPearl.VisitedNodeCriteria}(10,5), model, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic())== :FoundSolution
        @test length(toCall) == 4
        @test pop!(toCall).nodeLimitList[2] == 10
        @test pop!(toCall).nodeLimitList[3] == 20
        @test pop!(toCall).nodeLimitList[4] == 10
        @test pop!(toCall).nodeLimitList[5] == 10

        @test SeaPearl.generateLimitList(SeaPearl.lubyRBSearch{SeaPearl.VisitedNodeCriteria}(10,31))==[10,10,20,10,10,20,40,10,10,20,10,10,20,40,80,10,10,20,10,10,20,40,10,10,20,10,10,20,40,80,160]
    end

    @testset "Stopping Criteria" begin
        @testset "InfeasibleNodeCriteria" begin

            criteria = SeaPearl.staticRBSearch{SeaPearl.InfeasibleNodeCriteria}(1, 0)

            trailer = SeaPearl.Trailer()
            model = SeaPearl.CPModel(trailer)

            x = SeaPearl.IntVar(2, 2, "x", trailer)
            y = SeaPearl.IntVar(3, 3, "y", trailer)
            SeaPearl.addVariable!(model, x)
            SeaPearl.addVariable!(model, y)
            SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

            toCall = Stack{Function}()
            @test SeaPearl.expandRbs!(toCall, model, 1, criteria, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :Infeasible
            @test criteria(model,1) == false #one infeasible state has been reached, so the search stops.
        end
        @testset "VisitedNodeCriteria" begin

            criteria = SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(1, 0)

            trailer = SeaPearl.Trailer()
            model = SeaPearl.CPModel(trailer)

            x = SeaPearl.IntVar(2, 3, "x", trailer)
            y = SeaPearl.IntVar(2, 3, "y", trailer)
            SeaPearl.addVariable!(model, x)
            SeaPearl.addVariable!(model, y)
            SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

            toCall = Stack{Function}()
            @test SeaPearl.expandRbs!(toCall, model, 1, criteria, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :Feasible
            @test criteria(model,1) == true    #only one node has been visited
            pop!(toCall)(model) # :SavingState
            pop!(toCall)(model) # expandRbs!(...)
            @test criteria(model,1) == false
        end

        @testset "SolutionFoundCriteria" begin

            criteria = SeaPearl.staticRBSearch{SeaPearl.SolutionFoundCriteria}(2, 0)

            trailer = SeaPearl.Trailer()
            model = SeaPearl.CPModel(trailer)

            x = SeaPearl.IntVar(2, 3, "x", trailer)
            y = SeaPearl.IntVar(2, 3, "y", trailer)
            SeaPearl.addVariable!(model, x)
            SeaPearl.addVariable!(model, y)
            SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

            toCall = Stack{Function}()
            @test SeaPearl.expandRbs!(toCall, model, 2, criteria, SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :Feasible
            @test criteria(model,2) == true
            pop!(toCall)(model) # :SavingState
            @test pop!(toCall)(model) == :FoundSolution
            @test criteria(model,2) == true  #One solution found
            pop!(toCall)(model) # :Backtrack
            pop!(toCall)(model) # :SavingState
            @test pop!(toCall)(model) == :FoundSolution
            @test criteria(model,2) == false #Two solution found
        end
    end

    @testset "search!(::RBSearch)" begin
        ### Checking status ###
        # :LimitStop
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.numberOfNodes = 1
        @test SeaPearl.search!(model, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10,10), SeaPearl.MinDomainVariableSelection()) == :NodeLimitStop

        # :SolutionLimitStop
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)
        model.limit.numberOfSolutions = 0

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        @test SeaPearl.search!(model, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10,10), SeaPearl.MinDomainVariableSelection()) == :SolutionLimitStop

        # :Infeasible
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(3, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        @test SeaPearl.search!(model, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10,10), SeaPearl.MinDomainVariableSelection()) == :Infeasible

        # :Optimal
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(2, 2, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        push!(model.constraints, SeaPearl.Equal(x, y, trailer))
        @test SeaPearl.search!(model, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10,10), SeaPearl.MinDomainVariableSelection()) == :NonOptimal
        @test length(model.statistics.solutions) == 1

        ### Checking more complex solutions ###
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        @test SeaPearl.search!(model, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10,10), SeaPearl.MinDomainVariableSelection()) == :NonOptimal
        @test length(model.statistics.solutions) == 2
        @test model.statistics.solutions[1] == Dict("x" => 3,"y" => 3)
        @test model.statistics.solutions[2] == Dict("x" => 2,"y" => 2)

    end

    @testset "search!() with a BasicHeuristic" begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))

        @test SeaPearl.search!(model, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10,10), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic()) == :NonOptimal
        @test model.statistics.solutions[1] == Dict("x" => 3,"y" => 3)
        @test model.statistics.solutions[2] == Dict("x" => 2,"y" => 2)

        SeaPearl.empty!(model)

        x = SeaPearl.IntVar(2, 3, "x", model.trailer)
        y = SeaPearl.IntVar(2, 3, "y", model.trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, model.trailer))

        my_heuristic(x::SeaPearl.IntVar; cpmodel=nothing) = minimum(x.domain)
        @test SeaPearl.search!(model, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10,10), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(my_heuristic)) == :NonOptimal
        @test model.statistics.solutions[1] == Dict("x" => 2,"y" => 2)
        @test model.statistics.solutions[2] == Dict("x" => 3,"y" => 3)

    end

    @testset "search!() with a BasicHeuristic - out_solver being true" begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))
        @test SeaPearl.search!(model, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10,10), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(); out_solver=true) == :FoundSolution
        @test length(model.statistics.solutions) == 1
        @test model.statistics.solutions[1] == Dict("x" => 3,"y" => 3)

        SeaPearl.empty!(model)

        x = SeaPearl.IntVar(2, 3, "x", model.trailer)
        y = SeaPearl.IntVar(2, 3, "y", model.trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, model.trailer))

        my_heuristic(x::SeaPearl.IntVar; cpmodel=nothing) = minimum(x.domain)
        @test SeaPearl.search!(model, SeaPearl.staticRBSearch{SeaPearl.VisitedNodeCriteria}(10,10), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic(my_heuristic); out_solver=true) == :FoundSolution
        @test length(model.statistics.solutions) == 1
        @test model.statistics.solutions[1] == Dict("x" => 2,"y" => 2)

    end

    @testset "search!() with PPO/RBS" begin

    using Random
    numInFeatures = 3
    approximator_GNN = GeometricFlux.GraphConv(64 => 64, Flux.leakyrelu)
    target_approximator_GNN = GeometricFlux.GraphConv(64 => 64, Flux.leakyrelu)
    gnnlayers = 10
    seed = 123
    rng = MersenneTwister(seed)

    actor_approximator = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            GeometricFlux.GraphConv(numInFeatures => 64, Flux.leakyrelu),
            [approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu; initW = glorot_uniform(rng)),
            Flux.Dense(32, 32, Flux.leakyrelu; initW = glorot_uniform(rng)),
            Flux.Dense(32, 16, Flux.leakyrelu; initW = glorot_uniform(rng)),
        ),
        outputChain = Flux.Dense(16, 4; initW = glorot_uniform(rng)),
    )
    critic_approximator = SeaPearl.CPNN(
        graphChain = Flux.Chain(
            GeometricFlux.GraphConv(numInFeatures => 64, Flux.leakyrelu),
            [target_approximator_GNN for i = 1:gnnlayers]...
        ),
        nodeChain = Flux.Chain(
            Flux.Dense(64, 32, Flux.leakyrelu; initW = glorot_uniform(rng)),
            Flux.Dense(32, 32, Flux.leakyrelu; initW = glorot_uniform(rng)),
            Flux.Dense(32, 16, Flux.leakyrelu; initW = glorot_uniform(rng)),
        ),
        outputChain = Flux.Dense(16, 1; initW = glorot_uniform(rng)),
    )

    UPDATE_FREQ = 1

    # Agent definition
    agent = RL.Agent(
        policy = RL.QBasedPolicy(
            learner = RL.A2CLearner(
                approximator =RL.ActorCritic(
                    actor = actor_approximator,
                    critic = critic_approximator,
                    optimizer = ADAM(1e-3),
                ) |> cpu,
                γ = 0.99f0,
                actor_loss_weight = 1.0f0,
                critic_loss_weight = 0.5f0,
                entropy_loss_weight = 0.001f0,
                update_freq = UPDATE_FREQ,

            ),
            explorer = GumbelSoftmaxExplorer(),
            ),
        trajectory = RL.CircularArraySARTTrajectory(
            capacity = 8,
            state = SeaPearl.DefaultTrajectoryState => (),
            action = Vector{Int} => (1,),
            reward = Vector{Float32} => (1,),
            terminal = Vector{Bool} => (1,),
        )
    )

        # define the value selection
        valueSelection = SeaPearl.LearnedHeuristic(agent)

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x1 = SeaPearl.IntVar(1, 2, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 2, "x2", trailer)
        x3 = SeaPearl.IntVar(2, 3, "x3", trailer)
        x4 = SeaPearl.IntVar(1, 4, "x4", trailer)
        SeaPearl.addVariable!(model, x1)
        SeaPearl.addVariable!(model, x2)
        SeaPearl.addVariable!(model, x3)
        SeaPearl.addVariable!(model, x4)

        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x1, x2, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x2, x3, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x3, x4, trailer))

        # define the variable selection
        variableSelection = SeaPearl.MinDomainVariableSelection()

        # launch the search
        SeaPearl.search!(model, SeaPearl.geometricRBSearch{SeaPearl.InfeasibleNodeCriteria}(3,10,1.1), variableSelection, valueSelection)

        possible_solutions = [
            Dict("x1" => 1, "x2" => 2, "x3" => 3, "x4" => 1),
            Dict("x1" => 1, "x2" => 2, "x3" => 3, "x4" => 2),
            Dict("x1" => 1, "x2" => 2, "x3" => 3, "x4" => 4),
            Dict("x1" => 2, "x2" => 1, "x3" => 3, "x4" => 1),
            Dict("x1" => 2, "x2" => 1, "x3" => 3, "x4" => 2),
            Dict("x1" => 2, "x2" => 1, "x3" => 3, "x4" => 4),
            Dict("x1" => 2, "x2" => 1, "x3" => 2, "x4" => 1),
            Dict("x1" => 2, "x2" => 1, "x3" => 2, "x4" => 3),
            Dict("x1" => 2, "x2" => 1, "x3" => 2, "x4" => 4)
        ]

        for solution in model.statistics.solutions
            @test solution in possible_solutions
        end

    end

    @testset "search!() with a PPO/RBS" begin

        using Random
        numInFeatures = 3
        approximator_GNN = GeometricFlux.GraphConv(64 => 64, Flux.leakyrelu)
        target_approximator_GNN = GeometricFlux.GraphConv(64 => 64, Flux.leakyrelu)
        gnnlayers = 10
        seed = 123
        rng = MersenneTwister(seed)

        actor_approximator = SeaPearl.CPNN(
            graphChain = Flux.Chain(
                GeometricFlux.GraphConv(numInFeatures => 64, Flux.leakyrelu),
                [approximator_GNN for i = 1:gnnlayers]...
            ),
            nodeChain = Flux.Chain(
                Flux.Dense(64, 32, Flux.leakyrelu; initW = glorot_uniform(rng)),
                Flux.Dense(32, 32, Flux.leakyrelu; initW = glorot_uniform(rng)),
                Flux.Dense(32, 16, Flux.leakyrelu; initW = glorot_uniform(rng)),
            ),
            outputChain = Flux.Dense(16, 4; initW = glorot_uniform(rng)),
        )
        critic_approximator = SeaPearl.CPNN(
            graphChain = Flux.Chain(
                GeometricFlux.GraphConv(numInFeatures => 64, Flux.leakyrelu),
                [target_approximator_GNN for i = 1:gnnlayers]...
            ),
            nodeChain = Flux.Chain(
                Flux.Dense(64, 32, Flux.leakyrelu; initW = glorot_uniform(rng)),
                Flux.Dense(32, 32, Flux.leakyrelu; initW = glorot_uniform(rng)),
                Flux.Dense(32, 16, Flux.leakyrelu; initW = glorot_uniform(rng)),
            ),
            outputChain = Flux.Dense(16, 1; initW = glorot_uniform(rng)),
        )

        UPDATE_FREQ = 1

        # Agent definition   //     PPO
        agent = RL.Agent(
            policy = RL.QBasedPolicy(
                learner = RL.A2CLearner(
                    approximator =RL.ActorCritic(
                        actor = actor_approximator,
                        critic = critic_approximator,
                        optimizer = ADAM(1e-3),
                    ) |> cpu,
                    γ = 0.99f0,
                    actor_loss_weight = 1.0f0,
                    critic_loss_weight = 0.5f0,
                    entropy_loss_weight = 0.001f0,
                    update_freq = UPDATE_FREQ,

                ),
                explorer = GumbelSoftmaxExplorer(),
                ),
            trajectory = RL.CircularArraySARTTrajectory(
                capacity = 8,
                state = SeaPearl.DefaultTrajectoryState => (),
                action = Vector{Int} => (1,),
                reward = Vector{Float32} => (1,),
                terminal = Vector{Bool} => (1,),
            )
        )

        # define the value selection
        valueSelection = SeaPearl.LearnedHeuristic(agent)

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x1 = SeaPearl.IntVar(1, 2, "x1", trailer)
        x2 = SeaPearl.IntVar(1, 2, "x2", trailer)
        x3 = SeaPearl.IntVar(2, 3, "x3", trailer)
        x4 = SeaPearl.IntVar(1, 4, "x4", trailer)
        SeaPearl.addVariable!(model, x1)
        SeaPearl.addVariable!(model, x2)
        SeaPearl.addVariable!(model, x3)
        SeaPearl.addVariable!(model, x4)

        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x1, x2, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x2, x3, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x3, x4, trailer))

        # define the variable selection
        variableSelection = SeaPearl.MinDomainVariableSelection()

        # launch the search
        SeaPearl.search!(model, SeaPearl.geometricRBSearch{SeaPearl.InfeasibleNodeCriteria}(3,10,1.1), variableSelection, valueSelection; out_solver=true) #In simple example, we never get an infeasible state

        possible_solutions = [
            Dict("x1" => 1, "x2" => 2, "x3" => 3, "x4" => 1),
            Dict("x1" => 1, "x2" => 2, "x3" => 3, "x4" => 2),
            Dict("x1" => 1, "x2" => 2, "x3" => 3, "x4" => 4),
            Dict("x1" => 2, "x2" => 1, "x3" => 3, "x4" => 1),
            Dict("x1" => 2, "x2" => 1, "x3" => 3, "x4" => 2),
            Dict("x1" => 2, "x2" => 1, "x3" => 3, "x4" => 4),
            Dict("x1" => 2, "x2" => 1, "x3" => 2, "x4" => 1),
            Dict("x1" => 2, "x2" => 1, "x3" => 2, "x4" => 3),
            Dict("x1" => 2, "x2" => 1, "x3" => 2, "x4" => 4)
        ]

        for solution in model.statistics.solutions
            @test solution in possible_solutions
        end

    end

end


text:
# CSucVector
# ----------
#
# A variant of SucVector written in C. This is written to compare the
# performance of the Julia and C implementations.

const libbitvector = "deps/libbitvector"

type CSucVector <: AbstractIndexableBitVector
    ptr::Ptr{Void}

    function CSucVector(src::Union(BitVector,Vector{Bool}))
        ptr = ccall((:make_sucvector, libbitvector), Ptr{Void}, ())
        vec = new(ptr)
        finalizer(vec, vec -> ccall((:delete_sucvector, libbitvector), Void, (Ptr{Void},), vec.ptr))
        chunks = Uint64[]
        i = 1
        while i ≤ endof(src)
            chunk = IndexableBitVectors.read_chunk(src, i)
            push!(chunks, chunk)
            i += 64
        end
        @assert ccall((:read_chunks, libbitvector), Cint, (Ptr{Void}, Ptr{Uint64}, Csize_t), ptr, chunks, length(src)) == 0
        return vec
    end
end

function convert(::Type{CSucVector}, v::Union(BitVector,Vector{Bool}))
    return CSucVector(v)
end

function length(v::CSucVector)
    ccall((:length, libbitvector), Int64, (Ptr{Void},), v.ptr)
end

function getindex(v::CSucVector, i::Integer)
    if !(1 ≤ i ≤ endof(v))
        throw(BoundsError())
    end
    ccall((:access, libbitvector), Bool, (Ptr{Void}, Int64), v.ptr, i - 1)
end

function rank1(v::CSucVector, i::Int)
    if !(0 ≤ i ≤ endof(v))
        throw(BoundsError())
    end
    return unsafe_rank1(v, i)
end

function unsafe_rank1(v::CSucVector, i::Int)
    ccall((:rank1, libbitvector), Uint64, (Ptr{Void}, Int64), v.ptr, i - 1)
end

ext: jl


text:
#data processing
export fitStructures_full, fitStructures_random, getAllEvenRotations,
    sRAAR_bestfit,
    calculateSC,
    postprocess_run,
    get_qrange,
    calculate_cutoff,
    calculate_maximum_resolution,
    load_runs,
    analyse_ensemble,
    load_fitted_intensity,
    load_density,
    average_intensities,
    calculate_optimal,
    load_parameter_list,
    get_parameter_entry,
    get_optimal,
    load_sc_vs_triplets,
    substract_gamma,
    radial_difference_with_intensity,
    denoise_structure

"Fits structure `volume` to structure `reference` by exploring all possible SO(3) rotations"
function fitStructures_full(volume::SphericalHarmonicsVolume, reference::SphericalHarmonicsVolume, num::Int64, K_range::UnitRange{Int64}, L::Int64, ta::Float64 = 0.0, tb::Float64 = float(pi)/2, pa::Float64 =0.0, pb::Float64 = float(pi), ga::Float64 = 0.0, gb::Float64 = float(pi))
    println("Start with fitting")
    #check a single angle combination
    check_angles = function(theta::Float64,phi::Float64, gamma::Float64)
        rvolume = rotateStructure(volume, theta, phi, gamma, maximum(K_range), 2:2:L)
        sc = similarity(reference,rvolume, K_range)
        return (sc, theta, phi, gamma)
    end

    #generate all angle combinations
    angles = [ (theta, phi, gamma) for theta = ta:(tb-ta)/num:tb for phi = pa:(pb-pa)/num:pb for gamma = ga:(gb-ga)/num:gb]

    #check shell correlation for each angle combination
    results = pmap((angles)->check_angles(angles[1], angles[2], angles[3]), angles)

    #Retrieve the best result
    sort!(results, lt=(a,b)->a[1]>b[1])
    bestsc, bestt, bestp, bestg = results[1]
    println("Best intensity correlation: ", bestsc)
    return (rotateStructure(volume, bestt, bestp, bestg, maximum(K_range), 2:2:L), bestsc, bestt, bestp, bestg)
end

function fitStructures_full(volume::SphericalHarmonicsVolume, reference::SphericalHarmonicsVolume, num::Int64, K::Int64, L::Int64, ta::Float64 = 0.0, tb::Float64 = float(pi)/2, pa::Float64 =0.0, pb::Float64 = float(pi), ga::Float64 = 0.0, gb::Float64 = float(pi))
    fitStructures_full(volume, reference, num, 1:K, L, ta, tb, pa, pb, ga, gb)
end

"Fits structure `volume` to structure `reference` by exploring all possible SO(3) rotations"
function fitStructures_random(volume::Volume, reference::Volume, K_range::UnitRange{Int64}, L::Int64, stepsizefactor::Float64, repititions::Int64=4)
    stepsize = pi
    # results = []
    results = @parallel vcat for i = 1:repititions
        phi,theta,gamma = 2*pi*rand(), acos(rand()), 2*pi*rand()
        bestsc = 0.0

        while stepsize > pi/180.0
            new_theta = theta + randn()*stepsize/2
            new_phi = phi + randn()*stepsize
            new_gamma = gamma + randn()*stepsize
            rvolume = rotateStructure(volume, new_theta, new_phi, new_gamma, maximum(K_range), 2:2:L)
            sc = similarity(reference,rvolume, K_range)
            if sc > bestsc
                bestsc,theta,phi,gamma = sc, new_theta, new_phi, new_gamma
            end
            stepsize *= stepsizefactor
        end
        bestsc, theta, phi, gamma
    end
    sort!(results, lt=(a,b)->a[1]>b[1])
    bestsc,theta,phi,gamma = results[1]
    println("Best intensity correlation: ", bestsc)
    return (rotateStructure(volume, theta, phi, gamma, maximum(K_range), 2:2:L), bestsc, theta, phi, gamma)
end

function fitStructures_random(volume::Volume, reference::Volume, K::Int64, L::Int64, stepsizefactor::Float64, repititions::Int64=4)
    fitStructures_random(volume, reference, 1:K, L, stepsizefactor, repititions)
end

# "Calculates the power spectrum of a set of coefficients"
# function powerSpectrum(coeff)
#   sum = 0
#   for k = 1:KMAX
#     for l = 0:LMAX-1
#       # sum =  * sum
#       for m = -l:l
#         sum = sum + 1/(2*l+1)*abs(getc(coeff[k], l, m))^2
#       end
#     end
#   end
#   return sum
# end

"""Takes a structure and returns all 90 deg. rotations"""
function getAllEvenRotations(volume::CubeVolume)
    list = []
    da = pi
    for phi = 0:da:2.0*pi-da
        for theta = 0:da:2.0*pi-da
            let gamma = 0.0
                # for gamma = pi/2:pi/2:2.0*pi
                push!(list, rotateStructure(volume, phi, theta, gamma))
            end
        end
    end
    return list
end

"""Performs a spherical RAAR on the intensity and returns the best fit with the given densityCube (including mirror/rotated images)"""
function sRAAR_bestfit(referenceDensityCube::CubeVolume, intensity::SphericalHarmonicsVolume, iterations::Int64=1001, beta0::Float64 = 0.75, beta_max::Float64 = 0.90, tau::Float64 = 350.0; outputfile="", fourierCube::CubeVolume=CubeVolume(1, 1.0), cutoff_factor::Float64=0.5)

    #Do one phasing
    phased_density = sRAAR(intensity, iterations, beta0, beta_max, tau, cutoff_factor=cutoff_factor)

    #Get result as cube and all mirrored versions
    phased_densityCube = center_cube(getCube(phased_density))
    cubelist = [getAllEvenRotations(phased_densityCube); getAllEvenRotations(mirrorCube(phased_densityCube))]

    #Sort results by best similarity
    cube_cor_list = map((x)->(x,similarity(x,referenceDensityCube)),cubelist)
    sort!(cube_cor_list, lt=(a,b)->a[2]>b[2])

    #Choose best result and return
    best_cube,bestc = cube_cor_list[1]

    if outputfile != "" saveCube(best_cube, outputfile) end
    println("Correlation in realspace of fit: $bestc")#($([cube_cor_list[i][2] for i=1:length(cube_cor_list)]))

    if fourierCube.cubesize > 1
        best_resolution = calculate_maximum_resolution(FSC(best_cube, fourierCube), dr(fourierCube))
        println("Resolution of fit: ", best_resolution)
    end

    return best_cube
end

"""Calculates the FSC from a set of coefficients
Returning: the fitted density cube and the triple of ISC, FSC, ISC_nofitting"""
function calculateSC(volume::SphericalHarmonicsVolume, density::CubeVolume, fourier::CubeVolume, intensity::CubeVolume, num_phases::Int64=1, iterations::Int64=1001, beta0::Float64=0.75, beta1::Float64=0.90, tau::Float64=350.0)

    repetitions = nworkers() > 1 ? (nworkers() < num_phases ? num_phases : nworkers()) : num_phases
    #phasing tries on all available cores
    densities = pmap((vol)->sRAAR_bestfit(density, vol, iterations, beta0, beta1, tau, outputfile="", fourierCube=fourier), collect(Base.Iterators.repeated(volume,repetitions)))

    #Lets calculate the average density
    average_density = reduce(+, densities)

    #Calculate correlations for this result
    isc,fsc,isc_nofitting = ISC(average_density, intensity), FSC(average_density, fourier), shell_correlation_ISC(getCube(volume), intensity)

    println("Best resolution: ", calculate_maximum_resolution(fsc, dr(fourier)))
    println("Best density correlation: ", similarity(average_density,density))

    return average_density,isc,fsc,isc_nofitting
end

"""Performs a postprocessing on a run, including rotational fit in Fourier space and shell correlation calculations"""
function postprocess_run(params, state, reference_pdb_path::String, saving=false, num_points_sphere::Int64=35, sigma::Float64=0.0)

    #Load reference structures as spherical harmonics
    density,fourier,intensity = createSphericalHarmonicsStructure(reference_pdb_path, params["LMAX"], maximum(params["K2_range"]), qmax(maximum(params["K2_range"]), params["qmax"]))
    #Load reference structures as cube
    densCube,fourierCube,intensityCube = createCubicStructure(reference_pdb_path, 2*maximum(params["K2_range"])+1, qmax(maximum(params["K2_range"]), params["qmax"]))

    #intensity from runs
    input_intensity = deleteTerms(state["intensity"],maximum(params["K2_range"]),params["L"])

    #First check if we denoise?
    if sigma > 0.0
        input_intensity = denoise_structure(input_intensity, intensity, maximum(params["K2_range"]), sigma)[1]
    end

    if saving saveCube(input_intensity, "unfitted_intensity.mrc") end

    #Rotational fit in Fourier space
    # state["fittedIntensity"], bestsc, bestt, bestp, bestg = fitStructures(input_intensity, intensity, num_points_sphere, params["K"],params["L"], 0.0, float(pi), 0.0, 2.0*pi, 0.0, 2.0*pi)

    state["fittedIntensity"], bestsc, bestt, bestp, bestg = fitStructures_random(input_intensity, intensity, maximum(params["K2_range"]),params["L"] , 0.99)

    @everywhere gc()

    #Save the fittet intensity
    if saving saveCube(state["fittedIntensity"], "intensity.mrc") end

    #Calculate SC
    fittedDensityCube,isc,fsc,isc_nofitting = calculateSC(state["fittedIntensity"], densCube, fourierCube, intensityCube, 2*nworkers())

    #We are done
    state["sc"] = (isc,fsc,isc_nofitting)
    state["resolution"] = calculate_maximum_resolution(fsc, dr(intensityCube))
    state["state"] = "finished"

    #Save the shell correlations and the phased density
    if saving
        saveState(params, state)
        saveCube(fittedDensityCube,"density.mrc")
    end

    println("Done postprocessing.")
end

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

"""Calculates the corresponding q range for a dq and number of shells K"""
function get_qrange(K::Int64, dq::Float64)
    return collect(1:K) * dq
end

"Given a single fsc curve, calculates k_max for which 0.5=FSC(k_max)"
function calculate_cutoff(sc_curve::Array{Float64}, cutoff::Float64=0.5)
    for k = 1:length(sc_curve)-1
        x1 = k
        x2 = k+1
        y1 = sc_curve[x1]
        y2 = sc_curve[x2]
        if y1 > cutoff && y2 <= cutoff
            m = (y1-y2)/(x1-x2)
            b = y2-m*x2
            kcrossing = (cutoff-b)/m
            return kcrossing
        end
    end
    return length(sc_curve)
end

"Given a single fsc curve, calculates the resolution for which 0.5=FSC(k_max)"
function calculate_maximum_resolution(fsc::Array{Float64}, dq::Float64, cutoff::Float64=0.5)
    return 2.0*pi / ( calculate_cutoff(fsc, cutoff) * dq)
end

"""Alternative approach for resolution calculation"""
function calculate_maximum_resolution(fsc_list::Array{Array{Float64,1},1}, dq::Float64)
    if length(fsc_list) == 0 return 0.0, 0.0 end
    #Calculate the crossing of the mean
    resolutions = Float64[calculate_maximum_resolution(fsc, dq) for fsc in fsc_list]
    return mean(resolutions), std(resolutions)/sqrt(length(fsc_list))
end

"""Given a run in `dir`, returns the isc, fsc, isc_nofitting"""
function load_runs(dir::String, densityCube::CubeVolume, fourierCube::CubeVolume, intensityCube::CubeVolume; range=1000:1019)
    list = Any[]
    energylist = Float64[]
    density_list = CubeVolume[]
    intensity_list = CubeVolume[]
    for i in range
        try
            state = deserializeFromFile("$dir/$i/state.dat")
            push!(list, state["sc"])
            push!(energylist, state["E"])

            density = loadCube("$dir/$i/density.mrc")
            push!(density_list, density)

            intensity = loadCube("$dir/$i/intensity.mrc")
            push!(intensity_list, intensity)
        catch
            println("Couldn't load $dir/$i")
        end
    end

    return analyse_ensemble(density_list, intensity_list, densityCube, fourierCube, intensityCube)
end

# """Takes an ensemble of spherical correlations and calculates max, min, mean and std deviation"""
# function analyse_ensemble(sc_list::Array{Array{Float64}}, K::Int64)
#     if length(sc_list) > 0
#       maximum_sc = Float64[ maximum([length(sc_list[i]) > 1 ? sc_list[i][k] : 0.0 for i = 1:length(sc_list)]) for k=1:K]
#       minimum_sc = Float64[ minimum([length(sc_list[i]) > 1 ? sc_list[i][k] : 0.0 for i = 1:length(sc_list)]) for k=1:K]
#       mean_sc = Float64[ mean([length(sc_list[i]) > 1 ? sc_list[i][k] : 0.0 for i = 1:length(sc_list)]) for k=1:K]
#       stderr_sc = Float64[ std(Float64[length(sc_list[i]) > 1 ? sc_list[i][k] : 0.0 for i = 1:length(sc_list)]) for k=1:K]/sqrt(length(sc_list))
#       return Dict("max"=>maximum_sc,"min"=>minimum_sc,"mean"=>mean_sc,"std"=>stderr_sc)
#     else
#       return Dict("max"=>Float64[], "min"=>Float64[], "mean"=>Float64[], "std"=>Float64[])
#     end
#
# end

"""Takes an ensemble of densities/intensities and calculates the correlations"""
function analyse_ensemble(density_list::Array{CubeVolume}, intensity_list::Array{CubeVolume}, densityCube::CubeVolume, fourierCube::CubeVolume, intensityCube::CubeVolume)

    average_density = reduce(+, density_list)
    average_intensity = reduce(+, intensity_list)

    #Calculate averaged lines
    isc,fsc,isc_nofitting = ISC(average_density, intensityCube), FSC(average_density, fourierCube), shell_correlation_ISC(average_intensity, intensityCube)

    isc_list = map((density)-> ISC(density, intensityCube), density_list)
    fsc_list = map((density)-> FSC(density, fourierCube), density_list)
    isc_nofitting_list = map((intensity)->shell_correlation_ISC(intensity, intensityCube), intensity_list)

    res = calculate_maximum_resolution(fsc, dr(intensityCube))
    _,res_err = calculate_maximum_resolution(fsc_list, dr(intensityCube))

    get_err = (list) -> Float64[ std(Float64[list[i][k] for i = 1:length(list)]) for k=1:length(list[1])]/sqrt(length(list))

    isc_err = get_err(isc_list)
    fsc_err = get_err(fsc_list)
    isc_nofitting_err = get_err(isc_nofitting_list)
    return Dict("isc"=>isc, "isc_err"=>isc_err, "fsc"=>fsc, "fsc_err"=>fsc_err, "isc_nofitting"=>isc_nofitting, "isc_nofitting_err"=>isc_nofitting_err, "res"=>res, "res_err"=>res_err)

end

function load_fitted_intensity(path::String, energy_limit::Float64=1.0e10)
    statepath = "$path/state.dat"
    if isfile(statepath)
        res = deserializeFromFile(statepath)
        if res["E"] < energy_limit && haskey(res, "fittedIntensity")
            return res["fittedIntensity"]
        end
    end
end

function load_density(path::String)
    densitypath = "$path/density.mrc"
    if isfile(densitypath)
        return loadCube(densitypath)
    end
end

#e.g ../parallel/paper_res_vs_L_rotate_all_at_once_Bayes_N32_K16_I0.05_TD0.9999_SF1.01/P163840000_L14
function average_intensities(root::String, output_path::String, densityCube::CubeVolume, fourierCube::CubeVolume, intensityCube::CubeVolume; iterations::Int64=1001, beta0::Float64=0.75, beta1::Float64=0.90, tau::Float64=350.0, energy_limit::Float64=1.0e10, do_phasing::Bool=true)
    try mkpath(output_path) end
    #Average intensity
    intensity_list = [load_fitted_intensity("$root/$num", energy_limit) for num = 1000:1019]
    intensity_list = intensity_list[intensity_list.!=nothing]

    density_list = [load_density("$root/$num") for num = 1000:1019]
    density_list = density_list[density_list.!=nothing]

    average_intensity = getSphericalHarmonicsVolume(reduce(+, map(getSurfaceVolume, intensity_list)))
    average_density = reduce(+, density_list)
    average_density_res = calculate_maximum_resolution(FSC(average_density, fourierCube), dr(fourierCube))

    println("Averaged density has resoution of: $(average_density_res)")
    saveCube(average_intensity, "$(output_path)/intensity.mrc")
    saveCube(average_density, "$(output_path)/averaged_density.mrc")

    average_res = ()
    #Phase averaged intensity
    if do_phasing
        average_res = calculateSC(average_intensity, densityCube, fourierCube, intensityCube, 6, iterations, beta0, beta1, tau)
        saveCube(average_res[1], "$(output_path)/density.mrc")
    end
    serializeToFile("$(output_path)/state.dat", Dict("sc"=>average_res, "fittedIntensity"=>average_intensity))
end

# """Calculates the phasing of the optimal spherical harmonics structure with the given cutoff `K` and `L`"""
function calculate_optimal(intensity::SphericalHarmonicsVolume, densityCube::CubeVolume, fourierCube::CubeVolume, intensityCube::CubeVolume, K::Int64, L::Int64, dir::String, num_phases::Int64=8)
    res = calculateSC(deleteTerms(intensity, K, L), densityCube, fourierCube, intensityCube, num_phases, 1001, 0.75, 0.90, 350.0)
    saveCube(res[1], "$dir/density.mrc")
    saveCube(deleteTerms(intensity, K, L), "$dir/intensity.mrc")
    serializeToFile("$dir/optimal.dat", res)
end

function load_parameter_list(reference_intensity::SphericalHarmonicsVolume)
    parameter_list = Dict[]
    #K2 must be at least such that rmax > 35.0 Angstrom (Phasing)
    deltar = 1.0 #2 Angstrom resolution maximum
    qmax2 = pi/deltar
    K2_min = ceil(Int64, 35.0 / deltar)
    for K2 = K2_min:3:K2_min+12
        deltaq = qmax2 / K2
        for K3::Int64 = 20:2:min(K2,40)
            for L3::Int64 = 12:2:floor(Int64,(K2-1)/2)
                dir = "../parallel/optimal/optimal_K2_$(K2)_K3_$(K3)_L3_$(L3)"
                try
                    optimal = deserializeFromFile("$dir/optimal.dat")
                    retrieved = deserializeFromFile("$dir/retrieved.dat")
                    density = loadCube("$dir/density.mrc")
                    intensity = getCube(deleteTerms(reference_intensity,K3,L3))#loadCube("$dir/intensity.mrc")
                    res_optimal = calculate_maximum_resolution(optimal[3], deltaq)
                    res_retrieved = calculate_maximum_resolution(retrieved[3], deltaq)
                    println("K2=$(K2), K3=$(K3), L3=$(L3), res=$(res_optimal), retrieved=$(res_retrieved)")
                    #Note: K3*(K3+1)*(K3+2)/6*L3^4 is a rough estimation for the total number of operations needed to calculate the model, num_parameters is falsely named at this point
                    push!(parameter_list, Dict("density"=>density, "intensity"=>intensity, "K2"=>K2, "K3"=>K3, "L3"=>L3, "res_optimal"=>res_optimal, "res_retrieved"=>res_retrieved, "num_parameters"=>K3*(K3+1)*(K3+2)/6*L3^4, "optimal"=>optimal, "res_max"=>2*pi/(K3*deltaq)))
                catch
                    println("Failed loading $dir")
                end
            end
        end
    end
    return parameter_list
end

function get_parameter_entry(parameter_list::Array{Dict}, K::Int64, L::Int64)
    filter((x)-> x["K3"] == K && x["L3"] == L,parameter_list)[1]
end

function get_optimal(parameter_list::Array{Dict}, K::Int64, L::Int64)
    return get_parameter_entry(parameter_list, K,L)
end

function load_sc_vs_triplets(dict::Dict, densityCube::CubeVolume, fourierCube::CubeVolume, intensityCube::CubeVolume, K::Int64, L::Int64, dir::String, parameter_list::Array, correlation_list::Dict)
    if !haskey(dict,K) dict[K] = Dict() end
    dict[K][L] = Dict( name => merge(info, load_runs("$dir/$(name)", densityCube, fourierCube, intensityCube; range=1000:1019)) for (name, info) in correlation_list)

    opt = get_optimal(parameter_list, K, L)
    dict[K][L]["optimal"] = merge(
    Dict("shortlabel"=>"optimal", "label"=>"optimal", "linestyle"=>"--", "color"=>"grey"),
    analyse_ensemble([opt["density"]], [opt["intensity"]], densityCube, fourierCube, intensityCube)
    )
end

##################################################################
## noise
#################################################################

function substract_gamma(volume::SphericalHarmonicsVolume, intensity::SphericalHarmonicsVolume, factor::Float64, sigma::Float64, K::Int64)
    new_volume = deepcopy(volume)
    factor *= 1/sigma^3/(2.0*pi)^(3/2)
    for k = 1:K
        setc(new_volume, k, 0, 0, getc(volume,k,0,0) - factor*gaussian_distribution(k*dr(intensity), 0.0, sigma))
    end
    return new_volume
end

function radial_difference_with_intensity(volume::SphericalHarmonicsVolume, intensity::SphericalHarmonicsVolume, K::Int64)
    a = [real(getc(volume,k,0,0)) for k=1:K]
    a /= norm(a)
    b = [real(getc(intensity,k,0,0)) for k=1:K]
    b /= norm(b)
    return norm(a-b)
end

"""Denoises a single structure by making radial part of l=0,m=0 components fit the known radial distribution"""
function denoise_structure(volume::SphericalHarmonicsVolume, intensity::SphericalHarmonicsVolume, K::Int64, sigma::Float64)
    mydiff(x) = radial_difference_with_intensity(substract_gamma(volume, intensity, x, sigma, K), intensity, K)
    result = optimize(mydiff, 0.0, 1.0)
    gammamin = Optim.minimizer(result)
    return substract_gamma(volume, intensity, gammamin, sigma, K), gammamin
end

##################################################################
## Exp. data
#################################################################

export complete_core, phase_completed_intensity

"""Averages intensities and completes core via fitting"""
function complete_core(name::String, c1::C1, center_range::UnitRange{Int64}, range::UnitRange{Int64}; plotting::Bool=true)
    intensities = [deserializeFromFile("$name/$i/intensity.dat") for i = 1000:1019]
    average_intensity_surf = reduce(+, map(getSurfaceVolume, intensities))

    #Filter negativity
    average_intensity_surf.surf = map((x)-> max.(real(x), 0.0), average_intensity_surf.surf)
    average_intensity = getSphericalHarmonicsVolume(average_intensity_surf)
    saveCube(average_intensity, "intensity_averaged.mrc")

    reference_intensity = deepcopy(average_intensity)
    reference_surf = getSurfaceVolume(reference_intensity)
    c1_coliphage = [sum(abs, reference_surf.surf[k]) for k in 1:maximum(range)]

    # curve = deepcopy(c1_coliphage)
    fitting_range=4:10
    curve = log(deepcopy(c1))
    curve[1:3] = zeros(3)
    a = poly_fit(collect(fitting_range), curve[fitting_range], 2)
    func(x,a) = a[1] + x*a[2]+x^2*a[3]
    center_fit = [func(k,a) for k in 1:10]

    corrected_intensity = deepcopy(average_intensity)
    corrected_surf = getSurfaceVolume(corrected_intensity)
    shift = log.(sum(abs, corrected_surf.surf[minimum(range)]))-center_fit[minimum(range)]
    println("shfit = $shift")

    for k in center_range
        corrected_surf.surf[k] = exp(center_fit[k]+shift)*ones(length(corrected_surf.surf[k])) / length(corrected_surf.surf[k])
        # corrected_surf.surf[k] = shift*c1[k]*ones(length(corrected_surf.surf[k])) / length(corrected_surf.surf[k])
    end

    corrected_intensity = getSphericalHarmonicsVolume(corrected_surf)
    c1_coliphage_corrected = [sum(abs, corrected_surf.surf[k]) for k in 1:maximum(range)]
    saveCube(getSphericalHarmonicsVolume(corrected_surf), "intensity_averaged_corrected.mrc")

    if plotting == true
        plot(collect(1:10), curve[1:10], lw=5, label="orig")
        plot(collect(1:10), center_fit, label="fit 1")
        legend()
        savefig("fig1.pdf")

        figure()
        title("Radial sum of coliphage")
        plot(collect(1:maximum(range)),log(c1_coliphage), label="Surf coliphage", lw=3)
        # plot(collect(1:10), center_fit + shift, label="center fit")
        plot(collect(1:maximum(range)), log(c1_coliphage_corrected), label="surf coliphage corr.")
        ylabel("Sum of Radial Part")
        xlim(1,maximum(range))
        # ylim(1.0e-2, 5.0e3)
        xlabel("Shell")
        legend()
        savefig("fig2.pdf")
    end

    # Extend and add zeros at the end
    extended_corrected_intensity = deepcopy(corrected_intensity)
    extend_range = maximum(range)+1:40
    for i in extend_range
        push!(extended_corrected_intensity.coeff, zeros(Complex{Float64}, Base.size(extended_corrected_intensity.coeff[1])))
    end
    extended_corrected_intensity.KMAX += length(extend_range)
    return extended_corrected_intensity
end

function phase_completed_intensity(extended_corrected_intensity::SphericalHarmonicsVolume, num_tries::Int64=8, beta_end::Float64=0.90)

    reference_density_SH = sRAAR(extended_corrected_intensity, 1001, 0.75, beta_end, 350.0, cutoff_factor=0.3)
    reference_density = center_cube(getCube(reference_density_SH))

#     densities = pmap((vol)->sRAAR_bestfit(getCube(extended_corrected_intensity), vol,  1001, 0.75, 0.95, 350.0, cutoff_factor=0.3), collect(Base.Iterators.repeated(extended_corrected_intensity,num_tries-1)))
    densities = pmap((vol)->sRAAR_bestfit(reference_density, vol,  1001, 0.75, beta_end, 350.0, cutoff_factor=0.3), collect(Base.Iterators.repeated(extended_corrected_intensity,num_tries-1)))
    push!(densities, reference_density)
    map!((vol)->center_cube(vol), densities)
    for i = 1:length(densities) saveCube(densities[i], "density_$(i).mrc") end

    averaged_density = reduce(+, densities)

    serializeToFile("density_averaged.dat", averaged_density)

    saveCube(averaged_density, "density_averaged.mrc")
    return averaged_density
end


text:
module GeometricLagrangianRefFEsTests

using Test
using Gridap.Io
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Polynomials
using Gridap.ReferenceFEs
using JSON

orders = (2,3)
b = MonomialBasis(Float64,QUAD,orders)
r = [(0,0), (1,0), (2,0), (0,1), (1,1), (2,1), (0,2), (1,2), (2,2), (0,3), (1,3), (2,3)]
@test get_exponents(b) == r

orders = (1,1,2)
b = MonomialBasis(Float64,WEDGE,orders)
r = [(0,0,0), (1,0,0), (0,1,0), (0,0,1), (1,0,1), (0,1,1), (0,0,2), (1,0,2), (0,1,2)]
@test get_exponents(b) == r

orders = (1,1,1)
b = MonomialBasis(Float64,PYRAMID,orders)
r = [(0,0,0), (1,0,0), (0,1,0), (1,1, 0), (0,0,1)]
@test get_exponents(b) == r

orders = (1,1,1)
b = MonomialBasis(Float64,TET,orders)
r = [(0,0,0), (1,0,0), (0,1,0), (0,0,1)]
@test get_exponents(b) == r

orders = (2,2)
extrusion = Tuple(QUAD.extrusion)

dofs = LagrangianDofBasis(VectorValue{3,Float64},TET,1)
@test dofs.nodes == Point{3,Float64}[(0,0,0), (1,0,0), (0,1,0), (0,0,1)]
@test dofs.node_and_comp_to_dof == VectorValue{3,Int}[(1,5,9), (2,6,10), (3,7,11), (4,8,12)]

dofs = LagrangianDofBasis(Float64,WEDGE,(2,2,2))
r = Point{3,Float64}[
  (0.0,0.0,0.0),(1.0,0.0,0.0),(0.0,1.0,0.0),
  (0.0,0.0,1.0),(1.0,0.0,1.0),(0.0,1.0,1.0),
  (0.5,0.0,0.0),(0.5,0.0,1.0),(0.0,0.5,0.0),
  (0.5,0.5,0.0),(0.0,0.5,1.0),(0.5,0.5,1.0),
  (0.0,0.0,0.5),(1.0,0.0,0.5),(0.0,1.0,0.5),
  (0.5,0.0,0.5),(0.0,0.5,0.5),(0.5,0.5,0.5)]
@test dofs.nodes == r

dofs = LagrangianDofBasis(VectorValue{2,Int},VERTEX,())
@test dofs.node_and_comp_to_dof == VectorValue{2,Int}[(1,2)]

b = MonomialBasis(VectorValue{2,Int},VERTEX,())
@test length(b) == 2
@test evaluate(b,Point{0,Int}[(),()]) == VectorValue{2,Int}[(1, 0) (0, 1); (1, 0) (0, 1)]

reffe = LagrangianRefFE(VectorValue{2,Int},VERTEX,())
@test get_face_own_nodes(reffe) == [[1]]
@test get_face_own_dofs(reffe) == [[1,2]]
@test get_face_own_dofs_permutations(reffe) == [[[1, 2]]]
test_lagrangian_reference_fe(reffe)
@test ReferenceFE{0}(reffe,1) === reffe

reffe = LagrangianRefFE(VectorValue{2,Float64},SEGMENT,(2,))
@test get_face_own_dofs(reffe) == [[1, 4], [2, 5], [3, 6]]
test_lagrangian_reference_fe(reffe)

reffe = LagrangianRefFE(VectorValue{2,Float64},TRI,3)
@test get_face_own_dofs(reffe) == [[1, 11], [2, 12], [3, 13], [4, 5, 14, 15], [6, 7, 16, 17], [8, 9, 18, 19], [10, 20]]
test_lagrangian_reference_fe(reffe)

reffe = ReferenceFE(TRI,lagrangian,VectorValue{2,Float64},3)
@test get_face_own_dofs(reffe) == [[1, 11], [2, 12], [3, 13], [4, 5, 14, 15], [6, 7, 16, 17], [8, 9, 18, 19], [10, 20]]
test_lagrangian_reference_fe(reffe)

@test Lagrangian() == lagrangian

reffe = LagrangianRefFE(Float64,HEX,2)
test_lagrangian_reference_fe(reffe)

reffe = LagrangianRefFE(Float64,WEDGE,(1,1,2))
test_lagrangian_reference_fe(reffe)
refface = ReferenceFE{1}(reffe,3)
@test get_face_own_dofs(refface) == [[1], [2], []]
refface = ReferenceFE{1}(reffe,8)
@test get_face_own_dofs(refface) == [[1], [2], [3]]

orders = (4,)
reffe = LagrangianRefFE(VectorValue{2,Float64},SEGMENT,orders)
@test get_own_nodes_permutations(reffe) == [[1, 2, 3], [3, 2, 1]]
@test get_own_dofs_permutations(reffe) == [[1, 2, 3, 4, 5, 6], [3, 2, 1, 6, 5, 4]]

orders = (2,3)
reffe = LagrangianRefFE(VectorValue{2,Float64},QUAD,orders)
@test get_own_nodes_permutations(reffe) ==[[1, 2], [0, 0], [1, 2], [0, 0], [0, 0], [2, 1], [0, 0], [2, 1]]
@test get_own_dofs_permutations(reffe) == [
  [1, 2, 3, 4], [0, 0, 0, 0], [1, 2, 3, 4], [0, 0, 0, 0],
  [0, 0, 0, 0], [2, 1, 4, 3], [0, 0, 0, 0], [2, 1, 4, 3]]

reffe = LagrangianRefFE(VectorValue{2,Float64},QUAD,2)
@test get_node_and_comp_to_dof(reffe) == VectorValue{2,Int}[
  (1, 10), (2, 11), (3, 12), (4, 13), (5, 14), (6, 15), (7, 16), (8, 17), (9, 18)]

reffe = LagrangianRefFE(VectorValue{2,Float64},QUAD,2)
d = 1
@test get_face_own_dofs(reffe,d) == [[5, 14], [6, 15], [7, 16], [8, 17]]
@test get_face_own_nodes(reffe,d) == [[5], [6], [7], [8]]

# 0-order degenerated case

orders = (0,0)

b = compute_monomial_basis(Float64,QUAD,orders)
@test length(b.terms) == 1

own_nodes = compute_own_nodes(QUAD,orders)
@test own_nodes == Point{2,Float64}[(0.5, 0.5)]

face_orders = compute_face_orders(QUAD,SEGMENT,1,orders)
@test face_orders == (0,)

nodes, facenodeids = compute_nodes(QUAD,orders)
@test nodes == Point{2,Float64}[(0.5, 0.5)]
@test facenodeids == [Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], [1]]

nodeperms = compute_own_nodes_permutations(QUAD,own_nodes)
@test nodeperms == [[1], [1], [1], [1], [1], [1], [1], [1]]

reffaces = compute_lagrangian_reffaces(Float64,QUAD,orders)

reffe = LagrangianRefFE(Float64,QUAD,orders)
@test num_dofs(reffe) == 1
@test num_nodes(reffe) == 1
test_lagrangian_reference_fe(reffe)
@test Conformity(reffe) == L2Conformity()

@test get_face_own_nodes(reffe) == Vector{Int}[[], [], [], [], [], [], [], [], [1]]
#@test get_face_own_nodes_permutations(reffe) == Vector{Vector{Int}}[
#  [[]],[[]],[[]],[[]],[[],[]],[[],[]],[[],[]],[[],[]],[[1],[1],[1],[1],[1],[1],[1],[1]] ]
#@test get_own_nodes_permutations(reffe) == [[1], [1], [1], [1], [1], [1], [1], [1]]
@test get_face_own_nodes_permutations(reffe) == Vector{Vector{Int}}[
  [[]],[[]],[[]],[[]],[[]],[[]],[[]],[[]],[[1]] ]
@test get_own_nodes_permutations(reffe) == [[1]]

reffe = LagrangianRefFE(VectorValue{2,Float64},QUAD,orders)
@test num_dofs(reffe) == 2
@test num_nodes(reffe) == 1
test_lagrangian_reference_fe(reffe)

@test get_face_own_dofs(reffe) == Array{Int,1}[[], [], [], [], [], [], [], [], [1, 2]]
#@test get_face_own_dofs_permutations(reffe) == [
#  [[]],[[]],[[]],[[]],[[],[]],[[],[]],[[],[]],[[],[]],[[1,2],[1,2],[1,2],[1,2],[1,2],[1,2],[1,2],[1,2]]]
@test get_face_own_dofs_permutations(reffe) == [
  [[]],[[]],[[]],[[]],[[]],[[]],[[]],[[]],[[1,2]]]

own_nodes = compute_own_nodes(TRI,orders)
@test own_nodes == Point{2,Float64}[(1.0/3,1.0/3)]

# More API

reffe = LagrangianRefFE(Float64,WEDGE,1)

reffes = get_reffaces(ReferenceFE{2}, reffe)
iface_to_ftype = get_face_type(reffe,2)

@test length(reffes) == 2
@test iface_to_ftype == [1, 1, 2, 2, 2]

order = 4
reffe = LagrangianRefFE(Float64,TET,order)
@test get_order(reffe) == order
@test get_orders(reffe) == (order,order,order)
@test is_P(reffe) == true
@test is_Q(reffe) == false
@test is_S(reffe) == false

orders = (1,2)
reffe = LagrangianRefFE(Float64,QUAD,orders)
@test get_order(reffe) == 2
@test get_orders(reffe) == orders
@test is_P(reffe) == false
@test is_Q(reffe) == true
@test is_S(reffe) == false

reffe = LagrangianRefFE(Float64,QUAD,(2,3))
@test reffe == from_dict(LagrangianRefFE,to_dict(reffe))

reffe = LagrangianRefFE(Float64,QUAD,(2,3))

s = to_json(reffe)
@test reffe == from_json(LagrangianRefFE,s)

s = JSON.json(reffe)
@test reffe == from_json(LagrangianRefFE,s)

d = mktempdir()
f = joinpath(d,"reffe.jld2")

to_jld2_file(reffe,f)
@test reffe == from_jld2_file(typeof(reffe),f)

rm(d,recursive=true)

end # module


text:
"""
NETALIGNMR
----------
    solve the network alignment problem with Klau's algorithm.
    For details: [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-10-S1-S59]
    Usage
    ----
    Input:
    ----
    - `S`: the complete set of squares
    - `w`: the matching weights for all edges in the link graph L
    - `a`: the value of alpha in the netalign objective
    - `b`: the value of beta in the netalign objective
    - `li`: the start point of each edge in L
    - `lj`: the end point of each edge in L
    - `gamma`: starting step value (default 0.25)
    - `stepm`: number of non-decreasing iterations adjusting the step length (default 25)
    - `rtype`: the rounding type (default 1)
        * rtype = 1: only consider the current matching
        * rtype = 2: enrich the matching with info from other squares
    - `maxiter`: maximum number of iterations to take (default 1000)
    - `verbose`: output at each iterations (default false)
    Methods:
    -------
    xbest,st,status,hist = netalignmr(S,w,a,b,li,lj,gamma,stepm,rtype,maxiter)
    xbest,st,status,hist = netalignmr(S,w,a,b,li,lj,gamma,stepm,rtype)
    xbest,st,status,hist = netalignmr(S,w,a,b,li,lj,gamma,stepm)
    xbest,st,status,hist = netalignmr(S,w,a,b,li,lj,gamma)
    xbest,st,status,hist = netalignmr(S,w,a,b,li,lj)
    Example:
    -------
    S,w,li,lj = netalign_setup(A,B,L)
    xbest,st,status,hist = netalignmr(S,w,a,b,li,lj)
    ma,mb = edge_list(bipartite_matching(xbest,li,lj))
"""
function netalignmr(S::SparseMatrixCSC{T1,Int},w::Vector{T2},
                    a::Int,b::Int,li::Vector{Int},lj::Vector{Int},
                    gamma::Float64=0.25,stepm::Int=25,rtype::Int=1,maxiter::Int=1000,verbose::Bool=false) where {T1,T2}

  m = maximum(li)
  n = maximum(lj)
  Matching_Setup = bipartite_matching_setup(w,li,lj,m,n)
  tripi = Matching_Setup.tripi
  matm = Matching_Setup.m
  matn = Matching_Setup.n
  rp = Matching_Setup.rp
  ci = Matching_Setup.ci
  mperm = tripi[tripi.>0]

  U = spzeros(Float64,size(S,1),size(S,2))
  xbest = zeros(Float64,length(w))

  flower = 0.0                      #best lower bound on the solution
  fupper = Inf                      #best upper bound on the solution
  next_reduction_iteration = stepm  #reduction step

  hist = zeros(Float64,maxiter,7)
  iter = 1
  @inbounds for iter = 1:maxiter
    Q = (b/2)*S + U-copy(U')
    Qt = copy(Q')
    Qp = Qt.colptr
    Qr = Qt.rowval
    Qv = Qt.nzval
    M = size(Qt,1)
    N = size(Qt,2)
    nedges = length(li)

    all_matching = column_maxmatchsum(M,N,Qp,Qr,Qv,m,n,nedges,li,lj)
    q = all_matching.q
    qj = all_matching.mi
    qi = all_matching.mj
    medges = all_matching.medges

    SM = sparse(qi[1:medges],qj[1:medges],1,size(Q,1),size(Q,2))

    x = a*w + q
    ai = zeros(Float64,length(tripi))

    ai[tripi.>0] = x[mperm]

    M_output = bipartite_matching_primal_dual(rp,ci,ai,matm,matn)
    mi = MatrixNetworks.edge_indicator(M_output,li,lj)
    val = M_output.weight

    # compute stats
    matchval = dot(mi,w)
    overlap = dot(mi,(S*mi)/2)
    card = M_output.cardinality
    f = a*matchval + b*overlap

    if val < fupper
        fupper = val
        next_reduction_iteration = iter+stepm
    end

    if f > flower
        flower = f
        itermark = "*"
        xbest = convert(Vector{Float64},mi)
    else
        itermark = " "
    end

    if rtype == 1
        # no work
    elseif rtype==2
        mw = S*x
        mw = a*w + b/2*mw

        ai = zeros(Float64,length(tripi))
        ai[tripi.>0] = mw[mperm]

        M_output2 = MatrixNetworks.bipartite_matching_primal_dual(rp,ci,ai,matm,matn)
        mx = MatrixNetworks.edge_indicator(M_output2,li,lj)
        card = M_output2.cardinality
        matchval = dot(w,mx)
        overlap = dot(mx,(S*mx)/2)

        f = a*matchval + b*overlap

        if f > flower
            flower = f
            itermark = "**"
            mi = mx
            xbest = mw
        end
    end

    # report on current iter
    hist[iter,1:end] = [norm(nonzeros(U),1), flower, fupper, f, matchval, card, overlap]

    # the below if statement causes type instability due to @printf being type instable
    if verbose
        @printf("%5s   %4i   %8.1e   %7.2f %7.2f %7.2f  %7.2f %7.2f %7i %7i\n",
            itermark, iter, norm(nonzeros(U),1),
            flower, fupper, val,
            f, matchval, card, overlap)
    end
    if iter == next_reduction_iteration
        gamma = gamma*0.5
        # the below if statement causes type instability due to @printf being type instable
        if verbose
            @printf("%5s   %4s   reducing step to %f\n", "", "", gamma);
        end
        if gamma < 1e-24
            hist = hist[1:iter,:]
            break
        end
        next_reduction_iteration = iter+stepm
    end
    if (fupper-flower) < 1e-2
        hist = hist[1:iter,:]
        break
    end
    Wtemp = sparse(1:length(mi),1:length(mi),gamma*mi)
    U = U - Wtemp*triu(SM) + tril(SM)'*Wtemp
    Utemp1 = fillstored!(copy(U), 1)
    Utemp1 *= 0.5
    U = min.(U,Utemp1)
    Utemp1 *= -1
    U = max.(U,Utemp1)
  end
  status = zeros(Float64,2)

  st = (fupper-flower) < 1e-2
  status[1] = flower
  status[2] = fupper
  return xbest,st,status,hist

end

function netalignmr(S::SparseMatrixCSC{T,Int},L::SparseMatrixCSC{U,Int},
                    a::Int,b::Int,gamma::Float64=0.25,stepm::Int=25,
                    rtype::Int=1,maxiter::Int=1000,verbose::Bool=false) where {T,U}

     li,lj,w = findnz(L)
     netalignmr(S,w,a,b,li,lj,gamma,stepm,rtype,maxiter,verbose)
end


text:
using CrystalNets
using Test, Random
using PeriodicGraphs
using StaticArrays
using Graphs
using Combinatorics
import Base.Threads

CrystalNets.toggle_export(false)

function _finddirs()
    root = dirname(dirname(pathof(CrystalNets)))
    return joinpath(root, "test", "cif"), root
end

function capture_out(name)
    result = open(name, "w") do out
        redirect_stderr(devnull) do
            redirect_stdout(CrystalNets.julia_main, out)
        end
    end
    written = readlines(name)
    return result, written
end

const safeARCHIVE = deepcopy(CrystalNets.CRYSTAL_NETS_ARCHIVE)
const safeREVERSE = deepcopy(CrystalNets.REVERSE_CRYSTAL_NETS_ARCHIVE)
function __reset_archive!(safeARCHIVE, safeREVERSE)
    empty!(CrystalNets.CRYSTAL_NETS_ARCHIVE)
    empty!(CrystalNets.REVERSE_CRYSTAL_NETS_ARCHIVE)
    merge!(CrystalNets.CRYSTAL_NETS_ARCHIVE, safeARCHIVE)
    merge!(CrystalNets.REVERSE_CRYSTAL_NETS_ARCHIVE, safeREVERSE)
    nothing
end

import CrystalNets.Clustering: SingleNodes, AllNodes, Standard, PE, PEM

@testset "MOF examples" begin
    cifs, crystalnetsdir = _finddirs()
    kwargs = (; structure=StructureType.MOF, clusterings=[Clustering.Auto,Standard,PE,PEM])
    println(stderr, "The following warning about warnings altering performance is expected")
    mofdataset = determine_topology_dataset(joinpath(cifs, "MOFs"), false; kwargs...)

    @testset "Dataset analysis" begin
        @test length(mofdataset) == 14

        hkust1 = mofdataset["HKUST-1.cif"]
        @test hkust1[SingleNodes] == hkust1[AllNodes] == hkust1[Standard]
        @test hkust1[SingleNodes].name == "tbo"

        jxust1 = mofdataset["JXUST-1.cif"]
        @test jxust1[SingleNodes] == jxust1[AllNodes]
        @test PeriodicGraph(jxust1[SingleNodes]) == PeriodicGraph(REVERSE_CRYSTAL_NETS_ARCHIVE["pcu"])

        mil53 = mofdataset["MIL-53.cif"]
        @test string(mil53[SingleNodes]) == mil53[Standard].name == "bpq"
        @test mil53[AllNodes].name == "rna"
        @test mil53[PE].name == "sra, ABW"

        mil100 = mofdataset["MIL-100.cif"]
        @test mil100[SingleNodes].name == mil100[AllNodes].name == "moo"

        mil101 = mofdataset["MIL-101.cif"]
        @test mil101[SingleNodes].name == mil101[AllNodes].name == "mtn-e"
        @test mil101[PE].name == "mtn-e-a"

        mof5 = mofdataset["MOF-5.cif"]
        @test mof5[SingleNodes].name == mof5[AllNodes].name == mof5[Standard].name == "tbo"

        mof14 = mofdataset["MOF-14.cif/1"]
        @test mof14 == mofdataset["MOF-14.cif/2"]
        @test mof14[SingleNodes].name == mof14[AllNodes].name == mof14[Standard].name == "pto"
        @test mof14[PE] == parse(TopologicalGenome, "sqc11259")

        mof801 = mofdataset["MOF-801.cif"]
        @test mof801[SingleNodes].name == mof801[AllNodes].name == "fcu"
        @test mof801[Standard].name == "xbi"
        @test mof801[PE].name == "ubt"

        pcn700 = mofdataset["PCN-700.cif"]
        @test pcn700[SingleNodes].name == pcn700[AllNodes].name == "bcu"
        @test pcn700[PE].name == "pcb, ACO"

        uio66 = mofdataset["UiO-66.cif"]
        @test uio66[SingleNodes].name == uio66[AllNodes].name == "fcu"
        @test uio66[Standard].name == "xbi"
        @test uio66[PE].name == "ubt"

        zif8 = mofdataset["ZIF-8.cif"]
        @test mofdataset["ZIF-8.cif"] == parse(TopologyResult, "AllNodes, SingleNodes, Standard, PEM: sod, SOD | PE: sod-e")

        zif67 = mofdataset["ZIF-67.cif"]
        @test string(zif67) == "AllNodes, SingleNodes, Standard, PEM: sod, SOD\nPE: sod-e"
    end

    CrystalNets.toggle_warning(false)
    @test mofdataset["MIL-53.cif"] == determine_topology(joinpath(cifs, "MIL-53.cif"); kwargs...)
    @test_throws ArgumentError determine_topology(joinpath(cifs, "MIL-53.cif"); kwargs..., bonding=Bonding.Input)

    @test mofdataset["UiO-66.cif"] == determine_topology(joinpath(cifs, "UiO-66.cif"); kwargs...)

    juc101 = determine_topology(joinpath(cifs, "JUC-101.cif"); kwargs..., bonding=Bonding.Guess)
    @test juc101 == determine_topology(joinpath(cifs, "JUC-101.cif"); kwargs..., bonding=Bonding.Input)
    @test juc101[SingleNodes].name == "nia"
    @test juc101[AllNodes].name == "jjt"
    @test string(juc101[Standard].genome) == "3 1 2 0 0 0 1 2 0 1 1 1 3 0 0 0 1 3 0 1 0 1 4 0 0 0 1 4 0 0 1 1 5 0 0 0 1 5 0 1 1 1 6 0 0 0 1 6 0 1 0 1 7 0 0 0 1 7 0 0 1 2 3 0 0 -1 2 4 0 0 0 2 8 0 0 0 2 8 0 1 1 3 4 0 0 1 3 8 0 0 1 3 8 0 1 1 4 8 0 1 0 4 8 0 1 1 5 6 0 0 0 5 7 0 -1 0 5 8 1 0 0 5 8 1 1 1 6 7 0 -1 0 6 8 1 0 1 6 8 1 1 1 7 8 1 1 0 7 8 1 1 1"
    @test string(juc101[PE]) == "UNKNOWN 3 1 2 0 0 0 1 3 0 0 0 1 4 0 0 0 1 5 0 0 0 2 3 0 0 0 2 6 0 0 0 2 7 0 0 0 3 8 0 0 0 3 9 0 0 0 4 7 0 0 0 4 8 0 0 0 4 10 0 0 0 5 11 0 0 0 5 12 0 0 0 6 11 1 0 0 6 13 0 0 0 7 8 0 0 0 7 14 0 0 0 8 15 0 0 0 9 11 0 1 0 9 16 0 0 0 10 17 0 0 0 10 18 0 0 0 12 13 -1 1 0 12 16 -1 0 0 12 18 0 0 -1 13 16 0 -1 0 13 19 0 0 0 14 17 1 0 0 14 19 0 0 1 15 17 0 1 0 15 20 0 0 0 16 20 0 0 -1 18 19 -1 1 1 18 20 -1 0 0 19 20 0 -1 -1"
    @test juc101[PEM].genome == PeriodicGraph("3 1 2 0 0 0 1 3 0 0 0 1 4 0 0 0 1 5 0 0 0 1 6 0 0 0 1 7 0 0 0 2 4 0 0 0 2 8 0 0 0 3 4 0 0 0 3 9 0 0 0 4 6 0 0 0 4 10 0 0 0 4 11 0 0 0 5 6 0 0 0 5 12 0 0 0 6 7 0 0 0 6 10 0 0 0 6 11 0 0 0 7 13 0 0 0 8 14 0 0 0 8 15 0 0 0 9 16 0 0 0 9 17 0 0 0 10 18 0 0 0 11 19 0 0 0 12 15 -1 0 0 12 20 0 0 0 13 16 -1 0 0 13 21 0 0 0 14 22 0 0 0 14 23 0 0 0 15 18 1 0 -1 16 19 1 0 -1 17 22 -1 1 0 17 23 -1 1 0 18 24 0 0 0 19 25 0 0 0 20 22 -1 0 1 20 26 0 0 0 21 22 -2 1 1 21 26 -1 1 0 22 23 0 0 0 22 26 1 0 -1 23 24 0 0 -1 23 25 1 -1 -1 23 26 1 0 -1 24 26 1 0 0 25 26 0 1 0")

    ewetuw = determine_topology(joinpath(cifs, "EWETUW_clean.cif"); kwargs...)
    @test allunique(ewetuw)
    @test ewetuw[SingleNodes].name == "fit"
    @test unique!(sort!(degree(ewetuw[PE].genome))) == [3, 5, 6]
    @assert allunique(unique!(sort!(degree(x.genome))) for x in ewetuw)

    wemfif = determine_topology(joinpath(cifs, "WEMFIF_clean.cif"); kwargs...)
    @test wemfif[AllNodes] == wemfif[SingleNodes] == wemfif[Standard] == wemfif[PEM]
    @test wemfif[AllNodes].name == "dia"
    @test wemfif[PE].name == "crs"
    CrystalNets.toggle_warning(true)
end

@testset "Archive" begin
    @info "Checking that all known topologies are recognized (this can take a few minutes)."
    if Threads.nthreads() == 1
        @info "Use multiple threads to reduce this time"
    end
    reverse_archive = collect(CrystalNets.CRYSTAL_NETS_ARCHIVE)
    failurelock = ReentrantLock()
    failures = 0
    Threads.@threads for (genome, id) in reverse_archive
        test = try
            topological_genome(CrystalNet(PeriodicGraph(genome))).name == id
        catch e
            CrystalNets.isinterrupt(e) && rethrow()
            false
        end
        if !test
            lock(failurelock) do
                failures += 1
                # The following will throw as non-boolean, hence printing the failing test
                @test "$id failed (Archive)"
            end
        end
    end
    Test.get_testset().n_passed += length(reverse_archive) - failures
end

@testset "Module" begin
    targets = ["pcu", "afy, AFY", "apc, APC", "bam", "bcf", "cdp", "cnd", "ecb", "fiv",
    "ftd", "ftj", "ins", "kgt", "mot", "moz", "muh", "pbz", "qom", "sig",
    "sma", "sod-f", "sod-h", "utj", "utp"]
    failurelock = ReentrantLock()
    failures = 0
    Threads.@threads for target in targets
        @info "Testing $target"
        graph = PeriodicGraph(CrystalNets.REVERSE_CRYSTAL_NETS_ARCHIVE[target])
        n = PeriodicGraphs.nv(graph)
        for k in 1:50
            r = randperm(n)
            offsets = [SVector{3,Int}([rand(-3:3) for _ in 1:3]) for _ in 1:n]
            graph = swap_axes!(offset_representatives!(graph[r], offsets), randperm(3))
            if topological_genome(CrystalNet(graph)).name != target
                lock(failurelock) do
                    failures += 1
                    # The following will throw as non-boolean, hence printing the failing test
                    @test "$target failed (Module) with g = $(string(graph))"
                end
            end
        end
    end

    Test.get_testset().n_passed += length(targets) - failures
    cifs, crystalnetsdir = _finddirs()
    @test topological_genome(CrystalNet(redirect_stderr(devnull) do;
            parse_chemfile(joinpath(cifs, "Moganite.cif"))
          end)).name == "mog"
end

# # The following testset is too long to be run on CI
# @testset "Full-randomization test 3D" begin
#     reverse_archive3D = Tuple{String,String}[(g, id) for (g, id) in CrystalNets.CRYSTAL_NETS_ARCHIVE if g[1] == '3']
#     failurelock = ReentrantLock()
#     failures = 0
#     Threads.@threads for (genome, id) in reverse_archive3D
#         graph = PeriodicGraph3D(genome)
#         n = PeriodicGraphs.nv(graph)
#         for k in 1:30
#             r = randperm(n)
#             offsets = [SVector{3,Int}(rand(-3:3, 3)) for _ in 1:n]
#             graph = swap_axes!(offset_representatives!(graph[r], offsets), randperm(3))
#             topresult = topological_genome(CrystalNet(graph))
#             if topresult.name != id
#                 lock(failurelock) do
#                     failures += 1
#                     # The following will throw as non-boolean, hence printing the failing test
#                     @error """$id failed (full-random$(isempty(topresult.error) ? "" : "; error: "*
#                             Base._truncate_at_width_or_chars(true, topresult.error, 100)))"""
#                 end
#                 break
#             end
#         end
#     end
#     @info "Passed: $(length(reverse_archive3D) - failures)"
# end

@testset "Executable" begin
    cifs, crystalnetsdir = _finddirs()
    safeARGS = deepcopy(ARGS)

    out = tempname()
    empty!(ARGS)

    push!(ARGS, "-k", "3   1 2  0 0 0   1 2  0 0 1   1 2  0 1 0   1 2  1 0 0")
    result, written = capture_out(out)
    @test result == 0
    @test written == ["dia"]

    empty!(ARGS)
    path = joinpath(cifs, "ABW.cif")
    push!(ARGS, path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["sra, ABW"]

    empty!(ARGS)
    path = joinpath(cifs, "ABW.cif")
    push!(ARGS, "-a", joinpath(CrystalNets.arc_location, "rcsr.arc"), path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["sra"]
    __reset_archive!(safeARCHIVE, safeREVERSE)

    empty!(ARGS)
    path = joinpath(cifs, "RRO.cif")
    push!(ARGS, path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["RRO"]

    empty!(ARGS)
    path = joinpath(cifs, "RRO.cif")
    push!(ARGS, "-a", joinpath(CrystalNets.arc_location, "rcsr.arc"), path)
    result, written = capture_out(out)
    @test result == 1
    @test startswith(only(written), "UNKNOWN")
    __reset_archive!(safeARCHIVE, safeREVERSE)

    empty!(ARGS)
    path = joinpath(cifs, "HKUST-1.cif")
    push!(ARGS, "-s", "mof", path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["AllNodes, SingleNodes: tbo"]

    empty!(ARGS)
    path = joinpath(cifs, "HKUST-1_sym.cif")
    push!(ARGS, "-s", "mof", path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["AllNodes, SingleNodes: tbo"]

    empty!(ARGS)
    path = joinpath(cifs, "Diamond.cif")
    push!(ARGS, "-c", "atom", path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["EachVertex: dia"]

    empty!(ARGS)
    path = joinpath(cifs, "Diamond.cif")
    push!(ARGS, path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["dia"]

    empty!(ARGS)
    path = joinpath(cifs, "ALPO-3.10.7.163.001.cif")
    push!(ARGS, "-s", "guess", "-b", "input", path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["AllNodes, SingleNodes: gme, GME"]

    empty!(ARGS)
    path = joinpath(cifs, "ALPO-3.10.7.163.001.cif")
    push!(ARGS, "-s", "guess", "-b", "guess", path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["AllNodes, SingleNodes: gme, GME"]

    empty!(ARGS)
    path = joinpath(cifs, "ALPO-3.1.1.37.001.cif")
    push!(ARGS, "-s", "guess", "-b", "input", path)
    result, written = capture_out(out)
    @test result == 1 # Unknown topology with the input bonds

    empty!(ARGS)
    path = joinpath(cifs, "ALPO-3.1.1.37.001.cif")
    push!(ARGS, "-b", "auto", path)
    result, written = capture_out(out)
    @test_broken result == 0
    @test_broken written == ["afi, AFI"]

    # Test automatic removal of solvent residues and sites with multiple atoms
    empty!(ARGS)
    path = joinpath(cifs, "ALPO-3.1.1.128.001.cif")
    push!(ARGS, path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["sas, SAS"]

    empty!(ARGS)
    path = joinpath(cifs, "ALPO-3.1.2.20.002.cif")
    push!(ARGS, path)
    result, written = capture_out(out)
    @test result == 0
    @test written == ["gis, GIS"]

    empty!(ARGS)
    push!(ARGS, "--help")
    result, written = capture_out(out)
    @test result == 0
    @test startswith(popfirst!(written), "usage: CrystalNets")
    @test !isempty(popfirst!(written))
    @test occursin("CRYSTAL_FILE", popfirst!(written))
    # @test occursin("Form B", popfirst!(written))
    # @test occursin("Form C", popfirst!(written))
    @test isempty(popfirst!(written))
    @test isempty(popfirst!(written))
    @test popfirst!(written) == "Automatic recognition of crystal net topologies."

    empty!(ARGS)
    append!(ARGS, safeARGS)
    if basename(@__DIR__) != "test" # if used with include("runtests.jl")
        CrystalNets._reset_archive!()
    end
end

@testset "Unstable nets" begin
    for n in 2:4
        for m in 0:div(n*(n-1), 2)
            @info "Testing unstable node canonicalization for n = $n; m = $m"
            seen = SimpleGraph[]
            for _ in 1:500
                g = SimpleGraph(n, m; seed=-1)
                any(==(g), seen) && continue
                seencolors = Set{Vector{Int}}()
                seensubnodes = Set{Vector{Vector{Int}}}()
                for _ in 1:300
                    colors = rand(1:n, n)
                    sort!(colors) # allowed since the order of vertices is random
                    colors ∈ seencolors && continue
                    push!(seencolors, colors)
                    subnodes = [Int[] for _ in 1:n]
                    for i in 1:n
                        push!(subnodes[colors[i]], i)
                    end
                    filter!(!isempty, subnodes)
                    length(subnodes) == n && continue
                    subnodes ∈ seensubnodes && continue
                    push!(seensubnodes, subnodes)
                    sig = g[CrystalNets._order_collision(g, subnodes)]
                    for r in permutations(1:n)
                        newcolors = colors[r]
                        newsubnodes = [Int[] for _ in 1:n]
                        for i in 1:n
                            push!(newsubnodes[newcolors[i]], i)
                        end
                        filter!(!isempty, newsubnodes)
                        g2 = g[r]
                        @test sig == g2[CrystalNets._order_collision(g2, newsubnodes)]
                    end
                end
            end
        end
    end
end