text:
abstract type GMC_Tuner end

mutable struct τ_Tuner <: GMC_Tuner
    τ::Float64
    τ_history::Vector{Float64}
    α::Float64
    β::Float64
    a::BitVector
    memory::Integer

    function τ_Tuner(τ1, e::GMC_NS_Ensemble)
        new(τ1,
        zeros(0),
        e.GMC_tune_α,
        e.GMC_tune_β,
        falses(0),
        e.GMC_tune_μ)
    end
end

mutable struct τ_PID <: GMC_Tuner
    τ::Float64
    τ_history::Vector{Float64}
    α::Float64
    β::Float64

    a::BitVector
    a_history::BitVector
    ℵ::Vector{Float64}
    memory::Integer
    Kp::Float64
    Ki::Float64
    Kd::Float64

    ϵ::Float64
    pterm::Float64
    iterm::Float64
    dterm::Float64

    function τ_PID(e::GMC_NS_Ensemble)
        new(e.GMC_init_τ,
        zeros(0),
        e.GMC_tune_α,
        1.,
        falses(0),
        falses(0),
        Vector{Float64}(),
        e.GMC_tune_μ,
        e.GMC_tune_PID...,
        0.,0.,0.,0.)
    end
end

function process_report!(t::τ_PID, report::Bool)
    push!(t.a_history, report)
    length(t.a_history)>t.memory ? (t.a=t.a_history[end-t.memory:end]) : (t.a=t.a_history)
    push!(t.ℵ,mean(t.a))
    τ_update!(t)
end

function τ_update!(t::τ_PID)
    push!(t.τ_history,t.τ)

    t.ϵ=t.ℵ[end] - t.α
    t.pterm=t.Kp * t.ϵ
    !(t.τ==eps()) && (t.iterm+=(t.Ki * t.ϵ))

    length(t.ℵ) > Int64(round(.1*t.memory)) ? (d = t.ℵ[end-Int64(round(.1*t.memory))] - t.ℵ[end]) : #d is positive if acceptance is declining- in this case want smaller beta and tau
        d=0.
    t.dterm=t.Kd * d

    t.β=max(eps(), 1. + t.pterm + t.iterm - t.dterm)
    t.τ=max(eps(), t.τ * t.β)
end

function init_tune(e::GMC_NS_Ensemble)
    @info "Performing initial GMC timestep tuning on ensemble..."
    τ=1.; tuned=false; τd=1.1; accept=1.; it=1
    while !tuned && it<e.GMC_tune_ι
        test_report=falses(0)
        println(stdout, "τ: $τ, It: $it, Accept: $accept")

        while length(test_report)<Int64(floor(length(e.models)*e.GMC_tune_ρ))
            m=deserialize(rand([m.path for m in e.models]))
            candidate=galilean_trajectory_sample!(m,e,τ)

            candidate.id==m.id ? push!(test_report,0) : push!(test_report,1)
            candidate=0
            GC.gc()
        end

        last_accept=accept
        accept=sum(test_report)/length(test_report)

        if accept == 1.
            τ*=τd
            last_accept == 1. && (τd*=1+e.GMC_tune_β)
            last_accept < e.GMC_tune_α && (τd=min(1.000001,1-e.GMC_tune_β*τd))
        elseif accept < e.GMC_tune_α
            τ/=τd
            last_accept == 1. && (τd=min(1.000001,1-e.GMC_tune_β*τd))
            last_accept < e.GMC_tune_α && (τd*=1+e.GMC_tune_β)
        else
            return τ
        end

        it+=1
        move_cursor_up_while_clearing_lines(stdout, 1)
    end
    println("Done. τ=$τ")
    return τ
end

function tune_τ!(t::τ_Tuner, step_report::BitVector)
    push!(t.τ_history,t.τ)
    t.a=vcat(t.a,step_report)
    ls=length(t.a)
    ls>t.memory && (t.a=t.a[ls-t.memory:end])

    accept=sum(t.a)/ls
    if accept==1.
        t.τ=t.τ*(1+t.β)
    elseif accept<t.α
        t.τ=t.τ*(1-t.β)
    end
end

function Base.show(io::IO, t::τ_PID; progress=false)
    try
        plot=lineplot(t.ℵ, title="Recent Unobstructed Moves", xlabel="Proposals", ylabel="Rate", color=:white, name="Free moves")

        lbls=plot.labels_left
        (arrow="")
        all(t.α .> [v for v in parse.(Float64,values(lbls))]) && (arrow="↑")
        all(t.α .< [v for v in parse.(Float64,values(lbls))]) && (arrow="↓")

        lineplot!(plot, [t.α for i in 1:length(t.ℵ)], name="α setpoint "*arrow, color=:red)
        show(io, plot)
        println()
        printstyled("τ:$(t.τ), α:$(t.α), β:$(t.β)",bold=true)
        println()
        printstyled("Kp:$(t.Kp), Ki:$(t.Ki), Kd:$(t.Kd)",color=:green)
        println()
        printstyled("ptrm:$(t.pterm), itrm:$(t.iterm), dtrm:$(t.dterm)",color=:magenta)
        println()

        progress && return nrows(plot.graphics)+9
    catch
        printstyled("τ_PID NOT AVAILABLE. STANDBY", bold=true)
        progress && return 1
    end
end

function get_tuner_dict(e)
    d=Dict{Int64,τ_PID}()
    for rec in e.models
        d[rec.trajectory]=τ_PID(e)
    end
    return d
end


text:
include(joinpath("Molenaar", "molenaar-basic.jl"))
import .MolenaarBasic

res = MolenaarBasic.molenaar_basic(2)

# Figures
results = []
for sext in [0.01:0.05:5;]
    push!(results, MolenaarBasic.molenaar_basic(sext))
end

# Test for the remake of the molenaar model
include(joinpath("Molenaar", "molenaar_remake2.jl"))
import .MolenaarRemake2

#res = MolenaarRemake2.molenaar_remake2(5.0)

results_test = []
for sext in [0.01:0.05:5;]
    push!(results_test, MolenaarRemake2.molenaar_remake2(sext))
end

#########

include(joinpath("Molenaar", "molenaar-Fig4A.jl"))
import .MolenaarAlt

res = MolenaarAlt.molenaar_alt(2)

# Figure 4A
stats = []
for sext in [0.01:0.05:5;]
    push!(stats, MolenaarAlt.molenaar_alt(sext))
end

#########

include(joinpath("Molenaar", "molenaar-Fig4B.jl"))
import .MolenaarAltB

res = MolenaarAltB.molenaar_altb(2)

# Figure 4B
statsb = []
for sext = 0.001:20
    push!(statsb, MolenaarAltB.molenaar_altb(sext))
end

ext: jl


text:
function theo_acf(phi)
    foo    = zeros(15)
    foo[1] = 1
        for  k in 2:15
           foo[k] = phi ^ k
        end
    plot(x=0:14, y = foo, Geom.line)
end

function theo_ar_process(phi)
    res    = zeros(100)
    res[1] = 0
    for i in 2:100
        res[i] = phi * res[i-1] + randn()
    end
    plot(x=1:100, y = res, Geom.line)
end


text:
"""
    function cv_milk(ts, cs, bv, G, A, ID)
---
Using data 2021-03-11 for cross validation again.
"""
function cv_milk(ts, vs, G, A)
    id = ts[:, :ix]
    yt = ts[:, :Milk]
    wt = ts[:, :Milk_ERC]
    nt = length(id)

    yv = vs[:, :Milk]
    append!(id, vs.ix)
    nv = length(yv)

    gi = inv(G[id, id])
    ai = inv(A[id, id])
    _, _ = compare(yt, yv, gi, ai)
    _, _ = compare(yt, wt, yv, gi, ai)
end

function cv_2021_03_11()
    @load "$dat_dir/jld/cv-setup.jld"     dcs gcs ncs
    @load "$dat_dir/jld/drp-training.jld" dts gts nts
    @load "$dat_dir/jld/ebv-all.jld"      dbv gbv nbv
    @load "$dat_dir/jld/GAID.jld"         G A ID
    ts, vs = mk_milk_dt(dts, dcs, dbv)
    cv_milk(ts, vs, G, A);
    ts, vs = mk_milk_dt(gts, gcs, gbv)
    cv_milk(ts, vs, G, A);
    ts, vs = mk_milk_dt(nts, ncs, nbv)
    println(size(ts), ' ', size(vs))
    cv_milk(ts, vs, G, A);
end

function cv_2020_07_02()
    @load "$dat_dir/jld/cv-setup.jld"     dcs gcs ncs
    @load "$dat_dir/jld/drp-training-2020-07-02.jld" dts gts nts
    @load "$dat_dir/jld/ebv-all.jld"      dbv gbv nbv
    @load "$dat_dir/jld/GAID.jld"         G A ID
    ts, vs = mk_milk_dt(dts, dcs, dbv)
    cv_milk(ts, vs, G, A);
    ts, vs = mk_milk_dt(gts, gcs, gbv)
    cv_milk(ts, vs, G, A);
    ts, vs = mk_milk_dt(nts, ncs, nbv)
    cv_milk(ts, vs, G, A);
end

"""
This function is important.
It found that I mistakenly included the missing data,
i.e., 999.0, for the calculation,
which led to very wrong results.
Below line was added after.
```julia
    dft = filter(row -> row.Milk<990, df1)
```
"""
function compare_data_2021_3()
    @load "$dat_dir/jld/cv-setup.jld"     dcs gcs ncs
    @load "$dat_dir/jld/ebv-all.jld"      dbv gbv nbv
    pre = pst = nothing
    #test 1
    begin
        @load "$dat_dir/jld/drp-training.jld" dts gts nts
        ts, vs = mk_milk_dt(dts, dcs, dbv)
        pre = copy(ts)
        println(first(ts, 10))
        println(first(vs, 10))
        println(describe(ts))
        println(describe(vs))
    end
    #test 2
    begin
        @load "$dat_dir/jld/drp-training-2020-07-02.jld" dts gts nts
        ts, vs = mk_milk_dt(dts, dcs, dbv)
        pst = copy(ts)
        println(first(ts, 10))
        println(first(vs, 10))
        println(describe(ts))
        println(describe(vs))
    end
    #test 3
    pre, pst
end


text:
export rrule
"""
    ChainRulesCore.rrule(itp::AbstractInterpolation, x...)
ChainRulesCore.jl `rrule` for integration with automatic differentiation libraries.
"""
function ChainRulesCore.rrule(itp::AbstractInterpolation, x...)
    y = itp(x...)
    function pullback(Δy)
        (ChainRulesCore.DoesNotExist(), Δy * Interpolations.gradient(itp, x...))
    end
    y, pullback
end


text:
using Gadfly

t = [0:0.2:2pi;]
plot(x=cos(t), y=sin(t), Geom.line(preserve_order=true))
plot(x=cos(t), y=sin(t), Geom.path)


text:
using Test
using PDFIO
using PDFIO.PD
using PDFIO.Cos
using PDFIO.Common

# Internal methods for testing only
using PDFIO.Cos: parse_indirect_ref, decode_ascii85, CosXString, parse_value

include("debugIO.jl")

@testset "PDFIO tests" begin
    @testset "Miscellaneous" begin
        @test_throws ErrorException skipv(IOBuffer([UInt8(65), UInt8(66)]),
                                          UInt8(66))
        @test CDTextString(CosXString([UInt8('0'), UInt8('0'),
                                       UInt8('4'),UInt8('1')]))=="A"
        @test CDTextString(CosXString([UInt8('4'), UInt8('2'),
                                       UInt8('4'),UInt8('1')]))=="BA"
        @test CosFloat(CosInt(1)) == CosFloat(1f0)
        @test [CosFloat(1f0), CosInt(2)] == [CosFloat(1f0), CosFloat(2f0)]
        @test CDRect(CosArray([CosInt(0),
                               CosInt(0),
                               CosInt(640),
                               CosInt(480)])) == CDRect(0, 0, 640, 480)
        @test parse_indirect_ref(IOBuffer(b"10 0 R\n")) ==
            CosIndirectObjectRef(10, 0)
        @test string(parse_value(IOBuffer("% This is a comment\r\n"))) ==
                     "% This is a comment"
    end

    @testset "CDDate" begin
        @test string(CDDate("D:199812231952-08'30 "))==
            "D:19981223195200-08'30"
        @test_throws ErrorException CDDate("not a date")
        @test_throws ErrorException CDDate("D:209")
        @test CDDate("D:2009") == CDDate("D:20090101000000Z")
        @test CDDate("D:200902") == CDDate("D:20090201000000+00")
        @test CDDate("D:20090202") == CDDate("D:20090202000000-00")
        @test CDDate("D:2009020201") == CDDate("D:20090202010000+00'00")
        @test CDDate("D:200902020102") == CDDate("D:20090202010200+00'00")
        @test CDDate("D:20090202010203") == CDDate("D:20090202010203+00'00")
        @test CDDate("D:20090202010203-00'01") < CDDate("D:20090202010202") < CDDate("D:20090202010203") < CDDate("D:20090202010203+00'01")
        @test CDDate("D:20090202+01'01") > CDDate("D:20090202+00'01") > CDDate("D:20090202-00'01") > CDDate("D:20090202-01'01")
        @test isless(CDDate("D:2009020208-06"), CDDate("D:2009020204-01"))
        @test isequal(CDDate("D:2009020208-06"), CDDate("D:2009020204-02"))
    end

    @testset "Test FlateDecode" begin
        @test begin
            filename="files/1.pdf"
            DEBUG && println(filename)
            doc = pdDocOpen(filename)
            DEBUG && println(pdDocGetCatalog(doc))
            cosDoc = pdDocGetCosDoc(doc)
            DEBUG && map(println, cosDoc.trailer)
            info = pdDocGetInfo(doc)
            @assert info["Producer"] == "LibreOffice 5.3" && info["Creator"] == "Writer"
            @assert pdDocGetPageCount(doc) == 2
            page = pdDocGetPage(doc, 1)
            @assert pdPageGetMediaBox(page) == pdPageGetCropBox(page)
            @assert pdPageIsEmpty(page) == false
            contents = pdPageGetContents(page)
            bufstm = get(contents)
            buf = read(bufstm)
            close(bufstm)
            @assert length(buf) == 18669
            @assert length(pdPageGetContentObjects(page).objs)==190
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end
    @testset "Document without Info" begin
        @test begin
            filename="files/1_noinfo.pdf"
            DEBUG && println(filename)
            doc = pdDocOpen(filename)
            DEBUG && println(pdDocGetCatalog(doc))
            cosDoc = pdDocGetCosDoc(doc)
            DEBUG && map(println, cosDoc.trailer)
            info = pdDocGetInfo(doc)
            @assert info === nothing
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Document with empty property" begin
        @test begin
            filename="files/empty_property.pdf"
            DEBUG && println(filename)
            doc = pdDocOpen(filename)
            DEBUG && println(pdDocGetCatalog(doc))
            cosDoc = pdDocGetCosDoc(doc)
            DEBUG && map(println, cosDoc.trailer)
            info = pdDocGetInfo(doc)
            @assert info == Dict(
                "Producer" => "Scribus PDF Library 1.3.3.13",
                "CreationDate" => CDDate("D:20090807192622"),
                "ModDate" => CDDate("D:20090807192622"),
                "Creator" => "Scribus 1.3.3.13",
                "Trapped" => cn"False")
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "LaTeX Document" begin
        @test begin
            filename="files/outline.pdf"
            DEBUG && println(filename)
            doc = pdDocOpen(filename)
            @assert doc !== nothing
            # Following two lines will be extended to more sensible test in the future
            # For now this is the test for issue #44 - without respective fix pdPageExtractText() crashes
            item_pg = pdDocGetPage(doc, 1)
            pdPageExtractText(IOBuffer(), item_pg)
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "PDF File with ObjectStreams" begin
        @test begin
            filename="files/pdf-17.pdf"
            DEBUG && println(filename)
            doc = pdDocOpen(filename)
            @assert pdDocGetPageCount(doc) == 1
            page = pdDocGetPage(doc, 1)
            rm = pdPageGetMediaBox(page)
            rc = pdPageGetCropBox(page)
            @assert rm == CDRect(0, 0, 612, 792)
            @assert rc == CDRect(0, 0, 612, 792)
            @assert pdPageIsEmpty(page) == false
            contents = pdPageGetContents(page)
            bufstm = get(contents)
            buf = read(bufstm)
            close(bufstm)
            @assert length(buf) == 1021
            @assert length(pdPageGetContentObjects(page).objs)==1
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "General File Opening 3" begin
        @test begin
            filename="3.pdf"
            DEBUG && println(filename)
            isfile(filename)||
                download("http://www.stillhq.com/pdfdb/000003/data.pdf",filename)
            doc = pdDocOpen(filename)
            @assert pdDocGetPageCount(doc) == 30
            page = pdDocGetPage(doc, 1)
            @assert pdPageIsEmpty(page) == false
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Hybrid x-ref" begin
        @test begin
            filename="A1947-15.pdf"
            DEBUG && println(filename)
            resfile, template, filename = local_testfiles(filename)
            doc = pdDocOpen(filename)
            io = util_open(resfile, "w")
            try
                extract_text(io, doc)
            finally
                util_close(io)
                pdDocClose(doc)
            end
            @assert files_equal(resfile, template)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Non-standard CMap" begin
        @test begin
            filename="16-969_o7jp.pdf"
            DEBUG && println(filename)
            resfile, template, filename = local_testfiles(filename)
            doc = pdDocOpen(filename)
            io = util_open(resfile, "w")
            try
                extract_text(io, doc)
            finally
                util_close(io)
                pdDocClose(doc)
            end
            @assert files_equal(resfile, template)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Corrupt File" begin
        @test begin
            filename="files/A1947-14.pdf"
            DEBUG && println(filename)
            doc = pdDocOpen(filename)
            try
                npage= pdDocGetPageCount(doc)
                for i=1:npage
                    page = pdDocGetPage(doc, i)
                    if pdPageIsEmpty(page) == false
                        pdPageGetContentObjects(page)
                        pdPageExtractText(IOBuffer(), page)
                    end
                end
            finally
                pdDocClose(doc)
            end
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Test RunLengthDecode" begin
        @test begin
            filename="582.pdf"
            DEBUG && println(filename)
            isfile(filename)||
                download("http://www.stillhq.com/pdfdb/000582/data.pdf",filename)
            doc = pdDocOpen(filename)
            info = pdDocGetInfo(doc)
            @assert info["Trapped"] == cn"False"
            @assert pdDocGetPageCount(doc) == 12
            obj=PDFIO.Cos.cosDocGetObject(doc.cosDoc,
                                          PDFIO.Cos.CosIndirectObjectRef(177, 0))
            stm=get(obj)
            data=read(stm)
            close(stm)
            @assert length(data) == 273
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Test ASCIIHexDecode" begin
        @test begin
        filename="325.pdf"
        DEBUG && println(filename)
        isfile(filename)||
            download("http://www.stillhq.com/pdfdb/000325/data.pdf",filename)
        doc = pdDocOpen(filename)
        @assert pdDocGetPageCount(doc) == 1
        obj = cosDocGetObject(doc.cosDoc, CosIndirectObjectRef(7, 0))
        stm=get(obj)
        data=read(stm)
        close(stm)
        @assert length(data) == 121203
        pdDocClose(doc)
        length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Test ASCII85Decode" begin
        @test take!(decode_ascii85(IOBuffer("zzz!!!~>"))) == fill(0x0, 14)
        @test begin
        filename="388.pdf"
        DEBUG && println(filename)
        isfile(filename)||
            download("http://www.stillhq.com/pdfdb/000388/data.pdf",filename)
        doc = pdDocOpen(filename)
        @assert pdDocGetPageCount(doc) == 1
        obj=PDFIO.Cos.cosDocGetObject(doc.cosDoc, PDFIO.Cos.CosIndirectObjectRef(9, 0))
        stm=get(obj)
        data=read(stm)
        close(stm)
        @assert length(data) == 38117
        pdDocClose(doc)
        length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "LZWDecode Filter" begin
        @test begin
            filename="589.pdf"
            DEBUG && println(filename)
            isfile(filename)||
                download("http://www.stillhq.com/pdfdb/000589/data.pdf",filename)
            doc = pdDocOpen(filename)
            obj=PDFIO.Cos.cosDocGetObject(doc.cosDoc, PDFIO.Cos.CosIndirectObjectRef(70, 0))
            stm=get(obj)
            data=read(stm)
            close(stm)
            @assert length(data) == 768
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
        @test begin
            filename="339.pdf"
            DEBUG && println(filename)
            isfile(filename)||
                download("http://www.stillhq.com/pdfdb/000339/data.pdf",filename)
            doc = pdDocOpen(filename)
            stm = get(cosDocGetObject(doc.cosDoc, CosIndirectObjectRef(4, 0)))
            buf = read(stm)
            util_close(stm)
            @assert length(buf) == 6636
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Content Array" begin
        @test begin
            filename="504.pdf"
            DEBUG && println(filename)
            isfile(filename)||
                download("http://www.stillhq.com/pdfdb/000504/data.pdf",filename)
            doc = pdDocOpen(filename)
            page = pdDocGetPage(doc, 1)
            contents = pdPageGetContents(page)
            stm = get(contents)
            buf = read(stm)
            util_close(stm)
            @assert length(buf) == 5574
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Test read_string" begin
        @test begin
            DEBUG && PDFIO.Cos.parse_data("files/page5.txt")
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Page label test" begin
        @test begin
            filename="431.pdf"
            DEBUG && println(filename)
            isfile(filename) ||
                download("http://www.stillhq.com/pdfdb/000431/data.pdf",filename)
            doc = pdDocOpen(filename)
            @assert pdDocGetPageCount(doc) == 54
            @assert pdDocHasPageLabels(doc)
            @assert PDFIO.Cos.cosDocGetPageNumbers(doc.cosDoc, doc.catalog, "title") ==
                range(1, length=1)
            @assert PDFIO.Cos.cosDocGetPageNumbers(doc.cosDoc, doc.catalog, "ii") ==
                range(3, length=1)
            @assert PDFIO.Cos.cosDocGetPageNumbers(doc.cosDoc, doc.catalog, "42") ==
                range(46, length=1)
            pdDocGetPageRange(doc, "iii")
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
        doc = pdDocOpen("files/1.pdf")
        @test pdDocHasPageLabels(doc) == false
        @test_throws ErrorException(E_INVALID_PAGE_LABEL) pdDocGetPageRange(doc, "1")
        pdDocClose(doc)
    end

    @testset "Symbol Fonts test" begin
        @test begin
            filename="431.pdf"
            result, template_file = local_testfiles(filename)
            DEBUG && println(filename)
            isfile(filename) ||
                download("http://www.stillhq.com/pdfdb/000431/data.pdf",filename)
            doc = pdDocOpen(filename)
            (npage = pdDocGetPageCount(doc)) == 54
            try
                open(result, "w") do io
                    for i=1:npage
                        page = pdDocGetPage(doc, i)
                        if pdPageIsEmpty(page) == false
                            pdPageGetContentObjects(page)
                            pdPageExtractText(io, page)
                        end
                    end
                end
                @assert files_equal(result, template_file)
            finally
                pdDocClose(doc)
            end
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Font Flags test" begin
        @test begin
            filename="files/pdf-sample.pdf"
            doc = pdDocOpen(filename)
            (npage = pdDocGetPageCount(doc)) == 1
            page = pdDocGetPage(doc, 1)
            fonts = pdPageGetFonts(page)
            d = Dict{CosName, Tuple}()
            for (k, v) in fonts
                d[k] = (pdFontIsAllCap(v), pdFontIsBold(v),
                        pdFontIsFixedW(v), pdFontIsItalic(v),
                        pdFontIsSmallCap(v))
            end
            @assert d[cn"F2"]  == (false, false, false, false, false)
            @assert d[cn"TT2"] == (false, true,  false, false, false)
            @assert d[cn"TT8"] == (false, true,  false, true,  false)
            @assert d[cn"TT6"] == (false, false, false, false, false)
            @assert d[cn"TT4"] == (false, false, false, false, false)
            pdDocClose(doc)
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Forms XObjects Test" begin
        @test begin
            filename="Graphics-wpf.pdf"
            result, template_file = testfiles(filename)
            DEBUG && println(filename)
            isfile(filename) ||
                download("http://www.pdfsharp.net/wiki/GetFile.aspx?"*
                         "File=%2fGraphics-sample%2fGraphics-wpf.pdf",
                         filename)
            doc = pdDocOpen(filename)
            @assert (npage = pdDocGetPageCount(doc)) == 5
            try
                open(result, "w") do io
                    for i=npage:npage
                        page = pdDocGetPage(doc, i)
                        if pdPageIsEmpty(page) == false
                            pdPageGetContentObjects(page)
                            pdPageExtractText(io, page)
                        end
                    end
                end
                @assert files_equal(result, template_file)
            finally
                 pdDocClose(doc)
            end
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Inline Image test" begin
        @test begin
            filename="Pratham-Sanskaran.pdf"
            result, template_file, src = local_testfiles(filename)
            DEBUG && println(src)
            doc = pdDocOpen(src)
            (npage = pdDocGetPageCount(doc)) == 3
            try
                open(result, "w") do io
                    for i=1:npage
                        page = pdDocGetPage(doc, i)
                        if pdPageIsEmpty(page) == false
                            pdPageGetContentObjects(page)
                            pdPageExtractText(io, page)
                        end
                    end
                end
                @assert files_equal(result, template_file)
            finally
                 pdDocClose(doc)
            end
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "MacRomanEncoding Fonts test" begin
        @test begin
            filename="spec-2.pdf"
            result, template_file, src = local_testfiles(filename)
            DEBUG && println(src)
            doc = pdDocOpen(src)
            @assert (npage = pdDocGetPageCount(doc)) == 1
            try
                open(result, "w") do io
                    for i=1:npage
                        page = pdDocGetPage(doc, i)
                        if pdPageIsEmpty(page) == false
                            pdPageGetContentObjects(page)
                            pdPageExtractText(io, page)
                        end
                    end
                end
                @assert files_equal(result, template_file)
            finally
                pdDocClose(doc)
            end
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Text before header test" begin
        @test begin
            filename="spec-2c.pdf"
            result, template_file, src = local_testfiles(filename)
            DEBUG && println(src)
            doc = pdDocOpen(src)
            @assert (npage = pdDocGetPageCount(doc)) == 1
            try
                open(result, "w") do io
                    for i=1:npage
                        page = pdDocGetPage(doc, i)
                        if pdPageIsEmpty(page) == false
                            pdPageGetContentObjects(page)
                            pdPageExtractText(io, page)
                        end
                    end
                end
                @assert files_equal(result, template_file)
            finally
                pdDocClose(doc)
            end
            length(utilPrintOpenFiles()) == 0
        end
    end

    @testset "Attachment PDF" begin
        filename = "fileAttachment.pdf"
        result, template_file, src = local_testfiles(filename)
        DEBUG && println(src)
        doc = pdDocOpen(src)
        val = string(pdDocGetNamesDict(doc))
        @test val == "\n33 0 obj\n<<\n\t/EmbeddedFiles\t34 0 R\n\t/JavaScript\t35 0 R\n>>\nendobj\n\n" ||
            val == "\n33 0 obj\n<<\n\t/JavaScript\t35 0 R\n\t/EmbeddedFiles\t34 0 R\n>>\nendobj\n\n"
        pdDocClose(doc)
        @test length(utilPrintOpenFiles()) == 0
    end
    files=readdir(get_tempdir())
    @assert length(files) == 0
end

if isfile("pvt/pvttests.jl")
    include("pvt/pvttests.jl")
end


text:
using SumOfSquares
using DynamicPolynomials
using Test

function sos_lower_bound(factory, sparsity::Sparsity.Pattern)
    @polyvar x y z
    K = @set 1 - x^2 - y^2 >= 0 && 1 - y^2 - z^2 >= 0
    p = 1 - x*y - y*z
    model = Model(factory)
    @variable(model, t)
    @objective(model, Max, t)
    @constraint(model, p - t in SOSCone(), domain = K, sparsity = sparsity)
    optimize!(model)
    return objective_value(model)
end

function sparse_vs_dense(optimizer_constructor)
    # use chordal sparsity
    sparse_value = sos_lower_bound(optimizer_constructor, Sparsity.Variable())

    # no chordal sparsity
    dense_value = sos_lower_bound(optimizer_constructor, Sparsity.NoPattern())

    @test abs(sparse_value - dense_value) < 1e-6
end

import CSDP
sparse_vs_dense(CSDP.Optimizer)

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
"""
   Kernels for GPs

Author: Eric Ford
Adapted from: https://github.com/eford/RvSpectraKitLearn.jl/blob/master/src/deriv_spectra_gp.jl

Could update to use KernelFunctions.jl, etc.
If upgrade to make use of other packages, please keep this so we ahve somet that gets the job done with minimal dependancies.
Author: Eric Ford
"""

function kernel_matern32(d::T; rho::T=one(T), sigmasq::T =one(T)) where { T<:Real }
  x = abs(sqrt(3)*d/rho)
  sigmasq * (1+x) * exp(-x)
end

function dkerneldx_matern32(d::T; rho::T=one(T), sigmasq::T =one(T)) where { T<:Real }
  x = abs(sqrt(3)*d/rho)
  -sigmasq * (3*d/rho^2) * exp(-x)
end

function d2kerneldx2_matern32(d::T; rho::T=one(T), sigmasq::T =one(T)) where { T<:Real }
  x = abs(sqrt(3)*d/rho)
  -sigmasq * (3/rho^2) * (1-x)*exp(-x)
end

function kernel_matern52(d::T; rho::T=one(T), sigmasq::T =one(T)) where { T<:Real }
  x = abs(sqrt(5)*d/rho)
  sigmasq * (1+x*(1+x/3)) * exp(-x)
end

function dkerneldx_matern52(d::T; rho::T=one(T), sigmasq::T =one(T)) where { T<:Real }
  x = abs(sqrt(5)*d/rho)
  -sigmasq/3 * (5*d/rho^2) *(1+x)*exp(-x)
end

function d2kerneldx2_matern52(d::T; rho::T=one(T), sigmasq::T =one(T)) where { T<:Real }
  x = abs(sqrt(5)*d/rho)
  -sigmasq/3 * (5/rho^2) * (1+x*(1-x))*exp(-x)
end

function nuttall_kernel(x::T; rho::T=one(T)) where { T<:Real }
      a0 = 0.355768
      a1 = 0.487396
      a2 = 0.144232
      a3 = 0.012604
   if abs(x)>rho return zero(T) end
   cpx = cospi(x/rho)
   spx = sinpi(x/rho)
   c2px = cpx*cpx-spx*spx
   return a0+a1*cpx+(a2+a3*cpx)*cpx*cpx-(a2+3*a3*cpx)*spx*spx
end

function nuttall_dkerneldx(x::T; rho::T=one(T)) where { T<:Real }
      a0 = 0.355768
      a1 = 0.487396
      a2 = 0.144232
      a3 = 0.012604
   if abs(x)>rho return zero(T) end
   cpx = cospi(x/rho)
   spx = sinpi(x/rho)
   c2px = cpx*cpx-spx*spx
   s2px = 2*cpx*spx
   s3px = spx*c2px+cpx*s2px
   return -pi/rho*(a1*spx+2*a2*s2px+3*a3*s3px)
end

function nuttall_d2kerneldx2(x::T; rho::T=one(T)) where { T<:Real }
      a0 = 0.355768
      a1 = 0.487396
      a2 = 0.144232
      a3 = 0.012604
   if abs(x)>rho return zero(T) end
   cpx = cospi(x/rho)
   spx = sinpi(x/rho)
   c2px = cpx*cpx-spx*spx
   c3px = cpx*(cpx*cpx-3*spx*spx)
   return -(pi/rho)^2*(a1*cpx+4*a2*c2px+9*a3*c3px)
end

function matern32_sparse_kernel(x::T; rho::T=one(T), bandwidth::T=4*one(T), sigmasq::T =one(T)) where { T<:Real }
  km = kernel_matern32(x,rho=rho,sigmasq=sigmasq)
  kn = nuttall_kernel(x,rho=rho*bandwidth)
  km*kn
end

function dkerneldx_matern32_sparse(x::T; rho::T=one(T), bandwidth::T=4*one(T), sigmasq::T =one(T)) where { T<:Real }
  km = kernel_matern32(x,rho=rho,sigmasq=sigmasq)
  kn = nuttall_kernel(x,rho=rho*bandwidth)
  dkmdx = dkerneldx_matern32(x,rho=rho,sigmasq=sigmasq)
  dkndx = nuttall_dkerneldx(x,rho=rho*bandwidth)
  km*dkndx+kn*dkmdx
end

function d2kerneldx2_matern32_sparse(x::T; rho::T=one(T), bandwidth::T=4*one(T), sigmasq::T =one(T)) where { T<:Real }
  km = kernel_matern32(x,rho=rho,sigmasq=sigmasq)
  kn = nuttall_kernel(x,rho=rho*bandwidth)
  dkmdx = dkerneldx_matern32(x,rho=rho,sigmasq=sigmasq)
  dkndx = nuttall_dkerneldx(x,rho=rho*bandwidth)
  d2kmdx2 = d2kerneldx2_matern32(x,rho=rho,sigmasq=sigmasq)
  d2kndx2 = nuttall_d2kerneldx2(x,rho=rho*bandwidth)
  km*d2kndx2+2*dkndx*dkmdx+kn*d2kmdx2
end

function matern52_sparse_kernel(x::T; rho::T=one(T), bandwidth::T=4*one(T), sigmasq::T =one(T)) where { T<:Real }
  km = kernel_matern52(x,rho=rho,sigmasq=sigmasq)
  kn = nuttall_kernel(x,rho=rho*bandwidth)
  km*kn
end

function dkerneldx_matern52_sparse(x::T; rho::T=one(T), bandwidth::T=4*one(T), sigmasq::T =one(T)) where { T<:Real }
  km = kernel_matern52(x,rho=rho,sigmasq=sigmasq)
  kn = nuttall_kernel(x,rho=rho*bandwidth)
  dkmdx = dkerneldx_matern52(x,rho=rho,sigmasq=sigmasq)
  dkndx = nuttall_dkerneldx(x,rho=rho*bandwidth)
  km*dkndx+kn*dkmdx
end

function d2kerneldx2_matern52_sparse(x::T; rho::T=one(T), bandwidth::T=4*one(T), sigmasq::T =one(T)) where { T<:Real }
  km = kernel_matern52(x,rho=rho,sigmasq=sigmasq)
  kn = nuttall_kernel(x,rho=rho*bandwidth)
  dkmdx = dkerneldx_matern52(x,rho=rho,sigmasq=sigmasq)
  dkndx = nuttall_dkerneldx(x,rho=rho*bandwidth)
  d2kmdx2 = d2kerneldx2_matern52(x,rho=rho,sigmasq=sigmasq)
  d2kndx2 = nuttall_d2kerneldx2(x,rho=rho*bandwidth)
  km*d2kndx2+2*dkndx*dkmdx+kn*d2kmdx2
end

#=
function make_kernel_matern32_data(x::AbstractArray{T,1};
			sigmasq_obs::AbstractArray{T,1} = zeros(length(x)),
			sigmasq_cor::T = 1.0, rho::T=one(T))  where { T<:Real }
  make_kernel_data(x, kernel_matern32, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
end

function make_kernel_matern32_obs_pred(xobs::AbstractArray{T,1}, xpred::AbstractArray{T,1};
			#sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(x)),
			sigmasq_cor::T = 1.0, rho::T=one(T))  where { T<:Real }
  make_kernel_obs_pred(xobs, xpred, kernel_matern32, sigmasq_cor=sigmasq_cor, rho=rho)
end

function make_kernel_dmatern32dx_obs_pred(xobs::AbstractArray{T,1}, xpred::AbstractArray{T,1};
			#sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(x)),
			sigmasq_cor::T = 1.0, rho::T=one(T))  where { T<:Real }
  make_kernel_obs_pred(xobs, xpred, dkerneldx_matern32, sigmasq_cor=sigmasq_cor, rho=rho)
end

function make_kernel_d2matern32dx2_obs_pred(xobs::AbstractArray{T,1}, xpred::AbstractArray{T,1};
			#sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(x)),
			sigmasq_cor::T = 1.0, rho::T=one(T))  where { T<:Real }
  make_kernel_obs_pred(xobs, xpred, d2kerneldx2_matern32, sigmasq_cor=sigmasq_cor, rho=rho)
end

function make_kernel_matern32_sparse_data(x::AbstractArray{T,1};
			sigmasq_obs::AbstractArray{T,1} = zeros(length(x)),
			sigmasq_cor::T = 1.0, rho::T=one(T))  where { T<:Real }
  make_kernel_data(x, matern32_sparse_kernel, sigmasq_obs=sigmasq_obs, sigmasq_cor=sigmasq_cor, rho=rho)
end

function make_kernel_matern32_sparse_obs_pred(xobs::AbstractArray{T,1}, xpred::AbstractArray{T,1};
			#sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(x)),
			sigmasq_cor::T = 1.0, rho::T=one(T))  where { T<:Real }
  make_kernel_obs_pred(xobs, xpred, matern32_sparse_kernel, sigmasq_cor=sigmasq_cor, rho=rho)
end

function make_kernel_dmatern32dx_sparse_obs_pred(xobs::AbstractArray{T,1}, xpred::AbstractArray{T,1};
			#sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(x)),
			sigmasq_cor::T = 1.0, rho::T=one(T))  where { T<:Real }
  make_kernel_obs_pred(xobs, xpred, dkerneldx_matern32_sparse, sigmasq_cor=sigmasq_cor, rho=rho)
end

function make_kernel_d2matern32dx2_sparse_obs_pred(xobs::AbstractArray{T,1}, xpred::AbstractArray{T,1};
			#sigmasq_obs::AbstractArray{T,1} = 1e-16*ones(length(x)),
			sigmasq_cor::T = 1.0, rho::T=one(T))  where { T<:Real }
  make_kernel_obs_pred(xobs, xpred, d2kerneldx2_matern32_sparse, sigmasq_cor=sigmasq_cor, rho=rho)
end

=#


text:
#=

    MainResult_Setup.jl

    Setup data, model, etc to produce the main result.

=#

using Inference
using Model

using CSV, DataFrames, DataFramesMeta
using Distributions
using StatsBase

##############################################################
## Data
##############################################################

# Load data
data = CSV.read("Data/CSV/DualMarker.csv",DataFrame)

# Autofluorescence distribution
E = [@subset(data,:Time .== 0.0, :Quenched .== false, :Temp .== 37.0).Signal_Cy5,
     @subset(data,:Time .== 0.0, :Quenched .== false, :Temp .== 37.0).Signal_BDP]

# Observation times
T = sort(unique(data.Time))[2:end]

# Format data to match `simulate_model` output
X = [[[@subset(data, :Time .== t, :Quenched .== qu, :Temp .== 37.0)[:,ch]
        for ch in [:Signal_Cy5,:Signal_BDP]]
        for qu in [true,false]]
        for t  in T]

# Obtain quenching efficiency
η = 1 - (@df @subset(data,:Time .== 120.0,:Temp .== 4.0,:Quenched) mean(:Signal_Cy5)) /
        (@df @subset(data,:Time .== 120.0,:Temp .== 4.0,(!).(:Quenched)) mean(:Signal_Cy5))

##############################################################
## Setup model
##############################################################

# Parameter distributions
function param_dist(ξ)
    μR,σR,μλ,σλ,ωλ,μβ,σβ,ωβ,ρRλ,ρRβ,ρ̄λβ = ξ

    # Marginal distributions
    R = LogNormalAlt(1.0,μR,σR)
    λ = GammaAlt(μλ,σλ,ωλ)
    β = GammaAlt(μβ,σβ,ωβ)

    # Construct copula from correlation matrix
    ρλβ = ρRλ * ρRβ + ρ̄λβ  * sqrt((1 - ρRλ^2) * (1 - ρRβ^2))
    C = GaussianCopula([ρRλ,ρRβ,ρλβ])

    # Return multivariate distribution (use interpolated quantile for λ and β)
    MvDependent(C,[R,λ,β],q=[false,true,true])
end

# Model
model = θ -> simulate_model(θ,T,param_dist,E;N=1000,η)

##############################################################
## Setup inference
##############################################################

# Pre-compute discrepancy metric
disc = makediscrepancy(X,w=[1.0,80.0])

# ABC distance metric
dist = θ -> (disc ∘ model)(θ)

# Prior
prior = Product([
    Uniform(5000,12000),# α₁
    Uniform(20.0,80.0), # α₂
    Uniform(0.0,10.0),   # σ₁
    Uniform(0.0,3.0),   # σ₂
    Uniform(0.0,1.0),   # μR
    Uniform(0.0,1.0),   # σR
    Uniform(0.0,0.5),   # μλ
    Uniform(0.0,0.3),   # σλ
    Uniform(-2.0,1.0),  # ωλ
    Uniform(0.0,0.20),  # μβ
    Uniform(0.0,0.10),  # σβ
    Uniform(-2.0,1.0),  # ωβ
    Uniform(-1.0,1.0),  # ρRλ
    Uniform(-1.0,1.0),  # ρRβ
    Uniform(-1.0,1.0),  # ρ̄λβ
    Uniform(0.0,0.15),  # p,
])

# Parameter names
param_names = ["α₁","α₂","σ₁","σ₂","μR","σR","μλ","σλ","ωλ","μβ","σβ","ωβ","ρRλ","ρRβ","ρ̄λβ","p"]

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
function popchart(; x::AbstractVector = Int[], y::AbstractVector = Int[], group::AbstractVector = [])

    v = VegaVisualization(width = 640, height = 400)
    add_data!(v, x = x, y = y, group = group)

    #Get unique table name
    table = v.data[1].name

    legend!(v)

    v.scales = Array(VegaScale, 3)
    v.scales[1] = VegaScale(name = "g", domain = [2,1], range = "width",  _type = "ordinal")
    v.scales[2] = VegaScale(name = "y", _type = "ordinal", range = "height", reverse = true, domain = VegaDataRef(data = table, field = "y"))
    v.scales[3] = VegaScale(name = "group", _type = "ordinal", domain = [1,2], range = ["#1f77b4", "#e377c2"])

    v.marks = Array(VegaMark,2)

    #Y-axis labels
    v.marks[2] = VegaMark(
                            _type = "text",
                            from = VegaMarkFrom(data = table,
                                                transform = [VegaTransform(Dict{Any, Any}("type" => "aggregate", "groupby" => ["y"]))]),
                            properties = VegaMarkProperties(
                                                            enter = VegaMarkPropertySet(
                                                                                x = VegaValueRef(field = Dict{Any, Any}("group" => "width"), mult = 0.5, offset = -15),
                                                                                y = VegaValueRef(scale = "y", field = "y", offset = 11),
                                                                                text = VegaValueRef(field = "y"),
                                                                                baseline = VegaValueRef(value = "middle"),
                                                                                align = VegaValueRef(value = "center"),
                                                                                fill = VegaValueRef(value = "#000")
                                                                                )
                                                            )
                        )

    v.marks[1] = VegaMark(
                            _type = "group",
                            from = VegaMarkFrom(data = table,
                                                transform = [VegaTransform(Dict{Any, Any}("type" => "facet", "groupby" => ["group"])),
                                                             VegaTransform(Dict{Any, Any}("type" => "formula", "field" => "reverse", "expr" => "datum.key==2"))]),
                            properties = VegaMarkProperties(update = VegaMarkPropertySet(
                                                                                x = VegaValueRef(scale = "g", field = "key"),
                                                                                y = VegaValueRef(value = 0),
                                                                                width = VegaValueRef(scale = "g", band = true, offset = -30),
                                                                                height = VegaValueRef(field = Dict{Any, Any}("group" => "height"))
                                                                                )
                            ),
                            scales = [VegaScale(
                                name = "x",
                                _type = "linear",
                                range = "width",
                                reverse = VegaDataRef(field = "reverse"),
                                nice = true,
                                domain = VegaDataRef(data = table, field = "x")
                                )

                            ],
                            axes = [VegaAxis(_type = "x", scale = "x", format = "s")],

                            #Actual bars for chart
                            marks = [VegaMark(
                            _type = "rect",
                            properties = VegaMarkProperties(enter = VegaMarkPropertySet(x = VegaValueRef(scale = "x", field = "x"),
                                                                                x2 = VegaValueRef(scale = "x", value = 0),
                                                                                y = VegaValueRef(scale = "y", field = "y"),
                                                                                height = VegaValueRef(scale = "y", band = true, offset = -1),
                                                                                fill = VegaValueRef(scale = "group", field = "group")

                                )

                                )

        )]

    )

    #Default to Paired color scale, 12
    colorscheme!(v; palette = ("Paired", 12))
    return v
end

max_stars_repo_name: jpthor/vega


text:
#PhysicalCommunications:
#-------------------------------------------------------------------------------

module PhysicalCommunications

include("functions.jl")
include("prbs.jl")
include("eyediag.jl")

#==Exported interface
===============================================================================#
export MaxLFSR #Create MaxLFSR iterator object.
export sequence #Builds sequences with LFSR objects
export sequence_detecterrors #Tests validity of bit sequence using sequence generator algorithm
export buildeye

#==Unexported interface
================================================================================
	.DataEye #Stores eye data
==#

end # module

#Last line


text:
module Gasdynamics1D

  #=
  Quasi-1d gas dynamics routines
  =#

  using Roots
  using Requires
  using Reexport

  @reexport using ThermofluidQuantities

  import ThermofluidQuantities: Unitful, Temperature, Density, Pressure, StagnationTemperature, StagnationDensity,
          StagnationPressure, SoundSpeed, Enthalpy, StagnationEnthalpy, InternalEnergy,
          StagnationInternalEnergy, Velocity, MachNumber, MassFlowRate, Entropy,
          HeatFlux, FrictionFactor, Area, FLOverD, AreaRatio, TemperatureRatio,
          PressureRatio, StagnationPressureRatio, DensityRatio, VelocityRatio

  include("plot_recipes.jl")

  export ThermodynamicProcess, Isentropic, NormalShock, FannoFlow, RayleighFlow
  export SubsonicMachNumber,SupersonicMachNumber
  export SubsonicPOverP0,SupersonicPOverP0
  export T0OverT,P0OverP,ρ0Overρ,AOverAStar,AStar
  export FLStarOverD,POverPStar,ρOverρStar,TOverTStar,P0OverP0Star
  export T0OverT0Star,VOverVStar

  function __init__()

    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin

      @reexport using LaTeXStrings

      Plots.default(markerstrokealpha = 0, legend = false,
      dpi = 100, size = (400, 300), grid = false)

    end
  end

  ###### THERMODYNAMIC PROCESSES #######

  abstract type ThermodynamicProcess end

  abstract type Isentropic <: ThermodynamicProcess end
  abstract type NormalShock <: ThermodynamicProcess end
  abstract type FannoFlow <: ThermodynamicProcess end
  abstract type RayleighFlow <: ThermodynamicProcess end

  include("thermodynamics.jl")
  include("isentropic.jl")
  include("normalshocks.jl")
  include("fanno.jl")
  include("rayleigh.jl")
  include("nozzle.jl")
  include("vectors.jl")

end


text:
using Plots
pgfplots()
include("../src/Plotting.jl")

x = read_HFSS_pattern("test/data/lens.csv")

sphericalpattern(x)
polarpattern(x[ϕ(At(90))],lims=(-20,20))

phi = 0:360
theta = 0:360

freq = 38u"GHz"
lens_diameter = 15.8u"mm"
zero = 0u"mm"
gap = 0.1u"mm"
positions = [(-(lens_diameter + gap),zero,zero),
             (-(lens_diameter + gap) / 2,-√(3) * (lens_diameter + gap) / 2,zero),
             ((lens_diameter + gap) / 2,-√(3) * (lens_diameter + gap) / 2,zero),
             (lens_diameter + gap,zero,zero),
             ((lens_diameter + gap) / 2,√(3) * (lens_diameter + gap) / 2,zero),
             (-(lens_diameter + gap) / 2,√(3) * (lens_diameter + gap) / 2,zero),
             (zero,zero,zero)]

af = ArrayFactor(positions,[1,1,1,1,1,1,1])

plot(af,legend=false)

pattern = [af(p,t,freq) for p ∈ phi, t ∈ theta]

buh = RadiationPattern(pattern,phi,theta)

r = ustrip(gain(buh[ϕ(At(0))]))
theta = deg2rad.(Array(buh.dims[2].val))

p = polarpattern(theta,r,lims=(-30,10))

r1 = cos.(theta)
r2 = 2*cos.(theta)

plot(theta,[r1 r2],proj=:polar)


text:
const data_gmm = [
    -1.6518947325446898,
    -0.9012309098240023,
    -0.443220750663359,
]
const s1_gmm = 5.0
const s2_gmm = 1.0

const data_hmm = [
    0.8207426252019316,
    2.081341502733774,
    1.0565299991122412,
]
const s1_hmm = 0.1
const s2_hmm = 0.5

# data for IMM from R. Neal paper:
const data_neal = [-1.48, -1.40, -1.16, -1.08, -1.02, 0.14, 0.51, 0.53, 0.78]
const α_neal = 10.0
const s1_imm = 5.0
const s2_imm = 1.0

####################################################################################################
####################################################################################################
@model function hierarchical_gaussian(x)
    λ ~ Gamma(2.0, inv(3.0))
    m ~ Normal(0, sqrt(1 / λ))
    x ~ Normal(m, sqrt(1 / λ))
end

@model function bernoulli_mixture(x)
    # Mixture prior.
    w ~ Dirichlet(2, 1/2)

    # Latent probability.
    π ~ DiscreteNonParametric([0.3, 0.7], w)

    # Observation.
    x ~ Bernoulli(π)
end

bernoulli_example(;x = false) = bernoulli_mixture(x)

@model function gmm(x, K)
    N = length(x)

    # Cluster association prior.
    w ~ Dirichlet(K, 1/K)

    # Cluster assignments.
    z ~ filldist(DiscreteNonParametric(1:K, w), N)

    # Cluster centers.
    μ ~ filldist(Normal(0.0, s1_gmm), K)

    # Observations.
    for n = 1:N
        x[n] ~ Normal(μ[z[n]], s2_gmm)
    end

    return x
end

gmm_example(;x = data_gmm, K = 2) = gmm(x, K)
gmm_generate(N; kwargs...) = gmm_generate(Random.GLOBAL_RNG, N; kwargs...)
gmm_generate(rng, N; K = 2) =
    :x => identity.(gmm(fill(missing, N), K)(rng, VarInfo(), SampleFromPrior(), DefaultContext()))

@model function gmm_loopy(x, K, ::Type{T}=Float64) where {T<:Real}
    N = length(x)

    # Cluster centers.
    μ = Vector{T}(undef, K)
    for k = 1:K
        μ[k] ~ Normal(0.0, s2_gmm)
    end

    # Cluster association prior.
    w ~ Dirichlet(K, 1/K)

    # Cluster assignments & observations.
    z = Vector{Int}(undef, N)
    for n = 1:N
        z[n] ~ DiscreteNonParametric(1:K, w)
        x[n] ~ Normal(μ[z[n]], s2_gmm)
    end

    return x
end

gmm_loopy_example(;x = data_gmm, K = 2) = gmm_loopy(x, K)

# same as gmm_loopy, but with an affine transformation on μ.
@model function gmm_shifted(x, K, ::Type{T}=Float64) where {T<:Real}
    N = length(x)

    # Cluster centers.
    μ = Vector{T}(undef, K)
    for k = 1:K
        μ[k] ~ Normal(0.0, s2_gmm)
    end

    # Cluster association prior.
    w ~ Dirichlet(K, 1/K)

    # Cluster assignments & observations.
    z = Vector{Int}(undef, N)
    for n = 1:N
        z[n] ~ DiscreteNonParametric(1:K, w)
        x[n] ~ Normal(4μ[z[n]] - 1, s2_gmm)
    end

    return x
end

gmm_shifted_example(;x = data_gmm, K = 2) = gmm_shifted(x, K)

@model function gmm_marginalized(x, K)
    N = length(x)

    # Cluster association prior.
    w ~ Dirichlet(K, 1/K)

    # Cluster centers.
    μ ~ filldist(Normal(0.0, s1_gmm), K)

    # Observations.
    for n = 1:N
        x[n] ~ MixtureModel(Normal, μ, w)
    end

    return x
end

gmm_marginalized_example(;x = data_gmm, K = 2) = gmm_marginalized(x, K)
gmm_marginalized_generate(N; kwargs...) = gmm_marginalized_generate(Random.GLOBAL_RNG, N; kwargs...)
gmm_marginalized_generate(rng, N; K = 2) =
    :x => identity.(gmm_marginalized(fill(missing, N), K)(rng, VarInfo(), SampleFromPrior(), DefaultContext()))

# K clusters, each one around i for i = 1:K with fixed variance
@model function hmm(x, K, ::Type{X}=Float64) where {X<:Real}
    N = length(x)

    # Transition matrix.
    T = Vector{Vector{X}}(undef, K)
    for i = 1:K
        T[i] ~ Dirichlet(K, 1/K)
    end

    # State sequence.
    s = zeros(Int, N)
    s[1] ~ Categorical(K)
    for i = 2:N
        s[i] ~ Categorical(T[s[i-1]])
    end

    # Emission matrix.
    m = Vector{X}(undef, K)
    for i = 1:K
        m[i] ~ Normal(i, s1_hmm)
    end

    # Observe each point of the input.
    x[1] ~ Normal(m[s[1]], s2_hmm)
    for i = 2:N
        x[i] ~ Normal(m[s[i]], s2_hmm)
    end

    return x
end

hmm_example(;x = data_hmm, K = 2) = hmm(x, K)
hmm_generate(N; kwargs...) = hmm_generate(Random.GLOBAL_RNG, N; kwargs...)
hmm_generate(rng, N; K = 2) =
    :x => identity.(hmm(fill(missing, N), K)(rng, VarInfo(), SampleFromPrior(), DefaultContext()))

@model function changepoint(y)
    N = length(y)
    α = 1 / mean(y)
    λ₁ ~ Exponential(α)
    λ₂ ~ Exponential(α)
    τ ~ DiscreteUniform(1, N)

    for n in 1:N
        y[n] ~ Poisson(τ > n ? λ₁ : λ₂)
    end

    return y
end

# @model function reverse_deps(x)
#     m = Vector{Float64}(undef, 2)
#     m[1] ~ Normal()
#     m[2] ~ Normal()
#     x ~ MvNormal(m)
# end

function stickbreak(v)
    K = length(v) + 1
    cumprod_one_minus_v = cumprod(1 .- v)

    return map(1:K) do k
        if k == 1
            v[1]
        elseif k == K
            cumprod_one_minus_v[K - 1]
        else
            v[k] * cumprod_one_minus_v[k - 1]
        end
    end
end

# from https://luiarthur.github.io/TuringBnpBenchmarks/dpsbgmm
@model function imm_stick(y, α, K)
    N = length(y)
    crm = DirichletProcess(α)

    v ~ filldist(StickBreakingProcess(crm), K - 1)
    w = stickbreak(v)

    # Cluster assignments
    z = zeros(Int, N)
    for n = 1:N
        z[n] ~ Categorical(w)
    end

    # Cluster centers
    L = identity(K)
    μ ~ filldist(Normal(0.0, s1_imm), L)

    # Observations
    for n = 1:N
        y[n] ~ Normal(μ[z[n]], s2_imm)
    end

    return y
end

imm_stick_example(;y = data_neal, α = α_neal, K = 10) = imm_stick(y, α, K)
imm_stick_generate(N; kwargs...) = imm_stick_generate(Random.GLOBAL_RNG, N; kwargs...)
imm_stick_generate(rng, N; α = α_neal, K = 10) =
    :y => identity.(imm_stick(fill(missing, N), α, K)(rng, VarInfo(), SampleFromPrior(), DefaultContext()))

####################################################################################################
####################################################################################################
@model function gmm_tarray(x, K)
    N = length(x)

    # Cluster association prior.
    w ~ Dirichlet(K, 1/K)

    # Cluster assignments.
    z = tzeros(Int, N)
    for n = 1:N
        z[n] ~ DiscreteNonParametric(1:K, w)
    end

    # Cluster centers.
    μ ~ filldist(Normal(0.0, s1_gmm), K)

    # Observations.
    for n = 1:N
        x[n] ~ Normal(μ[z[n]], s2_gmm)
    end

    return x
end

gmm_tarray_example(;x = data_gmm, K = 2) = gmm_tarray(x, K)

# K clusters, each one around i for i = 1:K with fixed variance
@model function hmm_tarray(x, K)
    N = length(x)

    # Transition matrix.
    T ~ filldist(Dirichlet(K, 1/K), K)

    # State sequence.
    s = tzeros(Int, N)
    s[1] ~ Categorical(K)
    for i = 2:N
        s[i] ~ Categorical(T[:, s[i-1]])
    end

    # Emission matrix.
    m ~ arraydist(Normal.(1:K, s1_hmm))

    # Observe each point of the input.
    x[1] ~ Normal(m[s[1]], s2_hmm)
    for i = 2:N
        x[i] ~ Normal(m[s[i]], s2_hmm)
    end

    return x
end

hmm_tarray_example(;x = data_hmm, K = 2) = hmm_tarray(x, K)

function stickbreak(v)
    K = length(v) + 1
    cumprod_one_minus_v = cumprod(1 .- v)

    return map(1:K) do k
        if k == 1
            v[1]
        elseif k == K
            cumprod_one_minus_v[K - 1]
        else
            v[k] * cumprod_one_minus_v[k - 1]
        end
    end
end

@model function imm_stick_tarray(y, α, K)
    N = length(y)
    crm = DirichletProcess(α)

    v ~ filldist(StickBreakingProcess(crm), K - 1)
    w = stickbreak(v)

    # Cluster assignments
    z = tzeros(Int, N)
    for n = 1:N
        z[n] ~ Categorical(w)
    end

    # Cluster centers
    L = identity(K)
    μ ~ filldist(Normal(0.0, s1_imm), L)

    # Observations
    for n = 1:N
        y[n] ~ Normal(μ[z[n]], s2_imm)
    end

    return y
end

imm_stick_tarray_example(;y = data_neal, α = α_neal, K = 10) = imm_stick_tarray(y, α, K)


text:
using Test
using Scientific_Programming_in_Julia
using FiniteDifferences

import Scientific_Programming_in_Julia.ReverseDiff
import Scientific_Programming_in_Julia.ReverseDiff.σ
import Scientific_Programming_in_Julia.ReverseDiff.data

@testset "ReverseDiff" begin
    # test hcat rule
    x = rand(3,1)
    y = rand(3,1)
    z = hcat(x,y)

    xt = track(x)
    yt = track(y)
    zt = hcat(xt,yt)
    zt.grad = z
    accum!(xt)
    accum!(yt)
    @test xt.data ≈ z[:,1]
    @test yt.data ≈ z[:,2]

    # Simple Neural Network
    W1, b1 = track(rand(3,4)), track(rand(3))
    W2, b2 = track(rand(2,3)), track(rand(2))
    xs = [rand(4) |> track for _ in 1:10]
    ys = [rand(2) |> track for _ in 1:10]

    function forward(x, W1, b1, W2, b2)
        y1 = σ(W1*x + b1)
        y2 = W2*y1 + b2
    end

    function loss(xs, ys, model...)
        errs = hcat([forward(x,model...)-y for (x,y) in zip(xs,ys)]...)
        errs |> abs2 |> sum
    end

    dr = ReverseDiff.gradient(b1->loss(xs,ys,W1,b1,W2,b2), b1)[1]
    Base.abs2(x::Matrix) = abs2.(x)
    df = grad(central_fdm(5,1), b1->loss(data.(xs),data.(ys),W1.data,b1,W2.data,b2.data), b1.data)[1]
    @test dr ≈ df

    # make sure that reset works
    ReverseDiff.reset!.([W1,b1,W2,b2])
    xs = [rand(4) |> track for _ in 1:10]
    ys = [rand(2) |> track for _ in 1:10]
    df = grad(central_fdm(5,1), b1->loss(data.(xs),data.(ys),W1.data,b1,W2.data,b2.data), b1.data)[1]
    dr = ReverseDiff.gradient(b1->loss(xs,ys,W1,b1,W2,b2), b1)[1]
    @test dr ≈ df
end


text:
## --- Get Crust 1.0

	# Have to download crust 1 before testing it!
	@test get_crust1() == 0

	# Note: Available layers in Crust 1:
	# 1) water
	# 2) ice
	# 3) upper sediments   (VP, VS, rho not defined in all cells)
	# 4) middle sediments  "
	# 5) lower sediments   "
	# 6) upper crystalline crust
	# 7) middle crystalline crust
	# 8) lower crystalline crust

	const lats = [43.70, 39.2508]
	const lons = [-72.29, -106.2925]

## -- Test seismic data

	# Test single lat-lon-depth point from crust1
	@test [x[1] for x in find_crust1_seismic(50,50,8)] == [7.1, 4.05, 3000.0]

	# Test fetching lat-lon array from crust1
	@test find_crust1_seismic(lats,lons,6) == ([6.3, 6.0], [3.63, 3.5], [2790.0, 2720.0])
	@test find_crust1_seismic(lats,lons,7) == ([6.6, 6.6], [3.80, 3.8], [2860.0, 2860.0])
	@test find_crust1_seismic(lats,lons,8) == ([7.0, 7.2], [3.99, 4.1], [2950.0, 3030.0])

## --  Test layer thickness data

	@test find_crust1_thickness(lats,lons,6) ≈ [15.4, 17.73]
	@test find_crust1_thickness(lats,lons,7) ≈ [13.57, 16.21]
	@test find_crust1_thickness(lats,lons,8) ≈ [7.70, 16.72]

## -- Test cumulative thickness data

	@test find_crust1_base(lats,lons,6) ≈ [-14.99, -14.91]
	@test find_crust1_base(lats,lons,7) ≈ [-28.56, -31.12]
	@test find_crust1_base(lats,lons,8) ≈ [-36.26, -47.84]

## -- Test complete data

	@test all(find_crust1_layer(lats,lons,6) .≈ ([6.3, 6.0], [3.63, 3.5], [2.79, 2.72], [15.4, 17.73]))
	@test all(find_crust1_layer(lats,lons,7) .≈ ([6.6, 6.6], [3.8, 3.8], [2.86, 2.86], [13.57, 16.21]))
	@test all(find_crust1_layer(lats,lons,8) .≈ ([7.0, 7.2], [3.99, 4.1], [2.95, 3.03], [7.70, 16.72]))

## ---


text:
@autodiff struct PlanarQuadrotor <: ContinuousDynamics
    mass::Float64  # mass (kg)
    g::Float64     # gravity (m/s²)
    ℓ::Float64     # tip to tip distance (m)
    J::Float64
    function PlanarQuadrotor(mass::Real, g::Real, ℓ::Real)
        J = 0.2*mass*ℓ^2
        new(mass,g,ℓ,J)
    end
end
PlanarQuadrotor(; mass=1.0, g=9.81, ℓ=0.3) = PlanarQuadrotor(mass, g, ℓ)

state_dim(::PlanarQuadrotor) = 6
control_dim(::PlanarQuadrotor) = 2

begin
    body = quote

        mass,g,ℓ,J = model.mass, model.g, model.ℓ, model.J

        θ = x[3]
        s,c = sincos(θ)
        ẍ = (1/mass)*(u[1] + u[2])*s
        ÿ = (1/mass)*(u[1] + u[2])*c - g
        θddot = (1/J)*(ℓ/2)*(u[2] - u[1])
    end
    @eval function dynamics(model::PlanarQuadrotor, x, u)
        $body
        return SA[x[4], x[5], x[6], ẍ, ÿ, θddot]
    end
    @eval function dynamics!(model::PlanarQuadrotor, xdot, x, u)
        $body
        xdot[1:3] .= @view x[4:6]
        xdot[4] = ẍ
        xdot[5] = ÿ
        xdot[6] = θddot
        return nothing
    end
end


text:
"""
# module parameters

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-06-07

Laminart parameters saved as named tuple, in float32 or float64

Values from:

R. D. S. Raizada and S. Grossberg, ‘Context-sensitive binding by the laminar circuits of V1 and V2: A unified model of perceptual grouping, attention, and orientation contrast’, Visual Cognition, vol. 8, no. 3–5, p. 431, Jun. 2001, doi: 10/d8phhc.

unless stated otherwise.

Functions para_var and para_var_k to replace parameters for testing
"""
module Parameters

p_temp_f32 = (σ_1 = 1.0f0, σ_2 = 0.5f0, H_σ_x = 3.0f0, H_σ_y = 0.8f0)

parameters_f32 = (
    K = 2,
    δ_v = 1.25f0,
    δ_c = 0.25f0,
    δ_m = 0.01875f0,
    δ_z = 0.125f0,
    δ_s = 2.5f0,
    σ_1 = p_temp_f32.σ_1,
    C_1 = 1.5f0,
    C_2 = 0.075f0,
    σ_2 = p_temp_f32.σ_2,
    γ = 10.0f0,
    α = 0.5f0,
    ϕ = 2.0f0,
    Γ = 0.2f0,
    v_21 = 0.0f0, #1.0
    μ = 2.0f0,
    ν = 1.1f0,
    n = 6.0f0,
    att = 0.0f0,# p25
    η_p = 2.1f0,
    η_m = 1.5f0,
    λ = 1.5f0,
    a_23_ex = 3.0f0,
    a_23_in = 0.5f0,
    v12_6 = 1.0f0,
    v12_4 = 5.0f0,
    ψ = 0.5f0,
    C_AB_l = 4 * ceil(Int, p_temp_f32.σ_2) + 1,
    H_σ_x = p_temp_f32.H_σ_x,
    H_σ_y = p_temp_f32.H_σ_y,
    H_fact = 13.0f0,
    # H_l =  4*ceil(Int, max(p_temp.H_σ_x, p_temp.H_σ_y))+1,
    H_l = 19,
    T_fact = [0.87f0, 0.13f0],       #avg TP same orient, other orient
    T_p_m = 0.302f0,    #avg TM/TP
    T_v2_fact = 0.625f0,     #T in V2 = T*
    # W, H sizes
    H_size = 19,
    W_size = 19,
    # W σ, estimated
    W_p_σ_x_same_a = 6f0,
    W_p_σ_x_same_b = 3f0,
    W_p_σ_y_same_a = 2f0,
    W_p_σ_y_same_b = 4f0,
    W_p_σ_opp_a = 4f0,
    W_p_σ_x_opp_b = 3f0,
    W_p_σ_x_opp_c = 1.2f0,
    W_p_σ_y_opp_b = 1.2f0,
    W_p_σ_y_opp_c = 3f0,
    W_p_same_fact = 39f0,
    W_p_opp_fact_a = 70f0,
    W_p_opp_fact_b = 50f0,
    W_m_σ_x_same_a = 6f0,
    W_m_σ_x_same_b = 3f0,
    W_m_σ_y_same_a = 2.5f0,
    W_m_σ_y_same_b = 4f0,
    W_m_σ_opp_a = 4f0,
    W_m_σ_x_opp_b = 3f0,
    W_m_σ_x_opp_c = 1.2f0,
    W_m_σ_y_opp_b = 1.2f0,
    W_m_σ_y_opp_c = 3f0,
    W_m_same_fact = 330f0,
    W_m_opp_fact_a = 70f0,
    W_m_opp_fact_b = 50f0,
    # customn parameters for controling feedback, lgn equlibrum
    lgn_equ_u = 1.0f0,
    lgn_equ_A = 0.0f0,
    lgn_equ_B = 0.0f0,
    # filling = "circular",
    # filling = Fill(0f0),
    W_l = 19,
)

function para_var(va::NamedTuple)
    """
    Takes in named tuple with C_1, C_2, ϕ, Γ, η_p, η_m, λ, ψ
    Returns parameters named tuple with these values replaced
    """
    p_temp_f32 = (σ_1 = 1.0f0, σ_2 = 0.5f0, H_σ_x = 3.0f0, H_σ_y = 0.8f0)

    p = (
        K = 2,
        δ_v = 1.25f0,
        δ_c = 0.25f0,
        δ_m = 0.01875f0,
        δ_z = 0.125f0,
        δ_s = 2.5f0,
        σ_1 = p_temp_f32.σ_1,
        C_1 = va.C_1,
        C_2 = va.C_2,
        σ_2 = p_temp_f32.σ_2,
        γ = 10.0f0,
        α = 0.5f0,
        ϕ = va.ϕ,
        Γ = va.Γ,
        v_21 = 0.0f0, #1.0
        μ = 2.0f0,
        ν = 1.1f0,
        n = 6.0f0,
        att = 0.0f0,# p25
        η_p = va.η_p,
        η_m = va.η_m,
        λ = va.λ,
        a_23_ex = 3.0f0,
        a_23_in = 0.5f0,
        v12_6 = 1.0f0,
        v12_4 = 5.0f0,
        ψ = va.ψ,
        C_AB_l = 4 * ceil(Int, p_temp_f32.σ_2) + 1,
        H_σ_x = p_temp_f32.H_σ_x,
        H_σ_y = p_temp_f32.H_σ_y,
        H_fact = 13.0f0,
        # H_l =  4*ceil(Int, max(p_temp.H_σ_x, p_temp.H_σ_y))+1,
        H_l = 19,
        T_fact = [0.87f0, 0.13f0],       #avg TP same orient, other orient
        T_p_m = 0.302f0,    #avg TM/TP
        T_v2_fact = 0.625f0,     #T in V2 = T*
        # W, H sizes
        H_size = 19,
        W_size = 19,
        # W σ
        W_p_σ_x_same_a = 6f0,
        W_p_σ_x_same_b = 3f0,
        W_p_σ_y_same_a = 2f0,
        W_p_σ_y_same_b = 4f0,
        W_p_σ_opp_a = 4f0,
        W_p_σ_x_opp_b = 3f0,
        W_p_σ_x_opp_c = 1.2f0,
        W_p_σ_y_opp_b = 1.2f0,
        W_p_σ_y_opp_c = 3f0,
        W_p_same_fact = 39f0,
        W_p_opp_fact_a = 70f0,
        W_p_opp_fact_b = 50f0,
        W_m_σ_x_same_a = 6f0,
        W_m_σ_x_same_b = 3f0,
        W_m_σ_y_same_a = 2.5f0,
        W_m_σ_y_same_b = 4f0,
        W_m_σ_opp_a = 4f0,
        W_m_σ_x_opp_b = 3f0,
        W_m_σ_x_opp_c = 1.2f0,
        W_m_σ_y_opp_b = 1.2f0,
        W_m_σ_y_opp_c = 3f0,
        W_m_same_fact = 330f0,
        W_m_opp_fact_a = 70f0,
        W_m_opp_fact_b = 50f0,
        # customn parameters for controling feedback, lgn equlibrum
        lgn_equ_u = 1.0f0,
        lgn_equ_A = 0.0f0,
        lgn_equ_B = 0.0f0,
        # filling = "circular",
        # filling = Fill(0f0),
        W_l = 19,
    )
    return p
end

function para_var_k(va::NamedTuple)
    """
    Takes in named tuple with H_fact, T_fact, T_p_m, W_p_same_fact, W_m_same_fact, H_l, W_l
    Returns parameters named tuple with these values replaced
    """
    p_temp_f32 = (σ_1 = 1.0f0, σ_2 = 0.5f0, H_σ_x = 3.0f0, H_σ_y = 0.8f0)

    p_temp_f32 = (σ_1 = 1.0f0, σ_2 = 0.5f0, H_σ_x = 3.0f0, H_σ_y = 0.8f0)

    p = (
        K = 2,
        δ_v = 1.25f0,
        δ_c = 0.25f0,
        δ_m = 0.01875f0,
        δ_z = 0.125f0,
        δ_s = 2.5f0,
        σ_1 = p_temp_f32.σ_1,
        C_1 = 1.5f0,
        C_2 = 0.075f0,
        σ_2 = p_temp_f32.σ_2,
        γ = 10.0f0,
        α = 0.5f0,
        ϕ = 2.0f0,
        Γ = 0.2f0,
        v_21 = 0.0f0, #1.0
        μ = 2.0f0,
        ν = 1.1f0,
        n = 6.0f0,
        att = 0.0f0,# p25
        η_p = 2.1f0,
        η_m = 1.5f0,
        λ = 1.5f0,
        a_23_ex = 3.0f0,
        a_23_in = 0.5f0,
        v12_6 = 1.0f0,
        v12_4 = 5.0f0,
        ψ = 0.5f0,
        C_AB_l = 4 * ceil(Int, p_temp_f32.σ_2) + 1,
        H_σ_x = p_temp_f32.H_σ_x,
        H_σ_y = p_temp_f32.H_σ_y,
        H_fact = va.H_fact,
        # H_l =  4*ceil(Int, max(p_temp.H_σ_x, p_temp.H_σ_y))+1,
        H_l = va.H_l,
        T_fact = va.T_fact,       #avg TP same orient, other orient
        T_p_m = va.T_p_m,    #avg TM/TP
        T_v2_fact = 0.625f0,     #T in V2 = T*
        # W, H sizes
        H_size = 19,
        W_size = 19,
        # W σ
        W_p_σ_x_same_a = 6f0,
        W_p_σ_x_same_b = 3f0,
        W_p_σ_y_same_a = 2f0,
        W_p_σ_y_same_b = 4f0,
        W_p_σ_opp_a = 4f0,
        W_p_σ_x_opp_b = 3f0,
        W_p_σ_x_opp_c = 1.2f0,
        W_p_σ_y_opp_b = 1.2f0,
        W_p_σ_y_opp_c = 3f0,
        W_p_same_fact = va.W_p_same_fact,
        W_p_opp_fact_a = 70f0,
        W_p_opp_fact_b = 50f0,
        W_m_σ_x_same_a = 6f0,
        W_m_σ_x_same_b = 3f0,
        W_m_σ_y_same_a = 2.5f0,
        W_m_σ_y_same_b = 4f0,
        W_m_σ_opp_a = 4f0,
        W_m_σ_x_opp_b = 3f0,
        W_m_σ_x_opp_c = 1.2f0,
        W_m_σ_y_opp_b = 1.2f0,
        W_m_σ_y_opp_c = 3f0,
        W_m_same_fact = va.W_m_same_fact,
        W_m_opp_fact_a = 70f0,
        W_m_opp_fact_b = 50f0,
        # customn parameters for controling feedback, lgn equlibrum
        lgn_equ_u = 1.0f0,
        lgn_equ_A = 0.0f0,
        lgn_equ_B = 0.0f0,
        # filling = "circular",
        # filling = Fill(0f0),
        W_l = va.W_l,
    )
    return p
end

p_temp_f64 = (σ_1 = 1.0, σ_2 = 0.5, H_σ_x = 3.0, H_σ_y = 4.0)

parameters_f64 = (
    K = 2,
    δ_v = 1.25,
    δ_c = 0.25,
    δ_m = 0.01875,
    δ_z = 0.125,
    δ_s = 2.5,
    σ_1 = p_temp_f64.σ_1,
    C_1 = 1.5,
    C_2 = 0.075,
    σ_2 = p_temp_f64.σ_2,
    γ = 10.0,
    α = 0.5,
    ϕ = 2.0,
    Γ = 0.2,
    v_21 = 0.0, #1.0
    μ = 2.0,
    ν = 1.1,
    n = 6.0,
    att = 0.0,# p25
    η_p = 2.1,
    η_m = 1.5,
    λ = 1.5,
    a_23_ex = 3.0,
    a_23_in = 0.5,
    v12_6 = 1.0,
    v12_4 = 5.0,
    ψ = 0.5,
    C_AB_l = 4 * ceil(Int, p_temp_f64.σ_2) + 1,
    H_σ_x = p_temp_f64.H_σ_x,
    H_σ_y = p_temp_f64.H_σ_y,
    H_fact = 5.0,
    # H_l =  4*ceil(Int, max(p_temp.H_σ_x, p_temp.H_σ_y))+1,
    H_l = 19,
    T_fact = [0.87, 0.13],       #avg TP same orient, other orient
    T_p_m = 0.302,    #avg TM/TP
    T_v2_fact = 0.625,     #T in V2 = T*
    # W, H sizes
    H_size = 19,
    W_size = 19,
    # customn parameters for controling feedback, lgn equlibrum
    lgn_equ_u = 1.0,
    lgn_equ_A = 0.0,
    lgn_equ_B = 0.0,
    # filling = "circular",
    # filling = Fill(0.0),
    W_l = 19,
)

end

# # T_P_11 = 0.9032
# # T_P_21 = 0.1384
# # T_P_12 = 0.1282
# # T_P_22 = 0.8443
# # T_M_11 = 0.2719
# # T_M_21 = 0.0428
# # T_M_12 = 0.0388
# # T_M_22 = 0.2506
#
# # T_P in V2 0.625x T_P in V1  ??

# T_m ~ T_p /3.2 calculated..

#
# const δ_v = 1.25
# const δ_c = 0.25
# const δ_m = 0.01875
# const δ_z = 0.125
# const δ_s = 2.5
# const δ_c = 0.25
#
# const σ_1 = 1
#
# const C_1 = 1.5
# const C_2 = 0.075
#
#
# const σ_2 = 0.5
# const K = 2
# const γ = 10
#
# const α = 0.5
#
# const ϕ = 2.0
#
# const Γ = 0.2
#
# const v_21 = 0 #1.0
#
# const μ = 2
# const ν = 1.1
# const n = 6
#
# const att = 0
# # p25
#
# const η_p = 2.1
#
# const η_m = 1.5
#
# const λ = 1.5
#
# const a_23_ex = 3
# const a_23_in = 0.5
#
# const v12_6 = 1
# const v12_4 = 5
#
# const ψ = 0.5


text:
module MochaModelFactories
    using Mocha
    export make_dense_nn, make_dense_nn_regression, lenet_5_nn, alex_cifar_tutorial_nn
    include("MochaModelFactories/dense_nn_factory.jl")
    include("MochaModelFactories/dense_nn_regression.jl")
    include("MochaModelFactories/lenet-5.jl")
    include("MochaModelFactories/alex_cifar_tutorial.jl")

    # TODO: Include any extra model factories here
end


text:
using Test

@testset "Compare to SALib.py" begin
    include("test_sobolsequence.jl")
    include("test_uniform.jl")
    include("test_nonuniform.jl")
end

@testset "Unit Testing" begin
    include("unit_tests.jl")
end


text:
# FWI on the 2D Overthrust model with L-BFGS from the NLopt library
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using Statistics, Random, LinearAlgebra
using JUDI, NLopt, HDF5, SegyIO, PyPlot

# Load starting model
if ~isfile("$(JUDI.JUDI_DATA)/overthrust_model_2D.h5")
    ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/2DFWI/overthrust_model_2D.h5")
end
n, d, o, m0 = read(h5open("$(JUDI.JUDI_DATA)/overthrust_model_2D.h5", "r"), "n", "d", "o", "m0")
model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)

# Bound constraints
vmin = 1.4f0
vmax = 6.5f0

# Load data and create data vector
if ~isfile("$(JUDI.JUDI_DATA)/overthrust_2D.segy")
    ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/2DFWI/overthrust_2D.segy")
end
block = segy_read("$(JUDI.JUDI_DATA)/overthrust_2D.segy")
d_obs = judiVector(block)

# Set up wavelet and source vector
src_geometry = Geometry(block; key = "source", segy_depth_key = "SourceDepth")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.008f0) # 8 Hz wavelet
q = judiVector(src_geometry, wavelet)

############################################## FWI #################################################

# optimization parameters
fevals = 10    # allow 10 function evalutations
batchsize = 40
count = 0
fhistory = zeros(Float32, fevals+1)

# NLopt objective function
println("No.  ", "fval         ", "norm(gradient)")
function f!(x, grad)

    # Update model
    model0.m .= x

    # Select batch and calculate gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    fval, gradient = fwi_objective(model0, q[i], d_obs[i])

    # Reset gradient in water column to zero
    gradient = reshape(gradient, model0.n); gradient[:, 1:21] .= 0f0
    grad[1:end] .= gradient[1:end]

    global count; count += 1
    println(count, "    ", fval, "    ", norm(grad))
    fhistory[count] = fval
    return convert(Float64, fval)
end

# Optimization parameters
opt = Opt(:LD_LBFGS, prod(model0.n))
min_objective!(opt, f!)
mmin = (1f0 ./ vmax)^2   # bound constraints on slowness squared
mmax = (1f0 ./ vmin)^2
lower_bounds!(opt, mmin); upper_bounds!(opt, mmax)
maxeval!(opt, fevals)
(minf, minx, ret) = optimize(opt, copy(model0.m))

# Save results, function values and elapsed time
h5open("result_2D_overthrust_lbfgs.h5", "w") do file
    write(file, "x", sqrt.(1f0 ./ reshape(minx, model0.n)), "fhistory", fhistory)
end

# Plot convergence and final result
figure(); plot(fhistory/norm(fhistory, Inf));
xlabel("Iteration no."); ylabel("Normalized residual"); title("Convergence of FWI w/ L-BFGS")
figure(); imshow(sqrt.(1f0 ./ adjoint(model0.m)), extent=(0, 20.0, 5.15, 0)); title("FWI with L-BFGS");
xlabel("Lateral position [km]"); ylabel("Depth [km]")


text:
function templateTestTforeach(arr::AbstractArray)::Bool
	serial = 0;
	parallel = Threads.Atomic{Int}(0);
	foreach(x->serial+=x, arr);
	tforeach(x->Threads.atomic_add!(parallel, x), arr);
	return serial == parallel[];
end

function testTforeach()::Bool
	return templateTestTforeach(randArray());
end

function testTforeachOffset()::Bool
	return templateTestTforeach(offsetArray());
end

function testTforeachEmpty()::Bool
	return templateTestTforeach([]);
end

#####

function templateTestTforeachMultiple(arr1::AbstractArray, arr2::AbstractArray)::Bool
	serial1 = 0;
	serial2 = 0;
	parallel1 = Threads.Atomic{Int}(0);
	parallel2 = Threads.Atomic{Int}(0);
	foreach((x, y)->begin serial1 += x; serial2 += y; end, arr1, arr2);
	foreach((x, y)->begin Threads.atomic_add!(parallel1, x); Threads.atomic_add!(parallel2, y); end, arr1, arr2);
	return serial1 == parallel1[] && serial2 == parallel2[];
end

function testTforeachMultiple()::Bool
	return templateTestTforeachMultiple(randArray(), randArray());
end

function testTforeachMultipleOffset()::Bool
	return templateTestTforeachMultiple(offsetArray(), offsetArray());
end

function testTforeachMultipleEmpty()::Bool
	return templateTestTforeachMultiple([], []);
end


text:
# trades are filled state

function transact()
    for i in 1:6
        produce(print_with_color(:green, "trade has been filled"))
    end
end

t = Task(transact)

# prime numbers are the transition to the next state, where nothing is filled

function statemachine()
    j = 1
      while j < 5
        for k in 13:51
            if isprime(k)
                consume(t)
                println("")
                sleep(1)
                j+=1
            else
                print_with_color(:red, @sprintf("%d ",k))
                println("")
            end
        end
    end
end


text:
# This file contains types and structures for DEA peers
"""
    AbstractDEAPeersDMU
An abstract type representing the DEA Peers of a DMU.
"""
abstract type AbstractDEAPeersDMU end

"""
    DEAPeersDMU
An data structure representing the DEA Peers of a DMU.
"""
struct DEAPeersDMU <: AbstractDEAPeersDMU
    i::Int64
    J::Vector{Int64}
    V::Vector{Float64}
    dmuname::String
    dmunamesref::Vector{String}
end

Base.size(P::DEAPeersDMU) = (length(P.J),)

Base.length(P::DEAPeersDMU) = length(P.J)

Base.eltype(P::DEAPeersDMU) = Tuple{Tuple{Int64,String},Float64}

Base.firstindex(P::DEAPeersDMU) = 1

Base.lastindex(P::DEAPeersDMU) = length(P)

Base.getindex(P::DEAPeersDMU, i::Int) = return (P.J[i], P.dmunamesref[i]), P.V[i]

Base.iterate(P::DEAPeersDMU, state=1) = state > length(P) ? nothing : (P[state], state + 1)

function Base.show(io::IO, x::DEAPeersDMU)
    compact = get(io, :compact, false)

    J = x.J
    V = x.V
    dmuname = x.dmuname
    dmunamesref = x.dmunamesref

    print(io, dmuname, ": ")

    for p in x
        print(io, p[1][2], " ( ", p[2], " ) ")
    end

end

"""
    AbstractDEAPeers
An abstract type representing a DEA Peers.
"""
abstract type AbstractDEAPeers end

"""
    DEAPeers
An data structure representing a DEA Peers.
"""
struct DEAPeers <: AbstractDEAPeers
    n::Int64
    nref::Int64
    lambda::SparseMatrixCSC{Float64, Int64}
    I::Vector{Int64}
    J::Vector{Int64}
    V::Vector{Float64}
    dmunames::Vector{String}
    dmunamesref::Vector{String}

    function DEAPeers(x::AbstractDEAModel; atol::Float64 = 1e-10, namesref::Union{Vector{String},Nothing} = nothing)
        if ! isdefined(x, :lambda)
            throw(ArgumentError("Model does not have info on peers"));
        end

        n = nobs(x)
        lambda = deepcopy(x.lambda)
        nref = size(lambda, 2)
        dmunames = names(x)
        dmunamesref = String[]

        # Remove very small numbers
        droptol!(lambda, atol)

        # If peers matrix is square
        if namesref === nothing
            if n == nref
                dmunamesref = dmunames
            else
                dmunamesref = ["Ref$i" for i in 1:nref]
            end
        elseif length(namesref) != nref
            throw(DimensionMismatch("length of `namesref` and number of reference DMUs are not equal"));
        else
            dmunamesref = namesref
        end

        # Extract I, J and V vectors from SparseArray and order by I and J
        I, J, V = findnz(lambda)
        npeers = length(I)
        IJV = [I J V]
        IJV = sortslices(IJV, dims = 1, by = x -> (x[1], x[2]))

        I = convert.(Int, IJV[:,1])
        J = convert.(Int, IJV[:,2])
        V = IJV[:,3]

        new(n, nref, lambda, I, J, V, dmunames, dmunamesref)

    end

end

Base.size(P::DEAPeers) = (P.n,)

Base.length(P::DEAPeers) = P.n

Base.eltype(P::DEAPeers) = DEAPeersDMU

Base.firstindex(P::DEAPeers) = 1

Base.lastindex(P::DEAPeers) = length(P)

function Base.getindex(P::DEAPeers, i::Int)::DEAPeersDMU

    I = P.I
    J = P.J
    V = P.V
    dmunames = P.dmunames
    dmunamesref = P.dmunamesref

    DEAPeersDMU(i, J[I .== i], V[I .== i], dmunames[i], dmunamesref[unique(J[I .== i])])

end

Base.iterate(P::DEAPeers, state = 1) = state > length(P) ? nothing : (P[state], state+1)

function Base.show(io::IO, x::DEAPeers)
    compact = get(io, :compact, false)

    n = length(x)

    println(io, "DEA Peers")

    for p in x
        println(io, p)
    end

end

"""
    peers(model::AbstractDEAModel)
Return peers of a DEA model.

# Optional Arguments
- `atol=1e-10`: tolerance for zero values.
- `namesref`: a vector of strings with the names of the decision making units in the reference set.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> peers(deaio)
DEA Peers
1: 1 ( 1.0 )
2: 4 ( 0.424978317432784 ) 7 ( 0.10928013876843023 )
3: 1 ( 1.134321653189578 ) 4 ( 0.43800539083557943 )
4: 4 ( 1.0 )
5: 4 ( 0.25738077214231636 ) 7 ( 0.04844814534443607 )
6: 7 ( 0.3333333333333333 )
7: 7 ( 1.0 )
8: 4 ( 1.0348650979425895 ) 7 ( 0.11457435012935832 )
9: 7 ( 1.1481481481481481 )
10: 4 ( 0.49056603773584906 ) 7 ( 0.4905660377358491 )
11: 11 ( 1.0 )
```
"""
peers(model::AbstractDEAModel; atol::Float64 = 1e-10, namesref::Union{Vector{String},Nothing} = nothing) = DEAPeers(model, atol = atol, namesref = namesref) ;

"""
    peersmatrix(model::AbstractDEAModel)
Return peers matrix of a DEA model.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> deaio = dea(X, Y);

julia> peersmatrix(deaio)
11×11 SparseMatrixCSC{Float64, Int64} with 17 stored entries:
 1.0       ⋅    ⋅    ⋅         ⋅    ⋅    ⋅            ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅    ⋅   0.424978   ⋅    ⋅   0.10928       ⋅    ⋅    ⋅    ⋅
 1.13432   ⋅    ⋅   0.438005   ⋅    ⋅    ⋅            ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅    ⋅   1.0        ⋅    ⋅   1.48479e-16   ⋅    ⋅    ⋅    ⋅
 ⋮                                 ⋮                                ⋮
  ⋅        ⋅    ⋅   1.03487    ⋅    ⋅   0.114574      ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅    ⋅    ⋅         ⋅    ⋅   1.14815       ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅    ⋅   0.490566   ⋅    ⋅   0.490566      ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅    ⋅    ⋅         ⋅    ⋅    ⋅            ⋅    ⋅    ⋅   1.0
```
"""
peersmatrix(model::AbstractDEAModel) = model.lambda ;

"""
    ispeer(P::DEAPeers, i::Int64, j::Int64)
Return `true` if `j` is peer of decision making unit `i`.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> P = peers(dea(X, Y));

julia> ispeer(P, 2, 4)
true
```
"""
function ispeer(P::DEAPeers, i::Int64, j::Int64)
    return P.lambda[i, j] != 0
end

"""
    ispeer(P::DEAPeers, i::String, j::String)
Return `true` if `j` is peer of decision making unit `i`.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> firms = ["A"; "B"; "C"; "D"; "E"; "F"; "G"; "H"; "I"; "J"; "K"]

julia> P = peers(dea(X, Y, names = firms));

julia> ispeer(P, "B", "D")
true
```
"""
function ispeer(P::DEAPeers, i::String, j::String)
    # Get corresponding id's of the names
    ival = findall(x -> x == i, P.dmunames)
    jval = findall(x -> x == j, P.dmunames)

    if length(ival) == 0
        throw(ArgumentError("Name $i does not exists"));
    end
    if length(jval) == 0
        throw(ArgumentError("Name $j does not exists"));
    end

    if length(ival) > 1
        throw(ArgumentError("Name $i is repeated. Search by name not possible."));
    end
    if length(jval) > 1
        throw(ArgumentError("Name $j is repeated. Search by name not possible."));
    end

    ival = ival[1]
    jval = jval[1]

    return ispeer(P, ival, jval)
end

"""
    ispeer(P::DEAPeersDMU, j::Int64)
Return `true` if `j` is peer.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> P = peers(dea(X, Y));

julia> P2 = P[2];

julia> ispeer(P2, 4)
true
```
"""
function ispeer(p::DEAPeersDMU, j::Int64)::Bool
    if j <= 0
        throw(BoundsError(p, j))
    end

    return any(p.J .== j)
end

"""
    ispeer(P::DEAPeers, j::String)
Return `true` if `j` is peer.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> firms = ["A"; "B"; "C"; "D"; "E"; "F"; "G"; "H"; "I"; "J"; "K"]

julia> P = peers(dea(X, Y, names = firms));

julia> P2 = P[2]

julia> ispeer(P2, "D")
true
```
"""
function ispeer(p::DEAPeersDMU, j::String)::Bool
    return any(p.dmunamesref .== j)
end

"""
    sum(P::DEAPeers; dims)

Sum elements of the peers matrix over the given dimension.

# Examples
```jldoctest
julia> X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17];

julia> Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12];

julia> P = peers(dea(X, Y));

julia> sum(P, dims = 2)
11×1 Matrix{Float64}:
 1.0
 0.5342584562012143
 1.5723270440251573
 1.0
 ⋮
 1.1494394480719479
 1.1481481481481481
 0.9811320754716981
 1.0
```
"""
Base.sum(P::DEAPeers; dims) = sum(P.lambda; dims = dims)

# CONVERT FUNCTIONS

Base.convert(::Type{Matrix}, x::DEAPeers) = convert(Matrix, x.lambda);

Base.convert(::Type{SparseMatrixCSC}, x::DEAPeers) = x.lambda;


text:
struct StraightRoadway
    length::Float64
end

function mod_position_to_roadway(s::Float64, roadway::StraightRoadway)
    while s > roadway.length
        s -= roadway.length
    end
    while s < 0
        s += roadway.length
    end
    return s
end

function get_headway(s_rear::Float64, s_fore::Float64, roadway::StraightRoadway)
    while s_fore < s_rear
        s_fore += roadway.length
    end
    while s_fore > s_rear + roadway.length
        s_fore -= roadway.length
    end
    return s_fore - s_rear # positive distance
end

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
abstract type AbstractDevice end

struct CPU <: AbstractDevice
end

struct GPU <: AbstractDevice
    id::Int
end

GPU() = GPU(1)

"""
Check if the argument is of type CuArray. Doesn't require CuArrays.jl to be loaded
"""
is_cuarray(x) = startswith(string(typeof(x)), "CuArray")

# function has_cuda_inputs(tape::Tape)
#     res = false
#     for op in tape
#         if op isa Input && op.val isa CuArray
#             res = true
#             break
#         end
#     end
#     return res
# end

# currently GPU's ID is just a placeholder
guess_device(args) = any(is_cuarray, args) ? GPU(1) : CPU()

"""
Retrieve function compatible with specified device

See also: to_device(device, f)
"""
device_function(device::CPU, f) = f

"""
Convert object to a compatible with the specified device.

For CPU it's usually no-op. For GPU behavior differs between object types:

 * Arrays are converted to CuArrays
 * structs are converted recursively
 * functions are looked up using `device_function()` or transformed using tracer
 * all other objects are returned as is
"""
to_device(device::CPU, x) = x
to_device(device::CPU, f::Function, args) = f
# see also cuda.jl

(device::CPU)(x) = to_device(device, x)
(device::GPU)(x) = to_device(device, x)


text:
using FixedPointNumbers, FileIO, ColorTypes#, ColorVectorSpace, ImageCore

function saveframe(frame::Array{T, 2}, path) where T<:Color
    save(path, frame)
end

function saveframe(frame::Array{T, 2}, path) where T
    saveframe(Gray.(frame), path)
end

function saveframe(frame::Array{T, 1}, path) where T
    # save 1d images horizontally
    saveframe(reshape(frame, 1, length(frame)), path)
end

function saveframe(frame::Array{T, 3}, path) where T
    colors = mapslices(rgb -> RGB(rgb...), frame, 3)
    saveframe(squeeze(colors, 3), path)
end

snapshot(frame, tag) = saveframe(frame, "snapshots/$tag.png")
snap(frame) = snapshot(frame, "latest")

function saveframes(frames)
    length(frames) == 0 && return
    (sy, sx) = size(first(frames))
    try
        open(`ffmpeg -y -f rawvideo -pix_fmt gray -s:v $(sy)x$(sx) -r 25 -i pipe:0 movie.mkv`, "w") do out
            for frame in frames
                write(out, reinterpret.(UInt8, N0f8.(frame)))
            end
        end
    catch err
        println("export to video failed: $(err)")
        println("you may need to install the `ffmpeg` command-line program.")
    end
end


text:
export loadImage, saveImage

function loadImage(filename::String)
  ni = niread(filename)
  pixelspacing = voxel_size(ni.header)
  return makeAxisArray(ni.raw, pixelspacing)
end

function saveImage(filename::String, im::AbstractArray{T,5}, makeAbs::Bool=false) where T
  voxel_size = ustrip.(uconvert.(u"m", pixelspacing(im)[1:3]))

  orientation = zeros(Float32, 3, 4)
  orientation[1,1] = orientation[2,2] = orientation[3,3] = 1.0

  ni = NIVolume(makeAbs ? abs.(im) : im; voxel_size=voxel_size,
                              orientation = orientation)
  # TODO add  dim_info=dim_info)
  niwrite(filename, ni)
end


text:
using Base: @deprecate, @deprecate_binding, depwarn

@deprecate varm(v::RealArray, m::Real, wv::AbstractWeights) varm(v, wv, m)
@deprecate varm(A::RealArray, M::RealArray, wv::AbstractWeights, dim::Int) varm(v, wv, m, dim)
@deprecate stdm(v::RealArray, m::Real, wv::AbstractWeights) stdm(v, wv, m)
@deprecate stdm(v::RealArray, m::RealArray, wv::AbstractWeights, dim::Int) stdm(v, wv, m, dim)

@deprecate trimmean(x::RealArray, p::Real) mean(trim(x, p/2))

@deprecate _moment2(v::RealArray, m::Real, wv::AbstractWeights) _moment2(v, wv, m)
@deprecate _moment3(v::RealArray, m::Real, wv::AbstractWeights) _moment3(v, wv, m)
@deprecate _moment4(v::RealArray, m::Real, wv::AbstractWeights) _moment4(v, wv, m)
@deprecate _momentk(v::RealArray, k::Int, m::Real, wv::AbstractWeights) _momentk(v, k, wv, m)
@deprecate moment(v::RealArray, k::Int, m::Real, wv::AbstractWeights) moment(v, k, wv, m)

@deprecate AIC(obj::StatisticalModel) aic(obj)
@deprecate AICc(obj::StatisticalModel) aicc(obj)
@deprecate BIC(obj::StatisticalModel) bic(obj)

if !isdefined(Base, :stderr)
    @deprecate stderr(obj::StatisticalModel) stderror(obj)
else
    function (io::typeof(stderr))(obj::StatisticalModel)
        Base.depwarn("stderr(obj::StatisticalModel) is deprecated, use stderror(obj) instead", :stderr)
        io === stderr ? stderror(obj) : throw(MethodErrror(io, (obj,)))
    end
end

@deprecate R2(obj::StatisticalModel, variant::Symbol) r2(obj, variant)
@deprecate R²(obj::StatisticalModel, variant::Symbol) r²(obj, variant)
@deprecate adjR2(obj::StatisticalModel, variant::Symbol) adjr2(obj, variant)
@deprecate adjR²(obj::StatisticalModel, variant::Symbol) adjr²(obj, variant)
@deprecate model_response(obj::StatisticalModel) response(obj)

@deprecate norepeats(a::AbstractArray) allunique(a)

function findat!(r::IntegerArray, a::AbstractArray{T}, b::AbstractArray{T}) where T
    Base.depwarn("findat! is deprecated, use indexin instead", :findat!)
    length(r) == length(b) || raise_dimerror()
    d = indexmap(a)
    @inbounds for i = 1 : length(b)
        r[i] = get(d, b[i], 0)
    end
    return r
end

"""
    findat(a, b)

For each element in `b`, find its first index in `a`. If the value does
not occur in `a`, the corresponding index is 0.
"""
findat(a::AbstractArray, b::AbstractArray) = findat!(Array{Int}(undef, size(b)), a, b)

@deprecate df(obj::StatisticalModel) dof(obj)
@deprecate df_residual(obj::StatisticalModel) dof_residual(obj)

@deprecate_binding WeightVec Weights

struct RandIntSampler  # for generating Int samples in [0, K-1]
    a::Int
    Ku::UInt
    U::UInt

    function RandIntSampler(K::Int)
        Base.depwarn("RandIntSampler is deprecated, use Base.Random.RangeGenerator instead",
                     :RandIntSampler)
        Ku = UInt(K)
        new(1, Ku, div(typemax(UInt), Ku) * Ku)
    end
    function RandIntSampler(a::Int, b::Int)
        Base.depwarn("RandIntSampler is deprecated, use Base.Random.RangeGenerator instead",
                     :RandIntSampler)
        Ku = UInt(b-a+1)
        new(a, Ku, div(typemax(UInt), Ku) * Ku)
    end
end

function rand(rng::AbstractRNG, s::RandIntSampler)
    x = rand(rng, UInt)
    while x >= s.U
        x = rand(rng, UInt)
    end
    s.a + Int(rem(x, s.Ku))
end
rand(s::RandIntSampler) = rand(Random.GLOBAL_RNG, s)

@deprecate randi(rng::AbstractRNG, K::Int) rand(rng, 1:K)
@deprecate randi(K::Int) rand(1:K)
@deprecate randi(rng::AbstractRNG, a::Int, b::Int) rand(rng, a:b)
@deprecate randi(a::Int, b::Int) rand(a:b)

@deprecate(mad!(v::AbstractArray{<:Real}, center;
                constant::Real = BigFloat("1.482602218505601860547076529360423431326703202590312896536266275245674447622701")),
           mad!(v, center=center, constant=constant))


text:
module CodeInfoTools

using Core: CodeInfo

import Base: iterate, push!, pushfirst!, insert!, delete!, getindex, lastindex, setindex!, display, +, length, identity
import Base: show

#####
##### Exports
#####

export var, Variable, Canvas, renumber, code_info, finish, get_slot

#####
##### Utilities
#####

const Variable = Core.SSAValue
var(id::Int) = Variable(id)

Base.:(+)(v::Variable, id::Int) = Variable(v.id + id)
Base.:(+)(id::Int, v::Variable) = Variable(v.id + id)

function code_info(f, tt; generated=true, debuginfo=:default)
    ir = code_lowered(f, tt; generated=generated, debuginfo=:default)
    isempty(ir) && return nothing
    return ir[1]
end

function get_slot(ci::CodeInfo, s::Symbol)
    ind = findfirst(el -> el == s, ci.slotnames)
    ind == nothing && return
    return Core.Compiler.SlotNumber(ind)
end

@doc(
"""
    code_info(f, tt; generate = true, debuginfo = :default)

Return lowered code for function `f` with tuple type `tt`. Equivalent to `InteractiveUtils.@code_lowered` -- but a function call and requires a tuple type `tt` as input.
""", code_info)

walk(fn, x) = fn(x)
walk(fn, x::Variable) = fn(x)
walk(fn, x::Core.ReturnNode) = Core.ReturnNode(walk(fn, x.val))
walk(fn, x::Core.GotoNode) = Core.GotoNode(walk(fn, x.label))
walk(fn, x::Core.GotoIfNot) = Core.GotoIfNot(walk(fn, x.cond), walk(fn, x.dest))
walk(fn, x::Expr) = Expr(x.head, map(a -> walk(fn, a), x.args)...)
function walk(fn, x::Vector)
    map(x) do el
        walk(fn, el)
    end
end

resolve(x) = x
resolve(gr::GlobalRef) = getproperty(gr.mod, gr.name)

#####
##### Canvas
#####

struct Canvas
    defs::Vector{Tuple{Int, Int}}
    code::Vector{Any}
    codelocs::Vector{Int32}
end
Canvas() = Canvas(Tuple{Int, Int}[], Any[], Int32[])

@doc(
"""
```julia
struct Canvas
    defs::Vector{Int}
    code::Vector{Any}
    codelocs::Vector{Int32}
end
Canvas() = Canvas(Int[], Any[], Int32[])
```

A `Vector`-like abstraction for `Core` code nodes.

Properties to keep in mind:

1. Insertion anywhere is slow.
2. Pushing to beginning is slow.
2. Pushing to end is fast.
3. Deletion is fast.
4. Accessing elements is fast.
5. Setting elements is fast.
6. Calling `renumber` must walk the entire `Canvas` instance to update SSA values -- slow.

Thus, if you build up a `Canvas` instance incrementally, everything should be fast.
""", Canvas)

length(c::Canvas) = length(filter(x -> x[2] > 0, c.defs))

function getindex(c::Canvas, idx::Int)
    r, ind = c.defs[idx]
    @assert ind > 0
    getindex(c.code, r)
end
getindex(c::Canvas, v::Variable) = getindex(c, v.id)

function push!(c::Canvas, stmt)
    push!(c.code, stmt)
    push!(c.codelocs, Int32(1))
    l = length(c.defs) + 1
    push!(c.defs, (l, l))
    return Variable(length(c.defs))
end

function insert!(c::Canvas, idx::Int, x)
    r, ind = c.defs[idx]
    @assert(ind > 0)
    push!(c.code, x)
    push!(c.codelocs, Int32(1))
    for i in 1 : length(c.defs)
        r, k = c.defs[i]
        if k > 0 && k >= ind
            c.defs[i] = (r, k + 1)
        end
    end
    push!(c.defs, (length(c.defs) + 1, ind))
    return Variable(length(c.defs))
end
insert!(c::Canvas, v::Variable, x) = insert!(c, v.id, x)

pushfirst!(c::Canvas, x) = insert!(c, 1, x)

setindex!(c::Canvas, x, v::Int) = setindex!(c.code, x, v)
setindex!(c::Canvas, x, v::Variable) = setindex!(c, x, v.id)

function delete!(c::Canvas, idx::Int)
    c.code[idx] = nothing
    c.defs[idx] = (idx, -1)
end
delete!(c::Canvas, v::Variable) = delete!(c, v.id)

_get(d::Dict, c, k) = c
_get(d::Dict, c::Variable, k) = haskey(d, c.id) ? Variable(getindex(d, c.id)) : nothing

function renumber(c::Canvas)
    s = sort(filter(v -> v[2] > 0, c.defs); by = x -> x[2])
    d = Dict((s[i][1], i) for  i in 1 : length(s))
    ind = first.(s)
    swap = walk(k -> _get(d, k, k), c.code)
    return Canvas(Tuple{Int, Int}[(i, i) for i in 1 : length(s)],
        getindex(swap, ind), getindex(c.codelocs, ind))
end

#####
##### Pretty printing
#####

print_stmt(io::IO, ex) = print(io, ex)
print_stmt(io::IO, ex::Expr) = print_stmt(io::IO, Val(ex.head), ex)

const tab = "  "

function show(io::IO, c::Canvas)
    indent = get(io, :indent, 0)
    bs = get(io, :bindings, Dict())
    for (r, ind) in sort(c.defs; by = x -> x[2])
        ind > 0 || continue
        println(io)
        print(io, tab^indent, "  ")
        print(io, string("%", r), " = ")
        ex = get(c.code, r, nothing)
        ex == nothing ? print(io, "nothing") : print_stmt(io, ex)
    end
end

print_stmt(io::IO, ::Val, ex) = print(io, ex)

function print_stmt(io::IO, ::Val{:enter}, ex)
    print(io, "try (outer %$(ex.args[1]))")
end

function print_stmt(io::IO, ::Val{:leave}, ex)
    print(io, "end try (start %$(ex.args[1]))")
end

function print_stmt(io::IO, ::Val{:pop_exception}, ex)
    print(io, "pop exception $(ex.args[1])")
end

#####
##### Pipe
#####

struct NewVariable
    id::Int
end

mutable struct Pipe
    from::CodeInfo
    to::Canvas
    map::Dict{Any, Any}
    var::Int
end

function Pipe(ci::CodeInfo)
    canv = Canvas()
    p = Pipe(ci, canv, Dict(), 0)
    return p
end

@doc(
"""
    Pipe(ir)

A wrapper around a `Canvas` object. Call [`finish`](@ref) when done to produce a new `CodeInfo` instance.
""", Pipe)

get_slot(p::Pipe, s::Symbol) = get_slot(p.from, s)

# This is used to handle NewVariable instances.
substitute!(p::Pipe, x, y) = (p.map[x] = y; x)
substitute(p::Pipe, x) = get(p.map, x, x)
substitute(p::Pipe, x::Expr) = Expr(x.head, substitute.((p, ), x.args)...)
substitute(p::Pipe, x::Core.GotoNode) = Core.GotoNode(substitute(p, x.label))
substitute(p::Pipe, x::Core.GotoIfNot) = Core.GotoIfNot(substitute(p, x.cond), substitute(p, x.dest))
substitute(p::Pipe, x::Core.ReturnNode) = Core.ReturnNode(substitute(p, x.val))

length(p::Pipe) = length(p.to)

getindex(p::Pipe, v) = getindex(p.to, v)
function getindex(p::Pipe, v::Union{Variable, NewVariable})
    tg = substitute(p, v)
    return getindex(p.to, tg)
end

lastindex(p::Pipe) = length(p.to)

function pipestate(ci::CodeInfo)
    ks = sort([Variable(i) => v for (i, v) in enumerate(ci.code)], by = x -> x[1].id)
    return first.(ks)
end

function iterate(p::Pipe, (ks, i) = (pipestate(p.from), 1))
    i > length(ks) && return
    v = ks[i]
    st = walk(resolve, p.from.code[v.id])
    substitute!(p, v, push!(p.to, substitute(p, st)))
    return ((v, st), (ks, i + 1))
end

var!(p::Pipe) = NewVariable(p.var += 1)

function Base.push!(p::Pipe, x)
    tmp = var!(p)
    v = push!(p.to, substitute(p, x))
    substitute!(p, tmp, v)
    return tmp
end

function Base.pushfirst!(p::Pipe, x)
    tmp = var!(p)
    v = pushfirst!(p.to, substitute(p, x))
    substitute!(p, tmp, v)
    return tmp
end

function setindex!(p::Pipe, x, v::Union{Variable, NewVariable})
    k = substitute(p, v)
    setindex!(p.to, substitute(p, x), k)
end

function insert!(p::Pipe, v::Union{Variable, NewVariable}, x; after = false)
    v′ = substitute(p, v).id
    x = substitute(p, x)
    tmp = var!(p)
    substitute!(p, tmp, insert!(p.to, v′ + after, x))
    return tmp
end

function Base.delete!(p::Pipe, v::Union{Variable, NewVariable})
    v′ = substitute(p, v)
    delete!(p.to, v′)
end

function finish(p::Pipe)
    new_ci = copy(p.from)
    c = renumber(p.to)
    new_ci.code = c.code
    new_ci.codelocs = c.codelocs
    new_ci.slotnames = p.from.slotnames
    new_ci.slotflags = [0x00 for _ in new_ci.slotnames]
    new_ci.inferred = false
    new_ci.inlineable = p.from.inlineable
    new_ci.ssavaluetypes = length(p.to)
    return new_ci
end

@doc(
"""
    finish(p::Pipe)

Create a new `CodeInfo` instance from a [`Pipe`](@ref). Renumbers the wrapped `Canvas` in-place -- then copies information from the original `CodeInfo` instance and inserts modifications from the wrapped `Canvas`.
""", finish)

Base.display(p::Pipe) = display(p.to)
function Base.identity(p::Pipe)
    for (v, st) in p
    end
    return p
end

end # module


text:
insts = OrderedDict()
insts["minimal"] = [
    # tr neglog
    ((3, MatLogdetCone()),),
    ((3, MatNegLog()),),
    ((3, MatNegLogEigOrd()),),
    ((3, MatNegLogDirect()),),
    # tr negentropy
    ((3, MatNegEntropy()),),
    ((3, MatNegEntropyEigOrd()),),
    # tr negsqrt
    ((3, MatNegSqrt()),),
    ((3, MatNegSqrtEigOrd()),),
    # tr negpower01
    ((3, MatNegPower01(0.7)),),
    ((3, MatNegPower01(0.7)),),
    # tr power12
    ((3, MatPower12(1.3)),),
    ((3, MatPower12EigOrd(1.3)),),
    ]
insts["fast"] = [
    ((40, MatNegGeom()),),
    ((15, MatNegGeomEFExp()),),
    ((15, MatNegGeomEFPow()),),
    ((40, MatNegSqrtConj()),),
    ((8, MatNegSqrtConjEigOrd()),),
    ((20, MatNegSqrtConjDirect()),),
    ((50, MatLogdetCone()),),
    ((30, MatNegLog()),),
    ((6, MatNegLogEigOrd()),),
    ((15, MatNegLogDirect()),),
    ((30, MatNegEntropy()),),
    ((8, MatNegEntropyEigOrd()),),
    ((30, MatPower12(1.5)),),
    ((7, MatPower12EigOrd(1.5)),),
    ]
insts["various"] = [
    ((100, MatNegGeom()),),
    ((50, MatNegGeomEFExp()),),
    ((30, MatNegGeomEFPow()),),
    ((100, MatNegSqrtConj()),),
    ((12, MatNegSqrtConjEigOrd()),),
    ((50, MatNegSqrtConjDirect()),),
    ((150, MatLogdetCone()),),
    ((100, MatNegLog()),),
    ((18, MatNegLogEigOrd()), nothing, (default_tol_relax = 1000,)),
    ((80, MatNegLogDirect()),),
    ((75, MatNegEntropy()),),
    ((14, MatNegEntropyEigOrd()),),
    ((30, MatPower12(1.5)),),
    ((6, MatPower12EigOrd(1.5)),),
    ]
return (CovarianceEstJuMP, insts)