text:
# This file is part of the TaylorSeries.jl Julia package, MIT license

# Printing of TaylorSeries objects

# subscriptify is taken from the ValidatedNumerics.jl package, licensed under MIT "Expat".
# superscriptify is a small variation

const subscript_digits = [c for c in "₀₁₂₃₄₅₆₇₈₉"]
const superscript_digits = [c for c in "⁰¹²³⁴⁵⁶⁷⁸⁹"]

function subscriptify(n::Int)
    dig = reverse(digits(n))
    join([subscript_digits[i+1] for i in dig])
end

function superscriptify(n::Int)
    dig = reverse(digits(n))
    join([superscript_digits[i+1] for i in dig])
end

#  Fallback
function pretty_print(a::Taylor1)
    # z = zero(a[0])
    var = _params_Taylor1_.var_name
    space = string(" ")
    bigO = bigOnotation[end] ?
        string("+ 𝒪(", var, superscriptify(a.order+1), ")") :
        string("")
    # iszero(a) && return string(space, z, space, bigO)
    strout::String = space
    ifirst = true
    for i in eachindex(a)
        monom::String = i==0 ? string("") : i==1 ? string(" ", var) :
            string(" ", var, superscriptify(i))
        @inbounds c = a[i]
        # c == z && continue
        cadena = numbr2str(c, ifirst)
        strout = string(strout, cadena, monom, space)
        ifirst = false
    end
    strout = strout * bigO
    strout
end

function pretty_print(a::Taylor1{T}) where {T<:NumberNotSeries}
    z = zero(a[0])
    var = _params_Taylor1_.var_name
    space = string(" ")
    bigO = bigOnotation[end] ?
        string("+ 𝒪(", var, superscriptify(a.order+1), ")") :
        string("")
    iszero(a) && return string(space, z, space, bigO)
    strout::String = space
    ifirst = true
    for i in eachindex(a)
        monom::String = i==0 ? string("") : i==1 ? string(" ", var) :
            string(" ", var, superscriptify(i))
        @inbounds c = a[i]
        iszero(c) && continue
        cadena = numbr2str(c, ifirst)
        strout = string(strout, cadena, monom, space)
        ifirst = false
    end
    strout = strout * bigO
    strout
end

function pretty_print(a::Taylor1{T} where {T <: AbstractSeries{S}}) where {S<:Number}
    z = zero(a[0])
    var = _params_Taylor1_.var_name
    space = string(" ")
    bigO = bigOnotation[end] ?
        string("+ 𝒪(", var, superscriptify(a.order+1), ")") :
        string("")
    iszero(a) && return string(space, z, space, bigO)
    strout::String = space
    ifirst = true
    for i in eachindex(a)
        monom::String = i==0 ? string("") : i==1 ? string(" ", var) :
            string(" ", var, superscriptify(i))
        @inbounds c = a[i]
        iszero(c) && continue
        cadena = numbr2str(c, ifirst)
        ccad::String = i==0 ? cadena : ifirst ? string("(", cadena, ")") :
            string(cadena[1:2], "(", cadena[3:end], ")")
        strout = string(strout, ccad, monom, space)
        ifirst = false
    end
    strout = strout * bigO
    strout
end

function pretty_print(a::HomogeneousPolynomial{T}) where {T<:Number}
    z = zero(a[1])
    space = string(" ")
    iszero(a) && return string(space, z)
    strout::String = homogPol2str(a)
    strout
end

function pretty_print(a::TaylorN{T}) where {T<:Number}
    z = zero(a[0])
    space = string("")
    bigO::String = bigOnotation[end] ?
        string(" + 𝒪(‖x‖", superscriptify(a.order+1), ")") :
        string("")
    iszero(a) && return string(space, z, space, bigO)
    strout::String = space
    ifirst = true
    for ord in eachindex(a)
        pol = a[ord]
        iszero(pol) && continue
        cadena::String = homogPol2str( pol )
        strsgn = (ifirst || ord == 0 || cadena[2] == '-') ?
            string("") : string(" +")
        strout = string( strout, strsgn, cadena)
        ifirst = false
    end
    strout = strout * bigO
    strout
end

function homogPol2str(a::HomogeneousPolynomial{T}) where {T<:Number}
    numVars = get_numvars()
    order = a.order
    z = zero(T)
    space = string(" ")
    strout::String = space
    ifirst = true
    iIndices = zeros(Int, numVars)
    for pos = 1:size_table[order+1]
        monom::String = string("")
        @inbounds iIndices[:] = coeff_table[order+1][pos]
        for ivar = 1:numVars
            powivar = iIndices[ivar]
            if powivar == 1
                monom = string(monom, name_taylorNvar(ivar))
            elseif powivar > 1
                monom = string(monom, name_taylorNvar(ivar), superscriptify(powivar))
            end
        end
        @inbounds c = a[pos]
        iszero(c) && continue
        cadena = numbr2str(c, ifirst)
        strout = string(strout, cadena, monom, space)
        ifirst = false
    end
    return strout[1:prevind(strout, end)]
end

function homogPol2str(a::HomogeneousPolynomial{Taylor1{T}}) where {T<:Number}
    numVars = get_numvars()
    order = a.order
    z = zero(a[1])
    space = string(" ")
    strout::String = space
    ifirst = true
    iIndices = zeros(Int, numVars)
    for pos = 1:size_table[order+1]
        monom::String = string("")
        @inbounds iIndices[:] = coeff_table[order+1][pos]
        for ivar = 1:numVars
            powivar = iIndices[ivar]
            if powivar == 1
                monom = string(monom, name_taylorNvar(ivar))
            elseif powivar > 1
                monom = string(monom, name_taylorNvar(ivar),
                    superscriptify(powivar))
            end
        end
        @inbounds c = a[pos]
        iszero(c) && continue
        cadena = numbr2str(c, ifirst)
        ccad::String = (pos==1 || ifirst) ? string("(", cadena, ")") :
            string(cadena[1:2], "(", cadena[3:end], ")")
        strout = string(strout, ccad, monom, space)
        ifirst = false
    end
    return strout[1:prevind(strout, end)]
end

function numbr2str(zz, ifirst::Bool=false)
    plusmin = ifelse( ifirst, string(""), string("+ ") )
    return string(plusmin, zz)
end

function numbr2str(zz::T, ifirst::Bool=false) where
        {T<:Union{AbstractFloat,Integer,Rational}}
    iszero(zz) && return string( zz )
    plusmin = ifelse( zz < zero(T), string("- "),
                ifelse( ifirst, string(""), string("+ ")) )
    return string(plusmin, abs(zz))
end

function numbr2str(zz::T, ifirst::Bool=false) where {T<:Complex}
    zT = zero(zz.re)
    iszero(zz) && return string(zT)
    zre, zim = reim(zz)
    cadena = string("")
    if zre > zT
        if ifirst
            cadena = string(" ( ", abs(zre)," ")
        else
            cadena = string(" + ( ", abs(zre)," ")
        end
        if zim > zT
            cadena = string(cadena, "+ ", abs(zim), " im )")
        elseif zim < zT
            cadena = string(cadena, "- ", abs(zim), " im )")
        else
            cadena = string(cadena, ")")
        end
    elseif zre < zT
        cadena = string(" - ( ", abs(zre), " ")
        if zim > zT
            cadena = string(cadena, "- ", abs(zim), " im )")
        elseif zim < zT
            cadena = string(cadena, "+ ", abs(zim), " im )")
        else
            cadena = string(cadena, ")")
        end
    else
        if zim > zT
            if ifirst
                cadena = string("( ", abs(zim), " im )")
            else
                cadena = string("+ ( ", abs(zim), " im )")
            end
        else
            cadena = string("- ( ", abs(zim), " im )")
        end
    end
    return cadena
end

name_taylorNvar(i::Int) = string(" ", get_variable_names()[i])

# summary
summary(a::Taylor1{T}) where {T<:Number} =
    string(a.order, "-order ", typeof(a), ":")

function summary(a::Union{HomogeneousPolynomial{T}, TaylorN{T}}) where {T<:Number}
    string(a.order, "-order ", typeof(a), " in ", get_numvars(), " variables:")
end

# show
function show(io::IO, a::Union{Taylor1, HomogeneousPolynomial, TaylorN})
    if _show_default[end]
        return Base.show_default(IOContext(io, :compact => false), a)
    else
        return print(io, pretty_print(a))
    end
end


text:
# helper to create fake data from scores such that the detectors return the scores again
to_table(points) = reshape(repeat(points, length(points)), (length(points), length(points))) |> table
const base = [1., 2, 3]
const test_scores = [base, # train and test equal
                     [1.,2], # len(test) < len(train)
                     [-Inf, 0., 1, 2, 3, 4, Inf], # len(test) > len(train)
                     [1.,1,1],  # len(test) = len(train)
                     [2.,2,2],
                     [3.,3,3]]
const scores0, scores1, scores2, scores3, scores4, scores5 = map(t -> (base, t), test_scores)
const test0, test1, test2, test3, test4, test5 = map(to_table, test_scores)
const train = to_table(base)
const train_labels = to_categorical(fill(missing, length(base)))

@testset "normalization, combination and classification" begin
    # fake detectors
    supervised = MinimalSupervised()
    unsupervised = MinimalUnsupervised()

    # fake machines
    up = machine(ProbabilisticDetector(unsupervised), train) |> fit!
    sp = machine(ProbabilisticDetector(supervised), train, train_labels) |> fit!
    ud = machine(DeterministicDetector(unsupervised), train) |> fit!
    sd = machine(DeterministicDetector(supervised), train, train_labels) |> fit!

    minmax_test = score_tuple -> scale_minmax(score_tuple...)[2]
    unify_test = score_tuple -> scale_unify(score_tuple...)[2]
    classify = score_tuple -> classify_percentile(DEFAULT_THRESHOLD)(scale_minmax(score_tuple...)...)[2]

    @testset "scores" begin
        raw_predict(detector, data) = OutlierDetection.raw_scores.(predict(detector, data))
        to_scores(data) = [transform(up, data), transform(sp, data), raw_predict(up, data), raw_predict(sp, data)]
        test_scoring(scores, data, labels) = @test all((labels,) .== [minmax_test(scores), to_scores(data)...])

        # Tests expected results on simple score vectors
        test_scoring(scores0, test0, [0.0, 0.5, 1.0])
        test_scoring(scores1, test1, [0.0, 0.5])
        test_scoring(scores2, test2, [0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0])
        test_scoring(scores3, test3, [0.0, 0.0, 0.0])
        test_scoring(scores4, test4, [0.5, 0.5, 0.5])
        test_scoring(scores5, test5, [1.0, 1.0, 1.0])

        # Min-Max normalization and unify should not change sort order
        @test sortperm(minmax_test(scores0)) == sortperm(unify_test(scores0))
        @test sortperm(minmax_test(scores1)) == sortperm(unify_test(scores1))
        @test sortperm(minmax_test(scores2)) == sortperm(unify_test(scores2))
        @test sortperm(minmax_test(scores3)) == sortperm(unify_test(scores3))
        @test sortperm(minmax_test(scores4)) == sortperm(unify_test(scores4))
        @test sortperm(minmax_test(scores5)) == sortperm(unify_test(scores5))
    end

    @testset "classification" begin
        to_classes(data) = [transform(ud, data), transform(sd, data), predict(ud, data), predict(sd, data)]
        test_classification(scores, data, labels) = @test all((labels,) .== [classify(scores), to_classes(data)...])

        # Test expected classification results with given threshold
        test_classification(scores0, test0, ["normal", "normal", "outlier"])
        test_classification(scores1, test1, ["normal", "normal"])
        test_classification(scores2, test2, ["normal", "normal", "normal", "normal", "outlier", "outlier", "outlier"])
        test_classification(scores3, test3, ["normal", "normal", "normal"])
        test_classification(scores4, test4, ["normal", "normal", "normal"])
        test_classification(scores5, test5, ["outlier", "outlier", "outlier"])
    end

    @testset "combination of scores" begin
        @test combine_mean(minmax_test.([scores3, scores4])...) == [0.25, 0.25, 0.25]
        @test combine_mean(minmax_test.([scores3, scores5])...) == [0.5, 0.5, 0.5]
        @test combine_mean(minmax_test.([scores4, scores5])...) == [0.75, 0.75, 0.75]

        @test combine_median(minmax_test.([scores3, scores4])...) == [0.25, 0.25, 0.25]
        @test combine_median(minmax_test.([scores3, scores5])...) == [0.5, 0.5, 0.5]
        @test combine_median(minmax_test.([scores4, scores5])...) == [0.75, 0.75, 0.75]

        @test combine_max(minmax_test.([scores3, scores4])...) == [0.5, 0.5, 0.5]
        @test combine_max(minmax_test.([scores3, scores5])...) == [1.0, 1.0, 1.0]
        @test combine_max(minmax_test.([scores4, scores5])...) == [1.0, 1.0, 1.0]
    end
end

ext: jl


text:
function FieldDateDetector(arg0::FieldAbsoluteDate)
    return FieldDateDetector((FieldAbsoluteDate,), arg0)
end

function FieldDateDetector(arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{FieldTimeStamped})
    return FieldDateDetector((RealFieldElement, RealFieldElement, Vector{FieldTimeStamped}), arg0, arg1, arg2)
end

function add_event_date(obj::FieldDateDetector, arg0::FieldAbsoluteDate)
    return jcall(obj, "addEventDate", void, (FieldAbsoluteDate,), arg0)
end

function g(obj::FieldDateDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "g", RealFieldElement, (FieldSpacecraftState,), arg0)
end

function get_date(obj::FieldDateDetector)
    return jcall(obj, "getDate", FieldAbsoluteDate, ())
end


text:
using StatFiles
using IteratorInterfaceExtensions
using TableTraits
using DataValues
using Test

@testset "StatFiles" begin

ar = load("types.dta") |> IteratorInterfaceExtensions.getiterator |> collect

@test length(ar) == 3
@test ar[1] == (vfloat=DataValue{Float32}(3.14), vdouble=DataValue(3.14), vlong=DataValue{Int32}(2), vint=DataValue{Int16}(2), vbyte=DataValue{Int8}(2), vstring=DataValue("2"))
@test ar[2] == (vfloat=DataValue{Float32}(7.), vdouble=DataValue(7.), vlong=DataValue{Int32}(7), vint=DataValue{Int16}(7), vbyte=DataValue{Int8}(7), vstring=DataValue("7"))
@test ar[3] == (vfloat=DataValue{Float32}(), vdouble=DataValue{Float64}(), vlong=DataValue{Int32}(), vint=DataValue{Int16}(), vbyte=DataValue{Int8}(), vstring=DataValue(""))
@test TableTraits.isiterabletable(load("types.dta")) == true
@test IteratorInterfaceExtensions.isiterable(load("types.dta")) == true

ar = load("types.sas7bdat") |> IteratorInterfaceExtensions.getiterator |> collect

@test length(ar) == 3
@test ar[1] == (vfloat=DataValue{Float32}(3.14), vdouble=DataValue(3.14), vlong=DataValue{Int32}(2), vint=DataValue{Int16}(2), vbyte=DataValue{Int8}(2), vstring=DataValue("2"))
@test ar[2] == (vfloat=DataValue{Float32}(7.), vdouble=DataValue(7.), vlong=DataValue{Int32}(7), vint=DataValue{Int16}(7), vbyte=DataValue{Int8}(7), vstring=DataValue("7"))
@test ar[3] == (vfloat=DataValue{Float32}(), vdouble=DataValue{Float64}(), vlong=DataValue{Int32}(), vint=DataValue{Int16}(), vbyte=DataValue{Int8}(), vstring=DataValue(""))
@test TableTraits.isiterabletable(load("types.sas7bdat")) == true
@test IteratorInterfaceExtensions.isiterable(load("types.sas7bdat")) == true

ar = load("types.sav") |> IteratorInterfaceExtensions.getiterator |> collect

@test length(ar) == 3
@test ar[1] == (vfloat=DataValue{Float32}(3.14), vdouble=DataValue(3.14), vlong=DataValue{Int32}(2), vint=DataValue{Int16}(2), vbyte=DataValue{Int8}(2), vstring=DataValue("2"))
@test ar[2] == (vfloat=DataValue{Float32}(7.), vdouble=DataValue(7.), vlong=DataValue{Int32}(7), vint=DataValue{Int16}(7), vbyte=DataValue{Int8}(7), vstring=DataValue("7"))
@test ar[3] == (vfloat=DataValue{Float32}(), vdouble=DataValue{Float64}(), vlong=DataValue{Int32}(), vint=DataValue{Int16}(), vbyte=DataValue{Int8}(), vstring=DataValue(""))
@test TableTraits.isiterabletable(load("types.sav")) == true
@test IteratorInterfaceExtensions.isiterable(load("types.sav")) == true

end

@testset "Show" begin

ar = load("types.dta")

@test sprint((stream,data)->show(stream, "text/html", data), ar) == "

vfloat	vdouble	vlong	vint	vbyte	vstring

3.14	3.14	2	2	2	"2"

7.0	7.0	7	7	7	"7"

#NA	#NA	#NA	#NA	#NA	""

"
@test sprint((stream,data)->show(stream, "application/vnd.dataresource+json", data), ar) ==
"{\"schema\":{\"fields\":[{\"name\":\"vfloat\",\"type\":\"number\"},{\"name\":\"vdouble\",\"type\":\"number\"},{\"name\":\"vlong\",\"type\":\"integer\"},{\"name\":\"vint\",\"type\":\"integer\"},{\"name\":\"vbyte\",\"type\":\"integer\"},{\"name\":\"vstring\",\"type\":\"string\"}]},\"data\":[{\"vfloat\":3.14,\"vdouble\":3.14,\"vlong\":2,\"vint\":2,\"vbyte\":2,\"vstring\":\"2\"},{\"vfloat\":7.0,\"vdouble\":7.0,\"vlong\":7,\"vint\":7,\"vbyte\":7,\"vstring\":\"7\"},{\"vfloat\":null,\"vdouble\":null,\"vlong\":null,\"vint\":null,\"vbyte\":null,\"vstring\":\"\"}]}"
@test sprint(show, ar) == "3x6 STATA file\nvfloat │ vdouble │ vlong │ vint │ vbyte │ vstring\n───────┼─────────┼───────┼──────┼───────┼────────\n3.14   │ 3.14    │ 2     │ 2    │ 2     │ \"2\"    \n7.0    │ 7.0     │ 7     │ 7    │ 7     │ \"7\"    \n#NA    │ #NA     │ #NA   │ #NA  │ #NA   │ \"\"     "

@test showable("text/html", ar) == true
@test showable("application/vnd.dataresource+json", ar) == true

ar = load("types.sas7bdat")

@test sprint((stream,data)->show(stream, "text/html", data), ar) == "

vfloat	vdouble	vlong	vint	vbyte	vstring

3.14	3.14	2.0	2.0	2.0	"2"

7.0	7.0	7.0	7.0	7.0	"7"

#NA	#NA	#NA	#NA	#NA	""

"
@test sprint((stream,data)->show(stream, "application/vnd.dataresource+json", data), ar) ==
"{\"schema\":{\"fields\":[{\"name\":\"vfloat\",\"type\":\"number\"},{\"name\":\"vdouble\",\"type\":\"number\"},{\"name\":\"vlong\",\"type\":\"number\"},{\"name\":\"vint\",\"type\":\"number\"},{\"name\":\"vbyte\",\"type\":\"number\"},{\"name\":\"vstring\",\"type\":\"string\"}]},\"data\":[{\"vfloat\":3.140000104904175,\"vdouble\":3.14,\"vlong\":2.0,\"vint\":2.0,\"vbyte\":2.0,\"vstring\":\"2\"},{\"vfloat\":7.0,\"vdouble\":7.0,\"vlong\":7.0,\"vint\":7.0,\"vbyte\":7.0,\"vstring\":\"7\"},{\"vfloat\":null,\"vdouble\":null,\"vlong\":null,\"vint\":null,\"vbyte\":null,\"vstring\":\"\"}]}"
@test sprint(show, ar) == "3x6 SAS file\nvfloat │ vdouble │ vlong │ vint │ vbyte │ vstring\n───────┼─────────┼───────┼──────┼───────┼────────\n3.14   │ 3.14    │ 2.0   │ 2.0  │ 2.0   │ \"2\"    \n7.0    │ 7.0     │ 7.0   │ 7.0  │ 7.0   │ \"7\"    \n#NA    │ #NA     │ #NA   │ #NA  │ #NA   │ \"\"     "

@test showable("text/html", ar) == true
@test showable("application/vnd.dataresource+json", ar) == true

ar = load("types.sav")

@test sprint((stream,data)->show(stream, "text/html", data), ar) == "

vfloat	vdouble	vlong	vint	vbyte	vstring

3.14	3.14	2.0	2.0	2.0	"2"

7.0	7.0	7.0	7.0	7.0	"7"

#NA	#NA	#NA	#NA	#NA	""

"
@test sprint((stream,data)->show(stream, "application/vnd.dataresource+json", data), ar) ==
"{\"schema\":{\"fields\":[{\"name\":\"vfloat\",\"type\":\"number\"},{\"name\":\"vdouble\",\"type\":\"number\"},{\"name\":\"vlong\",\"type\":\"number\"},{\"name\":\"vint\",\"type\":\"number\"},{\"name\":\"vbyte\",\"type\":\"number\"},{\"name\":\"vstring\",\"type\":\"string\"}]},\"data\":[{\"vfloat\":3.140000104904175,\"vdouble\":3.14,\"vlong\":2.0,\"vint\":2.0,\"vbyte\":2.0,\"vstring\":\"2\"},{\"vfloat\":7.0,\"vdouble\":7.0,\"vlong\":7.0,\"vint\":7.0,\"vbyte\":7.0,\"vstring\":\"7\"},{\"vfloat\":null,\"vdouble\":null,\"vlong\":null,\"vint\":null,\"vbyte\":null,\"vstring\":\"\"}]}"
@test sprint(show, ar) == "3x6 SPSS file\nvfloat │ vdouble │ vlong │ vint │ vbyte │ vstring\n───────┼─────────┼───────┼──────┼───────┼────────\n3.14   │ 3.14    │ 2.0   │ 2.0  │ 2.0   │ \"2\"    \n7.0    │ 7.0     │ 7.0   │ 7.0  │ 7.0   │ \"7\"    \n#NA    │ #NA     │ #NA   │ #NA  │ #NA   │ \"\"     "

@test showable("text/html", ar) == true
@test showable("application/vnd.dataresource+json", ar) == true

end


text:
function plot_circle!(plt, x, y, r; points=100, color=:black)
    rng = range(0; stop=2pi, length=points)
    plot!(plt, θ -> x + r*cos(θ), θ -> y + r*sin(θ), rng; color=color)
end

function plot_winding(winding, awg, maxturns, fill=1.0)

    coil=winding.coil

    width = (coil.od-coil.id)/2 |> ustrip
    h = ustrip(coil.len)
    plt=plot([0, h,h, 0, 0], [0, 0, width, width, 0]; aspect_ratio=:equal,legend=:none)

    d = awg.d_insul |>  unit(winding.coil.id) |> ustrip
    r = d/2
    d /= sqrt(fill)
    turns = 0
    foreach(enumerate(range(d/2;step=sqrt(3)/2*d,length=winding.winding.layers))) do (layer,y)
        x_offset = isodd(layer) ? d/2 : d
        layerturns = isodd(layer) ? winding.winding.odd_layer : winding.winding.even_layer
        foreach(range(x_offset;step=d,length=layerturns)) do x
            if turns < maxturns
                plot_circle!(plt,x,y,r)
                turns += 1
            else

            end
        end
    end
    plt
end

max_stars_repo_path: scripts/ploting utils.jl

max_issues_repo_path: scripts/ploting utils.jl

max_forks_repo_path: scripts/ploting utils.jl


text:
include("rates.jl")
include("radiation.jl")
include("broadening.jl")
include("populations.jl")
include("voronoi_utils.jl")
include("characteristics.jl")
include("irregular_ray_tracing.jl")

function J_λ_regular(S_λ::Array{<:UnitsIntensity_λ, 4},
                     α_cont::Array{<:PerLength, 3},
                     populations::Array{<:NumberDensity, 4},
                     atmos::Atmosphere,
                     line::HydrogenicLine,
                     quadrature::String)

    # Ω = (θ, φ), space angle
    weights, θ_array, ϕ_array, n_angles = read_quadrature(quadrature)

    J_λ = zero(S_λ)

    γ = γ_constant(line,
                   atmos.temperature,
                   (populations[:, :, :, 1] .+ populations[:, :, :, 2]),
                   atmos.electron_density)

    damping_λ = Array{Float64, 4}(undef, size(S_λ))
    Threads.@threads for l in eachindex(line.λ)
        damping_λ[l, :, :, :] = damping.(γ, line.λ[l], line.ΔD)
    end

    # total extinction
    # α_tot = Array{Float64, 4}(undef, size(damping_λ))u"m^-1"

    for i in 1:n_angles
        θ = θ_array[i]
        ϕ = ϕ_array[i]
        k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]

        profile = compute_voigt_profile(line, atmos, damping_λ, k)

        Threads.@threads for l in eachindex(line.λ)

            α_tot = αline_λ(line,
                            profile[l, :, :, :],
                            populations[:, :, :, 2],
                            populations[:, :, :, 1]) + α_cont

            if θ > 90
                I_0 =  B_λ.(line.λ[l], atmos.temperature[1,:,:])
                J_λ[l,:,:,:] .+= weights[i].*short_characteristics_up(k,
                                                                      S_λ[l,:,:,:],
                                                                      I_0,
                                                                      α_tot,
                                                                      atmos;
                                                                      n_sweeps=3)
            elseif θ < 90
                I_0 = zero(S_λ[l, 1, :, :])
                J_λ[l,:,:,:] .+= weights[i].*short_characteristics_down(k,
                                                                        S_λ[l,:,:,:],
                                                                        I_0,
                                                                        α_tot,
                                                                        atmos;
                                                                        n_sweeps=3)
            end
        end
    end

    return J_λ, damping_λ
end

function J_λ_voronoi(S_λ::Matrix{<:UnitsIntensity_λ},
                     α_cont::Vector{<:PerLength},
                     populations::Matrix{<:NumberDensity},
                     sites::VoronoiSites,
                     line::HydrogenicLine,
                     quadrature::String)

    # Ω = (θ, φ), space angle
    weights, θ_array, ϕ_array, n_points = read_quadrature(quadrature)

    J_λ = zero(S_λ)

    γ = γ_constant(line,
                   sites.temperature,
                   (populations[:, 1] .+ populations[:, 2]),
                   sites.electron_density)

    damping_λ = Matrix{Float64}(undef, size(S_λ))
    Threads.@threads for l in eachindex(line.λ)
       damping_λ[l, :] = damping.(γ, line.λ[l], line.ΔD)
    end

    n_sweeps = 3

    for i in 1:n_points
        θ = θ_array[i]
        ϕ = ϕ_array[i]
        k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]

        profile = compute_voigt_profile(line, sites, damping_λ, k)

        Threads.@threads for l in eachindex(line.λ)

            α_tot = αline_λ(line,
                            profile[l, :],
                            populations[:, 2],
                            populations[:, 1]) + α_cont

            if θ_array[i] > 90
                bottom_layer = sites.layers_up[2] - 1
                bottom_layer_idx = sites.perm_up[1:bottom_layer]
                I_0 = B_λ.(line.λ[l], sites.temperature[bottom_layer_idx])
                J_λ[l,:] += weights[i]*Delaunay_upII(k, S_λ[l,:], I_0, α_tot, sites, n_sweeps)

            elseif θ_array[i] < 90
                top_layer = sites.layers_down[2] - 1
                I_0 = zeros(top_layer)u"kW*nm^-1*m^-2"
                J_λ[l,:] += weights[i]*Delaunay_downII(k, S_λ[l,:], I_0, α_tot, sites, n_sweeps)
            end

        end
    end
    return J_λ, damping_λ
end

function Λ_regular(ϵ::AbstractFloat,
                   maxiter::Integer,
                   atmos::Atmosphere,
                   line::HydrogenicLine,
                   quadrature::String,
                   DATA::String)

    # LTE populations
    LTE_pops = LTE_populations(line, atmos)

    # Initial populations
    populations = copy(LTE_pops)
    # zero_radiation_populations(line, atmos)

    # Find continuum extinction and absorption extinction (without Thomson and Rayleigh)
    α_cont = α_absorption.(line.λ0,
                           atmos.temperature,
                           atmos.electron_density*1.0,
                           LTE_pops[:,:,:,1].+LTE_pops[:,:,:,2],
                           LTE_pops[:,:,:,3]) .+
             α_scattering.(line.λ0,
                           atmos.electron_density,
                           LTE_pops[:,:,:,1])

    # destruction probability (Should I include line???)
    ελ = destruction(LTE_pops, atmos.electron_density, atmos.temperature, line)
    println("Minimum $(minimum(ελ)) destruction probability")

    # Start with the source function as the Planck function
    B_0 = Array{Float64, 4}(undef, (length(line.λ), size(α_cont)...))u"kW*m^-2*nm^-1"
    for l in eachindex(line.λ)
        B_0[l, :, :, :] = B_λ.(line.λ[l], atmos.temperature)
    end

    S_new = copy(B_0)

    S_old = zero(S_new)

    C = calculate_C(atmos, LTE_pops)

    i=0

    local J_new
    # check where ε < 5e-3, cut above heights
    while criterion(S_new, S_old, ϵ, i, maxiter, DATA)
        #############################
        # Calculate radiation field #
        #############################
        S_old = copy(S_new)
        J_new, damping_λ = J_λ_regular(S_old, α_cont, populations,
                                       atmos, line, quadrature)

        for l in eachindex(line.λ)
            S_new[l,:,:,:] = (1 .- ελ).*J_new[l,:,:,:] .+ ελ.*B_0[l,:,:,:]
        end

        #############################
        #      Calculate rates      #
        #############################
        R = calculate_R(atmos, line, J_new, damping_λ, LTE_pops)

        #############################
        #    Update populations     #
        #############################
        populations = get_revised_populations(R, C, atmos.hydrogen_populations*1.0)

        ############################
        #Write current data to file#
        ############################
        write_to_file(populations[:, 2:end-1, 2:end-1, :], DATA)
        write_to_file(S_new[:, :, 2:end-1, 2:end-1], DATA)

        ############################
        #     Iteration is done    #
        ############################
        i+=1
    end

    if i == maxiter
        println("Did not converge inside scope")
        return J_new, S_new, α_cont, populations
    end

    println("Converged in $i iterations")

    return J_new, S_new, α_cont, populations
end

function Λ_voronoi(ϵ::AbstractFloat,
                   maxiter::Integer,
                   sites::VoronoiSites,
                   line::HydrogenicLine,
                   quadrature::String,
                   DATA::String)
    println("---Iterating---")

    # Start in LTE
    LTE_pops = LTE_populations(line, sites)

    # Initial populations
    populations = copy(LTE_pops)
    # zero_radiation_populations(line, sites)

    # Find continuum extinction and absorption extinction (only with Thomson and Rayleigh)
    α_cont = α_absorption.(line.λ0,
                           sites.temperature,
                           sites.electron_density*1.0,
                           LTE_pops[:,1].+LTE_pops[:,2],
                           LTE_pops[:,3]) .+
             α_scattering.(line.λ0,
                           sites.electron_density,
                           LTE_pops[:,1])

    # Start with the source function as the Planck function
    B_0 = Array{Float64, 2}(undef, (length(line.λ), sites.n))u"kW*m^-2*nm^-1"

    for l in eachindex(line.λ)
        B_0[l, :] = B_λ.(line.λ[l], sites.temperature)
    end

    S_new = copy(B_0)

    S_old = zero(S_new)

    # destruction probability (Should I include line???)
    ελ = destruction(LTE_pops, sites.electron_density, sites.temperature, line)
    println("Minimum $(minimum(ελ)) destruction probability")

    C = calculate_C(sites, LTE_pops)

    i=0

    local J_new
    # check where ε < 1e-2, cut above heights
    while criterion(S_new, S_old, ϵ, i, maxiter, DATA)
        #############################
        # Calculate radiation field #
        #############################
        S_old = copy(S_new)
        J_new, damping_λ = J_λ_voronoi(S_old, α_cont, populations,
                                       sites, line, quadrature)

        for l in eachindex(line.λ)
            S_new[l,:] = (1 .- ελ).*J_new[l,:] .+ ελ.*B_0[l,:]
        end

        #############################
        #      Calculate rates      #
        #############################
        R = calculate_R(sites, line, J_new, damping_λ, LTE_pops)

        #############################
        #    Update populations     #
        #############################
        populations = get_revised_populations(R, C, sites.hydrogen_populations)

        ############################
        #Write current data to file#
        ############################
        write_to_file(populations, DATA)
        write_to_file(S_new, DATA)

        ############################
        #     Iteration is done    #
        ############################
        i+=1
    end

    if i == maxiter
        println("Did not converge inside scope")
        return J_new, S_new, α_cont, populations
    end

    println("Converged in $i iterations")

    return J_new, S_new, α_cont, populations
end

function criterion(S_new::Array{<:UnitsIntensity_λ, 4},
                   S_old::Array{<:UnitsIntensity_λ, 4},
                   ϵ::Float64,
                   i::Int,
                   maxiter::Int,
                   DATA::String)
    diff = 0
    nλ = size(S_new)[1]
    for l in 1:nλ
        l_diff = maximum(abs.(1 .- S_old[l,:,:,:]./S_new[l,:,:,:])) |> Unitful.NoUnits
        diff = max(diff, l_diff)
        if isnan(l_diff)
            println("NaN DIFF!, index $l")
        end
    end
    if i > 0
        println("   Rel. diff.: $diff")
    end
    println("Iteration $(i+1)...")

    # Write convergence history to file
    write_to_file(diff, i+1, DATA)

    diff > ϵ && i < maxiter
end

function criterion(S_new::Matrix{<:UnitsIntensity_λ},
                   S_old::Matrix{<:UnitsIntensity_λ},
                   ϵ::Float64,
                   i::Int,
                   maxiter::Int,
                   DATA::String)
    diff = 0
    nλ = size(S_new)[1]
    for l in 1:nλ
        l_diff = maximum(abs.(1 .- S_old[l, :]./S_new[l, :])) |> Unitful.NoUnits
        diff = max(diff, l_diff)
        if isnan(l_diff)
            println("NaN DIFF!, index $l")
        end
    end
    if i > 0
        println("   Rel. diff.: $diff")
    end
    println("Iteration $(i+1)...")

    # Write convergence history to file
    write_to_file(diff, i+1, DATA)

    diff > ϵ && i < maxiter
end

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
struct BACKSAFTParam <: EoSParam
    segment::SingleParam{Float64}
    sigma::PairParam{Float64}
    epsilon::PairParam{Float64}
    c::SingleParam{Float64}
    alpha::SingleParam{Float64}
end

abstract type BACKSAFTModel <: SAFTModel end
@newmodel BACKSAFT BACKSAFTModel BACKSAFTParam

export BACKSAFT
function BACKSAFT(components; idealmodel=BasicIdeal, userlocations=String[], ideal_userlocations=String[], verbose=false)
    params = getparams(components, ["SAFT/BACKSAFT","properties/molarmass.csv"]; userlocations=userlocations, verbose=verbose)
    segment = params["m"]
    c = params["c"]
    k = params["k"]
    alpha = params["alpha"]
    sigma = params["vol"]
    sigma.values .*= 6/N_A/1e6/π
    sigma.values .^= 1/3
    sigma = sigma_LorentzBerthelot(sigma)
    epsilon = epsilon_LorentzBerthelot(params["epsilon"], k)
    packagedparams = BACKSAFTParam(segment, sigma, epsilon, c, alpha)
    references = ["TODO BACKSAFT", "TODO BACKSAFT"]

    model = BACKSAFT(packagedparams, idealmodel; ideal_userlocations=ideal_userlocations, references=references, verbose=verbose)
    return model
end

function a_res(model::BACKSAFTModel ,V, T, z)
    a_hcb_ = @f(a_hcb)
    a_disp_ = @f(a_disp)
    a_chain_ = @f(a_chain)
    return  a_hcb_ + a_chain_ + (1.75*(a_chain_/a_hcb_)+1)*a_disp_
end

function a_hcb(model::BACKSAFTModel, V, T, z)
    α = model.params.alpha.values[1]
    m = model.params.segment.values[1]
    η = @f(ζ,3)
    return m*(α^2/(1-η)^2-(α^2-3α)/(1-η)-(1-α^2)*log(1-η)-3α)
end

function a_disp(model::BACKSAFTModel, V, T, z)
    m = model.params.segment.values[1]
    c = model.params.c.values[1]
    u = model.params.epsilon.values[1]*(1+c/T)
    η = @f(ζ,3)
    τ = 0.74048
    D1 = BACKSAFT_consts.D1
    D2 = BACKSAFT_consts.D2
    D3 = BACKSAFT_consts.D3
    D4 = BACKSAFT_consts.D4
    A1 = ∑(D1[j]*(u/T)*(η/τ)^j for j ∈ 1:6)
    A2 = ∑(D2[j]*(u/T)^2*(η/τ)^j for j ∈ 1:9)
    A3 = ∑(D3[j]*(u/T)^3*(η/τ)^j for j ∈ 1:5)
    A4 = ∑(D4[j]*(u/T)^4*(η/τ)^j for j ∈ 1:4)
    return m*(A1+A2+A3+A4)
end

function d(model::BACKSAFTModel, V, T, z, i)
    ϵ = model.params.epsilon.diagvalues
    σ = model.params.sigma.diagvalues
    return σ[i] * (1 - 0.12exp(-3ϵ[i]/T))
end

function ζ(model::BACKSAFTModel, V, T, z, n)
    ∑z = ∑(z)
    x = z/∑z
    m = model.params.segment.values
    return N_A*∑z*π/6/V * ∑(x[i]*m[i]*@f(d,i)^n for i ∈ @comps)
end

function a_chain(model::BACKSAFTModel, V, T, z)
    m = model.params.segment.values[1]
    return (1-m)*log(@f(g_hcb))
end

function g_hcb(model::BACKSAFTModel, V, T, z)
    α = model.params.alpha.values[1]
    η = @f(ζ,3)
    return 1/(1-η)+3*(1+α)*α*η/((1-η)^2*(1+3α))+3*η^2*α^2/((1-η)^3*(1+3α))
end

const BACKSAFT_consts = (
    D1 = [-8.8043,4.164627,-48.203555,140.4362,-195.23339,113.515],
    D2 = [2.9396,-6.0865383,40.137956,-76.230797,-133.70055,860.25349,-1535.3224,1221.4261,-409.10539],
    D3 = [-2.8225,4.7600148,11.257177,-66.382743,69.248785],
    D4 = [0.34,-3.1875014,12.231796,-12.110681],
)


text:
module SailingNavigation

export Position, lat, lon, GridPoint, TimeSlice, TimedPath, closestpoint, surround
export RoutingProblem, minimumtimeroute

using GeometryTypes
using SailingPolars

# NB: This uses latitude (often considered to be y) first then longitude (often considered to be x).
# This latitude, then longitude ordering is as per ISO 6709 (en.wikipedia.org/wiki/ISO_6709)

# Position is a Float32 2-tuple of latitude and longitude in degrees
Position = Point2f0

# latitude from Position
lat(p::Position) = p[1]

# longitude from Position
lon(p::Position) = p[2]

# A GridPoint is a Position with the SurfaceParameters of wind and current at the Position
mutable struct GridPoint
    pt::Position
    sp::SurfaceParameters
end

"""
    TimeSlice

A TimeSlice is a matrix of GridPoints, each Position point with their SurfaceParameters
A Vector of TimeSlice can give the surface characteristics for an ocean region over time.
"""
TimeSlice = Matrix{GridPoint}

"""
    mutable struct RoutingProblem

timeinterval: the seconds duration for each TimeSlice
timeframe: a vector of sequential timeslices for the ocean region
obstacleindices: the Cartesian indices in each timeslice for
    obstacles, such as land or shoals, where the ship may not go
startindex: the timeslice position for time of starting
start: starting location on grid of GridPoints
finish: destination / finish location on grid of GridPoints
allowrepeatvisits: whether the vessel may overlap its prior path, usually false
"""
mutable struct RoutingProblem
    timeinterval::Float64 # minutes between timeframe slices
    timeframe::Vector{TimeSlice}
    obstacleindices::Vector{Point2{Int}}
    startindex::Int
    start::Point2{Int}
    finish::Point2{Int}
    allowrepeatvisits::Bool
end

"""
    mutable struct TimedPath

duration: minutes total to travel the path
path: vector of Cartesian indices of points in grid for path to travel
"""
mutable struct TimedPath
    duration::Float64
    path::Vector{Point2{Int}}
end

"""
    closestpoint(p, mat)

Get the closest GridPoint in matrix mat to a given position p.
p: Cartesian indices of a Position (latitude, longitude in degrees) in grid of GridPoints
mat: matrix of Gridpoints
"""
closestpoint(p, mat) = findmin(gp -> haversine(p[1], p[2], gp.pt[1], gp.pt[2])[1], mat)[2]

function surround(p, mat, excluded)
    neighbors = Point2{Int}[(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]
    (xmax, ymax) = size(mat)
    return filter(q -> 1 <= q[1] <= xmax && 1 <= q[2] <= ymax && !(q in excluded),
        [x + p for x in neighbors])
end

"""
    function minimumtimeroute(rp::RoutingProblem, sp::SailingPolar, verbose=false)

Get the route (as a TimedPath) that minimizes time from start to finish for a given
RoutingProblem (sea parameters) given sailing polar data (ship parameters).
"""
function minimumtimeroute(rp::RoutingProblem, sp::SailingPolar, verbose=false)
    timedpaths = [TimedPath(0.0, [rp.start])]
    completed, mintime, minpath = false, 1000.0, TimedPath(1000.0, [])
    for i in 1:1000
        newpaths = TimedPath[]
        verbose && println("Checking ", length(timedpaths), " paths of length ",
            length(timedpaths[1].path) - 1)
        for tpath in timedpaths
            if tpath.path[end] == rp.finish
                completed = true
                push!(newpaths, tpath)
            else
                p1 = tpath.path[end]
                slice = rp.timeframe[div(Int(round(tpath.duration)),
                                     Int(round(rp.timeinterval))) %
                                     length(rp.timeframe) + 1]
                for p2 in surround([p1[1], p1[2]], slice, rp.obstacleindices)
                    !rp.allowrepeatvisits && p2 in tpath.path && continue
                    gp1, gp2 = slice[p1[1], p1[2]], slice[p2[1], p2[2]]
                    lat1, lon1, lat2, lon2 = gp1.pt[1], gp1.pt[2], gp2.pt[1], gp2.pt[2]
                    t = sailsegmenttime(sp, gp1.sp, lat1, lon1, lat2, lon2)
                    path = deepcopy(tpath.path)
                    push!(path, p2)
                    push!(newpaths, TimedPath(tpath.duration + t, path))
                end
            end
        end
        timedpaths = unique(newpaths)
        if completed
            mindur = minimum(map(x -> x.duration, timedpaths))
            finished = filter(x -> x.path[end] == rp.finish, timedpaths)
            minfindur, idx = findmin(map(x -> x.duration, finished))
            verbose && println("Current finished minimum: ", finished[idx], ", others $mindur")
            if mindur == minfindur
                minpath = finished[idx]
                break
            end
        end
    end
    return minpath
end

end # module


text:
@testset "Starting Area" begin

    # Test Starting Area
    x_nom = VehicleState(1.0, 1.0, 1.0, 1.0)
    x_min = VehicleState(0.0, 0.0, 0.0, 0.0)
    x_max = VehicleState(1.0, 2.0, 3.0, 4.0)
    start1 = StartingArea(x_nom, x_min, x_max)
    start2 = StartingArea(x_nom)
    start3 = StartingArea(x_min, x_max)

    @test start1.x_nom == x_nom
    @test start1.x_min == x_min
    @test start1.x_max == x_max

    @test start2.x_nom == x_nom
    @test start2.x_min == x_nom
    @test start2.x_max == x_nom

    @test start3.x_nom == VehicleState(0.5, 1.0, 1.5, 2.0)
    @test start3.x_min == x_min
    @test start3.x_max == x_max

    xr = randstate(start1)
    @test 0.0 <= xr.x <= 1.0
    @test 0.0 <= xr.y <= 2.0
    @test 0.0 <= xr.θ <= 3.0
    @test 0.0 <= xr.v <= 4.0

end


text:
function k=JacVec(J,v)
n1=size(v,1);
n2=size(v,2);
n=n1*n2;
rRho=reshape(permute(v(:,:,1),[2,1]),n,1);
rTh=reshape(permute(v(:,:,5),[2,1]),n,1);
rw=reshape(permute(v(:,:,4),[2,1]),n,1);
sRho=J.JRhoW*rw;
sTh=J.JThW*rw;
sw=J.JWRho*rRho+J.JWTh*rTh;
k=zeros(size(v));
k(:,:,1)=permute(reshape(sRho,n2,n1),[2,1]);
k(:,:,4)=permute(reshape(sw,n2,n1),[2,1]);
k(:,:,5)=permute(reshape(sTh,n2,n1),[2,1]);
end


text:
# assumes Choi matrix
apply_process(C::Operator,ψ::Ket) = apply_process(C,Operator(ψ))

apply_process(C::Operator,ρ::Operator{T,D}) where {T,D} =
    ptrace((transpose(ρ)⊗qeye(dims(ρ)))*C,ntuple(identity,D))

# assumes Kraus operators
apply_process(As::Vector{<:Operator},ρ::Operator) = sum(A->A*ρ*A',As)

function gate_fidelity_choi(C::Operator,U::Operator)
    # calculate the average gate fidelity given a Choi matrix and a unitary gate
    dims(C) != (dims(U)...,dims(U)...) && throw(DimensionMismatch())
    V = qeye(dims(U)) ⊗ U
    return gate_fidelity_choi(V'*C*V)
end

function gate_fidelity_choi(C::Operator)
    # calculate the average gate fidelity given a Choi matrix
    # Johnston, N. & Kribs, D. W. Quantum gate fidelity in terms of Choi
    # matrices. J. Phys. A Math. Theor. 44, 495303 (2011).
    d = check_choi(C)
    TrA² = let CC = Operator(data(C),(d,d)) # rewrap operator to use tensored indexing
        real(sum(CC[(i,i),(j,j)] for i=0:d-1, j=0:d-1))
    end
    return (d + TrA²)/(d^2 + d)
end

function gate_fidelity_kraus(As::Vector{Operator{T,D}}) where {T,D}
    # calculate the average gate fidelity given a list of Kraus operators
    d = prod(i->dims(first(As))[i],1:D)
    return (d + sum(abs2∘trace,As))/(d^2 + d)
end

function gate_fidelity_kraus(As::Vector{<:Operator{T1,D}},U::Operator{T2,D}) where {T1,T2,D}
    # calculate the average gate fidelity given a list of Kraus operators and
    # a unitary gate
    dimsmatch(As,U)
    d = prod(i->dims(U)[i],1:D)
    return (d + sum(abs2∘inner,product(As,(U,))))/(d^2 + d)
end

# conversion functions
"""
    operator_to_choi(O)

Compute the Choi matrix ``C_O`` of a unitary map represented by an operator `O`. We use the
convention of applying `O` to the second half of the extended Hilbert space, i.e.:
```math
C_O = ∑_{i,j} |i⟩⟨j| ⊗ O|i⟩⟨j|O^†
```
"""
operator_to_choi(O::Operator) = Operator(vec(O))

kraus_to_choi(As::Vector{<:Operator}) = sum(operator_to_choi,As)

kraus_to_natural(As::Vector{<:Operator}) = sum(A->conj(A)⊗A,As)

# go through Choi representation; the extra step is cheap compared to the eigendecomposition
natural_to_kraus(K::Operator) = choi_to_kraus(natural_to_choi(K))

# the conversion between the natural and Choi is a self-inverse
natural_to_choi(K::Operator) = choi_to_natural(K)

function choi_to_natural(C::Operator)
    # the "natural" representation is the same as Choi with shuffled dimensions
    d = check_choi(C)
    K = reshape(permutedims(reshape(data(C),(d,d,d,d)),[1,3,2,4]),(d^2,d^2))
    return Operator(K,dims(C))
end

function choi_to_kraus(C::Operator)
    check_choi(C)
    D,V = eigen(Hermitian(data(C)),1E-10,Inf)
    dm = _half_dims(dims(C))
    return [Operator(unvec(√(D[i])*V[:,i]),dm) for i = length(D):-1:1]
end

function choi_to_chi(C::Operator{T,D}) where {T,D}
    B = Operator(_choi_to_pauli_basis(_half_dims(dims(C))),dims(C))
    return (B*C*B')/2^(D÷2)
end

function chi_to_choi(χ::Operator{T,D}) where {T,D}
    B = Operator(_choi_to_pauli_basis(_half_dims(dims(χ))),dims(χ))
    return (B'*χ*B)/2^(D÷2)
end

function _choi_to_pauli_basis(D::Dims{N}) where N
    all(d->d==2,D) || throw(ArgumentError("only valid with qubits!"))
    B = Matrix{Float64}(undef,4^N,4^N)
    # change to I,X,Y,Z basis (Y = -i*σy)
    Ps = ([1 0; 0 1], [0 1; 1 0], [0 -1; 1 0], [1 0; 0 -1])
    for (i,ops) in enumerate(product(ntuple(_->Ps,N)...))
        B[:,i] = vec(reduce(⊗,ops))
    end
    return B
end

Base.vec(O::Operator) = Ket(vec(data(O)),(dims(O)...,dims(O)...))

unvec(v::Ket) = Operator(unvec(data(v)),_half_dims(dims(v)))

super(A::Operator) = conj(A) ⊗ A # represents A*ρ*A†
super(A::Operator,B::Operator) = transpose(B) ⊗ A # represents A*ρ*B
super(A::Operator,::UniformScaling) = qeye(dims(A)) ⊗ A # represents A*ρ
super(::UniformScaling,B::Operator) = transpose(B) ⊗ qeye(dims(B)) # represents ρ*B

# linear operator that represents the action AρB on a vectorized version of ρ
super(A::AbstractMatrix,B::AbstractMatrix=data(qeye(size(A,1)))) = copy(transpose(B)) ⊗ A

function unvec(vecA::AbstractVector)
    # unvectorize a vector into a square matrix
    d = isqrt(length(vecA))
    return reshape(vecA,(d,d))
end

function check_choi(C::Operator{T,D}) where {T,D}
    dm = dims(C)
    if D%2 != 0 || prod(i->dm[i],1:D÷2) != prod(i->dm[i],D÷2+1:D)
        throw(ArgumentError("not a valid Choi matrix!"))
    end
    return prod(i->dm[i],1:D÷2)
end

_half_dims(dm::Dims{D}) where{D} = ntuple(i->dm[i],D÷2)


text:
module TestCallbacks

using Gurobi
using Random
using Test

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

const MOI = Gurobi.MOI

const GRB_ENV =
    isdefined(Main, :GRB_ENV) ? Main.GRB_ENV : Gurobi.Env(output_flag = 0)

function callback_simple_model()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.RawOptimizerAttribute("OutputFlag"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("Cuts"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("Presolve"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("Heuristics"), 0)
    MOI.Utilities.loadfromstring!(
        model,
        """
    variables: x, y
    maxobjective: y
    c1: x in Integer()
    c2: y in Integer()
    c3: x in Interval(0.0, 2.5)
    c4: y in Interval(0.0, 2.5)
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    return model, x, y
end

function callback_knapsack_model()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.RawOptimizerAttribute("OutputFlag"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("Cuts"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("Presolve"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("PreCrush"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("Heuristics"), 0)
    N = 30
    x = MOI.add_variables(model, N)
    MOI.add_constraints(model, x, MOI.ZeroOne())
    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
    Random.seed!(1)
    item_weights, item_values = rand(N), rand(N)
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_weights, x), 0.0),
        MOI.LessThan(10.0),
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_values, x), 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    return model, x, item_weights
end

function test_lazy_constraint_callback()
    model, x, y = callback_simple_model()
    lazy_called = false
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> begin
            lazy_called = true
            x_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), x)
            y_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), y)
            status = MOI.get(
                model,
                MOI.CallbackNodeStatus(cb_data),
            )::MOI.CallbackNodeStatusCode
            int_sol = round.(Int, [x_val, y_val])
            if status == MOI.CALLBACK_NODE_STATUS_INTEGER
                @test isapprox(int_sol, [x_val, y_val], atol = 1e-6)
            end
            @test MOI.supports(model, MOI.LazyConstraint(cb_data))
            if y_val - x_val > 1 + 1e-6
                @test MOI.submit(
                    model,
                    MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([-1.0, 1.0], [x, y]),
                        0.0,
                    ),
                    MOI.LessThan{Float64}(1.0),
                ) === nothing
            elseif y_val + x_val > 3 + 1e-6
                @test MOI.submit(
                    model,
                    MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                        0.0,
                    ),
                    MOI.LessThan{Float64}(3.0),
                ) === nothing
            end
        end,
    )
    @test MOI.supports(model, MOI.LazyConstraintCallback())
    MOI.optimize!(model)
    @test lazy_called
    @test MOI.get(model, MOI.VariablePrimal(), x) == 1
    @test MOI.get(model, MOI.VariablePrimal(), y) == 2
end

function test_lazy_constraint_callback_fractional()
    # callback_simple_model() is not large enough to see MIPNODE callbacks.
    model, x, item_weights = callback_knapsack_model()
    lazy_called_integer = false
    lazy_called_fractional = false
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> begin
            status = MOI.get(
                model,
                MOI.CallbackNodeStatus(cb_data),
            )::MOI.CallbackNodeStatusCode
            if status == MOI.CALLBACK_NODE_STATUS_INTEGER
                lazy_called_integer = true
            elseif status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL
                lazy_called_fractional = true
            end
        end,
    )
    @test MOI.supports(model, MOI.LazyConstraintCallback())
    MOI.optimize!(model)
    @test lazy_called_integer
    @test lazy_called_fractional
end

function test_lazy_constraint_callback_OptimizeInProgress()
    model, x, y = callback_simple_model()
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> begin
            @test_throws(
                MOI.OptimizeInProgress(MOI.VariablePrimal()),
                MOI.get(model, MOI.VariablePrimal(), x)
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveValue()),
                MOI.get(model, MOI.ObjectiveValue())
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveBound()),
                MOI.get(model, MOI.ObjectiveBound())
            )
        end,
    )
    return MOI.optimize!(model)
end

function test_lazy_constraint_callback_UserCut()
    model, x, y = callback_simple_model()
    cb = nothing
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.UserCut(cb_data),
                MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
                MOI.LessThan(2.0),
            )
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(MOI.LazyConstraintCallback(), MOI.UserCut(cb)),
        MOI.optimize!(model)
    )
end

function test_lazy_constraint_callback_HeuristicSolution()
    model, x, y = callback_simple_model()
    cb = nothing
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(model, MOI.HeuristicSolution(cb_data), [x], [2.0])
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(
            MOI.LazyConstraintCallback(),
            MOI.HeuristicSolution(cb),
        ),
        MOI.optimize!(model)
    )
end

function test_user_cut_callback()
    model, x, item_weights = callback_knapsack_model()
    user_cut_submitted = false
    MOI.set(
        model,
        MOI.UserCutCallback(),
        cb_data -> begin
            terms = MOI.ScalarAffineTerm{Float64}[]
            accumulated = 0.0
            for (i, xi) in enumerate(x)
                if MOI.get(model, MOI.CallbackVariablePrimal(cb_data), xi) > 0.0
                    push!(terms, MOI.ScalarAffineTerm(1.0, xi))
                    accumulated += item_weights[i]
                end
            end
            @test MOI.supports(model, MOI.UserCut(cb_data))
            if accumulated > 10.0
                @test MOI.submit(
                    model,
                    MOI.UserCut(cb_data),
                    MOI.ScalarAffineFunction{Float64}(terms, 0.0),
                    MOI.LessThan{Float64}(length(terms) - 1),
                ) === nothing
                user_cut_submitted = true
            end
        end,
    )
    @test MOI.supports(model, MOI.UserCutCallback())
    MOI.optimize!(model)
    @test user_cut_submitted
end

function test_user_cut_callback_LazyConstraint()
    model, x, item_weights = callback_knapsack_model()
    cb = nothing
    MOI.set(
        model,
        MOI.UserCutCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.LazyConstraint(cb_data),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
                MOI.LessThan(5.0),
            )
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(MOI.UserCutCallback(), MOI.LazyConstraint(cb)),
        MOI.optimize!(model)
    )
end

function test_user_cut_callback_HeuristicSolution()
    model, x, item_weights = callback_knapsack_model()
    cb = nothing
    MOI.set(
        model,
        MOI.UserCutCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(model, MOI.HeuristicSolution(cb_data), [x[1]], [0.0])
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(
            MOI.UserCutCallback(),
            MOI.HeuristicSolution(cb),
        ),
        MOI.optimize!(model)
    )
end

function test_heuristic_callback()
    model, x, item_weights = callback_knapsack_model()
    solution_accepted = false
    solution_rejected = false
    MOI.set(
        model,
        MOI.HeuristicCallback(),
        cb_data -> begin
            x_vals = MOI.get.(model, MOI.CallbackVariablePrimal(cb_data), x)
            @test MOI.supports(model, MOI.HeuristicSolution(cb_data))
            status = MOI.get(
                model,
                MOI.CallbackNodeStatus(cb_data),
            )::MOI.CallbackNodeStatusCode
            int_sol = round.(Int, x_vals)
            if status == MOI.CALLBACK_NODE_STATUS_INTEGER
                @test isapprox(int_sol, x_vals, atol = 1e-6)
            end
            if MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                x,
                floor.(x_vals),
            ) == MOI.HEURISTIC_SOLUTION_ACCEPTED
                solution_accepted = true
            end
            if MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                x,
                ceil.(x_vals),
            ) == MOI.HEURISTIC_SOLUTION_REJECTED
                solution_rejected = true
            end
        end,
    )
    @test MOI.supports(model, MOI.HeuristicCallback())
    MOI.optimize!(model)
    @test solution_accepted
    @test solution_rejected
end

function test_heuristic_callback_LazyConstraint()
    model, x, item_weights = callback_knapsack_model()
    cb = nothing
    MOI.set(
        model,
        MOI.HeuristicCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.LazyConstraint(cb_data),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
                MOI.LessThan(5.0),
            )
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(
            MOI.HeuristicCallback(),
            MOI.LazyConstraint(cb),
        ),
        MOI.optimize!(model)
    )
end

function test_heuristic_callback_UserCut()
    model, x, item_weights = callback_knapsack_model()
    cb = nothing
    MOI.set(
        model,
        MOI.HeuristicCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.UserCut(cb_data),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
                MOI.LessThan(5.0),
            )
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(MOI.HeuristicCallback(), MOI.UserCut(cb)),
        MOI.optimize!(model)
    )
end

function test_CallbackFunction_callback_OptimizeInProgress()
    model, x, y = callback_simple_model()
    MOI.set(
        model,
        Gurobi.CallbackFunction(),
        (cb_data, cb_where) -> begin
            @test_throws(
                MOI.OptimizeInProgress(MOI.VariablePrimal()),
                MOI.get(model, MOI.VariablePrimal(), x)
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveValue()),
                MOI.get(model, MOI.ObjectiveValue())
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveBound()),
                MOI.get(model, MOI.ObjectiveBound())
            )
        end,
    )
    @test MOI.supports(model, Gurobi.CallbackFunction())
    return MOI.optimize!(model)
end

function test_CallbackFunction_callback_LazyConstraint()
    model, x, y = callback_simple_model()
    cb_calls = Int32[]
    function callback_function(cb_data::Gurobi.CallbackData, cb_where::Cint)
        push!(cb_calls, cb_where)
        if cb_where == Gurobi.GRB_CB_MIPSOL
            Gurobi.load_callback_variable_primal(cb_data, cb_where)
            x_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), x)
            y_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), y)
            if y_val - x_val > 1 + 1e-6
                MOI.submit(
                    model,
                    MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([-1.0, 1.0], [x, y]),
                        0.0,
                    ),
                    MOI.LessThan{Float64}(1.0),
                )
            elseif y_val + x_val > 3 + 1e-6
                MOI.submit(
                    model,
                    MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                        0.0,
                    ),
                    MOI.LessThan{Float64}(3.0),
                )
            end
        end
    end
    MOI.set(model, MOI.RawOptimizerAttribute("LazyConstraints"), 1)
    MOI.set(model, Gurobi.CallbackFunction(), callback_function)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), x) == 1
    @test MOI.get(model, MOI.VariablePrimal(), y) == 2
    @test length(cb_calls) > 0
    @test Gurobi.GRB_CB_MESSAGE in cb_calls
    @test Gurobi.GRB_CB_PRESOLVE in cb_calls
    @test Gurobi.GRB_CB_MIPSOL in cb_calls
end

function test_CallbackFunction_callback_UserCut()
    model, x, item_weights = callback_knapsack_model()
    user_cut_submitted = false
    cb_calls = Int32[]
    MOI.set(
        model,
        Gurobi.CallbackFunction(),
        (cb_data, cb_where) -> begin
            push!(cb_calls, cb_where)
            if cb_where != Gurobi.GRB_CB_MIPNODE
                return
            end
            status = Ref{Cint}()
            Gurobi.GRBcbget(
                cb_data,
                cb_where,
                Gurobi.GRB_CB_MIPNODE_STATUS,
                status,
            )
            if status[] != 2
                return  # Not optimal.
            end
            Gurobi.load_callback_variable_primal(cb_data, cb_where)
            terms = MOI.ScalarAffineTerm{Float64}[]
            accumulated = 0.0
            for (i, xi) in enumerate(x)
                if MOI.get(model, MOI.CallbackVariablePrimal(cb_data), xi) > 0.0
                    push!(terms, MOI.ScalarAffineTerm(1.0, xi))
                    accumulated += item_weights[i]
                end
            end
            if accumulated > 10.0
                MOI.submit(
                    model,
                    MOI.UserCut(cb_data),
                    MOI.ScalarAffineFunction{Float64}(terms, 0.0),
                    MOI.LessThan{Float64}(length(terms) - 1),
                )
                user_cut_submitted = true
            end
        end,
    )
    MOI.optimize!(model)
    @test user_cut_submitted
    @test Gurobi.GRB_CB_MIPNODE in cb_calls
end

function test_CallbackFunction_callback_HeuristicSolution()
    model, x, item_weights = callback_knapsack_model()
    solution_accepted = false
    solution_rejected = false
    cb_calls = Int32[]
    MOI.set(
        model,
        Gurobi.CallbackFunction(),
        (cb_data, cb_where) -> begin
            push!(cb_calls, cb_where)
            if cb_where != Gurobi.GRB_CB_MIPNODE
                return
            end
            status = Ref{Cint}()
            Gurobi.GRBcbget(
                cb_data,
                cb_where,
                Gurobi.GRB_CB_MIPNODE_STATUS,
                status,
            )
            if status[] != 2
                return  # Not optimal.
            end
            Gurobi.load_callback_variable_primal(cb_data, cb_where)
            x_vals =
                MOI.get.(model, MOI.CallbackVariablePrimal(cb_data), x)
            if MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                x,
                floor.(x_vals),
            ) == MOI.HEURISTIC_SOLUTION_ACCEPTED
                solution_accepted = true
            end
            if MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                x,
                ceil.(x_vals),
            ) == MOI.HEURISTIC_SOLUTION_REJECTED
                solution_rejected = true
            end
        end,
    )
    MOI.optimize!(model)
    @test solution_accepted
    @test solution_rejected
    @test Gurobi.GRB_CB_MIPNODE in cb_calls
end

function test_CallbackFunction_CallbackNodeStatus()
    model, x, item_weights = callback_knapsack_model()
    unknown_reached = false
    MOI.set(
        model,
        Gurobi.CallbackFunction(),
        (cb_data, cb_where) -> begin
            if MOI.get(model, MOI.CallbackNodeStatus(cb_data)) ==
               MOI.CALLBACK_NODE_STATUS_UNKNOWN
                unknown_reached = true
            end
        end,
    )
    MOI.optimize!(model)
    @test unknown_reached
end

function test_CallbackFunction_broadcast()
    model, x, _ = callback_knapsack_model()
    f(cb_data, x) = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), x)
    solutions = Vector{Float64}[]
    MOI.set(
        model,
        Gurobi.CallbackFunction(),
        (cb_data, cb_where) -> begin
            if cb_where == Gurobi.GRB_CB_MIPSOL
                Gurobi.load_callback_variable_primal(cb_data, cb_where)
                push!(solutions, f.(cb_data, x))
            end
        end,
    )
    MOI.optimize!(model)
    @test length(solutions) > 0
    @test length(solutions[1]) == length(x)
end

end  # module TestCallbacks

TestCallbacks.runtests()


text:
module SixelGraphics

export Screen, Bitmap, sixelplot, clear, putpixel, sixels

const NCOL = 16

# structs

struct Defaults
  ax1::Float64
  axy::Float64
  axy2::Float64
  ay1::Float64
  ayx::Float64
  tl::Float64
  function Defaults()
    new(0.2,0.1,0.9,0.15,0.15,0.01)
  end
end

mutable struct Scaler
  a::Float64
  b::Float64
  mn::Float64
  mx::Float64
  function Scaler(omn,omx,nmn,nmx)
    b = (omx-omn)/(nmx-nmn)
    a = omn - b*nmn
    new(a,b,nmn,nmx)
  end
end

mutable struct Color
  d::Bool
  r::Int
  g::Int
  b::Int
  function Color()
    new(false,0,0,0)
  end
end

struct Palette
  colors::Vector{Color}
  function Palette()
    new([Color() for i=1:NCOL])
  end
end

struct Font
  fnt::Vector{UInt8}
  function Font()
    fid = open(joinpath(dirname(pathof(SixelGraphics)),"Lat7-VGA8.raw"),"r")
    fnt = read(fid)
    close(fid)
    new(fnt)
  end
end

struct Bitmap
  nr::Int
  nc::Int
  planes::Array{UInt8,3}
  function Bitmap(nr::Int,nc::Int)
    new(nr,nc,zeros(UInt8, NCOL, nr, nc))
  end
end

struct Screen
  d::Defaults
  b::Bitmap
  sx::Scaler
  sy::Scaler
  p::Palette
  f::Font
  function Screen(nx::Int,ny::Int)
    new(Defaults(),Bitmap(div(ny,6),nx),Scaler(0,nx-1,0,1),Scaler(0,ny-1,0,1),Palette(),Font())
  end
end

# Scaler methods

function scale(s::Scaler, x)
  if x<s.mn
    return Int(floor(s.mn))
  elseif x>s.mx
    return Int(floor(s.mx))
  end
  Int(floor(s.a+s.b*x))
end

function set(s::Scaler,omn,omx,nmn,nmx)
  s.b = (omx-omn)/(nmx-nmn)
  s.a = omn - s.b*nmn
  s.mn = nmn
  s.mx = nmx
end

# Palette methods

function setcolor(p::Palette, c, r, g, b)
  p.colors[c].d = true
  p.colors[c].r = r
  p.colors[c].g = g
  p.colors[c].b = b
end

# Bitmap methods

function clear(b::Bitmap, v)
  b.planes[:,:,:] .= v
end

function putpixel(b::Bitmap, ic::Int, ix::Int, iy::Int)

  if ic<1 || ic>NCOL || ix<0 || ix>=b.nc || iy<0 || iy>=b.nr*6; return; end

  ix1 = ix+1
  ib1 = div(iy,6)
  ib2 = iy - ib1*6
  ib1 += 1
  ib3 = UInt8(32 >> ib2)

  b.planes[ic,ib1,ix1] = 63 + ((b.planes[ic,ib1,ix1] - 63) | ib3)
end

function idrawline(b::Bitmap, ic::Int, ix1::Int, iy1::Int, ix2::Int, iy2::Int)
  dx = ix2-ix1
  dy = iy2-iy1
  adx = abs(dx)
  ady = abs(dy)

  n = adx>ady ? adx : ady

  ax = dx/n
  ay = dy/n

  for i = 0:n
    ix = Int(round(ix1 + i*ax))
    iy = Int(round(iy1 + i*ay))
    putpixel(b,ic,ix,iy)
  end
end

function idrawchar(b::Bitmap, f::Font, ic::Int, ix::Int, iy::Int, c::Char)
  msk = UInt8(0)

  ptr = Int(c)*8+1
  iyj = iy+7

  for j = 0:7
    msk = 0x80
    ixi = ix
    for i = 0:7
      if (f.fnt[ptr] & msk) > 0; putpixel(b,ic,ixi,iyj); end
      msk >>= 1
      ixi += 1
    end
    ptr += 1
    iyj -= 1
  end
end

function idrawlchar(b::Bitmap, f::Font, ic::Int, ix::Int, iy::Int, c::Char)
  msk = UInt8(0)

  ptr = Int(c)*8+1
  ixj = ix

  for j = 0:7
    msk = 0x80
    iyi = iy
    for i = 0:7
      if (f.fnt[ptr] & msk) > 0; putpixel(b,ic,ixj,iyi); end
      msk >>= 1
      iyi += 1
    end
    ptr += 1
    ixj += 1
  end
end

function idrawstring(b::Bitmap, f::Font, ic::Int, ix::Int, iy::Int, a::String)
  ix1 = ix
  for i=1:length(a); idrawchar(b, f, ic, ix1, iy, a[i]); ix1+=8; end
end

function idrawlstring(b::Bitmap, f::Font, ic::Int, ix::Int, iy::Int, a::String)
  iy1 = iy
  for i=1:length(a); idrawlchar(b, f, ic, ix, iy1, a[i]); iy1+=8; end
end

function idrawmarker(b::Bitmap, f::Font, ic::Int, ix::Int, iy::Int, m::Int)
  if m==0
    putpixel(b, ic, ix,iy)
  elseif m==1
    idrawline(b, ic, ix-2,iy, ix+2, iy)
    idrawline(b, ic, ix,iy-2, ix, iy+2)
  elseif m==2
    idrawline(b, ic, ix-2,iy-2, ix+2, iy+2)
    idrawline(b, ic, ix-2,iy+2, ix+2, iy-2)
  elseif m==3
    idrawline(b, ic, ix-2,iy-2, ix-2, iy+2)
    idrawline(b, ic, ix-2,iy+2, ix+2, iy+2)
    idrawline(b, ic, ix+2,iy+2, ix+2, iy-2)
    idrawline(b, ic, ix+2,iy-2, ix-2, iy-2)
  else
    idrawchar(b, f, ic, ix-4, iy-4, Char(m))
  end
end

# Screen methods

function drawline(s::Screen, ic::Int, x1, y1, x2, y2)
  idrawline(s.b, ic, scale(s.sx,x1), scale(s.sy,y1), scale(s.sx,x2), scale(s.sy,y2))
end

function drawstring(s::Screen, ic::Int, x, y, a::String, or::Char)
  ix = scale(s.sx,x)
  iy = scale(s.sy,y)
  n = length(a)

  if or == 'c'
    idrawstring(s.b, s.f, ic, ix-n*4, iy-4, a)
  elseif or == 'u'
    idrawstring(s.b, s.f, ic, ix-n*4, iy, a)
  elseif or == 'd'
    idrawstring(s.b, s.f, ic, ix-n*4, iy-8, a)
  elseif or == 'r'
    idrawstring(s.b, s.f, ic, ix, iy-4, a)
  elseif or == 'l'
    idrawstring(s.b, s.f, ic, ix-n*8, iy-4, a)
  end
end

function drawlstring(s::Screen, ic::Int, x, y, a::String, or::Char)
  ix = scale(s.sx,x)
  iy = scale(s.sy,y)
  n = length(a)

  if or == 'r'
    idrawlstring(s.b, s.f, ic, ix, iy-n*4, a)
  elseif or == 'l'
    idrawlstring(s.b, s.f, ic, ix-8, iy-n*4, a)
  end
end

function drawplot(s::Screen, ic::Int, x::Array, y::Array, typ::Char, pch::Int)
  dl = typ=='l' || typ=='b'
  dp = typ=='p' || typ=='b'

  n = length(x)

  for i=2:n
    ix1 = scale(s.sx,x[i-1])
    iy1 = scale(s.sy,y[i-1])
    ix2 = scale(s.sx,x[i])
    iy2 = scale(s.sy,y[i])
    if dl; idrawline(s.b,ic,ix1,iy1,ix2,iy2); end
    if pch>=0
      if i==2 && dp; idrawmarker(s.b,s.f,ic,ix1,iy1,pch); end
      if dp; idrawmarker(s.b,s.f,ic,ix2,iy2,pch); end
    else
      if i==2 && dp; idrawmarker(s.b,s.f,1,ix1,iy1,96+i-1); end
      if dp; idrawmarker(s.b,s.f,1,ix2,iy2,96+i); end
    end
  end
end

function adjust(mn, mx)
  d = exp(log(10)*Int(floor(log10(mx-mn))))
  mn = floor(mn/d)
  mx = ceil(mx/d)
  nt = Int(round(mx-mn))
  if d>=1
    [Int(round(mn*d)) Int(round(mx*d)) nt]
  else
    [mn*d mx*d nt]
  end
end

function axes(s::Screen, xlab::String, autox::Bool, xmin, xmax,
                         ylab::String, autoy::Bool, ymin, ymax)
  drawline(s, 1, s.d.ax1, s.d.axy, s.d.axy2, s.d.axy)
  drawline(s, 1, s.d.ax1, s.d.axy, 0.2, s.d.axy+s.d.tl)
  drawline(s, 1, s.d.axy2, s.d.axy, s.d.axy2, s.d.axy+s.d.tl)

  drawline(s, 1, s.d.ayx, s.d.ay1, s.d.ayx, s.d.axy2)
  drawline(s, 1, s.d.ayx, s.d.ay1, s.d.ayx+s.d.tl, 0.15)
  drawline(s, 1, s.d.ayx, s.d.axy2, s.d.ayx+s.d.tl, s.d.axy2)

  drawstring(s, 1, (s.d.ax1+s.d.axy2)/2., s.d.axy-2*s.d.tl, xlab, 'd')
  drawlstring(s, 1, s.d.ayx-2*s.d.tl, (s.d.ay1+s.d.axy2)/2., ylab, 'l')

  if autox
    xp = adjust(xmin,xmax)
  else
    xp = [xmin xmax 1]
  end
  if autoy
    yp = adjust(ymin,ymax)
  else
    yp = [ymin ymax 1]
  end

  ntx = xp[3]
  nty = yp[3]

  for i = 0:ntx
    d = s.d.ax1 + (s.d.axy2-s.d.ax1)*i/ntx
    drawline(s, 1, d, s.d.axy, d, s.d.axy-s.d.tl)
  end

  for i = 0:nty
    d = s.d.ay1 + (s.d.axy2-s.d.ay1)*i/nty
    drawline(s, 1, s.d.ayx, d, s.d.ayx-s.d.tl, d)
  end

  drawstring(s, 1, s.d.ax1,s.d.axy-2*s.d.tl, repr(xp[1]), 'd')
  drawstring(s, 1, s.d.axy2,s.d.axy-2*s.d.tl, repr(xp[2]), 'd')

  drawstring(s, 1, s.d.ayx-2*s.d.tl,s.d.ay1, repr(yp[1]), 'l')
  drawstring(s, 1, s.d.ayx-2*s.d.tl,s.d.axy2, repr(yp[2]), 'l')

  set(s.sx, scale(s.sx,s.d.ax1), scale(s.sx,s.d.axy2), xp[1],xp[2])
  set(s.sy, scale(s.sy,s.d.ay1), scale(s.sy,s.d.axy2), yp[1],yp[2])
end

function frame(s::Screen)
  drawline(s, 1, 0.,0., 1.,0.)
  drawline(s, 1, 0.,1., 1.,1.)
  drawline(s, 1, 0.,0., 0.,1.)
  drawline(s, 1, 1.,0., 1.,1.)
end

# Sixel functions

function csixels(b::Bitmap, eol::Bool, lc::Int, ic::Int, i::Int)
  dc = true
  k = 1

  while k<=b.nc
    c = b.planes[ic,i,k]

    j = k + 1
    while j<=b.nc
      if b.planes[ic,i,j] != c; break; end
      j += 1
    end

    n = j-k
    if n==b.nc && c=='?'; return true; end

    if eol; println('$'); eol=false; end

    if dc
      if ic!=lc
        print('#',ic-1)
        lc = ic
      end
      dc = false
    end

    if n>3
      print('!',n,Char(c))
      k = j
    else
      print(Char(c))
      k += 1
    end
  end

  false
end

function sixels(s::Screen)
  lc = -1

  println("\n\x1bPq\n");

  for i = 1:NCOL
    if s.p.colors[i].d
      print('#',i-1,';',2,';',s.p.colors[i].r,';',s.p.colors[i].g,';',s.p.colors[i].b);
    end
  end

  for i = s.b.nr:-1:1
    eol = false
    for ic = 1:NCOL
      empty = csixels(s.b,eol,lc,ic,i);
      if !empty; eol = true; end
      if ic==NCOL; println('-'); end
    end
  end

  println("\x1b\\\n");
end

function julialogo(s::Screen)
  ix = scale(s.sx,0.9)-8*21
  iy = 100

  idrawstring(s.b, s.f, 4, ix,iy, "               _")
  iy-=8

  idrawstring(s.b, s.f, 2, ix,iy, "   _")
  idrawstring(s.b, s.f, 1, ix,iy, "           _")
  idrawstring(s.b, s.f, 4, ix,iy, "              (_) ")
  idrawstring(s.b, s.f, 3, ix,iy, "             _")
  idrawstring(s.b, s.f, 5, ix,iy, "                 _")
  iy-=8

  idrawstring(s.b, s.f, 2, ix,iy, "  (_)")
  idrawstring(s.b, s.f, 1, ix,iy, "          |")
  idrawstring(s.b, s.f, 3, ix,iy, "            (_)")
  idrawstring(s.b, s.f, 5, ix,iy, "                (_)")
  iy-=8

  idrawstring(s.b, s.f, 1, ix,iy, "   _ _   _| |_  __ _");
  iy-=8
  idrawstring(s.b, s.f, 1, ix,iy, "  | | | | | | |/ _` |")
  iy-=8
  idrawstring(s.b, s.f, 1, ix,iy, "  | | |_| | | | (_| |")
  iy-=8
  idrawstring(s.b, s.f, 1, ix,iy, " _/ |\\__'_|_|_|\\__'_|")
  iy-=8
  idrawstring(s.b, s.f, 1, ix,iy, "|__/")
end

function sixelplot(x=[], y=[]; title="", xlab="x", ylab="f(x)", xsize=384, ysize=288,
                   xlim=[0,0], ylim=[0,0],
                   typ='l', pch=1, clr=2, dclr=[-1,-1,-1],
                   showsixels=true, showframe=true, showlogo=false)

  if y==[]
    y = x
    n = length(x)
    x = Array{Int64,1}(1:n)
  end

  if length(x)!=length(y)
    error("input arrays have different lengths")
  end

  if rem(ysize,6)!=0
    error("ysize should be a multiple of 6")
  end

  if xlim[1]==xlim[2]
    xmin = minimum(x)
    xmax = maximum(x)
    autox = true
  else
    xmin = xlim[1]
    xmax = xlim[2]
    autox = false
  end

  if ylim[1]==ylim[2]
    ymin = minimum(y)
    ymax = maximum(y)
    autoy = true
  else
    ymin = ylim[1]
    ymax = ylim[2]
    autoy = false
  end

  if xmin>=xmax
    error("invalid values for xlim")
  end
  if ymin>=ymax
    error("invalid values for ylim")
  end

  s = Screen(xsize,ysize)
  clear(s.b, 63)

  if dclr[1]>=0
    setcolor(s.p, clr, dclr[1], dclr[2], dclr[3])
  end

  if showframe; frame(s); end
  if showlogo; julialogo(s); end

  if title!=""
    drawstring(s,1, 0.5,0.95, title, 'u')
  end

  axes(s, xlab,autox,xmin,xmax, ylab,autoy,ymin,ymax)

  drawplot(s, clr, x, y, typ, pch)

  if showsixels; sixels(s); end

  return s
end

function sixelplot(s::Screen, x=[], y=[]; typ='l', pch=2, clr=3, dclr=[-1,-1,-1], showsixels=true)
  if y==[]
    y = x
    n = length(x)
    x = Array{Int64,1}(1:n)
  end
  if dclr[1]>=0
    setcolor(s.p, clr, dclr[1], dclr[2], dclr[3])
  end
  drawplot(s, clr, x, y, typ, pch)
  if showsixels; sixels(s); end
end

end # module


text:
function validiban(iban::AbstractString)
    country2length = Dict(
        "AL" => 28, "AD" => 24, "AT" => 20, "AZ" => 28, "BE" => 16, "BH" => 22, "BA" => 20, "BR" => 29,
        "BG" => 22, "CR" => 21, "HR" => 21, "CY" => 28, "CZ" => 24, "DK" => 18, "DO" => 28, "EE" => 20,
        "FO" => 18, "FI" => 18, "FR" => 27, "GE" => 22, "DE" => 22, "GI" => 23, "GR" => 27, "GL" => 18,
        "GT" => 28, "HU" => 28, "IS" => 26, "IE" => 22, "IL" => 23, "IT" => 27, "KZ" => 20, "KW" => 30,
        "LV" => 21, "LB" => 28, "LI" => 21, "LT" => 20, "LU" => 20, "MK" => 19, "MT" => 31, "MR" => 27,
        "MU" => 30, "MC" => 27, "MD" => 24, "ME" => 22, "NL" => 18, "NO" => 15, "PK" => 24, "PS" => 29,
        "PL" => 28, "PT" => 25, "RO" => 24, "SM" => 27, "SA" => 24, "RS" => 22, "SK" => 24, "SI" => 19,
        "ES" => 24, "SE" => 24, "CH" => 21, "TN" => 24, "TR" => 26, "AE" => 23, "GB" => 22, "VG" => 24)

    # Ensure upper alphanumeric input.
    iban = replace(iban, r"\s", "")

    rst = ismatch(r"^[\dA-Z]+$", iban)
    # Validate country code against expected length.
    rst = rst && length(iban) == country2length[iban[1:2]]
    # Shift and convert.
    iban = iban[5:end] * iban[1:4]
    digs = parse(BigInt, join(parse(Int, ch, 36) for ch in iban))
    return rst && digs % 97 == 1
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ImageFiltering

# Part 1

const nhood_filter = centered( [ 1 1 1; 1 0 1; 1 1 1 ])

get_neighbors( I ) = imfilter( I , nhood_filter, Fill(0) )

function fill_seats( occupied, seats )
    neighbors = get_neighbors( occupied )
    occupied = ( neighbors .== 0 ) .* seats .| occupied
    neighbors = get_neighbors( occupied )
    occupied = ( neighbors .< 4 ) .* occupied
    return occupied
end

function fill_seats!( occupied, seats )
    neighbors = get_neighbors( occupied )
    occupied .= ( neighbors .== 0 ) .* seats .| occupied
    neighbors = get_neighbors( occupied )
    occupied .= ( neighbors .< 4 ) .* occupied
    return occupied
end

function read_input( filename )
    lines = readlines(filename)
    seats = hcat([Vector{Char}(l) .== 'L' for l in lines]...)
    return copy( transpose(seats) )
end

const seats = read_input("part_1\\input.txt")

# Part 2

const graph_neighbors = findall( nhood_filter .== 1 )

const ci_seats = CartesianIndices( seats )

const one_padded_seats = padarray( seats, Fill(1, ize(seats) ) )

const min_seat = ci_seats[1]
const max_seat = ci_seats[end]

function valid_neighbor( n )
    all( (n.I .>= min_seat.I) .& (n.I .<= max_seat.I) )
end
function filter_neighbors!( neighbors )
    filter!( valid_neighbor , neighbors )
end

function get_neighboring_seats( current, seats = one_padded_seats)
    n = 1
    neighbors = CartesianIndex{2}[]
    unit_neighbors = graph_neighbors
    current_neighbors = (current,) .+ graph_neighbors
    while( !isempty(current_neighbors) )
        neighbor_seats = seats[ current_neighbors ]
        append!( neighbors, current_neighbors[ neighbor_seats] )
        unit_neighbors = unit_neighbors[ .!neighbor_seats ]
        n = n + 1
        current_neighbors = (current,) .+ (unit_neighbors .* n)
    end
    filter_neighbors!( neighbors )
end

const seat_neighbors = [get_neighboring_seats( c, one_padded_seats ) for c in ci[ seats ]]

function get_neighbors_part_2( occupied )
    neighbors = zeros( Int64, size( occupied ) )
    neighbors[ seats ] = [sum( occupied[ s ] ) for s in seat_neighbors]
    return neighbors
end

function fill_seats_part_2!( occupied, seats )
    neighbors = get_neighbors_part_2( occupied )
    occupied .= ( neighbors .== 0 ) .* seats .| occupied
    neighbors = get_neighbors_part_2( occupied )
    occupied .= ( neighbors .< 5 ) .* occupied
    return occupied
end


text:
@testset "FPGroups" begin
    A = Alphabet([:a, :A, :b, :B, :c, :C], [2,1,4,3,6,5])

    @test FreeGroup(A) isa FreeGroup
    @test sprint(show, FreeGroup(A)) == "free group on 3 generators"

    F = FreeGroup([:a, :b, :c], A)
    @test sprint(show, F) == "free group on 3 generators"

    a,b,c = gens(F)
    @test c*b*a isa FPGroupElement

    # quotient of F:
    G = FPGroup(F, [a*b=>b*a, a*c=>c*a, b*c=>c*b])

    @test G isa FPGroup
    @test sprint(show, G) == "⟨ a  b  c | \n\t  a*b => b*a  a*c => c*a  b*c => c*b ⟩"
    @test rand(G) isa FPGroupElement

    f = a*c*b
    @test word(f) isa Word{UInt8}

    aG,bG,cG = gens(G)

    @test aG isa FPGroupElement
    @test_throws AssertionError aG == a
    @test word(aG) == word(a)

    g = aG*cG*bG

    @test_throws AssertionError f == g
    @test word(f) == word(g)
    @test word(g) == [1, 5, 3]
    Groups.normalform!(g)
    @test word(g) == [1, 3, 5]

    let g = aG*cG*bG
        # test that we normalize g before printing
        @test sprint(show, g) == "a*b*c"
    end

    # quotient of G
    H = FPGroup(G, [aG^2=>cG, bG*cG=>aG], maxrules=200)

    h = H(word(g))

    @test h isa FPGroupElement
    @test_throws AssertionError h == g
    @test_throws MethodError h*g

    H′ = FPGroup(G, [aG^2=>cG, bG*cG=>aG], maxrules=200)
    @test_throws AssertionError one(H) == one(H′)

    Groups.normalform!(h)
    @test h == H([5])

    @testset "GroupsCore conformance: H" begin
        test_Group_interface(H)
        test_GroupElement_interface(rand(H, 2)...)
    end

end


text:
function index_file(file::String)
    for word in _tokenize(open(readall, file))
        push!(get!(_INDEX, word, Set{String}()), file)
    end
end


text:
using DiffEqOperators

n = 100
x=0.0:0.005:2π
xprime = x[2:(end-1)]
dx=diff(x)
y = exp.(π*x)
y_im = exp.(π*im*x)
yim_ = y_im[2:(end-1)]
y_ = y[2:(end-1)]

@test_broken for dor in 1:6, aor in 1:6

    D1 = CenteredDifference(dor,aor,dx,length(x))

    #take derivative
    yprime1 = D1*y

    #test result
    @test_broken yprime1 ≈ (π^dor)*y_ # test operator with known derivative of exp(kx)

    #take derivatives
    y_imprime1 = D1*y_im

    #test result
    @test_broken y_imprime1 ≈ ((pi*im)^dor)*yim_ # test operator with known derivative of exp(jkx)

    #TODO: implement specific tests for the left and right boundary regions, waiting until after update
end


text:
using RecipesBase
using ColorTypes
import PlotUtils: cgrad

const mygreen = RGBA{Float64}(151/255,180/255,118/255,1)
const mygreen2 = RGBA{Float64}(113/255,161/255,103/255,1)
const myblue = RGBA{Float64}(74/255,144/255,226/255,1)

@userplot Streamlines

@recipe function f(s::Streamlines)
    elements = s.args[3]

    if (isa(elements, Array) || isa(elements, Tuple)) &&
        ConformalBody in typeof.(elements)

        b = elements[findfirst(isa.(elements,PotentialFlow.ConformalBody))]
        ζ = [r*exp(im*θ) for θ in s.args[2], r in s.args[1]]
        ψ = streamfunction(ζ,elements)
        Z = conftransform(ζ,b)

        # Determine an automatic range of contour levels
        vmin = minimum(ψ)
        vmax = maximum(ψ)
        vmean = vmax-vmin
        vfact = 1
        n = 31
        vmin = vmean+vfact*(vmin-vmean)
        vmax = vmean+vfact*(vmax-vmean)
        v = range(vmin,vmax,length=n)
        if (sign(vmin) != sign(vmax))
            v = v .- v[abs.(v).==minimum(abs.(v))];
        end

        @series begin
          seriestype --> :contour
            grid --> :none
            aspect_ratio --> 1
            linewidth --> 1
            legend --> :none
            seriescolor --> [:black, :black]
            levels --> v
            real.(Z), imag.(Z), ψ
        end

    else
      z = [x + im*y for y in s.args[2], x in s.args[1]]
      ψ = streamfunction(z, elements)

      @series begin
          seriestype --> :contour
          grid --> :none
            seriescolor --> cgrad([:grey, :grey])

            s.args[1], s.args[2], ψ
      end
    end
end

@recipe function plot(points::Array{P};
                      source_marker = :xcross,
                      vortex_marker = :circle) where {P <: Union{Points.Point, Blobs.Blob}}
    z = Elements.position(points)
    x = real.(z)
    y = imag.(z)

    sources = Int[]
    vortices = Int[]

    for (i, p) in enumerate(points)
        if abs2(real(p.S)) > 0
            push!(vortices, i)
        end

        if abs2(imag(p.S)) > 0
            push!(sources, i)
        end
    end

    if !isempty(vortices)
        @series begin
            seriestype --> :scatter
            markershape := vortex_marker
            marker_z --> Elements.circulation.(points[vortices])
            x[vortices], y[vortices]
        end
    end

    if !isempty(sources)
        @series begin
            seriestype --> :scatter
            markershape := source_marker
            marker_z --> Elements.flux.(points[sources])
            x[sources], y[sources]
        end
    end
end

@recipe function plot(s::Vortex.Sheet)
    z = s.zs
    Γ = Elements.circulation.(s.blobs)
    line_z --> 0.5(Γ[1:end-1] + Γ[2:end])./abs.(diff(z))
    x := real.(z)
    y := imag.(z)
    ()
end

@recipe function plot(p::Plate)
    z = [p.zs[1], p.zs[end]]
    linecolor --> :black
    x := real.(z)
    y := imag.(z)
    ()
end

@recipe function plot(b::ConformalBody)
    z = [b.zs; b.zs[1]]
    linecolor --> mygreen
    fillrange --> 0
    fillcolor --> mygreen
    aspect_ratio --> 1
    legend --> :none
    x := real.(z)
    y := imag.(z)
    ()
end

@recipe function plot(c::Corner{N};radius=2.0) where {N}

    θ0 = angle(c) - 0.5*N*π*c.ν
    θ1 = angle(c) + 0.5*N*π*c.ν
    θ = range(θ1-2π,θ0,length=200)
    xv = [0;2*radius*cos.(θ);0]
    yv = [0;2*radius*sin.(θ);0]
    z = xv+im*yv

    linecolor --> mygreen
    fillrange --> 0
    fillcolor --> mygreen
    aspect_ratio --> 1
    legend --> :none
    x := real.(z)
    y := imag.(z)
    ()
end

const VortexSystem = NTuple{N, Union{Element, Tuple, Array{V} where {V <: Element}}} where N
function RecipesBase.RecipesBase.apply_recipe(plotattributes::Dict{Symbol, Any}, sys::VortexSystem)
    series_list = RecipesBase.RecipeData[]
    for s in sys
        append!(series_list, RecipesBase.RecipesBase.apply_recipe(copy(plotattributes), s) )
    end
    series_list
end


text:
module rcmnl

using StatsFuns: logsumexp
using Base.Iterators: product
using FastGaussQuadrature: gausshermite
using StaticArrays
using LinearAlgebra: mul!, checksquare
using StatsFuns: sqrt2, sqrtπ, invsqrt2π, invπ

export simlogL

const AV = AbstractVector
const AM = AbstractMatrix

include("integration.jl")

"log lik of choices yᵢ conditional on mean inside good utilities uᵢ and shock eₖ"
function logLi!(tmp::Vector{T}, y::AV, u::AM, e::AV) where {T}
    length(y) == size(u,2) || throw(DimensionMismatch())
    logL = zero(T)
    tmp[1] = 0
    for j = 1:size(u,2)
        tmp[2:end] .= view(u,:,j) .+ e
        yit = y[j]
        logL += tmp[yit] - logsumexp(tmp)
    end
    return logL
end

"simulated logL of data using `npts` quadrature points"
function simlogL(y, X, θ::AbstractVector{T}; npts=15) where {T}

    # unpack parameters
    β = reshape(view(θ, 1:8), :, 2)
    Σchol = SMatrix{2,2}(θ[end-2], θ[end-1], 0, θ[end])

    # zip together vectors of nodes + weights for Multidim quad
    # and take cartesian product
    # https://gist.github.com/markvdw/f9ca12c99484cf2a881e84cb515b86c8
    vw = zip(gausshermite(npts)...)
    quadprod = product(vw, vw)

    # preallocate temp vars
    u = zeros(T, 2, size(y,1))
    tmp = zeros(T, 3)
    lse = zeros(T, size(quadprod))

    # loop over individuals
    # would be good to parallelize!
    sll = zero(T)
    for j in 1:size(y,2)

        # mean utilities of inside good
        @views mul!(u, β', X[:,j,:]')

        # loop over cartesian product of quadrature wts/nodes
        # (don't have to expand these)
        for (i, ((v1,w1), (v2,w2))) in enumerate(quadprod)

            # transform iid vₖ normal to correlatd eₖ
            # using StaticArrays b/c this allocates less
            # probably unecessary
            v = SVector(v1,v2)
            e = sqrt2 .* Σchol * v

            # quadrature weight wₖ * π*N/2
            w = w1*w2*π

            # log( wₖ * likᵢ conditional on eₖ)
            lse[i] = log(w) + logLi!(tmp, view(y,:,j), u, e)
        end

        # compute ∑ wtₖ * likₖ
        sll += logsumexp(lse)
    end
    return - invπ * sll
end

end # module


text:
# import .ImmerseCompose: Iterables, iterable, native
# import .DisplayGadfly: Figure, render_backend
import Compose: SVGClass, Form, Circle, Line, Backend
# using Gtk, Cairo, GtkUtilities

path2mask(backend::Backend, pathx, pathy) = path2mask(round(Int,width(Compose.surface(backend))), round(Int,height(Compose.surface(backend))), pathx, pathy)

function path2mask(w, h, pathx, pathy)
    length(pathx) == length(pathy) || error("pathx and pathy must have the same length")
    if isempty(pathx)
        return falses(w, h)
    end
    # Use Cairo to make the mask
    data = Array{UInt32}(undef, w, h)
    surf = Cairo.CairoImageSurface(data, Cairo.FORMAT_RGB24, flipxy=false)
    ctx = CairoContext(surf)
    set_source_rgb(ctx, 0, 0, 0)
    paint(ctx)
    move_to(ctx, pathx[end], pathy[end])
    for i = 1:length(pathx)
        line_to(ctx, pathx[i], pathy[i])
    end
    set_source_rgb(ctx, 1, 1, 1)
    fill(ctx)
    mask = reshape([x & 0xff > 0 for x in data], size(data))
    mask
end

const pathx = Float64[]
const pathy = Float64[]

function lasso_select_cb(f::Figure)
    c = f.canvas
    push!((c.mouse,:button1press), Gtk.@guarded (widget,event) -> begin
        empty!(pathx)
        empty!(pathy)
        push!(pathx, event.x)
        push!(pathy, event.y)
        push!((c.mouse,:button1motion), dragging)
    end)
    push!((c.mouse,:button1release), Gtk.@guarded (widget,event) -> begin
        # Restore the original mouse handlers
        pop!((c.mouse, :button1press))
        pop!((c.mouse, :button1motion))
        pop!((c.mouse, :button1release))
        # Toggle the lasso selection button
        lasso_button = guidata[c, :lasso_button]
        setproperty!(lasso_button, :active, false)
        # Redraw the canvas
        Gtk.draw(c)
        select_points(f, pathx, pathy)
    end)
end

Gtk.@guarded function dragging(widget,event)
    ctx = getgc(widget)
    move_to(ctx, pathx[end], pathy[end])
    line_to(ctx, event.x, event.y)
    set_source_rgb(ctx, 0, 0, 0)
    stroke(ctx)
    push!(pathx, event.x)
    push!(pathy, event.y)
    Gtk.reveal(widget, false)
end

function select_points(f::Figure, pathx, pathy)
    c = f.canvas
    # Check to see if the figure is empty (issue #56)
    isempty(f) && return nothing
    # Find the forms
    forms = find_panelforms(f.cc)
    # Create the mask
    backend = render_backend(c)
    mask = path2mask(backend, pathx, pathy)
    # Determine which points are in the mask
    coords = guidata[c,:panelcoords][1]
    selections = Dict()
    for form in forms
        selections[form] = find_inmask(backend, coords, form, mask)
    end
    # Run the callback
    lasso_button = guidata[c, :lasso_button]
    cb = guidata[lasso_button, :callback]
    cb(f.figno, selections)
end

const lasso_default = (figno, selections) -> export_selection(selections)

"""

`lasso_initialize(figno, cb)` sets the callback `cb` to run when the
user has selected points with the lasso.  When you have only one
"object" (Geom) in your plot, your callback can be of the form
```
    function my_simple_lasso_callback(figno, selections)
        println("You selected points ", first(values(selections)))
    end
```
For a `Geom.point` object, this would print a vector of integers
corresponding to the selected points.

More generally, the callback should have the syntax
```
    function my_lasso_callback(figno, selections)
        for (form, indexes) in selections
            if form.tag == :mydots
                println("In figure \$figno, from :mydots you selected ", indexes)
            end
        end
    end
```
`selections` is a `Dict` of `form=>indexes` pairs.  `form` is a
`Compose.Form`, the raw objects rendered by Gadfly; you may especially
want to query its `tag` to determine its identity (assuming you've
assigned a tag).  `indexes` is a vector describing the items selected;
for `Circle` forms, each element will be an `Int`, whereas for `Line`
Forms (which can hold multiple lines, perhaps drawn in different colors)
each element will be an `Tuple{Int,Int}` describing the line number
and vertex number.

To assign a tag to a Gadfly object, add the `tag` keyword argument to
the geometry, e.g. `Geom.point(tag = :mydots)`.
"""
function lasso_initialize(f::Figure, cb=lasso_default)
    c = f.canvas
    lasso_button = guidata[c, :lasso_button]
    guidata[lasso_button, :callback] = cb
    # Work around Gtk #161 & #185
    # Gtk.signal_connect(lasso_button, "clicked") do widget
    #     lasso_select_cb(f)
    # end
    #Gtk.signal_connect(lasso_wrapper, lasso_button, "clicked", Void, (), false, f)
end
lasso_initialize(i::Int, cb=lasso_default) = lasso_initialize(Figure(i), cb)

Gtk.@guarded function lasso_wrapper(widgetptr::Ptr, f)
    widget = convert(Gtk.GtkToggleToolButtonLeaf, widgetptr)
    if Gtk.get_gtk_property(widget, :active, Bool)
        lasso_select_cb(f)
    end
    nothing
end

function export_selection(selections)
    # Extract (tag,index) pairs. We don't use a Dict in case there are
    # multiple untagged objects.
    indexes = Array{Tuple{Symbol,Any}}(undef,length(selections))
    i = 0
    nonempty = false
    for (form,indx) in selections
        indexes[i+=1] = (form.tag,indx)
        nonempty |= !isempty(indx)
    end
    nonempty || return nothing
    resp, varname = Gtk.input_dialog("Pick variable name in Main for exporting selection", "selection", (("Cancel",0), ("Store",1)))
    if resp == 0
        return nothing
    end
    if !isempty(varname)
        sym = Symbol(varname)
        eval(Main, :($sym = $indexes))
    end
    nothing
end


text:
include("test_triangleints.jl")
include("test_contour0.jl")
include("test_contour1.jl")
include("test_contour2.jl")
include("scan1.jl")


text:
export GaussianNB

# X: Feature vector as a p x n matrix
# c: Vector of classes
# k: Number of distinct classes
import Distributions
import ScikitLearnBase

type GaussianNB
    # The model hyperparameters (not learned from data)
    bias::Float64

    # The parameters learned from data
    mu::Matrix{Float64} # Matrix of cluster centers: p x k
    sigma::Matrix{Float64} # Diagonals of covariance matrices for each class
	p::Vector{Float64} # Vector of probabilities for each class
	counts::Matrix{Int}

    # A constructor that accepts the hyperparameters as keyword arguments
    # with sensible defaults
    GaussianNB(; bias=0.0f0) = new(bias)
end

# This will define `clone`, `set_params!` and `get_params` for the model
ScikitLearnBase.@declare_hyperparameters(GaussianNB, [:bias])

# GaussianNB is a classifier
ScikitLearnBase.is_classifier(::GaussianNB) = true   # not required for transformers

function Distributions.rand(d::GaussianNB)
	p = size(d.mu, 1)
	c = Distributions.draw(d.drawtable)
	x = Array(Float64, p)
	for dim in 1:p
		x[dim] = rand(Normal(d.mu[dim, c], d.sigma[dim, c]))
	end
	x, c
end

function Distributions.rand(d::GaussianNB, n::Integer)
	p = size(d.mu, 1)
	X = Array(Float64, p, n)
	c = Array(Int, n)
	for obs in 1:n
		c[obs] = Distributions.draw(d.drawtable)
		for dim in 1:p
			X[dim, obs] = rand(Normal(d.mu[dim, c[obs]],
				                      d.sigma[dim, c[obs]]))
		end
	end
	X, c
end

function Distributions.mean(d::GaussianNB)
	p, c = size(d.mu)
	mx = zeros(Float64, p)
	for cl in 1:c
		mx += d.p[cl] * d.mu[:, c]
	end
	mx /= c
	mc = 0.0
	for cl in 1:c
		mc += d.p[cl] * cl
	end
	mx, mc
end

function ScikitLearnBase.fit!(model::GaussianNB, X::Matrix, c::Vector)
	p, n = size(X)
	nclasses = maximum(c)

	for i in 1:n
		model.mu[:, c[i]] += X[:, i]
		model.counts[c[i]] += 1
	end
	for cl in 1:nclasses
		model.mu[:, cl] /= model.counts[cl]
	end
	for i in 1:n
		model.sigma[:, c[i]] += (X[:, i] - model.mu[:, c[i]]).^2
	end
	for cl in 1:nclasses
		model.sigma[:, cl] = sqrt(model.sigma[:, cl] / (model.counts[cl] - 1)) + 1e-8
	end
	return model
end

function Distributions.logpdf(d::GaussianNB, x::Vector, c::Real)
	p = length(x)
	res = log(d.p[c])
	for dim in 1:p
		res += logpdf(Normal(d.mu[dim, c], d.sigma[dim, c]), x[dim])
	end
	return res
end

function Distributions.logpdf(d::GaussianNB, X::Matrix, c::Vector)
	p, n = size(X)
	res = zeros(Float64, n)
	for obs in 1:n
		res[obs] = log(d.p[c[obs]])
		for dim in 1:p
			res[obs] += logpdf(Normal(d.mu[dim, c[obs]],
				                      d.sigma[dim, c[obs]]),
			                   X[dim, obs])
		end
	end
	return res
end

function Distributions.loglikelihood(d::GaussianNB, X::Matrix, c::Vector)
	p, n = size(X)
	res = 0.0
	for obs in 1:n
		res += logpdf(d, X[:, obs], c[obs])
	end
	return res
end

function ScikitLearnBase.predict(d::GaussianNB, X::Matrix)
	nclasses = length(d.p)
	p, n = size(X)
	res = Array(Int, n)
	ll = Array(Float64, nclasses)
	for obs in 1:n
		for cl in 1:nclasses
			ll[cl] = logpdf(d, X[:, obs], cl)
		end
		res[obs] = indmax(ll)
	end
	return res
end


text:
using ArgParse
using CalibrateEDMF.TurbulenceConvectionUtils
using CalibrateEDMF.HelperFuncs
include("DiagnosticsTools.jl")

s = ArgParseSettings()
@add_arg_table s begin
    "--results_dir"
    help = "Path to the results directory (CEDMF output)"
    arg_type = String
    "--tc_output_dir"
    help = "Path to store TC output"
    arg_type = String
    "--run_set"
    help = "Set of cases to run TC. {reference, validation, test}."
    arg_type = String
    default = "validation"
    "--run_set_config"
    help = "Path to config if test is used for --run_set."
    arg_type = String
    default = nothing
    "--method"
    help = "Method for computing optimal parameters"
    arg_type = String
    default = "best_nn_particle_mean"
    "--metric"
    help = "mse metric to find the minimum of"
    arg_type = String
    default = "mse_full"
end
parsed_args = parse_args(ARGS, s)
include(joinpath(parsed_args["results_dir"], "config.jl"))
config = get_config()

"""
    run_TC_optimal(results_dir::String, tc_output_dir::String, config::Dict; method::String = "best_nn_particle_mean", metric::String = "mse_full" )
Given path to the results directory of completed calibration run, run TC on a given set of cases (`run_cases`).
    using optimal parameters and save resulting TC stats files.
Inputs:
 - results_dir - directory containing CEDMF `Diagnostics.nc` file.
 - tc_output_dir - directory to store TC output
 - namelist_args - additional arguments passed to the TurbulenceConvection namelist.
 - run_cases - Dictionary of cases to run. Must contain {"case_name","scm_suffix", "y_dir"}.
 - `method`: method for computing optimal parameters. Use parameters of:
        "best_particle" - particle with lowest mse in training (`metric` = "mse_full") or validation (`metric` = "val_mse_full") set.
        "best_nn_particle_mean" - particle nearest to ensemble mean for the iteration with lowest mse.
 - `metric` : mse metric to find the minimum of {"mse_full", "val_mse_full"}.
"""

function run_TC_optimal(
    results_dir::String,
    tc_output_dir::String,
    namelist_args::Union{Vector, Nothing} = nothing,
    run_cases::Union{Dict, Nothing} = nothing;
    method::String = "best_nn_particle_mean",
    metric::String = "mse_full",
)

    u_names, u = optimal_parameters(joinpath(results_dir, "Diagnostics.nc"); method = method, metric = metric)

    for i in 1:length(run_cases["case_name"])
        @info "Running TC suffix: $(run_cases["scm_suffix"][i])"
        les_path = if run_cases["case_name"][i] == "LES_driven_SCM"
            get_stats_path(run_cases["y_dir"][i])
        else
            nothing
        end
        run_SCM_handler(
            run_cases["case_name"][i],
            tc_output_dir;
            u = u,
            u_names = u_names,
            namelist_args = namelist_args,
            uuid = run_cases["scm_suffix"][i],
            les = les_path,
        )
    end
end

if parsed_args["run_set"] in ("reference", "validation")
    namelist_args = get_entry(config["scm"], "namelist_args", nothing)
    run_cases = get(config, parsed_args["run_set"], nothing)
    @info "Running TC on $(parsed_args["run_set"]) set"
    run_TC_optimal(
        parsed_args["results_dir"],
        parsed_args["tc_output_dir"],
        namelist_args,
        run_cases;
        method = parsed_args["method"],
        metric = parsed_args["metric"],
    )

elseif parsed_args["run_set"] == "test"

    if isnothing(parsed_args["run_set_config"])
        error("--run_set_config must be specified if --run_set is test")
    else
        include(parsed_args["run_set_config"])
        run_cases = get_reference_config(ScmTest())
    end
    namelist_args = get_entry(config["scm"], "namelist_args", nothing)
    @info "Running TC on $(parsed_args["run_set"]) set"
    run_TC_optimal(
        parsed_args["results_dir"],
        parsed_args["tc_output_dir"],
        namelist_args,
        run_cases;
        method = parsed_args["method"],
        metric = parsed_args["metric"],
    )

else
    error("--run_set must be one of {reference, validation, test}")
end