text:
#
# Tests depend on a running server at localhost:27017,
# and will create a database named "mongoc".
#
# Start a fresh MongoDB instance with:
#
# ```
# $ mkdir db
# $ mongod --dbpath ./db --smallfiles
# ```
#

import Mongoc

using Test
using Dates

const DB_NAME = "mongoc"

@testset "MongoDB" begin

    client = Mongoc.Client()

    @testset "Types" begin
        bson = Mongoc.BSON()
        @test_throws ErrorException Mongoc.Client("////invalid-url")
        @test client.uri == "mongodb://localhost:27017"
        Mongoc.set_appname!(client, "Runtests")
        db = client[DB_NAME]
        coll = db["new_collection"]

        io = IOBuffer()
        show(io, bson)
        show(io, client)
        show(io, db)
        show(io, coll)
        show(io, Mongoc.BSONCode("function() = 1"))
        show(io, Mongoc.QUERY_FLAG_TAILABLE_CURSOR)
    end

    @testset "Connection" begin

        @testset "ping" begin
            bson_ping_result = Mongoc.ping(client)
            @test haskey(bson_ping_result, "ok")
            @test Mongoc.as_json(Mongoc.ping(client)) == """{ "ok" : 1.0 }"""
        end

        @testset "Server Status" begin
            bson_server_status = Mongoc.command_simple(client["admin"], Mongoc.BSON("""{ "serverStatus" : 1 }"""))
            println("Server Mongo Version: ", bson_server_status["version"])
        end

        @testset "error print" begin
            error_happened = false
            try
                Mongoc.command_simple(client["hey"], Mongoc.BSON("""{ "you" : 1 }"""))
            catch e
                println(IOBuffer(), e)
                error_happened = true
            end

            @test error_happened
        end

        @testset "new_collection" begin
            coll = client[DB_NAME]["new_collection"]
            result = push!(coll, Mongoc.BSON("""{ "hello" : "world" }"""))
            @test Mongoc.as_json(result.reply) == """{ "insertedCount" : 1 }"""
            result = push!(coll, Mongoc.BSON("""{ "hey" : "you" }"""))
            @test Mongoc.as_json(result.reply) == """{ "insertedCount" : 1 }"""

            bson = Mongoc.BSON()
            bson["hey"] = "you"

            bson["zero_date"] = DateTime(0)
            bson["date_2018"] = DateTime(2018)

            result = push!(coll, bson)
            @test Mongoc.as_json(result.reply) == """{ "insertedCount" : 1 }"""

            i = 0
            for bson in coll
                @test haskey(bson, "hello") || haskey(bson, "hey")
                i += 1
            end
            @test i == length(coll)

            i = 0
            for bson in Mongoc.find(coll, Mongoc.BSON("""{ "hello" : "world" }"""))
                i += 1
            end
            @test i == 1

            Mongoc.command_simple(coll, Mongoc.BSON("""{ "collStats" : "new_collection" }"""))
            empty!(coll)
        end

        @testset "find_databases" begin
            found = false
            for obj in Mongoc.find_databases(client)
                if obj["name"] == DB_NAME
                    found = true
                end
            end
            @test found

            @test DB_NAME ∈ Mongoc.get_database_names(client)
        end

        @testset "catch cursor error" begin
            # issue #15
            invalid_client = Mongoc.Client("mongodb://invalid_url")
            collection = invalid_client["db_name"]["collection_name"]
            @test_throws ErrorException Mongoc.find_one(collection, Mongoc.BSON(""" { "a" : 1 } """))
        end

        @testset "Binary data" begin
            coll = client[DB_NAME]["new_collection"]
            bsonDoc = Mongoc.BSON()
            testdata = rand(UInt8, 100)
            bsonDoc["someId"] = "1234"
            bsonDoc["bindata"] = testdata
            result = push!(coll, bsonDoc)
            @test Mongoc.as_json(result.reply) == """{ "insertedCount" : 1 }"""

            # read data out and confirm
            selector = Mongoc.BSON("""{ "someId": "1234" }""")
            results = Mongoc.find_one(coll,  selector)

            @test results["bindata"] == testdata
        end

        @testset "find_collections" begin
            for obj in Mongoc.find_collections(client["local"])
                @test obj["name"] == "startup_log"
            end

            @test Mongoc.get_collection_names(client["local"]) == [ "startup_log" ]
        end

        @testset "BulkOperation" begin
            coll = client[DB_NAME]["new_collection"]
            bulk_operation = Mongoc.BulkOperation(coll)
            Mongoc.destroy!(bulk_operation)
            bulk_2 = Mongoc.BulkOperation(coll) # will be freed by GC
        end

        @testset "insert_many" begin
            collection = client[DB_NAME]["insert_many"]
            vector = Vector{Mongoc.BSON}()
            push!(vector, Mongoc.BSON("""{ "hey" : "you" }"""))
            push!(vector, Mongoc.BSON("""{ "out" : "there" }"""))
            push!(vector, Mongoc.BSON("""{ "count" : 3 }"""))

            append!(collection, vector)
            @test length(collection) == 3
            @test length(collect(collection)) == 3

            empty!(collection)
            @test isempty(collection)
        end

        @testset "delete_one" begin
            collection = client[DB_NAME]["delete_one"]
            doc = Mongoc.BSON("""{ "to" : "delete", "hey" : "you" }""")
            doc2 = Mongoc.BSON("""{ "to" : "keep", "out" : "there" }""")
            insert_result = push!(collection, doc)
            oid = Mongoc.BSONObjectId(insert_result.inserted_oid)
            push!(collection, doc2)

            selector = Mongoc.BSON()
            selector["_id"] = oid
            @test length(collection, selector) == 1
            result = Mongoc.delete_one(collection, selector)
            @test result["deletedCount"] == 1
            @test length(collection, selector) == 0

            empty!(collection)
        end

        @testset "delete_many" begin
            collection = client[DB_NAME]["delete_many"]
            append!(collection, [ Mongoc.BSON("""{ "first" : 1, "delete" : true }"""), Mongoc.BSON("""{ "second" : 2, "delete" : true }"""), Mongoc.BSON("""{ "third" : 3, "delete" : false }""") ])
            @test length(collection) == 3
            result = Mongoc.delete_many(collection, Mongoc.BSON("""{ "delete" : true }"""))
            @test result["deletedCount"] == 2
            @test length(collection) == 1
            result = Mongoc.delete_many(collection, Mongoc.BSON())
            @test result["deletedCount"] == 1
            @test isempty(collection)
        end

        @testset "update_one, update_many" begin
            collection = client[DB_NAME]["update_one"]
            append!(collection, [ Mongoc.BSON("""{ "first" : 1, "delete" : true }"""), Mongoc.BSON("""{ "second" : 2, "delete" : true }"""), Mongoc.BSON("""{ "third" : 3, "delete" : false }""") ])
            @test length(collection) == 3

            selector = Mongoc.BSON("""{ "delete" : false }""")
            update = Mongoc.BSON("""{ "\$set" : { "delete" : true, "new_field" : 1 } }""")
            result = Mongoc.update_one(collection, selector, update)

            @test result["modifiedCount"] == 1
            @test result["matchedCount"] == 1
            @test result["upsertedCount"] == 0

            updated_bson = Mongoc.find_one(collection, Mongoc.BSON("""{ "third" : 3 }"""))
            @test updated_bson != nothing
            @test updated_bson["delete"] == true
            @test updated_bson["new_field"] == 1

            selector = Mongoc.BSON("""{ "delete" : true }""")
            update = Mongoc.BSON("""{ "\$set" : { "delete" : false } }""")
            result = Mongoc.update_many(collection, selector, update)
            @test result["modifiedCount"] == 3
            @test result["matchedCount"] == 3
            @test result["upsertedCount"] == 0

            for doc in Mongoc.find(collection)
                @test doc["delete"] == false
            end

            @test Mongoc.find_one(collection, Mongoc.BSON("""{ "delete" : true }""")) == nothing

            empty!(collection)
        end

        @testset "aggregation, map_reduce" begin
            # reproducing the examples at https://docs.mongodb.com/manual/aggregation/
            docs = [
                Mongoc.BSON("""{ "cust_id" : "A123", "amount" : 500, "status" : "A" }"""),
                Mongoc.BSON("""{ "cust_id" : "A123", "amount" : 250, "status" : "A" }"""),
                Mongoc.BSON("""{ "cust_id" : "B212", "amount" : 200, "status" : "A" }"""),
                Mongoc.BSON("""{ "cust_id" : "A123", "amount" : 300, "status" : "D" }""")
            ]

            database = client[DB_NAME]
            collection = database["aggregation_example"]
            append!(collection, docs)
            @test length(collection) == 4

            # Aggregation
            let
                bson_pipeline = Mongoc.BSON("""
                    [
                        { "\$match" : { "status" : "A" } },
                        { "\$group" : { "_id" : "\$cust_id", "total" : { "\$sum" : "\$amount" } } }
                    ]""")

                # Response should be
                #   BSON("{ "_id" : "B212", "total" : 200 }")
                #   BSON("{ "_id" : "A123", "total" : 750 }")
                for doc in Mongoc.aggregate(collection, bson_pipeline)
                    if doc["_id"] == "A123"
                        @test doc["total"] == 750
                    elseif doc["_id"] == "B212"
                        @test doc["total"] == 200
                    else
                        # shouldn't get in here
                        @test false
                    end
                end
            end

            # map_reduce
            let
                input_collection_name = "aggregation_example"
                mapper = Mongoc.BSONCode(""" function() { emit( this.cust_id, this.amount ); } """)
                reducer = Mongoc.BSONCode(""" function(key, values) { return Array.sum( values ) } """)
                output_collection_name = "order_totals"
                query = Mongoc.BSON("""{ "status" : "A" }""")

                map_reduce_command = Mongoc.BSON()
                map_reduce_command["mapReduce"] = input_collection_name
                map_reduce_command["map"] = mapper
                map_reduce_command["reduce"] = reducer
                map_reduce_command["out"] = output_collection_name
                map_reduce_command["query"] = query

                result = Mongoc.command_simple(database, map_reduce_command)
                @test result["result"] == "order_totals"
                @test result["ok"] == 1.0

                for doc in Mongoc.find(database["order_totals"])
                   if doc["_id"] == "A123"
                        @test doc["value"] == 750
                    elseif doc["_id"] == "B212"
                        @test doc["value"] == 200
                    else
                        # shouldn't get in here
                        @test false
                    end
                end

               # Response should be
               # BSON("{ "_id" : "A123", "value" : 750.0 }")
               # BSON("{ "_id" : "B212", "value" : 200.0 }")
            end

            empty!(collection)
        end
    end

    @testset "Users" begin
        # creates admin user - https://docs.mongodb.com/manual/tutorial/enable-authentication/
        @test Mongoc.has_database(client, DB_NAME) # at this point, DB_NAME should exist
        database = client[DB_NAME]

        user_name = "myUser"
        pass = "abc123"
        roles = Mongoc.BSON()

        if Mongoc.has_user(database, user_name)
            Mongoc.remove_user(database, user_name)
        end

        Mongoc.add_user(database, user_name, pass, roles)
        Mongoc.remove_user(database, user_name)
        @test !Mongoc.has_user(database, user_name)
    end

    @testset "Session Options" begin
        opt = Mongoc.SessionOptions()
        @test Mongoc.get_casual_consistency(opt)
        Mongoc.set_casual_consistency!(opt, false)
        @test !Mongoc.get_casual_consistency(opt)
    end

    server_version = Mongoc.get_server_mongodb_version(client)

    if server_version < v"3.6"
        @warn("MongoDB server version $server_version does not support Sessions. Skipping tests.")
    else
        @testset "Session" begin
            session = Mongoc.Session(client)
            db = session[DB_NAME]
            collection = db["session_collection"]
            push!(collection, Mongoc.BSON("""{ "try-insert" : 1 }"""))
            empty!(collection)
        end
    end
end


text:
using Distributed
using Distributions
using PyPlot

function playWithData(;subj_id = 3
                        #, i_episode = 20
                        )

    data_path = "./mcmc_rl_fit/projects/fmri/data/SARSPEIC_all_fMRI.csv"
    max_nr_episodes = 50

    all_data, all_data_stats = RL_Fit.load_SARSPEI(data_path; max_episodes_per_subj = max_nr_episodes)

    subj_data = RL_Fit.filter_by_person(all_data, subj_id)

    all_puzzles = RL_Fit.get_all_puzzles(subj_data)
    i_puzzle = 1

    all_epi_in_puzzle = RL_Fit.get_all_episodes_in_puzzle(subj_data, i_puzzle)

    episode_data = []
    for i_episode in all_epi_in_puzzle
        episode_data = RL_Fit.get_data_by_puzzle_episode(subj_data, i_puzzle, i_episode)
        for i in 1:size(episode_data, 1)
            @show episode_data[i, :]
        end
        println("------")
    end
    # episode_data
end

# NOTE:
# ----
# Row coding of data file:
#SARSPE_COLS
#c_state=1
#c_action=2
#c_reward=3
#c_nextState=4
#c_puzzle=5
#c_episode=6
#c_subject=7
#c_surprise_trial=8
#c_RT_prctile=9
#c_action_value=10
#c_graph_version=11

function get_nr_episodes()
    data_path = "./mcmc_rl_fit/projects/fmri/data/SARSPEICZVG_all_fMRI.csv"
    max_episodes_per_subj = 70

    all_data, all_stats = RL_Fit.load_SARSPEI(data_path;max_episodes_per_subj=max_episodes_per_subj)
    all_subject_ids = all_stats.all_subject_ids

    for subj_id in all_subject_ids
        println("--------")
        @show subj_id

        subj_data = RL_Fit.filter_by_person(all_data, subj_id)
        graph_id = subj_data[1,Int(RL_Fit.c_graph_version)]
        @show graph_id
        nr_states = size(subj_data, 1)
        @show nr_states
        nr_surprise_trials = sum(subj_data[:, Int(RL_Fit.c_surprise_trial)])
        @show nr_surprise_trials
        @show nr_surprise_trials/nr_states
        all_episodes = unique(subj_data[:,Int(RL_Fit.c_episode)])
        @show maximum(all_episodes)
    end
end

function get_someinfonow()
    # Get some info about all subjects

    data_path = "./mcmc_rl_fit/projects/fmri/data/SARSPEICZVG_all_fMRI.csv"
    max_episodes_per_subj = 80
    (all_data, all_stats) = RL_Fit.load_SARSPEI(data_path;max_episodes_per_subj=max_episodes_per_subj)
    all_subject_ids = all_stats.all_subject_ids
    # println("all_subject_ids: $all_subject_ids")

    graph_ids = zeros(Int64, length(all_subject_ids))
    nr_episodes = zeros(Int64, length(all_subject_ids))
    nr_actions = zeros(Int64, length(all_subject_ids))

    for subj_id in all_subject_ids
        @show subj_id
        subj_data = RL_Fit.filter_by_person(all_data, subj_id)

        graph_id = subj_data[1,Int(RL_Fit.c_graph_version)]
        @show graph_id
        graph_ids[subj_id] = copy(graph_id)

        all_episodes = unique(subj_data[:,Int(RL_Fit.c_episode)])
        @show maximum(all_episodes)
        nr_episodes[subj_id] = copy(maximum(all_episodes))

        actions = subj_data[:,Int(RL_Fit.c_action)]
        @show length(actions)
        nr_actions[subj_id] = copy(length(actions))

        println("-------------------")
    end
    graph_ids, nr_episodes, nr_actions
end

function singleSubjectData(;id_subject=3)
    # Using an algorithm (edit: learner_function_V below) and a set of parameters
    # (edit params below) go through the actions of one subject (id_subject)
    # and obtain the log-likelihood and the corresponding learning signals
    # (eg RPE, SPE...)
    data_path = "./mcmc_rl_fit/projects/fmri/data/SARSPEIC_all_fMRI.csv"

    data_single_subj, data_stats = RL_Fit.load_SARSPEI(data_path; subject_id = id_subject)
    # -------
    # Surprise Actor-critic
    # ---------------------
    learner_function_V = RL_Fit.getLL_ActorCritic_TDLambda_modulatedActor_V
    params = [0.79724,  #alpha
    0.00433, #surpriseprobability
    2.78163, #stochasticity
    3.24411e-5, #eta
    0.99664, #gamma
    0.79392,  #lambda_actor
    0.30497,  #lambda_critic
    0.39196, # temp
    15, #nparticles
    ]
    # -------
    # SARSA-λ
    # -------
    # learner_function_V = RL_Fit.getLL_SarsaLambda_V
    # params = [0.2, #alpha
    # 0.9999, #gamma
    # 0.61, #lambda
    # 0.25,# temp
    # ]
    # -------
    # REINFORCE
    # ---------
    # learner_function_V = RL_Fit.getLL_REINFORCE_V
    # params = [0.96, #alpha
    # 0.76, #gamma
    # 0.95,# temp
    # ]
    # # -------
    # Actor-critic
    # ------------
    # learner_function_V = RL_Fit.getLL_ActorCritic_TDLambda_V
    # params = [0.5009131276460307, #alpha
    # 4.6834146251373146e-5, #eta
    # 0.8205988691213357, #gamma
    # 0.9599330699486914,#lambda actor
    # 0.5782935181795414,#lambda critic
    # 0.24363820949361994# temp
    # ]
    # -------
    # -------
    println("hi")
    @show learner_function_V
    LL, signals = learner_function_V(data_single_subj, params)
    signals=signals[id_subject]

    LL, signals
end

function lookatpolicyimbalances(signalsArray;
                                labels = ["REINFORCE", "Actor-critic", "SARSA"],
                                temperature_index = [3, 6, 4],
                                # save_path = "./"
                                subjID = 3
                                )
    # Plot policy imbalances (similar to Fig. A.3)
    # signalsArray = array with signals of different algorithms
    #                run on the data of one subject (subjID).
    #                (ie output of some learner_function_V,
    #                for example RL_Fit.getLL_ActorCritic_TDLambda_V)
    # temperature_index = index of temperature parameter in the list of parameters
    #                   of the algorithm (look at order in eg RL_Fit.getLL_ActorCritic_TDLambda_V)

    save_path = makehomesavepath("test_policyImbalance_")
    mkdir(save_path)

    fig, ax = subplots(nrows=3, ncols=1, figsize=(13, 19))
    ax_idx = 0

    # ---- f(st,at_chosen) - f(st,at_notchosen)
    ax_idx += 1
    for i in 1:length(signalsArray)
        ax[ax_idx].scatter(1:length(signalsArray[i].results_dict["PolicyParametersDiffChosenVsUnchosen"]),
                    signalsArray[i].results_dict["PolicyParametersDiffChosenVsUnchosen"],
                    label=labels[i])
    end
    ax[ax_idx].set_title(L"$f(s_t, a_{chosen}) - f(s_t, a_{not chosen})$")
    ax[ax_idx].set_xlabel("actions")
    ax[ax_idx].grid()
    ax[ax_idx].legend()

    # ---- f(st,at_chosen) - f(st,at_notchosen) / τ
    ax_idx += 1
    for i in 1:length(signalsArray)
        temperature = signalsArray[i].param_vector[temperature_index[i]]
        # @show temperature
        ax[ax_idx].scatter(1:length(signalsArray[i].results_dict["PolicyParametersDiffChosenVsUnchosen"]),
                    signalsArray[i].results_dict["PolicyParametersDiffChosenVsUnchosen"] ./ temperature,
                    label=labels[i])
    end
    ax[ax_idx].set_title(L"$(f(s_t, a_{chosen}) - f(s_t, a_{not chosen})) / \tau $")
    ax[ax_idx].set_xlabel("actions")
    ax[ax_idx].grid()
    ax[ax_idx].legend()

    # ---- π(st,at_chosen) - π(st,at_notchosen)
    # fig = figure(); ax = gca()
    ax_idx += 1
    for i in 1:length(signalsArray)
        ax[ax_idx].scatter(1:length(signalsArray[i].results_dict["PolicyDiffChosenVsUnchosen"]),
                    signalsArray[i].results_dict["PolicyDiffChosenVsUnchosen"],
                    label=labels[i])
    end
    ax[ax_idx].set_title(L"$\pi(s_t, a_{chosen}) - \pi(s_t, a_{not chosen})$")
    ax[ax_idx].set_xlabel("actions")
    ax[ax_idx].grid()
    ax[ax_idx].legend()

    fig.tight_layout()
    fig.savefig(joinpath(save_path, "subj_"*string(subjID)*"PolImbalances_.pdf"))
    plt.close()
end

function plotscattersignals(signalsArray;
                                labels = ["REINFORCE",
                                        "Actor-critic",
                                        "SARSA"],
                                save_path = "./",
                                subjID = 3
                                )
    # Plot some signals.
    # signalsArray = array with signals of different algorithms
    # run on the data of one subject (subjID).
    # (ie output of some learner_function_V)
    # NOTE: this function is just for illustration and playing around!
    # Need to edit which signals you would like to plot.

    # ---- SPE
    fig = figure(); ax = gca()
    for i in 1:length(signalsArray)
        ax.scatter(1:length(signalsArray[i].SPE),
                    signalsArray[i].SPE,
                    label=labels[i])
    end
    ax.set_title("SPE")
    ax.set_xlabel("actions")
    ax.grid()
    ax.legend()

    # --- Surprise
    fig2 = figure(); ax2 = gca()
    for i in 1:length(signalsArray)
        if  signalsArray[i].algorithm_name == "hybrid_actorcritic_FWD"
            ax2.scatter(1:length(signalsArray[i].SPE),
                        signalsArray[i].SPE,
                        label=labels[i])
        else
            ax2.scatter(1:length(signalsArray[i].results_dict["Surprise"]),
                        signalsArray[i].results_dict["Surprise"],
                        label=labels[i])
        end
    end
    ax2.set_title(L"$S_{BF}$")
    ax2.set_xlabel("actions")
    ax2.grid()
    ax2.legend()

end

ext: jl


text:
import Base.merge,Base.==

export merge, ==

"""
    merge(cons::Array{SeisCon,1})

Merge the elements of `con`, a vector of SeisCon objects, into one SeisCon object.
"""
function merge(cons::Array{SeisCon,1})

    # Check similar metadata
    ns = get_confield(cons, :ns)
    dsf = get_confield(cons, :dsf)

    if all(ns.==ns[1]) && all(dsf.==dsf[1])
        d = [cons[i].blocks for i in 1:length(cons)]
        return SeisCon(ns[1], dsf[1], vcat(d...), cons[1].fileheader)
    else
        @error "Dissimilar metadata, cannot merge"
    end
end

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

Merge two SeisCon objects together
"""
merge(a::SeisCon, b::SeisCon) = merge([a; b])

"""
    merge(a::SeisBlock, b::SeisBlock; kwargs...)

Merge two SeisBlocks into one.
"""
function merge(a::SeisBlock{DT}, b::SeisBlock{DT};
               force::Bool = false, consume::Bool = false) where {DT<:Union{IBMFloat32, Float32}}
   merge([a; b], force = force, consume = consume)
end

"""
    merge(a::SeisBlock{DT}, b::SeisBlock{DT};
               force::Bool = false, consume::Bool = false) where {DT<:Union{IBMFloat32, Float32}}

Merge `blocks`, a vector of `SeisBlock` objects into one `SeisBlock` object.

By default, `merge` checks to ensure that each `SeisBlock` has matching fileheaders. This ensures
that the returned `SeisBlock` is consistant with it's fileheader. This check can be turned off
using the `force` keyword. If set to false, `merge` will attempt to combine `blocks` without
checking fileheader metadata. The first fileheader in `blocks` will be used for the  returned
`SeisBlock`'s fileheader.

By default, `merge` references traceheaders, but copies data from the elements of `blocks`.
If memory use is a concern, the `consume` keyword will clear the traceheader and data field of each
element of `block` after it has been copied. This will prevent duplicating the data in memory,
at the cost of clearing the elements of `blocks`.

# Examples

    julia> s = segy_scan(joinpath(dirname(pathof(SegyIO)),"data/"), "overthrust", ["GroupX"; "GroupY"], verbosity = 0);

    julia> a = s[1:4]; b = s[5:8];

    julia> c = merge([a; b]);

    julia> c.traceheaders[1] === a.traceheaders[1]
    true

    julia> c = merge([a; b], consume = true);

    julia> a.traceheaders
    0-element Array{SegyIO.BinaryTraceHeader,1}

"""
function merge(blocks::Vector{SeisBlock{DT}};
               force::Bool = false, consume::Bool = false) where {DT<:Union{IBMFloat32, Float32}}

    # Check common file header
    n = length(blocks)
    if !force && !all([blocks[i].fileheader == blocks[i+1].fileheader for i in 1:n-1])
        @error "Fileheaders do not match for all blocks, use kw force=true to ignore"
    end

    fh = blocks[1].fileheader
    bth = Vector{BinaryTraceHeader}()
    data = Vector{DT}()
    for ii in 1:n
        if consume
            append!(bth, blocks[ii].traceheaders[:])
            append!(data, blocks[ii].data[:])
            (blocks[ii].data = Array{DT}(undef,0,0))
            (blocks[ii].traceheaders = Array{BinaryTraceHeader}(undef,0))
        else
            append!(bth, @view blocks[ii].traceheaders[:])
            append!(data, blocks[ii].data[:])
        end
    end

    out = SeisBlock(fh, bth, reshape(data, Int64(fh.bfh.ns), :))
end

function ==(a::FileHeader, b::FileHeader)
    match_th = a.th == b.th
    match_bfh = [getfield(a.bfh, field) == getfield(b.bfh,field) for field in fieldnames(BinaryFileHeader)]
    return match_th && all(match_bfh)
end


text:
using Distributions
using Random: shuffle
import Base: iterate

export select,
    SelectionOperator,
    SelectionResult,
    setup!

export FitnessProportionateSelection,
    L1Selection,
    PairWithBestSelection,
    SoftmaxSelection,
    TournamentSelection

abstract type SelectionOperator{G, P, K} end

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

"""
    selection(population, operator[, generation]) -> selection

Select parts of population of a population to be used in breeding.  Should compare fitnesses using
`isless`, if that is relevant.

The result is an iterable of `Family{G, K}`.
"""
selection(population::Population{G}, operator::SelectionOperator{G, P, K},
          generation::Integer) where {G, P, K} =
    selection(population, operator)

mutable struct PairWithBestSelection{G, P} <: SelectionOperator{G, P, 2}
    model::AbstractPopulationModel

    PairWithBestSelection{G, P}() where {G, P} = new{G, P}()
end

function setup!(operator::PairWithBestSelection, model::AbstractPopulationModel)
    operator.model = model
    operator
end

function selection(population::Population{G}, operator::PairWithBestSelection{G, P}) where {G, P}
    # simple naive groupings that pair the best entitiy with every other
    PairWithBestSelectionIterator{P}(population, operator.model.fittest)
end

struct PairWithBestSelectionIterator{P, M, G}
    population::Population{G}
    fittest::Individual{G}

    function PairWithBestSelectionIterator{P}(population::Population{G},
                                              fittest::Individual{G}) where {P, G}
        M = length(population) ÷ P
        new{P, M, G}(population, fittest)
    end
end

function iterate(itr::PairWithBestSelectionIterator{P, M}, state = 0) where {P, M}
    if state ≥ M
        return nothing
    else
        (itr.fittest, rand(itr.population)), state + 1
    end
end

# function iterate(itr::PairWithBestSelectionIterator{1})
#     next = iterate(itr.population)
#     if next == nothing
#         return nothing
#     else
#         individual, state = next
#         (itr.fittest, individual), state
#     end
# end

# function iterate(itr::PairWithBestSelectionIterator{1}, state)
#     next = iterate(itr.population, state)
#     if next == nothing
#         return nothing
#     else
#         individual, state = next
#         (itr.fittest, individual), state
#     end
# end

struct FitnessProportionateSelection{G, P, K, F} <: SelectionOperator{G, P, K}
    transform::F
    temperature::Rate
end

function selection(population::Population{G},
                   operator::FitnessProportionateSelection{G, P, K},
                   generation::Integer) where {G, P, K}
    function transform(f)
        operator.transform(f, operator.temperature(generation))
    end
    FitnessProportionateSelectionIterator{P, K}(population, transform)
end

struct FitnessProportionateSelectionIterator{M, P, K, G}
    population::Population{G}
    transform

    function FitnessProportionateSelectionIterator{P, K}(population::Population{G},
                                                         transform) where {P, K, G}
        M = length(population) ÷ P
        new{M, P, K, G}(population, transform)
    end
end

function iterate(itr::FitnessProportionateSelectionIterator{M, P, K}, state = 0) where {M, P, K}
    if state ≥ M
        return nothing
    else
        dist = Categorical(itr.transform(fitness.(itr.population)))
        indices = eachindex(itr.population)
        ntuple(i -> itr.population[rand(dist)], Val{K}()), state + 1
    end
end

function l1normalize(f, _)
    y = f .- minimum(f)
    s = sum(y)
    s == 0 ? fill(1/length(f), size(f)) : y ./ s
end

function softmax(f, θ = 1)
    y = float(f) ./ θ
    y .= exp.(y .- maximum(y))
    y ./ sum(y)
end

const SoftmaxSelection{G, P, K} = FitnessProportionateSelection{G, P, K, typeof(softmax)}
(::Type{SoftmaxSelection{G, P, K}})(rate::Rate = ConstantRate(1.0)) where {G, P, K} =
    SoftmaxSelection{G, P, K}(softmax, rate)

const L1Selection{G, P, K} = FitnessProportionateSelection{G, P, K, typeof(l1normalize)}
(::Type{L1Selection{G, P, K}})() where {G, P, K} =
    L1Selection{G, P, K}(l1normalize, ConstantRate(1.0))

struct TournamentSelection{G, S, P, K} <: SelectionOperator{G, P, K} end

selection(population::Population{G}, operator::TournamentSelection{G, S, P, K}) where {G, S, P, K} =
    TournamentSelectionIterator{S, P, K}(population)

struct TournamentSelectionIterator{M, S, P, K, G}
    population::Population{G}

    function TournamentSelectionIterator{S, P, K}(population::Population{G}) where {S, P, K, G}
        M = length(population) ÷ P
        new{M, S, P, K, G}(population)
    end
end

function iterate(itr::TournamentSelectionIterator{M, S, P, K}, state = 0) where {M, S, P, K}
    if state ≥ M
        return nothing
    else
        ntuple(i -> maximumby(fitness, randview(itr.population, S)), Val{K}()), state + 1
    end
end

struct TruncationSelection{G, P} <: SelectionOperator{G, P, 1} end

function selection(population::Population{G}, operator::TruncationSelection{G, P}) where {G, P}
    μ = length(population) ÷ P
    partialsortperm(model.population, 1:μ, by = fitness, rev = true)
end

# struct SUSSelection{G, P, K} <: SelectionOperator{G, P, K} end

# function selection(population::Population{G}, operator::SUSSelection{S, P, K})
#     fitness_cdf = fitness.(population)
#     fitness_cdf .-= minimum(fitness_cdf)
#     sum(fitness_cdf) == 0 && fitness_cdf .= 1
#     cumsum!(fitness_cdf, fitness_cdf)
#     range = fitness_cdf[end] / length(fitness_cdf)

#     SUSSelectionIterator{P, K}(population, shuffle(eachindex(population)), fitness_cdf, range)
# end

# struct SUSSelectionIterator{M, P, K, G}
#     population::Population{G}
#     permutation::Vector{Int}
#     fitness_cdf::Vector{Float64}
#     range::Float64

#     function SUSSelectionIterator{P, K}(population::Population{G}, permutation::Vector{Int}
#                                         fitness_cdf::Vector{Float64}, position::Float64) where {P, K, G}
#         M = length(population) ÷ P
#         new{M, P, K, G}(population, permutation, fitness_cdf, position)
#     end
# end

# function iterate(itr::SUSSelectionIterator{M, P, K}, state = rand() * itr.range) where {P, K}
#     if state > f[end]
#         return nothing
#     else
#         index = findfirst(f -> f ≥ itr.position, view(itr.fitness_cdf, itr.permutation))
#         return view(itr.population, itr.permutation)[index], state + itr.range
#     end
# end


text:
using Documenter, DependenciesParser

makedocs(modules = [DependenciesParser], sitename = "DependenciesParser.jl")
deploydocs(repo = "github.com/Nosferican/DependenciesParser.jl.git")


text:
"""-----------------------------------------------------------------------------
  This function runs the HOM algorithm [1] on an instance of a synthetic
  experiments in Park et al.[2] from the framework introduced in [3]. The
  routine takes in a number of points n, and builds source and target points
  using

                         t_i = scale*(s_i + σ*randn(d)),

  where s\\_i is the ith, d=2 dimensional random normal generated. The target
  points are then randomly shuffled and the HOM procedure is run.

  Inputs
  ------
  * n - (Int):
    The number of pairs of points to generate in the source and target.
  * sigma - (Float):
    The standard deviation of the random noise to perturb the points by.
  * outliers - (Int):
    Additional random normal points to add into the correspondence problem.
    Generated points are normal points centered at the origin with
    covariance = I. Default is 0.
  * scale - (Float):
    Desired scaling to the correspondence problem.
  * seed - (Any):
    Any desired seed, other experiment drivers use UInts.

  Ouputs
  ------
  * Graduated_Assignment_res - (Array{Float64,2}):
    The converged iterate X from the graduated assignment procedure. Run
    build_assignment to get the final matching.
  * runtime - (Float64):
    The runtime of the graduated assignment procedure run on the generated
    tensor.
  * p - (Array{Float,1}):
    The ground truth permutation used to shuffle the target points.

  Citation
  --------
  [1] - M. Chertok and Y. Keller, "Efficient high order matching," IEEE
        Transactions on Pattern Analysis and Machine Intelligence, vol. 32
        no. 12, pp. 2205-2215, 2010.
  [2] - S. Park, S.-K. Park, and M. Herbert, "Fast and scalable approximate
        spectral matching for higher order grpah matching," IEEE transactions
        on pattern analysis and machine intelligence, vol. 36, no. 3,
        pp. 479-492,2013.
  [3] - R. Zass and A. Shashua, "Probabalistic Graph and Hypergraph Matching"
        in 2008 IEEE Conference on Computer Vision and Pattern Recognition.
        IEEE, 2008, pp. 1-8.
-----------------------------------------------------------------------------"""
function synthetic_HOM(n::Int,sigma::Float64,outliers::Int=0,scale::Float64=1.0;seed=nothing)

    if seed !== nothing
        Random.seed!(seed)
    end

    d = 2
    source_points = randn(Float64,n+outliers,d)
    target_points = Array{Float64,2}(undef,n+outliers,d)
    for i::Int =1:n
        target_points[i,:] = source_points[i,:] + randn(d)*sigma
        target_points[i,:] *= scale
    end

    for i::Int=1:outliers
        target_points[n+i,:] = randn(d)
        source_points[n+i,:] = randn(d)
    end

    n += outliers
    m = n

    #shuffle the target points
    p = shuffle(1:n)
    target_points = target_points[p,:]

    #find all the triangles in between the points
    source_triangles = brute_force_triangles(source_points)
    target_triangles = brute_force_triangles(target_points)

    indices, vals = produce_HOM_tensor(source_triangles,target_triangles,n,tol=1e-12)

    marginalized_tensor = sym_mode1_marginalization(indices,vals,n^2)
    #marg_result = reshape(HOM_graduated_assignment(marginalized_tensor),n,n)
    #lam, vec, _ = eigs(marginalized_tensor,nev=1)

    Graduated_Assignment_res, runtime = @timed Array(reshape(HOM_graduated_assignment(marginalized_tensor,m),n,n)')

    return Graduated_Assignment_res, runtime , p
    #ten_result = reshape(HOM_graduated_assignment(ssten.COOTen(indices,vals)),n,n)
#    return marginalized_tensor,reshape(vec,n,n), marg_result, ten_result, p

end

"""-------------------------------------------------------------------------
     produces the angular difference between all pairs of triangles that
   meet the sparsification tolerance. Final output tensor is a symmetric
   tensor which computes the elements with

     C_{ii',jj',kk'} =
          exp((|θ_i - θ_i'|^2 + |θ_j - θ_j'|^2 + |θ_k - θ_k'|^2)/ϵ^2)

    Final output is a symmetric COO formatted tensor.

    Inputs:
    -------
    * image[1-2]_triangles - (list of (3-tuples,3-tuples) pairs):
      The set of triangles in each image, keys are the indices sorted in
      increasing order of the tensors linked to the corresponding angles in
      their triangle.
    * m - (Int):
      The number of points in image2, needed to format the final indices.
    * tol - (Float64)
      The sparsification tolerance as to whether or not to include that pair
      of triangles in the final tensor.

    Output:
    -------
     * HOM_indices - (Array{Int,2}):
       All the indices of the output tensor up to permutations (entries in
       sorted order).
     * HOM_vals  - (Array{Float64,1}):
       All the values of the corresponding hyper edges.
-------------------------------------------------------------------------"""
function produce_HOM_tensor(image1_triangles::Array{Tuple{Tuple{Int,Int,Int},Tuple{F,F,F}},1},
                            image2_triangles::Array{Tuple{Tuple{Int,Int,Int},Tuple{F,F,F}},1},
                            m::Int; tol::Float64=1e-5,test_mode=false) where {F <: AbstractFloat}
    max_tris  = size(image1_triangles,1)*size(image2_triangles,1)*6

    HOM_indices = Array{Int64,2}(undef,max_tris,3)
    HOM_vals = Array{Float64,1}(undef,max_tris)

    hyperedge_index = 1
    for ((i1,j1,k1),img1_angles) in image1_triangles
        for (indices,img2_angles) in image2_triangles
            for ((i2,j2,k2),angle_perm) in zip(permutations(indices),permutations(img2_angles))

                val = HOM_angle_diff(img1_angles,tuple(angle_perm...),1.0,test_mode=test_mode)
                if val > tol
                    HOM_indices[hyperedge_index,1] = (i1 -1)*m + i2
                    HOM_indices[hyperedge_index,2] = (j1 -1)*m + j2
                    HOM_indices[hyperedge_index,3] = (k1 -1)*m + k2
                    HOM_vals[hyperedge_index] = val
                    hyperedge_index += 1
                end
            end
        end
    end

    hyperedge_index -= 1 #correct last index increment
    return HOM_indices[1:hyperedge_index,:], HOM_vals[1:hyperedge_index]

end

function HOM_angle_diff(angles1::Tuple{F,F,F}, angles2::Tuple{F,F,F}, epsilon::F;test_mode=false) where {F <: AbstractFloat}
    a1,a2,a3 = angles1
    ap1,ap2,ap3 = angles2
    if test_mode
        return ((a1 - ap1)^2+(a2 - ap2)^2+(a3 - ap3)^2) / epsilon^2
    else
        return exp(-((a1 - ap1)^2+(a2 - ap2)^2+(a3 - ap3)^2) / epsilon^2)
    end
end


text:
function with_temp_dir(f::Function)
    original_directory = pwd()
    tmp_dir = mktempdir()
    atexit(() -> rm(tmp_dir; force = true, recursive = true))
    cd(tmp_dir)
    result = f(tmp_dir)
    cd(original_directory)
    rm(tmp_dir; force = true, recursive = true)
    return result
end


text:
export getNodalLaplacianMatrix, multOpNeumann!, Lap2DStencil,dxxMat,BC
function BC()
	return 2.0; # one is for first order Neumann, and 2 is for 2nd order. Note that with 2, the matrix is not symmetric!!!
end

function dxxMat(n::Int64,h::Float64)
O1 = -ones(n-1);
O1[n-1] = -BC();
O2 = 2.0*ones(n);
O2[1] = BC();
O2[n] = BC();
O3 = -ones(n-1);
O3[1] = -BC();
dxx = spdiagm((O1/(h^2),O2/(h^2),O3/(h^2)),[-1,0,1],n,n);
return dxx;
end

function getNodalLaplacianMatrix(Msh::RegularMesh)
nodes = Msh.n+1;
I1 = speye(nodes[1]);
Dxx1 = dxxMat(nodes[1],Msh.h[1]);
I2 = speye(nodes[2]);
Dxx2 = dxxMat(nodes[2],Msh.h[2]);
if Msh.dim==2
	L = kron(I2,Dxx1) + kron(Dxx2,I1);
else
	I3 = speye(nodes[3]);
	Dxx3 = dxxMat(nodes[3],Msh.h[3]);
	L = kron(I3,kron(I2,Dxx1) + kron(Dxx2,I1)) + kron(Dxx3,kron(I2,I1));
end
return L;
end

##
## [0  x2  0]
## [x4 x1 x5]
## [0  x3 0]

function Lap2DStencil(x1::Float64,x2::Float64,x3::Float64,x4::Float64,x5::Float64,h1invsq::Float64,h2invsq::Float64)
	return (2*h1invsq + 2*h2invsq)*x1 - h1invsq*(x2 + x3) - h2invsq*(x4 + x5);
end

function multOpNeumann!(M::RegularMesh,x::Array,y::Array,op::Function)
n = M.n;

if M.dim==2
	n1 = n[1]+1;
	n2 = n[2]+1;
	h1invsq = 1./(M.h[1]^2);
	h2invsq = 1./(M.h[2]^2);
	@inbounds y[1] = op(x[1],x[1],x[2],x[1],x[n1+1],h1invsq,h2invsq);
	@simd for i=2:n1-1
		@inbounds y[i] = op(x[i],x[i-1],x[i+1],x[i],x[i+n1],h1invsq,h2invsq);
	end
	@inbounds y[n1] = op(x[n1],x[n1-1],x[n1],x[n1],x[n1 + n1],h1invsq,h2invsq);
	for j=2:n2-1
		colShift = n1*(j-1);
		i = 1 + colShift;
		@inbounds y[i] = op(x[i],x[i],x[i+1],x[i-n1],x[i+n1],h1invsq,h2invsq);
		@simd for i = (2 + colShift):(n1-1 + colShift)
			@inbounds y[i] = op(x[i],x[i-1],x[i+1],x[i-n1],x[i+n1],h1invsq,h2invsq);
		end
		i = n1 + colShift;
		@inbounds y[i] = op(x[i],x[i-1],x[i],x[i-n1],x[i+n1],h1invsq,h2invsq);
	end
	colShift = n1*(n2-1);
	i = 1 + colShift;
	@inbounds y[i] = op(x[i],x[i],x[i+1],x[i-n1],x[i],h1invsq,h2invsq);
	@simd for i = (2 + colShift):(n1-1 + colShift)
		@inbounds y[i] = op(x[i],x[i-1],x[i+1],x[i-n1],x[i],h1invsq,h2invsq);
	end
	i = n1 + colShift;
	@inbounds y[i] = op(x[i],x[i-1],x[i],x[i-n1],x[i],h1invsq,h2invsq);
else
end
end

function multOpDirichlet!(M::RegularMesh,x::Array,y::Array,op::Function)
n = M.n;

if M.dim==2
	n1 = n[1]+1;
	n2 = n[2]+1;
	h1invsq = 1./(M.h[1]^2);
	h2invsq = 1./(M.h[2]^2);
	@inbounds y[1] = op(x[1],0.0,x[2],0.0,x[n1+1],h1invsq,h2invsq);
	@simd for i=2:n1-1
		@inbounds y[i] = op(x[i],x[i-1],x[i+1],0.0,x[i+n1],h1invsq,h2invsq);
	end
	@inbounds y[n1] = op(x[n1],x[n1-1],0.0,0.0,x[n1 + n1],h1invsq,h2invsq);
	for j=2:n2-1
		colShift = n1*(j-1);
		i = 1 + colShift;
		@inbounds y[i] = op(x[i],0.0,x[i+1],x[i-n1],x[i+n1],h1invsq,h2invsq);
		@simd for i = (2 + colShift):(n1-1 + colShift)
			@inbounds y[i] = op(x[i],x[i-1],x[i+1],x[i-n1],x[i+n1],h1invsq,h2invsq);
		end
		i = n1 + colShift;
		@inbounds y[i] = op(x[i],x[i-1],0.0,x[i-n1],x[i+n1],h1invsq,h2invsq);
	end
	colShift = n1*(n2-1);
	i = 1 + colShift;
	@inbounds y[i] = op(x[i],0.0,x[i+1],x[i-n1],0.0,h1invsq,h2invsq);
	@simd for i = (2 + colShift):(n1-1 + colShift)
		@inbounds y[i] = op(x[i],x[i-1],x[i+1],x[i-n1],0.0,h1invsq,h2invsq);
	end
	i = n1 + colShift;
	@inbounds y[i] = op(x[i],x[i-1],0.0,x[i-n1],0.0,h1invsq,h2invsq);
else
end
end


text:
@testset "Attrition" begin

@testset "Constructor" begin
    attrition = Attrition()
    @test all( [attrition.name == "default", attrition.rates == [0.0],
        attrition.curvePoints == [0.0], attrition.period == 1.0] )
    attrition = Attrition( "Attr1" )
    @test attrition.name == "Attr1"
end  # @testset "Constructor"

attrition = Attrition( "Attrition" )

@testset "function setAttritionRate!" begin
    setAttritionRate!( attrition, 0.04 )
    @test attrition.rates == [0.04]
    setAttritionRate!( attrition, 0.025 )
    @test all( [attrition.rates == [0.025],
        attrition.lambdas == [- log( 0.975 )],
        attrition.gammas == [1.0]] )
    @test !setAttritionRate!( attrition, -0.5 )
    @test !setAttritionRate!( attrition, 2.5 )
end  # @testset "function setAttritionRate!"

@testset "function setAttritionPeriod!" begin
    setAttritionPeriod!( attrition, 12 )
    @test ( attrition.period == 12 ) &&
        ( attrition.lambdas == [- log( 0.975 ) / 12] )
    setAttritionPeriod!( attrition, 6 )
    @test attrition.period == 6 &&
        ( attrition.lambdas == [- log( 0.975 ) / 6] )
    @test !setAttritionPeriod!( attrition, -6 )
end  # @testset "function setAttritionRate!"

@testset "function setAttritionCurve!" begin
    setAttritionCurve!( attrition, Dict( 0.0 => 0.05, 24.0 => 0.03 ) )
    @test ( attrition.curvePoints == [0.0, 24.0] ) &&
        ( attrition.rates == [0.05, 0.03] )
    setAttritionCurve!( attrition, [0.0 0.04; 30.0 0.025] )
    @test ( attrition.curvePoints == [0.0, 30.0] ) &&
        ( attrition.rates == [0.04, 0.025] )
    @test setAttritionCurve!( attrition, [6.0 0.01; 30.0 0.02] )
    @test ( attrition.curvePoints == [0.0, 6.0, 30.0] ) &&
        ( attrition.rates == [0.0, 0.01, 0.02] )
    @test setAttritionCurve!( attrition,
        Dict( -6.0 => 0.02, 12.0 => 0.03, 30.0 => 0.01 ) )
    @test all( [attrition.curvePoints == [0.0, 12.0, 30.0],
        attrition.rates == [0.02, 0.03, 0.01],
        attrition.lambdas == - log.( 1 .- [0.02, 0.03, 0.01] ) / 6,
        attrition.gammas == [1.0, 0.98^2, 0.98^2 * 0.97^3]] )
    @test setAttritionCurve!( attrition,
        Dict( -12.0 => 0.02, -6.0 => 0.03, 24.0 => 0.01 ) )
    @test ( attrition.curvePoints == [0.0, 24.0] ) &&
        ( attrition.rates == [0.03, 0.01] )
    @test !setAttritionCurve!( attrition, Dict( 0.0 => -0.02, 30.0 => 1.5 ) )
    @test !setAttritionCurve!( attrition, [0.0 0.05 0.03; 30.0 0.02 0.05] )
    @test !setAttritionCurve!( attrition, [0.0 0.05; 18.0 0.03; 18.0 0.02] )
end  # @testset "function setAttritionCurve!"

end  # @testset "Attrition"


text:
using Fingerprints
const FP = Fingerprints
using Test

abstract type AAbstractType end
struct AStruct end
mutable struct AMutableStruct end
primitive type APrimitiveType 64 end

@testset "namestub" begin
    @test FP.namestub(AAbstractType)   === :AAbstractType
    @test FP.namestub(AStruct)         === :AStruct
    @test FP.namestub(AMutableStruct)  === :AMutableStruct
    @test FP.namestub(APrimitiveType)  === :APrimitiveType
    @test FP.namestub(Vector{Float64}) === :Array
    @test FP.namestub(Vector) === :Array
    @test FP.namestub(Array) === :Array
end

@testset "Kind" begin
    @test FP.Kind(AMutableStruct)      === FP.KindMutable()
    @test FP.Kind(AStruct)             === FP.KindStruct()
    @test FP.Kind(AAbstractType)       === FP.KindAbstract()
    @test FP.Kind(APrimitiveType)      === FP.KindPrimitive()
    @test FP.Kind(Vector)              === FP.KindUnionAll()
    @test FP.Kind(Union{Vector, Int})  === FP.KindUnion()
    @test FP.Kind(AbstractArray)       === FP.KindUnionAll()
    @test FP.Kind(AbstractVector{Int}) === FP.KindAbstract()
    @test FP.Kind(Vector{Int})         === FP.KindBuiltin()
    @test FP.Kind(Matrix{Int})         === FP.KindBuiltin()
    @test FP.Kind(Symbol)              === FP.KindBuiltin()
    @test FP.Kind(String)              === FP.KindBuiltin()
end

module A
    struct S
        x
        y
    end
end

module B
    struct S
        x
        y
    end
    struct T
        x
        y
    end
end

module C
    struct S
        x
        z
    end
end

@testset "Module name invariance" begin
    @test fingerprint(A.S(1,2)) === fingerprint(B.S(1,2))
    @test fingerprint(A.S(1,2)) !== fingerprint(B.T(1,2))
    @test fingerprint(A.S(1,2)) !== fingerprint(C.S(1,2))
end

struct E1 end
struct E2 end

@testset "empty structs" begin
    @test fingerprint(E1())  !== fingerprint(E2())
    @test fingerprint(E1())  !== fingerprint(E1)
    @test fingerprint(E1())  === fingerprint(E1())
    @test fingerprint("hi")  !== fingerprint(:hi)
end

@testset "mutable structs" begin
    mutable struct M1
        x
    end
    @test fingerprint(M1(1)) === fingerprint(M1(1))
    @test fingerprint(M1(1)) !== fingerprint(M1(2))
end

struct Container{T}
    inner::T
end
@testset "persistence" begin
    @test fingerprint(Container(Container("1"))) === 0x182a028725dce688
end

@test fingerprint([1,2]) != fingerprint((1,2))
@test fingerprint((1,2)) != fingerprint((a=1,b=2))
@test fingerprint((a=2,b=1)) != fingerprint((a=1,b=2))


text:
macro strided(ex)
    ex = macroexpand(__module__, ex)
    esc(_strided(ex))
end

function _strided(ex::Expr)
    if ex.head == :. && ex.args[2] isa QuoteNode # field access: wrap whole expression in maybestrided
        return Expr(:call, :(Strided.maybestrided), ex)
    elseif ex.head == :call && ex.args[1] isa Symbol
        if ex.args[1] == :reshape
            return Expr(:call, :(Strided.sreshape), map(_strided, ex.args[2:end])...)
        elseif ex.args[1] == :view
            return Expr(:call, :(Strided.sview), map(_strided, ex.args[2:end])...)
        else
            return Expr(:call, ex.args[1], map(_strided, ex.args[2:end])...)
        end
    elseif (ex.head == :(=) || ex.head == :(kw)) && ex.args[1] isa Symbol
        return Expr(ex.head, ex.args[1],
                    Expr(:call, :(Strided.maybeunstrided), _strided(ex.args[2])))
    elseif (ex.head == :(->))
        return Expr(ex.head, ex.args[1],
                    Expr(:call, :(Strided.maybeunstrided), _strided(ex.args[2])))
    else
        return Expr(ex.head, map(_strided, ex.args)...)
    end
end
const exclusionlist = Symbol[:(:)]
_strided(ex::Symbol) =  ex in exclusionlist ? ex : Expr(:call, :(Strided.maybestrided), ex)
_strided(ex) = ex

maybestrided(A::AbstractStridedView) = A
maybestrided(A::AbstractArray) = StridedView(A)
maybestrided(A) = A
maybeunstrided(A::StridedView) = copy(A).parent
maybeunstrided(A) = A

macro unsafe_strided(args...)
    syms = args[1:end-1]
    ex = macroexpand(__module__, args[end]) #_strided(args[end])
    all(isa(s, Symbol) for s in syms) ||
        error("The first arguments to `@unsafe_strided` must be variable names")
    ex = Expr(:let, Expr(:block, [:($s = Strided.UnsafeStridedView($s)) for s in syms]...), ex)
    return esc(:(GC.@preserve $(syms...) $ex))
end

# macro sfor(args...)
#     syms = args[1:end-1]
#     all(isa(s, Symbol) for s in syms) || error("The first arguments to `@sfor` must be variable names that will be usek")
#     ex = macroexpand(__module__, args[end])
#     ex = _sfor(syms, ex)
# end


text:
"""
    AtomPosition{D} <: AbstractRealSpaceData{D}

Contains information about an atom in a crystal, including atomic number, the name of an atom,
and its position. Atomic names may be set arbitrarily, but will default to the atomic symbol if not
provided explicitly.

An atomic number of 0 indicates a dummy atom - this does not represent a real atom, but a position
of note within the structure. If no name is given for a dummy atom, it will be given the empty
string.
"""
struct AtomPosition{D} <: AbstractRealSpaceData{D}
    name::String
    num::Int
    pos::SVector{D,Float64}
    function AtomPosition{D}(
        name::AbstractString,
        num::Integer,
        pos::AbstractVector{<:Real}
    ) where D
        # Don't allow empty names by default
        name = isempty(name) && num > 0 ? ELEMENTS[num] : name
        return new{D}(name, num, pos)
    end
end

function AtomPosition{D}(num::Integer, pos::AbstractVector{<:Real}) where D
    return AtomPosition{D}(ELEMENTS[num], num, pos)
end

function AtomPosition{D}(name::AbstractString, pos::AbstractVector{<:Real}) where D
    return AtomPosition{D}(name, get(ELEMENT_LOOKUP, name, 0), pos)
end

"""
    atomname(a::AtomPosition) -> String

Returns the name of the atom. Note that this returns the name that was provided during
construction, which is not guaranteed to be the same as the name of the element that has the given
atomic number.
"""
atomname(a::AtomPosition) = a.name

"""
    atomicno(a::AtomPosition) -> Int

Gets the atomic number of an atom in an atomic position.
"""
atomicno(a::AtomPosition) = a.num

"""
    coord(a::AtomPosition{D}) -> SVector{D,Float64}

Returns the coordinate associated with an atomic position.
"""
coord(a::AtomPosition) = a.pos

Base.getindex(a::AtomPosition, ind) = a.pos[ind]

"""
    AtomList{D} <: AbstractRealSpaceData{D}

A list of atomic positions with an associated basis. Atomic coordinates are given in terms of the
specified basis. If the basis is a zero matrix, the coordinates are assumed to be given in
angstroms.
"""
struct AtomList{D} <: AbstractRealSpaceData{D}
    basis::BasisVectors{D}
    coord::Vector{AtomPosition{D}}
    function AtomList{D}(
        basis::BasisVectors{D},
        coord::AbstractVector{AtomPosition{D}}
    ) where D
        # Remove any duplicate atoms if they come up
        return new{D}(basis, unique(coord))
    end
end

function AtomList{D}(
    basis::AbstractLattice{D},
    coord::AbstractVector{AtomPosition{D}};
    prim=false
) where D
    # Use either the primitive or conventional vectors depending on user input
    if prim
        return AtomList{D}(basis.prim, unique(coord))
    else
        return AtomList{D}(basis.conv, unique(coord))
    end
end

# Do it without adding a lattice
function AtomList{D}(coord::AbstractVector{AtomPosition{D}}) where D
    return AtomList{D}(zeros(BasisVectors{D}), coord)
end

# Check if an AtomList{D} is empty
Base.isempty(l::AtomList) = isempty(l.coord)
Base.getindex(l::AtomList, ind) = l.coord[ind]

# Iterate through an AtomList
Base.iterate(l::AtomList) = iterate(l.coord)
Base.iterate(l::AtomList, n) = iterate(l.coord, n)

natom(l::AtomList) = length(l.coord)
basis(l::AtomList) = l.basis

"""
    cartesian(l::AtomList{D}) -> AtomList{D}

Converts an `AtomList` from reduced coordinates (relative to some crystal basis) to Cartesian
coordinates in space.
"""
function cartesian(l::AtomList{D}) where D
    # If there are no basis vectors specified, just return the original list
    # If the basis vectors are zero we're assuming Cartesian coordinates
    basis(l) == zeros(BasisVectors{D}) && return l
    newlist = map(a -> cartesian(l.basis, a), l.coord)
    return AtomList{D}(zeros(BasisVectors{D}), newlist)
end

# TODO: clean up these method definitions into something that makes more sense
"""
    cartesian(b::AbstractMatrix{<:Real}, a::AtomPosition{D}) -> AtomPosition{D}

Converts an `AtomPosition` defined in terms of basis `b` to a new `AtomPosition` defined in
Cartesian coordinates.
"""
function cartesian(b::AbstractMatrix{<:Real}, a::AtomPosition{D}) where D
    newpos = b * a.pos
    return AtomPosition{D}(a.name, a.num, newpos)
end

cartesian(b::BasisVectors{D}, a::AtomPosition{D}) where D = cartesian(matrix(b), a)

"""
    reduce_coords(basis::AbstractMatrix{<:Real}, a::AtomPosition; incell=false)

Convert a coordinate from a Cartesian basis to the crystal basis. If `incell` is true, the position
vector components will be truncated so they lie within the cell bounds (between 0 and 1).
"""
function reduce_coords(
    basis::AbstractMatrix{<:Real},
    a::AtomPosition{D};
    incell::Bool=false
) where D
    v = basis\coord(a)
    # If it needs to be in a cell, truncate the integer portion of the elements
    v = v - incell * floor.(v)
    return AtomPosition{D}(name(a), atomicno(a), v)
end

"""
    reduce_coords(basis::AbstractMatrix{<:Real}, a::AbstractVector{AtomPosition}; incell=false)

Convert a vector of atomic coordinates from a Cartesian basis to the crystal basis. If `incell` is
true, the position vector components will be truncated so they lie within the cell bounds (between
0 and 1).
"""
function reduce_coords(
    basis::AbstractMatrix{<:Real},
    va::AbstractVector{AtomPosition{D}};
    incell::Bool = false
) where D
    va_new = map(va) do a
        reduce_coords(basis, a, incell = incell)
    end
    return AtomPosition{D}(basis, va_new)
end


text:
SRBench["iterating"] = BenchmarkGroup();

itr_OneTo() = for i in Base.OneTo(1000); i; end
itr_OneToSRange() = for i in OneToSRange(1000); i; end
itr_OneToSRangeType() = for i in OneToSRange{1000}; i; end

SRBench["iterating"]["OneTo"] = BenchmarkGroup();
SRBench["iterating"]["OneTo"]["<:AbstractRange"] = @benchmarkable itr_OneTo()
SRBench["iterating"]["OneTo"]["OneToSRange"] = @benchmarkable itr_OneToSRange()
SRBench["iterating"]["OneTo"]["Type{OneToSRange}"] = @benchmarkable itr_OneToSRange()

# Indexing into OneTo iterators
SRBench["iterating"]["OneTo[i]"] = BenchmarkGroup();

itr_idx_OneTo() = for i in 1:100; Base.OneTo(1000)[i]; end
itr_idx_OneToSRange() = for i in 1:100; OneToSRange(1000)[i]; end
itr_idx_OneToSRangeType() = for i in 1:100; OneToSRange{1000}[i]; end

SRBench["iterating"]["OneTo[i]"]["<:AbstractRange"] = @benchmarkable itr_idx_OneTo()
SRBench["iterating"]["OneTo[i]"]["OneToSRange"] = @benchmarkable itr_OneToSRange()
SRBench["iterating"]["OneTo[i]"]["Type{OneToSRange}"] = @benchmarkable itr_OneToSRangeType()

# Iterating by steps
SRBench["iterating"]["1:2:1000"] = BenchmarkGroup();

itr_step_r() = for i in 1:2:1000; i; end
itr_step_sr() = for i in srange(1,1000,step=2); i; end
itr_step_srtype() = for i in StaticRange{Int,1,1000,2,500}; i; end

SRBench["iterating"]["1:2:1000"]["<:AbstractRange"] = @benchmarkable itr_step_r()
SRBench["iterating"]["1:2:1000"]["StaticRange"] = @benchmarkable itr_step_sr()
SRBench["iterating"]["1:2:1000"]["Type{StaticRange}"] = @benchmarkable itr_step_sr()

# Iterating into multidimensional linear indices
SRBench["iterating"]["LinearIndices((5,6))"] = BenchmarkGroup();

itr_li() = for i in LinearIndices((5,6)); i; end
itr_sli() = for i in LinearSIndices((5,6)); i; end

SRBench["iterating"]["LinearIndices((5,6))"]["LinearIndices"] = @benchmarkable itr_li()
SRBench["iterating"]["LinearIndices((5,6))"]["SIndices"] = @benchmarkable itr_sli()

### Sub Indexing
# TODO Figure out how to best test SubIndices
SRBench["iterating"]["LinearIndices((5,6))[2:3,1:3]"] = BenchmarkGroup();

# Indexing into
function idx_li()
    a = LinearIndices((5,6))[2:3,1:3]
    a[2,1]
end

function idx_li()
    a = LinearSIndices((5,6))[2:3,1:3]
    a[2,1]
end

function itr_sub_li()
    for i in
        i
    end
end

function itr_sub_sli()
    for i in LinearSIndices((5,6))[2:3,1:3]
        i
    end
end

function itr_sub_subli()
    a = SubLinearIndices(LinearSIndices((5,6)), (srange(2:3), srange(1:3)))
    for i in
        i
    end
end

function while_r()
    r = 1:1000
    i, state = iterate(r)
    while i < 1000
        i, state = iterate(r, state)
    end
end
SRBench["iterating"]["while_r"] = @benchmarkable while_r()

function while_sr()
    r = srange(1,1000,step=1)
    i, state = iterate(r)
    while i < 1000
        i, state = iterate(r, state)
    end
end
SRBench["iterating"]["while_sr"] = @benchmarkable while_sr()

function while_step_r()
    r = 1:2:1000
    i, state = iterate(r)
    while i < lastindex(r)
        i, state = iterate(r, state)
    end
end
SRBench["iterating"]["while_step_r"] = @benchmarkable while_step_r()

function while_step_sr()
    r = srange(1,1000,step=2)
    i, state = iterate(r)
    while i < lastindex(r)
        i, state = iterate(r, state)
    end
end

function itr_subindices()
    r = SubIndices((10,10), srange(1,10,step=1), srange(1,10,step=1))
    i, state = iterate(r)
    for i in 1:100
        r[i]
    end
end

function itr_linearindices()
    r = LinearIndices((10,10))
    i, state = iterate(r)
    while i < lastindex(r)
        i, state = iterate(r, state)
    end
end

function itr_staticindices()
    r = StaticIndices((10,10))
    i, state = iterate(r)
    while i < lastindex(r)
        i, state = iterate(r, state)
    end
end

SRBench["iterating"]["while_step_sr"] = @benchmarkable while_step_sr()


text:
using Test

function runtests_mpirun()
    nprocs = min(4, Sys.CPU_THREADS)
    testdir = dirname(@__FILE__)
    testfiles = ["lav.jl", "lavdense.jl", "matrix.jl", "distmatrix.jl", "props.jl", "generic.jl", "spectral.jl"]
    nfail = 0
    @info "Running Elemental.jl tests"
    for f in testfiles
        try
            run(`mpirun -np $nprocs $(Base.julia_cmd()) $(joinpath(testdir, f))`)
            Base.with_output_color(:green,stdout) do io
                println(io,"\tSUCCESS: $f")
            end
        catch ex
            Base.with_output_color(:red,stderr) do io
                println(io,"\tError: $f")
                showerror(io,ex,backtrace())
            end
            nfail += 1
        end
    end
    return nfail
end

function runtests_repl()
    nprocs = min(4, Sys.CPU_THREADS)
    exename = joinpath(Sys.BINDIR, Base.julia_exename())
    testdir = dirname(@__FILE__)
    testfiles = ["darray.jl"]
    nfail = 0
    @info "Running Elemental.jl tests"
    for f in testfiles
        try
            # FixMe! We temporarily run Finalize() explictly on the workers because the atexit hook
            # doesn't seem to be correctly triggered on workers as of 31 October 2018.
            cmdstr = "using Distributed, MPI; man = MPIManager(np = $nprocs); addprocs(man); include(\"$(joinpath(testdir, f))\"); asyncmap(p -> remotecall_fetch(() -> Elemental.Finalize(), p), workers())"
            run(`$exename -e $cmdstr`)
            Base.with_output_color(:green,stdout) do io
                println(io,"\tSUCCESS: $f")
            end
        catch ex
            Base.with_output_color(:red,stderr) do io
                println(io,"\tError: $f")
                showerror(io,ex,backtrace())
            end
            nfail += 1
        end
    end
    return nfail
end

function runtests()
    return runtests_mpirun() & runtests_repl()
end

exit(runtests())


text:
using MuZero
using Grids

# cartesian_coords_to_action(x,y) = (y-1)*3 + x

# function legal_actions(state::TicTacToeState)
#     game_is_won(state) ? [] : findall(iszero, calc_cell_states(state))
# end

# function dynamics(state::TicTacToeState, action::Integer)
#     player = to_play(state)
#     next_state = TicTacToeState(bits(state) | action_mask(player, action))
#     reward = last_reward(next_state)
#     reward, next_state
# end

# function to_image(state::TicTacToeState)
#     cell_states = calc_cell_states(state)
#     board = reverse(reshape(cell_states, 3, 3)', dims=1)
#     player_multiplier = to_play(state) == PLAYER_1 ? Float32(1) : Float32(-1)
#     player1_channel = Float32(1)*(board .== PLAYER_1)
#     player2_channel = Float32(-1)*(board .== PLAYER_2)
#     to_play_channel = player_multiplier*ones(Float32, 3, 3)
#     channels = cat(player1_channel, player2_channel, to_play_channel; dims=3)
# end

# const PLAYER_1 = 1
# const PLAYER_2 = 2

# game_is_won(state::TicTacToeState) = get_winner(state) > 0
# action_mask(player, action) = UInt32(player) << ((action-1)*2)
# just_played(arg) = to_play(arg) == PLAYER_1 ? PLAYER_2 : PLAYER_1
# bits(x::TicTacToeState) = x.bits

# function calc_cell_states(state::TicTacToeState)
#     cell_mask = UInt32(0b11)
#     function cell_state(cell_idx)
#         res = (bits(state) & (cell_mask << UInt32(2*(cell_idx-1)))) >>> UInt32(2*(cell_idx-1))
#         Int8(res)
#     end
#     map(cell_state, collect(1:9))
# end

# function calc_cell_state_counts(state::TicTacToeState)
#     cell_states = calc_cell_states(state)
#     not_played,p1,p2 = 0,0,0
#     for cell_state in cell_states
#         if cell_state == 0
#             not_played += 1
#         elseif cell_state == PLAYER_1
#             p1 += 1
#         elseif cell_state == PLAYER_2
#             p2 += 1
#         end
#     end
#     not_played, p1, p2
# end

# function to_play(state::TicTacToeState)
#     _, p1_actions, p2_actions = calc_cell_state_counts(state)
#     p1_actions == p2_actions ? PLAYER_1 : PLAYER_2
# end

# function last_reward(state::TicTacToeState)
#     winner = get_winner(state)
#     if winner == 0
#         # Either the game isn't over, or the game is over but stalemated.
#         return Float32(0)
#     end
#     winner == just_played(state) ? Float32(1) : Float32(-1)
# end

const PLAYER_1 = 1
const PLAYER_2 = 2

struct TicTacToeState <: AbstractState
    bits::UInt32
end

MuZero.key(x::TicTacToeState) = x.bits

function TicTacToeState(state_idx=0)
    TicTacToeState(UInt32(state_idx))
end

action_mask(player, action) = UInt32(player) << ((action-1)*2)

function MuZero.legal_actions(state::TicTacToeState)
    if game_is_won(state)
        return []
    end
    findall(iszero, calc_cell_states(state))
end

MuZero.hex_str(state::TicTacToeState) = "0x$(string(key(state), base=16, pad=8))"
Base.copy(state::TicTacToeState) = TicTacToeState(key(state))
is_terminal(state::TicTacToeState) = length(legal_actions(state)) == 0
game_is_won(state::TicTacToeState) = get_winner(state) > 0
islegal(state::TicTacToeState, action::Integer) = in(action, legal_actions(state))

MuZero.next_state(state::TicTacToeState, action::Integer) = next_state_bits(state, action) |> TicTacToeState
MuZero.next_state(state::TicTacToeState, actions::Vector{<:Integer}) = next_state_bits(state, actions) .|> TicTacToeState

next_state_bits(state::TicTacToeState, actions::Vector{<:Integer}) = [next_state_bits(state, a) for a in actions]

function MuZero.get_reward(state::TicTacToeState, action::Integer)
    next_state(state, action) |> last_reward
end

function MuZero.last_reward(state::TicTacToeState)
    winner = get_winner(state)
    if winner == 0
        # Either the game isn't over, or the game is over but stalemated.
        return Float32(0)
    end
    winner == just_played(state) ? Float32(1) : Float32(-1)
end

function next_state_bits(state::TicTacToeState, action::Integer)
    player = to_play(state)
    state.bits | action_mask(player, action)
end

function calc_cell_states(state::TicTacToeState)
    cell_mask = UInt32(0b11)
    function cell_state(cell_idx)
        res = (state.bits & (cell_mask << UInt32(2*(cell_idx-1)))) >>> UInt32(2*(cell_idx-1))
        Int8(res)
    end
    map(cell_state, collect(1:9))
end

function calc_cell_state_counts(state::TicTacToeState)
    cell_states = calc_cell_states(state)
    not_played,p1,p2 = 0,0,0
    for cell_state in cell_states
        if cell_state == 0
            not_played += 1
        elseif cell_state == PLAYER_1
            p1 += 1
        elseif cell_state == PLAYER_2
            p2 += 1
        end
    end
    not_played, p1, p2
end

function to_play(state::TicTacToeState)
    _, p1_actions, p2_actions = calc_cell_state_counts(state)
    p1_actions == p2_actions ? PLAYER_1 : PLAYER_2
end

function get_winner(state::TicTacToeState)
    winner_masks = [
        (0b111111000000000000, 0b010101000000000000, 0b101010000000000000), # Horizonal top
        (0b000000111111000000, 0b000000010101000000, 0b000000101010000000), # Horizonal middle
        (0b000000000000111111, 0b000000000000010101, 0b000000000000101010), # Horizonal bottom
        (0b110000110000110000, 0b010000010000010000, 0b100000100000100000), # Vertical left
        (0b001100001100001100, 0b000100000100000100, 0b001000001000001000), # Vertical center
        (0b000011000011000011, 0b000001000001000001, 0b000010000010000010), # Vertical right
        (0b110000001100000011, 0b010000000100000001, 0b100000001000000010), # Diagonal up
        (0b000011001100110000, 0b000001000100010000, 0b000010001000100000), # Diagonal down
    ]
    for (mask, player1, player2) in winner_masks
        res = state.bits & mask
        if res == player1
            return PLAYER_1
        elseif res == player2
            return PLAYER_2
        end
    end
    0
end

function MuZero.dynamics(state::TicTacToeState, action::Integer)
    player = to_play(state)
    next_state = TicTacToeState(state.bits | action_mask(player, action))
    reward = last_reward(next_state)
    reward, next_state
end

function valid_idx(state::TicTacToeState)
    cell_states = calc_cell_states(state)
    invalid_cell_mask = Int8(0b11)
    if in(invalid_cell_mask, cell_states)
        return false
    end
    _, p1_actions, p2_actions = calc_cell_state_counts(state)
    (p1_actions == p2_actions) || (p1_actions == (p2_actions+1))
end

function make_image(state::TicTacToeState)
    cell_states = calc_cell_states(state)
    board = reverse(reshape(cell_states, 3, 3)', dims=1)
    to_play_multiplier = to_play(state) == PLAYER_1 ? Float32(1) : Float32(-1)
    player1_channel = Float32(1)*(board .== PLAYER_1)
    player2_channel = Float32(-1)*(board .== PLAYER_2)
    to_play_channel = to_play_multiplier*ones(Float32, 3, 3)
    channels = cat(player1_channel, player2_channel, to_play_channel; dims=3)
    Image(channels)
end

function perfect_policy(curr_state::TicTacToeState; show_v_next=false)
    v_next = map(1:9) do action
        if !islegal(curr_state, action)
            return -Inf
        end
        state_next = next_state(curr_state, action)
        return last_reward(state_next) + -1*value_fn[key(state_next)]
    end
    if show_v_next
        @show v_next
    end
    max_v_next = maximum(v_next)
    num_matching = v_next[v_next .== max_v_next] |> length
    max_pr = 1.0f0 / num_matching
    map(x->x == max_v_next ? max_pr : 0.0f0, v_next)
end

function print_board(state::TicTacToeState)
    function cell_fn(x,y)
        action = cartesian_coords_to_action(x,y)
        cell_player = calc_cell_states(state)[action]
        if cell_player == PLAYER_1
            return "X"
        elseif cell_player == PLAYER_2
            return "O"
        end
        nothing
    end
    show(stdout, "text/plain", Grid(3,3)(cell_fn))
end

function init_grid_cell_fn(s::TicTacToeState)
    return (x,y) -> begin
        action = cartesian_coords_to_action(x,y)
        cell_player = calc_cell_states(s)[action]
        if cell_player == PLAYER_1
            return "X"
        elseif cell_player == PLAYER_2
            return "O"
        end
        nothing
    end
end

function Base.show(io::IO, m::MIME"text/plain", s::TicTacToeState)
    cell_fn = init_grid_cell_fn(s)
    g = Grid(3,3)
    show(io, m, g(cell_fn))
end

struct TicTacToeStateSVG
    s::TicTacToeState
end

svg(s::TicTacToeState) = TicTacToeStateSVG(s)

function Base.show(io::IO, m::MIME"image/svg+xml", state_svg::TicTacToeStateSVG)
    cell_fn = init_grid_cell_fn(state_svg.s)
    g = Grid(3,3) |> Grids.svg
    show(io, m, g(cell_fn))
end


text:
"Null configuration"
type NullConfiguration <: Configuration
    name::AbstractString
    source::AbstractString
    state::LifeCycle.State
    root::LoggerConfig

    NullConfiguration() = new("Null", "", LifeCycle.INITIALIZED, LoggerConfig(Level.OFF))
end
appender(cfg::NullConfiguration, name::AbstractString) = nothing
appenders(cfg::NullConfiguration) = APPENDERS()
appender!(cfg::NullConfiguration, nm::AbstractString, apnd::Appender) = nothing
logger(cfg::NullConfiguration, name::AbstractString) = cfg.root
loggers(cfg::NullConfiguration) = LOGCONFIGS()
logger!(cfg::NullConfiguration, nm::AbstractString, lc::LoggerConfig) = nothing
state(cfg::NullConfiguration) = LifeCycle.STOPPED
configure(cfg::NullConfiguration) = nothing
filter(cfg::NullConfiguration) = FILTER()

# Register configuration type
LOG4JL_CONFIG_TYPES[:NULL] = NullConfiguration

"Default Log4jl configuration"
type DefaultConfiguration <: Configuration
    name::AbstractString
    source::AbstractString
    state::LifeCycle.State
    root::LoggerConfig

    properties::PROPERTIES
    appenders::APPENDERS
    loggers::LOGCONFIGS
    filter::FILTER

    DefaultConfiguration() =
        new("Default", "", LifeCycle.INITIALIZED, LoggerConfig(), PROPERTIES(), APPENDERS(), LOGCONFIGS(), FILTER())
end
appender(cfg::DefaultConfiguration, name::AbstractString) = get(cfg.appenders, name, nothing)
appenders(cfg::DefaultConfiguration) = cfg.appenders
loggers(cfg::DefaultConfiguration) = cfg.loggers

function configure(cfg::DefaultConfiguration)
    # Add basic console appender
    appender!(cfg, Default = Appenders.Console(layout = Layouts.BasicLayout()))

    # Reference appender to root configuration
    reference!(cfg, ROOT_LOGGER_NAME => "Default")
    return
end

# Register configuration type
LOG4JL_CONFIG_TYPES[:DEFAULT] = DefaultConfiguration


text:
abstract type ParallelType end
struct Parallel{T1<:AbstractVector, T2<:AbstractVector, T3 <: Symbol, T4<: Real} <: ParallelType
    real::T1
    simulated::T1
    simulated_bounds::T2
    resampled::T1
    resampled_bounds::T2
    iter::T4
    reduction_method::T3
    correlation_method::T3
end

struct ParallelSets{T <: ParallelType}
    FA::T
    PCA::T
end

"""
    parallel(data, niter, f = fa)
Parallel Analysis.

# Arguments
- `data`
- `niter` is a size of simulation that generate (normal) random data matrix, calculate R, correlation matrix and its eigen values.
- `f` Function for dimension reduction. Default is `fa`.

# Values
- `real` Eigen values from read data.
- `simulated`, `resampled` Mean of eigen values from simulated data.
- `simulated_bound`, `resampled_bounds` 2.5th and 97.5th perceintile rank values of simulated eigen values.
- `iter` Number of iterations (of the simulation).
- `reduncion_method`
- `correlation_method`

# Example
```julia
julia> using ParallelAnalysis, Random, StatsPlots
julia> begin Random.seed!(1234)
           a = rand(Uniform(0.5, 2.0), 30)
           b = [sort(rand(Uniform(-3, 3), 4); rev = false) for i in 1:30]
           θ = rand(Normal(0, 1), 3000)
           resp = generate_response(θ, a, b)
       end;

julia> par_fit1 = parallel(resp, 10, x -> fa(x; cor_method = :Polychoric))
Progress: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:02
Parallel Analysis:
    Suggested that the number of factors is 1 (based on resampling)
    Suggested that the number of componets is 1 (based on resampling)

julia> plot(par_fit1) # visualize a result of parallel analysis

julia> plot(par_fit1, markershape = :none) # Don't show noisy makers.

```

"""
function parallel(data, niter, f = fa)
    n, j = size(data)
    fit = f(data)
    # println("Analyze read data.")
    if fit.cor == :Pearson
        V = cor(Matrix(data))
    elseif fit.cor == :Polychoric
        V = polycor(data)
    end
    eig_real_pca = sort!(eigvals(Symmetric(V)); rev = true)
    V[diagind(V)] = communalities(fit)
    eig_real_fa = sort!(eigvals(Symmetric(V)); rev = true)

    # println("Prepare matrix.")
    eig_sim_fa = Matrix{eltype(eig_real_fa)}(undef, niter, length(eig_real_fa))
    eig_sim_pca = Matrix{eltype(eig_real_pca)}(undef, niter, length(eig_real_pca))
    rsm_sim_fa = similar(eig_sim_fa)
    rsm_sim_pca = similar(eig_sim_pca)

    prg = Progress(niter)
    # println("Start simulation!")
    Threads.@threads for i in 1:niter

        W = random_sample(data)
        Z = randn(n, j)

        # resample data
        rsm_fit = f(W)
        M = rsm_fit.mat
        rsm_sim_pca[i, :] = sort!(eigvals(Symmetric(M)); rev = true)
        M[diagind(M)] = communalities(rsm_fit)
        rsm_sim_fa[i, :] = sort!(eigvals(Symmetric(M)); rev = true)

        # simulated data
        sim_fit = fa_pearson(Z)
        R = sim_fit.mat
        eig_sim_pca[i, :] = sort!(eigvals(Symmetric(R)); rev = true)
        R[diagind(R)] = communalities(sim_fit)
        eig_sim_fa[i, :] = sort!(eigvals(Symmetric(R)); rev = true)

        next!(prg)
    end

    # summarize
    eig_sim_fa_mean = mean(eig_sim_fa, dims = 1)[:]
    eig_sim_fa_bounds = map(x -> quantile(x, [0.025, 0.975]), eachcol(eig_sim_fa))
    rsm_sim_fa_mean = mean(rsm_sim_fa, dims = 1)[:]
    rsm_sim_fa_bounds = map(x -> quantile(x, [0.025, 0.975]), eachcol(rsm_sim_fa))

    eig_sim_pca_mean = mean(eig_sim_pca, dims = 1)[:]
    eig_sim_pca_bounds = map(x -> quantile(x, [0.025, 0.975]), eachcol(eig_sim_pca))
    rsm_sim_pca_mean = mean(rsm_sim_pca, dims = 1)[:]
    rsm_sim_pca_bounds = map(x -> quantile(x, [0.025, 0.975]), eachcol(rsm_sim_pca))

    ps = ParallelSets(
        Parallel(eig_real_fa, eig_sim_fa_mean, eig_sim_fa_bounds, rsm_sim_fa_mean, rsm_sim_fa_bounds, niter, :FA, Symbol(fit.cor)),
        Parallel(eig_real_pca, eig_sim_pca_mean, eig_sim_pca_bounds, rsm_sim_pca_mean, rsm_sim_pca_bounds, niter, :PCA, Symbol(fit.cor))
    )

    return ps
end

function findnfactors(x, y)
    for i in axes(x, 1)
        comp = x[i:end] .> y[i:end]
        findfirst(comp) != 1 && return i - 1
    end
    return 0
end

function Base.show(io::IO, x::Parallel)
    println(io, "Parallel Analysis: \nDimension reduction: $(x.reduction_method)")
    println(io, "Correlation method: $(x.correlation_method)")
    println(io, "Simulation size: $(x.iter)")
    println(io, "Suggested the number of factors is $(findnfactors(x.real, x.resampled)) (based on resampling)")
end

function Base.show(io::IO, x::ParallelSets)
    println(io, "Parallel Analysis:")
    println(io, "   Suggested that the number of factors is $(findnfactors(x.FA.real, x.FA.resampled)) (based on resampling)")
    println(io, "   Suggested that the number of componets is $(findnfactors(x.PCA.real, x.PCA.resampled)) (based on resampling)")
end

@recipe function f(pa::Parallel)

    @series begin
        # simulated mean
        label --> "$(pa.reduction_method) - Simulated mean"
        linecolor --> :red
        fillcolor --> :red
        linestyle --> :dash
        bds = hcat(pa.simulated_bounds...)
        l = bds[1, :] .- pa.simulated
        u = pa.simulated .- bds[2, :]
        ribbon := (l, u)
        pa.simulated
    end

    @series begin
        # simulated mean
        label --> "$(pa.reduction_method) - Resampled mean"
        linecolor --> :blue
        fillcolor --> :blue
        linestyle --> :dash
        bds = hcat(pa.resampled_bounds...)
        l = bds[1, :] .- pa.resampled
        u = pa.resampled .- bds[2, :]
        ribbon := (l, u)
        pa.resampled
    end

    @series begin
        label := ""
        linecolor --> :black
        linestyle --> :dot
        seriestype := :hline
        [1]
    end

    label --> "$(pa.reduction_method) - Real data"
    linecolor --> :black
    yguide --> "Eigen values"
    xguide --> "The number of components"
    title --> "Parallel Analysis over $(pa.iter) simulation."
    ylims --> (-Inf, Inf)
    pa.real
end

@recipe function f(x::ParallelSets)
    @series begin
        markershape --> :dtriangle
        markeralpha --> 0.5
        markercolor --> :white
        x.PCA
    end
    @series begin
        markershape --> :utriangle
        markeralpha --> 0.5
        markercolor --> :white
        x.FA
    end
    nothing
end


text:
_convert_dims(N::Integer, d::Integer) = N-d
_convert_dims(N::Integer, d::Dims)    = N.-d #map(i->_invert_dim(N, i), d)
_convert_dims(N::Integer, d::Integer...) = _convert_dims(N, Dims(d))
_convert_dims(N::Integer, d::Colon) = nothing

_convert_dims(::Dims{N}, d) where N = _convert_dims(N, d)
_convert_dims(::AbstractArray{T,N}, d) where {T,N} = _convert_dims(N, d)

_autokeep(dims::Colon) = false
_autokeep(dims) = true


text:
"""
Package: ForecastPlots

    candel(x, t, args...; colors = (:green, :red), kw...)

Plot a candlestick for stock prices

# Arguments
- `x`: Matrix with four columns for Open, High, Low and Close values.
- `t`: Vector with time units, it defaults to 1:size(x,1)
- `args...`: `plot` arguments
- `colors`: Tuple with the up and down color symbols, it defaults to (:green,:red)
- `kw...`: `plot` keyword arguments

# Returns
Candle plot

# Examples
```julia-repl
using MarketData

x = yahoo(:INTC)
last_month_intc = values(x)[end-30:end,1:4]
candle(last_month_intc)
candle(last_month_intc,timestamp(x)[end-30:end])
candle(last_month_intc,timestamp(x)[end-30:end]; colors=(:white,:black), title = "INTC")
```
"""
function candle(x::Matrix{<:Real},t=1:size(x,1); colors=(:green,:red), kw...)

    @assert length(colors) == 2 "Two colors required"
    @assert eltype(colors) == Symbol "'colors' must be symbols e.g. (:white,:black)"
    @assert colors[1] != colors[2] "'colors' must be different"

    plot(t,repeat([missing],size(x,1));labels=nothing,kw...)

    ttype = eltype(t) <: TimeType
    t = ttype ? Dates.value.(t) : t
    woc = ttype ? (t[2]-t[1])/2  : 1/2
    whl = woc/10

    rect(w, h, x, y) = Shape(x .+ [0, w, w, 0, 0], y .+ [0, 0, h, h, 0])

    for (i,ohlc) in zip(t,eachrow(x))

        # High / Low
        plot!(rect(whl,ohlc[2]-ohlc[3],i,ohlc[3]);
              seriescolor = :black,
              linewidth = 0,
              label = nothing)

        # Open / Close
        plot!(rect(woc,ohlc[4]-ohlc[1],i-woc/2,ohlc[1]);
              seriescolor = ohlc[1]<ohlc[4] ? colors[1] : colors[2],
              linewidth = 0,
              label = nothing)

    end
    current()
end


text:
import Singular: fetch, imap

function invariants(cforms::ClosedFormSystem)
    I, loopvars, expvars, lc = preprocess([cforms], singleloop = true)[1]
    elim = collect(Iterators.drop(Singular.gens(base_ring(I)), length(loopvars)*2))
    Singular.eliminate(I, elim...)
end

function invariants(loop::SingleLoop)
    invariants(closed_forms(loop))
end

function invariants(loop::MultiLoop)
    invariants(closed_forms(loop))
end

function invariants(loops::Array{ClosedFormSystem,1})
    if length(loops) == 0
        warn("No loops given!")
        return
    end
    loopvars = string.(loops[1].vars)
    inivars = ["$(v)_0" for v in loopvars]
    R, _ = PolynomialRing(QQ, [inivars; loopvars])

    preprocessed = preprocess(loops)

    index = 0
    I_new = initial_ideal(loopvars)
    I_o = 1
    I_n = nothing
    while I_n != I_o
        I_o = I_n
        for sys in preprocessed
            I_new = invariants(I_new, sys..., index)
            # println("New ideal: ", I_new)
            index += 1
        end

        idxvars = ["$(v)_$(index+1)" for v in loopvars]
        S, _ = PolynomialRing(QQ, [inivars; loopvars; idxvars])
        # println("Container: ", S)
        # rename variables to have v instead of v_index
        b, elim = renaming_polys(loopvars, index)
        # println("Basis (before map): ", b)
        b = [imap(g, S) for g in b]
        # println("Basis (after map): ", b)
        B = Singular.Ideal(S, b)

        elim = imap.(elim, Ref(S))

        I_n = Singular.eliminate(imap(I_new, S) + B, elim...)
        I_n = groebner(imap(I_n, R))
        # println("Final ideal: ", I_n)
        if index > (length(loopvars)+1) * length(loops)
            warn("something went wrong")
        end
    end
    I_n
end

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

function preprocess(loops::Array{ClosedFormSystem,1}; singleloop::Bool = false)
    # assume loops[i].vars == loops[j].vars for all i,j

    # cfslist = [rec_solve(sl.body) for sl in loops]
    # cfslist = loops.cforms
    # println("Closed forms:", cfslist)

    preprocessed = []
    for (i, loop) in enumerate(loops)
        cfs = loop.cforms
        exp = union(exponentials.(cfs)...)
        exp = filter(x -> x!=Sym(1), exp)

        # collect all variables
        lc = string(loops[i].lc)
        loopvars = string.(loops[i].vars)
        expvars = ["vvv_$(i)" for i in 1:length(exp)]

        vars, _ = var_order(loopvars, expvars, lc, 0)
        if singleloop
            vars, _ = var_order_single(loopvars, expvars, lc, 0)
        end

        R, svars = PolynomialRing(QQ, vars)
        varmap = Dict(zip(vars, svars))

        # replace exponentials in closed forms with variables
        if !isempty(expvars)
            expsym = Sym.(expvars)
            expmap = Dict(zip(exp, expsym))
            for cf in cfs
                expvars!(cf, expmap)
            end
        end

        # generate basis of closed forms
        function polyfn(cf)
            if singleloop
                expr = Sym("$(Sym(cf.f.x))") - replace_functions(polynomial(cf), 0)
            else
                expr = Sym("$(Sym(cf.f.x))_2") - replace_functions(polynomial(cf), 1)
            end
            # workaround as libsingular crashes when "real" rationals are contained
            return clear_denom(expr)
        end
        basis = polyfn.(cfs)

        sbasis = [sym2spoly(p, varmap) for p in basis]
        I = Singular.Ideal(R, sbasis)

        # dependencies
        if !isempty(exp)
            # println("Something with dependencies?????")
            A = AlgebraicDependencies.dependencies(sideal, exp, variables=expsym)
            if A != nothing
                # println("Something with imap?????")
                I += imap(A, R)
            else
                @debug "No algebraic dependencies among exponentials"
            end
        end

        res = (I, loopvars, expvars, lc)
        push!(preprocessed, res)
    end
    # println(preprocessed)
    preprocessed
end

function invariants(I_old::sideal, I_new::sideal, loopvars::Array{String,1}, expvars::Array{String,1}, lc::String, index::Int)

    vars, elimcnt = var_order(loopvars, expvars, lc, index)
    R, svars = PolynomialRing(QQ, vars)
    # map variables via order
    I_new = fetch(I_new, R)
    # map variables via name
    I_old = imap(I_old, R)

    elimvars = collect(Iterators.drop(svars, elimcnt))
    Singular.eliminate(I_old + I_new, elimvars...)
end

function var_order(loopvars::Array{String,1}, expvars::Array{String,1}, lc::String, index::Int)
    inivars = ["$(v)_0" for v in loopvars]
    finvars = ["$(v)_$(index+2)" for v in loopvars]
    midvars = ["$(v)_$(index+1)" for v in loopvars]
    auxvars = [expvars; lc]

    vars = [inivars; finvars; midvars; auxvars], length(inivars) + length(finvars) # also return number of variables not to eliminate
end

function var_order_single(loopvars::Array{String,1}, expvars::Array{String,1}, lc::String, index::Int)
    inivars = ["$(v)_0" for v in loopvars]
    finvars = ["$(v)" for v in loopvars]
    auxvars = [expvars; lc]

    vars = [inivars; finvars; auxvars], length(inivars) + length(finvars) # also return number of variables not to eliminate
end

function renaming_polys(loopvars::Array{String,1}, index::Int)
    midvars = ["$(v)_$(index+1)" for v in loopvars]
    # finvars = ["$(v)" for v in loopvars]

    _, svars = PolynomialRing(QQ, [midvars; loopvars])

    midcnt = length(midvars)
    midsvars = collect(Iterators.take(svars, midcnt))
    finsvars = collect(Iterators.drop(svars, midcnt))

    [ x - y for (x, y) in zip(midsvars, finsvars)], midsvars # also return variables to eliminate
end

function initial_ideal(loopvars::Array{String,1})
    inivars = ["$(v)_0" for v in loopvars]
    finvars = ["$(v)_1" for v in loopvars]

    R, svars = PolynomialRing(QQ, [inivars; finvars])

    inicnt = length(inivars)
    inisvars = collect(Iterators.take(svars, inicnt))
    finsvars = collect(Iterators.drop(svars, inicnt))

    basis = [ x - y for (x, y) in zip(inisvars, finsvars)]
    Singular.Ideal(R, basis)
end

function replace_functions(expr::Sym, index::Int)
    w0   = Wild("w0")
    fns  = Sym.(collect(expr.atoms(AppliedUndef)))
    vars = [Sym("$(string(Sym(fn.func)))_$(convert(Int, fn.args[1])+index)") for fn in fns]
    dict = Dict(zip(fns, vars))
    return subs(expr, dict)
end


text:
using Workspace
using Test
using Configurations
using StructTypes
using JSON

d = Dict(
    "name" => "test_workspace",
    "uuid" => "94f78cb8-c85a-11eb-0dc6-f7f6bd9facf1",
    "compiler" => Dict(
        "project" => Dict(
            "deps" => Dict(
                "JSON" => "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
            ),
            "compat" => Dict(
                "JSON" => "0.20.0",
            )
        )
    )
)

s = JSON.json(d)
ins = from_dict(WorkspaceInstance, JSON.parse(s))
table = WorkspaceTable()
add_workspace_process!(table, ins)
eval_in_workspace(table, ins, quote
    using Pkg
    Pkg.status()
end)

using Distributed
Distributed.remotecall_eval(Main, 3, Expr(:toplevel, quote
using JSON
end))