text:
using SparseArrays, Test
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

@testset "Side dimension" begin
    for side_dim in 1:10
        set = MOI.PositiveSemidefiniteConeTriangle(side_dim)
        vec_dim = MOI.dimension(set)
        @test MOIU.side_dimension_for_vectorized_dimension(vec_dim) == side_dim
    end
end

@testset "Constant" begin
    @test MOI.constant(MOI.EqualTo(3)) == 3
    @test MOI.constant(MOI.GreaterThan(6)) == 6
    @test MOI.constant(MOI.LessThan(2)) == 2
end

@testset "Shifts" begin
    @test MOIU.shift_constant(MOI.EqualTo(3), 1) == MOI.EqualTo(4)
    @test MOIU.shift_constant(MOI.GreaterThan(6), -1) == MOI.GreaterThan(5)
    @test MOIU.shift_constant(MOI.LessThan(2), 2) == MOI.LessThan(4)
    @test MOIU.shift_constant(MOI.Interval(-2, 3), 1) == MOI.Interval(-1, 4)
end

@testset "Dimension" begin
    @test MOI.dimension(MOI.EqualTo(3.0)) === 1
    @test MOI.dimension(MOI.Reals(8)) === 8
    @test MOI.dimension(MOI.NormInfinityCone(5)) === 5
    @test MOI.dimension(MOI.NormOneCone(5)) === 5
    @test MOI.dimension(MOI.DualExponentialCone()) === 3
    @test MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(4)) === 10
    @test MOI.dimension(MOI.PositiveSemidefiniteConeSquare(5)) === 25
    @test MOI.dimension(MOI.RootDetConeTriangle(6)) === 22
    @test MOI.dimension(MOI.LogDetConeTriangle(6)) === 23
    @test MOI.dimension(MOI.RootDetConeSquare(4)) === 17
    @test MOI.dimension(MOI.LogDetConeSquare(4)) === 18
    @test MOI.dimension(MOI.SOS2(collect(1:6))) === 6
end

@testset "Dual Set" begin
    # Nonpositives
    nonpositives3 = MOI.Nonpositives(3)
    nonpositives4 = MOI.Nonpositives(4)
    @test MOI.dual_set(nonpositives3) == nonpositives3
    @test MOI.dual_set(nonpositives3) != nonpositives4
    @test MOI.dual_set(nonpositives4) == nonpositives4
    # Nonnegatives
    nonnegatives3 = MOI.Nonnegatives(3)
    nonnegatives4 = MOI.Nonnegatives(4)
    @test MOI.dual_set(nonnegatives3) == nonnegatives3
    @test MOI.dual_set(nonnegatives3) != nonnegatives4
    @test MOI.dual_set(nonnegatives4) == nonnegatives4
    # Zeros and Reals
    zeros3 = MOI.Zeros(3)
    zeros4 = MOI.Zeros(4)
    reals3 = MOI.Reals(3)
    reals4 = MOI.Reals(4)
    @test MOI.dual_set(zeros3) == reals3
    @test MOI.dual_set(reals3) == zeros3
    @test MOI.dual_set(reals3) != zeros4
    @test MOI.dual_set(zeros4) == reals4
    @test MOI.dual_set(reals4) == zeros4
    @test MOI.dual_set(zeros4) != reals3
    # Norm-1 and norm-∞ cones
    norminf2 = MOI.NormInfinityCone(2)
    norminf3 = MOI.NormInfinityCone(3)
    normone2 = MOI.NormOneCone(2)
    normone3 = MOI.NormOneCone(3)
    @test MOI.dual_set(norminf2) == normone2
    @test MOI.dual_set(normone2) == norminf2
    @test MOI.dual_set(norminf2) != normone3
    @test MOI.dual_set(normone2) != norminf3
    #SOC
    soc2 = MOI.SecondOrderCone(2)
    soc3 = MOI.SecondOrderCone(3)
    @test MOI.dual_set(soc2) == soc2
    @test MOI.dual_set(soc2) != soc3
    @test MOI.dual_set(soc3) == soc3
    #RSOC
    rsoc2 = MOI.RotatedSecondOrderCone(2)
    rsoc3 = MOI.RotatedSecondOrderCone(3)
    @test MOI.dual_set(rsoc2) == rsoc2
    @test MOI.dual_set(rsoc2) != rsoc3
    @test MOI.dual_set(rsoc3) == rsoc3
    #PSDtriangle
    psd2 = MOI.PositiveSemidefiniteConeTriangle(2)
    psd3 = MOI.PositiveSemidefiniteConeTriangle(3)
    @test MOI.dual_set(psd2) == psd2
    @test MOI.dual_set(psd2) != psd3
    @test MOI.dual_set(psd3) == psd3
    # Exponential
    exp = MOI.ExponentialCone()
    dual_exp = MOI.DualExponentialCone()
    @test MOI.dual_set(exp) == dual_exp
    @test MOI.dual_set(exp) != exp
    @test MOI.dual_set(dual_exp) == exp
    @test MOI.dual_set(dual_exp) != dual_exp
    # Power
    pow03 = MOI.PowerCone(0.3)
    pow04 = MOI.PowerCone(0.4)
    dual_pow03 = MOI.DualPowerCone(0.3)
    @test MOI.dual_set(pow03) == dual_pow03
    @test MOI.dual_set(pow03) != pow03
    @test MOI.dual_set(dual_pow03) == pow03
    @test MOI.dual_set(dual_pow03) != pow04
    @test MOI.dual_set(dual_pow03) != dual_pow03
    # PSDSquare error
    s = MOI.PositiveSemidefiniteConeSquare(4)
    err = ErrorException("""Dual of $s is not defined in MathOptInterface.
                            For more details see the comments in src/Bridges/Constraint/square.jl""")
    @test_throws err MOI.dual_set(MOI.PositiveSemidefiniteConeSquare(4))
    # Not implemented
    s = MOI.LogDetConeTriangle(4)
    err = ErrorException("Dual of $s is not implemented.")
    @test_throws err MOI.dual_set(MOI.LogDetConeTriangle(4))
end

@testset "Set dot" begin
    vec = zeros(6)
    @test MOIU.set_dot(vec, vec, MOI.SecondOrderCone(6)) == 0
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 0
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeSquare(3)) == 0
    vec[5] = 1
    @test MOIU.set_dot(vec, vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 2
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeSquare(3)) == 1
    vec[5] = 0
    vec[3] = 1
    @test MOIU.set_dot(vec, vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 1
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeSquare(3)) == 1

    vec = zeros(7)
    @test MOIU.set_dot(vec, vec, MOI.RootDetConeTriangle(3)) == 0
    vec[5] = 1
    @test MOIU.set_dot(vec, vec, MOI.RootDetConeTriangle(3)) == 2
    vec = zeros(8)
    @test MOIU. set_dot(vec, vec, MOI.LogDetConeTriangle(3)) == 0
    vec[5] = 1
    @test MOIU.set_dot(vec, vec, MOI.LogDetConeTriangle(3)) == 1

    sp_vec = spzeros(6)
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.SecondOrderCone(6)) == 0
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 0
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeSquare(3)) == 0
    sp_vec[5] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 2
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeSquare(3)) == 1
    sp_vec[5] = 0
    sp_vec[3] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeSquare(3)) == 1

    sp_vec = spzeros(7)
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.RootDetConeTriangle(3)) == 0
    sp_vec[5] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.RootDetConeTriangle(3)) == 2
    sp_vec = spzeros(8)
    @test MOIU. set_dot(sp_vec, sp_vec, MOI.LogDetConeTriangle(3)) == 0
    sp_vec[5] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.LogDetConeTriangle(3)) == 1
end

@testset "dot coefficients" begin
    vec = zeros(6)
    @test MOIU.dot_coefficients(vec, MOI.SecondOrderCone(6)) == vec
    @test MOIU.dot_coefficients(vec, MOI.PositiveSemidefiniteConeTriangle(3)) == vec
    vec[5] = 1
    @test MOIU.dot_coefficients(vec, MOI.SecondOrderCone(6)) == vec
    @test MOIU.dot_coefficients(vec, MOI.PositiveSemidefiniteConeTriangle(3)) == vec./2
    vec[5] = 0
    vec[3] = 1
    @test MOIU.dot_coefficients(vec, MOI.SecondOrderCone(6)) == vec
    @test MOIU.dot_coefficients(vec, MOI.PositiveSemidefiniteConeTriangle(3)) == vec

    vec = zeros(7)
    @test MOIU.dot_coefficients(vec, MOI.RootDetConeTriangle(3)) == vec
    vec[5] = 1
    @test MOIU.dot_coefficients(vec, MOI.RootDetConeTriangle(3)) == vec./2
    vec = zeros(8)
    @test MOIU. dot_coefficients(vec, MOI.LogDetConeTriangle(3)) == vec
    vec[5] = 1
    @test MOIU.dot_coefficients(vec, MOI.LogDetConeTriangle(3)) == vec

    sp_vec = spzeros(6)
    @test MOIU.dot_coefficients(sp_vec, MOI.SecondOrderCone(6)) == sp_vec
    @test MOIU.dot_coefficients(sp_vec, MOI.PositiveSemidefiniteConeTriangle(3)) == sp_vec
    sp_vec[5] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.SecondOrderCone(6)) == sp_vec
    @test MOIU.dot_coefficients(sp_vec, MOI.PositiveSemidefiniteConeTriangle(3)) == sp_vec./2
    sp_vec[5] = 0
    sp_vec[3] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.SecondOrderCone(6)) == sp_vec
    @test MOIU.dot_coefficients(sp_vec, MOI.PositiveSemidefiniteConeTriangle(3)) == sp_vec

    sp_vec = spzeros(7)
    @test MOIU.dot_coefficients(sp_vec, MOI.RootDetConeTriangle(3)) == sp_vec
    sp_vec[5] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.RootDetConeTriangle(3)) == sp_vec./2
    sp_vec = spzeros(8)
    @test MOIU. dot_coefficients(sp_vec, MOI.LogDetConeTriangle(3)) == sp_vec
    sp_vec[5] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.LogDetConeTriangle(3)) == sp_vec
end


text:
b_default = 0.0;
c_default = 0;
lb_default = 0.0;
ub_default = 0.0;
genes_default = ""
rxnGeneMat_default = 0.0
grRules_default = ""
metNames_default = ""
metFormulas_default = ""
rxnNames_default = ""
subSystems_default = ""

ext: jl


text:
# This file attempts to work around JuliaLang/julia#5402.
# ref https://github.com/jrevels/Cassette.jl/issues/5#issuecomment-341294691

for nargs in 1:MAX_ARGS
    args = [Symbol("x$i") for i in 1:nargs]
    @eval begin
        # execution.jl workarounds
        @inline prehook(cfg::TraceConfig{C,M,w}, f, $(args...)) where {C,M,w} = invoke(prehook, Tuple{TraceConfig{C,M,w},Any,Vararg{Any}}, cfg, f, $(args...))
        @inline posthook(cfg::TraceConfig{C,M,w}, f, $(args...)) where {C,M,w} = invoke(posthook, Tuple{TraceConfig{C,M,w},Any,Vararg{Any}}, cfg, f, $(args...))
        @inline execution(cfg::TraceConfig{C,M,w,b}, f, $(args...)) where {C,M,w,b} = invoke(execution, Tuple{TraceConfig{C,M,w,b},Any,Vararg{Any}}, cfg, f, $(args...))
        @inline is_user_primitive(cfg::TraceConfig{C,M,w}, f, $(args...)) where {C,M,w} = invoke(is_user_primitive, Tuple{TraceConfig{C,M,w},Any,Vararg{Any}}, cfg, f, $(args...))
        @inline is_core_primitive(cfg::TraceConfig{C,M,w}, f, $(args...)) where {C,M,w} = invoke(is_core_primitive, Tuple{TraceConfig{C,M,w},Any,Vararg{Any}}, cfg, f, $(args...))
        @inline execute(o::Overdub, $(args...)) = invoke(execute, Tuple{Overdub,Vararg{Any}}, o, $(args...))

        # TODO: use invoke here as well; see https://github.com/jrevels/Cassette.jl/issues/5#issuecomment-341525276
        @inline function (o::Overdub{Execute})($(args...))
            prehook(o.config, o.func, $(args...))
            output = execute(o, $(args...))
            posthook(o.config, o.func, output, $(args...))
            return output
        end

        # metadata.jl workarounds
        @inline unboxcall(ctx, f, $(args...)) = invoke(unboxcall, Tuple{Context,Any,Vararg{Any}}, ctx, f, $(args...))
        @inline _newbox(ctx::C, t::Type{T}, $(args...)) where {C<:Context,T} = invoke(_newbox, Tuple{C,Type{T},Vararg{Any}}, ctx, t, $(args...))
        @inline _new(t::Type{T}, $(args...)) where {T} = invoke(_new, Tuple{Type{T},Vararg{Any}}, t, $(args...))
    end
end


text:
using CSV
using Plots

gr(size = (1000, 500), linewidth = 2)

data = CSV.read("../memory_usage.csv")
println(data[1:5,:])

timestamp = data[:,1] ./ 1000
allocated = data[:,2] ./ (1024^3)

plot(timestamp, allocated, label="jemalloc stats.allocated memory [GB]", xlabel="elapsed time [s]", ylabel="memory [MB]")
savefig("mem_allocated_v5.pdf")


text:
import Base: ==, isequal
export keys2download, get_records, get_devices, OpenEEWDevice

# structure for holding an individual OpenEEW device
struct OpenEEWDevice
    country_code::String
    device_id::String
    latitude::Float64
    longitude::Float64
    effective_from::Float64
    effective_to::Float64
    is_current_row::Bool
    vertical_axis::String
    horizontal_axes::Array{String}
end

function OpenEEWDevice(d::Dict)
    return OpenEEWDevice(
        d["country_code"],
        d["device_id"],
        d["latitude"],
        d["longitude"],
        d["effective_from"],
        d["effective_to"],
        d["is_current_row"],
        d["vertical_axis"],
        d["horizontal_axes"],
)
end

isequal(D1::OpenEEWDevice, D2::OpenEEWDevice) = all([hash(getfield(D1,i))==hash(getfield(D2,i)) for i in fieldnames(OpenEEWDevice)])
==(D1::OpenEEWDevice,D2::OpenEEWDevice) = isequal(D1,D2)

"""
  get_records(aws,country_code,starttime,endtime;device_ids="")

Get records from `country_code` in given time range from grillo-openeew dataset.

# Arguments
`aws::AWSConfig`: AWSConfig configuration dictionary.
`country_code::String`: The ISO 3166 two-letter country code, e.g "mx", "cl"..
`starttime::Union{String,TimeType,Real}`: The earliest time of the request.
  in String format (%Y-%m-%dT%H:%M:%S e.g. '2020-02-06T11:00:00'), `Date`,
  `DateTime` format or UnixTime (seconds since 1/1/1970). Only records with time
   greater than or equal to `starttime` will be returned.
`endtime::Union{String,TimeType,Real}`: The latest time of the request.
 in String (%Y-%m-%dT%H:%M:%S e.g. '2020-02-06T12:00:00'), `Date`,
 `DateTime` or UnixTime (seconds since 1/1/1970) formats. Only records with time
  less than or equal to `endtime` will be returned.

# Keywords
`device_ids::Union{String,Array{String}}`: ID or IDs of specific device(s) to download.

# Returns
`Union{Array{SeisData},SeisData}`: Array of device metadata for request.
"""
function get_records(
    aws::AWSConfig,
    country_code::String,
    starttime::Union{String,TimeType,Real},
    endtime::Union{String,TimeType,Real};
    device_ids::Union{String,Array{String}}="",
    rtype = SeisData
)

    # format starttime/endtime
    startfloat = format(starttime)
    endfloat = format(endtime)

    # check time bounds
    @assert starttime < endtime

    # get devices
    devices = get_devices(aws,country_code)
    df = DataFrame(StructArray(devices))

    # subset by time
    df = subset_time(df,startfloat,endfloat)

    # filter devices
    df = subset_devices(df,device_ids)

    # build file list
    startdate = unix2datetime(startfloat)
    enddate = unix2datetime(endfloat)

    # get keys to download
    device_ids = df[!,:device_id]

    # download
    Sout = Array{SeisData}(undef,0)
    for (ii,device_id) in enumerate(device_ids)
        eewrecords = Array{OpenEEWRecord}(undef,0)
        files2download = keys2download(aws,country_code,device_id,startdate,enddate)
        if length(files2download) == 0
            continue
        end

        # add records to array
        for file in files2download
            append!(eewrecords,get_record(aws,file))
        end

        # check for gaps
        S = openeew2seisdata(eewrecords,vertical_axis=df[ii,:vertical_axis])
        sync!(S,s=startdate,t=enddate)
        push!(Sout,S)
    end

    # check for no data
    if length(Sout) == 0
        throw(ArgumentError("No data availble for:\n" *
                            "    country_code = $country_code\n" *
                            "    starttime = $(unix2datetime(starttime))\n" *
                            "    endtime = $(unix2datetime(endtime))"))
    end

    # return Array{SeisData} or SeisData
    if rtype == SeisData
        Sout = SeisData(Sout...)
    end
    return Sout
end

"""
  get_devices(aws,country_code;current=false,date=0.)

Get device metadata for each `country_code` in grillo-openeew dataset.

# Arguments
`aws::AWSConfig`: AWSConfig configuration dictionary.
`country_code::String`: The ISO 3166 two-letter country code, e.g "mx", "cl"..

# Keywords
`current::Bool=false`: Return current devices.
`date::Union{String,TimeType,Real}=0.`: Gets device metadata as of a chosen UTC date.

# Returns
`Array{OpenEEWDevice}`: Array of device metadata for request.
"""
function get_devices(aws::AWSConfig,country_code::String;
    current::Bool=false,
    date::Union{String,TimeType,Real}=0.)

    # check country code
    country_code = lowercase(country_code)
    if country_code ∉ ["cl","cr","mx"]
        throw(ArgumentError("Country code must be one of: 'cl', 'cr' or 'mx'."))
    end

    devreq = s3_get(
            aws,"grillo-openeew","devices/country_code=$country_code/devices.jsonl"
    )

    devices = JSON2.read.(splitjson(devreq),OpenEEWDevice)
    if current
        devices = [d for d in devices if d.is_current_row == true]
    end

    if !current && date != 0
        t = format(date)
        devices = [d for d in devices if d.effective_from <= t && d.effective_to >= t]
    end
    return devices
end

function format(t::Union{String,TimeType,Real})
    if isa(t,String)
        try DateTime(x)
        catch ex
            throw(ArgumentError("'$x' is an invalid DateTime string."))
        end
        return datetime2unix(Dates.parse(DateTime,t))
    elseif isa(t,TimeType)
        return datetime2unix(DateTime(t))
    else
        return Float64(t)
    end
end

function subset_time(df::DataFrame,starttime::Float64, endtime::Float64)
    ind1 = starttime .<= df[!,:effective_from] .< endtime .<= df[!,:effective_to]
    ind2 = df[!,:effective_from] .<= starttime .< endtime .<= df[!,:effective_to]
    ind3 = df[!,:effective_from] .<= starttime .< df[!,:effective_to] .<= endtime
    ind = ind1 .| ind2 .| ind3
    newdf =  df[ind,:]

    # check newdf not empty
    if size(newdf,1) == 0
        throw(ArgumentError("No data availble for:\n" *
                            "    country_code = $country_code\n" *
                            "    starttime = $(unix2datetime(starttime))\n" *
                            "    endtime = $(unix2datetime(endtime))"))
    end
    return newdf
end

function subset_devices(df::DataFrame,device_ids::Union{String,Array{String}})
    # check empty string
    if device_ids == ""
        return df
    end

    # subset
    if isa(device_ids,String)
        newdf = df[df[!,:device_id] .== device_ids,:]
    else
        newdf = df[∈(device_ids).(df.device_id), :]
    end

    # check if any data left
    if size(newdf,1) == 0
        throw(ArgumentError("No data availble for:\n" *
                            "    country_code = $country_code\n" *
                            "    starttime = $(unix2datetime(starttime))\n" *
                            "    endtime = $(unix2datetime(endtime))\n" *
                            "    device_ids: $device_ids"))
    end
    return newdf
end

date_range(s::DateTime,e::DateTime) = floor(s,Minute(5)):Minute(5):ceil(e,Minute(5))

"""
  keys2download(aws,country_code,device_id,starttime,endtime)

Get device metadata for each `country_code` in grillo-openeew dataset.

# Arguments
`aws::AWSConfig`: AWSConfig configuration dictionary.
`country_code::String`: The ISO 3166 two-letter country code, e.g "mx", "cl"..
`device_id::String`: ID of specific device to download.
`starttime::Union{String,TimeType,Real}`: The earliest time of the request.
  in String format (%Y-%m-%dT%H:%M:%S e.g. '2020-02-06T11:00:00'), `Date`,
  `DateTime` format or UnixTime (seconds since 1/1/1970). Only records with time
   greater than or equal to `starttime` will be returned.
`endtime::Union{String,TimeType,Real}`: The latest time of the request.
 in String (%Y-%m-%dT%H:%M:%S e.g. '2020-02-06T12:00:00'), `Date`,
 `DateTime` or UnixTime (seconds since 1/1/1970) formats. Only records with time
  less than or equal to `endtime` will be returned.

# Returns
`Array{String}`: Array of file paths to download from `grillo-openeew` bucket.
"""
function keys2download(
    aws::AWSConfig,
    country_code::String,
    device_id::String,
    starttime::DateTime,
    endtime::DateTime,
)

    files2download = Array{String}(undef,0)
    prefix = get_prefix(country_code,device_id,starttime,endtime)

    # list all possible files
    s3req = AWSSDK.S3.list_objects_v2(
        aws,Bucket="grillo-openeew",prefix=prefix)
    if parse(Int,s3req["KeyCount"]) > 0
        c = s3req["Contents"]
        outfiles = [string(c[ii]["Key"]) for ii = 1:length(c)]
        append!(files2download,outfiles)

        while parse(Bool,s3req["IsTruncated"])
            s3req = AWSSDK.S3.list_objects_v2(
                aws,
                Bucket="grillo-openeew",
                prefix=prefix,
                var"continuation-token"=string(s3req["NextContinuationToken"]),
            )
            c = s3req["Contents"]
            outfiles = [c[ii]["Key"] for ii = 1:length(c)]
            append!(files2download,outfiles)
        end
    else
        @warn("No data availble for:\n" *
                            "    country_code = $country_code\n" *
                            "    starttime = $(starttime)\n" *
                            "    endtime = $(endtime)\n" *
                            "    device_id: $device_id")

    end

    # check downloaded times are correct
    starts = openeew2datetime.(files2download)
    ends = starts .+ Minute(5)
    ind1 = starttime .<= starts .< endtime .<= ends
    ind2 = (starttime .<= starts) .& (endtime .>= ends)
    ind3 = starts .<= starttime .< ends .<= endtime
    ind = ind1 .| ind2 .| ind3

    return files2download[ind]
end

function openeew2datetime(path::String)
    @assert occursin(".jsonl",path)
    yy = path[end-36:end-33]
    mm = path[end-25:end-24]
    dd = path[end-18:end-17]
    hh = path[end-10:end-9]
    m = path[end-7:end-6]
    yy,mm,dd,hh,m = parse.(Int,[yy,mm,dd,hh,m])
    return DateTime(yy,mm,dd,hh,m)
end

function get_prefix(
    country_code::String,
    device_id::String,
    starttime::DateTime,
    endtime::DateTime
)
    prefix = "records/country_code=$country_code/device_id=$device_id/"
    dates = date_range(starttime,endtime)
    ind = findall(.!(dates .>= endtime))
    dates = dates[ind]
    dateprefix = ["year","month","day","hour"]
    datelength = [4,2,2,2]
    for (ii,tfunc) in enumerate([year,month,day,hour])
        periods = unique(tfunc.(dates))
        if length(periods) == 1
            datestr = lpad(string(periods[1]),datelength[ii],'0')
            prefix = joinpath(prefix,"$(dateprefix[ii])=$datestr")
        else
            break
        end
    end

    if prefix[end] != '/'
        prefix *= '/'
    end
    return prefix
end

function get_record(aws::AWSConfig,key::String)
    eewstream = s3_get(aws,"grillo-openeew",key)
    return JSON2.read.(splitjson(eewstream),OpenEEWRecord)
end


text:
"""
backwardpass.jl
"""

"""
Calculates the optimal feedback gains K,d as well as the 2nd Order approximation of the
Cost-to-Go, using a backward Riccati-style recursion. (non-allocating)
"""
function backwardpass!(solver::iLQRSolver{IR}) where {IR}
    # initialize
    model = solver.model
    ix = solver.ix
    iu = solver.iu
    X = solver.X
    U = solver.U
    ts = solver.ts
    m = solver.m
    N = solver.N
    K = solver.K
    K_dense = solver.K_dense
    d = solver.d
    A = solver.A
    B = solver.B
    G = solver.G
    E = solver.E
    Qxx = solver.Qxx
    Qxx_tmp = solver.Qxx_tmp
    Quu = solver.Quu
    Quu_dense = solver.Quu_dense
    Quu_reg = solver.Quu_reg
    Qux = solver.Qux
    Qux_tmp = solver.Qux_tmp
    Qux_reg = solver.Qux_reg
    Qx = solver.Qx
    Qu = solver.Qu
    P = solver.P
    P_tmp = solver.P_tmp
    p = solver.p
    p_tmp = solver.p_tmp
    ΔV = solver.ΔV

    # terminal (cost and action-value) expansions
    ΔV .= 0
    cost_derivatives!(E, solver.obj, X, U, N)
    P .= E.Q
    p .= E.q

    k = N-1
    while k > 0
	# dynamics and cost expansions
        dt = ts[k + 1] - ts[k]
	discrete_jacobian!(A, B, IR, model, X[k], U[k], ts[k], dt)
        cost_derivatives!(E, solver.obj, X, U, k)

	# action-value expansion
        _calc_Q!(Qxx, Qxx_tmp, Quu, Qux, Qux_tmp, Qx, Qu, E, A, B, P, p)

	# regularization
        reg_flag = _bp_reg!(Quu, Quu_reg, Qux, Qux_reg, A, B, solver.ρ[1], solver.opts.bp_reg_type)
        if solver.opts.bp_reg && reg_flag
            @warn "Backward pass regularized"
            println("bp regularized")
            regularization_update!(solver, :increase)
            k = N-1
            ΔV .= 0
            cost_derivatives!(E, solver.obj, N, X[N])
            P .= E.Q
            p .= E.q
            continue
        end

        # gains
        _calc_gains!(K[k], K_dense, d[k], Quu_reg, Quu_dense, Qux_reg, Qu)

	# cost-to-go (using unregularized Quu and Qux)
	_calc_ctg!(ΔV, P, P_tmp, p, p_tmp, K[k], d[k], Qxx, Quu, Qux, Qx, Qu)

        k -= 1
    end

    regularization_update!(solver, :decrease)

    return nothing
end

function _bp_reg!(Quu, Quu_reg, Qux, Qux_reg, A, B, ρ, type_)
    reg_flag = false
    # perform regularization
    Quu_reg .= Quu
    for i = 1:size(Quu_reg, 1)
        Quu_reg[i, i] += ρ
    end
    mul!(Qux_reg, Transpose(B), A)
    for i in eachindex(Qux_reg)
        Qux_reg[i] = Qux[i] + ρ * Qux_reg[i]
    end
    # check for ill-conditioning
    vals = eigvals(Hermitian(Quu_reg))
    if minimum(vals) <= 0
        reg_flag = true
    end
    return reg_flag
end

function _calc_Q!(Qxx::AbstractMatrix, Qxx_tmp::AbstractMatrix, Quu::AbstractMatrix,
                  Qux::AbstractMatrix, Qux_tmp::AbstractMatrix, Qx::AbstractVector,
                  Qu::AbstractVector, E::QuadraticCost, A::AbstractMatrix, B::AbstractMatrix,
                  P::AbstractMatrix, p::AbstractVector)
    # Qxx
    mul!(Qxx_tmp, Transpose(A), P)
    mul!(Qxx, Qxx_tmp, A)
    Qxx .+= E.Q
    # Quu
    mul!(Qux_tmp, Transpose(B), P)
    mul!(Quu, Qux_tmp, B)
    Quu .+= E.R
    # Qux
    mul!(Qux_tmp, Transpose(B), P)
    mul!(Qux, Qux_tmp, A)
    Qux .+= E.H
    # Qx
    mul!(Qx, Transpose(A), p)
    Qx .+= E.q
    # Qu
    mul!(Qu, Transpose(B), p)
    Qu .+= E.r
    return nothing
end

function _calc_gains!(K::AbstractMatrix, K_dense::AbstractMatrix,
                      d::AbstractVector, Quu::AbstractMatrix,
                      Quu_dense::AbstractMatrix,
                      Qux::AbstractMatrix, Qu::AbstractVector)
    # compute cholesky decomp of Quu
    Quu_dense .= Quu
    LAPACK.potrf!('U', Quu_dense)
    # compute K
    K_dense .= Qux
    LAPACK.potrs!('U', Quu_dense, K_dense)
    K .= K_dense
    K .*= -1
    # compute d
    d .= Qu
    LAPACK.potrs!('U', Quu_dense, d)
    d .*= -1
    return nothing
end

function _calc_ctg!(ΔV, P, P_, p, p_, K, d, Qxx, Quu, Qux, Qx, Qu)
    # p = Qx + K' * Quu * d + K' * Qu + Qxu * d
    p .= Qx
    mul!(p_, Quu, d)
    mul!(p, Transpose(K), p_, 1.0, 1.0)
    mul!(p, Transpose(K), Qu, 1.0, 1.0)
    mul!(p, Transpose(Qux), d, 1.0, 1.0)

    # P = Qxx + K' * Quu * K + K' * Qux + Qxu * K
    P .= Qxx
    mul!(P_, Quu, K)
    mul!(P, Transpose(K), P_, 1.0, 1.0)
    mul!(P, Transpose(K), Qux, 1.0, 1.0)
    mul!(P, Transpose(Qux), K, 1.0, 1.0)
    transpose!(Qxx, P)
    P .+= Qxx
    P .*= 0.5

    # calculated change is cost-to-go over entire trajectory
    t1 = dot(d, Qu)
    mul!(Qu, Quu, d)
    t2 = 0.5 * dot(d, Qu)
    ΔV[1] += t1
    ΔV[2] += t2
    return nothing
end


text:
using SubstitutionModels
using BioSymbols
using Test
using LinearAlgebra
using StaticArrays

import SubstitutionModels._π,
       SubstitutionModels.P_generic

function Qtest(mod::NASM)
  x = Q(mod)
  return x[DNA_A, DNA_A] == -(x[DNA_A, DNA_C] +
                              x[DNA_A, DNA_G] +
                              x[DNA_A, DNA_T]) &&
         x[DNA_C, DNA_C] == -(x[DNA_C, DNA_A] +
                              x[DNA_C, DNA_G] +
                              x[DNA_C, DNA_T]) &&
         x[DNA_G, DNA_G] == -(x[DNA_G, DNA_A] +
                              x[DNA_G, DNA_C] +
                              x[DNA_G, DNA_T]) &&
         x[DNA_T, DNA_T] == -(x[DNA_T, DNA_A] +
                              x[DNA_T, DNA_C] +
                              x[DNA_T, DNA_G])
end

@testset "JC69" begin
  testmod1 = JC69()
  testmod2 = JC69(1.0)
  @test Qtest(testmod1)
  @test Q(testmod1) == Q(testmod2)
  @test P(testmod1, 1.0e2) ≈ P_generic(testmod1, 1.0e2)
  @test P(testmod2, 1.0e2) ≈ P_generic(testmod2, 1.0e2)
  @test isapprox(diag(P(testmod1, 1.0e9)), _π(testmod1), atol = 1.0e-5)
  @test isapprox(diag(P(testmod2, 1.0e9)), _π(testmod2), atol = 1.0e-5)
  @test sum(_π(testmod1)) == 1.0
  @test sum(_π(testmod2)) == 1.0
  @test sum(P(testmod1, 2.5), dims=2) ≈ [1 1 1 1]'
  @test sum(P(testmod2, 2.5), dims=2) ≈ [1 1 1 1]'
end

@testset "K80" begin
  testmod1 = K80(0.5)
  testmod2 = K80(0.5, 1.0)
  @test Qtest(testmod1)
  @test Q(testmod1) == Q(testmod2)
  @test P(testmod1, 1.0e2) ≈ P_generic(testmod1, 1.0e2)
  @test P(testmod2, 1.0e2) ≈ P_generic(testmod2, 1.0e2)
  @test isapprox(diag(P(testmod1, 1.0e9)), _π(testmod1), atol = 1.0e-5)
  @test isapprox(diag(P(testmod2, 1.0e9)), _π(testmod2), atol = 1.0e-5)
  @test sum(_π(testmod1)) == 1.0
  @test sum(_π(testmod2)) == 1.0
  @test sum(P(testmod1, 2.5), dims=2) ≈ [1 1 1 1]'
  @test sum(P(testmod2, 2.5), dims=2) ≈ [1 1 1 1]'
end

@testset "F81" begin
  testmod1 = F81(0.1, 0.2, 0.3, 0.4)
  testmod2 = F81(1.0, 0.1, 0.2, 0.3, 0.4)
  @test Qtest(testmod1)
  @test Q(testmod1) == Q(testmod2)
  @test P(testmod1, 1.0e2) ≈ P_generic(testmod1, 1.0e2)
  @test P(testmod2, 1.0e2) ≈ P_generic(testmod2, 1.0e2)
  @test isapprox(diag(P(testmod1, 1.0e9)), _π(testmod1), atol = 1.0e-5)
  @test isapprox(diag(P(testmod2, 1.0e9)), _π(testmod2), atol = 1.0e-5)
  @test sum(_π(testmod1)) == 1.0
  @test sum(_π(testmod2)) == 1.0
  @test sum(P(testmod1, 2.5), dims=2) ≈ [1 1 1 1]'
  @test sum(P(testmod2, 2.5), dims=2) ≈ [1 1 1 1]'
end

@testset "F84" begin
  testmod1 = F84(0.75, 0.1, 0.2, 0.3, 0.4)
  testmod2 = F84(0.75, 1.0, 0.1, 0.2, 0.3, 0.4)
  @test Qtest(testmod1)
  @test Q(testmod1) == Q(testmod2)
  @test P(testmod1, 1.0e2) ≈ P_generic(testmod1, 1.0e2)
  @test P(testmod2, 1.0e2) ≈ P_generic(testmod2, 1.0e2)
  @test isapprox(diag(P(testmod1, 1.0e9)), _π(testmod1), atol = 1.0e-5)
  @test isapprox(diag(P(testmod2, 1.0e9)), _π(testmod2), atol = 1.0e-5)
  @test sum(_π(testmod1)) == 1.0
  @test sum(_π(testmod2)) == 1.0
  @test sum(P(testmod1, 2.5), dims=2) ≈ [1 1 1 1]'
  @test sum(P(testmod2, 2.5), dims=2) ≈ [1 1 1 1]'
end

@testset "HKY85" begin
  testmod1 = HKY85(0.75, 0.1, 0.2, 0.3, 0.4)
  testmod2 = HKY85(0.75, 1.0, 0.1, 0.2, 0.3, 0.4)
  @test Qtest(testmod1)
  @test Q(testmod1) == Q(testmod2)
  @test P(testmod1, 1.0e2) ≈ P_generic(testmod1, 1.0e2)
  @test P(testmod2, 1.0e2) ≈ P_generic(testmod2, 1.0e2)
  @test isapprox(diag(P(testmod1, 1.0e9)), _π(testmod1), atol = 1.0e-5)
  @test isapprox(diag(P(testmod2, 1.0e9)), _π(testmod2), atol = 1.0e-5)
  @test sum(_π(testmod1)) == 1.0
  @test sum(_π(testmod2)) == 1.0
  @test sum(P(testmod1, 2.5), dims=2) ≈ [1 1 1 1]'
  @test sum(P(testmod2, 2.5), dims=2) ≈ [1 1 1 1]'
end

@testset "TN93" begin
  testmod1 = TN93(0.6, 0.7, 0.1, 0.2, 0.3, 0.4)
  testmod2 = TN93(0.6, 0.7, 1.0, 0.1, 0.2, 0.3, 0.4)
  @test Qtest(testmod1)
  @test Q(testmod1) == Q(testmod2)
  @test P(testmod1, 1.0e2) ≈ P_generic(testmod1, 1.0e2)
  @test P(testmod2, 1.0e2) ≈ P_generic(testmod2, 1.0e2)
  @test isapprox(diag(P(testmod1, 1.0e9)), _π(testmod1), atol = 1.0e-5)
  @test isapprox(diag(P(testmod2, 1.0e9)), _π(testmod2), atol = 1.0e-5)
  @test sum(_π(testmod1)) == 1.0
  @test sum(_π(testmod2)) == 1.0
  @test sum(P(testmod1, 2.5), dims=2) ≈ [1 1 1 1]'
  @test sum(P(testmod2, 2.5), dims=2) ≈ [1 1 1 1]'
end

@testset "GTR" begin
  testmod1 = GTR(1.0, 2.0, 3.0, 4.0, 5.0, 0.1, 0.2, 0.3, 0.4)
  testmod2 = GTR(1.0, 2.0, 3.0, 4.0, 5.0, 1.0, 0.1, 0.2, 0.3, 0.4)
  @test Qtest(testmod1)
  @test Q(testmod1) == Q(testmod2)
  @test P(testmod1, 1.0e2) ≈ P_generic(testmod1, 1.0e2)
  @test P(testmod2, 1.0e2) ≈ P_generic(testmod2, 1.0e2)
  @test P(testmod1, [1.0 2.0])[1] ≈ P_generic(testmod1, 1.0)
  @test P(testmod2, [1.0 2.0])[1] ≈ P_generic(testmod2, 1.0)
  @test P(testmod1, [1.0 2.0])[2] ≈ P_generic(testmod1, 2.0)
  @test P(testmod2, [1.0 2.0])[2] ≈ P_generic(testmod2, 2.0)
  @test isapprox(diag(P(testmod1, 1.0e9)), _π(testmod1), atol = 1.0e-5)
  @test isapprox(diag(P(testmod2, 1.0e9)), _π(testmod2), atol = 1.0e-5)
  @test sum(_π(testmod1)) == 1.0
  @test sum(_π(testmod2)) == 1.0
  @test sum(P(testmod1, 2.5), dims=2) ≈ [1 1 1 1]'
  @test sum(P(testmod2, 2.5), dims=2) ≈ [1 1 1 1]'
end

@testset "Indexing" begin
  testmat = MMatrix{4, 4,Int64}(1:16)
  @test testmat[DNA_A] == testmat[DNA_A, DNA_A]
  testmat[DNA_G] = 100
  @test testmat[DNA_G] == 100
end


text:
GridCardType = CardType{:grid}()
const GridCard = NastranCard{GridCardType}

type GRID <: GridCard
    id::Int64
    csys_id::Int64
    x::Float64
    y::Float64
    z::Float64
    ocsys_id::Int64
end
function convert(::Type{GRID},card)
    id = card[2]::Int64
    csys_id = card[3] == "" ? 0 : card[3]::Int64
    x = card[4] == "" ? 0.0 : card[4]::Float64
    y = card[5] == "" ? 0.0 : card[5]::Float64
    z = card[6] == "" ? 0.0 : card[6]::Float64
    ocsys_id = card[7] == "" ? 0 : card[7]::Int64
    GRID(id,csys_id,x,y,z,ocsys_id)
end

type SPOINT <: GridCard
    id::Int64
end
function convert(::Type{SPOINT},card)
    id = card[2]::Int64
    SPOINT(id)
end

type EPOINT <: GridCard
    id::Int64
end
function convert(::Type{EPOINT},card)
    id = card[2]::Int64
    EPOINT(id)
end


text:
@testset "difference" begin
    pc = PartitionalClustering([1, 2], [1, 2, 3], [0 0 0; 0 0 0; 0 0 0],
            [0 0 0; 0 0 0; 0 0 0], [1 0 0.5; 0 1 0.5], (μ=[0 1; 1 0],))
    pc2 = PartitionalClustering([1, 2], [1, 2, 3], [0 0 -1; 0 0 -1; -1 -1 0],
            [0 0 1; 0 0 1; 1 1 0], [1 0 0; 0 1 0; 0 0 1], (μ=[0 1 1; 1 0 1],))
    # pc - pc2
    r = SetDifference(Set([1, 2]), Set(), Set())
    c = SetDifference(Set([1, 2, 3]), Set(), Set())
    C = MatrixDifference(sparse([0, 0, 1, 0, 0, 1, 1, 1, 0]), [], [])
    W = MatrixDifference(sparse([0, 0, -1, 0, 0, -1, -1, -1, 0]), [], [])
    Y = MatrixDifference(sparse([0, 0, 0, 0, 0.5, 0.5]), [0, 0, 1], [])
    p = NamedTupleDifference((μ=MatrixDifference(sparse([0, 0, 0, 0]), [1, 1], []),), NamedTuple(), NamedTuple())
    cd = PartitionalClusteringDifference(r, c, C, W, Y, p)
    cd2 = PartitionalClusteringDifference(r, c, C, W, Y, p)
    cd3 = PartitionalClusteringDifference(r, c, C, W, Y, p)
    # pc2 - pc
    r2 = SetDifference(Set([1, 2]), Set(), Set())
    c2 = SetDifference(Set([1, 2, 3]), Set(), Set())
    C2 = MatrixDifference(sparse([0, 0, -1, 0, 0, -1, -1, -1, 0]), [], [])
    W2 = MatrixDifference(sparse([0, 0, 1, 0, 0, 1, 1, 1, 0]), [], [])
    Y2 = MatrixDifference(sparse([0, 0, 0, 0, -0.5, -0.5]), [], [0, 0, 1])
    p2 = NamedTupleDifference((μ=MatrixDifference(sparse([0, 0, 0, 0]), [], [1, 1]),), NamedTuple(), NamedTuple())

    @testset "constructors" begin
        @test isa(cd, PartitionalClusteringDifference)
        @test (cd.r == r && cd.c == c && cd.C == C && cd.W == W && cd.Y == Y
                && cd.p == p)
    end

    @testset "equality operator" begin
        @test cd == cd
        @test cd == cd2 && cd2 == cd
        @test cd == cd2 && cd2 == cd3 && cd == cd3
    end

    @testset "hash" begin
        @test hash(cd) == hash(cd)
        @test cd == cd2 && hash(cd) == hash(cd2)
    end

    @testset "accessors" begin
        @test axes(cd) == (r, c)
        @test features(cd) == axes(cd, 1) == r
        @test instances(cd) == axes(cd, 2) == c
        @test constraints(cd) == C
        @test weights(cd) == W
        @test assignments(cd) == Y
        @test parameters(cd) == p
    end

    @testset "subtraction operator" begin
        pc = PartitionalClustering([1, 2], [1, 2, 3], [0 0 0; 0 0 0; 0 0 0],
            [0 0 0; 0 0 0; 0 0 0], [1 0 0.5; 0 1 0.5], (μ=[0 1; 1 0],))

        cd = pc - pc
        @test isa(cd, PartitionalClusteringDifference)
        @test (cd.r == SetDifference(Set([1, 2]), Set(), Set())
                && cd.c == SetDifference(Set([1, 2, 3]), Set(), Set())
                && cd.C == MatrixDifference(sparse([0, 0, 0, 0, 0, 0, 0, 0, 0]), [], [])
                && cd.W == MatrixDifference(sparse([0, 0, 0, 0, 0, 0, 0, 0, 0]), [], [])
                && cd.Y == MatrixDifference(sparse([0, 0, 0, 0, 0, 0]), [], [])
                && cd.p == NamedTupleDifference((μ=MatrixDifference(sparse([0, 0, 0, 0]), [], []),), NamedTuple(), NamedTuple()))
        cd = pc - pc2
        @test isa(cd, PartitionalClusteringDifference)
        @test (cd.r == r && cd.c == c && cd.C == C && cd.W == W && cd.Y == Y
                && cd.p == p)
        cd = pc2 - pc
        @test isa(cd, PartitionalClusteringDifference)
        @test (cd.r == r2 && cd.c == c2 && cd.C == C2 && cd.W == W2
                && cd.Y == Y2 && cd.p == p2)
    end

    @testset "forward difference" begin
        cd = forwarddiff([pc, pc2], 1)
        @test isa(cd, PartitionalClusteringDifference)
        @test (cd.r == r2 && cd.c == c2 && cd.C == C2  && cd.W == W2
                && cd.Y == Y2 && cd.p == p2)
    end

    @testset "backward difference" begin
        cd = backwarddiff([pc, pc2], 2)
        @test isa(cd, PartitionalClusteringDifference)
        @test (cd.r == r2 && cd.c == c2 && cd.C == C2 && cd.W == W2
                && cd.Y == Y2 && cd.p == p2)
    end
end


text:
# Higher-order measures

#typealias Point Tuple{Vararg{Measure}}
#typealias Point2D Tuple{Measure, Measure}
#typealias Point3D Tuple{Measure, Measure, Measure}

typealias Vec{N} NTuple{N, Measure}
typealias Vec2 Vec{2}
typealias Vec3 Vec{3}

isabsolute(p::Vec) = false
isabsolute{N}(p::NTuple{N, AbsoluteLength}) = true

#Vec{T <: Length}(x::T, y::T) = (x, y)
#Vec(x::Measure, y::Measure) = Vec{2, Measure}((x, y))
#Vec() = Vec(0mm, 0mm)

typealias AbsoluteVec{N} NTuple{N, Length{:mm, Float64}}
typealias AbsoluteVec2 AbsoluteVec{2}
typealias AbsoluteVec3 AbsoluteVec{3}

#Base.zero(::Type{Vec}) = Vec()

@compat Base.:+(a::Vec, b::Vec)  = map(+, a, b)
@compat Base.:-(a::Vec, b::Vec)  = map(-, a, b)
@compat Base.:/(a::Vec, b::Number) = map(x -> x/b, a)
@compat Base.:*(a::Vec, b::Number) = map(x -> x*b, a)
@compat Base.:*(a::Number, b::Vec) = b*a

@compat @generated function Base.:+{N}(a::Vec{N}, b::Vec{N})
    Expr(:tuple, [:(a[$i] + b[$i]) for i=1:N]...)
end


text:
using RAFF
using Printf
using Random

using Logging
using Base.CoreLogging

"""

    calc_ratio( model_str::String, np::Int, p::Int,
        sol::Vector{Float64}, ntests::Int=10,
        initguess::Vector{Float64}=nothing, maxms::Int=1,
        fp::IOStream=stdout)

Run a sequence of `ntests` tests for model `model_str` using a dataset
of `np` points and `p` trusted points (i. e. `np - p` outliers). The
points are generated by perturbations around *solution* `sol`. See
function [`generate_noisy_data`](@ref) for more details about the
generation of random points and outliers.

If provided, `initguess` is used as a starting points for
[`raff`](@ref), `maxms` is the number of random initial points used
and `fp` is the output.

"""
function calc_ratio( model_str::String, np::Int, p::Int,
    sol::Vector{Float64}, ntests::Int=10,
    initguess=nothing, maxms::Int=1,fp=stdout; cluster=nothing)

    # Set Logging
    global_logger(ConsoleLogger(stdout, Logging.Error))

    n, model, = RAFF.model_list[model_str]

    tmpsol = Vector{Float64}(undef, n)

    large_number = 179424673

    # Tests
    tot_tim = 0.0
    n_match = zeros(Int, np + 1)
    n_exact = 0
    n_out = 0
    # Number of correct outliers found
    n_cout = 0
    # Number of false positives found
    n_fp = 0

    for i = 1:ntests

        # Define seed for this run
        Random.seed!(large_number + i)

        if initguess == nothing

            x = zeros(Float64, n)
        else

            x = copy(initguess)

        end

        tmpsol .= sol

        data, tmpsol, tout = if cluster == nothing

            generate_noisy_data(model, n, np, p, tmpsol, (1.0, 30.0))

        else

            generate_clustered_noisy_data(model, n, np, p, tmpsol,
                (1.0, 30.0), cluster)

        end

        rsol, t, = @timed praff(model, data[:, 1:end - 1], n;
                                MAXMS=maxms, initguess=x, ftrusted=0.7)

        cnt = 0

        for k in rsol.outliers

            (k in tout) && (cnt += 1)

        end

        n_out += length(rsol.outliers)

        n_cout += cnt

        n_fp += length(rsol.outliers) - cnt

        n_match[cnt + 1] += 1

        (cnt == np - p) && (length(rsol.outliers) == np - p) && (n_exact += 1)

        tot_tim += t

    end

    @printf(fp, "%10s %5d %5d %10.3f %10.3f %10.3f %10.3f %10.2f %5d %5d %5d %10.3f\n",
            model_str, np, p, n_match[np - p + 1] /
            (1.0 * ntests), n_exact / (1.0 * ntests), n_cout / ntests,
            n_fp / ntests, n_out / ntests, n_match[1], n_match[2],
            n_match[3], tot_tim)

    return n_match

end

"""

    run_calc_ratio(filename="/tmp/table.txt")

Perform a sequence of tests for different models and a combination of
different number of outliers.

Saves the results in `filename`.

"""
function run_calc_ratio(filename="/tmp/table.txt")

    for (model_str, sol) in [ ("linear", [-200.0, 1000.0]), ("cubic", [0.5, -20.0, 300.0, 1000.0]),
                              ("expon", [5000.0, 4000.0, 0.2]),
                              ("logistic", [6000.0, -5000, -0.2, -3.7]) ]

        for (np, p) in [(10, 9), (10, 8), (100, 99), (100, 90)]

            for maxms in [1, 10, 100, 1000]

                open(filename, "a") do fp

                    calc_ratio(model_str, np, p, sol, 1000, nothing, maxms, fp);

                end

            end

        end

    end

end

"""

    run_calc_ratio_clustered(filename="/tmp/table.txt")

Perform a sequence of tests for different models for solving problems
with clustered outliers.

Saves the results in `filename`.

"""
function run_calc_ratio_clustered(filename="/tmp/table.txt")

    np = 100

    p  = 90

    maxms = 100

    for (model_str, sol) in [ ("linear", [-200.0, 1000.0]), ("cubic", [0.5, -20.0, 300.0, 1000.0]),
                              ("expon", [5000.0, 4000.0, 0.2]),
                              ("logistic", [6000.0, -5000, -0.2, -3.7]) ]

        open(filename, "a") do fp

            calc_ratio(model_str, np, p, sol, 1000, nothing, maxms, fp, cluster=(5.0, 10.0));

        end

    end

end


text:
module PlotExamples

using Plots
using Colors
using Compat

const DOCDIR = Pkg.dir("Plots") * "/docs"
const IMGDIR = Pkg.dir("Plots") * "/img"

"""
Holds all data needed for a documentation example... header, description, and plotting expression (Expr)
"""
type PlotExample
  header::@compat(AbstractString)
  desc::@compat(AbstractString)
  exprs::Vector{Expr}
end

# the examples we'll run for each
const examples = PlotExample[
  PlotExample("Lines",
              "A simple line plot of the columns.",
              [
                :(plot(Plots.fakedata(50,5), w=3))
              ]),
  PlotExample("Functions, adding data, and animations",
              "Plot multiple functions.  You can also put the function first, or use the form `plot(f, xmin, xmax)` where f is a Function or AbstractVector{Function}.\n\nGet series data: `x, y = plt[i]`.  Set series data: `plt[i] = (x,y)`. Add to the series with `push!`/`append!`.\n\nEasily build animations.  (`convert` or `ffmpeg` must be available to generate the animation.)  Use command `gif(anim, filename, fps=15)` to save the animation.",
              [
                :(p = plot([sin,cos], zeros(0), leg=false)),
                :(anim = Animation()),
                :(for x in linspace(0, 10π, 100)
                    push!(p, x, Float64[sin(x), cos(x)])
                    frame(anim)
                  end)
              ]),
  PlotExample("Parametric plots",
              "Plot function pair (x(u), y(u)).",
              [
                :(plot(sin, x->sin(2x), 0, 2π, line=4, leg=false, fill=(0,:orange)))
              ]),
  PlotExample("Colors",
              "Access predefined palettes (or build your own with the `colorscheme` method).  Line/marker colors are auto-generated from the plot's palette, unless overridden.  Set the `z` argument to turn on series gradients.",
              [
                :(y = rand(100)),
                :(plot(0:10:100,rand(11,4),lab="lines",w=3, palette=:grays, fill=(0.5,:auto))),
                :(scatter!(y, zcolor=abs(y-.5), m=(:heat,0.8,stroke(1,:green)), ms=10*abs(y-0.5)+4, lab="grad"))
              ]),
  PlotExample("Global",
              "Change the guides/background/limits/ticks.  Convenience args `xaxis` and `yaxis` allow you to pass a tuple or value which will be mapped to the relevant args automatically.  The `xaxis` below will be replaced with `xlabel` and `xlims` args automatically during the preprocessing step. You can also use shorthand functions: `title!`, `xaxis!`, `yaxis!`, `xlabel!`, `ylabel!`, `xlims!`, `ylims!`, `xticks!`, `yticks!`",
              [
                :(y = rand(20,3)),
                :(plot(y, xaxis=("XLABEL",(-5,30),0:2:20,:flip), background_color = RGB(0.2,0.2,0.2), leg=false)),
                :(hline!(mean(y,1)+rand(1,3), line=(4,:dash,0.6,[:lightgreen :green :darkgreen]))),
                :(vline!([5,10])),
                :(title!("TITLE")),
                :(yaxis!("YLABEL", :log10))
              ]),
  PlotExample("Two-axis",
              "Use the `axis` arguments.\n\nNote: Currently only supported with Qwt and PyPlot",
              [
                :(plot(Vector[randn(100), randn(100)*100], axis = [:l :r], ylabel="LEFT", yrightlabel="RIGHT", xlabel="X", title="TITLE"))
              ]),
  PlotExample("Arguments",
              "Plot multiple series with different numbers of points.  Mix arguments that apply to all series (marker/markersize) with arguments unique to each series (colors).  Special arguments `line`, `marker`, and `fill` will automatically figure out what arguments to set (for example, we are setting the `linestyle`, `linewidth`, and `color` arguments with `line`.)  Note that we pass a matrix of colors, and this applies the colors to each series.",
              [
                :(ys = Vector[rand(10), rand(20)]),
                :(plot(ys, line=(:dot,4,[:black :orange]), marker=([:hex :d],12,0.8,stroke(3,:gray))))
              ]),
  PlotExample("Build plot in pieces",
              "Start with a base plot...",
              [
                :(plot(rand(100)/3, reg=true, fill=(0,:green)))
              ]),
  PlotExample("",
              "and add to it later.",
              [
                :(scatter!(rand(100), markersize=6, c=:orange))
              ]),
  PlotExample("Heatmaps",
              "",
              [
                :(heatmap(randn(10000),randn(10000), nbins=20))
              ]),
  PlotExample("Line types",
              "",
              [
                :(types = intersect(supportedTypes(), [:line, :path, :steppre, :steppost, :sticks, :scatter])'),
                :(n = length(types)),
                :(x = Vector[sort(rand(20)) for i in 1:n]),
                :(y = rand(20,n)),
                :(plot(x, y, line=(types,3), lab=map(string,types), ms=15))
              ]),
  PlotExample("Line styles",
              "",
              [
                :(styles = setdiff(supportedStyles(), [:auto])'),
                :(plot(cumsum(randn(20,length(styles)),1), style=:auto, label=map(string,styles), w=5))
              ]),
  PlotExample("Marker types",
              "",
              [
                :(markers = setdiff(supportedMarkers(), [:none,:auto,Shape])'),
                :(n = length(markers)),
                :(x = linspace(0,10,n+2)[2:end-1]),
                :(y = repmat(reverse(x)', n, 1)),
                :(scatter(x, y, m=(8,:auto), lab=map(string,markers), bg=:linen))
              ]),
  PlotExample("Bar",
              "x is the midpoint of the bar. (todo: allow passing of edges instead of midpoints)",
              [
                :(bar(randn(999)))
              ]),
  PlotExample("Histogram",
              "",
              [
                :(histogram(randn(1000), nbins=20))
              ]),
  PlotExample("Subplots",
              """
                subplot and subplot! are distinct commands which create many plots and add series to them in a circular fashion.
                You can define the layout with keyword params... either set the number of plots `n` (and optionally number of rows `nr` or
                number of columns `nc`), or you can set the layout directly with `layout`.
              """,
              [
                :(subplot(randn(100,5), layout=[1,1,3], t=[:line :hist :scatter :step :bar], nbins=10, leg=false))
              ]),
  PlotExample("Adding to subplots",
              "Note here the automatic grid layout, as well as the order in which new series are added to the plots.",
              [
                :(subplot(Plots.fakedata(100,10), n=4, palette=[:grays :blues :heat :lightrainbow], bg=[:orange :pink :darkblue :black]))
              ]),
  PlotExample("",
              "",
              [
                :(subplot!(Plots.fakedata(100,10)))
              ]),
  PlotExample("Open/High/Low/Close",
              "Create an OHLC chart.  Pass in a vector of OHLC objects as your `y` argument.  Adjust the tick width with arg `markersize`.",
              [
                :(n=20),
                :(hgt=rand(n)+1),
                :(bot=randn(n)),
                :(openpct=rand(n)),
                :(closepct=rand(n)),
                :(y = [OHLC(openpct[i]*hgt[i]+bot[i], bot[i]+hgt[i], bot[i], closepct[i]*hgt[i]+bot[i]) for i in 1:n]),
                :(ohlc(y; markersize=8))
              ]),
  PlotExample("Annotations",
              "Currently only text annotations are supported.  Pass in a tuple or vector-of-tuples: (x,y,text).  `annotate!(ann)` is shorthand for `plot!(; annotation=ann)`",
              [
                :(y = rand(10)),
                :(plot(y, ann=(3,y[3],text("this is #3",:left)))),
                :(annotate!([(5,y[5],text("this is #5",16,:red,:center)),
                             (10,y[10],text("this is #10",:right,20,"courier"))]))
              ]),
  PlotExample("Custom Markers",
              "A `Plots.Shape` is a light wrapper around vertices of a polygon.  For supported backends, pass arbitrary polygons as the marker shapes.  Note: The center is (0,0) and the size is expected to be rougly the area of the unit circle.",
              [
                :(verts = [(-1.0,1.0),(-1.28,0.6),(-0.2,-1.4),(0.2,-1.4),(1.28,0.6),(1.0,1.0),
                           (-1.0,1.0),(-0.2,-0.6),(0.0,-0.2),(-0.4,0.6),(1.28,0.6),(0.2,-1.4),
                           (-0.2,-1.4),(0.6,0.2),(-0.2,0.2),(0.0,-0.2),(0.2,0.2),(-0.2,-0.6)])
                :(plot(0.1:0.2:0.9, 0.7rand(5)+0.15,
                       l=(3,:dash,:lightblue),
                       m=(Shape(verts),30,RGBA(0,0,0,0.2)),
                       bg=:pink, fg=:darkblue,
                       xlim = (0,1), ylim=(0,1), leg=false))
              ]),

  PlotExample("Contours",
              "",
              [
                :(x = 1:0.3:20),
                :(y = x),
                :(f(x,y) = sin(x)+cos(y)),
                :(contour(x, y, f, fill=true))
              ]),

  PlotExample("Pie",
              "",
              [
                :(x = ["Nerds", "Hackers", "Scientists"]),
                :(y = [0.4, 0.35, 0.25]),
                :(pie(x, y, title="The Julia Community", l=0.5))
              ]),

  PlotExample("3D",
              "",
              [
                :(n = 100),
                :(ts = linspace(0,8π,n)),
                :(x = ts .* map(cos,ts)),
                :(y = 0.1ts .* map(sin,ts)),
                :(z = 1:n),
                :(plot(x, y, z, zcolor=reverse(z), m=(10,0.8,:blues,stroke(0)), leg=false, w=5)),
                :(plot!(zeros(n),zeros(n),1:n, w=10))
              ])

]

function createStringOfMarkDownCodeValues(arr, prefix = "")
  string("`", prefix, join(sort(map(string, arr)), "`, `$prefix"), "`")
end
createStringOfMarkDownSymbols(arr) = isempty(arr) ? "" : createStringOfMarkDownCodeValues(arr, ":")

function generate_markdown(pkgname::Symbol)

  # set up the backend, and don't show the plots by default
  pkg = backend(pkgname)
  # default(:show, false)

  # mkdir if necessary
  try
    mkdir("$IMGDIR/$pkgname")
  end

  # open the markdown file
  md = open("$DOCDIR/$(pkgname)_examples.md", "w")

  write(md, "## Examples for backend: $pkgname\n\n")

  write(md, "### Initialize\n\n```julia\nusing Plots\n$(pkgname)()\n```\n\n")

  for (i,example) in enumerate(examples)

    try

      # we want to always produce consistent results
      srand(1234)

      # run the code
      map(eval, example.exprs)

      # # save the png
      # imgname = "$(pkgname)_example_$i.png"

      # NOTE: uncomment this to overwrite the images as well
      if i == 2
        imgname = "$(pkgname)_example_$i.gif"
        gif(anim, "$IMGDIR/$pkgname/$imgname", fps=15)
      else
        imgname = "$(pkgname)_example_$i.png"
        png("$IMGDIR/$pkgname/$imgname")
      end

      # write out the header, description, code block, and image link
      write(md, "### $(example.header)\n\n")
      write(md, "$(example.desc)\n\n")
      write(md, "```julia\n$(join(map(string, example.exprs), "\n"))\n```\n\n")
      write(md, "![](../img/$pkgname/$imgname)\n\n")

    catch ex
      # TODO: put error info into markdown?
      warn("Example $pkgname:$i failed with: $ex")
    end

    #
  end

  write(md, "- Supported arguments: $(createStringOfMarkDownCodeValues(supportedArgs(pkg)))\n")
  write(md, "- Supported values for axis: $(createStringOfMarkDownSymbols(supportedAxes(pkg)))\n")
  write(md, "- Supported values for linetype: $(createStringOfMarkDownSymbols(supportedTypes(pkg)))\n")
  write(md, "- Supported values for linestyle: $(createStringOfMarkDownSymbols(supportedStyles(pkg)))\n")
  write(md, "- Supported values for marker: $(createStringOfMarkDownSymbols(supportedMarkers(pkg)))\n")
  write(md, "- Is `subplot`/`subplot!` supported? $(subplotSupported(pkg) ? "Yes" : "No")\n\n")

  write(md, "(Automatically generated: $(now()))")
  close(md)

end

# make and display one plot
function test_examples(pkgname::Symbol, idx::Int; debug = true)
  Plots._debugMode.on = debug
  println("Testing plot: $pkgname:$idx:$(examples[idx].header)")
  backend(pkgname)
  backend()
  map(eval, examples[idx].exprs)
  plt = current()
  gui(plt)
  plt
end

# generate all plots and create a dict mapping idx --> plt
function test_examples(pkgname::Symbol; debug = false)
  Plots._debugMode.on = debug
  plts = Dict()
  for i in 1:length(examples)
    # if examples[i].header == "Subplots" && !subplotSupported()
    #   break
    # end

    try
      plt = test_examples(pkgname, i, debug=debug)
      plts[i] = plt
    catch ex
      # TODO: put error info into markdown?
      warn("Example $pkgname:$i:$(examples[i].header) failed with: $ex")
    end
  end
  plts
end

# axis            # :left or :right
# color           # can be a string ("red") or a symbol (:red) or a ColorsTypes.jl
#                 #   Colorant (RGB(1,0,0)) or :auto (which lets the package pick)
# label           # string or symbol, applies to that line, may go in a legend
# width           # width of a line
# linetype        # :line, :step, :stepinverted, :sticks, :scatter, :none, :heatmap, :hexbin, :hist, :bar
# linestyle       # :solid, :dash, :dot, :dashdot, :dashdotdot
# marker          # :none, :ellipse, :rect, :diamond, :utriangle, :dtriangle,
#                 #   :cross, :xcross, :star1, :star2, :hexagon
# markercolor     # same choices as `color`, or :match will set the color to be the same as `color`
# markersize      # size of the marker
# nbins           # number of bins for heatmap/hexbin and histograms
# heatmap_c       # color cutoffs for Qwt heatmaps
# fill          # fill value for area plots
# title           # string or symbol, title of the plot
# xlabel          # string or symbol, label on the bottom (x) axis
# ylabel          # string or symbol, label on the left (y) axis
# yrightlabel     # string or symbol, label on the right (y) axis
# reg             # true or false, add a regression line for each line
# size            # (Int,Int), resize the enclosing window
# pos             # (Int,Int), move the enclosing window to this position
# windowtitle     # string or symbol, set the title of the enclosing windowtitle
# screen          # Integer, move enclosing window to this screen number (for multiscreen desktops)

@compat const _ltdesc = Dict(
    :none       => "No line",
    :line       => "Lines with sorted x-axis",
    :path       => "Lines",
    :steppre    => "Step plot (vertical then horizontal)",
    :steppost   => "Step plot (horizontal then vertical)",
    :sticks     => "Vertical lines",
    :scatter    => "Points, no lines",
    :heatmap    => "Colored regions by density",
    :hexbin     => "Similar to heatmap",
    :hist       => "Histogram (doesn't use x)",
    :bar        => "Bar plot (centered on x values)",
    :hline      => "Horizontal line (doesn't use x)",
    :vline      => "Vertical line (doesn't use x)",
    :ohlc       => "Open/High/Low/Close chart (expects y is AbstractVector{Plots.OHLC})",
    :contour    => "Contour lines (uses z)",
    :path3d     => "3D path (uses z)",
    :scatter3d  => "3D scatter plot (uses z)",
  )

function buildReadme()
  readme = readall("$DOCDIR/readme_template.md")

  # build keyword arg table
  table = "Keyword | Default | Type | Aliases \n---- | ---- | ---- | ----\n"
  allseries = merge(Plots._seriesDefaults, @compat(Dict(:line=>nothing, :marker=>nothing, :fill=>nothing)))
  allplots = merge(Plots._plotDefaults, @compat(Dict(:xaxis=>nothing, :yaxis=>nothing)))
  alldefs = merge(allseries, allplots)
  for k in Plots.sortedkeys(alldefs)
  # for d in (Plots._seriesDefaults, Plots._plotDefaults)
  #   for k in Plots.sortedkeys(d)
    aliasstr = createStringOfMarkDownSymbols(aliases(Plots._keyAliases, k))
    table = string(table, "`:$k` | `$(alldefs[k])` | $(haskey(allseries,k) ? "Series" : "Plot") | $aliasstr  \n")
    # end
  end
  readme = replace(readme, "[[KEYWORD_ARGS_TABLE]]", table)

  # build linetypes table
  table = "Type | Desc | Aliases\n---- | ---- | ----\n"
  for lt in Plots._allTypes
    aliasstr = createStringOfMarkDownSymbols(aliases(Plots._typeAliases, lt))
    table = string(table, "`:$lt` | $(_ltdesc[lt]) | $aliasstr  \n")
  end
  readme = replace(readme, "[[LINETYPES_TABLE]]", table)

  # build linestyles table
  table = "Type | Aliases\n---- | ----\n"
  for s in Plots._allStyles
    aliasstr = createStringOfMarkDownSymbols(aliases(Plots._styleAliases, s))
    table = string(table, "`:$s` | $aliasstr  \n")
  end
  readme = replace(readme, "[[LINESTYLES_TABLE]]", table)

  # build markers table
  table = "Type | Aliases\n---- | ----\n"
  for s in Plots._allMarkers
    aliasstr = createStringOfMarkDownSymbols(aliases(Plots._markerAliases, s))
    table = string(table, "`:$s` | $aliasstr  \n")
  end
  readme = replace(readme, "[[MARKERS_TABLE]]", table)

  readme_fn = Pkg.dir("Plots") * "/README.md"
  f = open(readme_fn, "w")
  write(f, readme)
  close(f)

  gadfly()
  Plots.dumpSupportGraphs()
end

default(size=(500,300))

# run it!
# note: generate separately so it's easy to comment out
# @osx_only generate_markdown(:unicodeplots)
# generate_markdown(:qwt)
# generate_markdown(:gadfly)
# generate_markdown(:pyplot)
# generate_markdown(:immerse)
# generate_markdown(:winston)

end # module


text:
# Where are we?
_BASE_PATH = pwd()
_PATH_TO_DATA = joinpath(_BASE_PATH, "data")

# Package that we need -
using DataFrames
using CSV
using Distributions
using Statistics
using LinearAlgebra

# load the actual data -
path_to_experimental_data = joinpath(_PATH_TO_DATA, "Training-Thrombin-TF.csv")
full_df = CSV.read(path_to_experimental_data, DataFrame)

# what are the col names? (exclude subject id and visit id)
data_col_symbols = Symbol.(names(full_df)[3:end])
data_col_names_tuple = Tuple(Symbol.(names(full_df)[3:end]))
number_of_fields = length(data_col_symbols)

# what visit are we looking at?
# visit_df = filter(:visitid => x->(x==2 || x==3), full_df)

# sample -
number_of_synthetic_samples = 1000
data_array = Matrix(full_df[!, 3:end])
μ = mean(data_array,dims=1)
Σ = cov(data_array)
D = MvNormal(reshape(μ,16), Σ)
synthetic_sampled_data = rand(D, number_of_synthetic_samples)

# create a new data frame -
synthetic_data_frame = DataFrame()
for i ∈ 1:number_of_synthetic_samples

    # grab a col of data (a random sample)
    data_col = max.(0.0, synthetic_sampled_data[:,i])

    # build a data tuple to add to df -
    data_tuple = NamedTuple{data_col_names_tuple}(data_col)
    push!(synthetic_data_frame, data_tuple)
end

# setup scale factor dictionary -
SF = 1e9
scale_factor_dict = Dict()
scale_factor_dict[:II] = (1.4e-6)*SF
scale_factor_dict[:V] = (2e-8)*SF
scale_factor_dict[:VII] = (1e-8)*SF
scale_factor_dict[:VIII] = (7e-10)*SF
scale_factor_dict[:IX] = (9e-8)*SF
scale_factor_dict[:X] = (1.6e-7)*SF
scale_factor_dict[:XI] = (1e-8)*SF
scale_factor_dict[:XII] = (1e-8)*SF
scale_factor_dict[:AT] = (3.4e-6)*SF
scale_factor_dict[:PC] = (6.3e-8)*SF
scale_factor_dict[:TFPI] = 1.0

# ok, finally - let's convert these data into the correct units - TFPI is already in nM, the rest
# are in percentages of nominal -
for i ∈ 1:number_of_synthetic_samples

    for j ∈ 1:number_of_fields

        # look up the scale factor -
        fld_symbol = data_col_symbols[j]
        if (fld_symbol != :TFPI && haskey(scale_factor_dict,fld_symbol) == true) # TFPI is already in the correct units -
            SF_val = scale_factor_dict[fld_symbol]
            old_value = synthetic_data_frame[i,fld_symbol]
            new_value = SF_val*(old_value/100)
            synthetic_data_frame[i,fld_symbol] = new_value
        end
    end
end

# dump sample data to disk -
path_to_synthetic_data = joinpath(_PATH_TO_DATA, "Training-Synthetic-Thrombin-TF.csv")
CSV.write(path_to_synthetic_data, synthetic_data_frame)

max_stars_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_issues_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_forks_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation


text:
# Copyright 2020-2021 John T. Foster
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
using TaD
using Test
using SafeTestsets

@testset "TaD Tests" begin

    @safetestset "BSpline Tests" begin
        include("test_bspline.jl")
    end

end


text:
@deprecate ringGraph(n) ring_graph(n)
@deprecate pathGraph(n) path_graph(n)
@deprecate completeGraph(n) complete_graph(n)
@deprecate generalizedRing(n) generalized_ring(n)
@deprecate randGenRing(n::Int64, k::Integer; verbose=false) rand_gen_ring(n, k; verbose=verbose)
@deprecate hyperCube(n::Int64) hypercube(n)
@deprecate completeBinaryTree(n::Int64) complete_binary_tree(n)

@deprecate productGraph(a0::SparseMatrixCSC, a1::SparseMatrixCSC) product_graph(a0, a1)
#@deprecate wGrid2(n::Integer; weightGen::Function=rand) wgrid2(n, weightGen)
@deprecate randMatching(n::Integer) rand_matching(n)
@deprecate randRegular(n::Integer, k::Integer) rand_regular(n, k)
@deprecate grownGraph(n::Integer, k::Integer) grown_graph(n, k)
@deprecate grownGraphD(n::Integer, k::Integer) grown_graph_d(n, k)
@deprecate prefAttach(n::Integer, k::Integer, p::Float64) pref_attach(n, k, p)

@deprecate randWeight(a) rand_weight(a)
@deprecate randWeightSub(a) rand_weight_sub(a)
@deprecate wtedChimera(n::Integer, k::Integer; verbose=false) wted_chimera(n, k, verbose=verbose)
@deprecate wtedChimera(n::Integer; verbose=false) wted_chimera(n, verbose=verbose)
@deprecate semiWtedChimera(n::Integer; verbose=false, prefix="") semi_wted_chimera(n, verbose=verbose, prefix=prefix)
@deprecate semiWtedChimera(n::Integer, k::Integer; verbose=false, prefix="") semiwted_chimera(n, k, verbose=verbose, prefix=prefix)

@deprecate thickenOnce(a) thicken_once(a)

@deprecate twoLift(a, flip::AbstractArray{Bool,1}) two_lift(a, flip)
@deprecate twoLift(a) two_lift(a)
@deprecate twoLift(a, k::Integer) two_lift(a, k)

@deprecate plotGraph(gr,x,y;color=[0,0,1],dots=true,setaxis=true,number=false) plot_graph(gr,x,y;color=color,dots=dots,setaxis=setaxis,number=number)
@deprecate plotGraph(gr,x,y,z;color=[0,0,1],dots=true,setaxis=true,number=false) plot_graph(gr,x,y,z;color=color,dots=dots,setaxis=setaxis,number=number)

@deprecate approxCholLap(a; args...) approxchol_lap(a; args...)
@deprecate approxCholSDDM(a; args...) approxchol_sddm(a; args...)
@deprecate cholLap(a; args...) chol_lap(a; args...)
@deprecate cholSDDM(a; args...) chol_sddm(a; args...)

@deprecate spectralDrawing(a) spectral_drawing(a)
@deprecate spectralCoords(a) spectral_coords(a)
@deprecate compStretches(tree, mat) comp_stretches(tree, mat)

@deprecate readIJ(fn) read_graph(fn)
@deprecate readIJ(fn, sep) read_graph(fn)
@deprecate readIJV(fn) read_graph(fn)


text:
export encodevideo, encode!, prepareencoder, appendencode!, finishencode!, mux

# A struct collecting encoder objects for easy passing
mutable struct VideoEncoder
    codec_name::String
    apCodecContext::Array{Ptr{AVCodecContext}}
    apFrame::Array{Ptr{AVFrame}}
    apPacket::Array{Ptr{AVPacket}}
    pix_fmt::Integer
end

"""
    encode(encoder::VideoEncoder, io::IO)

Encode frame in memory
"""
function encode!(encoder::VideoEncoder, io::IO; flush=false)
    if flush
        ret = avcodec_send_frame(encoder.apCodecContext[1],C_NULL)
    else
        ret = avcodec_send_frame(encoder.apCodecContext[1],encoder.apFrame[1])
    end
    ret < 0 && error("Error $ret sending a frame for encoding")

    while ret >= 0
        ret = avcodec_receive_packet(encoder.apCodecContext[1], encoder.apPacket[1])
        if (ret == -35 || ret == -11 || ret == -541478725) # -35=EAGAIN or -11.  -541478725=AVERROR_EOF
             return
        elseif (ret < 0)
            error("Error $ret during encoding")
        end
        pkt = unsafe_load(encoder.apPacket[1])
        @debug println("Write packet $(pkt.pts) (size=$(pkt.size))")
        data = unsafe_wrap(Array,pkt.data,pkt.size)
        write(io,data)
        av_packet_unref(encoder.apPacket[1])
    end
end

const AVCodecContextPropertiesDefault =[:priv_data => ("crf"=>"22","preset"=>"medium")]

function islossless(prop)
    for p in prop
        if p[1] == :priv_data
            for subp in p[2]
                if (subp[1] == "crf") && (subp[2] == "0")
                    return true
                end
            end
        end
    end
    return false
end

isfullcolorrange(props) = (findfirst(map(x->x == Pair(:color_range,2),props)) != nothing)

"""
    prepareencoder(firstimg;framerate=30,AVCodecContextProperties=[:priv_data => ("crf"=>"22","preset"=>"medium")],codec_name::String="libx264")

Prepare encoder and return AV objects.
"""
function prepareencoder(firstimg; framerate=30, AVCodecContextProperties=[:priv_data => ("crf"=>"22","preset"=>"medium")], codec_name::String="libx264")
    height = size(firstimg,1)
    width = size(firstimg,2)
    ((width % 2 !=0) || (height % 2 !=0)) && error("Encoding error: Image dims must be a multiple of two")

    codec = avcodec_find_encoder_by_name(codec_name)
    codec == C_NULL && error("Codec '$codec_name' not found")

    apCodecContext = Ptr{AVCodecContext}[avcodec_alloc_context3(codec)]
    apCodecContext == [C_NULL] && error("Could not allocate video codec context")

    apPacket = Ptr{AVPacket}[av_packet_alloc()]
    apPacket == [C_NULL] && error("av_packet_alloc() error")

    if eltype(firstimg) == UInt8 && (codec_name in ["libx264", "h264_nvenc"])
        if islossless(AVCodecContextProperties)
            if !isfullcolorrange(AVCodecContextProperties)
                @warn """Encoding output not lossless.
                Lossless libx264 encoding of UInt8 requires
                :color_range=>2 within AVCodecContextProperties, to represent
                full 8-bit pixel value range."""
            end
        end
        pix_fmt = AV_PIX_FMT_GRAY8
    elseif eltype(firstimg) == Gray{N0f8} && (codec_name in ["libx264", "h264_nvenc"])
        if islossless(AVCodecContextProperties)
            if !isfullcolorrange(AVCodecContextProperties)
                @warn """Encoding output not lossless.
                Lossless libx264 encoding of Gray{N0f8} requires
                :color_range=>2 within AVCodecContextProperties, to represent
                full 8-bit pixel value range."""
            end
        end
        pix_fmt = AV_PIX_FMT_GRAY8
    elseif eltype(firstimg) == Gray{N0f8} && (codec_name == "libx264rgb")
        pix_fmt = AV_PIX_FMT_RGB24
    elseif eltype(firstimg) == RGB{N0f8} && (codec_name == "libx264rgb")
        if !islossless(AVCodecContextProperties)
            @warn """Codec libx264rgb has limited playback support. Given that
            selected encoding settings are not lossless (crf!=0), codec_name="libx264"
            will give better playback results"""
        end
        pix_fmt = AV_PIX_FMT_RGB24
    elseif eltype(firstimg) == RGB{N0f8} && (codec_name in ["libx264", "h264_nvenc"])
        if islossless(AVCodecContextProperties)
            @warn """Encoding output not lossless.
            libx264 does not support lossless RGB planes. RGB will be downsampled
            to lossy YUV420P. To encode lossless RGB use codec_name="libx264rgb" """
        end
        pix_fmt = AV_PIX_FMT_YUV420P
    else
        error("VideoIO: Encoding image element type $(eltype(firstimg)) with
        codec $codec_name not currently supported")
    end

    av_setfield(apCodecContext[1],:width,width)
    av_setfield(apCodecContext[1],:height,height)
    framerate_rat = Rational(framerate)
    av_setfield(apCodecContext[1],:time_base,AVRational(1/framerate_rat))
    av_setfield(apCodecContext[1],:framerate,AVRational(framerate_rat))
    av_setfield(apCodecContext[1],:pix_fmt,pix_fmt)

    codecContext = unsafe_load(apCodecContext[1])
    for prop in AVCodecContextProperties
        if prop[1] == :priv_data
            for pd in prop[2]
                av_opt_set(codecContext.priv_data, string(pd[1]), string(pd[2]), AV_OPT_SEARCH_CHILDREN)
            end
        else
            av_setfield(apCodecContext[1],prop[1],prop[2])
        end
    end

    # open it
    ret = avcodec_open2(apCodecContext[1], codec, C_NULL)
    ret < 0 && error("Could not open codec: Return code $(ret)")

    apFrame = Ptr{VideoIO.AVFrame}[VideoIO.av_frame_alloc()]
    apFrame == [C_NULL] && error("Could not allocate video frame")

    av_setfield(apFrame[1],:format,pix_fmt)
    av_setfield(apFrame[1],:width,width)
    av_setfield(apFrame[1],:height,height)

    ret = av_frame_get_buffer(apFrame[1], 32)
    ret < 0 && error("Could not allocate the video frame data")

    return VideoEncoder(codec_name, apCodecContext, apFrame, apPacket, pix_fmt)
end

"""
    appendencode(encoder::VideoEncoder, io::IO, img, index::Integer)

Send image object to ffmpeg encoder and encode
"""
function appendencode!(encoder::VideoEncoder, io::IO, img, index::Integer)

    flush(stdout)
    ret = av_frame_make_writable(encoder.apFrame[1])
    ret < 0 && error("av_frame_make_writable() error")

    frame = unsafe_load(encoder.apFrame[1]) #grab data from c memory
    img_eltype = eltype(img)

    if (img_eltype == UInt8) && (encoder.pix_fmt == AV_PIX_FMT_GRAY8)
        for r = 1:frame.height, c = 1:frame.width
            unsafe_store!(frame.data[1], img[r,c], ((r-1)*frame.linesize[1])+c)
        end
    elseif (img_eltype == Gray{N0f8}) && (encoder.pix_fmt == AV_PIX_FMT_GRAY8)
        img_uint8 = rawview(channelview(img))
        for r = 1:frame.height, c = 1:frame.width
            unsafe_store!(frame.data[1], img_uint8[r,c], ((r-1)*frame.linesize[1])+c)
        end
    elseif (img_eltype == Gray{N0f8}) && (encoder.pix_fmt == AV_PIX_FMT_RGB24)
        img_uint8 = repeat(PermutedDimsArray(rawview(channelview(img)),(2,1))[:],inner=3)
        unsafe_copyto!(frame.data[1], pointer(img_uint8), length(img_uint8))
    elseif (img_eltype  == RGB{N0f8}) && (encoder.pix_fmt == AV_PIX_FMT_RGB24)
        img_uint8 = PermutedDimsArray(rawview(channelview(img)),(1,3,2))[:]
        unsafe_copyto!(frame.data[1], pointer(img_uint8), length(img_uint8))
    elseif (img_eltype  == RGB{N0f8}) && (encoder.pix_fmt == AV_PIX_FMT_YUV420P)
        img_YCbCr = convert.(YCbCr{Float64}, img)
        img_YCbCr_half = convert.(YCbCr{Float64}, restrict(img))
        for r = 1:frame.height, c = 1:frame.width
            unsafe_store!(frame.data[1], round(UInt8,img_YCbCr[r,c].y), ((r-1)*frame.linesize[1])+c)
        end
        for r = 1:Int64(frame.height/2), c = 1:Int64(frame.width/2)
            unsafe_store!(frame.data[2], round(UInt8,img_YCbCr_half[r,c].cb), ((r-1)*frame.linesize[2])+c)
            unsafe_store!(frame.data[3], round(UInt8,img_YCbCr_half[r,c].cr), ((r-1)*frame.linesize[3])+c)
        end
    else
        error("Array element type not supported")
    end
    av_setfield(encoder.apFrame[1],:pts,index)
    encode!(encoder, io)
end

"""
    finishencode(encoder::VideoEncoder, io::IO)

End encoding by sending endencode package to ffmpeg, and close objects.
"""
function finishencode!(encoder::VideoEncoder, io::IO)
    encode!(encoder, io, flush=true) # flush the encoder
    write(io,UInt8[0, 0, 1, 0xb7]) # add sequence end code to have a real MPEG file
    avcodec_free_context(encoder.apCodecContext)
    av_frame_free(encoder.apFrame)
    av_packet_free(encoder.apPacket)
end

ffmpeg_framerate_string(fr::Real) = string(fr)
ffmpeg_framerate_string(fr::String) = fr
ffmpeg_framerate_string(fr::Rational) = "$(numerator(fr))/$(denominator(fr))"
ffmpeg_framerate_string(fr) = error("Framerate type not valid. Mux framerate should be a subtype of Real (Integer, Float64, Rational etc.), or String")

"""
    mux(srcfilename,destfilename,framerate;silent=false,deletestream=true)

Multiplex stream file into video container. Deletes stream file by default.
"""
function mux(srcfilename, destfilename, framerate; silent=false, deletestream=true)
    fr_str = ffmpeg_framerate_string(framerate)
    muxout = FFMPEG.exe(`-y -framerate $fr_str -i $srcfilename -c copy $destfilename`, collect=true)
    filter!(x->!occursin.("Timestamps are unset in a packet for stream 0.",x),muxout) #known non-bug issue with h264
    if occursin("ffmpeg version ",muxout[1]) && occursin("video:",muxout[end])
        deletestream && rm("$srcfilename")
        !silent && (@info "Video file saved: $destfilename")
        !silent && (@info muxout[end-1])
        !silent && (@info muxout[end])
    else
        @warn "Stream Muxing may have failed: $srcfilename into $destfilename"
        println.(muxout)
    end
end

"""
    encodevideo(filename::String,imgstack::Array;
        AVCodecContextProperties = AVCodecContextPropertiesDefault,
        codec_name = "libx264",
        framerate = 24)

Encode image stack to video file and return filepath.

"""
function encodevideo(filename::String,imgstack::Array;
    AVCodecContextProperties = AVCodecContextPropertiesDefault,
    codec_name = "libx264",
    framerate = 24,
    silent=false)

    io = Base.open("temp.stream","w")
    encoder = prepareencoder(imgstack[1],codec_name=codec_name,framerate=framerate,AVCodecContextProperties=AVCodecContextProperties)
    !silent && (p = Progress(length(imgstack), 1))   # minimum update interval: 1 second
    index = 1
    for img in imgstack
        appendencode!(encoder, io, img, index)
        !silent && next!(p)
        index += 1
    end
    finishencode!(encoder, io)
    close(io)
    mux("temp.stream",filename,framerate,silent=silent)
    return filename
end


text:
module MountainCarExperiment

using DeepRL
using Flux
using Random
using ProgressMeter
using TensorBoardLogger
using Logging
using Statistics

function construct_agent(s, num_actions)

    ϵ=0.1
    γ=1.0f0
    batch_size=32

    tn_update_freq=100
    update_freq=1
    min_mem_size=1000
    er_size = 10000
    hist = 1

    model = Chain(Dense(length(s)*hist, 64, Flux.relu),
                  Dense(64, 64, Flux.relu),
                  Dense(64, num_actions))

    # er = DeepRL.ExperienceReplayDef(er_size, 1, Int, DeepRL.StateBuffer{Float32}(er_size, length(s)))
    target_network = deepcopy(model)

    return DQNAgent(
        model,
        target_network,
        QLearning(γ, Flux.mse),
        DeepRL.RMSPropTFCentered(0.001),
        ϵGreedy(ϵ, num_actions),
        er_size,
        hist,
        s,
        batch_size,
        tn_update_freq,
        update_freq,
        min_mem_size,
        hist_squeeze = Val{true}()
    )
end

function main_experiment(seed, max_num_steps)

    mc = MountainCar(0.0, 0.0, true)
    Random.seed!(Random.GLOBAL_RNG, seed)
    rng = Random.GLOBAL_RNG

    s = MinimalRLCore.get_state(mc)

    agent = construct_agent(s, length(MinimalRLCore.get_actions(mc)))

    total_rews = Float32[]
    steps = Int[]

    max_episode_length = 5000

    front = ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇']
    p = ProgressMeter.Progress(
        max_num_steps;
        dt=0.1,
        desc="Step:",
        barglyphs=ProgressMeter.BarGlyphs('|','█',front,' ','|'))

    eps = 1
    while sum(steps) < max_num_steps
        cur_step = 0
        episode_cut_off = min(max_episode_length, max_num_steps - sum(steps))
        tr, stp =
            run_episode!(mc, agent, episode_cut_off, rng) do (s, a, s′, r)
                cur_step+=1
                # next!(p, showvalues=[(:step, sum(steps)+cur_step), (:episode, eps)])
                next!(p)
            end
        push!(total_rews, tr)
        push!(steps, stp)

        eps += 1
    end

    return agent.model, total_rews, steps

end

end


text:
const PyFraction_Type__ref = Ref(PyNULL)
PyFraction_Type() = begin
    ptr = PyFraction_Type__ref[]
    if isnull(ptr)
        m = PyImport_ImportModule("fractions")
        isnull(m) && return ptr
        ptr = PyObject_GetAttrString(m, "Fraction")
        Py_DecRef(m)
        isnull(m) && return ptr
        PyFraction_Type__ref[] = ptr
    end
    ptr
end

PyFraction_From(x::Union{Rational,Integer}) = begin
    t = PyFraction_Type()
    isnull(t) && return PyNULL
    a = PyTuple_New(2)
    isnull(a) && return PyNULL
    b = PyLong_From(numerator(x))
    isnull(b) && (Py_DecRef(a); return PyNULL)
    err = PyTuple_SetItem(a, 0, b)
    ism1(err) && (Py_DecRef(a); return PyNULL)
    b = PyLong_From(denominator(x))
    isnull(b) && (Py_DecRef(a); return PyNULL)
    err = PyTuple_SetItem(a, 1, b)
    ism1(err) && (Py_DecRef(a); return PyNULL)
    r = PyObject_CallObject(t, a)
    Py_DecRef(a)
    return r
end


text:
using Celeste

const read_with_strategy = function(strategy)
    rcf = Celeste.RunCamcolField(4263, 5, 119)
    Celeste.SDSSIO.load_field_images(strategy, [rcf])
    Celeste.SDSSIO.read_photoobj_files(strategy, [rcf])
end

function test_config(config)
    @show config
    strategy = Celeste.read_settings_file(joinpath("ioconfigs", config))
    read_with_strategy(strategy)
end

for config in ["plain.yml", "slurp.yml"]
    test_config(config)
end

rm("data/dr12", force = true, recursive=true)
cd("data") do
    # Create the original SDSS directory layout for testing purposes
    mkpath("dr12/boss/photo/redux/301/4263/objcs/5")
    mkpath("dr12/boss/photoObj/frames/301/4263/5")
    mkpath("dr12/boss/photoObj/301/4263/5")
    for band in ['z','g','r','i','u']
        cp("4263/5/119/frame-$band-004263-5-0119.fits", "dr12/boss/photoObj/frames/301/4263/5/frame-$band-004263-5-0119.fits")
        cp("4263/5/119/fpM-004263-$(band)5-0119.fit", "dr12/boss/photo/redux/301/4263/objcs/5/fpM-004263-$(band)5-0119.fit")
    end
    cp("4263/5/119/psField-004263-5-0119.fit", "dr12/boss/photo/redux/301/4263/objcs/5/psField-004263-5-0119.fit")
    cp("4263/5/119/photoObj-004263-5-0119.fits", "dr12/boss/photoObj/301/4263/5/photoObj-004263-5-0119.fits")
    cp("4263/5/photoField-004263-5.fits", "dr12/boss/photoObj/301/4263/photoField-004263-5.fits")
end

test_config("sdsslayout.yml")

rm("data/dr12", recursive=true)

worker = addprocs(1)[]
fetch(@spawnat worker Base.require(:Celeste))

# Test fetching over the network
# The basedir is relative, if the worker doesn't try to fetch via the network after this, it'll fail
fetch(@spawnat worker cd("/tmp"))
remotecall_fetch(read_with_strategy, worker, Celeste.SDSSIO.MasterRPCStrategy(Celeste.read_settings_file(joinpath("ioconfigs", "plain.yml"))))
rmprocs([worker])


text:
# """
#     (This method is not efficient for intence computations)
#     In Braunstein et al paper (see README)
#     KK = αF'F
#     Σ^-1 = (KK + D)
#     return Σ, covariance matrix of Q, the full multivariate Gaussian
# """
# function Q_sigma(model::MetNet, epout::EPOut, α)

#     # in epout sol d is scale using maxflux, we'll scale it back
#     maxflux = max(maximum(abs.(lb)), maximum(abs.(ub)))
#     a = epout.sol.a * maxflux;
#     d = epout.sol.b * maxflux^2;
#     #= a, d are the mean and variance of the univariate Gaussians
#     used as site approximations in the EP =#
#     S = model.S
#     D = Matrix(Diagonal(1.0 ./ d))
#     invΣ = α*S'*S + D
#     Σ = inv(invΣ)
# end


text:
"""
	grdpaste(cmd0::String="", arg1=nothing, arg2=nothing, kwargs...)

Combine grids ``grid1`` and ``grid2`` into ``grid3`` by pasting them together along their common edge.
Both grids must have the same dx, dy and have one edge in common.

Full option list at [`grdpaste`]($(GMTdoc)grdpaste.html)

Parameters
----------

- **G** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdpaste(....) form.
    ($(GMTdoc)grdpaste.html#g)
- $(GMT.opt_V)
- $(GMT.opt_f)
"""
function grdpaste(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:V_params :f :params])
	cmd = add_opt(cmd, 'G', d, [:G :outgrid])

	cmd, got_fname, arg1, arg2 = find_data(d, cmd0, cmd, arg1, arg2)
	return common_grd(d, "grdpaste " * cmd, arg1, arg2)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdpaste(arg1, arg2=nothing, cmd0::String=""; kw...) = grdpaste(cmd0, arg1, arg2; kw...)


text:
using ALifeBenchmark
import SharedArrays.SharedArray
using Plots
using Random
using Serialization
import Base.Threads.@threads
using StringDistances: Levenshtein
using Measures
using LightGraphs

function plot_result(times, values, name, trial_id)
    plot(times, values,
            title = "",
            label = "",
            xguide = "Time",
            yguide = name,
            seriestype = :scatter,
            markersize = 1.5,
            markerstrokewidth = 0,
            size = (900, 600),
            margin = 10mm,
            dpi = 1000)
    savefig("$name$trial_id.svg")
end

function plot_result(times_1, times_2, values, name, trial_id)
    plot(times_1, times_2,
            marker_z  = values,
            title = "",
            label = "",
            xguide = "Time",
            yguide = name,
            seriestype = :scatter,
            markersize = 3,
            markerstrokewidth = 0,
            size = (900, 600),
            margin = 10mm,
            dpi = 1000)
    savefig("$name$trial_id")
end

function save_result(data, name, trial_id)
    serialize("$name$trial_id", data)
end

function collect_basic_statistic(abbrevation, trial_id, get_statistic)
    collect_compare_to_end_statistic(abbrevation, trial_id, (x, _) -> get_statistic(x), identity)
end

function collect_compare_to_end_statistic(abbrevation, trial_id, get_statistic, prepare)
    @info "Start $abbrevation"

    snapshot_ids = get_snapshot_ids(trial_id)
    num_snapshots = length(snapshot_ids)

    prepared_last = last(snapshot_ids) |> get_snapshot |> prepare

    values = SharedArray{Float64}(num_snapshots)
    times = SharedArray{UInt64}(num_snapshots)

    done = Threads.Atomic{Int}(0)

    @threads for (i, snapshot_id) in (shuffle(snapshot_ids) |> enumerate |> unique)
        snapshot = get_snapshot(snapshot_id)
        prepared = prepare(snapshot)

        current_value = get_statistic(prepared, prepared_last)
        current_time = get_time(snapshot)

        values[i] = current_value
        times[i] = current_time

        Threads.atomic_add!(done, 1)
        @info "$abbrevation \t $(done[] / num_snapshots) \t $current_time \t $current_value"

        if done[] % 100 == 0
            # plot_result(Array(times), Array(values), abbrevation, trial_id)
            # save_result((Array(times), Array(values)), abbrevation, trial_id)
        end
    end

    return (Array(times), Array(values))
end

function collect_cross_statistic(abbrevation, trial_id, get_statistic, prepare)
    @info "Start $abbrevation"

    snapshot_ids = get_snapshot_ids(trial_id)
    num_snapshots = length(snapshot_ids)

    values = SharedArray{Float64}(num_snapshots^2)
    times_1 = SharedArray{UInt64}(num_snapshots^2)
    times_2 = SharedArray{UInt64}(num_snapshots^2)

    done = Threads.Atomic{Int}(0)

    @threads for (i, snapshot_id_1) in (shuffle(snapshot_ids[1:5:end]) |> enumerate |> unique)
        snapshot_1 = get_snapshot(snapshot_id_1)
        prepared_1 = prepare(snapshot_1)
        time_1 = get_time(snapshot_1)

        for (j, snapshot_id_2) in (shuffle(snapshot_ids[1:5:end]) |> enumerate |> unique)
            snapshot_2 = get_snapshot(snapshot_id_2)
            prepared_2 = prepare(snapshot_2)
            time_2 = get_time(snapshot_2)

            current_value = get_statistic(prepared_1, prepared_2)

            index = (i-1)*num_snapshots + j

            values[index] = current_value
            times_1[index] = time_1
            times_2[index] = time_2

            Threads.atomic_add!(done, 1)

            if done[] % 1000 == 0
                @info "$abbrevation \t $(done[] / num_snapshots / num_snapshots) \t $time_1 \t $time_2 \t $current_value"
            end

            if done[] % 100_000 == 0
                plot_result(Array(times_1), Array(times_2), Array(values), abbrevation, trial_id)
                save_result((Array(times_1), Array(times_2), Array(values)), abbrevation, trial_id)
            end
        end
    end

    return (Array(times_1), Array(times_2), Array(values))
end

"""
Build Cache
"""
phenotype_cache = Dict()
graph_data = load_graph_data()

for edge in edges(graph_data.phenotype_graph)
    u, v = src(edge), dst(edge)
    genotype_u, genotype_v = graph_data.genotype_vertex_mapping(u), graph_data.genotype_vertex_mapping(v)

    phenotype_similarity = edge.weight
    phenotype_cache[(genotype_u, genotype_v)] = phenotype_similarity
    phenotype_cache[(genotype_v, genotype_u)] = phenotype_similarity
end

"""
COLLECT STATISTIC
"""

if length(ARGS) == 2
    trial_id, type_of_analysis = ARGS
    threshold = 0.001
    min_samples, max_samples = 50, 500
elseif length(ARGS) == 5
    trial_id, type_of_analysis, threshold, min_samples, max_samples = ARGS
    threshold = parse(Float64, threshold)
    min_samples = parse(Int, min_samples)
    max_samples = parse(Int, max_samples)
else
    trial_id = "12433992799852588"
    type_of_analysis = "LA"
    threshold = 0.001
    min_samples, max_samples = 50, 500
end

@info ARGS

if type_of_analysis == "LA"
    result = collect_compare_to_end_statistic(
        type_of_analysis,
        trial_id,
        (snapshot, last_snapshot) -> get_adaption_of_snapshot(last_snapshot, snapshot, threshold, min_samples, max_samples),
        identity
    )

elseif type_of_analysis == "RF"
    result = collect_basic_statistic(type_of_analysis, trial_id, snapshot -> get_reachable_fitness(snapshot, threshold, min_samples, max_samples))

elseif type_of_analysis == "RD"
    result = collect_basic_statistic(type_of_analysis, trial_id, snapshot -> get_reachable_diversity(snapshot, threshold, min_samples, max_samples))

elseif type_of_analysis == "PD"
    result = collect_compare_to_end_statistic(
        type_of_analysis,
        trial_id,
        (distribution, last_distribution) -> _wasserstein(last_distribution, distribution, Levenshtein()),
        get_genotype_distribution
    )

elseif type_of_analysis == "EP"
    result = collect_basic_statistic(type_of_analysis, trial_id, snapshot -> get_evolutionary_potential(snapshot, 600_000, threshold, min_samples, max_samples))

elseif type_of_analysis == "CPD"
    result = collect_cross_statistic(
        type_of_analysis,
        trial_id,
        (distribution_1, distribution_2) -> _wasserstein(distribution_1, distribution_2, Levenshtein()),
        get_genotype_distribution
    )

elseif type_of_analysis == "CLA"
    result = collect_cross_statistic(
        type_of_analysis,
        trial_id,
        (snapshot_1, snapshot_2) -> get_adaption_of_snapshot(snapshot_1, snapshot_2, threshold, min_samples, max_samples),
        identity
    )

elseif type_of_analysis == "ENT"
    result = collect_basic_statistic(type_of_analysis, trial_id, get_entropy)

elseif type_of_analysis == "GD"
    result = collect_basic_statistic(type_of_analysis, trial_id, snapshot -> get_genotype_diversity(snapshot, Levenshtein(), threshold, min_samples, max_samples))

elseif type_of_analysis == "PHD"
    result = collect_basic_statistic(type_of_analysis, trial_id, snapshot -> get_phenotype_diversity(snapshot, threshold, min_samples, max_samples))

elseif type_of_analysis == "NN"
    result = collect_basic_statistic(type_of_analysis, trial_id, snapshot -> get_neutrality(snapshot, 1e-5, phenotype_cache))

elseif type_of_analysis == "NullNN"
    result = collect_basic_statistic(type_of_analysis, trial_id, snapshot -> get_neutrality_null_model(snapshot, 1e-5, phenotype_cache))

end

plot_result(result..., type_of_analysis, trial_id)
save_result(result, type_of_analysis, trial_id)

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
"""
Return the principal axes of the rigid body
defined by the input coordinates.
The order of the principals axes is sorted
by the eigenvalues of the moment of inertia matrix
(from low to high).
Therefore the first column of the returned
principal axes matrix is the direction
along which the object is the most elongated.

Arguments
-----------
:Type{PrincipalAxes}
    must be type `PrincipalAxes`

coordinates:AbstractArray
    must be a 2D matrix and must have exactly 3 columns

center=[]:AbstractArray
    (keyword) center of rotation
"""
function gage(::Type{PrincipalAxes},coordinates::AbstractArray; center::AbstractArray=[])
    if length(coordinates) == 0
        return Array{AbstractFloat,1}()
    end

    # If input coordinates is an array of array
    # then convert it to a 2D matrix
    if issubtype(typeof(coordinates[1]), Array)
        coordinates = vcat([x' for x in coordinates]...)
    end

    # change rotation center
    if length(center) == 3
        coordinates = coordinates .- transpose(center)
    end

    shape = size(coordinates)
    @debug @assert length(shape) == 2
    @debug @assert shape[2] == 3

    I = zeros(3,3) # moment of inertia tensor
    I[1,1] = sum(coordinates[:,2].^2 + coordinates[:,3].^2)
    I[2,2] = sum(coordinates[:,1].^2 + coordinates[:,3].^2)
    I[3,3] = sum(coordinates[:,1].^2 + coordinates[:,2].^2)
    I[1,2] = -sum(coordinates[:,1] .* coordinates[:,2])
    I[1,3] = -sum(coordinates[:,1] .* coordinates[:,3])
    I[2,3] = -sum(coordinates[:,2] .* coordinates[:,3])
    I[2,1] = I[1,2]
    I[3,1] = I[1,3]
    I[3,2] = I[2,3]
    D,V = eig(I)
    return V[:,sortperm(D)] # sort the eigenvectors by the eigenvalues
end


text:
function from_geneatlas(parsed_args)
    bridge = CSV.File(parsed_args["bridge"], header=["SAMPLE_ID", "geneatlas_id"]) |> DataFrame
    # Loading in pure datafame objects seems to be fine with 16GB of memory
    phenotypes  = CSV.File(parsed_args["phenotypes-file"], drop=[:FID]) |> DataFrame

    merged = innerjoin(bridge,
                       phenotypes,
                       on = :geneatlas_id => :IID)

    select!(merged, Not(:geneatlas_id))

    if ~(parsed_args["phenotypes-list"] isa Nothing)
        filter_list = readlines(open(parsed_args["phenotypes-list"]))
        phen_list = filter(x-> x ∈ filter_list, names(phenotypes))
        merged = merged[!, vcat(["SAMPLE_ID"], phen_list)]
    end

    if ~(parsed_args["withdrawal-list"] isa Nothing)
        withdrawn = CSV.File(parsed_args["withdrawal-list"], header=["SAMPLE_ID"]) |> DataFrame
        merged = filter(:SAMPLE_ID => ∉(withdrawn.SAMPLE_ID), merged)
    end

    CSV.write(parsed_args["output"], merged)
end

function prepare_phenotypes(parsed_args)
    from_geneatlas(parsed_args)
end

function write_phenotype_batch(phenotypes, batch, batch_start, batch_end)
    batch_file = batch_filename(batch)
    open(batch_file, "w") do io
        for i in batch_start:batch_end
            println(io, phenotypes[i])
        end
    end
end

batch_filename(batch) = string("phenotypes_batch_", batch, ".csv")

function tmle_phenotypes_batches(parsed_args)
    phenotypes = filter(
        !=("SAMPLE_ID"),
        split(readline(open(parsed_args["phenotypes-file"])), ",")
    )
    n_phenotypes = size(phenotypes, 1)
    batch_size = parsed_args["batch-size"]
    batch_size = batch_size == "max" ? size(phenotypes, 1) : parse(Int, batch_size)

    batch = 1
    batch_start = 1
    batch_end = batch_start + batch_size - 1
    while batch_end <= n_phenotypes
        write_phenotype_batch(phenotypes, batch, batch_start, batch_end)
        batch += 1
        batch_start = batch_end + 1
        batch_end = batch_start + batch_size - 1
    end

    if batch_start <= n_phenotypes
        write_phenotype_batch(phenotypes, batch, batch_start, n_phenotypes)
    end

end


text:
# Tests for Von-Mises Fisher distribution

using Distributions
using Base.Test

D = 3
mu = randn(D)
mu = mu / norm(mu)
kappa = 100.0
d = VonMisesFisher(mu, kappa)

# Basics

@test length(d) == D
@test d.kappa == kappa
@test_approx_eq d.mu mean(d)
@test_approx_eq norm(d.mu) 1.0

# MLE

x = rand(d, 10_000)
dmle = fit_mle(VonMisesFisher, x')
@test all(abs(mean(d) - mean(dmle)) .< .01)
@test_approx_eq norm(dmle.mu) 1.0
#@test abs(scale(dmle) - scale(d)) < .01 * scale(d) # within 1%? not always...

# Density

# TODO: Check against R's movMF. (Currently I'm a bit suspicious about their code.)
# > set.seed(1)
# > mu=c(1,0,0)
# > kappa=1.
# > x = rmovMF(1, mu, kappa)
# > x
#           [,1]       [,2]     [,3]
# [1,] 0.1772372 -0.4566632 0.871806
# > dmovMF(x, mu, kappa)
# [1] 1.015923
# WEIRD:
# > dmovMF(x, mu, 100)
# [1] 1.015923
# > dmovMF(x, mu, 1)
# [1] 1.015923

# mu = [1.0, 0., 0.]
# kappa = 1.0
# x = [0.1772372, -0.4566632, 0.871806]
# d = VonMisesFisher(mu, kappa)
#@test abs(logpdf(d, x) - 1.015923) < .00001


text:
using HorizonSideRobots
include("Functions.jl")

function draw_skew_cross(robot::Robot)
    sides_tuple = [(Nord, Ost), (Nord, West), (Sud, West), (Sud, Ost)]
    for sides in sides_tuple
        while (!isborder(robot, sides[1]) && !isborder(robot, sides[2]))
            move!(robot, sides[1])
            move!(robot, sides[2])
            putmarker!(robot)
        end

        rsides = reverse_side.(sides)

        while ismarker(robot)
            move!(robot, rsides[1])
            move!(robot, rsides[2])
        end
    end
    putmarker!(robot)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
################################################################################
#
#   HONEYCOMB LATTICE
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcellHoneycomb(
            unitcell_type  :: Type{U},
            implementation :: Val{1}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2}, U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # Bravais lattice vectors
        Vector{Float64}[
            Float64[sqrt(3.0)/2, -0.5],
            Float64[sqrt(3.0)/2, +0.5]
        ],
        # Sites
        S[
            newSite(S, Float64[0.0, 0.0],         getDefaultLabelN(LS,1)),
            newSite(S, Float64[1/sqrt(3.0), 0.0], getDefaultLabelN(LS,2))
        ],
        # Bonds
        B[
            newBond(B, 1, 2, getDefaultLabel(LB), (0, 0)),
            newBond(B, 1, 2, getDefaultLabel(LB), (-1, 0)),
            newBond(B, 1, 2, getDefaultLabel(LB), (0, -1)),
            newBond(B, 2, 1, getDefaultLabel(LB), (0, 0)),
            newBond(B, 2, 1, getDefaultLabel(LB), (1, 0)),
            newBond(B, 2, 1, getDefaultLabel(LB), (0, 1))
        ]
    )
end

# Implementation
# - implementation 4 (Kitaev)
# - labels <: Any
function getUnitcellHoneycomb(
            unitcell_type  :: Type{U},
            implementation :: Val{4}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2}, U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # Bravais lattice vectors
        Vector{Float64}[
            Float64[sqrt(3.0)/2, -0.5],
            Float64[sqrt(3.0)/2, +0.5]
        ],
        # Sites
        S[
            newSite(S, Float64[0.0, 0.0],         getDefaultLabelN(LS,1)),
            newSite(S, Float64[1/sqrt(3.0), 0.0], getDefaultLabelN(LS,2))
        ],
        # Bonds
        B[
            newBond(B, 1, 2, getDefaultLabelX(LB), (0, 0)),
            newBond(B, 1, 2, getDefaultLabelY(LB), (-1, 0)),
            newBond(B, 1, 2, getDefaultLabelZ(LB), (0, -1)),
            newBond(B, 2, 1, getDefaultLabelX(LB), (0, 0)),
            newBond(B, 2, 1, getDefaultLabelY(LB), (1, 0)),
            newBond(B, 2, 1, getDefaultLabelZ(LB), (0, 1))
        ]
    )
end