text:
function place_first_grid_order!(current, ticker, config, open_orders, bot_path, client, usr)
    @info("NEW order found!")
    id_str = open_orders["NEW"]["trigger"] * "--" * string(now())
    quant = ceil(config["order_volume"] / current, digits=config["precision"])
    if config["precision"] < 1
        quant = string(Int(quant))
    end
    order = try client.order_limit_buy(
        symbol=config["symbol"],
        quantity=quant,
        price=ticker);
    catch
        @info("order_limit_buy failed")
        nothing
    end
    if order !== nothing
        open_orders["BUY"]["EXTRA-0_" * id_str] = order
        open_orders["NEW"] = Dict()
        open(joinpath(bot_path, "open_orders.json"), "w") do f
            JSON.print(f, open_orders, 4)
        end
        p = Dict()
        p["ticker"] = ticker
        create_notification("new cycle", config, usr, params = p)
    end
    return open_orders
end
export place_first_grid_order!

function place_buy_order_command(k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees)
    quant = round(v["origQty"], digits=config["precision"])
    new_id = k
    if v["price"] != "market"
        new_price = Pinguu.price2string(parse(Float64, v["price"]), config["sell_precision"])
        order = try client.order_limit_buy(
            symbol= config["symbol"],
            quantity=quant,
            price=new_price);
            catch
            nothing
        end
    else
        order = try client.order_market_buy(
            symbol= config["symbol"],
            quantity=quant);
            catch
            nothing
        end
    end
    if order !== nothing
        order["time"] = string(now())
        open_orders["BUY"][new_id] = order
    else
        open_orders["BUY"][new_id] = Dict()
        open_orders["BUY"][new_id]["origQty"] = quant
        open_orders["BUY"][new_id]["price"]   = new_price
        open_orders["BUY"][new_id]["status"]  = "place order"
        open_orders["BUY"][new_id]["time"]    = string(now())
    end
    open(joinpath(bot_path, "open_orders.json"), "w") do f
        JSON.print(f, open_orders, 4)
    end
    open(joinpath(bot_path, "trigger.json"), "w") do f
        JSON.print(f, Dict(), 4)
    end
end
export place_buy_order_command

function buy_executed(tmp, k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees)
    #
    ## Cancel order, if not done yet
    if tmp["status"] != client.ORDER_STATUS_FILLED && tmp["status"] != client.ORDER_STATUS_CANCELED
        result = try client.cancel_order(
            symbol=config["symbol"],
            orderId=tmp["orderId"])
        catch
            nothing
        end
        if result == nothing
            return nothing
        end
    end
    extra_iteration = parse(Int64, split(split(k, "_")[1], "-")[2])
    #
    ## Move to closed_orders
    !(id_str in keys(closed_orders)) ? closed_orders[id_str] = Dict() : ""
    tmp["time"] = string(now())
    closed_orders[id_str][split(k, "_")[1]] = deepcopy(tmp);
    open(joinpath(bot_path, "closed_orders.json"), "w") do f
        JSON.print(f, closed_orders, 4)
    end
    #
    ## Remove from open_orders
    delete!(open_orders["BUY"], k)
    open(joinpath(bot_path, "open_orders.json"), "w") do f
        JSON.print(f, open_orders, 4)
    end
    #
    ## Place EXTRA
    new_id = "EXTRA-" * string(extra_iteration + 1) * "_" * id_str

    last_price = parse(Float64, tmp["price"]) == 0.0 ? v["price"] : tmp["price"]

    new_price = Pinguu.price2string(parse(Float64, last_price) * (1 - config["extra_order_step"]), config["sell_precision"], round_up=true)
    quant = floor(config["order_volume"] * (config["extra_order_martingale"])^(extra_iteration + 1) / parse(Float64, new_price), digits=config["precision"])
    if extra_iteration < parse(Float64, config["extra_order_count"])
        open_orders["BUY"][new_id] = Dict()
        open_orders["BUY"][new_id]["origQty"] = quant
        open_orders["BUY"][new_id]["price"]   = new_price
        open_orders["BUY"][new_id]["status"]  = "place order"
        open_orders["BUY"][new_id]["time"]    = string(now())
        open(joinpath(bot_path, "open_orders.json"), "w") do f
            JSON.print(f, open_orders, 4)
        end
        open(joinpath(bot_path, "trigger.json"), "w") do f
            JSON.print(f, Dict(), 4)
        end
    end
    #
    ## Cancel current take profit orders
    if length(open_orders["SELL"]) > 0
        for (k2, v2) in open_orders["SELL"]
            if split(k2, "_")[end] == id_str
                tmp_sell = try client.get_order(
                    symbol=config["symbol"],
                    orderId=v2["orderId"])
                    catch
                    try sleep(0.5)
                        client.get_order(
                            symbol=config["symbol"],
                            orderId=v2["orderId"])
                        catch nothing
                    end
                end
                if !(tmp_sell["status"] in [client.ORDER_STATUS_CANCELED, client.ORDER_STATUS_FILLED]) || tmp_sell == nothing
                    sleep(0.5)
                    tmp_sell = try client.cancel_order(
                        symbol=config["symbol"],
                        orderId=v2["orderId"])
                    catch
                        nothing
                    end
                end
                if tmp_sell !== nothing
                    if parse(Float64, v2["executedQty"]) > 0
                        random_id_str = "SELL-" * lpad(rand(1:1:10000), 5, "0")
                        while haskey(closed_orders[id_str], random_id_str)
                            random_id_str = "SELL-" * lpad(rand(1:1:10000), 5, "0")
                        end
                        closed_orders[id_str][random_id_str] = result
                        open(joinpath(bot_path, "closed_orders.json"), "w") do f
                            JSON.print(f, closed_orders, 4)
                        end
                    end
                    delete!(open_orders["SELL"], k2)
                else
                    #
                    ## Cancel failed
                    delete!(open_orders["SELL"], k2)
                    v2["status"] = "cancel order"
                    #
                    ## Give different ID to be able to create a new "main" order later
                    random_id_str = "SELL-" * lpad(rand(1:1:10000), 5, "0")
                    open_orders[random_id_str] = v2
                end
            end
        end
    end

    #
    ## Place command for take profit order
    open_orders["SELL"]["SELL_" * id_str] = Dict()
    open_orders["SELL"]["SELL_" * id_str]["status"] = "place order"
    open_orders["SELL"]["SELL_" * id_str]["id_str"] = id_str

    open(joinpath(bot_path, "open_orders.json"), "w") do f
        JSON.print(f, open_orders, 4)
    end
    open(joinpath(bot_path, "trigger.json"), "w") do f
        JSON.print(f, Dict(), 4)
    end

    #
    ## Create notification
    p = Dict()
    p["extra_iteration"] = extra_iteration
    p["executedQty"] = tmp["executedQty"]
    create_notification("buy exec", config, usr, params = p)
    if extra_iteration >= parse(Float64, config["extra_order_count"])
        create_notification("no volume", config, usr)
    end
end
export buy_executed

function cancel_buy(k, v, id_str, config, open_orders, closed_orders, bot_path, client)
    @info("Cancel Order")
    if haskey(v, "orderId")
        result = try client.cancel_order(
            symbol=config["symbol"],
            orderId=v["orderId"])
        catch
            nothing
        end
        if result !== nothing
            delete!(open_orders["BUY"], k)
            open(joinpath(bot_path, "open_orders.json"), "w") do f
                JSON.print(f, open_orders, 4)
            end
            if parse(Float64, result["executedQty"]) > 0
                random_id_str = "EXTRA-" * lpad(rand(100:1:10000), 5, "0")
                while haskey(closed_orders[id_str], random_id_str)
                    random_id_str = "EXTRA-" * lpad(rand(100:1:10000), 5, "0")
                end
                closed_orders[id_str][random_id_str] = result
                open(joinpath(bot_path, "closed_orders.json"), "w") do f
                    JSON.print(f, closed_orders, 4)
                end
            end
        end
    else
        @info("Order has not been placed")
        #
        ## Order not placed, yet
        delete!(open_orders["BUY"], k)
        open(joinpath(bot_path, "open_orders.json"), "w") do f
            JSON.print(f, open_orders, 4)
        end

    end
    open(joinpath(bot_path, "trigger.json"), "w") do f
        JSON.print(f, Dict(), 4)
    end
end
export cancel_buy

function buy_now(tmp, k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees)
    extra_iteration = parse(Int64, split(split(k, "_")[1], "-")[2])
    if extra_iteration < parse(Float64, config["extra_order_count"])
        #
        # Cancel current BUY order
        if tmp["status"] != client.ORDER_STATUS_CANCELED
            result = try client.cancel_order(
                symbol=config["symbol"],
                orderId=v["orderId"])
            catch
                nothing
            end
            if result == nothing
                return nothing
            end
        else
            result = tmp
        end

        #
        # Check if that order was partially filled -> if so: move cancelled order to closed_orders
        if parse(Float64, result["executedQty"]) > 0
            result["time"] = string(now())
            !haskey(closed_orders, id_str) ? closed_orders[id_str] = Dict() : ""
            closed_orders[id_str][split(k, "_")[1]] = result
            extra_iteration += 1
            open(joinpath(bot_path, "closed_orders.json"), "w") do f
                JSON.print(f, closed_orders, 4)
            end
        end
        #
        # Remove from open_orders
        delete!(open_orders["BUY"], k)
        #
        # Get new order details
        ticker = client.get_ticker(symbol=config["symbol"])["lastPrice"]
        quant = floor(parse(Float64, v["origQty"]), digits=config["precision"])
        new_id = "EXTRA-" * string(extra_iteration) * "_" * id_str
        #
        # Place new order
        order = try client.order_limit_buy(
            symbol=config["symbol"],
            quantity=quant,
            price=ticker)
        catch
            nothing
        end

        if order !== nothing
            order["time"] = string(now())
#                 if parse(Float64, order["price"]) == 0
#                     order["price"] = Pinguu.price2string(parse(Float64, order["cummulativeQuoteQty"]) / parse(Float64, order["origQty"]), config["sell_precision"])
#                 end
            open_orders["BUY"][new_id] = order
        else
            open_orders["BUY"][new_id] = Dict()
            open_orders["BUY"][new_id]["origQty"] = quant
            open_orders["BUY"][new_id]["price"]   = ticker
            open_orders["BUY"][new_id]["status"]  = "place order"
            open_orders["BUY"][new_id]["time"]    = string(now())
        end
        open(joinpath(bot_path, "open_orders.json"), "w") do f
            JSON.print(f, open_orders, 4)
        end
    end
    open(joinpath(bot_path, "trigger.json"), "w") do f
        JSON.print(f, Dict(), 4)
    end
end
export buy_now

function partially_order(tmp, k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees, minutes)
    if v["partial_time"] < string(now() - Minute(minutes))
        result = try client.cancel_order(
            symbol=config["symbol"],
            orderId=tmp["orderId"])
        catch
            nothing
        end
        if result !== nothing
            tmp["partial_time"] = v["partial_time"]
            v = deepcopy(tmp)
            v["status"] = "continue"
            open_orders["BUY"][k] = v
            open(joinpath(bot_path, "open_orders.json"), "w") do f
                JSON.print(f, open_orders, 4)
            end
            p = Dict()
            p["extra_iteration"] = extra_iteration
            p["executedQty"] = tmp["executedQty"]
            create_notification("partial buy exec", config, usr, params = p)
        end
        if tmp["status"] == client.ORDER_STATUS_CANCELED
            v["status"] = "continue"
            open_orders["BUY"][k] = v
            open(joinpath(bot_path, "open_orders.json"), "w") do f
                JSON.print(f, open_orders, 4)
            end
        end
    end
end
export partially_order

function first_order_takes_too_long(tmp, k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees, minutes)
    t = DateTime(split(id_str, "--")[end])
    if now() - t > Minute(minutes)
        result = try client.cancel_order(
            symbol=config["symbol"],
            orderId=v["orderId"])
        catch
            nothing
        end
        if result !== nothing
            delete!(open_orders["BUY"], k)
            open(joinpath(bot_path, "open_orders.json"), "w") do f
                JSON.print(f, open_orders, 4)
            end
            Pinguu.create_notification("cycle closed", config, usr)
        end
    end
end
export first_order_takes_too_long

function add_fees_to_closed_order(order, config, bot_path, client)
    for t in order["trades"]
        sleep(0.5)
        commission = parse(Float64, t["commission"])
        commission_coin = parse(Float64, t["commission"])
        if commission != 0
            from = split(string(unix2datetime(t["time"]/1000)), ".")[1]
            to = split(string(unix2datetime(t["time"]/1000) + Minute(1)), ".")[1]
            #
            ## for commission
            if t["commissionAsset"] == config["farm_coin"]
                commission *= parse(Float64, t["price"])
            elseif t["commissionAsset"] != config["base_coin"]
                translator = string(t["commissionAsset"], config["base_coin"])
                klines = client.get_historical_klines(translator, client.KLINE_INTERVAL_1MINUTE, from, to)
                tt = Pinguu.klines2table(klines, time_warp=0)
                commission *= tt[1].close
                tmp_price_for_conversion = tt[1].close
            end
            #
            ## for commission_coin
            if t["commissionAsset"] != config["farm_coin"]
                translator = string(t["commissionAsset"], config["farm_coin"])
                # Translator exists
                klines = try client.get_historical_klines(translator, client.KLINE_INTERVAL_1MINUTE, from, to)
                    catch
                    false
                end
                if klines == false
                    # Translator exists in other order
                    translator = string(config["farm_coin"], t["commissionAsset"])
                    klines = try client.get_historical_klines(translator, client.KLINE_INTERVAL_1MINUTE, from, to)
                        catch
                        false
                    end
                else
                    tt = Pinguu.klines2table(klines, time_warp=0)
                    weightedAvgPrice = 1/tt[1].close
                end
                if klines == false
                    # Translator does not exist e.g. DOGEBNB
                    weightedAvgPrice = tmp_price_for_conversion / parse(Float64, t["price"])
                else
                    tt = Pinguu.klines2table(klines, time_warp=0)
                    weightedAvgPrice = 1/tt[1].close
                end

                if typeof(weightedAvgPrice) == String
                    weightedAvgPrice = parse(Float64, weightedAvgPrice)
                end
                commission_coin *= weightedAvgPrice
            end
            t["com"] = commission
            t["com_coin"] = commission_coin
        else
            t["com"] = 0
            t["com_coin"] = 0
        end
    end
    return order
end
export add_fees_to_closed_order

function check_ongoing_closed(id_str, config, closed_orders, bot_path, client)
    store = false
    for (id, order) in closed_orders[id_str]
        if !haskey(order, "trades")
            trades = client.get_my_trades(orderId=order["orderId"], symbol=order["symbol"])
            order["trades"] = trades
            order = try add_fees_to_closed_order(order, config, bot_path, client)
                catch
                @info("Getting fees failed.")
                return
            end
            store = true
        end
    end
    if store
        open(joinpath(bot_path, "closed_orders.json"), "w") do f
            JSON.print(f, closed_orders, 4)
        end
    end
end
export check_ongoing_closed


text:
#=
This file contains the helper functions necessary for file import/export
of various file formats.
=#

"""
    determine_marker(geno_parse::T, digits::Int) where T<:AbstractDataFrame
Return either `Int8` or `Int16` depending on largest allelic value in all genotypes
in the first 10 samples of an input DataFrame (or all the samples if less than 10 samples).
If the largest allele is 20 or greater, the marker will be considered a Microsatellite
and coded in `PopData` as `Int16`, and the opposite is true for SNPs. There's no
specific reason 100 was chosen other than it being a reasonable buffer for edge
cases since SNP data <= 4, and haplotyped data could be a bit higher. Even if the
microsatellite markers are coded incorrectly, there will be zero impact to performance,
and considering how few microsatellite markers are used in typical studies, the
effect on in-memory size should be negligible (as compared to SNPs).
"""
function determine_marker(geno_parse::T, digits::Int) where T<:AbstractDataFrame
    # get the total # columns
    total_col = size(geno_parse,2)
    # find the 25% cutoff
    if total_col > 11
        num_test_loc = (total_col - 2) ÷ 8
    else
        num_test_loc = total_col - 2
    end
    # remove everything else
    test_df = @view geno_parse[:, 3:num_test_loc]

    # isolate the largest allele value
    max_allele = map(eachcol(test_df)) do i
        phase.(i, Int16, digits)  |>
        skipmissing |> Base.Iterators.flatten |> maximum
    end |> maximum

    if max_allele <= 100
        return Int8
    else
        return Int16
    end
end

"""
    find_ploidy(genotypes::T where T<:SubArray)
Used internally in the `genepop` and `delimited` file parsers to scan the genotypes
of a sample and return the ploidy of the first non-missing locus.
"""
@inline function find_ploidy(genotypes::T) where T<:GenoArray
    return Int8(length(first(skipmissing(genotypes))))
end

"""
    phase(loc::Union{String, Int}, type::DataType, digit::Int)
Takes a String of numbers or Integers and returns a typed locus appropriate for PopGen.jl as used in the
`genepop` and `csv` file parsers. Use `type` to specify output type (`Int8` or `Int16`),
and `digit` to specify the number of digits/characters used per allele in a locus.

## Examples
```
ph_locus = phase("128114", Int16, 3)
map(i -> phase(i, Int16, 3), ["112131", "211112", "001003", "516500"])
# or #
[phase(i, Int8, 2) for i in ["0101", "0103", "0202", "0103"]]
```
"""
@inline function phase(loc::T, type::DataType, digit::Int) where T<:AbstractString
    loc == "-9" || iszero(parse(Int, loc)) && return missing
    phased = map(i -> parse(type, join(i)), Iterators.partition(loc, digit))
    sort!(phased)
    Tuple(phased)
end

phase(loc::Missing, type::DataType, digit::Int) = missing

@inline function phase(loc::T, type::DataType, digit::Int) where T<:Integer
    loc == -9 || iszero(loc) && return missing
    units = 10^digit
    allele1 = loc ÷ units |> type
    allele2 = loc % units |> type
    return [allele1, allele2] |> sort |> Tuple
end

"""
    unphase(geno::T; digits::Int = 3, ploidy::Int = 2, miss::Int = 0) where T <: Genotype
Takes a `Genotype` (e.g. `(131, 94)`) and returns a string of concatenated
alleles padded with *n* number of zeroes, where *n* is given by `digits = `.
`missing` values are returned as either a string of 'digits × ploidy' zeroes (`miss = 0`)
or `"-9"` (`miss = -9`). The `ploidy` flag is only relevant for unphasing `missing` genotypes
and not used otherwise.

### Example
```
unphase((1,2,3,4), digits = 3)
"001002003004"

unphase(missing, digits = 2, ploidy = 2, miss = -9)
"-9"

unphase(missing, digits = 2, ploidy = 2, miss = 0)
"0000"
```
"""
function unphase(geno::T; digits::Int = 3, ploidy::Int = 2, miss::Int = 0) where T <: Genotype
    join(map(i -> lpad(i, digits, "0"), geno))
end

function unphase(geno::Missing; digits::Int = 3, ploidy::Int, miss::Int = 0)
    if miss == 0
        return "0"^(digits * ploidy)
    else
        return "$miss"
    end
end

ext: jl


text:
function generate(rng::AbstractRNG,
                  ::Type{Transmission},
                  tr::TransmissionRates,
                  event::Event{T}) where T <: DiseaseStateSequence
  id = event.individual
  if _new_transmission(event)
    external_or_internal = Weights([tr.external[id]; sum(tr.internal[:,id])])
    if sum(external_or_internal) == 0.0
      @error "All transmission rates = 0.0, No transmission can be generated"
      return NoTransmission()
    elseif sample(rng, [true; false], external_or_internal)
      @debug "Exogenous tranmission generated"
      return ExogenousTransmission(id)
    else
      source = sample(rng, 1:individuals(tr), Weights(tr.internal[:, id]))
      @debug "Endogenous transmission generated (source id = $source)"
      return EndogenousTransmission(id, source)
    end
  else
    @debug "No transmission generated"
    return NoTransmission()
  end
end

function generate(rng::AbstractRNG,
                  ::Type{Transmission},
                  tr::TransmissionRates,
                  event::NoEvent)
  return NoTransmission()
end


text:
using Dates
using Distributions

#start_time = now()
#end_time = now()
@testset "Black Scholes with Banyan with size=$size" for size in [128_000_000, 256_000_000, 512_000_000]
    use_job_for_testing(sample_rate = 128) do

        # NOTE: 64M works but 128M doesn't (and also doesn't batch the I/O)
        price = BanyanArrays.fill(4.0, size)
        strike = BanyanArrays.fill(4.0, size)
        t = BanyanArrays.fill(4.0, size)
        rate = BanyanArrays.fill(4.0, size)
        vol = BanyanArrays.fill(4.0, size)

        d1 = map(
            (p, s, t, r, v) -> ((log(p / s) + (r + v^2 * 0.5) * t) / (v * sqrt(t))),
            price,
            strike,
            t,
            rate,
            vol,
        )

        d2 = map((d1, v, t) -> (d1 - (v * sqrt(t))), d1, vol, t)

        call = map(
            (d1, d2, p, s, t, r) ->
                ((cdf(Normal(), d1) * p) - (cdf(Normal(), d2) * s * exp(-r * t))),
            d1,
            d2,
            price,
            strike,
            t,
            rate,
        )

        call_sum = sum(call)
        call = nothing
        res = collect(call_sum)
        # println(res)
        @show res
        #res = collect(res)
        #end = now()
        #@test typeof(res) == Base.Vector{Float64}
        #@test all(v->v==3.999999985812889, res)
    end
end
#end_time = now()
#println(end_time - start_time)

@testset "Black Scholes without Banyan with size=$size" for size in [128_000, 128_000_000, 512_000_000]
    use_job_for_testing(sample_rate = 128) do

        price = Base.fill(4.0, size)
        strike = Base.fill(4.0, size)
        t = Base.fill(4.0, size)
        rate = Base.fill(4.0, size)
        vol = Base.fill(4.0, size)

        d1 = Base.map(
            (p, s, t, r, v) -> ((log(p / s) + (r + v^2 * 0.5) * t) / (v * sqrt(t))),
            price,
            strike,
            t,
            rate,
            vol,
        )

        d2 = Base.map((d1, v, t) -> (d1 - (v * sqrt(t))), d1, vol, t)

        call = Base.map(
            (d1, d2, p, s, t, r) ->
                ((cdf(Normal(), d1) * p) - (cdf(Normal(), d2) * s * exp(-r * t))),
            d1,
            d2,
            price,
            strike,
            t,
            rate,
        )

        call_sum = Base.sum(call)
        res = Base.collect(call_sum)
        println(res)
        # @test typeof(res) == Base.Vector{Float64}
        # @test all(v->v==3.999999985812889, res)
    end
end

@testset "Black Scholes without Banyan vectorized with size=$size" for size in [
    128_000,
    128_000_000,
    512_000_000,
]
    use_job_for_testing(sample_rate = 128) do

        price = Base.fill(4.0, size)
        strike = Base.fill(4.0, size)
        t = Base.fill(4.0, size)
        rate = Base.fill(4.0, size)
        vol = Base.fill(4.0, size)

        d1 = (log.(price ./ strike) .+ (rate .+ vol .^ 2 .* 0.5) .* t) ./ (vol .* sqrt.(t))

        d2 = d1 .- (vol .* sqrt.(t))

        call =
            (cdf.(Normal(), d1) .* price) - (cdf.(Normal(), d2) .* strike .* exp.(-rate .* t))

        call_sum = Base.sum(call)
        res = Base.collect(call_sum)
        println(res)
    end
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
# This file calculates the first application in Y.-S. Chan, A. C. Fannjiang, and G. H. Paulino,
# Integral equations with hypersingular kernels -- theory and applications to fracture mechanics,
# Int. J. Eng. Sci., 41:683--720, 2003.

using ApproxFun, SingularIntegralEquations

x = Fun(identity)
w = 1/sqrt(1-x^2)
d = domain(x)
d2 = d^2
B = dirichlet(d)
H2 = Hilbert(d,2)
Σ = DefiniteIntegral(d)
ϵ = 2.0
K = LowRankFun((x,y)->-1./(x+y+2ϵ).^2+12(x+ϵ)./(x+y+2ϵ).^3-12(x+ϵ).^2./(x+y+2ϵ).^4,d2)
L = H2[w] + Σ[K*(w/π)]
f = -Fun(one)

uSIE = [B;L]\[zeros(2);f]
@time uSIE = [B;L]\[zeros(2);f]
println("The rank of K is: ",rank(K))
println("The length of uSIE is: ",length(uSIE))
println("The extrema of uSIE are: ",extrema(uSIE))
println("The normalized Stress Intensify Factors are: ",(uSIE/(1-x^2))(-1),"  ",(uSIE/(1-x^2))(1))


text:
include("../src/utilities.jl")

@testset "Plots request massaging" begin
    wanted = "all"
    needed = Set([:fit, :residuals, :normal_checks, :cooksd, :leverage, :homoscedasticity])
    @test needed == get_needed_plots(wanted)

    wanted = []
    needed = Set([])
    @test needed == get_needed_plots(wanted)

    wanted = :none
    needed = Set()
    @test needed == get_needed_plots(wanted)

    wanted = Set()
    needed = Set()
    @test needed == get_needed_plots(wanted)

    wanted = ["fit", "cooksd", "homoscedasticity"]
    needed = Set([:fit, :cooksd, :homoscedasticity])
    @test needed == get_needed_plots(wanted)

    wanted = [:residuals, :normal_checks]
    needed = Set([:residuals, :normal_checks])
    @test needed == get_needed_plots(wanted)

    wanted = :leverage
    needed = Set([:leverage])
    @test needed == get_needed_plots(wanted)

end

@testset "Covariance estimator stats massaging" begin
    wanted = [:white, :white , :bogus , :nw]
    needed = ([:white], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = [:hc0, :hc1, :nw]
    needed = ([:hc0, :hc1], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = ["White", "HC0" , "nw" , "nothing"]
    needed = ([:white, :hc0], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = ["Hc1", "HC0" , "nothing"]
    needed = ([:hc0 , :hc1], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = "Nw"
    needed = ([], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = :hc0
    needed = ([:hc0], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = :hc0
    needed = ([:hc0], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = []
    needed = ([], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = :none
    needed = ([], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = :all
    needed = ([:white, :hc0, :hc1, :hc2, :hc3], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

end

@testset "model stats massaging" begin
    wanted = ["default"]
    needed = Set([:coefs, :sse, :mse, :sst, :rmse, :sigma, :t_statistic, :r2, :adjr2, :stderror, :t_values, :p_values, :ci])
    @test needed == get_needed_model_stats(wanted)
    wanted = [:default, :vif]
    needed = Set([:coefs, :sse, :vif, :mse, :sst, :rmse, :sigma, :t_statistic, :r2, :adjr2, :stderror, :t_values, :p_values, :ci])
    @test needed == get_needed_model_stats(wanted)
    wanted = [:default, :diag_ks]
    needed = Set([:coefs, :sse, :diag_ks, :mse, :sst, :rmse, :sigma, :t_statistic, :r2, :adjr2, :stderror, :t_values, :p_values, :ci])
    @test needed == get_needed_model_stats(wanted)
    wanted = [:default, :diag_normality]
    needed = Set([:coefs, :sse, :diag_ks, :diag_ad, :diag_jb, :mse, :sst, :rmse, :sigma, :t_statistic, :r2, :adjr2, :stderror, :t_values, :p_values, :ci])
    @test needed == get_needed_model_stats(wanted)
    wanted = [:default, :diag_heteroskedasticity]
    needed = Set([:coefs, :sse, :diag_white, :diag_bp, :mse, :sst, :rmse, :sigma, :t_statistic, :r2, :adjr2, :stderror, :t_values, :p_values, :ci])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["r2", "rmse"]
    needed = Set([:r2, :rmse, :coefs, :sse, :mse, :sst])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["aic"]
    needed = Set([:coefs, :mse, :sse, :aic])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["Ci"]
    needed = Set([:coefs, :mse, :sse, :ci, :sigma, :stderror, :t_statistic])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["Adjr2"]
    needed = Set([:coefs, :mse, :sse, :r2, :adjr2, :sst])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["stderror"]
    needed = Set([:coefs, :mse, :sse, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["t_values"]
    needed = Set([:coefs, :mse, :sse, :t_values, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["p_values"]
    needed = Set([:coefs, :mse, :sse, :p_values, :t_values, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["ci", "P_values"]
    needed = Set([:coefs, :mse, :sse, :p_values, :t_values, :stderror, :sigma, :ci, :t_statistic])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["none"]
    needed = Set([:coefs, :mse, :sse])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["all"]
    needed = Set([:coefs, :sse, :mse, :sst, :rmse, :aic, :sigma, :t_statistic, :vif, :r2, :adjr2, :stderror, :t_values, :p_values, :ci,
                    :diag_normality, :diag_ks, :diag_ad, :diag_jb, :diag_heteroskedasticity, :diag_white, :diag_bp, :press,
                    :t1ss, :t2ss, :pcorr1, :pcorr2, :scorr1, :scorr2])
    @test needed == get_needed_model_stats(wanted)
    wanted = [ ]
    needed = Set([:coefs, :mse, :sse])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["press" ]
    needed = Set([:coefs, :mse, :sse, :press])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["t1ss" ]
    needed = Set([:coefs, :mse, :sse, :t1ss])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["t2ss", "t1ss" ]
    needed = Set([:coefs, :mse, :sse, :t1ss, :t2ss])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["pcorr1" ]
    needed = Set([:coefs, :mse, :sse, :t1ss, :pcorr1])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["pcorr1", "pcorr2" ]
    needed = Set([:coefs, :mse, :sse, :t1ss, :t2ss, :pcorr1, :pcorr2])
    @test needed == get_needed_model_stats(wanted)
    wanted = [:scorr1 ]
    needed = Set([:coefs, :mse, :sse, :sst, :t1ss, :scorr1])
    @test needed == get_needed_model_stats(wanted)
    wanted = [:scorr1, :scorr2 ]
    needed = Set([:coefs, :mse, :sse, :sst, :t1ss, :t2ss, :scorr1, :scorr2])
    @test needed == get_needed_model_stats(wanted)
    wanted = [ "bogus"]
    needed = Set([:coefs, :mse, :sse])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["stderror", "Bogus"]
    needed = Set([:coefs, :mse, :sse, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = Set([:stderror, :Bogus])
    needed = Set([:coefs, :mse, :sse, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = Set()
    needed = Set([:coefs, :mse, :sse])
    @test needed == get_needed_model_stats(wanted)
end

@testset "prediction stats massaging" begin
    wanted = ["none"]
    target_need = Set([:predicted])
    target_present = Set([:predicted])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["all"]
    target_need = Set([:predicted, :residuals, :leverage, :stdp, :stdi, :stdr, :student, :rstudent, :lcli, :ucli, :lclp, :uclp, :press, :cooksd])
    target_present = Set([:predicted, :residuals, :leverage, :stdp, :stdi, :stdr, :student, :rstudent, :lcli, :ucli, :lclp, :uclp, :press, :cooksd])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["CooksD"]
    target_need = Set([:predicted, :residuals, :leverage, :stdp, :stdr, :student, :cooksd])
    target_present = Set([:predicted, :cooksd])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["CooksD", "lcli"]
    target_need = Set([:predicted, :residuals, :leverage, :stdp, :stdr, :student, :cooksd, :stdi, :lcli])
    target_present = Set([:predicted, :cooksd, :lcli])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["CooksD", "Bogus", "lcli"]
    target_need = Set([:predicted, :residuals, :leverage, :stdp, :stdr, :student, :cooksd, :stdi, :lcli])
    target_present = Set([:predicted, :cooksd, :lcli])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["PRESS"]
    target_need = Set([:predicted, :residuals, :leverage, :press])
    target_present = Set([:predicted, :press])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["lclp"]
    target_need = Set([:predicted, :stdp, :leverage, :lclp])
    target_present = Set([:predicted, :lclp])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["uclp"]
    target_need = Set([:predicted, :stdp, :leverage, :uclp])
    target_present = Set([:predicted, :uclp])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["ucli"]
    target_need = Set([:predicted, :stdi, :leverage, :ucli])
    target_present = Set([:predicted, :ucli])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["ucli", "LCLP"]
    target_need = Set([:predicted, :stdi, :leverage, :ucli, :lclp, :stdp])
    target_present = Set([:predicted, :ucli, :lclp])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["student"]
    target_need = Set([:predicted, :residuals, :leverage, :stdr, :student])
    target_present = Set([:predicted, :student])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = [:rstudent]
    target_need = Set([:predicted, :residuals, :leverage, :stdr, :student, :rstudent])
    target_present = Set([:predicted, :rstudent])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = [:stdp]
    target_need = Set([:predicted, :leverage, :stdp])
    target_present = Set([:predicted, :stdp])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = [:stdp, :stdi]
    target_need = Set([:predicted, :leverage, :stdp, :stdi])
    target_present = Set([:predicted, :stdp, :stdi])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = [:stdpp, :stdr]
    target_need = Set([:predicted, :leverage, :stdr])
    target_present = Set([:predicted, :stdr])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = []
    target_need = Set([:predicted])
    target_present = Set([:predicted])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = Set()
    target_need = Set([:predicted])
    target_present = Set([:predicted])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

end


text:
getpvtestresult(learner, 300) |> svg
getpolicy(learner, getstate(Tac([2,8,6]))) |> display
getvalue(learner, getstate(Tac([2,8,6])), 1) |> display
getvalue(learner, getstate(Tac([2,8,6])), 3) |> display
getvalue(learner, getstate(Tac([2,8,6,9]))) |> display
getvalue(learner, getstate(Tac([2,8,6])), 9) |> display
getstate(Tac([2,8,6]))
mzg = play_game(learner, [2,8,6])
ts1 = TrainingSample(mzg, 1, 9, 9)

# STOPPED HERE
# The following two are not similar! (and
# this is the main problem). UCB score is being badly
# effected by the bad recurrent_inference value results
# via the inferred value in backpropagate.
getvalue(learner, getstate(Tac([2,8,6,9,7,1,5, 4])))
getvalue(learner, getstate(Tac([2,8,6,9,7,1,5]), 4))

getstate(Tac([2,8,6])) |> display

value_fn = getvaluefn(Tac)
optimalp = getoptimalpolicyfn(Tac)

###########
s1 = getstate(Tac([2,8,6,1]))
value_fn[s1]
p1 = getpolicy(learner, s1)
optimalp(s1)

s3 = getstate(Tac([2,8,6,3]))
value_fn[s3]
p3 = getpolicy(learner, s3)
optimalp(s3)
###########

s4 = getstate(Tac([2,8,6,4]))
s4 |> display
value_fn[s4]
p4 = getpolicy(learner, s4)
optimalp(s4)

s5 = getstate(Tac([2,8,6,5]))
s5 |> display
value_fn[s5]
p5 = getpolicy(learner, s5)
optimalp(s5)

s7 = getstate(Tac([2,8,6,7]))
s7 |> display
value_fn[s7]
p7 = getpolicy(learner, s7)
optimalp(s7)

s9 = getstate(Tac([2,8,6,9]))
s9 |> display
value_fn[s9]
p9 = getpolicy(learner, s9)
optimalp(s9)

import Pkg
Pkg.rm(["Grids", "GridGames"])
Pkg.develop(["Grids", "GridGames"])

using MuZero
using GridGames
using Base.Threads: @spawn
cd(@__DIR__)
include("game.jl")
include("muzero_config.jl")
config = getconfig(Tac)
include("inference_tests.jl")
pv_tests = makepolicyvaluetests(Tac)
pv_tests.period = 100
learner = Learner(
    config;
    uuid="435c256d-1c60-4e07-ac03-a2430ae7b7bd",
    dir="~/.julia/dev/MuZeroModels/tictactoe/",
    tests=[pv_tests],
    schedule=tacschedule!
)
getpvtestresult(learner, 300)
getpvtestresult(learner, 1405) |> svg
length(learner.replay_buffers)
mzgame = learner.replay_buffers[1]
mzgame.game.states[end]
target = getstate(Tac([1,4,8,2])) #6,3,6
found = findmuzerogames(learner, target)

length(learner.replay_buffers)
learner.stats.policy_test_fails
#40 fails: 6,8,25,31,47,51,57,78,116,138,139,142,143,180,185,186,187,201,202,218,225,259,284,300,323,334,418,512,514,538,628,725,778,1098,1114,1188,1413,2120,2760,300
showresultsfull(learner; mean=100)
get_steps(learner)
acc, fails, value_mse = runtests(pv_tests, training_network(learner))

using Flux
test_net = new_test_network(learner.config)
testmode!.(get_nets(test_net), true)
updatetestparams!(learner.train, test_net)
testmode!.(get_nets(test_net), true)
acc2, fails2, value_mse2 = runtests(pv_tests, test_net)

#run all tests one by one! (instead of in batch mode)
fails = []
for pvtest in pv_tests
    out = initial_inference(test_net, pvtest.state) |> NetworkOutput
    append!(fails, getpolicyfails(pvtest, out))
end

fit!(learner; steps=500)
get_steps(learner)
learner.stats.policy_test_fails
# 45 fails: 6,25,31,47,51,57,78,116,121,128,131,138,139,142,172,180,185,187,201,218,225,259,284,323,334,418,512,514,538,628,725,761,778,839,1098,1114,1188,1318,1413,1606,2066,2120,2760,3254,300

getpvtestresult(learner, 3254) |> svg
target2 = getstate(Tac([1,5,4,6,2,9]))
found = findmuzerogames(
    learner,
    getpolicyvaluetests(learner)[3254].state
)

fit!(learner; steps=90)
# 47 fails: 6,25,31,47,51,57,78,116,121,128,131,139,142,167,172,180,185,187,201,218,225,238,259,284,285,323,334,418,512,514,538,628,778,834,839,1098,1114,1188,1318,1413,1606,1861,2120,2760,3254,3286,300
found = findmuzerogames(
    learner,
    getpolicyvaluetests(learner)[3286].state
)
getpolicy(
    learner,
    getpolicyvaluetests(learner)[1405].state
) |> display
getpvtestresult(learner, 3254) |> svg


text:
"""
**Gaussian Likelihood**

Classical Gaussian noise : ``p(y|f) = \\mathcal{N}(y|f,\\epsilon)``

```julia
GaussianLikelihood(ϵ::T=1e-3) #ϵ is the variance
```

There is no augmentation needed for this likelihood which is already conjugate
"""
struct GaussianLikelihood{T<:Real} <: RegressionLikelihood{T}
    ϵ::LatentArray{T}
    θ::LatentArray{Vector{T}}
    function GaussianLikelihood{T}(ϵ::AbstractVector{T}) where {T<:Real}
        new{T}(ϵ)
    end
    function GaussianLikelihood{T}(ϵ::AbstractVector{T},θ::AbstractVector{<:AbstractVector{T}}) where {T<:Real}
        new{T}(ϵ,θ)
    end
end

function GaussianLikelihood(ϵ::T=1e-3) where {T<:Real}
    GaussianLikelihood{T}([ϵ])
end

function GaussianLikelihood(ϵ::AbstractVector{T}) where {T<:Real}
    GaussianLikelihood{T}(ϵ)
end

function pdf(l::GaussianLikelihood,y::Real,f::Real)
    pdf(Normal(y,l.ϵ[1]),f) #WARNING multioutput invalid
end

function logpdf(l::GaussianLikelihood,y::Real,f::Real)
    logpdf(Normal(y,l.ϵ[1]),f) #WARNING multioutput invalid
end

function Base.show(io::IO,model::GaussianLikelihood{T}) where T
    print(io,"Gaussian likelihood")
end

function compute_proba(l::GaussianLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    return μ,max.(σ²,0.0).+ l.ϵ[1]
end

function init_likelihood(likelihood::GaussianLikelihood{T},inference::Inference{T},nLatent::Int,nSamplesUsed::Int,nFeatures::Int) where {T<:Real}
    if length(likelihood.ϵ) ==1 && length(likelihood.ϵ) != nLatent
        return GaussianLikelihood{T}([likelihood.ϵ[1] for _ in 1:nLatent],[fill(inv(likelihood.ϵ[1]),nSamplesUsed) for _ in 1:nLatent])
    elseif length(likelihood.ϵ) != nLatent
        @warn "Wrong dimension of ϵ : $(length(likelihood.ϵ)), using first value only"
        return GaussianLikelihood{T}([likelihood.ϵ[1] for _ in 1:nLatent])
    else
        return GaussianLikelihood{T}(likelihood.ϵ,[fill(likelihood.ϵ[i],nSamplesUsed) for i in 1:nLatent])
    end
end

function local_updates!(model::GP{T,<:GaussianLikelihood}) where {T}
    if model.opt_noise
        model.likelihood.ϵ .= inv(model.nFeatures)*broadcast((y,μ,Σ)->sum(abs2,y)-2*dot(y,μ)+sum(abs2,μ)+tr(Σ),model.y,model.Knn.*model.invKnn.*model.y,model.Knn.*([I].-model.invKnn.*model.Knn))
    end
end

function local_updates!(model::SVGP{T,<:GaussianLikelihood}) where {T}
    if model.inference.Stochastic
        #TODO make it a moving average
        ρ = inv(sqrt(1+model.inference.nIter))
        model.likelihood.ϵ .= (1-ρ)*model.likelihood.ϵ + ρ/model.inference.nSamplesUsed *broadcast((y,κ,μ,Σ,K̃)->sum(abs2.(y-κ*μ))+opt_trace(κ*Σ,κ)+sum(K̃),model.inference.y,model.κ,model.μ,model.Σ,model.K̃)
        model.likelihood.ϵ
    else
        model.likelihood.ϵ .= 1.0/model.inference.nSamplesUsed *broadcast((y,κ,μ,Σ,K̃)->sum(abs2.(y-κ*μ))+opt_trace(κ*Σ,κ)+sum(K̃),model.y,model.κ,model.μ,model.Σ,model.K̃)
    end
    model.likelihood.θ .= broadcast(ϵ->fill(inv(ϵ),model.inference.nSamplesUsed),model.likelihood.ϵ)
end

function local_updates!(model::VStP{T,<:GaussianLikelihood}) where {T}
    model.likelihood.ϵ .= 1.0/model.inference.nSamples *broadcast((y,μ,Σ)->sum(abs2,y-μ)+Σ,model.y,model.μ,tr.(model.Σ))
    model.likelihood.θ .= broadcast(ϵ->fill(inv(ϵ),model.inference.nSamplesUsed),model.likelihood.ϵ)
end

@inline ∇E_μ(model::AbstractGP{T,<:GaussianLikelihood,<:AnalyticVI}) where {T} = model.inference.y./model.likelihood.ϵ
@inline ∇E_μ(model::AbstractGP{T,<:GaussianLikelihood,<:AnalyticVI},i::Int) where {T} = model.inference.y[i]./model.likelihood.ϵ[i]
@inline ∇E_Σ(model::AbstractGP{T,<:GaussianLikelihood,<:AnalyticVI}) where {T} = 0.5*model.likelihood.θ
@inline ∇E_Σ(model::AbstractGP{T,<:GaussianLikelihood,<:AnalyticVI},i::Int) where {T} = 0.5*model.likelihood.θ[i]

function predict_f(model::GP{T,GaussianLikelihood{T},Analytic{T}},X_test::AbstractMatrix{T};covf::Bool=true,fullcov::Bool=false) where {T}
    k_star = kernelmatrix.([X_test],[model.X],model.kernel)
    μf = k_star.*model.invKnn.*model.y
    if !covf
        return model.nLatent == 1 ? μf[1] : μf
    end
    if fullcov
        Σf = Symmetric.(kernelmatrix.([X_test],model.kernel) .- k_star.*model.invKnn.*transpose.(k_star))
        i = 0
        ϵ = 1e-16
        while count(isposdef.(Σf))!=model.nLatent
            Σf .= ifelse.(isposdef.(Σf),Σf,Σf.+ϵ.*[I])
            if i > 100
                println("DAMN")
                break;
            end
            ϵ *= 2
            i += 1
        end
        @assert count(isposdef.(Σf))==model.nLatent
        return model.nLatent == 1 ? (μf[1],Σf[1]) : (μf,Σf)
    else
        σ²f = kerneldiagmatrix.([X_test],model.kernel) .- opt_diag.(k_star.*model.invKnn,k_star)
        return model.nLatent == 1 ? (μf[1],σ²f[1]) : (μf,σ²f)
    end
end

function proba_y(model::GP{T,GaussianLikelihood{T},Analytic{T}},X_test::AbstractMatrix{T}) where {T}
    μf, σ²f = predict_f(model,X_test,covf=true)
    σ²f .+= model.likelihood.ϵ
    return μf,σ²f
end

function proba_y(model::SVGP{T,GaussianLikelihood{T},AnalyticVI{T}},X_test::AbstractMatrix{T}) where {T}
    μf, σ²f = predict_f(model,X_test,covf=true)
    σ²f .+= model.likelihood.ϵ
    return μf,σ²f
end

### Special case where the ELBO is equal to the marginal likelihood
function ELBO(model::GP{T,GaussianLikelihood{T}}) where {T}
    return -0.5*sum(broadcast((y,invK)->dot(y,invK*y) - logdet(invK)+ model.nFeatures*log(twoπ),model.y,model.invKnn))
end

function ELBO(model::SVGP{T,GaussianLikelihood{T}}) where {T}
    return expecLogLikelihood(model) - GaussianKL(model)
end

function expecLogLikelihood(model::SVGP{T,GaussianLikelihood{T}}) where {T}
    return -0.5*model.inference.ρ*sum(broadcast((y,ϵ,κ,Σ,μ,K̃)->1.0/ϵ*(sum(abs2.(y[model.inference.MBIndices]-κ*μ))+sum(K̃)+opt_trace(κ*Σ,κ))+model.inference.nSamplesUsed*(log(twoπ)+log(ϵ)),model.y,model.likelihood.ϵ,model.κ,model.Σ,model.μ,model.K̃))
end

function hyperparameter_gradient_function(model::GP{T,GaussianLikelihood{T}}) where {T}
    A = ([I].-model.invKnn.*(model.y.*transpose.(model.y))).*model.invKnn
    if model.IndependentPriors
        return (function(Jnn,index)
                    return -0.5*hyperparameter_KL_gradient(Jnn,A[index])
                end,
                function(kernel,index)
                    return -0.5/getvariance(kernel)*opt_trace(model.Knn[index],A[index])
                end,
                function(index)
                    return -model.invKnn[index]*(model.μ₀[index]-model.y[index])
                end)
    else
        return (function(Jnn,index)
            return -0.5*sum(hyperparameter_KL_gradient(Jnn,A[i]) for i in 1:model.nLatent)
                end,
                function(kernel,index)
                    return -0.5/getvariance(kernel)*sum(opt_trace(model.Knn[1],A[i]) for i in 1:model.nLatent)
                end,
                function(index)
                    return -sum(model.invKnn.*(model.μ₀.-model.μ))
                end)
    end
end


text:
"""
    struct ComplexSpace <: EuclideanSpace{ℂ}

A standard complex vector space ℂ^d with Euclidean inner product and no additional
structure. It is completely characterised by its dimension and whether its the normal space
or its dual (which is canonically isomorphic to the conjugate space).
"""
struct ComplexSpace <: EuclideanSpace{ℂ}
  d::Int
  dual::Bool
end
ComplexSpace(d::Int) = ComplexSpace(d, false)

# convenience constructor
Base.:^(::ComplexNumbers, d::Int) = ComplexSpace(d)
Base.getindex(::ComplexNumbers) = ComplexSpace
Base.getindex(::ComplexNumbers, d::Int) = ComplexSpace(d)

# Corresponding methods:
#------------------------
dim(V::ComplexSpace) = V.d
Base.axes(V::ComplexSpace) = Base.OneTo(dim(V))

Base.conj(V::ComplexSpace) = ComplexSpace(V.d, !V.dual)
isdual(V::ComplexSpace) = V.dual

Base.oneunit(::Type{ComplexSpace}) = ComplexSpace(1)
⊕(V1::ComplexSpace, V2::ComplexSpace) = (V1.dual==V2.dual ?
    ComplexSpace(V1.d+V2.d, V1.dual) :
    throw(SpaceMismatch("Direct sum of a vector space and its dual does not exist")))
fuse(V1::ComplexSpace, V2::ComplexSpace) = ComplexSpace(V1.d*V2.d)
flip(V::ComplexSpace) = dual(V)

Base.min(V1::ComplexSpace, V2::ComplexSpace) = V1.dual == V2.dual ?
    ComplexSpace(min(V1.d, V2.d), V1.dual) :
    throw(SpaceMismatch("V1 and V2 should both be normal or dual spaces"))
Base.max(V1::ComplexSpace, V2::ComplexSpace) = V1.dual == V2.dual ?
    ComplexSpace(max(V1.d, V2.d), V1.dual) :
    throw(SpaceMismatch("V1 and V2 should both be normal or dual spaces"))

Base.show(io::IO, V::ComplexSpace) = print(io, V.dual ? "(ℂ^$(V.d))'" : "ℂ^$(V.d)")
Base.show(io::IO, ::Type{ComplexSpace}) = print(io, "ComplexSpace")


text:
using tele
using Test

@testset "tele.jl" begin
    # Write your tests here.
end


text:
"""
    logo(cmd0::String=""; kwargs...)

Plots the GMT logo on a map. By default, the GMT logo is 5 cm wide and 2.5 cm high and will be
positioned relative to the current plot origin. Use various options to change this and to place
a transparent or opaque rectangular map panel behind the GMT logo.

Full option list at [`gmtlogo`](http://gmt.soest.hawaii.edu/doc/latest/gmtlogo.html)

Parameters
----------

- **D** : **pos** : **position** : -- Str --

	Sets reference point on the map for the image using one of four coordinate systems.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/gmtlogo.html#d)
- **F** : **box** : -- Str --

	Without further options, draws a rectangular border around the GMT logo using MAP_FRAME_PEN.
    or map rose (T)
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/gmtlogo.html#f)
- **julia** : -- Number --

	Create the Julia instead of the GMT logo. Provide circle diameter in centimeters
- **GMTjulia** : -- Number --

    Create the GMT Julia GMT logo. Provide circle diameter in centimeters
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- $(GMT.opt_P)
- $(GMT.opt_R)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_t)

- Example, make a GMT Julia logo with circles of 1 cm: logo(GMTjulia=1, show=true)
"""
function logo(cmd0::String=""; K=false, O=false, first=true, kwargs...)

	length(kwargs) == 0 && return monolitic("gmtlogo", cmd0, arg1)	# Speedy mode
	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

    cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, "")
	cmd, = parse_JZ(cmd, d)
	cmd = parse_UVXY(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd = parse_params(cmd, d)

	cmd, K, O = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	cmd = add_opt(cmd, 'D', d, [:D :pos :position])
	cmd = add_opt(cmd, 'F', d, [:F :box])

	cmd = finish_PS(d, cmd, output, K, O)

	do_julia    = haskey(d, :julia)
	do_GMTjulia = haskey(d, :GMTjulia)
	if (do_julia || do_GMTjulia)
		if (do_julia)	r = d[:julia]	else	r = d[:GMTjulia]	end
		c,t,r2 = jlogo(r)			# r2 is the diameter of the inner circle
		if (!occursin("-R", cmd))  cmd = @sprintf("-R0/%f/0/%f ", 2r, 2r) * cmd  end
		if (!occursin("-J", cmd))  cmd = " -Jx1 " * cmd  end
		do_show = false
		if (do_GMTjulia && haskey(d, :show))	# Too soon to make the display
			delete!(d, :show);	do_show = true
		end
		fmt = fname_ext
		if (do_GMTjulia && haskey(d, :fmt))		# Too soon to set the format. Need to finish the PS first
			fmt = d[:fmt];	delete!(d, :fmt);
			fname_ext = "ps"
		end
		finish_PS_module(d, c * cmd, "", output, fname_ext, opt_T, K, "psxy", t)
		if (haskey(d, :GMTjulia))
			letter_height = 0.75 * r2 / 2.54 * 72 		# Make the letters 75% of the cicle's diameter
			opt_F = @sprintf("+f%d,NewCenturySchlbk-Italic",letter_height)
			text!(text_record(t[1:3,1:2], ["G", "T", "M"]), R=[], J=[], F=opt_F, fmt=fmt, show=do_show)
		end
	else
		if (!occursin("-D", cmd))  cmd = " -Dx0/0+w5c " * cmd	end
		return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "gmtlogo")
	end
end

# ---------------------------------------------------------------------------------------------------
logo!(cmd0::String=""; K=true, O=true, first=false, kw...) = logo(cmd0; K=K, O=O, first=first, kw...)
logo(; K=false, O=false, first=true, kw...) = logo(""; K=K, O=O, first=first, kw...)
logo!(; K=true, O=true, first=false, kw...) = logo(""; K=K, O=O, first=first, kw...)

# -------------------------------------------------------------------------
function jlogo(L=5)
	# Create the Julia "Terminator" 3 colored circles triangle
	# L is the length of the equilateral triangl
	W = 2 * L 					# Region width
	H = L * sind(60) 			# Triangle height
	s_size = 0.8 * L 			# Circle diameter
	l_thick = s_size * 0.06 	# Line thickness

	s1 = s_size					# Outer circle diameter. To simulate a line.
	if (GMTver < 6)  s1 = s1 / 2.54  end	# There is a bug in older versions
	s2 = s1 * (1 - 0.06)		# Inner circle diameter. The one that will be filled.
	t = [L/2 L/2 0 s1; L+L/2 L/2 1 s1; L L/2+H 2 s1; L/2 L/2 3 s2; L+L/2 L/2 4 s2; L L/2+H 5 s2]
	cmd = " -Sc -C171/43/33,130/83/171,81/143/24,191/101/95,158/122/190,128/171/93 "
	return cmd, t, s2
end


text:
using DelimitedFiles

"""
ULMFiT - LANGUAGE MODEL

    LanguageModel(load_pretrained::Bool=false, vocabpath::String="vocabs/lm_vocab.csv";kw...)

The Language model structure for ULMFit is defined by 'LanguageModel' struct.
It contains has two fields:
    vocab   : vocabulary, which will be used for language modelling
    layers  : embedding, RNN and dropout layers of the whole model
In this language model, the embedding matrix used in the embedding layer
is same for the softmax layer, following Weight-tying technique.
The field 'layers' also includes the Variational Dropout layers.
It takes several dropout probabilities for different dropout for different layers.

[Usage and arguments are discussed in the docs]
There are several keyword argunments to set the dropout probabilities
of the layers of model checkout those arguments in the docs.
# Example:

julia> lm = LanguageModel()
"""
mutable struct LanguageModel
    vocab :: Vector
    layers :: Flux.Chain
end

function LanguageModel(load_pretrained::Bool=false, vocabpath::String=joinpath(@__DIR__,"vocabs/lm_vocab.csv");embedding_size::Integer=400, hid_lstm_sz::Integer=1150, out_lstm_sz::Integer=embedding_size,
    embed_drop_prob::Float64 = 0.05, in_drop_prob::Float64 = 0.4, hid_drop_prob::Float64 = 0.5, layer_drop_prob::Float64 = 0.3, final_drop_prob::Float64 = 0.3)
    vocab = (string.(readdlm(vocabpath, ',')))[:, 1]
    de = gpu(DroppedEmbeddings(length(vocab), embedding_size, embed_drop_prob; init = (dims...) -> init_weights(0.1, dims...)))
    lm = LanguageModel(
        vocab,
        Chain(
            de,
            VarDrop(in_drop_prob),
            gpu(AWD_LSTM(embedding_size, hid_lstm_sz, hid_drop_prob; init = (dims...) -> init_weights(1/hid_lstm_sz, dims...))),
            VarDrop(layer_drop_prob),
            gpu(AWD_LSTM(hid_lstm_sz, hid_lstm_sz, hid_drop_prob; init = (dims...) -> init_weights(1/hid_lstm_sz, dims...))),
            VarDrop(layer_drop_prob),
            gpu(AWD_LSTM(hid_lstm_sz, out_lstm_sz, hid_drop_prob; init = (dims...) -> init_weights(1/hid_lstm_sz, dims...))),
            VarDrop(final_drop_prob),
            x -> de(x, true),
            softmax
        )
    )
    load_pretrained && load_model!(lm, datadep"Pretrained ULMFiT Language Model/ulmfit_lm_en.bson")
    return lm
end

Flux.@functor LanguageModel

"""
    test_lm(lm::LanguageModel, data_gen, num_of_iters::Integer; unknown_token::String="_unk_")

This function is used to test the `LanguageModel` on the given data given by data_gen.
`num_of_iters` refers to number of batches for which the model has to be tested.
It returns loss, accuracy, precsion, recall and F1 score.

# Example:

julia> test_lm(lm, data_gen, 200, "<unk")
"""
function test_lm(lm::LanguageModel, data_gen, num_of_iters::Integer; unknown_token::String="_unk_")
    model_layers = lm.layers
    testmode!(model_layers)
    loss = 0
    len = length(vocab)
    TP, FP, FN, TN = zeros(len, 1), zeros(len, 1), zeros(len, 1), zeros(len, 1)
    for iter=1:num_of_iters
        X, Y = take!(gen), take!(gen)
        H = broadcast(w -> indices(w, lm.vocab, unknown_token), X)
        H = model_layers.(H)
        Y = broadcast(x -> gpu(Flux.onehotbatch(x, lm.vocab, unknown_token)), Y)
        l = sum(crossentropy.(H, Y))
        Flux.reset!(model_layers)
        for (h, y) in zip(H, Y)
            TP .+= sum(h .* y, dims=2)
            FN .+= sum(((-1 .* h) .+ 1) .* y, dims=2)
            FP .+= sum(h .* ((-1 .* y) .+ 1), dims=2)
            TN .+= sum(((-1 .* h) .+ 1) .* ((-1 .* y) .+ 1), dims=2)
        end
        loss += l/length(X[1])
    end
    precision = *((TP./(TP .+ FP))...)
    recall = *((TP./(TP .+ FN))...)
    F1 = (2*precision*recall)/(precision+recall)
    accuracy = (TP[1] + TN[1])/(TP[1] + TN[1] + FP[1] + FN[1])
    return (loss/num_of_iters), accuarcy, precision, recall, F1
end

# computes the forward pass while training
function forward(lm, batch)
    batch = map(x -> indices(x, lm.vocab, "_unk_"), batch)
    batch = gpu(batch)
    batch = lm.layers.(batch)
    return batch
end

# loss funciton - Calculates crossentropy loss
function loss(lm, gen)
    H = forward(lm, take!(gen))
    Y = broadcast(x -> gpu(Flux.onehotbatch(x, lm.vocab, "_unk_")), take!(gen))
    l = sum(Flux.crossentropy.(H, Y))
    reset!(lm.layers)
    return l
end

# Backpropagation step while training
function backward!(layers, lm, gen, opt)
    # Calulating gradients and weights updation
    p = get_trainable_params(layers)
    grads = Zygote.gradient(() -> loss(lm, gen), p)
    Flux.Optimise.update!(opt, p, grads)
    return
end

"""
pretrain_lm!

This funciton contains main training loops for pretrainin the Language model
including averaging step for the 'AWD_LSTM' layers.

Usage and arguments are explained in the docs of ULMFiT
"""
function pretrain_lm!(lm::LanguageModel=LanguageModel(), data_loader::Channel=load_wikitext_103;
    base_lr=0.004, epochs::Integer=1, checkpoint_iter::Integer=5000)

    # Initializations
    opt = ADAM(base_lr, (0.7, 0.99))    # ADAM Optimizer

    # Pre-Training loops
    for epoch=1:epochs
        println("\nEpoch: $epoch")
        gen = data_loader()
        num_of_batches = take!(gen) # Number of mini-batches
        T = num_of_iters-Int(floor((num_of_iters*2)/100))   # Averaging Trigger
        set_trigger!.(T, lm.layers)  # Setting triggers for AWD_LSTM layers
        for i=1:num_of_batches

            # REVERSE PASS
            backward!(lm.layers, lm, gen, opt)

            # ASGD Step, works after Triggering
            asgd_step!.(i, lm.layers)

            # Resets dropout masks for all the layers with Varitional DropOut or DropConnect masks
            reset_masks!.(lm.layers)

            # Saving checkpoints
            if i == checkpoint_iter save_model!(lm) end
        end
    end
end

# To save model
function save_model!(m::LanguageModel, filepath::String)
    weights = cpu.(params(m))
    BSON.@save filepath weights
end

# To load model
function load_model!(lm::LanguageModel, filepath::String)
    BSON.@load filepath weights
    # reshape saved weights to match Recurr (h, c) shape
    layers = [5, 6, 10, 11, 15, 16]
    for l in layers
        weights[l] = reshape(weights[l], length(weights[l]), 1)
    end
    Flux.loadparams!(lm, weights)
end

"""
sample(starting_text::AbstractDocument, lm::LanguageModel)

Prints sampling results taking `starting_text` as initial tokens for the sampling for LanguageModel.

# Example:

julia> sampling("computer science", lm)
SAMPLING...
, is fast growing field ......

"""
function sample(starting_text::AbstractDocument, lm::LanguageModel)
    testmode!(lm.layers)
    model_layers = lm.layers
    tokens = tokens(starting_text)
    word_indices = map(x -> indices([x], lm.vocab, "_unk_"), tokens)
    h = (model_layers.(word_indices))[end]
    prediction = lm.vocab[argmax(h)[1]]
    println("SAMPLING...")
    print(prediction, ' ')
    while true
        h = indices([prediction], lm.vocab, "_unk_")
        h = model_layers(h)
        prediction = lm.vocab[argmax(h)[1]]
        print(prediction, ' ')
        prediction == "_pad_" && break
    end
end


text:
## ONE FACTOR MODELS ##
type OneFactorShortRateTree{S <: ShortRateDynamics} <: ShortRateTree
  tree::TrinomialTree
  dynamics::S
  tg::TimeGrid
  treeLattice::TreeLattice1D

  function OneFactorShortRateTree{S}(tree::TrinomialTree, dynamics::S, tg::TimeGrid)
    oneFactorTree = new(tree, dynamics, tg)
    oneFactorTree.treeLattice = TreeLattice1D(tg, get_size(tree, 2), oneFactorTree)

    return oneFactorTree
  end
end

get_size(tr::OneFactorShortRateTree, i::Int) = get_size(tr.tree, i)

function discount(tr::OneFactorShortRateTree, i::Int, idx::Int)
  x = get_underlying(tr.tree, i, idx)
  r = short_rate(tr.dynamics, tr.tg.times[i], x)
  return exp(-r * tr.tg.dt[i])
end

descendant(tr::OneFactorShortRateTree, i::Int, idx::Int, branch::Int) = descendant(tr.tree, i, idx, branch)
probability(tr::OneFactorShortRateTree, i::Int, idx::Int, branch::Int) = probability(tr.tree, i, idx, branch)

get_params(m::OneFactorModel) = Float64[get_a(m), get_sigma(m)]

type RStarFinder{M <: ShortRateModel}
  model::M
  strike::Float64
  maturity::Float64
  valueTime::Float64
  fixedPayTimes::Vector{Float64}
  amounts::Vector{Float64}
end

function operator(rsf::RStarFinder)
  function _inner(x::Float64)
    _value = rsf.strike
    _B = discount_bond(rsf.model, rsf.maturity, rsf.valueTime, x)
    sz = length(rsf.fixedPayTimes)
    for i = 1:sz
      dbVal = discount_bond(rsf.model, rsf.maturity, rsf.fixedPayTimes[i], x) / _B
      _value -= rsf.amounts[i] * dbVal
    end

    return _value
  end

  return _inner
end

discount_bond(model::OneFactorModel, tNow::Float64, maturity::Float64, factors::Vector{Float64}) = discount_bond(model, tNow, maturity, factors[1])
discount_bond{M <: OneFactorModel}(model::M, tNow::Float64, maturity::Float64, _rate::Float64) = A(model, tNow, maturity) * exp(-B(model, tNow, maturity) * _rate)


text:
abstract type InterpretableModel end

notdef(T::DataType, method) =
    error("Method $(string(method)) not defined for type $(T)")
notdef(o, method) = notdef(typeof(o), method)
istimeseriesmodel(m::InterpretableModel) = false
isclassificationmodel(m::InterpretableModel) = notdef(m, :isclassificationmodel)
isregressionmodel(m::InterpretableModel) = notdef(m, :isregressionmodel)

abstract type InterpretableRegressionModel <: InterpretableModel end
isclassificationmodel(m::InterpretableRegressionModel) = false
isregressionmodel(m::InterpretableRegressionModel) = true

abstract type InterpretableClassificationModel <: InterpretableModel end
isclassificationmodel(m::InterpretableClassificationModel) = true
isregressionmodel(m::InterpretableClassificationModel) = false

abstract type InterpretableTimeseriesModel <: InterpretableModel end
istimeseriesmodel(m::InterpretableTimeseriesModel) = true

"""
    horizon(model)

The number of steps into the future that a time series model can predict.
Returns Inf if it can predict a flexible number of steps into the future.
"""
horizon(m::InterpretableTimeseriesModel) = notdef(m, :horizon)

abstract type InterpretableTimeseriesRegressionModel <: InterpretableTimeseriesModel end
isclassificationmodel(m::InterpretableTimeseriesRegressionModel) = false
isregressionmodel(m::InterpretableTimeseriesRegressionModel) = true


text:
module Loading

using TOML
using Pkg3.Operations: user_depot, load_manifest

abstract type Loader end

struct LoadInstalled <: Loader
    depot::String
end

function Base.load_hook(loader::LoadInstalled, name::String, found)
    # load and parse the package manifest
    manifest = load_manifest()
    # pull the package UUID and verion hash out of the manifest
    package  = haskey(manifest, name)       ? manifest[name]       : return found
    uuid     = haskey(package, "uuid")      ? package["uuid"]      : return found
    hash     = haskey(package, "hash-sha1") ? package["hash-sha1"] : return found
    # relative path for package and version
    path = joinpath(uuid, hash, "src", "$name.jl")
    # see if previously found matches expected suffix
    found isa String && endswith(found, path) && return found
    # see if package and version is installed in the current depot
    path = joinpath(loader.depot, "packages", path)
    # if it is, return it, otherwise pass through
    return isfile(path) ? path : found
end

empty!(LOAD_PATH)
ENV["JULIA_PKGDIR"] = tempname()
@eval Base module Loading; DEPOTS = []; end
push!(LOAD_PATH, LoadInstalled(user_depot()))
push!(Base.Loading.DEPOTS, user_depot())

end # module


text:
# primary features:
# rollout sequence of controls
# df/dx, df/du derivatives

# features to add:
# address mjdata, mjmodel fields by string or symbol
#   replicate the get_sensors stuff
__precompile__()

module mjWrap

using MuJoCo
using JLD2

using UnsafeArrays
using LinearAlgebra.BLAS
using Distributed

struct mjSet
    m::mj.jlModel
    d::mj.jlData
    datas::Vector{mj.jlData}

    # helpers
    sensor_names::Dict{Symbol, AbstractRange}
    name::String
    skip::Integer
    mode::Integer

    nq::Integer
    nv::Integer
    nu::Integer
    ns::Integer
    dt::Float64

    #function mjSet(model_name::String) # TODO?
    #end
end

struct TrajSamples
    T::Integer
    K::Integer
    s0::Vector{Float64}      # nqnv
    ctrl::Array{Float64,3}   # nu x T x K
    state::Array{Float64,3}  # nqnv x T x K
    obs::Array{Float64,3}    # nu x T x K
    reward::Array{Float64,2} # T x K
end

struct Trajectory
    ctrl::Matrix{Float64}   # nu x T x K
    state::Matrix{Float64}  # nqnv x T x K
    obs::Matrix{Float64}    # nu x T x K
    reward::Vector{Float64} # T x K
end

#=
type TrajDerivatives
    T::Integer
    K::Integer
    d_x::Array{Float64}
    d_u::Array{Float64}
    s_x::Array{Float64}
    s_u::Array{Float64}

    #c::Array{Float64}
    #c_x::Array{Float64}
    #c_u::Array{Float64}
    #c_xx::Array{Float64}
    #c_uu::Array{Float64}
    #c_ux::Array{Float64}
end
=#

const mjw = mjWrap
export TrajSamples, mjSet, mjw
const GDD_POS = 1
const GDD_ACC = 2
const GDD_Q = 3
const GDD_COST = 4

function finish()
    mj.deactivate()
    return 0;
end

function modelparams(x::mjSet)
    return x.nq, x.nv, x.nu, x.ns
end

function clear_model(x::mjSet)
    mj.deleteModel(x.m.m)
    mj.deleteData(x.d.d)
    for i=1:length(x.datas)
        mj.deleteData(x.datas[i])
    end
end

function load_model(model_name::String, skip::Integer, mode::String="normal", ns::Integer=0)
    pm = mj_loadXML(model_name, C_NULL)
    if pm == nothing
        return
    end
    pd = mj_makeData(pm)

    m, d = mj.mapmujoco(pm, pd)
    ndata = Threads.nthreads()
    datas = Vector{jlData}(undef, ndata)
    print("Making mjDatas: 1 ")
    datas[1] = d
    for i=2:ndata
        print("$i ")
        datas[i] = mj.mapdata(pm, mj_makeData(pm))
    end
    println()

    if m.m[].nsensor > 0
        sensors = mj.name2range(m, m.m[].nsensor,
                                m.name_sensoradr, m.sensor_adr, m.sensor_dim)
    else
        sensors = Dict(:none=>0:0)
    end

    nq = Int64(m.m[].nq)
    nv = Int64(m.m[].nv)
    nu = Int64(m.m[].nu)
    if ns==0 ns = Int64(m.m[].nsensordata) end
    dt = m.m[].opt.timestep #mj.get(m, :opt, :timestep)

    return mjSet(m, d, datas, sensors, model_name, skip, 0,
                 nq, nv, nu, ns, dt)
end

############ UTILS ############
function wraptopi(ang::Float64)
    return mod(ang+pi, 2pi) - pi
end

function unitfit(ang::Float64, val::Float64)
    clamp.(ang, -val, val)./val
end

function inertiaweight(m::mj.jlModel, d::mj.jlData,
                       v::AbstractVector{Float64})
    nv = length(d.qvel)
    mark = mj.MARKSTACK(d)

    My = @view d.stack[(mark+1):(mark+nv)]
    mj_solveM(m.m, d.d, My, v, 1)
    force = v'*My

    mj.FREESTACK(d, Int32(mark))
    return force/mj_getTotalmass(m.m)
end

#function stackAlloc(m::mj.jlModel, d::mj.jlData)
#   mark = mj.MARKSTACK(d)
#end

function limitcontrols(x::mjSet, s::TrajSamples)
    for u=1:x.nu
        clamp!(uview(s.ctrl,u,:,:),
               x.m.actuator_ctrlrange[1,u],
               x.m.actuator_ctrlrange[2,u])
    end
end

# modify ctrl vector with the noise stored in ctrl
function inertiascale!(m::mj.jlModel, d::mj.jlData,
                       input::AbstractVector{Float64},
                       output::AbstractVector{Float64})
    mj.solveM2(m.m, d.d, output, input, 1)
end

function state2data(d::mj.jlData, s::AbstractVector{Float64})
    nq = length(d.qpos)
    nv = length(d.qvel)
    copyto!(d.qpos, 1, s, 1, nq)
    copyto!(d.qvel, 1, s, nq+1, nv)
end

function data2state(s::AbstractVector{Float64}, d::mj.jlData)
    nq = length(d.qpos)
    nv = length(d.qvel)
    copyto!(s, d.qpos)
    copyto!(s, nq+1, d.qvel, 1, nv)
end

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

function reset(m::mj.jlModel, d::mj.jlData,
               s0::AbstractVector{Float64},
               ctrl0::AbstractVector{Float64},
               obs::AbstractVector{Float64})
    # reset
    d.qacc .= 0.0
    state2data(d, s0)
    copyto!(d.ctrl, ctrl0)
    d.d[].time = 0.0 #mj.set(d, :time, 0.0)

    # out
    if length(d.sensordata) > 0
        mj_forward(m, d)
        copyto!(obs, d.sensordata)
    else
        mj_forwardSkip(m.m, d.d, Integer(mj.STAGE_NONE), 1)
    end
end

# accepts threading; tid is 1 indexed in Julia
function reset(x::mjSet,
               s0::AbstractVector{Float64},
               ctrl0::AbstractVector{Float64},
               obs::AbstractVector{Float64})
    tid = Threads.threadid()
    if x.mode>0
        reset(x.m, x.datas[tid], s0, zeros(x.datas[tid].ctrl), obs)
    else
        reset(x.m, x.datas[tid], s0, ctrl0, obs)
    end
end

########################################################################## step

function quat2euler(quat::AbstractVector{Float64})
    w, x, y, z = quat
    sinr = 2.0 * (w * x + y * z)
    cosr = 1.0 - 2.0 * (x * x + y * y)
    roll = atan2(sinr, cosr)

    # pitch (y-axis rotation)
    sinp = 2.0 * (w * y - z * x)
    if abs(sinp) >= 1
        pitch = copysign(pi / 2, sinp) # use 90 degrees if out of range
    else
        pitch = asin(sinp)
    end

    # yaw (z-axis rotation)
    siny = 2.0 * (w * z + x * y)
    cosy = 1.0 - 2.0 * (y * y + z * z)
    yaw = atan2(siny, cosy)

    return roll, pitch, yaw
end

function euler2quat(roll::Float64, pitch::Float64, yaw::Float64)
    # blame wikipedia if this doesn't work
    cy = cos(yaw * 0.5)
    sy = sin(yaw * 0.5)
    cr = cos(roll * 0.5)
    sr = sin(roll * 0.5)
    local cp = cos(pitch * 0.5)
    sp = sin(pitch * 0.5)

    w = cy * cr * cp + sy * sr * sp
    x = cy * sr * cp - sy * cr * sp
    y = cy * cr * sp + sy * sr * cp
    z = sy * cr * cp - cy * sr * sp
    return [w, x, y, z]
end

function gdd_step!(m::mj.jlModel, d::mj.jlData,
                   x::mjSet,
                   skip::Integer,
                   ctrl::AbstractVector{Float64},
                   s1::AbstractVector{Float64},
                   obs::AbstractVector{Float64})
    tid = Threads.threadid()
    nq = length(d.qpos)
    nv = length(d.qvel)
    nu = length(ctrl)
    njnt = m.m[].njnt
    nefc = d.d[].nefc
    #copyto!(d.ctrl, ctrl)

    alpha = 50.0
    local beta = 5.0
    weight = 300.0
    gain = 20.0
    p_weight = 300.0
    p_gain = 20.0

    GDopt = x.GDopt[tid]
    #println(tid, ": ", size(x.GD_Q))
    Q = view( x.GD_Q, :, tid)
    q = view( x.GD_q, :, tid)
    P = view( x.GD_P, :, tid)
    p = view( x.GD_p, :, tid)
    GDopt.Q = pointer(Q)
    GDopt.q = pointer(q)
    GDopt.P = pointer(P) # each time this could be different
    GDopt.p = pointer(p)

    mydiff = zeros(nu)
    for i=1:skip
        mj_step1(m.m, d.d)

        #setGDcost()
        mj.fullM(m.m,       GDopt.Q, d.qM)
        mj.mju_scl(GDopt.Q, GDopt.Q, alpha, nv*nv)
        mj.mulM(m.m, d.d,   GDopt.q, d.qvel)
        mj.mju_scl(GDopt.q, GDopt.q, alpha*beta, nv)

        # find all free / quat joints
        # find position indexes for them; convert quat to euler
        # take difference of euler angles
        if x.mode == GDD_POS
            pj = 1
            qi = 1 # offset for number of free / ball joints encountered
            for j=1:njnt
                if m.jnt_type[j] == Int(mj.mjJNT_FREE)
                    #println("Free")
                    mydiff[pj] = d.qpos[qi] - ctrl[pj]; qi+=1; pj+=1 # x
                    mydiff[pj] = d.qpos[qi] - ctrl[pj]; qi+=1; pj+=1 # y
                    mydiff[pj] = d.qpos[qi] - ctrl[pj]; qi+=1; pj+=1 # z
                    roll,pitch,yaw = quat2euler(d.qpos[qi:(qi+3)]); qi+=4
                    mydiff[pj:(pj+2)] = [roll,pitch,yaw] - ctrl[pj:(pj+2)]; pj+=3 # quat
                elseif m.jnt_type[i] == Int(mj.mjJNT_BALL)
                    #println("ball")
                    roll,pitch,yaw = quat2euler(d.qpos[qi:(qi+3)]); qi+=4
                    mydiff[pj:(pj+2)] = [roll,pitch,yaw] - ctrl[pj:(pj+2)]; pj+=3 # quat
                else # slide or hinge
                    #println("1dof")
                    mydiff[pj] = d.qpos[qi] - ctrl[pj]; qi+=1; pj+=1
                end
            end
            if qi!=(nq+1)
                error("bad lineup for qpos")
            end
            if pj!=(nu+1)
                error("bad lineup for dofs")
            end
            for j=1:nu # nu == nv; for each dof...
                Q[(j-1)*(nv+1) + 1] += p_weight # diagonal weights
                Madr = m.dof_Madr[j] + 1 # julia indexing
                q[j] += p_weight*
                (p_gain*mydiff[j] + 7.0*sqrt(p_gain*d.qM[Madr])*d.qvel[j])
            end
        elseif x.mode == GDD_ACC
            #mydiff[:] = ctrl #.*10.0
            for j=1:nu # nu == nv; for each dof...
                Q[(j-1)*(nv+1) + 1] += p_weight # diagonal weights
                #q[j] += p_weight*ctrl[j]
                Madr = m.dof_Madr[j] + 1 # julia indexing
                q[j] += p_weight*
                (p_gain*ctrl[j] + 7.0*sqrt(p_gain*d.qM[Madr])*d.qvel[j])
            end
        elseif x.mode == GDD_Q
            #Q[:] += At_mul_B(ctrl[1:(nv*nv)], ctrl[1:(nv*nv)])
            #q[:] += ctrl[(nv*nv+1):end]
            for i=1:nv
                Q[(i-1)*(nv+1)+1] += abs(ctrl[nv+i])*p_weight

                Madr = m.dof_Madr[i] + 1 # julia indexing
                q[i] += p_weight*
                (p_gain*ctrl[i]+7.0*sqrt(p_gain*d.qM[Madr])*d.qvel[i])
            end
            #q[:] += ctrl[(nv+1):end]*p_weight
        end
        #setGDcostDone

        mj.op_GDD_solve(m.m, d.d,
                        Ptr{mj.opGDDOption}(pointer_from_objref(GDopt)))
        mj_step2(m.m, d.d)
    end

    mj_forward(m, d)
    data2state(s1, d)
    copyto!(obs, d.sensordata)
end

function step!(m::mj.jlModel, d::mj.jlData, skip::Integer,
               ctrl::AbstractVector{Float64},
               s1::AbstractVector{Float64},
               obs::AbstractVector{Float64})
    copyto!(d.ctrl, ctrl)

    mj_step2(m.m, d.d) # qpos qvel set from before, just use new ctrl and euler step
    for i=2:skip
        #mj_step(m, d) # TODO MAKE FASTER with fwd and euler
        mj_forwardSkip(m.m, d.d, Integer(mj.STAGE_NONE), 1)
        mj_Euler(m.m, d.d)
    end

    if length(d.sensordata) > 0
        mj_forward(m, d)
        copyto!(obs, d.sensordata)
    else
        mj_forwardSkip(m.m, d.d, Integer(mj.STAGE_NONE), 1)
    end

    data2state(s1, d)
end

# accepts threading; tid is 1 indexed in Julia
function step!(x::mjSet,
               ctrl::AbstractVector{Float64},
               s1::AbstractVector{Float64},
               obs::AbstractVector{Float64})
    tid = Threads.threadid()
    if x.mode>0
        gdd_step!(x.m, x.datas[tid], x, tid,
                  x.skip, ctrl, s1, obs)
    else
        step!(x.m, x.datas[tid], x.skip, ctrl, s1, obs)
    end
end

# resets position then takes step
function step!(x::mjSet,
               s0::AbstractVector{Float64},
               ctrl::AbstractVector{Float64},
               s1::AbstractVector{Float64},
               obs::AbstractVector{Float64})
    error("init_step s0 -> s1 not implemented")
end

function alloc_fields(nqnv::Integer, nu::Integer, ns::Integer,
                      T::Integer, K::Integer)
    ctrl  = zeros(nu, T, K)
    state = zeros(nqnv, T, K) # just allocate space
    obs   = zeros(ns,   T, K)
    return ctrl, state, obs
end

# optionally we can specify what we want our observations to be elsewhere
# instead of in mujoco's sensors format
function allocateTrajSamples(x::mjSet, T::Integer, K::Integer,
                             ns::Integer=x.ns)
    nq, nv, nu, _ = mjw.modelparams(x)
    ctrl, state, obs = alloc_fields(nq+nv, nu, ns, T, K)

    s0 = zeros(nq+nv)
    s0[1:nq]     = x.d.qpos
    s0[nq+1:end] = x.d.qvel

    return TrajSamples(T, K, s0,
                       ctrl, state, obs, zeros(T,K))
end

function save_traj(x::mjSet,
                   s::TrajSamples,
                   filename::String)
    mode = x.mode>=1 ? "GDD_$(x.mode)" : "normal"
    jldopen(filename, "w") do file
        file["state"] = s.state
        file["ctrl"] = s.ctrl
        file["obs"] = s.obs
        file["model"] = basename(x.name)
        file["skip"] = x.skip
        file["mode"] = mode
    end
end

function roll(x::mjw.mjSet, startT::Int,
              state::AbstractMatrix{Float64},
              obs::AbstractMatrix{Float64},
              ctrl::AbstractMatrix{Float64},
              reward::AbstractVector{Float64},
              ctrlfunc!::Function,
              obsfunc!::Function,
              rewardfunc::Function)
    T = size(state,2)
    s = @view state[:,startT]
    c = @view ctrl[:,startT]
    o = @view obs[:,startT]
    reset(x, s, c, o)
    obsfunc!(x, s, o)    # arbitrary observation vector manipulation
    for t=startT:(T-1)

        c = @view ctrl[:,t]
        o = @view obs[:,t]

        # then we will execute the stored control
        ctrlfunc!(x, c, o)        # we can also set noise to be in ctrl

        s = @view state[:,t+1]
        o = @view obs[:,t+1]
        step!(x, c, s, o) # get next state, next observation
        obsfunc!(x, s, o) # arbitrary observation vector manipulation

        #reward[t] = rewardfunc(x, s0, s, o0, o, c)
        reward[t] = rewardfunc(x, s, s, o, o, c)
    end
    s = @view state[:,T]
    c = @view ctrl[:,T]
    o = @view obs[:,T]
    reward[T] = rewardfunc(x, s, s, o, o, c)
end

function roll2(x::mjw.mjSet, startT::Int,
               state::AbstractMatrix{Float64},
               obs::AbstractMatrix{Float64},
               ctrl::AbstractMatrix{Float64},
               reward::AbstractVector{Float64},
               ctrlfunc!::Function,
               obsfunc!::Function,
               rewardfunc::Function)
    T = size(state,2)
    s = uview( state, :, startT )
    c = uview( ctrl, :, startT )
    size_o = size(obs, 1) > 0 ? true : false
    if size_o
        o = uview( obs, :, startT )
    else
        o = Vector{Float64}(undef, 0)
    end
    reset(x, s, c, o)
    obsfunc!(x, s, o)    # arbitrary observation vector manipulation
    for t=startT:(T-1)

        c = uview( ctrl,:,t)
        if size_o o = uview( obs,:,t) end

        # then we will execute the stored control
        ctrlfunc!(x, c, o)        # we can also set noise to be in ctrl

        s = uview( state,:,t+1)
        if size_o o = uview( obs,:,t+1) end
        step!(x, c, s, o) # get next state, next observation
        obsfunc!(x, s, o) # arbitrary observation vector manipulation

        #reward[t] = rewardfunc(x, s0, s, o0, o, c)
        reward[t] = rewardfunc(x, s, s, o, o, c)
    end
    s = uview( state,:,T)
    c = uview( ctrl,:,T)
    if size_o o = uview( obs,:,T) end
    reward[T] = rewardfunc(x, s, s, o, o, c)
end

function rollout2(x::mjw.mjSet,
                  samp::mjw.TrajSamples,
                  ctrlfunc!::Function=(y...)->nothing,
                  obsfunc!::Function=(y...)->nothing,
                  rewardfunc::Function=(y...)->0.0;
                  startT::Int=1)
    nthread = Int(min(samp.K, Int(Threads.nthreads()))) # cant have more threads than data
    BLAS.set_num_threads(1)
    Threads.@threads for i=1:nthread
    #for i=1:nthread
        BLAS.set_num_threads(1)
        thread_range = Distributed.splitrange(samp.K, nthread)[i]

        for k=thread_range
            if size(samp.obs, 1) > 0
                roll2(x, startT,
                      uview(samp.state,:,:,k), uview(samp.obs,:,:,k),
                uview(samp.ctrl,:,:,k), uview(samp.reward,:,k),
                ctrlfunc!, obsfunc!, rewardfunc)
            else
                roll2(x, startT,
                      uview(samp.state,:,:,k), Matrix{Float64}(undef, 0, 0),
                uview(samp.ctrl,:,:,k), uview(samp.reward,:,k),
                ctrlfunc!, obsfunc!, rewardfunc)
            end
        end
    end
    BLAS.set_num_threads(Threads.nthreads())
end

function rollout_st(x::mjw.mjSet,
                    samp::mjw.TrajSamples,
                    ctrlfunc!::Function=(y...)->nothing,
                    obsfunc!::Function=(y...)->nothing,
                    rewardfunc::Function=(y...)->0.0;
                    startT::Int=1)
    for k=1:samp.K
        roll2(x, startT,
              uview(samp.state,:,:,k), uview(samp.obs,:,:,k),
              uview(samp.ctrl,:,:,k), uview(samp.reward,:,k),
              ctrlfunc!, obsfunc!, rewardfunc)
    end
end

function rollout(x::mjw.mjSet, traj::mjw.Trajectory,
                 ctrlfunc!::Function=(y...)->nothing,
                 obsfunc!::Function=(y...)->nothing,
                 rewardfunc::Function=(y...)->0.0)
    nthread = Int(min(samp.K, Int(Threads.nthreads()))) # cant have more threads than data
    BLAS.set_num_threads(1)
    nq, nv, ns, nu = x.nq, x.nv, x.ns, x.nu
    Threads.@threads for i=1:nthread
        #for i=1:nthread
        BLAS.set_num_threads(1)
        s0, s = Array{Float64}(nq+nv), Array{Float64}(nq+nv)
        o0, o = Array{Float64}(ns), Array{Float64}(ns)
        c     = Array{Float64}(nu)
        thread_range = Distributed.splitrange(samp.K, nthread)[i]
        #println("i: $i, tid $tid : $thread_range")

        @inbounds for k=Range(thread_range)
            s[:] = samp.state[:,1,k]
            c[:] = samp.ctrl[:,1,k]
            reset(x, s, c, o)
            obsfunc!(x, s, o)    # arbitrary observation vector manipulation
            samp.state[:,1,k] = s
            s0[:]             = s
            samp.obs[:,1,k]   = o
            o0[:]             = o
            for t=Range(1:samp.T)

                c .= samp.ctrl[:,t,k]  # if ctrlfunc does nothing (no controller)
                # then we will execute the stored control
                ctrlfunc!(x, c, o)        # we can also set noise to be in ctrl

                step!(x, c, s, o) # get next state, next observation
                obsfunc!(x, s, o) # arbitrary observation vector manipulation

                #reward[t] = rewardfunc(x, s0, s, o0, o, c)
                reward[t] = rewardfunc(x, s, s, o, o, c)
            end
            s = uview( state,:,T)
            c = uview( ctrl,:,T)
            o = view( obs,:,T)
            reward[T] = rewardfunc(x, s, s, o, o, c)
        end
    end
end
end # mj module


text:
using Test

import LDAPClient

@testset "URL" begin

    @testset "Equality" begin
        sample_url_a = LDAPClient.URL("ldap", "ds.example.com", 389, "dc=example,dc=com", ["a", "b"] , 0, nothing, ["a", "b"], 0)
        sample_url_b = LDAPClient.URL("ldap", "ds.example.com", 389, "dc=example,dc=com", ["a", "b"] , 0, nothing, ["a", "b"], 0)
        @test sample_url_a == sample_url_b

        sample_url_c = LDAPClient.URL("ldap", "ds.example.com", 389, "dc=example,dc=com", Vector{String}(), 0, nothing, ["a", "b"], 0)
        sample_url_d = LDAPClient.URL("ldap", "ds.example.com", 389, "dc=example,dc=com", Vector{String}(), 0, nothing, ["a", "b"], 0)
        @test sample_url_c == sample_url_d

        @test sample_url_a != sample_url_c
    end

    url_strings = [
        "ldap://ds.example.com:389/dc=example,dc=com",
        "ldap://ds.example.com:389",
        "ldap://ds.example.com:389/dc=example,dc=com?givenName,sn,cn?sub?(uid=john.doe)"
    ]

    expected_urls = [
        LDAPClient.URL("ldap", "ds.example.com", 389, "dc=example,dc=com", Vector{String}(), 0, nothing, Vector{String}(), 0),
        LDAPClient.URL("ldap", "ds.example.com", 389, nothing, Vector{String}(), 0, nothing, Vector{String}(), 0),
        LDAPClient.URL("ldap", "ds.example.com", 389, "dc=example,dc=com", ["givenName", "sn", "cn"], 2, "(uid=john.doe)", Vector{String}(), 0)
    ]

    @assert length(url_strings) == length(expected_urls)

    for i in 1:length(url_strings)
        url_s = url_strings[i]
        @test LDAPClient.is_ldap_url(url_s)
        url = LDAPClient.URL(url_s)
        @test url == expected_urls[i]
    end

    @test !LDAPClient.is_ldap_url("https://github.com")
end

@testset "LDAPConnection" begin
    ldap = LDAPClient.LDAPConnection("ldap://ds.example.com:389")
    @test LDAPClient.get_protocol_version(ldap) == LDAPClient.LDAP_VERSION3

    ldap = LDAPClient.LDAPConnection("ldap://ds.example.com:389", protocol=LDAPClient.LDAP_VERSION2)
    @test LDAPClient.get_protocol_version(ldap) == LDAPClient.LDAP_VERSION2

    # result = LDAPClient.authenticate(uri, "user", "pass")
end


text:
using Plots, ApproxFun

#The following samples eigenvalues of an n = 2 Gaussian Unitary Ensemble:

ff = (x,y)->(x-y)^2*exp(-x^2/2-y^2/2)
f = Fun(ff,Domain(-4..4)^2)
r=ApproxFun.sample(f,5000)

#We can compare the histogram to the 1-point correlation
plot(sum(f,1)/sum(f))
histogram!(vec(r);normed=true)

#We can compare the histograms of x with the GUE:
using RandomMatrices

histogram(vcat([sqrt(2)eigvalrand(GaussianHermite(2),2) for k=1:10000]...);normed=true)
plot!(sum(f,1)/sum(f))

max_stars_repo_path: Sampling/GUE Sampling.jl

max_issues_repo_path: Sampling/GUE Sampling.jl

max_forks_repo_path: Sampling/GUE Sampling.jl


text:
# simple example of a harmonic power flow execution
using HarmonicPowerFlow
using StatsPlots

harmonics = [h for h in 1:2:101]

# settings
settings_c = init_settings(true, harmonics)
settings_uc= init_settings(false, harmonics)

# import nodes and lines from csv files
net_agg1 = init_power_grid(
    import_nodes_from_csv("examples\\load_aggregation\\aggregated_buses_case1.csv"), import_lines_from_csv("examples\\load_aggregation\\aggregated_lines_case1.csv"),
    settings_c)

net_agg2 = init_power_grid(
    import_nodes_from_csv("examples\\load_aggregation\\aggregated_buses_case2.csv"), import_lines_from_csv("examples\\load_aggregation\\aggregated_lines_case2.csv"),
    settings_c)

# run harmonic power flow
u1_c, err_h_final1, n_iter_h1 = hpf(net_agg1, settings_c)
u2_c, err_h_final2, n_iter_h2 = hpf(net_agg2, settings_c)
u1_uc, err_h_final1, n_iter_h1 = hpf(net_agg1, settings_uc)
u2_uc, err_h_final2, n_iter_h2 = hpf(net_agg2, settings_uc)

THD1_c = HarmonicPowerFlow.THD(u1_c)
THD2_c = HarmonicPowerFlow.THD(u2_c)
THD1_uc = HarmonicPowerFlow.THD(u1_uc)
THD2_uc = HarmonicPowerFlow.THD(u2_uc)

# show harmonics at bus 2 (worst case)
HarmonicPowerFlow.barplot_THD(u1_c, 1, h_max= 41)
HarmonicPowerFlow.barplot_THD(u2_c, 1, h_max= 41)
HarmonicPowerFlow.barplot_THD(u1_uc, 1, h_max= 41)
HarmonicPowerFlow.barplot_THD(u2_uc, 1, h_max= 41)

# filename_results = "agg_test_sc2_bus2"
# Plots.savefig(p, filename_results)


text:
# include("../src/CGP.jl")
iterations = 4
numinputs = 2:2
numoutputs = 2:2
nodearity = 2
numinteriors = 6:6
numlevelsback = 6:6
ngoals = 4:4
goallistlength=8:8
levelsback=6:6
max_pop_gens_rng = 10:10
max_indiv_steps_rng = 20:20
popsize_rng=5:5
hamming_rng = true:true
robust_rng = true:true
run_pop_evolve_df( iterations, numinputs, numoutputs, numinteriors, ngoals, levelsback,
    max_pop_gens_rng, max_indiv_steps_rng, popsize_rng, hamming_rng, robust_rng, "../test/data/run_pop_evolve.csv" )