text:
@testset "linear problems" begin
    @testset "lp1_test" begin
        #=
        primal
            min -4x_2 - 1
        s.t.
            x_1 >= 3         :y_2
            x_1 + 2x_2 <= 3  :y_3
        dual
            max 3y_2 + 3y_3 - 1
        s.t.
            y_2 >= 0
            y_3 <= 0
            y_2 + y_3 == 0    :x_1
            2y_3 == -4        :x_2
        =#
        primal_model = lp1_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 2
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent())
        @test Set(list_of_cons) == Set(
            [
                (VI, MOI.GreaterThan{Float64})
                (VI, MOI.LessThan{Float64})
                (SAF{Float64}, MOI.EqualTo{Float64})
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.GreaterThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.LessThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 2
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == -1.0
        @test MOI.coefficient.(obj.terms) == [3.0; 3.0]
    end

    @testset "lp7_test compact" begin
        #=
        primal
            min -4x1 -3x2 -1
        s.t.
            2x1 + x2 - 3 <= 0  :y_3
            x1 + 2x2 - 3 <= 0  :y_4
            x1 >= 1            :y_1 (behaves like affine since rhs!=0)
            x2 >= 0            :y_2 (does not appear in compact)
        standard dual
            max 3y_4 + 3y_3 + y_1 - 1
        s.a.
            y_1 + 2y_3 + y_4 = -4 :x_1
            y_2 + y_3 + 2y_4 = -3 :x_2
            y_1 >= 0
            y_2 >= 0
            y_3 <= 0
            y_4 <= 0
        compact dual
            max 3y_4 + 3y_3 + y_1 - 1
        s.a.
            y_1 + 2y_3 + y_4 = -4 :x_1
            - y_3 - 2y_4 >= 3 :x_2
            y_1 >= 0
            y_3 <= 0
            y_4 <= 0
        =#
        primal_model = lp7_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 3
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent())
        @test Set(list_of_cons) == Set(
            [
                (VI, MOI.GreaterThan{Float64})
                (SAF{Float64}, MOI.EqualTo{Float64})
                (SAF{Float64}, MOI.GreaterThan{Float64})
                (VVF, MOI.Nonpositives)
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.GreaterThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VVF,MOI.Nonpositives}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.GreaterThan{Float64}}(),
        ) == 1
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == -1.0
        @test Set(MOI.coefficient.(obj.terms)) == Set([3.0; 1.0; 3.0])

        ci = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.EqualTo{Float64}}(),
        )[]
        f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
        @test Set(MOI.coefficient.(f.terms)) == Set([1.0; 2.0; 1.0])
        @test MOI.constant(f) == 0.0
        @test MOI.get(dual_model, MOI.ConstraintSet(), ci) == MOI.EqualTo(-4.0)
        ci = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.GreaterThan{Float64}}(),
        )[]
        f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
        @test Set(MOI.coefficient.(f.terms)) == Set([-1.0; -2.0])
        @test MOI.constant(f) == 0.0
        @test MOI.get(dual_model, MOI.ConstraintSet(), ci) ==
              MOI.GreaterThan(3.0)
    end

    @testset "lp7_test standard" begin
        #=
        primal
            min -4x1 -3x2 -1
        s.t.
            2x1 + x2 - 3 <= 0  :y_3 (nopos1)
            x1 + 2x2 - 3 <= 0  :y_4 (nopos1)
            x1 >= 1            :y_1 (behaves like affine since rhs!=0)
            x2 >= 0            :y_2 (does not appear in compact)
        standard dual
            max 3y_4 + 3y_3 + y_1 - 1
        s.a.
            y_1 + 2y_3 + y_4 = -4 :x_1
            y_2 + y_3 + 2y_4 = -3 :x_2
            y_1 >= 0
            y_2 >= 0
            y_3 <= 0 (nopos1)
            y_4 <= 0 (nopos1)
        compact dual
            max 3y_4 + 3y_3 + y_1 - 1
        s.a.
            y_1 + 2y_3 + y_4 = -4 :x_1
            - y_3 - 2y_4 >= 3 :x_2
            y_1 >= 0
            y_3 <= 0
            y_4 <= 0
        =#
        primal_model = lp7_test()
        dual = dualize(primal_model, consider_constrained_variables = false)
        dual_model, primal_dual_map = dual.dual_model, dual.primal_dual_map
        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 4
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent())
        @test Set(list_of_cons) == Set(
            [
                (SAF{Float64}, MOI.EqualTo{Float64})
                (VI, MOI.GreaterThan{Float64})
                (VVF, MOI.Nonpositives)
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VVF,MOI.Nonpositives}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.GreaterThan{Float64}}(),
        ) == 2
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 2
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == -1.0
        @test Set(MOI.coefficient.(obj.terms)) == Set([3.0; 1.0; 3.0])

        cis = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{VI,MOI.GreaterThan{Float64}}(),
        )
        for ci in cis
            @test MOI.get(dual_model, MOI.ConstraintSet(), ci) ==
                  MOI.GreaterThan(0.0)
        end

        cis = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.EqualTo{Float64}}(),
        )
        @test Set(MOI.get.(dual_model, MOI.ConstraintSet(), cis)) ==
              Set([MOI.EqualTo(-4.0), MOI.EqualTo(-3.0)])
    end

    @testset "lp14_test_min" begin
        #=
            min -1 +2x1 -3x2 +4x3
        s.t.
            +5x1  -6x2  +7x3 <= 8   (y1)
            -9x1 +10x2 -11x3 >= -12 (y2)
            13x1 -14x2 +15x3 == 16  (y3)
            x1 >= 0 (y4)
            x2 <= 0 (y5)
            x3 free (y6)
        standard dual
            ...
        compact dual
            max 8y1 -12y2 +16y3 - 1
        s.a.
            -5y1  +9y2 -13y3 +2 >= 0 (x1)
            +6y1 -10y2 +14y3 -3 <= 0 (x2)
            -7y1 +11y2 -15y3 +4 == 0 (x3)
            y1 <= 0
            y2 >= 0
            y3 free
        =#
        primal_model = lp14_test(MOI.MIN_SENSE)
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 3
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent())
        @test Set(list_of_cons) == Set(
            [
                (VI, MOI.GreaterThan{Float64})
                (VI, MOI.LessThan{Float64})
                (SAF{Float64}, MOI.EqualTo{Float64})
                (SAF{Float64}, MOI.GreaterThan{Float64})
                (SAF{Float64}, MOI.LessThan{Float64})
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.GreaterThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.LessThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.GreaterThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.LessThan{Float64}}(),
        ) == 1

        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == -1.0
        @test Set(MOI.coefficient.(obj.terms)) == Set([8.0, -12.0, 16])

        ci = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.GreaterThan{Float64}}(),
        )[]
        f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
        @test Set(MOI.coefficient.(f.terms)) == Set([-5.0; 9.0; -13.0])
        @test MOI.constant(f) == 0.0
        @test MOI.get(dual_model, MOI.ConstraintSet(), ci) ==
              MOI.GreaterThan(-2.0)

        ci = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.LessThan{Float64}}(),
        )[]
        f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
        @test Set(MOI.coefficient.(f.terms)) == Set([6.0; -10.0; +14.0])
        @test MOI.constant(f) == 0.0
        @test MOI.get(dual_model, MOI.ConstraintSet(), ci) == MOI.LessThan(3.0)

        # TODO flip these
        # ci = MOI.get(dual_model, MOI.ListOfConstraintIndices{SAF{Float64},MOI.EqualTo{Float64}}())[]
        # f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
        # @test Set(MOI.coefficient.(f.terms)) == Set([-7.0; +11.0; -15.0])
        # @test MOI.constant(f) == 0.0
        # @test MOI.get(dual_model, MOI.ConstraintSet(), ci) == MOI.EqualTo(-4.0)

    end
    @testset "lp14_test_max compact" begin
        #=
            max -1 +2x1 -3x2 +4x3
        s.t.
            +5x1  -6x2  +7x3 <= 8   (y1)
            -9x1 +10x2 -11x3 >= -12 (y2)
            13x1 -14x2 +15x3 == 16  (y3)
            x1 >= 0 (y4)
            x2 <= 0 (y5)
            x3 free (y6)
        standard dual
            min -8y1 +12y2 -16y3 - 1
        s.a.
            -5y1  +9y2 -13y3 -y4 -2 == 0 (x1)
            +6y1 -10y2 +14y3 -y5 +3 == 0 (x2)
            -7y1 +11y2 -15y3 -y6 -4 == 0 (x3)
            y1 <= 0
            y2 >= 0
            y3 free
            y4 >= 0
            y5 <= 0
            y6 == 0
        compact dual
            min -8y1 +12y2 -16y3 - 1
        s.a.
            -5y1  +9y2 -13y3 -2 >= 0 (x1)
            +6y1 -10y2 +14y3 +3 <= 0 (x2)
            -7y1 +11y2 -15y3 -4 == 0 (x3)
            y1 <= 0
            y2 >= 0
            y3 free
        =#
        primal_model = lp14_test(MOI.MAX_SENSE)
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 3
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent())
        @test Set(list_of_cons) == Set(
            [
                (VI, MOI.GreaterThan{Float64})
                (VI, MOI.LessThan{Float64})
                (SAF{Float64}, MOI.EqualTo{Float64})
                (SAF{Float64}, MOI.GreaterThan{Float64})
                (SAF{Float64}, MOI.LessThan{Float64})
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.GreaterThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.LessThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.GreaterThan{Float64}}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.LessThan{Float64}}(),
        ) == 1

        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == -1.0
        @test Set(MOI.coefficient.(obj.terms)) == Set([-8.0, 12.0, -16])

        ci = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.GreaterThan{Float64}}(),
        )[]
        f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
        @test Set(MOI.coefficient.(f.terms)) == Set([-5.0; 9.0; -13.0])
        @test MOI.constant(f) == 0.0
        @test MOI.get(dual_model, MOI.ConstraintSet(), ci) ==
              MOI.GreaterThan(2.0)

        ci = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.LessThan{Float64}}(),
        )[]
        f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
        @test Set(MOI.coefficient.(f.terms)) == Set([6.0; -10.0; +14.0])
        @test MOI.constant(f) == 0.0
        @test MOI.get(dual_model, MOI.ConstraintSet(), ci) == MOI.LessThan(-3.0)

        # TODO flip these
        # ci = MOI.get(dual_model, MOI.ListOfConstraintIndices{SAF{Float64},MOI.EqualTo{Float64}}())[]
        # f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
        # @test Set(MOI.coefficient.(f.terms)) == Set([-7.0; +11.0; -15.0])
        # @test MOI.constant(f) == 0.0
        # @test MOI.get(dual_model, MOI.ConstraintSet(), ci) == MOI.EqualTo(4.0)

    end
    @testset "lp14_test_min standard" begin
        #=
            min -1 +2x1 -3x2 +4x3
        s.t.
            +5x1  -6x2  +7x3 <= 8   (y1)
            -9x1 +10x2 -11x3 >= -12 (y2)
            13x1 -14x2 +15x3 == 16  (y3)
            x1 >= 0 (y4)
            x2 <= 0 (y5)
            x3 free
        standard dual
            min 8y1 -12y2 +16y3 - 1
        s.a.
            -5y1  +9y2 -13y3 -y4 +2 == 0 (x1)
            +6y1 -10y2 +14y3 -y5 -3 == 0 (x2)
            -7y1 +11y2 -15y3     +4 == 0 (x3)
            y1 <= 0
            y2 >= 0
            y3 free
            y4 >= 0
            y5 <= 0
        compact dual
            max 8y1 -12y2 +16y3 - 1
        s.a.
            -5y1  +9y2 -13y3 +2 >= 0 (x1)
            +6y1 -10y2 +14y3 -3 <= 0 (x2)
            -7y1 +11y2 -15y3 +4 == 0 (x3)
            y1 <= 0
            y2 >= 0
            y3 free
        =#
        primal_model = lp14_test(MOI.MIN_SENSE)
        dual = dualize(primal_model, consider_constrained_variables = false)
        dual_model, primal_dual_map = dual.dual_model, dual.primal_dual_map

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 5
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent())
        @test Set(list_of_cons) == Set(
            [
                (SAF{Float64}, MOI.EqualTo{Float64})
                (VI, MOI.GreaterThan{Float64})
                (VI, MOI.LessThan{Float64})
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 3
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.GreaterThan{Float64}}(),
        ) == 2
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.LessThan{Float64}}(),
        ) == 2

        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == -1.0
        @test Set(MOI.coefficient.(obj.terms)) == Set([8.0, -12.0, 16])

        cis = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{VI,MOI.GreaterThan{Float64}}(),
        )
        for ci in cis
            @test MOI.get(dual_model, MOI.ConstraintSet(), ci) ==
                  MOI.GreaterThan(0.0)
        end
        cis = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{VI,MOI.LessThan{Float64}}(),
        )
        for ci in cis
            @test MOI.get(dual_model, MOI.ConstraintSet(), ci) ==
                  MOI.LessThan(0.0)
        end

        cis = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.EqualTo{Float64}}(),
        )
        # TODO: review signs here
        @test Set(MOI.get.(dual_model, MOI.ConstraintSet(), cis)) ==
              Set([MOI.EqualTo(2.0), MOI.EqualTo(-3.0), MOI.EqualTo(4.0)])
        for ci in cis
            set = MOI.get(dual_model, MOI.ConstraintSet(), ci)
            f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
            if set == MOI.EqualTo(2.0)
                @test Set(MOI.coefficient.(f.terms)) ==
                      Set(-[-5.0; 9.0; -13.0; -1])
                @test MOI.constant(f) == 0.0
            elseif set == MOI.EqualTo(-3.0)
                @test Set(MOI.coefficient.(f.terms)) ==
                      Set(-[6.0; -10.0; +14.0; -1])
                @test MOI.constant(f) == 0.0
            elseif set == MOI.EqualTo(4.0)
                @test Set(MOI.coefficient.(f.terms)) ==
                      Set(-[-7.0; +11.0; -15.0])
                @test MOI.constant(f) == 0.0
            else
                error("wrong ci")
            end
        end
    end
    @testset "lp14_test_max standard" begin
        #=
            max -1 +2x1 -3x2 +4x3
        s.t.
            +5x1  -6x2  +7x3 <= 8   (y1)
            -9x1 +10x2 -11x3 >= -12 (y2)
            13x1 -14x2 +15x3 == 16  (y3)
            x1 >= 0 (y4)
            x2 <= 0 (y5)
            x3 free
        standard dual
            min -8y1 +12y2 -16y3 - 1
        s.a.
            -5y1  +9y2 -13y3 -y4 -2 == 0 (x1)
            +6y1 -10y2 +14y3 -y5 +3 == 0 (x2)
            -7y1 +11y2 -15y3     -4 == 0 (x3)
            y1 <= 0
            y2 >= 0
            y3 free
            y4 >= 0
            y5 <= 0
        compact dual
            min -8y1 +12y2 -16y3 - 1
        s.a.
            -5y1  +9y2 -13y3 -2 >= 0 (x1)
            +6y1 -10y2 +14y3 +3 <= 0 (x2)
            -7y1 +11y2 -15y3 -4 == 0 (x3)
            y1 <= 0
            y2 >= 0
            y3 free
        =#
        primal_model = lp14_test(MOI.MAX_SENSE)
        dual = dualize(primal_model, consider_constrained_variables = false)
        dual_model, primal_dual_map = dual.dual_model, dual.primal_dual_map

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 5
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent())
        @test Set(list_of_cons) == Set(
            [
                (SAF{Float64}, MOI.EqualTo{Float64})
                (VI, MOI.GreaterThan{Float64})
                (VI, MOI.LessThan{Float64})
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 3
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.GreaterThan{Float64}}(),
        ) == 2
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VI,MOI.LessThan{Float64}}(),
        ) == 2

        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.constant(obj) == -1.0
        @test Set(MOI.coefficient.(obj.terms)) == Set([-8.0, 12.0, -16])

        cis = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{VI,MOI.GreaterThan{Float64}}(),
        )
        for ci in cis
            @test MOI.get(dual_model, MOI.ConstraintSet(), ci) ==
                  MOI.GreaterThan(0.0)
        end
        cis = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{VI,MOI.LessThan{Float64}}(),
        )
        for ci in cis
            @test MOI.get(dual_model, MOI.ConstraintSet(), ci) ==
                  MOI.LessThan(0.0)
        end

        cis = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.EqualTo{Float64}}(),
        )
        # TODO: review signs here
        @test Set(MOI.get.(dual_model, MOI.ConstraintSet(), cis)) ==
              Set([MOI.EqualTo(-2.0), MOI.EqualTo(3.0), MOI.EqualTo(-4.0)])
        for ci in cis
            set = MOI.get(dual_model, MOI.ConstraintSet(), ci)
            f = MOI.get(dual_model, MOI.ConstraintFunction(), ci)
            if set == MOI.EqualTo(-2.0)
                @test Set(MOI.coefficient.(f.terms)) ==
                      Set(-[-5.0; 9.0; -13.0; -1])
                @test MOI.constant(f) == 0.0
            elseif set == MOI.EqualTo(3.0)
                @test Set(MOI.coefficient.(f.terms)) ==
                      Set(-[6.0; -10.0; +14.0; -1])
                @test MOI.constant(f) == 0.0
            elseif set == MOI.EqualTo(-4.0)
                @test Set(MOI.coefficient.(f.terms)) ==
                      Set(-[-7.0; +11.0; -15.0])
                @test MOI.constant(f) == 0.0
            else
                error("wrong ci")
            end
        end
    end
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/nonconvex.jl
# Functions which determine when the upper bounding (local nlp) problem should
# be solved as well as routines used to call the upper bounding problem.
#############################################################################

"""
$(TYPEDSIGNATURES)

Checks thats the solution of the lower (relaxed problem) is integer feasible to
within tolerances specified by the parameters: `integer_abs_tol` (absolute tolerance)
and `integer_rel_tol` (relative tolerance).
"""
function is_integer_feasible_relaxed(m::GlobalOptimizer)
    bool = true
    atol = _integer_abs_tol(m)
    rtol = _integer_rel_tol(m)
    for i = 1:_variable_num(BranchVar(), m)
        if is_integer(BranchVar(), m, i)
            xsol = _lower_solution(BranchVar(), m, i)
            if isapprox(floor(xsol), xsol; atol = atol, rtol = rtol)
                continue
            elseif isapprox(ceil(xsol), xsol; atol = atol, rtol = rtol)
                continue
            else
                bool &= false
                break
            end
        end
    end
    return bool
end

"""
$(TYPEDSIGNATURES)

Default check to see if the upper bounding problem should be run. By default,
The upper bounding problem is run on every node up to depth `upper_bounding_depth`
and is triggered with a probability of `0.5^(depth - upper_bounding_depth)`
afterwards for continuous problems. For integral problems, the `upper_bounding_depth`
the above approach is used as well as running on every node up to depth
`upper_bounding_depth + cont_depth` and is triggered with a probability of
`0.5^(depth - upper_bounding_depth - cont_depth)`.
"""
function default_upper_heurestic(m::GlobalOptimizer)
    bool = false
    ubd_limit = m._parameters.upper_bounding_depth
    n = _current_node(m)
    if is_integer_feasible_relaxed(m)
        Δdepth = n.depth - n.cont_depth
        bool |= (Δdepth <= ubd_limit)
        bool |= (rand() < 0.5^(Δdepth - ubd_limit))
    end
    bool |= (n.depth <= ubd_limit)
    bool |= (rand() < 0.5^(n.depth - ubd_limit))
    return bool
end

"""
$(TYPEDSIGNATURES)

Default upper bounding problem which simply calls `solve_local_nlp!` to solve
the nlp locally.
"""
function upper_problem!(t::ExtensionType, m::GlobalOptimizer)
    if !default_upper_heurestic(m)
        m._upper_feasibility = false
        m._upper_objective_value = Inf
    else
        solve_local_nlp!(m)
    end
    return
end
upper_problem!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = upper_problem!(_ext(m), m)

ext: jl

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using MetaGraphs: MetaDiGraph, set_prop!, props

#####################
# general parser infrastructure for circuits
#####################

# The `ParserCombinator` library works correctly but is orders of magnitude too slow.
# Instead here we hardcode some simpler parsers to speed things up

"""
Load a probabilistic circuit from file.
Support circuit file formats:
 * ".psdd" for PSDD files
"""
function load_prob_circuit(file::String)::ProbΔ
    @assert endswith(file,".psdd")
    compile_prob(parse_psdd_file(file))
end

"""
Load a structured probabilistic circuit from file.
Support circuit file formats:
 * ".psdd" for PSDD files
Supported vtree file formats:
 * ".vtree" for VTree files
"""
function load_struct_prob_circuit(circuit_file::String, vtree_file::String)::Tuple{ProbΔ,PlainVtree}
    @assert endswith(circuit_file,".psdd")
    circuit_lines = parse_circuit_file(circuit_file)
    vtree_lines = parse_vtree_file(vtree_file)
    compile_struct_prob(circuit_lines, vtree_lines)
end

function load_logistic_circuit(circuit_file::String, classes::Int)::LogisticΔ
    @assert endswith(circuit_file,".circuit")
    circuit_lines = parse_circuit_file(circuit_file)
    compile_logistic(circuit_lines, classes)
end

#####################
# parse based on file extension
#####################

function parse_circuit_file(file::String)::CircuitFormatLines
    if endswith(file,".circuit")
        parse_lc_file(file)
    elseif endswith(file,".psdd")
        parse_psdd_file(file)
    elseif endswith(file,".sdd")
        parse_sdd_file(file)
    else
        throw("Cannot parse this file type as a circuit: $file")
    end
end

"Parse a clt from given file"
function parse_clt(filename::String)::MetaDiGraph
    f = open(filename)
    n = parse(Int32,readline(f))
    n_root = parse(Int32,readline(f))
    clt = MetaDiGraph(n)
    for i in 1 : n_root
        root, prob = split(readline(f), " ")
        root, prob = parse(Int32, root), parse(Float64, prob)
        set_prop!(clt, root, :parent, 0)
        set_prop!(clt, root, :cpt, Dict(1=>prob,0=>1-prob))
    end

    for i = 1 : n - n_root
        dst, src, prob1, prob0 = split(readline(f), " ")
        dst, src, prob1, prob0 = parse(Int32, dst), parse(Int32, src), parse(Float64, prob1), parse(Float64, prob0)
        add_edge!(clt, src,dst)
        set_prop!(clt, dst, :parent, src)
        set_prop!(clt, dst, :cpt, Dict((1,1)=>prob1, (0,1)=>1-prob1, (1,0)=>prob0, (0,0)=>1-prob0))
    end
    return clt
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Conduit_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Conduit")
JLLWrappers.@generate_main_file("Conduit", UUID("36ec8eac-713e-50f0-b9f1-e7516df2c55e"))
end  # module Conduit_jll


text:
export
FactorGraph,
currentGraph,
setCurrentGraph,
nodes,
edges,
Terminal,
ultimatePartner

"""
A factor graph consisting of factor nodes and edges.
"""
mutable struct FactorGraph
    nodes::Dict{Symbol, FactorNode}
    edges::Vector{Edge}
    variables::Dict{Symbol, Variable}
    counters::Dict{String, Int} # Counters for automatic node id assignments
    placeholders::Dict{Clamp, Tuple{Symbol, Int}}
end

"""
Return currently active `FactorGraph`.
Create one if there is none.
"""
function currentGraph()
    try
        return current_graph
    catch
        return FactorGraph()
    end
end

setCurrentGraph(graph::FactorGraph) = global current_graph = graph

function FactorGraph()

    global current_posterior_factorization = nothing

    setCurrentGraph(FactorGraph(Dict{Symbol, FactorNode}(),
                                            Edge[],
                                            Dict{Symbol, Variable}(),
                                            Dict{Type, Int}(),
                                            Dict{Clamp, Tuple{Symbol, Int}}()))
end
"""
Automatically generate a unique id based on the current counter value for the element type.
"""
function generateId(t::Type)
    graph = currentGraph()
    tname = lowercase(String(nameof(t))) # Remove module prefix from typename
    counter = haskey(graph.counters, tname) ? graph.counters[tname]+1 : 1
    id = Symbol("$(tname)_$(counter)")

    # Make sure we have a unique id (if we can check it)
    collection =
        if t <: FactorNode
            graph.nodes
        elseif t <: Edge
            graph.edges
        elseif t <: Variable
            graph.variables
        end
    if collection != nothing
        # Increase counter until we have a unique id
        while haskey(collection, id)
            counter += 1
            id = Symbol("$(tname)_$(counter)")
        end
    end

    graph.counters[tname] = counter # save counter

    return id
end

"""
Add a FactorNode to a FactorGraph
"""
function addNode!(graph::FactorGraph, nd::FactorNode)
    !haskey(graph.nodes, nd.id) || error("Graph already contains a FactorNode with id $(nd.id)")
    graph.nodes[nd.id] = nd
    return graph
end

"""
Add a Variable to a FactorGraph
"""
function addVariable!(graph::FactorGraph, var::Variable)
    !haskey(graph.variables, var.id) || error("Graph already contains a Variable with id $(var.id)")
    graph.variables[var.id] = var
    return graph
end

"""
`hasNode(graph, node)` checks if `node` is part of `graph`.
"""
hasNode(graph::FactorGraph, nd::FactorNode) = (haskey(graph.nodes, nd.id) && (graph.nodes[nd.id] === nd))

"""
`hasVariable(graph, var)` checks if `var` is part of `graph`.
"""
hasVariable(graph::FactorGraph, var::Variable) = (haskey(graph.variables, var.id) && (graph.variables[var.id] === var))

nodes(graph::FactorGraph = currentGraph()) = Set{FactorNode}(values(graph.nodes))

function nodes(edgeset::Set{Edge})
    # Return all nodes connected to edgeset
    connected_nodes = Set{FactorNode}()
    for edge in edgeset
        (edge.a == nothing) || push!(connected_nodes, edge.a.node)
        (edge.b == nothing) || push!(connected_nodes, edge.b.node)
    end

    return connected_nodes
end

edges(graph::FactorGraph=currentGraph()) = Set{Edge}(graph.edges)
edges(node::FactorNode) = Set{Edge}([intf.edge for intf in node.interfaces])
edges(nodeset::Set{FactorNode}) = union(Set((edges(node) for node=nodeset))...)

"""
Description:

    Terminal is a special type of node that is only used in the internal
    graph of a CompositeFactor. A Terminal is used to terminate an Edge in the
    internal graph that is linked to an interface of the CompositeFactor.

    A Terminal is linked to an interface of the
    CompositeFactor containing the Terminal.

Interfaces:

    1. out

Construction:

    Terminal(id=:some_id)
"""
mutable struct Terminal <: FactorNode
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}
    outer_interface::Interface # Interface of CompositeFactor linked to this Terminal

    function Terminal(out::Variable, outer_interface::Interface; id=generateId(Terminal))
        self = new(id, Array{Interface}(undef, 1), Dict{Symbol,Interface}(), outer_interface)
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)

        return self
    end
end

"""
ultimatePartner(interface) finds the 'ultimate partner' of interface.
If interface.partner does not belong to a Terminal, it simply returns
interface.partner. In case of a Terminal node, it finds the first
non-Terminal partner on a higher level factor graph.
"""
function ultimatePartner(interface::Interface)
    if (interface.partner != nothing) && isa(interface.partner.node, Terminal)
        return ultimatePartner(interface.partner.node.outer_interface)
    else
        return interface.partner
    end
end

"""
Deterministic edges are associated with a PointMass belief.
Find all edges in the graph that are deterministic.
"""
function deterministicEdges(fg::FactorGraph)
    deterministic_edges = Set{Edge}()
    for node in nodes(fg)
        if isDeltaConstraint(node)
            union!(deterministic_edges, deterministicEdgeSet(node.i[:out].edge))
        end
    end

    return deterministic_edges
end

"""
Extend the root_edge to a deterministically connected set, and determine
which edges in the extended set carry PointMass beliefs.
"""
function deterministicEdgeSet(root_edge::Edge)
    # Collect deterministically connected edges
    edge_set = extend(root_edge, terminate_at_soft_factors=true)

    # Collect all interfaces in edge set
    interfaces = Interface[]
    for edge in edge_set
        (edge.a == nothing) || push!(interfaces, edge.a)
        (edge.b == nothing) || push!(interfaces, edge.b)
    end

    # Identify a schedule that propagates through the entire edge set
    dg = summaryDependencyGraph(edge_set)
    schedule = children(interfaces, dg)

    # Build dictionary of deterministicness
    is_deterministic = Dict{Interface, Bool}()
    for interface in schedule
        is_deterministic[interface] = isDeterministic(interface, is_deterministic)
    end

    # Determine deterministic edges
    deterministic_edges = Set{Edge}()
    for edge in edge_set
        deterministic_a = (edge.a != nothing) && is_deterministic[edge.a]
        deterministic_b = (edge.b != nothing) && is_deterministic[edge.b]

        # An edge is deterministic if any of its interfaces are deterministic
        if deterministic_a || deterministic_b
            push!(deterministic_edges, edge)
        end
    end

    return deterministic_edges
end

"""
Determine whether interface sends a PointMass message, based on the connected
node, and the deterministic status of the inbound messages.
"""
function isDeterministic(interface::Interface, is_deterministic::Dict{Interface, Bool})
    # First handle immediate cases
    node = interface.node
    if !isa(node, DeltaFactor)
        # We assume a non-delta factor sends stochastic messages in all directions.
        # This includes the CompositeFactor, for which this is the conservative choice.
        return false
    elseif isDeltaConstraint(node)
        return true # Clamp nodes only send deterministic messages
    end

    # Collect deterministicness of inbounds
    inbounds_deterministic = Bool[]
    for iface in node.interfaces
        if iface != interface
            partner = ultimatePartner(iface)
            if partner == nothing # Dangling edge
                push!(inbounds_deterministic, false) # Unconstrained inbound is considered stochastic (uninformative)
            else
                push!(inbounds_deterministic, is_deterministic[partner])
            end
        end
    end

    # Determine deterministicness of outbound from inbounds
    if isa(node, Equality) && any(inbounds_deterministic)
        return true # Deterministic
    elseif any(.!inbounds_deterministic)
        return false # Delta factor with any stochastic inbounds sends a stochastic message
    elseif all(inbounds_deterministic)
        return true # Delta factor with all inbounds deterministic sends a deterministic message
    else # Fallback
        return false # Stochastic
    end
end


text:
@testset "Utilities" begin
  a, b, c = P2(0,0), P2(1,0), P2(0,1)
  @test signarea(a, b, c) == T(0.5)
  a, b, c = P2(0,0), P2(0,1), P2(1,0)
  @test signarea(a, b, c) == T(-0.5)

  p1, p2, p3 = P2(0,0), P2(1,1), P2(0.25,0.5)
  s = Segment(P2(0.5,0.0), P2(0.0,1.0))
  @test sideof(p1, s) == :LEFT
  @test sideof(p2, s) == :RIGHT
  @test sideof(p3, s) == :ON

  p1, p2, p3 = P2(0.5,0.5), P2(1.5,0.5), P2(1,1)
  c = Chain(P2[(0,0),(1,0),(1,1),(0,1),(0,0)])
  @test sideof(p1, c) == :INSIDE
  @test sideof(p2, c) == :OUTSIDE
  @test sideof(p3, c) == :INSIDE
end


text:
module COFILS

# package code goes here

using Persa
using DecisionTree

include("latentvariable.jl")
include("method.jl")

end # module


text:
using CSV
using Statistics
using SurNoR_2020
using Distributions
using Distributed
using DataFrames
using Random

# ------------------------------------------------------------------------------
# Reading Data and making input
# ------------------------------------------------------------------------------
BehavData = CSV.read("./data/BehavData_SwitchState.csv")

# ------------------------------------------------------------------------------
# Going through updated rules
# ------------------------------------------------------------------------------
for update_form = ["Leaky"]
        if update_form == "SMiLe"
            surprise_form="CC"
        else
            surprise_form="BF"
        end

        novelty_form="log-count"
        update_form_model_free = "Q_learning"
        gamma_form="fix"

        # ----------------------------------------------------------------------
        # Defining Hyper Parameters
        HyperParam = SurNoR_2020.Str_HyperParam(; gamma_form=gamma_form,
                                                surprise_form=surprise_form,
                                                novelty_form=novelty_form,
                                                update_form=update_form,
                                                update_form_model_free=update_form_model_free)

        # ----------------------------------------------------------------------
        # Defining Parameters
        Param1 = SurNoR_2020.Str_Param(; num_act=4, num_obs=11,
                                        count_decay=1, Q_R_0 = 0, Δα=0)

        Param2 = deepcopy(Param1)

        Param = [Param1, Param2]

        # ----------------------------------------------------------------------
        # Subject Information
        Block_Set = 1:2
        Fold_Set = [[1,2,3,4,5,6,7,8],
                    [1,2,3,4,9,11,12,13],
                    [5,6,7,8,9,11,12,13]]
        Epi_Set = 1:5

        # ----------------------------------------------------------------------
        # Agents
        T_Iter = 40
        for Chain = 1:100
            @distributed for Fold = 1:3
                    Random.seed!((Fold-1)*100 + Chain)

                    Sub_Set = Fold_Set[Fold]
                    LogL_Seq = zeros(T_Iter, 1)
                    Param_Seq = zeros(T_Iter, 17)                                                                            # WARNING: SHOULD CHANGE

                    Param_Seq[1,1] = 0.0 + 0.01*rand(1)[1]             # pc
                    Param_Seq[1,2] = 1e-6 + (9e-6)*rand(1)[1]          # ϵ
                    Param_Seq[1,3] = rand(Rayleigh(5) ,1)[1]           # β_R
                    Param_Seq[1,4] = rand(Rayleigh(0.05) ,1)[1]        # β_RN1
                    Param_Seq[1,5] = rand(Rayleigh(0.05) ,1)[1]        # β_RN2
                    Param_Seq[1,6] = 0.8 + 0.19*rand(1)[1]             # γ_R
                    Param_Seq[1,7] = 0.4 + 0.59*rand(1)[1]             # γ_N
                    Param_Seq[1,8] = 0.8 + 0.19*rand(1)[1]             # λ_R
                    Param_Seq[1,9] = 0.4 + 0.59*rand(1)[1]             # λ_N
                    Param_Seq[1,10] = 0.0 + 0.19*rand(1)[1]            # α
                    Param_Seq[1,11] = 0 + 4*rand(1)[1]                 # Q_N_0
                    Param_Seq[1,12] = 20                               # N_PS
                    Param_Seq[1,13] = 0 + 1*rand(1)[1]                 # w_2
                    Param_Seq[1,14] = 0 + 1*rand(1)[1]                 # w_11
                    Param_Seq[1,15] = 0 + 1*rand(1)[1]                 # w_12
                    Param_Seq[1,16] = 0 + 10*rand(1)[1]                # w_scaler
                    Param_Seq[1,17] = rand(Rayleigh(5) ,1)[1]          # β_R_2

                    Param1.p_change = Param_Seq[1,1]
                    Param1.ϵ = Param_Seq[1,2]
                    Param1.β_R = Param_Seq[1,3]
                    Param1.β_ratio_N2R = Param_Seq[1,4]
                    Param1.γ_RL_R = Param_Seq[1,6]
                    Param1.γ_RL_N = Param_Seq[1,7]
                    Param1.λ_elig_R = Param_Seq[1,8]
                    Param1.λ_elig_N = Param_Seq[1,9]
                    Param1.α = Param_Seq[1,10]
                    Param1.Q_N_0 = Param_Seq[1,11]
                    Param1.N_prior_sweep = Param_Seq[1,12]
                    Param1.w_MF_to_MB_2 = Param_Seq[1,13]
                    Param1.w_MF_to_MB_1 = Param_Seq[1,14]
                    Param1.w_MF_to_MB_Scaler = Param_Seq[1,16]

                    Param2 = deepcopy(Param1)
                    Param2.β_ratio_N2R = Param_Seq[1,5]
                    Param2.w_MF_to_MB_1 = Param_Seq[1,15]
                    Param2.β_R = Param_Seq[1,17]

                    Param = [Param1, Param2]

                    (Agents, LogL) = SurNoR_2020.Func_train_set_of_Agents(BehavData;
                                                                           Sub_Set = Sub_Set,
                                                                           Epi_Set=Epi_Set,
                                                                           Block_Set=Block_Set,
                                                                           Param=deepcopy(Param),
                                                                           HyperParam=HyperParam)
                    LogL_Seq[1] = LogL
                    start = time()
                    Flag = true
                    t = 2
                    while (t<=T_Iter)&&Flag
                            if mod(t,1)==0
                                    @show "-----------------------------"
                                    @show [update_form, Chain, Fold, t]
                                    elapsed = time() - start
                                    @show elapsed
                                    start = time()
                                    @show "-----------------------------"
                            end
                            Param_temp = Param_Seq[t-1,:]

                            Param_perm_order = shuffle(1:17)                                                                    # WARNING: SHOULD CHANGE
                            for i = 1:length(Param_perm_order)
                                    num_param = Param_perm_order[i]
                                    if num_param==1
                                            # P_change
                                            @show "p_change"
                                            x = 0.00:0.0025:0.3
                                    elseif num_param==2
                                            # Eps
                                            @show "eps"
                                            x = -12:0.2:0.2
                                            x = 10 .^ x
                                    elseif num_param==3
                                            # β_R
                                            @show "beta_R"
                                            x = 0.1:0.2:15
                                    elseif num_param==4
                                            # β_RN1
                                            @show "beta_RN1"
                                            x = 0.00:0.005:0.6
                                    elseif num_param==5
                                            # β_RN2
                                            @show "beta_RN2"
                                            x = 0.00:0.005:0.6
                                    elseif num_param==6
                                            # γ_R
                                            @show "γ_R"
                                            x = 0.85:0.01:0.99
                                    elseif num_param==7
                                            # γ_N
                                            @show "γ_N"
                                            x = 0.20:0.05:0.99
                                    elseif num_param==8
                                            # λ_R
                                            @show "λ_R"
                                            x = 0.85:0.01:1
                                    elseif num_param==9
                                            # λ_N
                                            @show "λ_N"
                                            x = 0.20:0.05:1
                                    elseif num_param==10
                                            # α
                                            @show "α"
                                            x = 0.0:0.02:0.2
                                    elseif num_param==11
                                            # Q_N_0
                                            @show "Q_N_0"
                                            x = 0.0:0.25:12
                                    elseif num_param==12
                                            # N_PS
                                            @show "N_PS"
                                            x = cat(0:14,15:5:40,dims=1)
                                    elseif num_param==13
                                            # w_2
                                            @show "w_2"
                                            x = 0:0.05:1
                                    elseif num_param==14
                                            # w_11
                                            @show "w_11"
                                            x = 0:0.05:1
                                    elseif num_param==15
                                            # w_12
                                            @show "w_12"
                                            x = 0:0.05:1
                                    elseif num_param==16
                                            # w_scaler
                                            @show "w_scaler"
                                            x = 0:0.1:6
                                    elseif num_param==17
                                            # β_R_2
                                            @show "beta_R_2"
                                            x = 0.1:0.2:15
                                    end

                                    LogL_range = zeros(size(x))
                                    for x_ind = 1:length(x)
                                            Param_temp[num_param] = x[x_ind]

                                            Param1.p_change = Param_temp[1]
                                            Param1.ϵ = Param_temp[2]
                                            Param1.β_R = Param_temp[3]
                                            Param1.β_ratio_N2R = Param_temp[4]
                                            Param1.γ_RL_R = Param_temp[6]
                                            Param1.γ_RL_N = Param_temp[7]
                                            Param1.λ_elig_R = Param_temp[8]
                                            Param1.λ_elig_N = Param_temp[9]
                                            Param1.α = Param_temp[10]
                                            Param1.Q_N_0 = Param_temp[11]
                                            Param1.N_prior_sweep = Param_temp[12]
                                            Param1.w_MF_to_MB_2 = Param_temp[13]
                                            Param1.w_MF_to_MB_1 = Param_temp[14]
                                            Param1.w_MF_to_MB_Scaler = Param_temp[16]

                                            Param2 = deepcopy(Param1)
                                            Param2.β_ratio_N2R = Param_temp[5]
                                            Param2.w_MF_to_MB_1 = Param_temp[15]
                                            Param2.β_R = Param_temp[17]

                                            Param = [Param1, Param2]

                                            (Agents, LogL) = SurNoR_2020.Func_train_set_of_Agents(BehavData;
                                                                                Sub_Set = Sub_Set,
                                                                                Epi_Set=Epi_Set,
                                                                                Block_Set=Block_Set,
                                                                                Param=deepcopy(Param),
                                                                                HyperParam=HyperParam)
                                            LogL_range[x_ind] = LogL
                                    end
                                    if sum(isnan.(LogL_range)) > 0
                                            @show LogL_range
                                            LogL_range[isnan.(LogL_range)] .= -Inf
                                    end
                                    (LogL, x_ind) = findmax(LogL_range)
                                    Param_temp[num_param] = x[x_ind]
                            end

                            Param_Seq[t,:] = Param_temp
                            LogL_Seq[t] = LogL
                            @show LogL

                            Flag = (!(LogL_Seq[t-1] == LogL))&&(!(LogL == -Inf))

                            if (t==T_Iter)||(mod(t,5)==0)||(!Flag)                                                              # WARNING: SHOULD CHANGE
                                    Path = string("src/Optimization/CV_3Fold/CoAscData_HybN_F", string(Fold), "_C",
                                                    string(Chain), ".csv")
                                    Saved_Data = DataFrame(LogL = LogL_Seq[:],
                                                           p_change = Param_Seq[:,1],
                                                           eps = Param_Seq[:,2],
                                                           beta_R_1 = Param_Seq[:,3],
                                                           beta_Ratio_1 = Param_Seq[:,4],
                                                           beta_Ratio_2 = Param_Seq[:,5],
                                                           gamma_R = Param_Seq[:,6],
                                                           gamma_N = Param_Seq[:,7],
                                                           lambda_R = Param_Seq[:,8],
                                                           lambda_N = Param_Seq[:,9],
                                                           alpha = Param_Seq[:,10],
                                                           Q_N_0 = Param_Seq[:,11],
                                                           N_PS = Param_Seq[:,12],
                                                           w2 = Param_Seq[:,13],
                                                           w11 = Param_Seq[:,14],
                                                           w12 = Param_Seq[:,15],
                                                           w_scaler = Param_Seq[:,16],
                                                           beta_R_2 = Param_Seq[:,17])
                                   CSV.write(Path, Saved_Data)
                            end
                            t += 1
                    end
            end
        end
end


text:
struct IdealGasCpEquation1{T<:Number}
    vals::Array{T,1}
end

function thermocalc(eq::IdealGasCpEquation1, T)
    A,B,C,D,E = eq.vals
    return A+B*T+C*T^2+D*T^3+E/(T^2)
end
function thermocalc(eq::IdealGasCpEquation1, P, T)
    return thermocalc(eq, T)
end


text:
#
# Object headers
#

const HM_NIL = 0x00
const HM_DATASPACE = 0x01
const HM_LINK_INFO = 0x02
const HM_DATATYPE = 0x03
const HM_FILL_VALUE_OLD = 0x04
const HM_FILL_VALUE = 0x05
const HM_LINK_MESSAGE = 0x06
const HM_EXTERNAL_FILE_LIST = 0x07
const HM_DATA_LAYOUT = 0x08
const HM_BOGUS = 0x09
const HM_GROUP_INFO = 0x0a
const HM_FILTER_PIPELINE = 0x0b
const HM_ATTRIBUTE = 0x0c
const HM_OBJECT_COMMENT = 0x0d
const HM_SHARED_MESSAGE_TABLE = 0x0f
const HM_OBJECT_HEADER_CONTINUATION = 0x10
const HM_SYMBOL_TABLE = 0x11
const HM_MODIFICATION_TIME = 0x12
const HM_BTREE_K_VALUES = 0x13
const HM_DRIVER_INFO = 0x14
const HM_ATTRIBUTE_INFO = 0x15
const HM_REFERENCE_COUNT = 0x16

const OH_ATTRIBUTE_CREATION_ORDER_TRACKED = 2^2
const OH_ATTRIBUTE_CREATION_ORDER_INDEXED = 2^3
const OH_ATTRIBUTE_PHASE_CHANGE_VALUES_STORED = 2^4
const OH_TIMES_STORED = 2^5

const OBJECT_HEADER_CONTINUATION_SIGNATURE = htol(0x4b48434f) # "OCHK"

struct ObjectStart
    signature::UInt32
    version::UInt8
    flags::UInt8
end
ObjectStart(flags::UInt8) = ObjectStart(OBJECT_HEADER_SIGNATURE, 2, flags)
define_packed(ObjectStart)

# Reads the start of an object including the signature, version, flags,
# and (payload) size. Returns the size.
function read_obj_start(io::IO)
    os = jlread(io, ObjectStart)
    os.signature == OBJECT_HEADER_SIGNATURE || throw(InvalidDataException())
    os.version == 2 || throw(UnsupportedVersionException())

    if (os.flags & OH_TIMES_STORED) != 0
        # Skip access, modification, change and birth times
        skip(io, 128)
    end
    if (os.flags & OH_ATTRIBUTE_PHASE_CHANGE_VALUES_STORED) != 0
        # Skip maximum # of attributes fields
        skip(io, 32)
    end

    read_size(io, os.flags)
end

struct HeaderMessage
    msg_type::UInt8
    size::UInt16
    flags::UInt8
end
define_packed(HeaderMessage)

function isgroup(f::JLDFile, roffset::RelOffset)
    io = f.io
    seek(io, fileoffset(f, roffset))

    sz = read_obj_start(io)
    pmax::Int64 = position(io) + sz
    while position(io) <= pmax-4
        msg = jlread(io, HeaderMessage)
        endpos = position(io) + msg.size
        if msg.msg_type == HM_LINK_INFO || msg.msg_type == HM_GROUP_INFO || msg.msg_type == HM_LINK_MESSAGE
            return true
        elseif msg.msg_type == HM_DATASPACE || msg.msg_type == HM_DATATYPE || msg.msg_type == HM_FILL_VALUE || msg.msg_type == HM_DATA_LAYOUT
            return false
        end
        seek(io, endpos)
    end
    return false
end


text:
using TensorNetworkEvolve, Random, Test
using Graphs, Yao, ForwardDiff
using TensorNetworkEvolve.TensorAD

@testset "sr" begin
    function rand_hamiltonian(g::SimpleGraph)
        blocks = []
        nbit = nv(g)
        for edge in edges(g)
            i,j = edge.src, edge.dst
            push!(blocks, randn()*kron(nbit,i=>X,j=>X) + randn()*kron(nbit, i=>Y,j=>Y) + randn()*kron(nbit, i=>Z, j=>Z))
        end
        for i in vertices(g)
            push!(blocks, put(nbit,(i,)=>rand()*X))
        end
        return +(blocks...)
    end
    g = SimpleGraph(5)
    for (i,j) in [(1,2), (1,3), (2,4), (2,5), (3,4), (3,5)]
        add_edge!(g, i, j)
    end

    h = rand_hamiltonian(g)
    p1 = rand_simplepeps(ComplexF64, g, 2; Dmax=4)
    vars = TensorAD.DiffTensor(variables(p1))
    gvars = TensorAD.gradient(vars->norm(load_variables(p1, vars)), vars)[1]
    @test gvars isa DiffTensor

    @test TensorNetworkEvolve.iloss2(h, p1, variables(p1)) isa Array{T,0} where T
    fvec = TensorNetworkEvolve.fvec(p1, h)
    @test fvec isa DiffTensor
    _complex(x::AbstractVector) = [Complex(x[2i-1], x[2i]) for i=1:length(x)÷2]
    fvec2 = ForwardDiff.gradient(x->TensorNetworkEvolve.iloss2(h, p1, _complex(x))[], reinterpret(Float64, variables(p1)))
    @test fvec.data ≈ _complex(fvec2) * -im
end


text:
using LinearAlgebra

mutable struct QLearning{T,TQ,OQ} <: AbstractAgent where {T<:Real,TQ<:AbstractPolicy,OQ<:AbstractOptimizer}
    qf::TQ
    qopt::OQ
    γ::T
    ϵ::T
    gp::Array{T, 2}
    action::Int

    function QLearning(qf::AbstractFuncApprox, qopt::AbstractOptimizer, γ::T, ϵ::T) where {T<:Real}
        tp = eltype(get_params(qf))
        sp = size(get_params(qf))
        qf2 = clone(qf)
        qopt2 = clone(qopt, qf2)
        new{tp,typeof(qf2),typeof(qopt2)}(qf2, qopt2, tp(γ), tp(ϵ), zeros(tp, sp), -1)
    end
end

function act!(agent::TA, env::TE, rng::AbstractRNG, learn::Bool) where {TA<:QLearning,TE<:AbstractEnvironment} # do not learn if learn flag is set
    qs = zeros(get_num_outputs(agent.qf))
    call!(qs, agent.qf, env.state)

    if rand(rng, 1)[1] < agent.ϵ
        agent.action = rand(rng, 1:length(qs))[1]
    else
        agent.action = argmax(qs)
    end

    r = step!(env, agent.action, rng)
end

function act!(agent::QLearning, env::AbstractEnvironment, rng::AbstractRNG)
    x = env.state
    fill!(agent.gp, 0.)
    q = 0.
    if agent.action != -1
        q = call_gradient!(agent.gp, agent.qf, env.state, agent.action)
    else
        qs = zeros(get_num_outputs(agent.qf))
        call_gradient!(qs, agent.gp, agent.qf, env.state)

        if rand(rng, 1)[1] < agent.ϵ
            agent.action = rand(rng, 1:length(qs))[1]
        else
            agent.action = argmax(qs)
        end
        q = qs[agent.action]
        agent.gp[:, 1:end .!= agent.action] .= 0.
    end

    r = step!(env, agent.action, rng)

    if is_terminal(env)
        target = r
    else
        qs = call(agent.qf, env.state)
        probs = softmax(qs)

        if rand(rng, 1)[1] < agent.ϵ
            agent.action = rand(rng, 1:length(qs))[1]
        else
            agent.action = argmax(qs)
        end

        #TD update
        q2 = maximum(qs)
        target = r + agent.γ * q2
    end

    delta = target - q

    update!(agent.qopt, delta, agent.gp)
    if is_terminal(env)
        new_episode!(agent, rng)
    end
end

function act!(agent::QLearning{T,TQ,TO}, env::AbstractEnvironment, rng::AbstractRNG) where {T,TQ,TO<:Union{Parl2AccumulatingTraceOptimizer{T,TQ},AutoTIDBDOptimizer}}
    x = env.state
    fill!(agent.gp, 0.)
    q = 0.
    if agent.action != -1
        q = call_gradient!(agent.gp, agent.qf, env.state, agent.action)
    else
        qs = zeros(get_num_outputs(agent.qf))
        call_gradient!(qs, agent.gp, agent.qf, env.state)

        if rand(rng, 1)[1] < agent.ϵ
            agent.action = rand(rng, 1:length(qs))[1]
        else
            agent.action = argmax(qs)
        end
        q = qs[agent.action]
        agent.gp[:, 1:end .!= agent.action] .= 0.
    end

    r = step!(env, agent.action, rng)
    gp2 = zeros(T, size(agent.gp))
    if is_terminal(env)
        target = r
    else
        qs = call(agent.qf, env.state)
        probs = softmax(qs)

        if rand(rng, 1)[1] < agent.ϵ
            agent.action = rand(rng, 1:length(qs))[1]
        else
            agent.action = argmax(qs)
        end
        gradient!(gp2, agent.qf, env.state, agent.action)
        #TD update
        q2 = maximum(qs)
        target = r + agent.γ * q2
    end

    delta = target - q

    update!(agent.qopt, delta, agent.gp, gp2)
    if is_terminal(env)
        new_episode!(agent, rng)
    end
end

function new_episode!(agent::QLearning, rng::AbstractRNG)
    new_episode!(agent.qopt)
    agent.action = -1
end

function clone(agent::QLearning)::QLearning
    a = QLearning(agent.qf, agent.qopt, agent.γ, agent.ϵ)
    return a
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
@testset "HyperRectangles" begin
    c = translate([4,5,6])Cuboid([1,2,3])
    @test HyperRectangle(c) ==
          HyperRectangle{3,Float64}(SVector{3}(4,5,6),SVector{3}(1,2,3))
    s = Sphere(3)
    @show s
end


text:
using Insolation
using Plots
using Dates
using Statistics
using Formatting

using CLIMAParameters: AbstractParameterSet
const APS = AbstractParameterSet

"""
    calc_day_lat_insolation(n_days::I,
                            n_lats::I,
                            param_set::APS) where {I<:Int}
"""
function calc_day_lat_insolation(n_days::I,
                                n_lats::I,
                                param_set::APS) where {I<:Int}
  d_arr = Array{I}(round.(collect(range(0, stop = 365, length = n_days))))
  l_arr = collect(range(-90, stop = 90, length = n_lats))
  F_arr = zeros(n_days, n_lats)
  # loop over days
  for (i, d) in enumerate(d_arr)
    for (j, lat) in enumerate(l_arr)
        date = Dates.DateTime(2000,1,1) + Dates.Day(d)
        θ, dist = daily_zenith_angle(date, lat, param_set, milankovitch=false)
        F_arr[i, j] = insolation(θ, dist, param_set)
    end
  end
  return d_arr, l_arr, F_arr
end

"""
    plot_day_lat_insolation(n_days::I,
                            n_lats::I,
                            F_arr::Array{FT},
                            scmap,
                            stitle,
                            file_name) where {FT<:AbstractFloat,I<:Int}
"""
function plot_day_lat_insolation(d_arr::Array{I},
                                 l_arr::Array{FT},
                                 F_arr::Array{FT},
                                 scmap,
                                 stitle,
                                 file_name) where {FT<:AbstractFloat,I<:Int}
  if scmap == "YlOrRd"
    cmap = :YlOrRd
    vmin, vmax = 0, ceil(max(F_arr...)/100)*100
  elseif scmap == "PRGn"
    cmap = :PRGn
    vmin, vmax = ceil(max(abs.(F_arr)...)/10)*-10, ceil(max(abs.(F_arr)...)/10)*10
  end

  p1 = contourf(d_arr, l_arr, F_arr', c=cmap, clims=(vmin,vmax),
    title=stitle, xlabel="Days since Jan 1 2000", ylabel="Latitude", colorbar_title="ToA Insolation [W/m2]")
  meanF = mean(F_arr, dims=1)[1,:]
  x1,x2 = [floor(minimum(meanF)), ceil(maximum(meanF))]
  p2 = plot(meanF, l_arr, xlabel="Annual-mean TOA \nInsolation [W/m2]",
    ylims=[-90,90], yticks=[-90,-60,-30,0,30,60,90], xlims=[x1,x2], legend = false)
  plot(p1, p2, size=(800,400), layout = grid(1,2, widths=(0.8,0.2)), dpi=250,
        left_margin=4Plots.mm, bottom_margin=6Plots.mm, right_margin=6Plots.mm)

  savefig(file_name)
end


text:
module TerminalExtensions

using Compat
#
# None of these functions are operating system specific because they
# might be connect via e.g. SSH to a different client operating system
#

import Base.Terminals: CSI

const DCS = "\eP"
const ST  = "\e\\"

function readDCS(io::IO)
    while nb_available(STDIN) >= 2
        c1 = read(io,UInt8)
        c1 == 0x90 && return true
        if c1 == UInt8('\e')
            read(io,UInt8) == UInt8('P') && return true
        end
    end
    return false
end

function readST(io::IO)
    c1 = read(io,UInt8)
    c1 == 0x90 && return true
    c1 != UInt8('\e') && return false
    read(io,UInt8) != UInt8('\\') && return false
    return true
end

module iTerm2

    import Base: display
    using Compat

    immutable InlineDisplay <: Display; end

    function set_mark()
        "\033]50;SetMark\007"
    end

    # Runs after interactively edited command but before execution
    function preexec()
        "\033]133;C\007"
    end

    function remotehost_and_currentdir()
        return string("\033]1337;RemoteHost=",ENV["USER"],"@",readstring(`hostname -f`),"\007","\033]1337;CurrentDir=",pwd(),"\007")
    end

    function prompt_prefix(last_success = true)
        return string("\033]133;D;$(convert(Int, last_success))\007",remotehost_and_currentdir(),"\033]133;A\007")
    end

    function prompt_suffix()
        return "\033]133;B\007"
    end

    function shell_version_number()
        return "\033]1337;ShellIntegrationVersion=1\007"
    end

    function prepare_display_file(;filename="Unnamed file", size=nothing, width=nothing, height=nothing, preserveAspectRation::Bool=true, inline::Bool=false)
        q = "\e]1337;File="
        options = ASCIIString[]
        filename != "Unnamed file" && push!(options,"name=" * base64encode(filename))
        size !== nothing && push!(options,"size=" * dec(size))
        height !== nothing && push!(options,"height=" * height)
        width !== nothing && push!(options,"width=" * width)
        preserveAspectRation !== true && push!(options,"preserveAspectRation=0")
        inline !== false && push!(options,"inline=1")
        q *= join(options,';')
        q *= ":"
        write(STDOUT,q)
    end

    function display_file(data::Vector{UInt8}; kwargs...)
        prepare_display_file(;kwargs...)
        write(STDOUT,base64encode(data))
        write(STDOUT,'\a')
    end

    # Incomplete list. Will be extended as necessity comes up
    const iterm2_mimes = ["image/png", "image/gif", "image/jpeg", "application/pdf", "application/eps"]

    for mime in iterm2_mimes
        @eval begin
            function display(d::InlineDisplay, m::MIME{Symbol($mime)}, x)
                prepare_display_file(;filename="image",inline=true)
                buf = IOBuffer()
                writemime(Base.Base64EncodePipe(buf),m,x)
                write(STDOUT, takebuf_array(buf))
                write(STDOUT,'\a')
            end
        end
    end

    function display(d::InlineDisplay,x)
        for m in iterm2_mimes
            if mimewritable(m,x)
                return display(d,m,x)
            end
        end
        throw(MethodError(display, (d,x)))
    end

end

function __init__()
    if !(isinteractive() && isdefined(Base, :active_repl))
        return
    end
    begin
        term = Base.Terminals.TTYTerminal("xterm",STDIN,STDOUT,STDERR)
        Base.Terminals.raw!(term,true)
        Base.start_reading(STDIN)

        # Detect iTerm support
        println("\e[1337n\e[5n")
        readuntil(STDIN, "\e")
        itermname = ""
        if read(STDIN, Char) != '0'
            itermname = readuntil(STDIN, "\e")[1:end-2]
        end
        # Read the rest of the \e[5n query
        read(STDIN, 3)

        if startswith(itermname, "ITERM2")
            pushdisplay(iTerm2.InlineDisplay())
            repl = Base.active_repl#REPL.LineEditREPL(Terminals.TTYTerminal("xterm",STDIN,STDOUT,STDERR))

            if !isdefined(repl,:interface)
                repl.interface = Base.REPL.setup_interface(repl)
            end

            let waserror = false
                prefix = repl.interface.modes[1].prompt_prefix
                repl.interface.modes[1].prompt_prefix = function ()
                    (TerminalExtensions.iTerm2.prompt_prefix(waserror) * (isa(prefix,Function) ? prefix() : prefix))
                end
                suffix = repl.interface.modes[1].prompt_suffix
                repl.interface.modes[1].prompt_suffix = function ()
                    ((isa(suffix,Function) ? suffix() : suffix) * TerminalExtensions.iTerm2.prompt_suffix())
                end
                for mode in repl.interface.modes
                    if isdefined(mode,:on_done)
                        of = mode.on_done
                        mode.on_done = function (args...)
                            print(STDOUT,TerminalExtensions.iTerm2.preexec())
                            of(args...)
                            waserror = repl.waserror
                        end
                    end
                end
            end
        end
    end
end

export queryTermcap, iTerm2

end # module


text:
using LinearAlgebra
using Random
using NNlib: softplus
using Roots # for inverse

################################################################################
#                            Planar and Radial Flows                           #
#             Ref: Variational Inference with Normalizing Flows,               #
#               D. Rezende, S. Mohamed(2015) arXiv:1505.05770                  #
################################################################################

###############
# PlanarLayer #
###############

# TODO: add docstring

mutable struct PlanarLayer{T1<:AbstractVector{<:Real}, T2<:Union{Real, AbstractVector{<:Real}}} <: Bijector{1}
    w::T1
    u::T1
    b::T2
end
function Base.:(==)(b1::PlanarLayer, b2::PlanarLayer)
    return b1.w == b2.w && b1.u == b2.u && b1.b == b2.b
end

function get_u_hat(u, w)
    # To preserve invertibility
    x = w' * u
    return u .+ (planar_flow_m(x) - x) .* w ./ sum(abs2, w)   # from A.1
end

function PlanarLayer(dims::Int, wrapper=identity)
    w = wrapper(randn(dims))
    u = wrapper(randn(dims))
    b = wrapper(randn(1))
    return PlanarLayer(w, u, b)
end

planar_flow_m(x) = -1 + softplus(x)   # for planar flow from A.1
ψ(z, w, b) = (1 .- tanh.(w' * z .+ b).^2) .* w    # for planar flow from eq(11)

# An internal version of transform that returns intermediate variables
function _transform(flow::PlanarLayer, z::AbstractVecOrMat)
    return _planar_transform(flow.u, flow.w, first(flow.b), z)
end
function _planar_transform(u, w, b, z)
    u_hat = get_u_hat(u, w)
    transformed = z .+ u_hat .* tanh.(w' * z .+ b) # from eq(10)
    return (transformed = transformed, u_hat = u_hat)
end

(b::PlanarLayer)(z) = _transform(b, z).transformed

function forward(flow::PlanarLayer, z::AbstractVecOrMat)
    transformed, u_hat = _transform(flow, z)
    # Compute log_det_jacobian
    psi = ψ(z, flow.w, first(flow.b)) .+ zero(eltype(u_hat))
    if psi isa AbstractVector
        T = eltype(psi)
    else
        T = typeof(vec(psi))
    end
    log_det_jacobian::T = log.(abs.(1.0 .+ psi' * u_hat)) # from eq(12)
    return (rv = transformed, logabsdetjac = log_det_jacobian)
end

function (ib::Inverse{<:PlanarLayer})(y::AbstractVector{<:Real})
    flow = ib.orig
    w = flow.w
    b = first(flow.b)
    u_hat = get_u_hat(flow.u, w)

    # Find the scalar ``alpha`` from A.1.
    wt_y = dot(w, y)
    wt_u_hat = dot(w, u_hat)
    alpha = find_alpha(y, wt_y, wt_u_hat, b)

    return y .- u_hat .* tanh(alpha * norm(w, 2) + b)
end

function (ib::Inverse{<:PlanarLayer})(y::AbstractMatrix{<:Real})
    flow = ib.orig
    w = flow.w
    b = first(flow.b)
    u_hat = get_u_hat(flow.u, flow.w)

    # Find the scalar ``alpha`` from A.1 for each column.
    wt_u_hat = dot(w, u_hat)
    alphas = mapvcat(eachcol(y)) do c
        find_alpha(c, dot(w, c), wt_u_hat, b)
    end

    return y .- u_hat .* tanh.(alphas' .* norm(w, 2) .+ b)
end

"""
    find_alpha(y::AbstractVector{<:Real}, wt_y, wt_u_hat, b)

Compute an (approximate) real-valued solution ``α`` to the equation
```math
wt_y = α + wt_u_hat tanh(α + b)
```

The uniqueness of the solution is guaranteed since ``wt_u_hat ≥ -1``.
For details see appendix A.1 of the reference.

# References

D. Rezende, S. Mohamed (2015): Variational Inference with Normalizing Flows.
arXiv:1505.05770
"""
function find_alpha(y::AbstractVector{<:Real}, wt_y, wt_u_hat, b)
    # Compute the initial bracket ((-Inf, 0) or (0, Inf))
    f0 = wt_u_hat * tanh(b) - wt_y
    zero_f0 = zero(f0)
    if f0 < zero_f0
        initial_bracket = (zero_f0, oftype(f0, Inf))
    else
        initial_bracket = (oftype(f0, -Inf), zero_f0)
    end
    alpha = find_zero(initial_bracket) do x
        x + wt_u_hat * tanh(x + b) - wt_y
    end

    return alpha
end

logabsdetjac(flow::PlanarLayer, x) = forward(flow, x).logabsdetjac
isclosedform(b::Inverse{<:PlanarLayer}) = false


text:
# Try to approximate squared-exponential kernel:
using PyPlot
using Optim
using Roots

function psd(abest)
nfreq = 1000
#omega = logspace(-1,1,nfreq)
omega = logspace(0,1,nfreq)
psd = zeros(nfreq)
# Compute the asymptotic value:
coeff_inf = 0.0
for k=1:nexp
  a1 = abest[1+(k-1)*4]
  a2 = abest[2+(k-1)*4]
  b  = abest[3+(k-1)*4]
  om = abest[4+(k-1)*4]
  psdt = (a2*om.*(b^2 + om^2 - omega.^2) + a1*b.*(b^2 + om^2 + omega.^2))./
         (b^4 + (om^2 - omega.^2).^2 + 2*b^2.*(om^2 + omega.^2))
  psd += psdt
end
return psd
end

#function gaussian_approx6()

ntime = 100
time = collect(linspace(0,10,ntime))
gaussian=exp(-time.^2/2.)

function gauss_model(x)
model = zeros(ntime)
for k=1:nexp
  model += x[1+(k-1)*4].*exp(-x[3+(k-1)*4].*time).*cos(x[4+(k-1)*4].*time)
  model += x[2+(k-1)*4].*exp(-x[3+(k-1)*4].*time).*sin(x[4+(k-1)*4].*time)
end
return model
end

function compute_chi(x)
  function dpsd_domega(omega)
# Let the kernel components be:
# a1*Exp[-b*t]*Cos[omega_0*t] + a2*Exp[-b*t]*Sin[omega_0*t]
# x = [{a1,a2,b,omega_0}_1; {a1,a2,b,omega_0}_2; {a1,a2,b,omega_0}_3; ... ]
#
  dpsd = zero(Real)
  for i=1:nexp
    a1 = x[(i-1)*4+1]
    a2 = x[(i-1)*4+2]
    b  = x[(i-1)*4+3]
    om0 = x[(i-1)*4+4]
    dpsd += (-2*sqrt(2/pi)*omega*(a1*b*(b^4 - 3*om0^4 + 2*om0^2*omega^2 + omega^4 - 2*b^2*(om0^2 - omega^2)) -
       a2*om0*(-3*b^4 + (om0^2 - omega^2)^2 - 2*b^2*(om0^2 + omega^2))))/
       ((b^2 + (om0 - omega)^2)^2*(b^2 + (om0 + omega)^2)^2)
  end
  return dpsd
  end

  function psd_func(omega)
  psd = 0.
  for k=1:nexp
    a1 = x[1+(k-1)*4]
    a2 = x[2+(k-1)*4]
    b  = x[3+(k-1)*4]
    om = x[4+(k-1)*4]
    psdt = (a2*om*(b^2 + om^2 - omega^2) + a1*b*(b^2 + om^2 + omega^2))/
         (b^4 + (om^2 - omega^2)^2 + 2*b^2*(om^2 + omega^2))
    psd += psdt
  end
  return psd
  end

model = gauss_model(x)
sum_coeff = zeros(x)
# Compute Taylor-series coefficients:
for k=1:nexp
  sum_coeff[1] += x[1+(k-1)*4]
  sum_coeff[2] +=-x[1+(k-1)*4]*x[3+(k-1)*4]+x[2+(k-1)*4]*x[4+(k-1)*4]
  sum_coeff[3] += x[1+(k-1)*4]*x[3+(k-1)*4]^2-2*x[2+(k-1)*4]*x[3+(k-1)*4]*x[4+(k-1)*4]-x[1+(k-1)*4]*x[4+(k-1)*4]^2
  sum_coeff[4] +=-x[1+(k-1)*4]*x[3+(k-1)*4]^3+3*x[2+(k-1)*4]*x[3+(k-1)*4]^2*x[4+(k-1)*4]+3*x[1+(k-1)*4]*x[3+(k-1)*4]*x[4+(k-1)*4]^2-x[2+(k-1)*4]*x[4+(k-1)*4]^3
end

#minpsd = minimum(psd(x))
#omega_min = fzero(dpsd_domega,[1e-3,1e3])
#minpsd = psd_func(omega_min)
result = optimize(psd_func,1.0,10.0,Brent())
minpsd = psd_func(result.minimum)
println("PSD minimum: ",minpsd," ",result.minimum)

return sum((gaussian-model).^2) +100.*((sum_coeff[1]-1.0)^2+sum_coeff[2]^2+(sum_coeff[3]+1.0)^2+sum_coeff[4]^2)+(minpsd < 0)*abs(minpsd)*1e5
end

# The number of exponentials is actually +1:
nexp = 4
npar = 4*nexp
#abest=zeros(npar)

#atrial = zeros(npar)
#for k=1:npar
#  atrial[k] = abest[k]
#end
#atrial[1:8]=0.5
#atrial = [-0.2967663878261003,0.25057478717710174,1.3827297351466428,1.702023384832311,-0.2967663878261003,0.25057478717710174,1.3827297351466428,1.702023384832311,3.0791800790217154,0.0,1.2400386725361678,0.0,-0.3561416933055925,0.0,2.801219489380058,0.0,-1.1290392201099397,0.0,1.0162289199237395,0.0]
#atrial = [-0.2967663878261003*2,0.25057478717710174*2,1.3827297351466428,1.702023384832311,3.0791800790217154,0.0,1.2400386725361678,0.0,-0.3561416933055925,0.0,2.801219489380058,0.0,-1.1290392201099397,0.0,1.0162289199237395,0.0]
#atrial = [-0.5489307813321151,0.3318064860491032,1.3836581492240345,1.7999142202401186,26.778627442628697,0.0,1.2786729767481086,0.0,-1.8044641240364967,0.0,2.018989606638419,0.0,-23.42487315195119,0.0,1.248306537682799,0.0] #0.0003011873652880382
#atrial = [-0.5336021091255497,0.2999401472397189,1.3833756425260124,1.8206931811172136,16.30785328812252,0.0,1.3870095791283663,0.0,-4.448314627596262,0.0,1.7985112282392832,0.0,-10.3255967541902,0.0,1.29143086997431,0.0]
#atrial = [-0.5297041429112943,0.2933240470063606,1.3830031311925879,1.8253942281094937,14.986091075429751,0.0,1.464640372400348,0.0,-7.931800351569044,0.0,1.717799667957279,0.0,-5.524251375777077,0.0,1.2773197102081726,0.0]
#atrial = [-0.5243881693053298,0.2836999647812343,1.3832697630116844,1.8321298873666765,49.904780921019984,0.0,1.4654597590489415,0.0,-22.280787592453734,0.0,1.6135419037436551,0.0,-26.099271853203692,0.0,1.3769569862282118,0.0]
#atrial = [-0.5251847078236657,0.2835380989822935,1.3837836249167428,1.8318749094338078,53.123080563224285,0.0,1.4718814839416798,0.0,-25.331145427490217,0.0,1.6066377692092113,0.0,-26.26641632953569,0.0,1.3799755188402958,0.0]
#atrial = [-0.5249565598434786,0.28333358471586606,1.383715102743298,1.8320654762189466,55.52393985773224,0.0,1.4749129722610685,0.0,-27.470786754029175,0.0,1.601868924110799,0.0,-26.523862614593764,0.0,1.3812751618641164,0.0]
#New minimum: 0.00023883842674643924 [-0.5247762751816804,0.2828915570634826,1.3837656415003445,1.8323467099030555,63.120917132667905,0.0,1.478130410600365,0.0,-32.17183668366479,0.0,1.5929285243364661,0.0,-29.423970266940803,0.0,1.38693874263686,0.0] 0.00047359417598702186 0.0012219351046597327
#atrial = [-0.5247762751816804,0.2828915570634826,1.3837656415003445,1.8323467099030555,63.120917132667905,0.0,1.478130410600365,0.0,-32.17183668366479,0.0,1.5929285243364661,0.0,-29.423970266940803,0.0,1.38693874263686,0.0] # 0.00047359417598702186 0.0012219351046597327
#New minimum: 0.00023883842674643924 [-0.5247762751816804,0.2828915570634826,1.3837656415003445,1.8323467099030555,63.120917132667905,0.0,1.478130410600365,0.0,-32.17183668366479,0.0,1.5929285243364661,0.0,-29.423970266940803,0.0,1.38693874263686,0.0] 0.00047359417598702186 0.0012219351046597327
#atrial = [-0.5237762751816804,0.2828915570634826,1.3837656415003445,1.8323467099030555,63.120917132667905,0.0,1.478130410600365,0.0,-32.17183668366479,0.0,1.5929285243364661,0.0,-29.423970266940803,0.0,1.38693874263686,0.0]

atrial = [-0.5185744409672999,0.2854238724264553,1.3848025999389613,1.8339730386273703,63.11883856322811,-0.002585613714221086,1.479593361171615,-0.001569459359505545,-32.174013349165506,-0.00439276758185438,1.593669981824972,0.00035151389750538633,-29.42595216593894,-0.0013496416338785663,1.3889097102817822,0.002608575323797371]

model = gauss_model(atrial)
clf()
plot(time,gaussian)
plot(time,model)
#chibest=sum((gaussian-model).^2)
chibest=compute_chi(atrial)

println("Initial chi-square: ",chibest)
read(STDIN,Char)

#result = optimize(compute_chi, atrial, BFGS(), OptimizationOptions(autodiff = true))
result = optimize(compute_chi, atrial)

println(Optim.minimizer(result),Optim.minimum(result))
abest= Optim.minimizer(result)
clf()
model = gauss_model(abest)
plot(time,gaussian)
plot(time,model)
plot(time,(gaussian-model))
chi=sum((gaussian-model).^2)
println("New minimum: ",chi," ",abest," ",std(gaussian-model)," ",maximum(abs(gaussian-model)))
read(STDIN,Char)
# Now, plot the power spectrum:
nfreq = 100000
#omega = logspace(-1,1,nfreq)
omega = logspace(-3,1,nfreq)
psd_grid = zeros(nfreq)
clf()
#plot(omega,psd_grid)
# Compute the asymptotic value:
coeff_inf = 0.0
for k=1:nexp
  a1 = abest[1+(k-1)*4]
  a2 = abest[2+(k-1)*4]
  b  = abest[3+(k-1)*4]
  om = abest[4+(k-1)*4]
  psdt = (a2*om.*(b^2 + om^2 - omega.^2) + a1*b.*(b^2 + om^2 + omega.^2))./
         (b^4 + (om^2 - omega.^2).^2 + 2*b^2.*(om^2 + omega.^2))
  psd_grid += psdt
  coeff_inf += a1*b-a2*om
#  plot(omega,psdt)
#  plot(omega,(a1*b-a2*om)./omega.^2,"r.")
#  read(STDIN,Char)
end
clf()
plot(omega,psd_grid)
plot(omega,exp(-omega.^2./2)*sqrt(pi/2))
#plot(omega,coeff_inf./omega.^2)
#return
#end

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess


text:
"""************************ initialization schemes **************************"""

"""
    finiteMPS(N, d, maxD)

constructs a random finite MPS for a  system with N sites with local physical dimsension
d and max. bond dimension maxD. Bond dimension will increas dynamically with
system size but wont exceed maxD.

# Arguments
- N: System size
- d: local physical dimension
- maxD: maximal bond Dimension

return:
	- MPSvec vector storing M for each site
	- Cvec vector storing singular Values between site and site+1
	  so the last entrie of Cvec is basically the norm of the sate
 	  and of dimension 1 x 1

"""
function randMPS(N::Int, d::Int, maxD::Int)
	MPSvec = []
	Cvec = []

	N_half = round(N/2, RoundDown)
	for i = 1:N_half

		α_dim = d^(i-1) > maxD ? maxD : Int(d^(i-1))
		β_dim = d^i > maxD ? maxD : Int(d^i)
		push!(MPSvec, randn(ComplexF64, (α_dim, d, β_dim)))
		push!(Cvec, randn(ComplexF64, (β_dim, β_dim)))
	end

	if N-2*N_half != 0
		bDim = d^N_half > maxD ? maxD : Int(d^N_half)
		push!(MPSvec, randn(ComplexF64, (bDim, d, bDim)))
		push!(Cvec, randn(ComplexF64, (bDim, bDim)))
	end

	for i = N_half:-1:1

		β_dim = d^(i-1) > maxD ? maxD : Int(d^(i-1))
		α_dim = d^i > maxD ? maxD : Int(d^i)

		push!(MPSvec, randn(ComplexF64, (α_dim, d, β_dim)))
		push!(Cvec, randn(ComplexF64, (β_dim, β_dim)))

	end

	return MPSvec, Cvec
end

function randomBloch(N::Int, maxD::Int)
		MPSvec = []
		C = []

		N_half = round(N/2, RoundDown)

		for i = 1:N_half

			ϕ = rand(Uniform(0.0,2.0*pi))
			θ = rand(Uniform(0.0,2.0*pi))
			α_dim = 2^(i-1) > maxD ? maxD : Int(2^(i-1))
			β_dim = 2^i > maxD ? maxD : Int(2^i)

			A = zeros(ComplexF64,α_dim, 2, β_dim)
			c = zeros(Float64, β_dim, β_dim)
			A[1,1,1] = cos(θ/2.0)
			A[1,2,1] = sin(θ/2.0)*exp(1.0im*ϕ)
			c[1,1] = 1.0

			push!(MPSvec, A)
			push!(C, c)
		end

		if N-2*N_half != 0
			bDim = 2^N_half > maxD ? maxD : Int(2^N_half)

			ϕ = rand(Uniform(0.0,2.0*pi))
			θ = rand(Uniform(0.0,2.0*pi))

			A = zeros(ComplexF64, bDim, 2, bDim)
			c = zeros(Float64, bDim, bDim)
			A[1,1,1] = cos(θ/2.0)
			A[1,2,1] = sin(θ/2.0)*exp(1.0im*ϕ)
			c[1,1] = 1.0

			push!(MPSvec, A)
			push!(C, c)
		end

		for i = N_half:-1:1

			β_dim = 2^(i-1) > maxD ? maxD : Int(2^(i-1))
			α_dim = 2^i > maxD ? maxD : Int(2^i)

			ϕ = rand(Uniform(0.0,2.0*pi))
			θ = rand(Uniform(0.0,2.0*pi))

			A = zeros(ComplexF64,α_dim, 2, β_dim)
			c = zeros(Float64, β_dim, β_dim)

			push!(MPSvec, A)
			push!(C, c)

		end

		return MPSvec, C
 end

# this should be essential turned into a convertion once MPS struct is defined
function Coef_to_MPS(N::Int, Coef::Vector{<:Number})
	MPSvec = []

	tmp = zeros(ComplexF64, 1,size(Coef)[1])
	tmp[1,:] = Coef
	for site in 1:N
		bout_pre, d_left = size(tmp)

		tmp =  reshape(tmp, bout_pre *2, Int(d_left/2))

		F = svd(tmp)

		mps = reshape(F.U, bout_pre, 2, size(F.U)[2])

		tmp = diagm(0=>F.S)*F.V'

		push!(MPSvec, mps)

	end

	return MPSvec
end

max_stars_repo_name: Gawatz/finiteMPS

max_issues_repo_name: Gawatz/finiteMPS

max_forks_repo_name: Gawatz/finiteMPS


text:
function SolveModel(A, obj, b, n)
# Solves the forward optimization problem

# Step 1: Create the model
    m1 = Model(solver=GurobiSolver(OutputFlag = 0))
    @variable(m1, x[1:n])
    @objective(m1, Max, dot(obj, x))
    for i = 1:length(A[:, 1])
        @constraint(m1, dot(A[i, :], x) <= b[i])
    end

# Step 2: Solve the model
    status = solve(m1; suppress_warnings = true)
# Step 3: return optimal value;
    val = getvalue(x);
end

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

function SolveFP(cone, iteration)
    # Feasibility Problem to be solved at every iteration of the decomposition algorithm
    epsilon = 1
    m2 = Model(solver=GurobiSolver(OutputFlag = 1))
    @variable(m2, y[1:n]>=1e-4)
    @variable(m2, gamma[1:iteration, 1:c] >= 0)
    @variable(m2, w, Bin)

    for i = 1:iteration
        for j = 1:n
            @constraint(m2, y[j] == sum(gamma[i,k]*cone[i][k,j] for k = 1:length(cone[i][:,1])))
        end
    end

    @constraint(m2, sum(y) == epsilon)

    @objective(m2, Min, 0)
    status = solve(m2; suppress_warnings = true)
    return status
end

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

function SolveP2(C_bar, cone, iteration)

    m3 = Model(solver=GurobiSolver(OutputFlag = 0))
    @variable(m3, y[1:n])
    @variable(m3, gamma[1:iteration, 1:c] >= 0)

    for i = 1:iteration
        for j = 1:n
            @constraint(m3, y[j] == sum(gamma[i,k]*cone[i][k,j] for k = 1:length(cone[i][:,1])))
        end
    end

    @objective(m3, Min, sum((y[j]-C_bar[j])^2 for j = 1:n))
    status = solve(m3; suppress_warnings = true)
    return getvalue(y), getobjectivevalue(m3)
end

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

function SolveP1(A, b, x, NoS, J)
    # NoS (number of samples, i.e., I is always 1 here.)
    n = length(A[1, :, 1])
    m4 = Model(solver = GurobiSolver(OutputFlag = 0, IntFeasTol = 1e-9, MIPGapAbs = 1e-5, MIPGap = 0))
    M = 50
    epsilon = 1

    @variable(m4, xbar[1:NoS, 1:n], start = 0)
    @variable(m4, lambda[1:NoS, 1:c] >= 0)
    @variable(m4, s[1:NoS, 1:c] >= 0)
    @variable(m4, z[1:NoS, 1:c], Bin)
    @variable(m4, C[1:n]>=1e-4)
    @variable(m4, t[1:J, 1:n, 1:NoS] >= 0)

    @constraint(m4, S[k = 1:NoS, j = 1:n], C[j] - sum(lambda[k, i]*A[i, j, k] for i = 1:c) == 0)
    @constraint(m4, PF[k = 1:NoS, i = 1:c], sum(A[i, j, k]*xbar[k, j] for j = 1:n) + s[k, i] == b[i])
    @constraint(m4, CS1[k = 1:NoS, i = 1:c], s[k, i] <= M*(1-z[k, i]))
    @constraint(m4, CS2[k = 1:NoS, i = 1:c], lambda[k, i] <= M*z[k, i])
    @constraint(m4, VC[k = 1:NoS] , sum(z[k,:]) >= n)
    @constraint(m4, linearization1[j = 1:J, i = 1:n, k = 1:NoS], t[j,i,k] >= xbar[k,i] - x[j,i,k])
    @constraint(m4, linearization2[j = 1:J, i = 1:n, k = 1:NoS], t[j,i,k] >= -(xbar[k,i] - x[j,i,k]))
    @constraint(m4, sum(C) == epsilon)

    @objective(m4, Min, sum(t))
    status = solve(m4; suppress_warnings = true)
    loss1 = getobjectivevalue(m4)

    return getvalue(xbar), loss1, getvalue(z)
end

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

function SolveP1m(A, b, x, NoS, J, instance, initx, initz)
    # Solves P1 with multiple (NoS > 1) experiments
    n = length(A[1, :, 1])
    initz = round.(initz)
    m5 = Model(solver = GurobiSolver(OutputFlag = 1, TimeLimit = 7200, ConcurrentMIP = 3,
    IntFeasTol = 1e-9, MIPGapAbs = 1e-5, MIPGap = 0))
    M = 20
    epsilon = 1

    @variable(m5, xbar[j = 1:NoS, k = 1:n])
    @variable(m5, lambda[1:NoS, 1:c] >= 0)
    @variable(m5, s[1:NoS, 1:c] >= 0)
    @variable(m5, z[j = 1:NoS, k = 1:c], Bin)
    for j = 1:NoS-1
        for k = 1:c
            setvalue(z[j,k], initz[j,k])
        end
    end
    @variable(m5, C[1:n]>=1e-4)
    @variable(m5, t[1:J, 1:n, 1:NoS] >= 0)

    @constraint(m5, S[k = 1:NoS, j = 1:n], C[j] - sum(lambda[k, i]*A[i, j, k] for i = 1:c) == 0)
    @constraint(m5, PF[k = 1:NoS, i = 1:c], sum(A[i, j, k]*xbar[k, j] for j = 1:n) + s[k, i] == b[k, i])
    @constraint(m5, CS1[k = 1:NoS, i = 1:c], s[k, i] <= M*(1-z[k, i]))
    @constraint(m5, CS2[k = 1:NoS, i = 1:c], lambda[k, i] <= M*z[k, i])
    @constraint(m5, VC[k = 1:NoS] , sum(z[k,:]) >= n)
    @constraint(m5, epsilon == sum(C))
    @constraint(m5, linearization1[j = 1:J, i = 1:n, k = 1:NoS], t[j,i,k] >= xbar[k,i] - x[j,i,k])
    @constraint(m5, linearization2[j = 1:J, i = 1:n, k = 1:NoS], t[j,i,k] >= -(xbar[k,i] - x[j,i,k]))

    @objective(m5, Min, sum(t))
    status = solve(m5; suppress_warnings = true)
    loss = getobjectivevalue(m5)
    return getvalue(xbar), loss, getvalue(z), status
end

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

function accuracy(x, x_hat)
NoS = length(x[:, 1])
n = length(x[1,:])
error = 0
    for i = 1:NoS
        if !isapprox(norm(x[i,:] - x_hat[i,:]), 0, atol = 1e-3)
            error += 1
        end
    end
    return error/NoS
end

max_stars_repo_name: ddolab/InvLP

max_issues_repo_name: ddolab/InvLP

max_forks_repo_name: ddolab/InvLP