text:
using UUIDs

@testset "AutoCommitTransactionTestCase" begin
    tc = setup()

    @testset "test_can_run_simple_statement" begin
        sess = session(tc.driver)
        result = run(sess, "RETURN 1 AS n")
        for record in result
            @test value(record, 1) == 1
            @test value(record, "n") == 1
            @test record["n"] == 1
            @test_throws KeyError record["x"]
            @test_throws KeyError record[(1,2)]
            @test length(record) == 1
        end
        close(sess)
    end

    @testset "test_can_run_simple_statement_with_params" begin
        sess = session(tc.driver)
        count = 0
        result = run(sess, "RETURN \$x AS n", Dict("x"=>Dict("abc"=>["d","e","f"])))
        for record in result
            @test value(record, 1) == Dict("abc"=>["d","e","f"])
            @test record["n"] == Dict("abc"=>["d","e","f"])
            @test length(record) == 1
            count+=1
        end
        close(sess)
        @test count == 1
    end

    @testset "test_autocommit_transactions_use_bookmarks" begin
        bookmarks = []

        sess = session(tc.driver)
        result = run(sess, "CREATE ()")
        consume(result)
        lb = last_bookmark(sess)
        @test lb != nothing
        push!(bookmarks, lb)
        close(sess)

        sess = session(tc.driver, bookmarks=bookmarks)
        @test sess.bookmarks_in == bookmarks
        result = run(sess, "CREATE ()")
        consume(result)
        lb = last_bookmark(sess)
        @test lb != nothing
        @test !in(lb, bookmarks)
        close(sess)
    end

    @testset "test_fails_on_bad_syntax" begin
        sess = session(tc.driver)
        @test_throws Neo4jBolt.JuliaBolt.CypherSyntaxError run(sess, "X")
    end

    @testset "test_fails_on_missing_parameter" begin
        sess = session(tc.driver)
        @test_throws Neo4jBolt.JuliaBolt.ClientError run(sess, "RETURN \$x")
    end

    @testset "test_can_run_simple_statement_from_bytes_string" begin
        sess = session(tc.driver)
        count = 0
        result = run(sess, b"RETURN 1 AS n")
        for record in result
            @test value(record, 1) == 1
            @test record["n"] == 1
            @test length(record) == 1
            count+=1
        end
        close(sess)
        @test count == 1
    end

    @testset "test_can_run_statement_that_returns_multiple_records" begin
        sess = session(tc.driver)
        count = 0
        result = run(sess, "unwind(range(1, 10)) AS z RETURN z")
        for record in result
            @test 1 <= value(record, 1) <= 10
            count+=1
        end
        close(sess)
        @test count == 10
    end

    @testset "test_can_use_with_to_auto_close_session" begin
        session(tc.driver) do sess
            record_list = collect(run(sess, "RETURN 1"))
            @test length(record_list) == 1
            for record in record_list
                @test value(record, 1) == 1
            end
        end
    end

    @testset "test_can_return_node" begin
        session(tc.driver) do sess
            record_list = collect(run(sess, "CREATE (a:Person {name:'Alice'}) RETURN a"))
            @test length(record_list) == 1
            for record in record_list
                alice = value(record, 1)
                @test isa(alice, Node)
                @test alice.labels == Set(["Person"])
                @test alice.properties == Dict("name"=>"Alice")
            end
        end
    end

    @testset "test_can_return_relationship" begin
        session(tc.driver) do sess
            record_list = collect(run(sess, "CREATE ()-[r:KNOWS {since:1999}]->() RETURN r"))
            @test length(record_list) == 1
            for record in record_list
                rel = value(record, 1)
                @test isa(rel, Relationship)
                @test rel.type == "KNOWS"
                @test rel.properties == Dict("since"=>1999)
            end
        end
    end

    @testset "test_can_return_path" begin
        session(tc.driver) do sess
            record_list = collect(run(sess, "MERGE p=({name:'Alice'})-[:KNOWS]->({name:'Bob'}) RETURN p"))
            @test length(record_list) == 1
            for record in record_list
                path = value(record, 1)
                @test isa(path, Path)
                @test start_node(path).properties["name"] == "Alice"
                @test end_node(path).properties["name"] == "Bob"
                @test path.relationships[1].type == "KNOWS"
                @test length(path.nodes) == 2
                @test length(path.relationships) == 1
            end
        end
    end

    @testset "test_can_handle_cypher_error" begin
        session(tc.driver) do sess
            @test_throws Neo4jBolt.JuliaBolt.CypherSyntaxError run(sess, "X")
        end
    end

    @testset "test_keys_are_available_before_and_after_stream" begin
        session(tc.driver) do sess
            result = run(sess, "UNWIND range(1, 10) AS n RETURN n")
            @test Neo4jBolt.keys(result) == ["n"]
            collect(result)
            @test Neo4jBolt.keys(result) == ["n"]
        end
    end

    @testset "test_keys_with_an_error" begin
        session(tc.driver) do sess
            @test_throws Neo4jBolt.JuliaBolt.CypherSyntaxError begin
                run(sess, "X")
                Neo4jBolt.keys(result)
            end
        end
    end

    @testset "test_should_not_allow_empty_statements" begin
        session(tc.driver) do sess
            @test_throws ErrorException run(sess, "")
        end
    end

    @testset "test_statement_object" begin
        session(tc.driver) do sess
            v = value(single(run(sess, "RETURN \$x", x=1)))
            @test v == 1
        end
    end

    @testset "test_statement_object" begin
        session(tc.driver) do sess
            v = value(single(run(sess, "RETURN \$x", x=1)))
            @test v == 1
        end
    end

    @testset "test_autocommit_transactions_should_support_timeout" begin
        session(tc.driver) do sess1
            consume(run(sess1, "CREATE (a:Node)"))
            session(tc.driver) do sess2
                tx1 = begin_transaction(sess1)
                consume(run(tx1, "MATCH (a:Node) SET a.property = 1"))
                @test_throws Neo4jBolt.JuliaBolt.TransientError consume(run(sess2, Statement("MATCH (a:Node) SET a.property = 2", timeout=0.25)))
            end
        end
    end

    teardown(tc.driver)
end

@testset "SummaryTestCase" begin
    tc = setup()

    @testset "test_can_obtain_summary_after_consuming_result" begin
        session(tc.driver) do sess
            result = run(sess, "CREATE (n) RETURN n")
            summ = Neo4jBolt.summary(result)
            @test summ.statement == "CREATE (n) RETURN n"
            @test summ.parameters == Dict()
            @test summ.statement_type == "rw"
            @test summ.counters.nodes_created == 1
        end
    end

    @testset "test_no_plan_info" begin
        session(tc.driver) do sess
            result = run(sess, "CREATE (n) RETURN n")
            summ = Neo4jBolt.summary(result)
            @test summ.plan == nothing
            @test summ.profile == nothing
        end
    end

    @testset "test_can_obtain_plan_info" begin
        session(tc.driver) do sess
            result = run(sess, "EXPLAIN CREATE (n) RETURN n")
            summ = Neo4jBolt.summary(result)
            plan = summ.plan
            @test plan.operator_type == "ProduceResults"
            @test plan.identifiers == ["n"]
            @test length(plan.children) == 1
        end
    end

    @testset "test_can_obtain_profile_info" begin
        session(tc.driver) do sess
            result = run(sess, "PROFILE CREATE (n) RETURN n")
            summ = Neo4jBolt.summary(result)
            profile = summ.profile
            @test profile.db_hits == 0
            @test profile.rows == 1
            @test profile.operator_type == "ProduceResults"
            @test profile.identifiers == ["n"]
            @test length(profile.children) == 1
        end
    end

    @testset "test_no_notification_info" begin
        session(tc.driver) do sess
            result = run(sess, "CREATE (n) RETURN n")
            summ = Neo4jBolt.summary(result)
            notifications = summ.notifications
            @test notifications == []
        end
    end

    @testset "test_can_obtain_notification_info" begin
        session(tc.driver) do sess
            result = run(sess, "EXPLAIN MATCH (n), (m) RETURN n, m")
            summ = Neo4jBolt.summary(result)
            notifications = summ.notifications

            @test length(notifications) == 1
            notification = notifications[1]
            @test notification.code == "Neo.ClientNotification.Statement.CartesianProductWarning"
            @test notification.title == "This query builds a cartesian product between " *
                                         "disconnected patterns."
            @test notification.severity == "WARNING"
            @test notification.description == "If a part of a query contains multiple " *
                                               "disconnected patterns, this will build a " *
                                               "cartesian product between all those parts. This " *
                                               "may produce a large amount of data and slow down " *
                                               "query processing. While occasionally intended, " *
                                               "it may often be possible to reformulate the " *
                                               "query that avoids the use of this cross product, " *
                                               "perhaps by adding a relationship between the " *
                                               "different parts or by using OPTIONAL MATCH " *
                                               "(identifier is: (m))"
            position = notification.position
            @test position != nothing
        end
    end

    @testset "test_contains_time_information" begin
        session(tc.driver) do sess
            summary = consume(run(sess, "UNWIND range(1,1000) AS n RETURN n AS number"))
            @test isa(summary.t_first, Integer)
            @test isa(summary.t_last, Integer)
        end
    end

    teardown(tc.driver)
end

@testset "ResetTestCase" begin
    tc = setup()

    @testset "test_automatic_reset_after_failure" begin
        session(tc.driver) do sess
            try
                consume(run(sess, "X"))
                throw(ErrorException("A Cypher error should have occurred"))
            catch ex
                if isa(ex, Neo4jBolt.JuliaBolt.CypherSyntaxError)
                    result = run(sess, "RETURN 1")
                    record = [r for r in result][1]
                    @test value(record, 1) == 1
                else
                    @test false
                end
            end
        end
    end

    teardown(tc.driver)
end

@testset "ExplicitTransactionTestCase" begin
    tc = setup()

    @testset "test_can_commit_transaction" begin
        session(tc.driver) do sess
            tx = begin_transaction(sess)
            result = run(tx, "CREATE (a) RETURN id(a)")
            record = [r for r in result][1]
            node_id = value(record, 1)
            @test isa(node_id, Integer)

            run(tx, "MATCH (a) WHERE id(a) = \$n SET a.foo = \$foo", Dict("n"=>node_id, "foo"=>"bar"))

            commit(tx)

            result = run(sess, "MATCH (a) WHERE id(a) = \$n RETURN a.foo", Dict("n"=>node_id))
            record = [r for r in result][1]
            v = value(record, 1)

            @test v == "bar"
        end
    end

    @testset "test_can_rollback_transaction" begin
        session(tc.driver) do sess
            tx = begin_transaction(sess)

            result = run(tx, "CREATE (a) RETURN id(a)")
            record = [r for r in result][1]
            node_id = value(record, 1)
            @test isa(node_id, Integer)

            run(tx, "MATCH (a) WHERE id(a) = \$n SET a.foo = \$foo", Dict("n"=>node_id, "foo"=>"bar"))

            rollback(tx)

            result = run(sess, "MATCH (a) WHERE id(a) = \$n RETURN a.foo", Dict("n"=>node_id))
            @test length([r for r in result]) == 0
        end
    end

    @testset "test_can_rollback_transaction_using_with_block" begin
        session(tc.driver) do sess
            node_id = 0
            begin_transaction(sess) do tx
                result = run(tx, "CREATE (a) RETURN id(a)")
                record = [r for r in result][1]
                node_id = value(record, 1)
                @test isa(node_id, Integer)

                run(tx, "MATCH (a) WHERE id(a) = \$n SET a.foo = \$foo", Dict("n"=>node_id, "foo"=>"bar"))

                tx.success = false
            end
            result = run(sess, "MATCH (a) WHERE id(a) = \$n RETURN a.foo", Dict("n"=>node_id))
            @test length([r for r in result]) == 0
        end
    end

    @testset "test_broken_transaction_should_not_break_session" begin
        session(tc.driver) do sess
            @test_throws Neo4jBolt.JuliaBolt.CypherSyntaxError begin
                begin_transaction(sess) do tx
                    run(tx, "X")
                end
            end
            begin_transaction(sess) do tx
                 run(tx, "RETURN 1")
            end
        end
    end

    # @testset "test_last_run_statement_should_be_cleared_on_failure" begin
    #     session(tc.driver) do sess
    #         tx = begin_transaction(sess)
    #         consume(run(tx, "RETURN 1"))
    #         connection_1 = sess.connection
    #         @test connection_1.last_run_statement == "RETURN 1"
    #         result = run(tx, "X")
    #         connection_2 = sess.connection
    #         @test_throws Neo4jBolt.JuliaBolt.CypherSyntaxError consume(result)
    #         @test connection_2 === connection_1
    #         @test connection_2.last_run_statement == nothing
    #         close(tx)
    #     end
    # end

    @testset "test_statement_object_not_supported" begin
        session(tc.driver) do sess
            begin_transaction(sess) do tx
                 @test_throws ErrorException run(tx, Statement("RETURN 1", timeout=0.25))
            end
        end
    end

    @testset "test_transaction_timeout" begin
        session(tc.driver) do sess1
            consume(run(sess1, "CREATE (a:Node)"))
            session(tc.driver) do sess2
                tx1 = begin_transaction(sess1)
                consume(run(tx1, "MATCH (a:Node) SET a.property = 1"))
                tx2 = begin_transaction(sess2, timeout=0.25)
                @test_throws Neo4jBolt.JuliaBolt.TransientError consume(run(tx2, "MATCH (a:Node) SET a.property = 2"))
            end
        end
    end

    @testset "test_exit_after_explicit_close_should_be_silent" begin
        session(tc.driver) do sess
            tx_ = nothing
            begin_transaction(sess) do tx
                tx_ = tx
                @test !tx.closed
                close(tx)
                @test tx.closed
            end
            @test tx_.closed
        end
    end

    teardown(tc.driver)
end

@testset "BookmarkingTestCase" begin
    tc = setup()

    @testset "test_can_obtain_bookmark_after_commit" begin
        session(tc.driver) do sess
            begin_transaction(sess) do tx
                run(tx, "RETURN 1")
            end
            @test last_bookmark(sess) != nothing
        end
    end

    @testset "test_can_pass_bookmark_into_next_transaction" begin
        unique_id = string(uuid4())
        bookmark = nothing

        session(tc.driver, access_mode=WRITE_ACCESS) do sess
            begin_transaction(sess) do tx
                run(tx, "CREATE (a:Thing {uuid:\$uuid})", uuid=unique_id)
            end
            bookmark = last_bookmark(sess)
        end

        @test bookmark != nothing

        session(tc.driver, access_mode=READ_ACCESS, bookmark=bookmark) do sess
            begin_transaction(sess) do tx
                result = run(tx, "MATCH (a:Thing {uuid:\$uuid}) RETURN a", uuid=unique_id)
                record_list = [r for r in result]
                @test length(record_list) == 1
                record = record_list[1]
                @test length(record) == 1
                thing = value(record, 1)
                @test isa(thing, Node)
                @test thing.properties["uuid"] == unique_id
            end
        end
    end

    @testset "test_bookmark_should_be_none_after_rollback" begin
        sess_ = nothing

        session(tc.driver, access_mode=WRITE_ACCESS) do sess
            sess_ = sess
            begin_transaction(sess) do tx
                run(tx, "CREATE (a)")
            end
        end

        @test last_bookmark(sess_) != nothing

        session(tc.driver, access_mode=WRITE_ACCESS) do sess
            sess_ = sess
            begin_transaction(sess) do tx
                run(tx, "CREATE (a)")
                tx.success = false
            end
        end

        @test last_bookmark(sess_) == nothing
    end

    teardown(tc.driver)
end

@testset "SessionCompletionTestCase" begin
    tc = setup()

    @testset "test_should_sync_after_commit" begin
        session(tc.driver) do sess
            tx = begin_transaction(sess)
            result = run(tx, "RETURN 1")
            commit(tx)
            @test length(result.records) == 1
            @test value(result.records[1], 1) == 1
        end
    end

    @testset "test_should_sync_after_rollback" begin
        session(tc.driver) do sess
            tx = begin_transaction(sess)
            result = run(tx, "RETURN 1")
            rollback(tx)
            @test length(result.records) == 1
            @test value(result.records[1], 1) == 1
        end
    end

    @testset "test_errors_on_write_transaction" begin
        sess = session(tc.driver)
        @test_throws MethodError write_transaction(sess, (tx, args, kwargs)->run(tx, "CREATE (a:Thing {uuid:\$uuid})", uuid=args[1]), UUIDs.uuid4())
        close(sess)
    end

    @testset "test_errors_on_run_transaction" begin
        sess = session(tc.driver)
        tx = begin_transaction(sess)
        @test_throws MethodError run(tx, "CREATE (a:Thing {uuid:\$uuid})", uuid=uuid4())
        rollback(tx)
        close(sess)
    end

    @testset "test_errors_on_run_session" begin
        sess = session(tc.driver)
        close(sess)
        @test_throws ErrorException run(sess, "RETURN 1")
    end

    @testset "test_errors_on_begin_transaction" begin
        sess = session(tc.driver)
        close(sess)
        @test_throws ErrorException begin_transaction(sess)
    end

    @testset "test_large_values" begin
        for i in 1:6
            sess = session(tc.driver)
            run(sess, "RETURN '$(repeat('A', 2^20))'")
            close(sess)
            @test true # no error thrown
        end
    end

    teardown(tc.driver)
end

@testset "TransactionCommittedTestCase" begin
    tc = setup()

    @testset "test_errors_on_run" begin
        sess = session(tc.driver)
        tx = begin_transaction(sess)
        run(tx, "RETURN 1")
        commit(tx)
        @test_throws ErrorException run(tx, "RETURN 1")
    end

    teardown(tc.driver)
end

@testset "TransactionFunctionTestCase" begin
    tc = setup()

    @testset "test_simple_read" begin
        function work(tx, args, kwargs)
            return value(single(run(tx, "RETURN 1")))
        end

        session(tc.driver) do sess
            v = read_transaction(sess, work)
            @test v == 1
        end
    end

     @testset "test_read_with_arg" begin
        function work(tx, args, kwargs)
            return value(single(run(tx, "RETURN \$x", x=kwargs[:x])))
        end

        session(tc.driver) do sess
            v = read_transaction(sess, work, x=7)
            @test v == 7
        end
    end

     @testset "test_simple_write" begin
        function work(tx, args, kwargs)
            return value(single(run(tx, "CREATE (a {x: 7}) RETURN a.x")))
        end

        session(tc.driver) do sess
            v = write_transaction(sess, work)
            @test v == 7
        end
    end

     @testset "test_write_with_arg" begin
        function work(tx, args, kwargs)
            return value(single(run(tx, "CREATE (a {x: \$x}) RETURN a.x", x=kwargs[:x])))
        end

        session(tc.driver) do sess
            v = write_transaction(sess, work, x=7)
            @test v == 7
        end
    end

     @testset "test_write_with_arg_and_metadata" begin

        work = UnitOfWork(
            function(tx, args, kwargs)
                return value(single(run(tx, "CREATE (a {x: \$x}) RETURN a.x", x=kwargs[:x])))
            end,
            Dict("foo"=>"bar"),
            25
        )

        session(tc.driver) do sess
            v = write_transaction(sess, work, x=7)
            @test v == 7
        end
    end

    teardown(tc.driver)
end


text:
# Weight functions
# ----------------
# A good overview of these can be found at:
# http://research.microsoft.com/en-us/um/people/zhang/INRIA/Publis/Tutorial-Estim/node24.html

"""
An m-estimator is a cost/loss function used in modified (weighted) least squares
problems of the form:
    min ∑ᵢ ρ(rᵢ)
"""
abstract type MEstimator
end

# Default implementations from definintion of cost-function rho
"The cost (a.k.a. loss) function ρ for the M-estimator"
function estimator_rho end
"The derivative of the cost (a.k.a. loss) function ψ for the M-estimator"
function estimator_psi end
# TODO m_estimator_psi(r,est::MEstimator) = first derivative of m_estimator_rho(r,est) w.r.t. r
"""
The weight function, w, for the M-estimator, to be used for modifying the normal
equations of a least-square problem
"""
estimator_weight(r,est::MEstimator) = m_estimator_psi(r,est) ./ r
"""
The square root of the weight function, sqrt(w), for the M-estimator, to be used
for modifying a least-squares problem
"""
estimator_sqrtweight(r,est::MEstimator) = sqrt.(m_estimator_weight(r,est))
isconvex(::MEstimator) = false

"The (convex) L2 M-estimator is that of the standard least squares problem."
struct L2Estimator <: MEstimator; end
L2Estimator(width) = L2Estimator()
estimator_rho(r,::L2Estimator) = 0.5*r.^2
estimator_psi(r,::L2Estimator) = r
estimator_weight(r,::L2Estimator) = ones(size(r))
estimator_sqrtweight(r,::L2Estimator) = ones(size(r))
isconvex(::L2Estimator) = true

"""
The standard L1 M-estimator takes the absolute value of the residual, and is
convex but non-smooth.
"""
struct L1Estimator <: MEstimator; end
L1Estimator(width) = L1Estimator()
estimator_rho(r,::L1Estimator) = abs.(r)
estimator_psi(r,::L1Estimator) = sign.(r)
estimator_weight(r,::L1Estimator) = 1.0 ./ abs.(r)
estimator_sqrtweight(r,::L1Estimator) = abs.(r) .^ (-0.5)
isconvex(::L1Estimator) = true

"""
The convex L1-L2 estimator interpolates smoothly between L2 behaviour for small
residuals and L1 for outliers.
"""
struct L1L2Estimator <: MEstimator; width::Float64; end
estimator_rho(r,::L1L2Estimator) = 2.0 * (sqrt.(1.0 .+ 0.5*r.*r) .- 1.0)
estimator_psi(r,::L1L2Estimator) = r ./ sqrt.(1.0 .+ 0.5*r.*r)
estimator_weight(r,::L1L2Estimator) = 1.0 ./ sqrt.(1 .+ 0.5*r.*r)
estimator_sqrtweight(r,::L1L2Estimator) = (1.0 .+ 0.5*r.*r) .^ (-1/4)
isconvex(::L1L2Estimator) = true

"""
The convex Huber estimator switches from between quadratic and linear cost/loss
function at a certain cutoff.
"""
struct HuberEstimator <: MEstimator; width::Float64; end
function estimator_rho(r,est::HuberEstimator)
    rho = 0.5*r.^2
    absr = abs.(r)
    rho[absr .> est.width] = est.width * (absr[absr .> est.width] .- 0.5*est.width)
    return rho
end
function estimator_psi(r,est::HuberEstimator)
    psi = copy(r)
    absr = abs.(r)
    psi[absr .> est.width] = est.width * sign.(r[absr .> est.width])
    return psi
end
function estimator_weight(r,est::HuberEstimator)
    w = ones(size(r))
    absr = abs.(r)
    w[absr .> est.width] = est.width ./ absr[absr .> est.width]
    return w
end
function estimator_sqrtweight(r,est::HuberEstimator)
    w = ones(size(r))
    absr = abs.(r)
    w[absr .> est.width] = sqrt.(est.width ./ absr[absr .> est.width])
    return w
end
isconvex(::HuberEstimator) = true

"""
The (convex) "fair" estimator switches from between quadratic and linear
cost/loss function at a certain cutoff, and is C3 but non-analytic.
"""
struct FairEstimator <: MEstimator; width::Float64; end
estimator_rho(r,est::FairEstimator) = est.width^2 * (abs.(r)/est.width - log.(1.0 .+ abs.(r)/est.width))
estimator_psi(r,est::FairEstimator) = r ./ (1.0 .+ abs.(r)/est.width)
estimator_weight(r,est::FairEstimator) = 1.0 ./ (1.0 .+ abs.(r)/est.width)
estimator_sqrtweight(r,est::FairEstimator) = 1.0 ./ sqrt.(1.0 .+ abs.(r)/est.width)
isconvex(::FairEstimator) = true

"""
The non-convex Cauchy estimator switches from between quadratic behaviour to
logarithmic tails. This rejects outliers but may result in mutliple minima.
"""
struct CauchyEstimator <: MEstimator; width::Float64; end
estimator_rho(r,est::CauchyEstimator) = 0.5*est.width^2 * log.(1.0 .+ r.*r/(est.width*est.width))
estimator_psi(r,est::CauchyEstimator) = r ./ (1.0 .+ r.*r/(est.width*est.width))
estimator_weight(r,est::CauchyEstimator) = 1.0 ./ (1.0 .+ r.*r/(est.width*est.width))
estimator_sqrtweight(r,est::CauchyEstimator) = 1.0 ./ sqrt.(1.0 .+ r.*r/(est.width*est.width))
isconvex(::CauchyEstimator) = false

"""
The non-convex Geman-McClure for strong supression of outliers and does not guarantee a unique solution
"""
struct GemanEstimator <: MEstimator; end
estimator_rho(r,est::GemanEstimator) = 0.5*r.*r ./(1.0 .+ r.*r)
estimator_psi(r,est::GemanEstimator) = r ./ (1.0 .+ r.*r).^2
estimator_weight(r,est::GemanEstimator) = 1.0 ./ (1.0 .+ r.*r).^2
estimator_sqrtweight(r,est::GemanEstimator) = 1.0 ./ (1.0 .+ r.*r)
isconvex(::GemanEstimator) = false

"""
The non-convex Welsch for strong supression of outliers and does not guarantee a unique solution
"""
struct WelschEstimator <: MEstimator; width::Float64; end
estimator_rho(r,est::WelschEstimator) = 0.5 * est.width^2 *( 1.0 .- exp.(-(r/est.width).^2) )
estimator_psi(r,est::WelschEstimator) = r.* exp.(-(r/est.width).^2)
estimator_weight(r,est::WelschEstimator) = exp.(-(r/est.width).^2)
estimator_sqrtweight(r,est::WelschEstimator) = exp.(-0.5*(r/est.width).^2)
isconvex(::WelschEstimator) = false

"""
The non-convex Tukey biweight estimator which completly suppress the outliers,
and does not guarantee a unique solution
"""
struct TukeyEstimator <: MEstimator; width::Float64; end
function estimator_rho(r,est::TukeyEstimator)
    rho = est.width^2/6.0 * (1.0 .- (1.0 .- (r/est.width).^2 ).^3)
    rho[est.width .< abs.(r)] .= est.width^2/6.0
    return rho
end
function estimator_psi(r,est::TukeyEstimator)
    psi = r.* (1.0 .- (r/est.width).^2).^2
    psi[est.width .< abs.(r)] .= 0.0
    return psi
end
function estimator_weight(r,est::TukeyEstimator)
    w = (1.0 .- (r/est.width).^2).^2
    w[est.width .< abs.(r)] .= 0.0
    return w
end
function estimator_sqrtweight(r,est::TukeyEstimator)
    w = (1.0 .- (r/est.width).^2)
    w[est.width .< abs.(r)] .= 0.0
    return w
end
isconvex(::TukeyEstimator) = true

"""
A custom M-Estimator which composes multiple estimators on different ranges.

Construct with syntax:
    MultiEstimator(MEstimator => Range, MEstimator => Range, ...)

Unfilled ranges will automatically use the standard L2-estimator.
"""
struct MultiEstimator{Estimators<:Tuple,Ranges<:Tuple} <: MEstimator
    est::Estimators
    rng::Ranges

    function MultiEstimator{Estimators, Ranges}(est::Estimators, rng::Ranges) where {Estimators<:Tuple, Ranges<:Tuple}
        if length(est) != length(rng)
            error("Must have same number of estimators and ranges")
        end
        for i = 1:length(est)
            if !isa(est[i],MEstimator)
                error("Expecting an M-estimator, got $(est[i]) of type $(typeof(est[i]))")
            end
            if !hasmethod(Base.getindex,(Vector{Float64},typeof(rng[i])))
                error("Expecting an indexible range, got $(rng[i]) of type $(typeof(rng[i]))")
            end
        end
        new{Estimators, Ranges}(est,rng)
    end
end

# Convenience constructors for up to 4 different ranges
MultiEstimator(x::Pair{Est1,T1}) where {Est1<:MEstimator,T1} = MultiEstimator{Tuple{Est1},Tuple{T1}}((x.first,), (x.second,))
MultiEstimator(x1::Pair{Est1,T1},x2::Pair{Est2,T2}) where {Est1<:MEstimator,T1,Est2<:MEstimator,T2} = MultiEstimator{Tuple{Est1,Est2},Tuple{T1,T2}}((x1.first,x2.first), (x1.second,x2.second))
MultiEstimator(x1::Pair{Est1,T1},x2::Pair{Est2,T2},x3::Pair{Est3,T3}) where {Est1<:MEstimator,T1,Est2<:MEstimator,T2,Est3<:MEstimator,T3} = MultiEstimator{Tuple{Est1,Est2,Est3},Tuple{T1,T2,T3}}((x1.first,x2.first,x3.first), (x1.second,x2.second,x3.second))
MultiEstimator(x1::Pair{Est1,T1},x2::Pair{Est2,T2},x3::Pair{Est3,T3},x4::Pair{Est4,T4}) where {Est1<:MEstimator,T1,Est2<:MEstimator,T2,Est3<:MEstimator,T3,Est4<:MEstimator,T4} = MultiEstimator{Tuple{Est1,Est2,Est3,Est4},Tuple{T1,T2,T3,T4}}((x1.first,x2.first,x3.first,x4.first), (x1.second,x2.second,x3.second,x4.second))
function estimator_rho(r,est::MultiEstimator)
    rho = 0.5*r.^2
    for i = 1:length(est.est)
        rho[est.rng[i]] = estimator_rho(r[est.rng[i]],est.est[i])
    end
    return rho
end
function estimator_psi(r,est::MultiEstimator)
    psi = r
    for i = 1:length(est.est)
        psi[est.rng[i]] = estimator_psi(r[est.rng[i]],est.est[i])
    end
    return psi
end
function estimator_weight(r,est::MultiEstimator)
    w = ones(size(r))
    for i = 1:length(est.est)
        w[est.rng[i]] = estimator_weight(r[est.rng[i]],est.est[i])
    end
    return w
end
function estimator_sqrtweight(r,est::MultiEstimator)
    w = ones(size(r))
    for i = 1:length(est.est)
        w[est.rng[i]] = estimator_sqrtweight(r[est.rng[i]],est.est[i])
    end
    return w
end
function isconvex(r,est::MultiEstimator)
    out = true
    for i = 1:length(est.est)
        out = out && isconvex(est.est[i])
    end
    return out
end

"""
Recreate an M-Estimator of the same type using the median absolute
deviation of the residual
"""
refit_estimator(Est::Union{T,Type{T}}, res, factor=1.0) where {T<:MEstimator} = T(factor*1.43*StatsBase.mad(res, normalize=false))
function refit_estimator(est::T, res, factor=1.0) where T<:MultiEstimator
    tmp = ntuple(i->(typeof(est.est[i])(factor*1.43*StatsBase.mad(res[est.rng[i]], normalize=false)) => est.rng[i]),length(est.est))
    return MultiEstimator(tmp...)
end

ext: jl


text:
#=

This file contains code for calling Matlab solvers from inside Julia.
To the extent possible, it uses the same interfaces as the solvers in Laplacians.

They were mainly written to allow speed testing of linear solvers.

=#

using MATLAB
test_var_ = 1
@mput test_var_
mat"""
disp('Connected to Matlab!')
"""

"""
    f = matlab_ichol(sdd)

Computes a zero-fill incomplete Cholesky factorization of sdd, and returns a function that solves systems of linear equations in it.  Uses the symrcm ordering of the matrix.  Requires that MATLAB.jl be included.
"""
function matlab_ichol(sdd)

    if ~isdefined(:MATLAB)
        error("Type using MATLAB before calling this routine.")
    end

    n = sdd.n

    @mput sdd

    mat"""
    p = symrcm(sdd);
    laperm = sdd(p,p);
    L = ichol(laperm)
    """

    @mget p
    @mget L

    p = convert(Array{Int,1},vec(p))
    pi = zeros(Int,n)
    pi[p] = collect(1:n)

    Ltri = LowerTriangular(L)
    Utri = Ltri'

    # Create the solver function
    f = function(b::Array{Float64,1})

        y = Utri \ (Ltri \ b[p])
        return y[pi]

    end

    return f
end

"""
    f = matlab_ichol_sddm(sddm; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)

Computes a zero-fill incomplete Cholesky factorization of sddm, and returns a function that solves systems of linear equations in it.  Uses the symrcm ordering of the matrix.  Requires that MATLAB.jl be included.
"""
function matlab_ichol_sddm(sddm::SparseMatrixCSC; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[])

    F = matlab_ichol(sddm)
    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

    f(b;tol=tol_,maxits=maxits_, maxtime=maxtime_, verbose=verbose_,pcgIts=pcgIts_) = pcg(sddm, b, F, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)

    return f

end

"""
    f = matlab_ichol_lap(A; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)

Wrapper for matlab_ichol_sddm.
"""
matlab_ichol_lap = lapWrapSDDM(matlab_ichol_sddm)

"""
    x = matlabCmgSolver(mat, b; tol::Real=1e-6, maxits=10000)

This runs Koutis's CMG solver.  You must have installed the solver, and it must be on Matlab's default path.  This routine does not implement all of our preferred interface.  Use the same solver for sddm and Laplacian matrices.
"""
function matlabCmgSolver(mat::SparseMatrixCSC, b; tol::Real=1e-6, maxits=10000, verbose=false, params...)

    maxits = min(10000,maxits)

    @mput mat
    @mput b
    @mput tol
    @mput maxits

    mat"""
    pfun = cmg_sdd(mat);

    [x,flag,relres,iter] = pcg(mat, b, tol, maxits, pfun);
    """

    @mget iter
    @mget x

    if verbose
        println("PCG stopped after $(iter) iterations.")
    end

    return x
end

"""
    solver = matlabCmgSolver(mat; tol::Real=1e-6, maxits=10000)

This runs Koutis's CMG solver.  You must have installed the solver, and it must be on Matlab's default path.  This routine does not implement all of our preferred interface.  Use the same solver for sddm and Laplacian matrices.

Note that this does not build the solver.  Rather, it just makes sure that the call to `solver` will.
"""
function matlabCmgSolver(mat::SparseMatrixCSC; tol::Real=1e-6, maxits=1000, params...)

    tol_=tol
    maxits_=maxits

    f(b; tol=tol_, maxits=maxits_, params...) = matlabCmgSolver(mat, b; tol=tol, maxits=maxits)

    return f
end

"""
    f = matlabCmgSolverHalf(mat; tol::Real=1e-6, maxits=10000, verbose=false)

This builds Koutis's CMG solver.
You must have installed the solver, and it must be on Matlab's default path.
This routine does not implement all of our preferred interface.
Use the same solver for sddm and Laplacian matrices.

This returns a function that you can use to solve systems in mat thereafter.
However, if you interact with Matlab before calling f, objects that it needs
in Matlab's namespace can be overwritten.  Use with caution!
"""
function matlabCmgSolverHalf(mat::SparseMatrixCSC; tol::Real=1e-6, maxits=10000,
  verbose=false, pcgIts=Int[], params...)

    # check if this is an ajacency matrix instead of a Laplacian
    if maximum(diag(mat)) == 0
      mat = lap(mat)
    end

    maxits = min(10000,maxits)

    @mput mat

    t1 = time()

    mat"""
    pfun = cmg_sdd(mat);
    """

    if verbose
      println("Build time ", time() - t1)
    end

    tol_=tol
    maxits_=maxits
    verbose_=verbose
    pcgIts_ =pcgIts

    f = function(b; tol=tol_, maxits=maxits_, verbose=verbose_, pcgIts=pcgIts_)
      @mput b
      @mput tol
      @mput maxits

      mat"""
      [x,flag,relres,iter] = pcg(mat, b, tol, maxits, pfun);
      """

      @mget iter
      @mget x

      if verbose
          println("PCG stopped after $(iter) iterations.")
      end
      if length(pcgIts) > 0
          pcgIts[1] = iter
      end

      return x
    end

    return f
end

"""
    solver = matlabCmgLap(a; tol::Real=1e-6, maxits=10000)

This runs Koutis's CMG solver.  You must have installed the solver, and it must be on Matlab's default path.  This routine does not implement all of our preferred interface.

The input `a` should be the adjacency matrix.  This will find connected components, and then call the sddm version of the solver.

Note that this does not build the solver.  Rather, it just makes sure that the call to `solver` will.
"""
matlabCmgLap = lapWrapSDDM(matlabCmgSolver)


text:
#########################################################################################################################################
############################# baseextension.jl features some functionality that ought to be in Julia.Base. ##############################
#########################################################################################################################################
#
# Functions:
#       +()             -  addition for tuples
#       allcomb()       -  take all possible combinations of the elements of a set of n vectors
#       ndgrid()        -  construct an n-dimensional cartesian grid using n input vectors with coordinates (ndim-generalization of meshgrid in MatLab)
#       blockmat()      -  construct an abstract matrix of type T with N abstract matrices of type T on its diagonal
#
#########################################################################################################################################
#
# File part of the Julia IsoGeometric Analysis toolbox
# R.R.Hiemstra
# 14-09-2015
#
#########################################################################################################################################

import Base.getindex, Base.setindex!, Base.size, Base.length, Base.similar, Base.*, Base.+, Base.-

# multiplication by a scalar for tuples
*(a::S,q::NTuple{N,T}) where {N,S<:Real,T<:Real} = ntuple(i -> a*q[i], N)
*(q::NTuple{N,T},a::S) where {N,S<:Real,T<:Real} = ntuple(i -> a*p[i], N)

# addition for tuples
+(p::NTuple{N,T},q::NTuple{N,T}) where {N,T<:Real} = ntuple(i -> p[i]+q[i], N)
-(p::NTuple{N,T},q::NTuple{N,T}) where {N,T<:Real} = ntuple(i -> p[i]-q[i], N)
# input:
#    p :: Tuple{N,T}  - ntuple of real numbers of type T
#    q :: Tuple{N,T}  - ntuple of real numbers of type T
# output:
#    r :: Tuple{N,T}  - ntuple of real numbers of type T

# adding a real number to a tuple
+(p::NTuple{N,S},q::T) where {N,S<:Real,T<:Real} = ntuple(i -> p[i]+q, N)
-(p::NTuple{N,S},q::T) where {N,S<:Real,T<:Real} = ntuple(i -> p[i]-q, N)
+(p::S,q::NTuple{N,T}) where {N,S<:Real,T<:Real} = ntuple(i -> p[i]+q, N)
-(p::S,q::NTuple{N,T}) where {N,S<:Real,T<:Real} = ntuple(i -> p[i]-q, N)

"create linearly spaced vector"
linearspace(a,b,n) = collect(range(a, stop=b, length=n))

"create linearly spaced Matrix"
function linearspace(v::Vector{T},k) where {T<:Real}
    m = length(v)-1
    V = zeros(k,m)
    [V[:,i] = linearspace(v[i],v[i+1],k) for i in 1:m]
    return V
end

"Take all possible combinations of the elements of a set of n vectors"
allcomb(v::AbstractVector) = copy(v)
allcomb(v1::AbstractVector{T}, v2::AbstractVector{T}) where {T} = [repeat(v1, length(v2), 1) repeat(reshape(v2,1,length(v2)), size(v1,1), 1)[:]]
allcomb(a::AbstractMatrix{T}, v::AbstractVector{T}) where {T} = [repeat(a, length(v), 1) repeat(reshape(v,1,length(v)), size(a,1), 1)[:]]

"""
take all possible combinations of the elements of a set of n vectors

## input:
    vs :: AbstractVector{T}...  -  set of n abstract vectors of type T

## output:
    out :: AbstractMatrix{T}    -  one output of abstract matrix of type T which contains all possible combinations of the
                                   elements of the n input vectors
"""
function allcomb(vs::AbstractVector{T}...) where {T}
    n = length(vs)
    out = allcomb(vs[1])
    for i in 2:n
        out = allcomb(out,vs[i])
    end
    return out
end

"Construct an n-dimensional cartesian grid using n input vectors with coordinates (ndim-generalization of meshgrid in MatLab)"
ndgrid(v::AbstractVector) = copy(v)

function ndgrid(v1::AbstractVector{T}, v2::AbstractVector{T}) where {T}
    m, n = length(v1), length(v2)
    v1 = reshape(v1, m, 1)
    v2 = reshape(v2, 1, n)
    (repeat(v1, 1, n), repeat(v2, m, 1))
end

function ndgrid_fill(a, v, s, snext)
    for j = 1:length(a)
        a[j] = v[div(rem(j-1, snext), s)+1]
    end
end

"""
Construct an n-dimensional cartesian grid using n input vectors with coordinates (ndim-generalization of meshgrid in MatLab)

## input:
    vs :: AbstractVector{T}...  -  set of n abstract vectors of type T representing univariate coordinates

## output:
    out :: AbstractMatrix{T}...    -  set of n abstract matrices of type T representing an n-dimensional grid
"""
function ndgrid(vs::AbstractVector{T}...) where {T}
    n = length(vs)
    sz = map(length, vs)
    out = ntuple(n, i->Array(T, sz))
    s = 1
    for i=1:n
        a = out[i]::Array
        v = vs[i]
        snext = s*size(a,i)
        ndgrid_fill(a, v, s, snext)
        s = snext
    end
    out
end

"""
Construct an abstract matrix of type T with N abstract matrices of type T on its diagonal

## input:
    A :: AbstractMatrix{T}...  -  set of n abstract matrices of type T

## output:
    B :: AbstractMatrix{T}     -  abstract matrix of type T with the N abstract matrices on its diagonal
"""
function blockmat(A::AbstractMatrix{T}...) where {T}
    # determine size of matrices
    n = [0,cumsum([size(A[k],1) for k in 1:length(A)])]
    m = [0,cumsum([size(A[k],2) for k in 1:length(A)])]

    # allocate space for blockmatrix
    B = similar(A[1],n[end],m[end]); B[:] = 0

    # fill blockmatrix
    for k in 1:length(A)
        B[n[k]+1:n[k+1],m[k]+1:m[k+1]] = A[k]
    end

    return B
end

max_stars_repo_name: michakraus/iga-demos

max_issues_repo_name: michakraus/iga-demos

max_forks_repo_name: michakraus/iga-demos


text:
using SententialDecisionDiagrams
const SDD = SententialDecisionDiagrams

var_count = 4
var_order = [2,1,4,3]
vtree_type = "balanced"

vtree = SDD.vtree(var_count, var_order, vtree_type)
manager = SDD.sdd_manager(vtree)

println("constructing SDD ... ")
a, b, c, d = [SDD.literal(i,manager) for i in 1:4]
α = (a & b) | (b & c) | (c & d)
println("done")

println("saving sdd and vtree ... ")
SDD.dot("$(@__DIR__)/output/sdd.dot",α)
SDD.dot("$(@__DIR__)/output/vtree.dot",vtree)
println("done")


text:
module Bodies
    import LinearAlgebra
    using NBodyLeapFrog: V3

    const mass_solar = 1.988544e30 # in kg
    const mass_mercury = 3.302e23
    const mass_venus = 48.685e23
    const mass_earth = 5.97219e24
    const mass_mars = 6.4185e23
    const mass_jupyter = 1898.13e24
    const mass_saturn = 5.68319e26
    const mass_neptune = 102.41e24
    const mass_uranus = 86.8103e24
    const mass_pluto= 1.307e22

    const mass_planetX = 3.2e24 # Planet refers to an additional planet I can add to the solar
    # system to test if it is maximally packed. (The solar system is not
    # maximally packed. There is room to add another planet and preserve the
    # stable orbit of all planets.)

    #  unit of time -> year, unit of space -> AU
    const year = 3.154e7 #year in seconds
    const AU = 1.496e11 #in m

    const G_standard = 6.67259e-11 # in m^3/(kg-s^2)
    const G_year_AU = G_standard*(1/AU)^3/(1/mass_solar*(1/year)^2)

    struct Body{T}
        r::V3{T}
        v::V3{T}
        m::T
    end

    const mercury = Body(
        V3(8.887985138765460E-02, -4.426150338141062E-01, -4.475716356484761E-02),
        # velocity in AU/day
        V3(2.190877912081542E-02, 7.161568136528000E-03, -1.425929443086507E-03),
        mass_mercury / mass_solar
    )

    const venus = Body(
        V3(4.043738093622098E-02, -7.239789211502183E-01, -1.241560658530024E-02),
        V3(2.005742309538389E-02, 1.141448268256643E-03, -1.142174441569258E-03),
        mass_venus/mass_solar
    )

    const earth =  Body(
        V3(-2.020844529756663E-02, -1.014332737790859E+00, -1.358267619371298E-05),
        V3(1.692836723212859E-02, -3.484006532982474E-04, 6.028542314557626E-07),
        mass_earth/mass_solar
    )

    # target 1998
    #  -3.466334931755365E-02 -1.013773181327570E+00  2.111689861662178E-04
    #   1.692129864984556E-02 -5.252811129268817E-04  3.987686870581435E-07

    const mars = Body(
        V3(7.462481663749645E-01, -1.181663652521456E+00, -4.321921404013512E-02),
        V3(1.235610918162121E-02, 8.680869489377649E-03 , -1.220500608452554E-04),
        mass_mars/mass_solar
    )

    const planetX = Body(
        V3(0.0, 2.06, 0.0),
        V3(1.235610918162121E-02, 0.0, 0.0),
        mass_planetX / mass_solar
    )

    const jupyter = Body(
        V3(3.384805319103406E+00, 3.658805636759595E+00, -9.100441946210819E-02),
        V3(-5.634671617093230E-03, 5.479180979634376E-03, 1.034981407898108E-04),
        mass_jupyter/mass_solar
    )

    const saturn = Body(
        V3(-1.083899692644216E-01, -1.003995196286016E+01, 1.793391553155583E-01),
        V3(5.278410787728323E-03, -7.712342079566598E-05, -2.084447335785041E-04),
        mass_saturn / mass_solar
    )

    const neptune = Body(
        V3(4.675566709791660E+00, -2.985428200863175E+01, 5.070034142531887E-01),
        V3(3.080716380724798E-03, 5.030733458293977E-04, -8.101711269674541E-05),
        mass_neptune / mass_solar
    )

    const uranus = Body(
        V3(-2.693448460292631E-01, -1.927606446869220E+01, -6.808868692550485E-02),
        V3(3.903100242621723E-03, -2.380111092360100E-04, -5.164025224695875E-05),
        mass_uranus / mass_solar
    )

    const pluto = Body(
        V3(-2.129074273328636E+01, -1.896633337434039E+01, 8.187955378677129E+00),
        V3(2.276295756013608E-03, -2.670481848836963E-03, -3.669545371032554E-04),
        mass_pluto / mass_solar
    )

    # target  -1.156541154581570E+01 -2.704864218000164E+01  6.239749761161465E+00
    #   2.964408290188142E-03 -1.722224413824548E-03 -6.839434010481107E-04

    const sun = Body(
        V3(-3.430031536367300E-03, 1.761881027012596E-03, 1.246691303879918E-05),
        V3(3.433119412673547E-06, -5.231300927361546E-06, -2.972974735550750E-08),
        1.0
    )

    # Arrays for modeling the solar system
    const set = [sun, mercury, venus, earth, mars, jupyter, saturn, uranus, neptune, pluto]
    const set_packed = [sun, mercury, venus, earth, mars, planetX, jupyter, saturn, uranus, neptune, pluto]

    const dayToYear = 365.25

    function chunit_day2year(body::Body)
        Body(body.r, body.v * dayToYear, body.m)
    end
end


text:
@testset "conversions.jl" begin
    gp = GP(SqExponentialKernel())
    x = rand(10)
    f = gp(x, 1e-9)
    mu, var = mean_and_var(f)
    z = rand(10)

    # `SymBand`
    @test convert_arguments(SymBand, z, f) == (map(Point2f, z, mu), map(sqrt, var))

    # `GPSample`
    @test convert_arguments(GPSample, z, f) == (z, f)

    # Default fallback
    for P in (Lines, LineSegments, Scatter)
        @test convert_arguments(P, z, f) == convert_arguments(P, z, mu)
    end

    # Fallback for `Band`
    lower = mu .- sqrt.(var)
    upper = mu .+ sqrt.(var)
    @test convert_arguments(Band, z, f) == (map(Point2f, z, lower), map(Point2f, z, upper))

    scale = rand()
    lower = mu .- scale .* sqrt.(var)
    upper = mu .+ scale .* sqrt.(var)
    @test convert_arguments(Band, z, f; bandscale=scale) ==
          (map(Point2f, z, lower), map(Point2f, z, upper))

    # Fallback for `AbstractGP` and `FiniteGP`
    for P in (Lines, LineSegments, Scatter, Band, SymBand, GPSample)
        @test convert_arguments(P, x, gp) == convert_arguments(P, x, f)
        @test convert_arguments(P, f) == convert_arguments(P, x, f)
    end
end


text:
using DataStructures: CircularBuffer
include("AlgorithmicEnv.jl")

"""
Task is to return every nth character from the input tape.
http://arxiv.org/abs/1511.07275
"""
mutable struct DuplicatedInputEnv <: TapeAlgorithmicEnv
    MIN_REWARD_SHORTFALL_FOR_PROMOTION::Float32
    base::Int
    episode_total_reward::Float32
    reward_shortfalls::CircularBuffer{Float32}
    charmap::Array{Char, 1}
    min_length::Int
    action_space::TupleSpace
    observation_space::Discrete

    MOVEMENTS::NTuple
    READ_HEAD_START::Int

    duplication::Int

    target::Array{Int8, 1}  # target tape
    input_data::Array{Int8, 1}
    time::Int8
    read_head_position::Int8  # current position of the read head.
    write_head_position::Int8
    last_action
    last_reward::Float32
end

DuplicatedInputEnv(base::Int=5, duplication::Int=2) =
    DuplicatedInputEnv(-1f0, # MIN_REWARD_SHORTFALL_FOR_PROMOTION
            base,
            0f0, # episode_total_reward
            CircularBuffer{Float32}(10),  # reward_shortfalls
            push!(['A'+i for i=0:base-1], ' '),  # charmap
            2,  # starting min_length
            TupleSpace([Discrete(2), Discrete(2), Discrete(base)]),  # action_space
            Discrete(base+1),  # observation_space
            (:left, :right), 1,  # MOVEMENTS and READ_HEAD_START
            duplication,
            Int8[],
            Int8[],
            0, 1, 1, nothing, 0.0)

function generate_input_data!(env::DuplicatedInputEnv, size_)
    res = []
    size_ < env.duplication &&
        (size_ = env.duplication)
    for i=1:div(size_, env.duplication)
        char = rand(1:env.base)
        for _ in 1:env.duplication
            push!(res, char)
        end
    end
    env.input_data = res
end

target_from_input_data!(env::DuplicatedInputEnv) =
    env.target = [env.input_data[i] for i=1:env.duplication:length(env.input_data)]


text:
using BenchmarkTools
BLAS.set_num_threads(4)

verb, samples, seconds = 0, 5, 60^2
suite = BenchmarkGroup()

suite["SparseDeconvolution"] = BenchmarkGroup()
include("SparseDeconvolution.jl")
suite["SparseDeconvolution"]["FullMatrix"] = SparseDeconvolution.benchmark(verb = verb, samples = samples, seconds = seconds)
suite["SparseDeconvolution"]["MatrixFree"] = SparseDeconvolution.benchmarkMatrixFree(verb = verb, samples = samples, seconds = seconds)

suite["LineSpectraEstimation"] = BenchmarkGroup()
include("LineSpectraEstimation.jl")
suite["LineSpectraEstimation"]["FullMatrix"] = LineSpectraEstimation.benchmark(verb = verb, samples = samples, seconds = seconds)
suite["LineSpectraEstimation"]["MatrixFree"] = LineSpectraEstimation.benchmarkMatrixFree(verb = verb, samples = samples, seconds = seconds)

suite["MatrixDecomposition"] = BenchmarkGroup()
include("MatrixDecomposition.jl")
suite["MatrixDecomposition"] = MatrixDecomposition.benchmark(verb = verb, samples = samples, seconds = seconds)

suite["DNN"] = BenchmarkGroup()
include("DNN.jl")
suite["DNN"] = DNN.benchmark(verb = verb, samples = samples, seconds = seconds)

suite["TotalVariation"] = BenchmarkGroup()
include("TotalVariation.jl")
suite["TotalVariation"] = TotalVariation.benchmark(verb = verb, samples = samples, seconds = seconds)

suite["AudioDeclipping"] = BenchmarkGroup()
include("AudioDeclipping.jl")
suite["AudioDeclipping"] = AudioDeclipping.benchmark(verb = verb, samples = 1, seconds = seconds)

println("\n")
showall(median(suite))
println("\n")
#showall(memory(suite))


text:
function norm_two_constraint(ps_m::CanonicalModel,
                            range_data::Array{Tuple{String,Float64},1},
                            time_range::UnitRange{Int64},
                            cons_name::Symbol,
                            var_names::Tuple{Symbol, Symbol})

    ps_m.constraints[cons_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, [r[1] for r in range_data], time_range)

    for t in time_range, r in range_data

            ps_m.constraints[cons_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, r[2]^2 == ps_m.variables[var_names[1]][r[1], t] + ps_m.variables[var_names[2]][r[1], t])

    end

    return

end


text:
struct CalciumDynamics{T<:Number} <: SpeciesDynamics{Calcium}
    τCa::T
    Ca∞::T
    calc_multiplier::T # f * area = 14.96 * 0.0628
end

function (p::CalciumDynamics)(n::Neuron, vars, varnames, currents)
    τCa, Ca∞, calc_multiplier, Cₘ = get_parameters(p)
    Ca = vars[findfirst(x -> x == Calcium, varnames)]
    D(Ca) ~ (1 / τCa) * (-Ca + Ca∞ + (calc_multiplier * currents / Cₘ))
end

get_parameters(::CalciumDynamics) = @parameters τCa, Ca∞, calc_multiplier, Cₘ
default_params(p::CalciumDynamics, n::Neuron, vars, varnames) = Dict(
    get_parameters(p) .=> (p.τCa, p.Ca∞, p.calc_multiplier, capacitance(n.geometry))
)

struct CaReversalDynamics <: SpeciesDynamics{Calcium} end

function (p::CaReversalDynamics)(n::Neuron, vars, varnames, currents)
    Ca = vars[findfirst(x -> x == Calcium, varnames)]
    ECa = vars[findfirst(x -> x == Reversal{Calcium}, varnames)]
    return ECa ~ (500.0) * (8.6174e-5) * (283.15) * (log(max((3000.0 / Ca), 0.001)))
end


text:
"""
    A Solver takes a System and solves it.
"""
abstract Solver

" Should be overridden by the Solver implementations to initialize their data "
function initialize!(solver::Solver)
end

" Must be overriden by Solver implementations to actually do the solve "
function solve!(solver::Solver)
    throw(MethodError(solve!, solver))
end

" Assemble the FE problem "
function assemble!(solver::Solver)
end

" Solvers that solve using a matrix and rhs "
abstract ImplicitSolver <: Solver

" Solvers that use dense matrices "
abstract DenseImplicitSolver <: ImplicitSolver
    # " The System the solve will be done for "
    # system::System

    # " The 'stiffness matrix' "
    # mat::Matrix{Float64}

    # " The RHS (force vector) "
    # rhs::Vector{Float64}

    # " The solution vector "
    # solution::Vector{Float64}

    # " The ghosted solution vector (for Assembly) "
    # ghosted_solution::Vector{Float64}

    # " Whether or not this Solver has been initialized "
    # initialized::Bool

" Initializes the matrix and rhs to the correct size "
function initialize!(solver::DenseImplicitSolver)
    @assert solver.system.initialized

    n_dofs = solver.system.n_dofs

    solver.mat = zeros(Float64, n_dofs, n_dofs)
    solver.rhs = zeros(Float64, n_dofs)
    solver.solution = zeros(Float64, n_dofs)
    solver.ghosted_solution = zeros(Float64, n_dofs)

    solver.initialized = true
end

" Update the ghosted_solution vector from the actual solution "
function updateGhostedSolution!(solver::DenseImplicitSolver)
    solver.ghosted_solution = solver.solution
end


text:
@testsuite "broadcasting" AT->begin
    broadcasting(AT)
    vec3(AT)
end

test_idx(idx, A::AbstractArray{T}) where T = A[idx] * T(2)

function testv3_1(a, b)
    x = a .* b
    y =  x .+ b
    @inbounds return y[1]
end

function testv3_2(a, b)
    x = a .* b
    y =  x .+ b
    return y
end

# more complex function for broadcast
function test_kernel(a::T, b) where T
    c = a
    for i=1:5
        c += T(b)
        c = a * T(2)
    end
    return c
end

function broadcasting(AT)
    for ET in supported_eltypes()
        N = 10
        @testset "broadcast $ET" begin
            @testset "RefValue" begin
                cidx = rand(1:Int(N), 2*N)
                gidx = AT(cidx)
                cy = rand(ET, 2*N)
                gy = AT(cy)
                cres = fill(zero(ET), size(cidx))
                gres = AT(cres)
                gres .= test_idx.(gidx, Base.RefValue(gy))
                cres .= test_idx.(cidx, Base.RefValue(cy))
                @test Array(gres) == cres
            end

            @testset "Tuple" begin
                @test compare(AT, rand(ET, 3, N), rand(ET, 3, N), rand(ET, N), rand(ET, N), rand(ET, N)) do out, arr, a, b, c
                    broadcast!(out, arr, (a, b, c)) do xx, yy
                        xx + first(yy)
                    end
                end
            end

            @testset "Adjoint and Transpose" begin
                A = AT(rand(ET, N))
                A' .= ET(2)
                @test all(isequal(ET(2)'), A)
                transpose(A) .= ET(1)
                @test all(isequal(ET(1)), A)
            end

            ############
            # issue #27
            @test compare((a, b)-> a .+ b, AT, rand(ET, 4, 5, 3), rand(ET, 1, 5, 3))
            @test compare((a, b)-> a .+ b, AT, rand(ET, 4, 5, 3), rand(ET, 1, 5, 1))

            ############
            # issue #22
            dim = (32, 32)
            @test compare(AT, rand(ET, dim), rand(ET, dim), rand(ET, dim)) do tmp, a1, a2
                tmp .=  a1 .+ a2 .* ET(2)
            end

            ############
            # issue #21
            if ET in (Float32, Float64)
                @test compare((a1, a2)-> muladd.(ET(2), a1, a2), AT, rand(ET, dim), rand(ET, dim))
                #########
                # issue #41
                # The first issue is likely https://github.com/JuliaLang/julia/issues/22255
                # since GPUArrays adds some arguments to the function, it becomes longer longer, hitting the 12
                # so this wont fix for now
                @test compare(AT, rand(ET, dim), rand(ET, dim), rand(ET, dim), rand(ET, dim), rand(ET, dim), rand(ET, dim)) do a1, a2, a3, a4, a5, a6
                    @. a1 = a2 + (1.2) *((1.3)*a3 + (1.4)*a4 + (1.5)*a5 + (1.6)*a6)
                end

                @test compare(AT, rand(ET, dim), rand(ET, dim), rand(ET, dim), rand(ET, dim)) do u, uprev, duprev, ku
                    fract = ET(1//2)
                    dt = ET(1.4)
                    dt2 = dt^2
                    @. u = uprev + dt*duprev + dt2*(fract*ku)
                end
                @test compare((x) -> (-).(x), AT, rand(ET, 2, 3))

                @test compare(AT, rand(ET, dim), rand(ET, dim), rand(ET, dim), rand(ET, dim), rand(ET, dim), rand(ET, dim)) do utilde, gA, k1, k2, k3, k4
                    btilde1 = ET(1)
                    btilde2 = ET(1)
                    btilde3 = ET(1)
                    btilde4 = ET(1)
                    dt = ET(1)
                    @. utilde = dt*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4)
                end
            end

            @test compare((x) -> fill!(x, 1), AT, rand(ET, 3,3))
            @test compare((x, y) -> map(+, x, y), AT, rand(ET, 2, 3), rand(ET, 2, 3))

            @test compare((x) -> 2x, AT, rand(ET, 2, 3))
            @test compare((x, y) -> x .+ y, AT, rand(ET, 2, 3), rand(ET, 1, 3))
            @test compare((z, x, y) -> z .= x .+ y, AT, rand(ET, 2, 3), rand(ET, 2, 3), rand(ET, 2))

            @test compare(A -> A .= identity.(ET(10)), AT, rand(ET, 40, 40))
            @test compare(A -> test_kernel.(A, ET(10)), AT, rand(ET, 40, 40))
            @test compare(A -> A .* ET(10), AT, rand(ET, 40, 40))
            @test compare((A, B) -> A .* B, AT, rand(ET, 40, 40), rand(ET, 40, 40))
            @test compare((A, B) -> A .* B .+ ET(10), AT, rand(ET, 40, 40), rand(ET, 40, 40))
        end

        @testset "map! $ET" begin
            @test compare(AT, rand(2,2), rand(2,2)) do x,y
                map!(+, x, y)
            end
            @test compare(AT, rand(2), rand(2,2)) do x,y
                map!(+, x, y)
            end
            @test compare(AT, rand(2,2), rand(2)) do x,y
                map!(+, x, y)
            end
        end

        @testset "map $ET" begin
            @test compare(AT, rand(2,2), rand(2,2)) do x,y
                map(+, x, y)
            end
            @test compare(AT, rand(2), rand(2,2)) do x,y
                map(+, x, y)
            end
            @test compare(AT, rand(2,2), rand(2)) do x,y
                map(+, x, y)
            end
        end
    end

    @testset "0D" begin
        x = AT{Float64}(undef)
        x .= 1
        @test collect(x)[] == 1
        x /= 2
        @test collect(x)[] == 0.5
    end

    @testset "Ref" begin
        # as first arg, 0d broadcast
        @test compare(x->getindex.(Ref(x),1), AT, [0])

        void_setindex!(args...) = (setindex!(args...); return)
        @test compare(x->(void_setindex!.(Ref(x),1); x), AT, [0])

        # regular broadcast
        a = AT(rand(10))
        b = AT(rand(10))
        cpy(i,a,b) = (a[i] = b[i]; return)
        cpy.(1:10, Ref(a), Ref(b))
        @test Array(a) == Array(b)
    end

    @testset "stackoverflow in copy(::Broadcast)" begin
        copy(Base.broadcasted(identity, AT(Int[])))
    end
end

function vec3(AT)
    @testset "vec 3" begin
        N = 20

        xc = map(x-> ntuple(i-> rand(Float32), Val(3)), 1:N)
        yc = map(x-> ntuple(i-> rand(Float32), Val(3)), 1:N)

        x = AT(xc)
        y = AT(yc)

        res1c = fill(0f0, N)
        res2c = similar(xc)

        res1 = AT(res1c)
        res2 = AT(res2c)

        res1 .= testv3_1.(x, y)
        res1c .= testv3_1.(xc, yc)
        @test Array(res1) ≈ res1c

        res2 .= testv3_2.(x, y)
        res2c .= testv3_2.(xc, yc)
        @test all(map((a,b)-> all((1,2,3) .≈ (1,2,3)), Array(res2), res2c))
    end
end


text:
"""
    update(x, measurement, H)

UPDATE step in Kalman filtering for linear dynamics models:
```math
K = P_{n+1}^P * H^T * S^{-1}
m_{n+1} = m_{n+1}^P + K * (0 - z)
P_{n+1} = P_{n+1}^P - K*S*K^T
```

This function provides a very simple UPDATE implementation.
In the solvers, we recommend to use the non-allocating [`update!`](@ref).
"""
function update(x::Gaussian, measurement::Gaussian, H::AbstractMatrix)
    m, C = x
    z, S = measurement

    K = C * H' * inv(S)
    m_new = m - K * z
    C_new = C - K * S * K'

    return Gaussian(m_new, C_new)
end
"""UPDATE step in Joseph-form, with square-root matrix inputs"""
function update(x::SRGaussian, measurement::SRGaussian, H::AbstractMatrix)
    m, C = x
    z, S = measurement

    K = C * H' * inv(S)
    m_new = m - K * z
    C_new = X_A_Xt(C, (I - K * H))

    return Gaussian(m_new, C_new)
end

"""
    update!(x_out, x_pred, measurement, H, R=0)

UPDATE step in Kalman filtering for linear dynamics models, given a measurement `Z=N(z, S)`.
In-place implementation of [`update`](@ref), saving the result in `x_out`.

```math
K = P_{n+1}^P * H^T * S^{-1}
m_{n+1} = m_{n+1}^P + K * (0 - z)
P_{n+1} = P_{n+1}^P - K*S*K^T
```

Implemented in Joseph Form.

See also: [`predict`](@ref)
"""
function update!(
    x_out::Gaussian,
    x_pred::Gaussian,
    measurement::Gaussian,
    H::AbstractMatrix,
    K_cache::AbstractMatrix,
    M_cache::AbstractMatrix,
    m_tmp,
)
    z, S = measurement.μ, copy!(m_tmp.Σ, measurement.Σ)
    m_p, P_p = x_pred.μ, x_pred.Σ
    D = length(m_p)

    # K = P_p * H' / S
    S_chol = cholesky!(S)
    K = _matmul!(K_cache, Matrix(P_p), H')
    rdiv!(K, S_chol)

    # x_out.μ .= m_p .+ K * (0 .- z)
    x_out.μ .= m_p .- _matmul!(x_out.μ, K, z)

    # M_cache .= I(D) .- mul!(M_cache, K, H)
    _matmul!(M_cache, K, H, -1.0, 0.0)
    @inbounds @simd ivdep for i in 1:D
        M_cache[i, i] += 1
    end

    X_A_Xt!(x_out.Σ, P_p, M_cache)

    return x_out
end


text:
# p14.jl - solve nonlinear BVP u_xx = exp(u), u(-1)=u(1)=0
#         (compare p13.jl)

N = 16;
(D,x) = cheb(N); D2 = D^2; D2 = D2[2:N,2:N];
u = zeros(N-1);
change = 1; it = 0;
while change > 1e-15                   # fixed-point iteration
    unew = D2\exp.(u);
    change = norm(unew-u,Inf);
    u = unew; it = it+1;
end
u = [0;u;0];
clf(); axes([.1,.4,.8,.5]);
plot(x,u,".",markersize=6);
xx = -1:.01:1;
uu = polyval(polyfit(x,u),xx);
plot(xx,uu), grid(true);
title("no. steps = $it      u(0) =$(u[Int(N/2+1)])")


text:
using AutomotivePOMDPs, POMDPs
using Base.Test
using ProgressMeter

function test_index(mdp::PedCarMDP, state_space::Vector{PedCarMDPState})
    @showprogress for (i, s) in enumerate(state_space)
        si = stateindex(mdp, s)
        if si != i
            println("bug for $i, index is $si")
            break
        end
    end
end

mdp = PedCarMDP(pos_res=3.0)

n_states(mdp)

@times state_space = states(mdp)

test_index(mdp, state_space)

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs