text:
#
# Tests depend on a running server at localhost:27017,
# and will create a database named "mongoc".
#

import Mongoc

if VERSION < v"0.7-"
    using Base.Test
else
    using Test
    using Dates
end

const DB_NAME = "mongoc"

#=
function gc_on_osx_v6()
    @static if VERSION < v"0.7-" && is_apple()
            gc()
    else
        nothing
    end
end
=#

@testset "BSON" begin

    @testset "as_json" begin
        @test_throws ErrorException Mongoc.BSON(""" { ajskdla sdjsafd } """)
        bson = Mongoc.BSON("""{"hey" : 1}""")
        @test bson["hey"] == 1
        @test Mongoc.as_json(bson) == """{ "hey" : 1 }"""
        @test Mongoc.as_json(bson, canonical=true) == """{ "hey" : { "\$numberInt" : "1" } }"""
    end

    @testset "oid conversion" begin
        local oid::Mongoc.BSONObjectId = "5b9fb22b3192e3fa155693a1"
        local oid_str::String = oid
    end

    @testset "oid compare" begin
        x = Mongoc.BSONObjectId()
        y = Mongoc.BSONObjectId()
        @test x != y
        @test x == x
        @test y == y
        @test hash(x) == hash(x)
        @test hash(y) == hash(y)
        @test hash(x) != hash(y)
        @test Mongoc.bson_oid_compare(x, y) < 0
        @test Mongoc.bson_oid_compare(y, x) > 0
        @test Mongoc.bson_oid_compare(x, x) == 0
        @test Mongoc.bson_oid_compare(y, y) == 0
    end

    @testset "oid string" begin
        x = Mongoc.BSONObjectId()
        y = Mongoc.BSONObjectId(string(x))
        @test x == y

        @test_throws ErrorException Mongoc.BSONObjectId("invalid_objectid")
    end

    @testset "oid time" begin
        Mongoc.get_time( Mongoc.BSONObjectId("5b9eaa2711c3dd0d6a46a5c4") ) == DateTime(2018, 9, 16, 19, 8, 23) # 2018-09-16T19:08:23
    end

    # https://github.com/JuliaLang/julia/issues/29193
    #=
    if VERSION < v"0.7-"
        @testset "BSONObjectId segfault issue" begin
            io = IOBuffer()
            v = Vector{Mongoc.BSONObjectId}()

            for i in 1:5
                push!(v, Mongoc.BSONObjectId())
            end
            show(io, v)
        end

        @testset "oid sort" begin
            v = Vector{Mongoc.BSONObjectId}()
            for i in 1:10_000
                push!(v, Mongoc.BSONObjectId())
            end
            @test length(v) == length(unique(v))

            v_sorted = sort(v, lt = (a,b) -> Mongoc.bson_oid_compare(a,b) < 0)
            @test length(v_sorted) == length(v)
            for i in 1:length(v_sorted)
                @test v_sorted[i] == v[i]
            end
        end
    end
    =#

    @testset "BSON Iterator" begin
        doc = Mongoc.BSON("""{ "a" : 1, "b" : 2.2, "str" : "my string", "bool_t" : true, "bool_f" : false, "array" : [1, 2, false, "inner_string"], "document" : { "a" : 1, "b" : "b_string"}  }""")

        new_id = Mongoc.BSONObjectId()
        doc["_id"] = new_id

        @test haskey(doc, "a")
        @test haskey(doc, "b")
        @test !haskey(doc, "c")
        @test haskey(doc, "str")
        @test haskey(doc, "_id")
        @test haskey(doc, "array")
        @test haskey(doc, "document")

        @test doc["a"] == 1
        @test doc["b"] == 2.2
        @test doc["str"] == "my string"
        @test doc["bool_t"]
        @test !doc["bool_f"]
        @test doc["_id"] == new_id
        @test doc["array"] == [1, 2, false, "inner_string"]
        @test doc["document"] == Dict("a"=>1, "b"=>"b_string")

        @test_throws ErrorException doc["invalid key"]

        doc_dict = Mongoc.as_dict(doc)
        @test doc_dict["a"] == 1
        @test doc_dict["b"] == 2.2
        @test doc_dict["str"] == "my string"
        @test doc_dict["bool_t"]
        @test !doc_dict["bool_f"]
        @test doc_dict["_id"] == new_id
        @test doc_dict["array"] == [1, 2, false, "inner_string"]
        @test doc_dict["document"] == Dict("a"=>1, "b"=>"b_string")
    end

    @testset "BSON write" begin
        bson = Mongoc.BSON()
        new_oid = Mongoc.BSONObjectId()
        bson["_id"] = new_oid
        bson["int32"] = Int32(10)
        bson["int64"] = Int64(20)
        bson["string"] = "hey you"
        bson["bool_true"] = true
        bson["bool_false"] = false
        bson["double"] = 2.3
        bson["datetime"] = DateTime(2018, 2, 1, 10, 20, 35, 10)
        bson["vector"] = collect(1:10)
        bson["source"] = Mongoc.BSONCode("function() = 1")
        @test_throws ErrorException bson["key"] = Date(2018, 9, 18)

        let
            sub_bson = Mongoc.BSON()
            sub_bson["hey"] = "you"
            sub_bson["num"] = 10
            bson["sub_document"] = sub_bson
        end

        @test bson["_id"] == new_oid
        @test bson["int32"] == Int32(10)
        @test bson["int64"] == Int64(20)
        @test bson["string"] == "hey you"
        @test bson["bool_true"]
        @test !bson["bool_false"]
        @test bson["double"] == 2.3
        @test bson["datetime"] == DateTime(2018, 2, 1, 10, 20, 35, 10)
        @test bson["sub_document"]["hey"] == "you"
        @test bson["sub_document"]["num"] == 10
        @test bson["vector"] == collect(1:10)
        @test bson["source"] == Mongoc.BSONCode("function() = 1")

        let
            sub_bson = bson["sub_document"]
            @test sub_bson["hey"] == "you"
            @test sub_bson["num"] == 10
        end
    end

    @testset "BSON to Dict conversion" begin
        dict = Dict("a" => 1, "b" => false, "c" => "string")
        doc = Mongoc.BSON(dict)
        @test doc["a"] == 1
        @test doc["b"] == false
        @test doc["c"] == "string"
        @test dict == Mongoc.as_dict(doc)
    end
end

@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 Mongoc.find(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

        #gc_on_osx_v6() # avoid segfault on Cursor destroy

        @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 "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

        #gc_on_osx_v6() # avoid segfault on Cursor destroy
    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

    #gc_on_osx_v6()

    @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"
        @static if VERSION < v"0.7"
            warn("MongoDB server version $server_version does not support Sessions. Skipping tests.")
        else
            @warn("MongoDB server version $server_version does not support Sessions. Skipping tests.")
        end
    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 DrWatson
@quickactivate
using ArgParse
using GenerativeAD
import StatsBase: fit!, predict
using StatsBase
using BSON
using Flux
using IPMeasures
using GenerativeModels
using Distributions

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
        arg_type = Int
        help = "seed"
        default = 1
    "dataset"
        default = "iris"
        arg_type = String
        help = "dataset"
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, contamination = parsed_args

#######################################################################################
################ THIS PART IS TO BE PROVIDED FOR EACH MODEL SEPARATELY ################
modelname = "wae"
# sample parameters, should return a Dict of model kwargs
"""
	sample_params()

Should return a named tuple that contains a sample of model parameters.
"""
function sample_params()
	par_vec = (2 .^(3:8), 2 .^(4:9), 10f0 .^(-4:-3), 2 .^ (5:7), ["relu", "swish", "tanh"], 3:4, 1:Int(1e8),
		["imq", "gauss", "rq"], 10f0 .^ (-3:0), 10f0 .^(-1:0))
	argnames = (:zdim, :hdim, :lr, :batchsize, :activation, :nlayers, :init_seed, :kernel, :sigma, :lambda)
	parameters = (;zip(argnames, map(x->sample(x, 1)[1], par_vec))...)
	# ensure that zdim < hdim
	while parameters.zdim >= parameters.hdim
		parameters = merge(parameters, (zdim = sample(par_vec[1])[1],))
	end
	return parameters
end
"""
	fit(data, parameters)

This is the most important function - returns `training_info` and a tuple or a vector of tuples `(score_fun, final_parameters)`.
`training_info` contains additional information on the training process that should be saved, the same for all anomaly score functions.
Each element of the return vector contains a specific anomaly score function - there can be multiple for each trained model.
Final parameters is a named tuple of names and parameter values that are used for creation of the savefile name.
"""
function fit(data, parameters)
	# construct model - constructor should only accept kwargs
	model = GenerativeAD.Models.vae_constructor(;idim=size(data[1][1],1), parameters...)

	# construct loss function
	if parameters.kernel == "imq"
		k = IMQKernel(parameters.sigma)
	elseif parameters.kernel == "gauss"
		k = GaussianKernel(parameters.sigma)
	elseif parameters.kernel == "rq"
		k = RQKernel(parameters.sigma)
	else
		error("given kernel not known")
	end
	loss(m::GenerativeModels.VAE,x) = parameters.lambda*mmd_mean(m, x, k) .- mean(logpdf(m.decoder, x, rand(m.encoder, x)))
	loss(m::GenerativeModels.VAE, x, batchsize::Int) =
		mean(map(y->loss(m,y), Flux.Data.DataLoader(x, batchsize=batchsize)))

	# fit train data
	try
		global info, fit_t, _, _, _ = @timed fit!(model, data, loss; max_train_time=82800/max_seed,
			patience=200, check_interval=10, parameters...)
	catch e
		# return an empty array if fit fails so nothing is computed
		@info "Failed training due to \n$e"
		return (fit_t = NaN, history=nothing, npars=nothing, model=nothing), []
	end

	# construct return information - put e.g. the model structure here for generative models
	training_info = (
		fit_t = fit_t,
		history = info.history,
		npars = info.npars,
		model = info.model
		)

	# now return the different scoring functions
	training_info, [
		(x -> GenerativeAD.Models.reconstruction_score(info.model, x), merge(parameters, (score = "reconstruction",))),
		(x -> GenerativeAD.Models.reconstruction_score_mean(info.model, x), merge(parameters, (score = "reconstruction-mean",))),
		(x -> GenerativeAD.Models.latent_score(info.model, x), merge(parameters, (score = "latent",))),
		(x -> GenerativeAD.Models.latent_score_mean(info.model, x), merge(parameters, (score = "latent-mean",))),
		]
end
function GenerativeAD.edit_params(data, parameters)
	idim = size(data[1][1],1)
	# put the largest possible zdim where zdim < idim, the model tends to converge poorly if the latent dim is larger than idim
	if parameters.zdim >= idim
		zdims = 2 .^(1:8)
		zdim_new = zdims[zdims .< idim][end]
		parameters = merge(parameters, (zdim=zdim_new,))
	end
	parameters
end

####################################################################
################ THIS PART IS COMMON FOR ALL MODELS ################
# only execute this if run directly - so it can be included in other files
if abspath(PROGRAM_FILE) == @__FILE__
	# set a maximum for parameter sampling retries
	try_counter = 0
	max_tries = 10*max_seed
	cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
	while try_counter < max_tries
	    parameters = sample_params()

	    for seed in 1:max_seed
			savepath = datadir("experiments/tabular$cont_string/$(modelname)/$(dataset)/seed=$(seed)")
			mkpath(savepath)

			# get data
			data = GenerativeAD.load_data(dataset, seed=seed, contamination=contamination)

			# edit parameters
			edited_parameters = GenerativeAD.edit_params(data, parameters)

			@info "Trying to fit $modelname on $dataset with parameters $(edited_parameters)..."
			@info "Train/validation/test splits: $(size(data[1][1], 2)) | $(size(data[2][1], 2)) | $(size(data[3][1], 2))"
			@info "Number of features: $(size(data[1][1], 1))"

			# check if a combination of parameters and seed alread exists
			if GenerativeAD.check_params(savepath, edited_parameters)
				# fit
				training_info, results = fit(data, edited_parameters)

				# save the model separately
				if training_info.model != nothing
					tagsave(joinpath(savepath, savename("model", edited_parameters, "bson", digits=5)),
						Dict("model"=>training_info.model,
							"fit_t"=>training_info.fit_t,
							"history"=>training_info.history,
							"parameters"=>edited_parameters
							), safe = true)
					training_info = merge(training_info, (model = nothing,))
				end

				# here define what additional info should be saved together with parameters, scores, labels and predict times
				save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, contamination = contamination))

				# now loop over all anomaly score funs
				for result in results
					GenerativeAD.experiment(result..., data, savepath; save_entries...)
				end
				global try_counter = max_tries + 1
			else
				@info "Model already present, trying new hyperparameters..."
				global try_counter += 1
			end
		end
	end
	(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing
end

ext: jl


text:
"""
    euler_method(f, x0, span, h=1.0e-2)

Calculate the solution to a differential equation using forward euler method.
"""
function euler_method(f, x0, span, h = 1.0e-2)
    s, e = span
    steps = Int(floor((e - s) / h))
    x = Vector{Float64}(undef, steps)
    t = Vector{Float64}(undef, steps)
    x[1] = x0
    t[1] = 0
    for i in 1:(steps-1)
        x[i+1] = x[i] + h * f(x[i], t[i])
        t[i+1] = t[i] + h
    end
    return x, t
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
# This file contains algorithms needed for the integratiorn of
# rational functions from chapter 2 of the book
#
#    Manuel Bronstein, Symbolic Integration I: Transcendental Functions,
#    2nd ed., Springer 2005.
#

function HermiteReduce_original(A::PolyElem{T}, D::PolyElem{T}) where T <: FieldElement
    # See Bronstein's book, Section 2.2, p. 40
    Ds = Squarefree(D)
    n = length(Ds)
    P, As = PartialFraction(A, [Ds[j]^j for j=1:n])
    g = zero(D)//one(D) # rational function with value 0
    h = P + As[1]//Ds[1] # rational function
    for k=2:n
        if degree(Ds[k])>0
            V = Ds[k]
            for j=k-1:-1:1
                B, C = gcdx(derivative(V), V, -divexact(As[k], j))
                g += B//V^j
                As[k] = -j*C - derivative(B)
            end
            h += As[k]//V
        end
    end
    g, h
end

function HermiteReduce_quadratic(A::PolyElem{T}, D::PolyElem{T}) where T <: FieldElement
    # See Bronstein's book, Section 2.2, p. 41
    g = zero(D)//one(D) # rational function with value 0
    Ds = Squarefree(D)
    m = length(Ds)
    for i=2:m
        if degree(Ds[i])>0
            V = Ds[i]
            U = divexact(D, V^i)
            for j=i-1:-1:1
                B, C = gcdx(U*derivative(V), V, -divexact(A, j))
                g += B//V^j
                A = -j*C - U*derivative(B)
            end
            D = U*V
        end
    end
    g, A//D
end

function HermiteReduce(A::PolyElem{T}, D::PolyElem{T}) where T <: FieldElement
    # See Bronstein's book, Section 2.2, p. 44
    g = zero(D)//one(D) # rational function with value 0
    Dminus = gcd(D, derivative(D))
    Dstar = divexact(D, Dminus)
    while degree(Dminus)>0
        Dminus2 = gcd(Dminus, derivative(Dminus))
        Dminusstar = divexact(Dminus, Dminus2)
        B, C = gcdx(-divexact(Dstar*derivative(Dminus), Dminus), Dminusstar, A)
        A = C - divexact(derivative(B)*Dstar, Dminusstar)
        g += B//Dminus
        Dminus = Dminus2
    end
    return g, A//Dstar
end

struct SumOfLogTerms{T<:FieldElement, P<:PolyElem{T}} <: Term
    R::P
    S::PolyElem{P}
end

function Base.show(io::IO, t::SumOfLogTerms)
    s = string(parent(t.R).S)
    print(io, "Σ_{$s: ")
    show(io, t.R)
    print(io, " = 0}$s*log(")
    show(io, t.S)
    print(io, ")")
end

function IntRationalLogPart(A::PolyElem{T}, D::PolyElem{T}; make_monic::Bool=false, symbol=:α) where T <: FieldElement
    # See Bronstein's book, Section 2.5, p. 51
    F = base_ring(A)
    Ft, t = PolynomialRing(F, symbol)
    FtX, X = PolynomialRing(Ft, symbols(parent(A))[1])
    R, Rs = SubResultant(D(X), A(X)-t*derivative(D)(X))
    Qs = Squarefree(R)
    ds = degree.(Qs)
    Qs = [f for f in Qs if degree(f)>0]
    Ss = typeof(X)[]
    ii = 1
    for i=1:length(ds)
        if ds[i]>0
            if i==degree(D)
                S = D(X)
            else
                m = findfirst(z->z==i, degree.(Rs))
                S =  Rs[m]
                As = Squarefree(leading_coefficient(S))
                for j=1:length(As)
                    S = divexact(S, gcd(As[j],Qs[ii])^j)
                end
            end
            if make_monic
                Qs[ii] = divexact(Qs[ii], leading_coefficient(Qs[ii]))
                h = invmod(leading_coefficient(S), Qs[ii])
                S = FtX([rem(h*c, Qs[ii]) for c in coefficients(S)])
            else
                S = FtX([rem(c, Qs[ii]) for c in coefficients(S)])
            end
            push!(Ss, S)
            ii+=1
        end
    end
    [SumOfLogTerms(Q, S) for (Q, S) in zip(Qs, Ss)]
end

function Complexify(R::PolyElem{T}; symbols=[:α, :β]) where T <: FieldElement
    F = base_ring(R)
    Fuv, uv = PolynomialRing(F, symbols)
    u = uv[1]
    v = uv[2]
    c = collect(coefficients(R))
    d = length(c)-1
    P = zero(Fuv)
    for n=0:d
        for k=0:2:n
            P += (-1)^div(k,2)*binomial(n,k)*c[n+1]*u^(n-k)*v^k
        end
    end
    Q = zero(Fuv)
    for n=1:d
        for k=1:2:n
            Q += (-1)^div(k,2)*binomial(n,k)*c[n+1]*u^(n-k)*v^k
        end
    end
    P, Q
end

function LogToAtan(A::PolyElem{T}, B::PolyElem{T}) where T <: FieldElement
    # See Bronstein's book, Section 2.8, p. 63
    Q, R = divrem(A, B)
    if iszero(R)
        return [FunctionTerm(atan, 2, Q)]
    end
    if degree(A)<degree(B)
        return LogToAtan(-B, A)
    end
    G, D, C = gcdx(B, -A)
    return vcat(FunctionTerm(atan, 2, divexact(A*D+B*C, G)), LogToAtan(D, C))
end

struct SumOfRealTerms <: Term
    R
    S
    P
    Q
    LT::Term # log term
    ATs::Vector{Term} #atan terms
end

function Base.show(io::IO, t::SumOfRealTerms)
    s = string.(parent(t.P).S)
    print(io, "Σ_{$(s[1]),$(s[2]): $(s[2])>0, ")
    show(io, t.P)
    print(io, " = ")
    show(io, t.Q)
    print(io, " = 0}")
    show(io, t.LT)
    print(io, " + $(s[2])*(")
    print(io, join(string.(t.ATs), "+"))
    print(io, "))")
    s0 = string(parent(t.R).S)
    print(io, " + Σ_{$s0: Im($s0)=0, ")
    show(io, t.R)
    print(io, " = 0}$s0*log(")
    show(io, t.S)
    print(io, ")")

end

function LogToReal(t::SumOfLogTerms; symbols=[:α, :β]) #{T, PP}) where {T<:FieldElement, PP<:PolyElem{T}}
    # See Bronstein's book, Section 2.8, p. 69
    F = base_ring(t.R)
    R, uv = PolynomialRing(F, symbols)
    K = FractionField(R)
    Kx, x = PolynomialRing(K, "x")
    P, Q = Complexify(t.R)
    cc =[Complexify(c) for c in coefficients(t.S)]
    A = Kx([c[1] for c in cc])
    B = Kx([c[2] for c in cc])
    SumOfRealTerms(t.R, t.S, P, Q, FunctionTerm(log, 1, A^2+B^2), LogToAtan(A, B))
end

function positive_constant_coefficient(f::PolyElem)
    if constant_coefficient(f)<0
        return -f
    else
        return f
    end
end

function rationalize_if_possible(x::qqbar) #Nemo algebraic number type
    if degree(x)==1
        return fmpq(x)
    else
        return x
    end
end

function rationalize_if_possible(f::PolyElem{qqbar})
    if maximum(degree.(coefficients(f)))==1
        return polynomial(Nemo.QQ, fmpq.(coefficients(f)))
    else
        return f
    end
end

function Eval(t::SumOfLogTerms; real_output::Bool=true)
    var = string(symbols(parent(t.S))[1])
    F = base_ring(t.R)
    as = roots(t.R)
    if length(as)==degree(t.R)
        return [FunctionTerm(log, a, positive_constant_coefficient(
            polynomial(F, [c(a) for c in coefficients(t.S)], var))) for a in as]
    end

    as = roots(t.R, QQBar)
    us = real.(as)
    vs = imag.(as)
    if iszero(vs) || !real_output
        return [FunctionTerm(log, rationalize_if_possible(a), rationalize_if_possible(positive_constant_coefficient(
            polynomial(QQBar, [c(a) for c in coefficients(t.S)], var)))) for a in as]
    end
    r = LogToReal(t)
    result = Term[]
    for i = 1:length(as)
        a = as[i]
        u, v = us[i], vs[i]
        if iszero(v)
            push!(result, FunctionTerm(log, rationalize_if_possible(u), rationalize_if_possible(positive_constant_coefficient(
            polynomial(QQBar, [c(u) for c in coefficients(r.S)], var)))))
        elseif v>0
            if !iszero(u)
                push!(result, FunctionTerm(log, rationalize_if_possible(r.LT.coeff*u), rationalize_if_possible(positive_constant_coefficient(
                polynomial(QQBar, [numerator(c)(u, v)//denominator(c)(u, v) for c in coefficients(r.LT.arg)], var)))))
            end
            for AT in r.ATs
                push!(result, FunctionTerm(atan, rationalize_if_possible(AT.coeff*v), rationalize_if_possible(
                polynomial(QQBar, [numerator(c)(u, v)//denominator(c)(u, v) for c in coefficients(AT.arg)], var))))
            end
        end
    end
    result
end

function IntegrateRationalFunction(f::FracElem{P}) where {T<:FieldElement, P<:PolyElem{T}}
    # See Bronstein's book, Section 2.5, p. 52
    g, h = HermiteReduce(numerator(f), denominator(f))
    Q, R = divrem(numerator(h), denominator(h))
    result = Term[]
    IQ = integral(Q)
    if !iszero(IQ)
        push!(result, IdTerm(IQ))
    end
    if !iszero(g)
        push!(result, IdTerm(g))
    end
    if !iszero(R)
        SL =  IntRationalLogPart(R, denominator(h), make_monic=true, symbol=:α)
        result = vcat(result, [Eval(h, real_output=true) for h in SL]...)
    end
    result
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using ACE, Test, Printf, LinearAlgebra, StaticArrays, BenchmarkTools,
      JuLIP, JuLIP.Testing

##
@testset "ACE.jl" begin
    # ------------------------------------------
    #   basic polynomial basis building blocks
    include("polynomials/test_ylm.jl")
    include("polynomials/test_rylm.jl")
    include("polynomials/test_transforms.jl")
    include("polynomials/test_orthpolys.jl")

    # --------------------------------------------
    # core permutation-invariant functionality
    include("test_1pbasis.jl")
    include("test_pibasis.jl")
    include("test_pipot.jl")

    # ------------------------
    #   rotation_invariance
    include("rpi/test_cg.jl")
    include("rpi/test_rpibasis.jl")

    # ----------------------
    #   pair potentials
    include("pair/test_pair_basis.jl")
    include("pair/test_pair_pot.jl")
    include("pair/test_repulsion.jl")

    # ----------------------
    #   miscallaneous ...
    include("compat/test_compat_v05.jl")
    include("compat/test_compat.jl")
    include("test_any.jl")

    # ----------------------------------
    #    old tests to be re-introduced
    # include("test_real.jl")
    # include("test_orth.jl")
    # include("test_descriptor.jl")
    # include("bonds/test_cylindrical.jl")
    # include("bonds/test_fourier.jl")
    # include("bonds/test_envpairbasis.jl")
end


text:
#######################################
#
# Plot spectrum of signal
#
#######################################

function plot_spectrum(signal::Vector, fs::Real;
titletext::S="", Fmin::Number=0, Fmax::Number=90, targetFreq::F=0.0, dBPlot::Bool=true,
noise_level::Number=0, signal_level::Number=0) where {S <: AbstractString, F <: AbstractFloat}

    # Determine fft frequencies
    signal_length = length(signal)
    frequencies = range(0, stop=1, length=Int(signal_length / 2 + 1))*fs/2

    # Calculate fft and convert to power
    fftSweep = 2 / signal_length * fft(signal)
    spectrum = abs.(fftSweep[1:div(signal_length, 2) + 1])  # Amplitude
    spectrum = spectrum.^2

    valid_idx = (frequencies .<= Fmax) .& (frequencies .>= Fmin)
    spectrum = spectrum[valid_idx]
    frequencies = frequencies[valid_idx]

    # Want a log plot?
    if dBPlot
        spectrum = 10*log10.( spectrum )
        ylabel="Response Power (dB)"
    else
        ylabel="Response Power (uV^2)"
    end

    # Plot signal
    Plots.plot(frequencies, spectrum, lab = "Spectrum", color = :black)
    xlims!(Fmin, Fmax)
    xlabel!("Frequency (Hz)")
    ylabel!(ylabel)
    p = title!(titletext)

    # Plot the noise level if requested
    if noise_level != 0
        if dBPlot
            noise_level = 10*log10.(noise_level)
        end
        p = plot!([Fmin, targetFreq+2], [noise_level, noise_level], lab = "Noise", color = :red)
    end

    # Plot the signal level if requested
    if signal_level != 0
        if dBPlot
            signal_level = 10*log10.(signal_level)
        end
        plot!([Fmin, targetFreq], [signal_level, signal_level], lab = "Signal", color = :green)

        targetFreqIdx =   something(findfirst(isequal(minimum(abs.(frequencies.-targetFreq))), abs.(frequencies.-targetFreq)), 0)
        targetFreq    = frequencies[targetFreqIdx]
        targetResults = spectrum[targetFreqIdx]
        #TODO remove label for circle rather than make it empty
        p = plot!([targetFreq], [targetResults], marker = (:circle, 5, 0.1, :green), markerstrokecolor = :green, lab = "")
    end
    return p
end

function plot_spectrum(eeg::SSR, chan::Int; targetFreq::Number=0)

    channel_name = channelnames(eeg)[1]

    # Check through the processing to see if we have done a statistical test at target frequency
    signal = nothing
    result_idx = find_keys_containing(eeg.processing, "statistics")

    for r = 1:length(result_idx)
        result = get(eeg.processing, collect(keys(eeg.processing))[result_idx[r]], 0)
        if result[:AnalysisFrequency][1] == targetFreq

            result_snr = result[:SNRdB][chan]
            signal = result[:SignalAmplitude][chan]^2
            noise  = result[:NoiseAmplitude][chan]^2
            title  = "Channel $(channel_name). SNR = $(round(result_snr, sigdigits=4)) dB"
        end
    end

    if signal == nothing
        title  = "Channel $(channel_name)"
        noise  = 0
        signal = 0
    end

    title = replace(title, "_", " ")

    avg_sweep = squeeze(mean(eeg.processing["sweeps"], 2), 2)
    avg_sweep = avg_sweep[:, chan]
    avg_sweep = convert(Array{Float64}, vec(avg_sweep ))

    p = plot_spectrum(avg_sweep, eeg.header["sampRate"][1];
                        titletext=title, targetFreq = targetFreq,
                        noise_level = noise, signal_level = signal)

    return p
end

function plot_spectrum(eeg::SSR, chan::AbstractString; targetFreq::Number=modulationrate(eeg))

    return plot_spectrum(eeg,  something(findfirst(isequal(chan), channelnames(eeg)), 0), targetFreq=targetFreq)
end

#######################################
#
# Plot time series
#
# Automatic plotting for single channel (with units) or multichannel (with channel names) time series
#
#######################################

"""
Plot a single channel time series

#### Input

* signal: Vector of data
* fs: Sample rate
* channels: Name of channel to plot
* plot_points: Number of points to plot, they will be equally spread. Used to speed up plotting
* Other optional arguements are passed to gadfly plot function

#### Output

Returns a figure

"""
function plot_single_channel_timeseries(signal::AbstractVector{T}, fs::Real;
        xlabel::S="Time (s)", ylabel::S="Amplitude (uV)", lab::S="", kwargs...) where {T <: Number, S <: AbstractString}

    @debug("Plotting single channel waveform of size $(size(signal))")

    time_s = collect(1:size(signal, 1))/fs   # Create time axis

    Plots.plot(time_s, signal, t=:line, c=:black, lab=lab, xlabel=xlabel, ylabel=ylabel)
end

"""
Plot a multi channel time series

#### Input

* signals: Array of data
* fs: Sample rate
* channels: Name of channels
* plot_points: Number of points to plot, they will be equally spread. Used to speed up plotting
* Other optional arguements are passed to gadfly plot function

#### Output

Returns a figure

"""
function plot_multi_channel_timeseries(signals::Array{T, 2}, fs::Number, channels::Array{S};
        xlabel::S="Time (s)", ylabel::S="Amplitude (uV)", kwargs...) where {T <: Number, S <: AbstractString}

    @debug("Plotting multi channel waveform of size $(size(signals))")

    time_s = collect(1:size(signals, 1))/fs                        # Create time axis

    variances = var(signals, dims = 1)          # Variance of each figure for rescaling
    mean_variance = mean(variances)     # Calculate for rescaling figures

    p = Plots.plot(t=:line, c=:black, xlabel=xlabel, ylabel=ylabel, ylim=(-0.5, size(signals, 2) - 0.5))

    for c in 1:size(signals, 2)                                  # Plot each channel
        signals[:, c] = signals[:, c] - mean(signals[:, c])      # Remove mean
        signals[:, c] = signals[:, c] / mean_variance .+ (c-1)   # Rescale and shift so all chanels are visible
        p = plot!(time_s, signals[:, c], c=:black, lab="")
    end
    p = plot!(yticks = (0:length(channels) - 1, channels))

    return p
end

#######################################
#
# Filter response
#
#######################################

# Plot filter response
function plot_filter_response(zpk_filter::FilterCoefficients, fs::Integer;
              lower::Number=1, upper::Number=30, sample_points::Int=1024)

    frequencies = range(lower, stop=upper, length=1024)
    h = freqz(zpk_filter, frequencies, fs)
    magnitude_dB = 20*log10.(convert(Array{Float64}, abs.(h)))
    phase_response = (360/(2*pi))*unwrap(convert(Array{Float64}, angle.(h)))

    p1 = plot(frequencies, magnitude_dB,   lab = "")
    p2 = plot(frequencies, phase_response, lab = "")

    p = plot(p1,p2, ylabel = ["Magnitude (dB)" "Phase (degrees)"], xlabel = "Frequency (Hz)", layout = @layout([a; b]))
end


text:
# Elias Kuepper, Holger Teichgraeber, 2019

 ######################
 # CapacityExpansion
 # Capacity Expansion Probelem formulation
 #
 #####################
using TimeSeriesClustering
using CSV
using YAML
using DataFrames
using StatsBase
using JLD2
using FileIO
using JuMP

include(joinpath("utils","datastructs.jl"))
include(joinpath("utils","optvariable.jl"))
include(joinpath("utils","utils.jl"))
include(joinpath("utils","load_data.jl"))
include(joinpath("optim_problems","run_opt.jl"))
include(joinpath("optim_problems","set.jl"))
include(joinpath("optim_problems","opt.jl"))


text:
using Plots
using LinearAlgebra
using Statistics
using Random
using Distributions
using PDMats
using PDMatsExtras
using ShiftedArrays
using Parameters
using Polynomials

using DelimitedFiles
using DataFrames

# Import data
data = readdlm("test/pca/real_data/data_deseasoned.csv", ',')
regstatematch = readdlm("test/otrok_whiteman/real_data/regstatematch.csv", ',')

# Data as matrix
datamat = Matrix(data[2:end, :])     # remove first row containing col names
datamat = Float64.(datamat)         # transform entires to Float64

for i in 1:size(datamat)[2]
    datamat[:, i] = (datamat[:, i] .- mean(datamat[:, i])) ./ std(datamat[:, i])
end

regstatematch = Matrix(regstatematch[2:end, :])
regassign = Int64.(regstatematch[:, 3])

# Estimate HDFM

nlevels = 2

nvar = size(datamat)[2]

nfactors = [1, max(regassign...)]

fassign = ones(Int, nvar, nlevels)
fassign[:, 2] = regassign

flags = [3, 3]

varlags = 3 * ones(Int, nvar)

hdfmpriors = HDFMStruct(nlevels = nlevels,
    nfactors = nfactors,
    factorassign = fassign,
    factorlags = flags,
    errorlags = varlags,
    ndraws = 1000,
    burnin = 50)

#results = PCA2LevelEstimator(datamat, hdfmpriors)
resultsKN = KN2LevelEstimator(datamat, hdfmpriors)
resultsOW = OW2LevelEstimator(datamat, hdfmpriors)
resultsPCA = PCA2LevelEstimator(datamat, hdfmpriors)

#Save state names
statenames = data[1, :]

medians = Any[]
quant33 = Any[]
quant66 = Any[]
stds = Any[]

j = 1
for i in 1:size(results.F)[1]
    push!(stds, std(results.F[i, j, :]))
    push!(quant33, quantile(results.F[i, j, :], 0.33))
    push!(quant66, quantile(results.F[i, j, :], 0.66))
    push!(medians, median(results.F[i, j, :]))
end

plot(results.means.F[:, j][1:end])
#plot!(results.means.F[:, 1])
plot(medians)
plot!(quant33)
plot!(quant66)

vardecomp = vardecomp2level(datamat, results.means.F, reshape(results.means.B, 3, 50)', fassign)
#vardecomp2 = vardecomp2level(datamat, results2.means.F, reshape(results2.means.B, 3, 50)', fassign)
vardecomp = [vardecomp zeros(nvar)]
vardecomp[:, 3] = ones(nvar) - vardecomp[:, 1] - vardecomp[:, 2]

groupedbar(vardecomp,
    bar_position = :stack,
    bar_width = 0.9,
    xticks = (1:nvar, 1:nvar),
    label = ["National" "Regional" "Idiosyncratic"],
    title = "Variance Decompositions",
    ylabel = "Variance Contribution",
    xlabel = "State")

plot(vardecomp[:, 1])
plot!(vardecomp[:, 2])

plot(ones(size(datamat)[2]) - vardecomp[:, 1] - vardecomp[:, 2])

histogram(vardecomp[:, 1], normalize = :probability)
histogram(vardecomp[:, 2], normalize = :probability)


text:
for orientation in (:row, :col)

    row_oriented = orientation == :row
    dim_obs, dim_param = row_oriented ? (1, 2) : (2, 1)
    NT, TN = row_oriented ? ('N', 'T') : ('T', 'N')

    @eval begin

        function dotvectors!(
                ::Val{$(Meta.quot(orientation))},
                xᵀx::Vector{T},
                X::Matrix{T}
            ) where {T<:AbstractFloat}
            if !(size(X, $dim_obs) == length(xᵀx))
                errorstring = string("Dimension mismatch on dimension ", $dim_obs)
                throw(DimensionMismatch(errorstring))
            end
            fill!(xᵀx, zero(T))
            for I in Base.Cartesian.CartesianIndices(size(X))
                xᵀx[I.I[$dim_obs]] += X[I]^2
            end
            xᵀx
        end

        @inline function dotvectors(σ::Val{$(Meta.quot(orientation))}, X::Matrix{T}) where {T<:AbstractFloat}
            dotvectors!(σ, Array{T}(undef, size(X,$dim_obs)), X)
        end

        function gramian!(
                ::Val{$(Meta.quot(orientation))},
                G::Matrix{T},
                X::Matrix{T},
                symmetrize::Bool
            ) where {T<:LinearAlgebra.BLAS.BlasReal}
            LinearAlgebra.BLAS.syrk!('U', $NT, one(T), X, zero(T), G)
            symmetrize ? LinearAlgebra.copytri!(G, 'U') : G
        end

        @inline function gramian!(
                ::Val{$(Meta.quot(orientation))},
                G::Matrix{T},
                X::Matrix{T},
                Y::Matrix{T}
            ) where {T<:LinearAlgebra.BLAS.BlasReal}
            LinearAlgebra.BLAS.gemm!($NT, $TN, one(T), X, Y, zero(T), G)
        end
    end
end


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

JLLWrappers.@generate_main_file_header("Libbpf")
JLLWrappers.@generate_main_file("Libbpf", UUID("02f9a84d-9555-5f1a-8c3c-42027521038d"))
end  # module Libbpf_jll


text:
export TinkerforgeDevice
abstract type TinkerforgeDevice end

# returned by get_connection_state
export TinkerforgeIPConConnectionState, CONNECTION_STATE_DISCONNECTED, CONNECTION_STATE_CONNECTED, CONNECTION_STATE_PENDING
@enum TinkerforgeIPConConnectionState begin
    CONNECTION_STATE_DISCONNECTED = 0
    CONNECTION_STATE_CONNECTED = 1
    CONNECTION_STATE_PENDING = 2 # auto-reconnect in process
end

export IPConnection
mutable struct IPConnection
	ipconInternal::PyObject

  function IPConnection()
    package = pyimport("tinkerforge.ip_connection")
    ipconInternal = package.IPConnection()

    return new(ipconInternal)
  end
end

export connect
"""
    $(SIGNATURES)

Creates a TCP/IP connection to the given *host* and *port*. The host
and port can point to a Brick Daemon or to a WIFI/Ethernet Extension.
Devices can only be controlled when the connection was established
successfully.
Blocks until the connection is established and throws an exception if
there is no Brick Daemon or WIFI/Ethernet Extension listening at the
given host and port.
"""
connect(ipcon::IPConnection, host::IPAddr, port::Integer) = connect(ipcon, string(host), port)
connect(ipcon::IPConnection, host::String, port::Integer) = ipcon.ipconInternal.connect(host, port)
connect(ipcon::IPConnection) = connect(ipcon, "localhost", 4223)

export disconnect
"""
    $(SIGNATURES)

Disconnects the TCP/IP connection from the Brick Daemon or the
WIFI/Ethernet Extension.
"""
disconnect(ipcon::IPConnection) = ipcon.ipconInternal.disconnect()

export authenticate
"""
    $(SIGNATURES)

Performs an authentication handshake with the connected Brick Daemon or
WIFI/Ethernet Extension. If the handshake succeeds the connection switches
from non-authenticated to authenticated state and communication can
continue as normal. If the handshake fails then the connection gets closed.
Authentication can fail if the wrong secret was used or if authentication
is not enabled at all on the Brick Daemon or the WIFI/Ethernet Extension.
For more information about authentication see
https://www.tinkerforge.com/en/doc/Tutorials/Tutorial_Authentication/Tutorial.html
"""
authenticate(ipcon::IPConnection, secret::String) = ipcon.ipconInternal.authenticate(secret)

export get_connection_state
"""
    $(SIGNATURES)

Can return the following states:
- CONNECTION_STATE_DISCONNECTED: No connection is established.
- CONNECTION_STATE_CONNECTED: A connection to the Brick Daemon or
  the WIFI/Ethernet Extension is established.
- CONNECTION_STATE_PENDING: IP Connection is currently trying to
  connect.
"""
get_connection_state(ipcon::IPConnection) = TinkerforgeIPConConnectionState(ipcon.ipconInternal.get_connection_state())

export set_auto_reconnect
"""
    $(SIGNATURES)

Enables or disables auto-reconnect. If auto-reconnect is enabled,
the IP Connection will try to reconnect to the previously given
host and port, if the connection is lost.
Default value is *True*.
"""
set_auto_reconnect(ipcon::IPConnection, auto_reconnect::Bool) = ipcon.ipconInternal.set_auto_reconnect(auto_reconnect)

export get_auto_reconnect
"""
    $(SIGNATURES)

Returns *true* if auto-reconnect is enabled, *false* otherwise.
"""
get_auto_reconnect(ipcon::IPConnection) = ipcon.ipconInternal.get_auto_reconnect()

export set_timeout
"""
    $(SIGNATURES)

Sets the timeout in seconds for getters and for setters for which the
response expected flag is activated.
Default timeout is 2.5.
"""
set_timeout(ipcon::IPConnection, timeout::Real) = ipcon.ipconInternal.set_timeout(timeout)

export get_timeout
"""
    $(SIGNATURES)

Returns the timeout as set by set_timeout.
"""
get_timeout(ipcon::IPConnection) = ipcon.ipconInternal.get_timeout()

export enumerate
"""
    $(SIGNATURES)

Broadcasts an enumerate request. All devices will respond with an
enumerate callback.
"""
enumerate(ipcon::IPConnection) = ipcon.ipconInternal.enumerate()

export wait
"""
    $(SIGNATURES)

Stops the current thread until unwait is called.
This is useful if you rely solely on callbacks for events, if you want
to wait for a specific callback or if the IP Connection was created in
a thread.
Wait and unwait act in the same way as "acquire" and "release" of a
semaphore.
"""
wait(ipcon::IPConnection) =  ipcon.ipconInternal.wait()

export unwait
"""
    $(SIGNATURES)

Unwaits the thread previously stopped by wait.
Wait and unwait act in the same way as "acquire" and "release" of
a semaphore.
"""
unwait(ipcon::IPConnection) = ipcon.ipconInternal.unwait()

export register_callback
"""
    $(SIGNATURES)

Registers the given *function* with the given *callback_id*.
"""
register_callback(ipcon::IPConnection, callback_id::Integer, func::Function) = ipcon.ipconInternal.register_callback(callback_id, func)


text:
module BatchedRoutines

function reshaped_batch_size(n::Int, A::AbstractArray)
    return prod(k->size(A, k+n), ndims(A) - n)
end

function reshaped_batch_size(n::Int, A::AbstractArray{T, 3}) where T
    n == 2 && return size(A, 3)
    return prod(k->size(A, k+n), ndims(A) - n)
end

macro iterate_batch(elty, Xs, dims, f)
    ex = Expr(:block)
    for X in Xs.args
        push!(ex.args, :(BLAS.chkstride1($(esc(X)))))
    end
    # check batch size
    batch_check = Expr(:comparison)
    for k in 1:length(Xs.args)-1
        X = Xs.args[k]; n = dims.args[k]
        push!(batch_check.args, :(reshaped_batch_size($n, $(esc(X)))))
        push!(batch_check.args, :(==))
    end
    push!(ex.args, :($(esc(:batch_n)) = reshaped_batch_size($(dims.args[end]), $(esc(Xs.args[end])))))
    push!(batch_check.args, esc(:batch_n))
    push!(ex.args, Expr(:(||), batch_check, :(error("batch size mismatch"))))

    # assign batch size

    # generated ptrs
    for_block = Expr(:block, esc(f))

    for (X, n) in zip(Xs.args, dims.args)
        ptrX = Symbol(:ptr, X)
        push!(ex.args, :($(esc(ptrX)) = Base.unsafe_convert(Ptr{$(esc(elty))}, $(esc(X)))))
        push!(for_block.args, :($(esc(ptrX)) += stride($(esc(X)), $(n + 1)) * sizeof($(esc(elty)))))
    end

    forloop = Expr(:for, :($(esc(:batch_k)) = 1:$(esc(:batch_n))), for_block)
    push!(ex.args, forloop)
    return ex
end

include("blas.jl")
include("lapack.jl")
include("linalg.jl")

end # module


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct NonMaximumSuppression_PickTop <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function NonMaximumSuppression_PickTop(; kwargs...)
        obj = new(meta(NonMaximumSuppression_PickTop), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct NonMaximumSuppression_PickTop
const __meta_NonMaximumSuppression_PickTop = Ref{ProtoMeta}()
function meta(::Type{NonMaximumSuppression_PickTop})
    ProtoBuf.metalock() do
        if !isassigned(__meta_NonMaximumSuppression_PickTop)
            __meta_NonMaximumSuppression_PickTop[] = target = ProtoMeta(NonMaximumSuppression_PickTop)
            allflds = Pair{Symbol,Union{Type,String}}[:perClass => Bool]
            meta(target, NonMaximumSuppression_PickTop, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_NonMaximumSuppression_PickTop[]
    end
end
function Base.getproperty(obj::NonMaximumSuppression_PickTop, name::Symbol)
    if name === :perClass
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

mutable struct NonMaximumSuppression <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function NonMaximumSuppression(; kwargs...)
        obj = new(meta(NonMaximumSuppression), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct NonMaximumSuppression
const __meta_NonMaximumSuppression = Ref{ProtoMeta}()
function meta(::Type{NonMaximumSuppression})
    ProtoBuf.metalock() do
        if !isassigned(__meta_NonMaximumSuppression)
            __meta_NonMaximumSuppression[] = target = ProtoMeta(NonMaximumSuppression)
            fnum = Int[1,100,101,110,111,200,201,202,203,210,211]
            allflds = Pair{Symbol,Union{Type,String}}[:pickTop => NonMaximumSuppression_PickTop, :stringClassLabels => StringVector, :int64ClassLabels => Int64Vector, :iouThreshold => Float64, :confidenceThreshold => Float64, :confidenceInputFeatureName => AbstractString, :coordinatesInputFeatureName => AbstractString, :iouThresholdInputFeatureName => AbstractString, :confidenceThresholdInputFeatureName => AbstractString, :confidenceOutputFeatureName => AbstractString, :coordinatesOutputFeatureName => AbstractString]
            oneofs = Int[1,2,2,0,0,0,0,0,0,0,0]
            oneof_names = Symbol[Symbol("SuppressionMethod"),Symbol("ClassLabels")]
            meta(target, NonMaximumSuppression, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_NonMaximumSuppression[]
    end
end
function Base.getproperty(obj::NonMaximumSuppression, name::Symbol)
    if name === :pickTop
        return (obj.__protobuf_jl_internal_values[name])::NonMaximumSuppression_PickTop
    elseif name === :stringClassLabels
        return (obj.__protobuf_jl_internal_values[name])::StringVector
    elseif name === :int64ClassLabels
        return (obj.__protobuf_jl_internal_values[name])::Int64Vector
    elseif name === :iouThreshold
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :confidenceThreshold
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :confidenceInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :coordinatesInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :iouThresholdInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :confidenceThresholdInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :confidenceOutputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :coordinatesOutputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

export NonMaximumSuppression_PickTop, NonMaximumSuppression


text:
__precompile__(true)

module QPnorm

using LinearAlgebra
using Polynomials
using GeneralQP
using TRS
using DataFrames
using Printf
using SparseArrays
using CSV

include("algebra.jl")
include("qp_norm_boundary.jl")
include("printing.jl")
include("qp_norm.jl")
include("find_feasible_point.jl")
export solve, find_feasible_point

end


text:
export adaboost

import ..AdaBoostModel

using mlpack._Internal.io

import mlpack_jll
const adaboostLibrary = mlpack_jll.libmlpack_julia_adaboost

# Call the C binding of the mlpack adaboost binding.
function adaboost_mlpackMain()
  success = ccall((:adaboost, adaboostLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module adaboost_internal
  import ..adaboostLibrary

import ...AdaBoostModel

# Get the value of a model pointer parameter of type AdaBoostModel.
function IOGetParamAdaBoostModel(paramName::String)::AdaBoostModel
  AdaBoostModel(ccall((:IO_GetParamAdaBoostModelPtr, adaboostLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type AdaBoostModel.
function IOSetParamAdaBoostModel(paramName::String, model::AdaBoostModel)
  ccall((:IO_SetParamAdaBoostModelPtr, adaboostLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeAdaBoostModel(stream::IO, model::AdaBoostModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeAdaBoostModelPtr, adaboostLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeAdaBoostModel(stream::IO)::AdaBoostModel
  buffer = read(stream)
  AdaBoostModel(ccall((:DeserializeAdaBoostModelPtr, adaboostLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    adaboost(; [input_model, iterations, labels, test, tolerance, training, verbose, weak_learner])

This program implements the AdaBoost (or Adaptive Boosting) algorithm. The
variant of AdaBoost implemented here is AdaBoost.MH. It uses a weak learner,
either decision stumps or perceptrons, and over many iterations, creates a
strong learner that is a weighted ensemble of weak learners. It runs these
iterations until a tolerance value is crossed for change in the value of the
weighted training error.

For more information about the algorithm, see the paper "Improved Boosting
Algorithms Using Confidence-Rated Predictions", by R.E. Schapire and Y. Singer.

This program allows training of an AdaBoost model, and then application of that
model to a test dataset.  To train a model, a dataset must be passed with the
`training` option.  Labels can be given with the `labels` option; if no labels
are specified, the labels will be assumed to be the last column of the input
dataset.  Alternately, an AdaBoost model may be loaded with the `input_model`
option.

Once a model is trained or loaded, it may be used to provide class predictions
for a given test dataset.  A test dataset may be specified with the `test`
parameter.  The predicted classes for each point in the test dataset are output
to the `predictions` output parameter.  The AdaBoost model itself is output to
the `output_model` output parameter.

Note: the following parameter is deprecated and will be removed in mlpack 4.0.0:
`output`.
Use `predictions` instead of `output`.

For example, to run AdaBoost on an input dataset `data` with labels `labels`and
perceptrons as the weak learner type, storing the trained model in `model`, one
could use the following command:

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> labels = CSV.read("labels.csv"; type=Int)
julia> _, model, _, _ = adaboost(labels=labels, training=data,
            weak_learner="perceptron")
```

Similarly, an already-trained model in `model` can be used to provide class
predictions from test data `test_data` and store the output in `predictions`
with the following command:

```julia
julia> using CSV
julia> test_data = CSV.read("test_data.csv")
julia> _, _, predictions, _ = adaboost(input_model=model,
            test=test_data)
```

# Arguments

 - `input_model::AdaBoostModel`: Input AdaBoost model.
 - `iterations::Int`: The maximum number of boosting iterations to be run
      (0 will run until convergence.)  Default value `1000`.

 - `labels::Array{Int, 1}`: Labels for the training set.
 - `test::Array{Float64, 2}`: Test dataset.
 - `tolerance::Float64`: The tolerance for change in values of the
      weighted error during training.  Default value `1e-10`.

 - `training::Array{Float64, 2}`: Dataset for training AdaBoost.
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

 - `weak_learner::String`: The type of weak learner to use:
      'decision_stump', or 'perceptron'.  Default value `decision_stump`.

# Return values

 - `output::Array{Int, 1}`: Predicted labels for the test set.
 - `output_model::AdaBoostModel`: Output trained AdaBoost model.
 - `predictions::Array{Int, 1}`: Predicted labels for the test set.
 - `probabilities::Array{Float64, 2}`: Predicted class probabilities for
      each point in the test set.

"""
function adaboost(;
                  input_model::Union{AdaBoostModel, Missing} = missing,
                  iterations::Union{Int, Missing} = missing,
                  labels = missing,
                  test = missing,
                  tolerance::Union{Float64, Missing} = missing,
                  training = missing,
                  verbose::Union{Bool, Missing} = missing,
                  weak_learner::Union{String, Missing} = missing,
                  points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, adaboostLibrary), Nothing, ());

  IORestoreSettings("AdaBoost")

  # Process each input argument before calling mlpackMain().
  if !ismissing(input_model)
    adaboost_internal.IOSetParamAdaBoostModel("input_model", convert(AdaBoostModel, input_model))
  end
  if !ismissing(iterations)
    IOSetParam("iterations", convert(Int, iterations))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(tolerance)
    IOSetParam("tolerance", convert(Float64, tolerance))
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if !ismissing(weak_learner)
    IOSetParam("weak_learner", convert(String, weak_learner))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  IOSetPassed("output_model")
  IOSetPassed("predictions")
  IOSetPassed("probabilities")
  # Call the program.
  adaboost_mlpackMain()

  return IOGetParamURow("output"),
         adaboost_internal.IOGetParamAdaBoostModel("output_model"),
         IOGetParamURow("predictions"),
         IOGetParamMat("probabilities", points_are_rows)
end


text:
module FusionDirect

using OpenGene
using OpenGene.Algorithm
using OpenGene.Reference

# package code goes here

# make it compatible for different version of Julia
include("compat.jl")

export detect

import Base: -,
    abs

include("index/index.jl")
include("detect/detect.jl")

end # module


text:
function request_handler(req::HTTP.Request)
    target = String(HTTP.Messages.getfield(req, :target))
    uri = HTTP.URI(target)
    query = uri.query
    urisstring = unescapeuri(query)
    params = HTTP.queryparams(urisstring) # create URI out of target String
    wewant = ["access_token", "token_type", "state", "expires_in"]
    ipa = filter(pa-> in(pa[1], wewant), params)
    if length(ipa) > 3
        if get(ipa, "state", "") == "987"
            if get(ipa, "token_type", "") == "Bearer"
                c = spotcred()::SpotifyCredentials
                c.ig_access_token = get(ipa, "access_token", "")
                c.token_type = get(ipa, "token_type", "")
                expi_in = parse(Int, get(ipa, "expires_in", 0))
                c.expires_at = string(Dates.now() + Dates.Second(expi_in))
                SPOTCRED[] = c
                # Feedback in 'wait_for_ig_access'
            else
                @error "Unexpected token_type"
            end
        else
            @error "Unexpected state"
        end
        HTTP.Response(202, repr(ipa))
    elseif length(ipa) > 0
        @error "Unexpected query string: $query"
        HTTP.Response(202, repr(ipa))
    else
        if target == "/favicon.ico"
            HTTP.Response(404)
        else
            fpath = joinpath(@__DIR__, "standardize_uri.html")
            HTTP.Response( 200, body = read(fpath ) )
        end
    end
end

function has_ig_access_token()
    c = spotcred()::SpotifyCredentials
    t = c.ig_access_token
    t != ""
end
function wait_for_ig_access(; stopwait::DateTime = now() + Dates.Second(20))
    while !has_ig_access_token() && now() < stopwait
        sleep(0.1)
    end
    now() >= stopwait  ? println("\n\tTimeout server") : println("\n\tImplicit grant token expires at $(SPOTCRED[].expires_at) - closing server")
end

function start_a_listening_server()
    stored = spotcred().redirect
    stripped = strip(stored, '/')
    u = parse(URI, stripped)
    host() = parse(IPAddr, u.host)
    port() = parse(Int, u.port)
    path() = u.path
    server = Sockets.listen(host(), port())
    close_server_when_ready_task = @async begin
        wait_for_ig_access()
        close(server)
    end
    println("\tListening for authorization on $(host()):$(port()) and path $(path())")
    listen_task = @async HTTP.serve(request_handler, host(), port(); server, stream =false)
    listen_task, close_server_when_ready_task
end

function launch_a_browser_that_asks_for_implicit_grant()
    scopes = join(DEFAULT_IMPLICIT_GRANT, "%20")
    uri = OAUTH_AUTHORIZE_URL * "?" *
            "client_id=" * spotcred().client_id *
            "&redirect_uri=" * replace(spotcred().redirect, "/" => "%2F") *
            "&scope=" * scopes *
            "&show_dialog=true" *
            "&response_type=token" *
            "&state=987"
    println("\tLaunching a browser at: $uri")

    # Try opening a browser window based on OS type
    if Sys.isapple()
        run(`open $uri`)
    end
    if Sys.iswindows()
        run(`explorer $uri`)
    end
    if Sys.islinux()
        run(`xdg-open $uri`)
    end

    # Added delay so that user can authorize via browser window while running unit tests
    @info "Waiting for 10 seconds"
    sleep(10)
    @info "Try running that command again, you should be authorized now"
end

"""
If the user gives us access, that will be in the future.
Is it worth waiting for? If so, call `wait_for_ig_access()`
"""
function get_implicit_grant()
    listen_task, close_server_when_ready_task = start_a_listening_server()
    launch_a_browser_that_asks_for_implicit_grant()
    listen_task, close_server_when_ready_task
end


text:
using Lazy
import StatsBase

abstract type AbstractParameterFunction end

function Base.get(apf::AbstractParameterFunction, state_t, action_t, state_tp1, action_tp1, preds_tilde) end

function call(apf::AbstractParameterFunction, state_t, action_t, state_tp1, action_tp1, preds_tilde)
    Base.get(apf::AbstractParameterFunction, state_t, action_t, state_tp1, action_tp1, preds_tilde)
end

"""
Cumulants
"""

abstract type AbstractCumulant <: AbstractParameterFunction end

function Base.get(cumulant::AbstractCumulant, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    throw(DomainError("Base.get(CumulantType, args...) not defined!"))
end

"""
    FeatureCumulant
    - Basic Cumulant which takes the value c_t = s_tp1[idx] for 1<=idx<=length(s_tp1)
"""
struct FeatureCumulant <: AbstractCumulant
    idx::Int
end

Base.get(cumulant::FeatureCumulant, state_t, action_t, state_tp1, action_tp1, preds_tp1) = state_tp1[cumulant.idx]

struct PredictionCumulant <: AbstractCumulant
    idx::Int
end

Base.get(cumulant::PredictionCumulant, state_t, action_t, state_tp1, action_tp1, preds_tp1) = preds_tp1[cumulant.idx]

struct FunctionalCumulant{F} <: AbstractCumulant
    func::F
end

Base.get(cumulant::FunctionalCumulant, state_t, action_t, state_tp1, action_tp1, preds_tp1) =
    cumulant.func(state_t, action_t, state_tp1, action_tp1, preds_tp1)

"""
Discounting
"""
abstract type AbstractDiscount <: AbstractParameterFunction end

function Base.get(γ::AbstractDiscount, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    throw(DomainError("Base.get(DiscountType, args...) not defined!"))
end

struct ConstantDiscount{T} <: AbstractDiscount
    γ::T
end

Base.get(cd::ConstantDiscount, state_t, action_t, state_tp1, action_tp1, preds_tp1) = cd.γ

struct StateTerminationDiscount{T<:Number, F} <: AbstractDiscount
    γ::T
    condition::F
    terminal::T
    StateTerminationDiscount(γ, condition) = new{typeof(γ), typeof(condition)}(γ, condition, convert(typeof(γ), 0.0))
end

Base.get(td::StateTerminationDiscount, state_t, action_t, state_tp1, action_tp1, preds_tp1) =
    td.condition(state_t, action_t, state_tp1) ? td.terminal : td.γ

"""
Policies
"""
abstract type AbstractPolicy <: AbstractParameterFunction end

function Base.get(π::AbstractPolicy, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    throw(DomainError("Base.get(PolicyType, args...) not defined!"))
end

StatsBase.sample(rng::Random.AbstractRNG, π::AbstractPolicy, state) = nothing
StatsBase.sample(π::AbstractPolicy, state) = StatsBase.sample(Random.GLOBAL_RNG, π, state)

struct NullPolicy <: AbstractPolicy
end
Base.get(π::NullPolicy, state_t, action_t, state_tp1, action_tp1, preds_tp1) = 1.0

struct PersistentPolicy <: AbstractPolicy
    action::Int64
end

Base.get(π::PersistentPolicy, state_t, action_t, state_tp1, action_tp1, preds_tp1) = π.action == action_t ? 1.0 : 0.0
StatsBase.sample(rng::Random.AbstractRNG, π::PersistentPolicy, state=nothing) = π.action

struct RandomPolicy{T<:AbstractFloat} <: AbstractPolicy
    probabilities::Array{T,1}
    weight_vec::StatsBase.Weights{T, T, Array{T, 1}}
    RandomPolicy(probabilities::Array{T,1}) where {T<:AbstractFloat} = new{T}(probabilities, StatsBase.Weights(probabilities))
end

UniformRandomPolicy(num_actions::Integer) = RandomPolicy(fill(1.0/num_actions, num_actions))
FavoredRandomPolicy(num_actions::Integer, favored_action::Integer, favored_prob::T) where {T<:AbstractFloat} =
    RandomPolicy([ act == favored_action ? favored_prob : (T(1.0) - favored_prob)/(num_actions-1) for act in 1:num_actions])

Base.get(π::RandomPolicy, state_t, action_t, state_tp1, action_tp1, preds_tp1) = π.probabilities[action_t]

StatsBase.sample(rng::Random.AbstractRNG, π::RandomPolicy, state=nothing) = StatsBase.sample(rng, 1:length(π.weight_vec), π.weight_vec)

struct InterpolationPolicy{P<:AbstractPolicy, M<:AbstractPolicy} <: AbstractPolicy
    beta::Float64
    policy_1::P
    policy_2::M
end

Base.get(π::InterpolationPolicy, state_t, action_t, state_tp1, action_tp1, preds_tp1) =
    π.beta*Base.get(π.policy_1, state_t, action_t, state_tp1, action_tp1, preds_tp1) +
    (1.0-π.beta)*Base.get(π.policy_2, state_t, action_t, state_tp1, action_tp1, preds_tp1)

struct FunctionalPolicy{F} <: AbstractPolicy
    func::F
end

Base.get(π::FunctionalPolicy, state_t, action_t, state_tp1, action_tp1, preds_tp1) =
    π.func(state_t, action_t, state_tp1, action_tp1, preds_tp1)

abstract type AbstractGVF end

function Base.get(gvf::AbstractGVF, state_t, action_t, state_tp1, action_tp1, preds_tp1) end

Base.get(gvf::AbstractGVF, state_t, action_t, state_tp1, preds_tp1) =
    Base.get(gvf::AbstractGVF, state_t, action_t, state_tp1, nothing, preds_tp1)

Base.get(gvf::AbstractGVF, state_t, action_t, state_tp1) =
    Base.get(gvf::AbstractGVF, state_t, action_t, state_tp1, nothing, nothing)

function cumulant(gvf::AbstractGVF) end
function discount(gvf::AbstractGVF) end
function policy(gvf::AbstractGVF) end

struct GVF{C<:AbstractCumulant, D<:AbstractDiscount, P<:AbstractPolicy} <: AbstractGVF
    cumulant::C
    discount::D
    policy::P
end

cumulant(gvf::GVF) = gvf.cumulant
discount(gvf::GVF) = gvf.discount
policy(gvf::GVF) = gvf.policy

function Base.get(gvf::GVF, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    c = Base.get(gvf.cumulant, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    γ = Base.get(gvf.discount, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    π_prob = Base.get(gvf.policy, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    return c, γ, π_prob
end

abstract type AbstractHorde end

struct Horde{T<:AbstractGVF} <: AbstractHorde
    gvfs::Vector{T}
end

function Base.get(gvfh::Horde, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    C = map(gvf -> Base.get(cumulant(gvf), state_t, action_t, state_tp1, action_tp1, preds_tp1), gvfh.gvfs)
    Γ = map(gvf -> Base.get(discount(gvf), state_t, action_t, state_tp1, action_tp1, preds_tp1), gvfh.gvfs)
    Π_probs = map(gvf -> Base.get(policy(gvf), state_t, action_t, state_tp1, action_tp1, preds_tp1), gvfh.gvfs)
    return C, Γ, Π_probs
end

function Base.get!(C::Array{T, 1}, Γ::Array{T, 1}, Π_probs::Array{T, 1}, gvfh::Horde, state_t, action_t, state_tp1, action_tp1, preds_tp1) where {T<:AbstractFloat}
    C .= map(gvf -> Base.get(cumulant(gvf), state_t, action_t, state_tp1, action_tp1, preds_tp1), gvfh.gvfs)
    Γ .= map(gvf -> Base.get(discount(gvf), state_t, action_t, state_tp1, action_tp1, preds_tp1), gvfh.gvfs)
    Π_probs .= map(gvf -> Base.get(policy(gvf), state_t, action_t, state_tp1, action_tp1, preds_tp1), gvfh.gvfs)
    return C, Γ, Π_probs
end

Base.get(gvfh::Horde, state_tp1, preds_tp1) = Base.get(gvfh::Horde, nothing, nothing, state_tp1, nothing, preds_tp1)
Base.get(gvfh::Horde, action_t, state_tp1, preds_tp1) = Base.get(gvfh::Horde, nothing, action_t, state_tp1, nothing, preds_tp1)
Base.get(gvfh::Horde, state_t, action_t, state_tp1, preds_tp1) = Base.get(gvfh::Horde, state_t, action_t, state_tp1, nothing, preds_tp1)

Base.get!(C, Γ, Π_probs,gvfh::Horde, action_t, state_tp1, preds_tp1) = Base.get!(C, Γ, Π_probs, gvfh::Horde, nothing, action_t, state_tp1, nothing, preds_tp1)

@forward Horde.gvfs Base.length


text:
@testset "PSSE Parsing " begin
    files = readdir(PSSE_RAW_DIR)
    if length(files) == 0
        error("No test files in the folder")
    end

    for f in files[1:1]
        @info "Parsing $f ..."
        pm_data = PowerSystems.PowerModelsData(joinpath(PSSE_RAW_DIR, f))
        @info "Successfully parsed $f to PowerModelsData"
        sys = System(pm_data)
        for g in get_components(Generator, sys)
            @test haskey(get_ext(g), "z_source")
        end
        @info "Successfully parsed $f to System struct"
    end

    # Test bad input
    pm_data = PowerSystems.PowerModelsData(joinpath(PSSE_RAW_DIR, files[1]))
    pm_data.data["bus"] = Dict{String, Any}()
    @test_throws PowerSystems.DataFormatError System(pm_data)
end

@testset "PSSE PowerModel Parsing " begin
    PM_PSSE_PATH = joinpath(DATA_DIR, "pm_data", "pti")
    files = readdir(PM_PSSE_PATH)
    if length(files) == 0
        error("No test files in the folder")
    end

    for f in files[1:1]
        @info "Parsing $f ..."
        sys = System(joinpath(PM_PSSE_PATH, f))
        for g in get_components(Generator, sys)
            @test haskey(get_ext(g), "z_source")
        end
        @info "Successfully parsed $f to System struct"
    end
end


text:
module TestBootstrap
    using StreamStats
    using Distributions
    using Base.Test

    # CI's for mean of uniform draws
    for n in rand(1:10_000, 10)
        xs = rand(n)
        inner_stat = StreamStats.Mean()
        stat1 = StreamStats.BernoulliBootstrap(inner_stat)
        stat2 = StreamStats.PoissonBootstrap(inner_stat)
        for x in xs
            update!(stat1, x)
            update!(stat2, x)
        end

        # test using quantiles to generate bootstrap CIs
        lower1, upper1 = StreamStats.ci(stat1, 0.05)
        lower2, upper2 = StreamStats.ci(stat2, 0.05)
        m, sem = mean(xs), std(xs) / sqrt(n)
        @test m - 3 * sem <= lower1 <= m - 1 * sem
        @test m + 1 * sem <= upper1 <= m + 3 * sem
        @test m - 3 * sem <= lower2 <= m - 1 * sem
        @test m + 1 * sem <= upper2 <= m + 3 * sem

        # test normal approximation of the bootstrap CIs
        lower1, upper1 = StreamStats.ci(stat1, 0.05, :normal)
        lower2, upper2 = StreamStats.ci(stat2, 0.05, :normal)
        m, sem = mean(xs), std(xs) / sqrt(n)
        @test m - 3 * sem <= lower1 <= m - 1 * sem
        @test m + 1 * sem <= upper1 <= m + 3 * sem
        @test m - 3 * sem <= lower2 <= m - 1 * sem
        @test m + 1 * sem <= upper2 <= m + 3 * sem
    end

    # CI's for variance of uniform draws
    for n in rand(1:10_000, 10)
        xs = rand(n)
        inner_stat = StreamStats.Variance()
        stat1 = StreamStats.BernoulliBootstrap(inner_stat)
        stat2 = StreamStats.PoissonBootstrap(inner_stat)
        for x in xs
            update!(stat1, x)
            update!(stat2, x)
        end
        lower1, upper1 = StreamStats.ci(stat1, 0.05, :quantile)
        lower2, upper2 = StreamStats.ci(stat2, 0.05, :quantile)
        @test lower1 <= var(xs) <= upper1
        @test lower2 <= var(xs) <= upper2

        n_lower1, n_upper1 = StreamStats.ci(stat1, 0.05, :normal)
        n_lower2, n_upper2 = StreamStats.ci(stat2 , 0.05, :normal)
        @test n_lower1 <= var(xs) <= n_upper1
        @test n_lower2 <= var(xs) <= n_upper2
    end

    # test rand() method
    n = 1_000
    b1 = StreamStats.BernoulliBootstrap(StreamStats.Mean(), 10_000)
    b2 = StreamStats.BernoulliBootstrap(StreamStats.Mean(), 10_000)
    xs = rand(n)
    for i in xs
        update!(b1, i)
        update!(b2, i / 2.0)
    end
    xb = [rand(b1) for i in 1:1_000]

    # test -() method
    # TODO: rename Boostrap.cached_state() to replicates(), so that
    # below can simply be mean(replicates(b2-b1))
    # @test mean(b2 - b1) < 0

end


text:
struct LazyArrayStyle{N} <: AbstractArrayStyle{N} end
LazyArrayStyle(::Val{N}) where N = LazyArrayStyle{N}()
LazyArrayStyle{M}(::Val{N}) where {N,M} = LazyArrayStyle{N}()

struct BroadcastArray{T, N, F, Args} <: LazyArray{T, N}
    f::F
    args::Args
end

const BroadcastVector{T,F,Args} = BroadcastArray{T,1,F,Args}
const BroadcastMatrix{T,F,Args} = BroadcastArray{T,2,F,Args}

LazyArray(bc::Broadcasted) = BroadcastArray(bc)

BroadcastArray{T,N,F,Args}(bc::Broadcasted) where {T,N,F,Args} = BroadcastArray{T,N,F,Args}(bc.f,bc.args)
BroadcastArray{T,N}(bc::Broadcasted{Style,Axes,F,Args}) where {T,N,Style,Axes,F,Args} = BroadcastArray{T,N,F,Args}(bc.f,bc.args)
BroadcastArray{T}(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Tuple{Vararg{Any,N}},<:Any,<:Tuple}) where {T,N} =
    BroadcastArray{T,N}(bc)

BroadcastVector(bc::Broadcasted) = BroadcastVector{combine_eltypes(bc.f, bc.args)}(bc)
BroadcastMatrix(bc::Broadcasted) = BroadcastMatrix{combine_eltypes(bc.f, bc.args)}(bc)

_broadcast2broadcastarray() = ()
_broadcast2broadcastarray(a, b...) = tuple(a, _broadcast2broadcastarray(b...)...)
_broadcast2broadcastarray(a::Broadcasted{DefaultArrayStyle{0}}, b...) = tuple(materialize(a), _broadcast2broadcastarray(b...)...)
_broadcast2broadcastarray(a::Broadcasted, b...) = tuple(BroadcastArray(a), _broadcast2broadcastarray(b...)...)

_BroadcastArray(bc::Broadcasted) = BroadcastArray{combine_eltypes(bc.f, bc.args)}(bc)
BroadcastArray(bc::Broadcasted{S}) where S =
    _BroadcastArray(instantiate(Broadcasted{S}(bc.f, _broadcast2broadcastarray(bc.args...))))

BroadcastArray(f, A, As...) = BroadcastArray(broadcasted(f, A, As...))
BroadcastMatrix(f, A...) = BroadcastMatrix(broadcasted(f, A...))
BroadcastVector(f, A...) = BroadcastVector(broadcasted(f, A...))

BroadcastArray{T,N}(f, A...) where {T,N} = BroadcastArray{T,N,typeof(f),typeof(A)}(f, A)

BroadcastArray(b::BroadcastArray) = b
BroadcastVector(A::BroadcastVector) = A
BroadcastMatrix(A::BroadcastMatrix) = A

Broadcasted(A::BroadcastArray) = instantiate(broadcasted(A.f, A.args...))

@inline BroadcastArray(A::AbstractArray) = BroadcastArray(call(A), arguments(A)...)

axes(A::BroadcastArray) = axes(Broadcasted(A))
size(A::BroadcastArray) = map(length, axes(A))

@propagate_inbounds getindex(A::BroadcastArray{<:Any,N}, kj::Vararg{Int,N}) where N = Broadcasted(A)[kj...]

@propagate_inbounds _broadcast_getindex_range(A::Union{Ref,AbstractArray{<:Any,0},Number}, I) = A[] # Scalar-likes can just ignore all indices
# Everything else falls back to dynamically dropping broadcasted indices based upon its axes
@propagate_inbounds _broadcast_getindex_range(A, I) = A[I]

getindex(B::BroadcastArray{<:Any,1}, kr::AbstractVector{<:Integer}) =
    BroadcastArray(Broadcasted(B).f, map(a -> _broadcast_getindex_range(a,kr), Broadcasted(B).args)...)
getindex(B::BroadcastArray{<:Any,1}, kr::AbstractUnitRange{<:Integer}) =
    BroadcastArray(Broadcasted(B).f, map(a -> _broadcast_getindex_range(a,kr), Broadcasted(B).args)...)

copy(bc::Broadcasted{<:LazyArrayStyle}) = BroadcastArray(bc)

copyto!(dest::AbstractArray{<:Any,N}, bc::BroadcastArray{<:Any,N}) where N =
    copyto!(dest, Broadcasted(bc))

# Replacement for #18.
# Could extend this to other similar reductions in Base... or apply at lower level?
# for (fname, op) in [(:sum, :add_sum), (:prod, :mul_prod),
#                     (:maximum, :max), (:minimum, :min),
#                     (:all, :&),       (:any, :|)]
function Base._sum(f, A::BroadcastArray, ::Colon)
    bc = Broadcasted(A)
    T = Broadcast.combine_eltypes(f ∘ bc.f, bc.args)
    out = zero(T)
    @simd for I in eachindex(bc)
        @inbounds out += f(bc[I])
    end
    out
end
function Base._prod(f, A::BroadcastArray, ::Colon)
    bc = Broadcasted(A)
    T = Broadcast.combine_eltypes(f ∘ bc.f, bc.args)
    out = one(T)
    @simd for I in eachindex(bc)
        @inbounds out *= f(bc[I])
    end
    out
end

BroadcastStyle(::Type{<:LazyArray{<:Any,N}}) where N = LazyArrayStyle{N}()
BroadcastStyle(::Type{<:Adjoint{<:Any,<:LazyVector{<:Any}}}) where N = LazyArrayStyle{2}()
BroadcastStyle(::Type{<:Transpose{<:Any,<:LazyVector{<:Any}}}) where N = LazyArrayStyle{2}()
BroadcastStyle(::Type{<:Adjoint{<:Any,<:LazyMatrix{<:Any}}}) where N = LazyArrayStyle{2}()
BroadcastStyle(::Type{<:Transpose{<:Any,<:LazyMatrix{<:Any}}}) where N = LazyArrayStyle{2}()
BroadcastStyle(L::LazyArrayStyle{N}, ::StaticArrayStyle{N}) where N = L
BroadcastStyle(::StaticArrayStyle{N}, L::LazyArrayStyle{N})  where N = L

"""
    BroadcastLayout{F}()

is returned by `MemoryLayout(A)` if a matrix `A` is a `BroadcastArray`.
`F` is the typeof function that broadcast operation is applied.
"""
struct BroadcastLayout{F} <: MemoryLayout end

tuple_type_memorylayouts(::Type{I}) where I<:Tuple = MemoryLayout.(I.parameters)
tuple_type_memorylayouts(::Type{Tuple{A}}) where {A} = (MemoryLayout(A),)
tuple_type_memorylayouts(::Type{Tuple{A,B}}) where {A,B} = (MemoryLayout(A),MemoryLayout(B))
tuple_type_memorylayouts(::Type{Tuple{A,B,C}}) where {A,B,C} = (MemoryLayout(A),MemoryLayout(B),MemoryLayout(C))

broadcastlayout(::Type{F}, _...) where F = BroadcastLayout{F}()
broadcastlayout(::Type, ::LazyLayout...) = LazyLayout()
broadcastlayout(::Type, _, ::LazyLayout) = LazyLayout()
broadcastlayout(::Type, _, _, ::LazyLayout) = LazyLayout()
broadcastlayout(::Type, _, _, _, ::LazyLayout) = LazyLayout()
MemoryLayout(::Type{BroadcastArray{T,N,F,Args}}) where {T,N,F,Args} =
    broadcastlayout(F, tuple_type_memorylayouts(Args)...)
## scalar-range broadcast operations ##
# Ranges already support smart broadcasting
for op in (+, -, big)
    @eval begin
        broadcasted(::LazyArrayStyle{1}, ::typeof($op), r::AbstractRange) =
            broadcast(DefaultArrayStyle{1}(), $op, r)
    end
end

for op in (-, +, *, /)
    @eval broadcasted(::LazyArrayStyle{1}, ::typeof($op), r::AbstractRange, x::Real) =
        broadcast(DefaultArrayStyle{1}(), $op, r, x)
end

for op in (-, +, *, \)
    @eval broadcasted(::LazyArrayStyle{1}, ::typeof($op), x::Real, r::AbstractRange) =
        broadcast(DefaultArrayStyle{1}(), $op, x, r)
end

broadcasted(::LazyArrayStyle{N}, op, r::AbstractFill{T,N}) where {T,N} =
    broadcast(DefaultArrayStyle{N}(), op, r)
broadcasted(::LazyArrayStyle{N}, op, r::AbstractFill{T,N}, x::Number) where {T,N} =
    broadcast(DefaultArrayStyle{N}(), op, r, x)
broadcasted(::LazyArrayStyle{N}, op, x::Number, r::AbstractFill{T,N}) where {T,N} =
    broadcast(DefaultArrayStyle{N}(), op, x, r)
broadcasted(::LazyArrayStyle{N}, op, r::AbstractFill{T,N}, x::Ref) where {T,N} =
    broadcast(DefaultArrayStyle{N}(), op, r, x)
broadcasted(::LazyArrayStyle{N}, op, x::Ref, r::AbstractFill{T,N}) where {T,N} =
    broadcast(DefaultArrayStyle{N}(), op, x, r)
broadcasted(::LazyArrayStyle{N}, op, r1::AbstractFill{T,N}, r2::AbstractFill{V,N}) where {T,V,N} =
    broadcast(DefaultArrayStyle{N}(), op, r1, r2)

broadcasted(::LazyArrayStyle{N}, ::typeof(*), a::Zeros{T,N}, b::Zeros{V,N}) where {T,V,N} =
    broadcast(DefaultArrayStyle{N}(), *, a, b)
broadcasted(::LazyArrayStyle{N}, ::typeof(*), a::AbstractArray{T,N}, b::Zeros{V,N}) where {T,V,N} =
    broadcast(DefaultArrayStyle{N}(), *, a, b)
broadcasted(::LazyArrayStyle{N}, ::typeof(*), a::Zeros{T,N}, b::AbstractArray{V,N}) where {T,V,N} =
    broadcast(DefaultArrayStyle{N}(), *, a, b)

###
# support
###

_broadcast_colsupport(sz, A::Number, j) = OneTo(sz[1])
_broadcast_colsupport(sz, A::AbstractVector, j) = colsupport(A,j)
_broadcast_colsupport(sz, A::AbstractMatrix, j) = size(A,1) == 1 ? OneTo(sz[1]) : colsupport(A,j)
_broadcast_rowsupport(sz, A::Number, j) = OneTo(sz[2])
_broadcast_rowsupport(sz, A::AbstractVector, j) = OneTo(sz[2])
_broadcast_rowsupport(sz, A::AbstractMatrix, j) = size(A,2) == 1 ? OneTo(sz[2]) : rowsupport(A,j)

colsupport(::BroadcastLayout{typeof(*)}, A, j) = intersect(_broadcast_colsupport.(Ref(size(A)), A.args, j)...)
rowsupport(::BroadcastLayout{typeof(*)}, A, j) = intersect(_broadcast_rowsupport.(Ref(size(A)), A.args, j)...)

for op in (:+, :-)
    @eval begin
        colsupport(::BroadcastLayout{typeof($op)}, A, j) = convexunion(_broadcast_colsupport.(Ref(size(A)), A.args, j)...)
        rowsupport(::BroadcastLayout{typeof($op)}, A, j) = convexunion(_broadcast_rowsupport.(Ref(size(A)), A.args, j)...)
    end
end

###
# SubArray
###

call(b::BroadcastLayout, a::SubArray) = call(b, parent(a))

sublayout(b::BroadcastLayout, _) = b

_broadcastviewinds(::Tuple{}, inds) = ()
_broadcastviewinds(sz, inds) =
    tuple(isone(sz[1]) ? OneTo(sz[1]) : inds[1], _broadcastviewinds(tail(sz), tail(inds))...)

_broadcastview(a, inds) = view(a, _broadcastviewinds(size(a), inds)...)

function arguments(b::BroadcastLayout, V::SubArray)
    args = arguments(parent(V))
    _broadcastview.(args, Ref(parentindices(V)))
end

###
# Transpose
###

call(b::BroadcastLayout, a::AdjOrTrans) = call(b, parent(a))

transposelayout(b::BroadcastLayout) = b
arguments(b::BroadcastLayout, A::Adjoint) = map(adjoint, arguments(b, parent(A)))
arguments(b::BroadcastLayout, A::Transpose) = map(transpose, arguments(b, parent(A)))


text:
"""
    EasyIO

user-friendly I/O methods. It covers:
1. timeTag(), generating time taggings
2. readcsv() & writecsv(), recoverage of two deprecated functions since Julia 0.6

"""
module EasyIO
   import LinearAlgebra: diagind, diag  # std lib, diagind() gets index range for an assigned offline parameter, used in diagwrite()
   import Dates  # std lib, for time taggings
   import DelimitedFiles # std lib, for delimiter-separated file I/O
   export timeTag, readcsv, writecsv, diagw, diagw!, expand, expand!

# -------------
"""
   timeTag()

generates a time tagging; the tag can be used to mark time or name files.
returns a String.
"""
function timeTag()
   local tagstr = replace(string(Dates.now()), ":" => "_" )
   # tagstr = replace( tagstr, "T" => "_" )
   tagstr = replace( tagstr, "." => "_" )
   return tagstr::String
end
# -------------
# Two IO functions from Julia 0.6
# NOTE: well ... the readcsv() & writecsv() before Julai 0.6 are SO convenient in some ways ...
#       I re-define the two functions here; (directly copied from base Julia 0.6)
readcsv(io; opts...)          = DelimitedFiles.readdlm(io, ','; opts...)
readcsv(io, T::Type; opts...) = DelimitedFiles.readdlm(io, ',', T; opts...)
"""
    writecsv(filename, A; opts)
Equivalent to [`writedlm`](@ref) with `delim` set to comma.
"""
writecsv(io, a; opts...) = DelimitedFiles.writedlm(io, a, ','; opts...)
# -------------
"""
    diagw( MAT::Matrix, VEC::Union{Vector,Tuple,NamedTuple,AbstractRange} ; offset::Int = 0 )

Writes a vector "VEC" into a matrix "MAT" diagonally, according to offset= (value starts from 0).
Returns a copy of MAT;

Depends:
1. LinearAlgebra.diagind()

Note:
1. diagind(MAT) (the index of positions to fill) can have different length from vector VEC
   1. if length(diagind(MAT)) < length(VEC): extra elements of VEC will be neglected
   2. if length(diagind(MAT)) > length(VEC): extra elements of MAT[diagind(MAT)] will be kept
"""
function diagw( MAT::Matrix, VEC::Union{Vector,Tuple,NamedTuple,AbstractRange} ; offset::Int = 0 )
   # declare
   local Res::Matrix = copy(MAT)
   local tmpIdxRange = diagind(Res, offset) # get indexes range
   # fill elements
   for s in 1:min( length(VEC), length(tmpIdxRange) ); Res[ tmpIdxRange[s] ] = VEC[s]; end
   return Res::Matrix
end # end function diagw
# ----------------
"""
    diagw!( MAT::Matrix, VEC::Union{Vector,Tuple,NamedTuple,AbstractRange} ; offset::Int = 0 )

Writes a vector "VEC" into a matrix "MAT" diagonally (in-place), according to offset= (value starts from 0).
Returns nothing;

Depends:
1. LinearAlgebra.diagind()

Note:
1. diagind(MAT) (the index of positions to fill) can have different length from vector VEC
   1. if length(diagind(MAT)) < length(VEC): extra elements of VEC will be neglected
   2. if length(diagind(MAT)) > length(VEC): extra elements of MAT[diagind(MAT)] will be kept
"""
function diagw!( MAT::Matrix, VEC::Union{Vector,Tuple,NamedTuple,AbstractRange} ; offset::Int = 0 )
   # declare
   local tmpIdxRange = diagind(MAT, offset) # get indexes range
   # fill elements
   for s in 1:min( length(VEC), length(tmpIdxRange) ); MAT[ tmpIdxRange[s] ] = VEC[s]; end
   return nothing::Nothing
end # end function diagw
# ----------------
"""
   expand( VEC::Vector, LEN::Int )

expands VEC to specific length LEN;
if LEN > length(VEC), fill extra positions with the last element of VEC;
if LEN < length(VEC), drop the ending elements;
returns a Vector/1-D Array of the same type
"""
function expand( VEC::Vector, LEN::Int )
   LEN >= 0  ||  "expects a non-negative length to expand but received: $(LEN)"
   tmpLen = length(VEC); tmpVec = copy(VEC)
   if LEN > tmpLen
      tmpVal = tmpVec[end]
      for s in range( tmpLen+1, length= LEN-tmpLen )
         push!( tmpVec, tmpVal )
      end
      return tmpVec::Vector
   else
      return tmpVec[1:LEN]::Vector
   end
end # end function expand
function expand!( VEC::Vector, LEN::Int )
   LEN >= 0  ||  "expects a non-negative length to expand but received: $(LEN)"
   tmpLen = length(VEC)
   if LEN > tmpLen
      tmpVal = VEC[end]
      for s in range( tmpLen + 1, length= LEN - tmpLen )
         push!( VEC, tmpVal )
      end
      return nothing::Nothing
   else
      splice!( VEC, (LEN + 1):tmpLen ) # similar to pop!() but on a range
      return nothing::Nothing
   end
end # end function expand!
# ----------------

end # end module EasyIO

max_stars_repo_name: Clpr/NumToolLibs

max_issues_repo_name: Clpr/NumToolLibs

max_forks_repo_name: Clpr/NumToolLibs


text:
using Knet, Base.Test

# x* = f(1, · · · , 1)
# f(x*) = 0

rosenbrock(x) = sum((1-x[1:end-1]).^2 + 100*(x[2:end]-x[1:end-1].^2).^2)

grads = grad(rosenbrock)
srand(123456789)

function test_base(w, params)
	prev = 0
	i = 1

	current = rosenbrock(w)
	while i <= 1000000 && abs(current - prev) > 1e-10 && current > 1e-3
		g = grads(w)
		w, params = update!(w, g, params)
		prev = current
		current = rosenbrock(w)

		i += 1
	end

	@test current <= 1e-3
	info("$(typeof(params)) Passed!\nConverged to $(current) at epoch $(i-1)")
end

function test_sgd(w)
	params = Sgd(;lr=0.0005)
	test_base(w, params)
end

function test_momentum(w)
	params = Momentum(w;lr=0.00005, gamma=0.95, velocity=zeros(w))
	test_base(w, params)
end

function test_adam(w)
	params = Adam(w; lr=0.005, beta1=0.9, beta2=0.95, t=1, eps=1e-8, fstm=zeros(w), scndm=zeros(w))
	test_base(w, params)
end

function test_adagrad(w)
	params = Adagrad(w; lr=0.35, eps=1e-6, G=zeros(w))
	test_base(w, params)
end

function test_adadelta(w)
	params = Adadelta(w; lr=0.001, rho=0.9, eps=1e-6, G=zeros(w), delta=zeros(w))
	test_base(w, params)
end

function test_rmsprop(w)
	params = Rmsprop(w; lr=0.0002, rho=0.9, eps=1e-6, G=zeros(w))
	test_base(w, params)
end

#GPU Tests
info("GPU Tests")
raw = randn(100, 1)
w = convert(KnetArray, raw)
@time test_sgd(copy(w))
@time test_momentum(copy(w))
@time test_adam(copy(w))
@time test_adagrad(copy(w))
@time test_adadelta(copy(w))
@time test_rmsprop(copy(w))

#CPU Tests
gpu(-1)
info("CPU Tests")
w = randn(100, 1)
@time test_sgd(copy(w))
@time test_momentum(copy(raw))
@time test_adam(copy(raw))
@time test_adagrad(copy(raw))
@time test_adadelta(copy(raw))
@time test_rmsprop(copy(raw))


text:
mutable struct Moments
    μ
    Σ
    Moments(ndims) = new(zeros(Float64, ndims), diagm(fill(1.0, ndims)))
end

function updateDN!(pind, nind, post, gh, pmoments)
    G = size(pind, 1)
    for g in axes(pind, 1)
        N = 0
        d = @views @fastmath @inbounds sum([post[i, :] .* nind[i] for i in pind[g]])
        d = d ./ sum(d)
        m = gh[g].n'd
        deviance = gh[g].n .- m'
        s = d' * deviance.^2
        pmoments[g].μ = m
        pmoments[g].Σ = s
        if g > 1
            gh[g].w = map(x -> pdf(MvNormal(m, sqrt.(s)), x), eachrow(gh[g].n))
            gh[g].w ./= sum(gh[g].w)
        end
    end
end


text:
module TKMPMeasures

export integrate
export ∫
export TKMPMeasure

include("tkmp_measure.jl")

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

end # module


text:
# The SteinMaternPowerKernel is given by
#
# k(x, y) = ((1 + ||x||^alpha)*(1 + ||y||^alpha))^-1 * km(x,y)
#
# where km is the matern kernel.

mutable struct SteinMaternPowerKernel <: SteinMaternWeightedKernel
    # alpha is the exponent on the norms
    alpha::Float64
    # gamma is the rate parameter for the matern kernel
    gamma::Float64
end

# default is use gamma = sqrt{3} and alpha = 2.5
SteinMaternPowerKernel() = SteinMaternPowerKernel(2.5, sqrt(3.0))

# default is use gamma sqrt{3}
SteinMaternPowerKernel(alpha::Float64) = SteinMaternPowerKernel(alpha, sqrt(3.0))

# The term (1 + ||z||^alpha)^-1
function _Qz(ker::SteinMaternPowerKernel, z::Array{Float64, 1})
    Qz = 1.0 + norm(z,2)^ker.alpha
    1.0 / Qz
end

function Q(ker::SteinMaternPowerKernel, x::Array{Float64,1}, y::Array{Float64,1})
    _Qz(ker, x) * _Qz(ker, y)
end

function gradxQ(ker::SteinMaternPowerKernel, x::Array{Float64,1}, y::Array{Float64,1})
    fx = _Qz(ker, x)
    fy = _Qz(ker, y)
    mx = norm(x, 2)^(ker.alpha - 2)
    -x * ker.alpha * fx^2 * fy * mx
end

function gradxyQ(ker::SteinMaternPowerKernel, x::Array{Float64,1}, y::Array{Float64,1})
    mx = norm(x, 2)^(ker.alpha - 2.0)
    my = norm(y, 2)^(ker.alpha - 2.0)
    dot(x, y) * ker.alpha^2 * mx * my * Q(ker, x, y)^2
end