text:
module LevelSetAdvection

export  get_cellcenter_coordinates, get_cell_diameter,
        upwind2d_step, ReinitHJ2d_update, Signum

##############################################
using Reexport
@reexport using Gridap
@reexport using Gridap.Geometry # get_grid_topology
@reexport using Gridap.ReferenceFEs # get_faces
@reexport using Gridap.TensorValues
##############################################

# Auxiliar functions
Signum(x;β=2.0) = tanh.(β.*x)
SigNum(x;ϵ=0.1) = x/√(x^2 + ϵ^2)

Sₑ(ϕ;ϵ=0.1) = 0.5/ϵ.*(1 .+ cos.(π*ϕ/ϵ)).*(abs.(ϕ).<=ϵ)
Hₑ(ϕ;ϵ=0.1) = 1.0.*(ϕ.>ϵ) .+ 0.5*(1 .+ ϕ/ϵ .+ 1/ϵ*sin.(π*ϕ/ϵ)).*(abs.(ϕ).<=ϵ)

"""
Godunov Hamiltonian 2D:
Hg(sgnϕ,a,b,c,d) ≈ |∇ϕ|

"""
function Hg(sgnϕ,a,b,c,d)
    h = 0

    aM = a*(a>=0)
    am = a - aM
    bM = b*(b>=0)
    bm = b - bM
    cM = c*(c>=0)
    cm = c - cM
    dM = d*(d>=0)
    dm = d - dM

    if sgnϕ>0
        h = √(max(am^2,bM^2) + max(cm^2,dM^2) )
    else
        h = √(max(aM^2,bm^2) + max(cM^2,dm^2) )
    end

    return h
end

"""
Get cell centered coordinates
x,y = get_center_coord(Ω::Triangulation)
"""
function get_cellcenter_coordinates(Ω::Triangulation)
    #=
    # This is slower
    xy = get_cell_points(Ω) # all vertices in cells
    xc = [sum(xy.cell_phys_point[i])/length(xy.cell_phys_point[i]) for i in 1:length(xy.cell_phys_point)] # center at cell
    =#
    # this is faster
    xy = get_cell_coordinates(Ω) # same as get_cell_points(topo)
    midpoint(xs) = sum(xs)/length(xs)
    xc = lazy_map(midpoint,xy)
    return xc
end

"""
Get the cell diameters
(perhaps there is a better way)

diam::Vector = get_cell_diameter(Ω::Triangulation)
"""
function get_cell_diameter(Ω::Triangulation)
    xy = get_cell_coordinates(Ω)
    #diam(x) = maximum(norm.(x[][:] - x[][end:-1:1])) #! not all combinations! [1.2.3.4] - [4,3,2,1] ... but seems to be enough?!
    #return lazy_map(diam,xy) # not working
    #? not optimal
    d = zeros(num_cells(Ω))
    for (k,x) in enumerate(xy)
        d_ij = 0.0
        for i in eachindex(x)
            for j in i:length(x)
                d_ij = max(d_ij,norm(x[i] - x[j]))
            end
        end
        d[k] = d_ij
    end

    return d
end

"""
Local Finite Difference derivatives across faces

returns ϕx⁺, ϕx⁻, ϕy⁺, ϕy⁻
"""
function local_FD_derivatives(ϕ,S,N,W,E,xc)
    #? Trick: We repeat to BoundsError when 'e' is on the boundary
    S = repeat(S,2)
    N = repeat(N,2)
    E = repeat(E,2)
    W = repeat(W,2)

    # Finite differences: According to the ordering in Gridap #!CHECK for unstructured!
    # * Here we avoid division by 0 on the boundary
    ϕyM = (ϕ[N[2]] - ϕ[N[1]])/(xc[N[2]][2]-xc[N[1]][2] + 1e-30)  # Dy⁺(ϕ)  point upwards
    ϕym = (ϕ[S[2]] - ϕ[S[1]])/(xc[S[2]][2]-xc[S[1]][2] + 1e-30)  # Dy⁻(ϕ)  point downwards
    ϕxM = (ϕ[E[2]] - ϕ[E[1]])/(xc[E[2]][1]-xc[E[1]][1] + 1e-30)  # Dx⁺(ϕ)  point to the right
    ϕxm = (ϕ[W[2]] - ϕ[W[1]])/(xc[W[2]][1]-xc[W[1]][1] + 1e-30)  # Dx⁻(ϕ)  point to the left

    return ϕxM,ϕxm,ϕyM,ϕym
end

"""
simple Upwind 2D step

We want to solve

∂ϕ/∂t + V|∇ϕ| = 0

       ϕ(t=0) = ϕ₀

Via Godunov FD scheme:
`Gϕ = upwind2d_step(G::GridTopology,xc::AbstractArray,ϕ::AbstractArray,V::AbstractArray)`

where
Gϕ ≈ V|∇ϕ|

Then, update your function as
ϕ(n+1) = ϕ(n) - Δt⋅Gϕ(n)
"""
function upwind2d_step(topo::GridTopology,xc,ϕ::AbstractArray,V::AbstractArray)

    # Recover mesh connectivity
    nelem = num_cells(topo) #length(topo.cell_type)

    if nelem != length(ϕ)
        error("Vector ϕ has different length (",length(ϕ),") to number of cells in the mesh (",nelem,")")
    end

    dim = num_cell_dims(topo)
    cell_to_faces  = get_faces(topo,dim,dim-1)
    face_to_cells  = get_faces(topo,dim-1,dim)

    if dim !=2
        error("Sorry: 3D is not implemented yet in `ls_advection.jl`")
    end

    g = Array{Float64,1}(undef,length(ϕ))

    for e in 1:nelem

        # Recovering neighbourhood around element 'e'
        S,N,W,E = face_to_cells[cell_to_faces[e]] # should have 4 terms in 2D-structured grid
        ϕxM,ϕxm,ϕyM,ϕym = local_FD_derivatives(ϕ,S,N,W,E,xc)

        # Godunov:
        g⁺ = √(max(min(ϕxM,0)^2,max(ϕxm,0)^2) + max(min(ϕyM,0)^2,max(ϕym,0)^2))
        g⁻ = √(max(max(ϕxM,0)^2,min(ϕxm,0)^2) + max(max(ϕyM,0)^2,min(ϕym,0)^2))

        V⁺ = max(V[e],0)  # upwind with V magnitude
        V⁻ = min(V[e],0)

        g[e] = V⁺ * g⁺ + V⁻ * g⁻
    end

    return g
end

"""
solves some time steps for

    ∂ϕ/∂t + sign(ϕ₀)(|∇ϕ|-1) = 0     in the Grid

                      ϕ(t=0) = ϕ₀
"""
function ReinitHJ2d_update(topo::GridTopology,xc,ϕ::AbstractArray,niter::Int;vmax::AbstractFloat=1.0,scheme="RK2",show_errormsg::Bool=false,Δt::Real=0.1)
    # Euler & RK2 copyed from
    # "*On reinitilizing level set functions*", C. Min, J Comput.Phys 299:8 (2010) DOI:10.1016/j.jcp.2009.12.032

    #TODO: Check if we are on QUADs (and structured!?)

    # Recover mesh connectivity
    nelem = length(topo.cell_type)

    if nelem != length(ϕ)
        error("Vector ϕ has different length (",length(ϕ),") to number of cells in the mesh (",nelem,")")
    end

    dim = num_cell_dims(topo)
    cell_to_faces  = get_faces(topo,dim,dim-1)
    face_to_cells  = get_faces(topo,dim-1,dim)

    if dim !=2
        error("Sorry: 3D is not implemented yet in `ls_advection.jl`")
    end

    #TODO: Choose the signum function: TEST required!
    #signψ₀ = sign.(ϕ)
    signψ₀ = Signum(ϕ)
    #signψ₀ = SigNum.(ϕ,ϵ=dx)
    #signψ₀ = lazy_map(SigNum(x,ϵ=1.2*dx)->x,ϕ) # not working

    CFL = 0.3 # later: Let the user choose it

    if lowercase(scheme) == "upwind"
        #? Upwind scheme
        # slow convergence and show oscillations near the boundaries, even with smooth signum function and Laplacian regularization

        # init
        δ = Array{Float64,1}(undef,length(ϕ)) ###

        res0 = 1e+10

        for iter in 1:niter
            for e in 1:nelem

                # Recovering neighbourhood around element 'e'
                S,N,W,E = face_to_cells[cell_to_faces[e]] # should be 4 in 2D!
                ϕxM,ϕxm,ϕyM,ϕym = local_FD_derivatives(ϕ,S,N,W,E,xc)

                if signψ₀[e]<0
                    δ[e] = √(max(min(ϕxM,0)^2,max(ϕxm,0)^2) + max(min(ϕyM,0)^2,max(ϕym,0)^2)) - 1.0 # g⁻ + 1
                else
                    δ[e] = √(max(max(ϕxM,0)^2,min(ϕxm,0)^2) + max(max(ϕyM,0)^2,min(ϕym,0)^2)) - 1.0 # g⁺ - 1
                end
                δ[e] = δ[e]*signψ₀[e] #- ν*Δϕ  # when applies the smoothed signum

            end

            # This is the scheme with Laplacian regularization:
            # ϕ(n+1) = ϕ(n) - Δt⋅sign(ϕ₀)⋅(|∇ϕ(n)| - 1) + Δt⋅ν⋅Δϕ
            # ϕ(n+1) ≈ ϕ(n) - Δt[ -(ϕ₀<0)⋅g⁻ + (ϕ₀>0)⋅g⁺ - sign(ϕ₀) ] + Δt⋅ν⋅Δϕ
            # ϕ(n+1) ≈ ϕ(n) - Δt[ -(ϕ₀<0)⋅(g⁻-1) + (ϕ₀>0)⋅(g⁺-1) ] + Δt⋅ν⋅Δϕ
            # ϕ(n+1) ≈ ϕ(n) - Δt⋅[sign(ϕ₀)⋅δ(n) - ν⋅Δϕ(n)]

            res1 = norm(δ)

            if res0<res1
                Δt *= 0.9
                if show_errormsg
                    print("R.")
                end
            end

            ϕ = ϕ - Δt .* δ

            if show_errormsg
                println("residue[$(iter)]: ",res1*Δt)
            end
            res0=res1
        end

    elseif lowercase(scheme) == "euler"
        #? Forward Euler method using one-sided ENO Finite Differences
        # better than upwind

        ν = 0.001  # smoothing ... just in case

        ϕ1 = copy(ϕ)
        for iter in 1:niter
            ϕold = ϕ

            for e in 1:nelem

                # Recovering neighbourhood around element 'e'
                S,N,W,E = face_to_cells[cell_to_faces[e]] # should be 4 in 2D!
                ϕxM,ϕxm,ϕyM,ϕym = local_FD_derivatives(ϕ,S,N,W,E,xc)

                ϕ[e] = ϕ[e] - Δt*signψ₀[e]*(Hg(signψ₀[e],ϕxM,ϕxm,ϕyM,ϕym) - 1.0) #+ ν*Δϕ
            end

            if show_errormsg
                error_ = norm(ϕ - ϕold)/norm(ϕold)
                println("error[$(iter)] = ",error_)
            end
        end

    elseif lowercase(scheme) == "rk2"
        #? TVD Runge–Kutta method
        ϕ1 = ϕ
        ϕ2 = ϕ

        for iter in 1:niter
            ϕold = ϕ

            # step 1 in RK2
            for e in 1:nelem

                # Recovering neighbourhood around element 'e'
                S,N,W,E = face_to_cells[cell_to_faces[e]] # should be 4 in 2D!
                ϕxM,ϕxm,ϕyM,ϕym = local_FD_derivatives(ϕ,S,N,W,E,xc)

                ϕ1[e] = ϕ[e] - Δt*signψ₀[e]*(Hg(signψ₀[e],ϕxM,ϕxm,ϕyM,ϕym) - 1.0)
            end
            # step 2 in RK2
            for e in 1:nelem

                # Recovering neighbourhood around element 'e'
                S,N,W,E = face_to_cells[cell_to_faces[e]] # should be 4 in 2D!
                ϕxM,ϕxm,ϕyM,ϕym = local_FD_derivatives(ϕ,S,N,W,E,xc)

                ϕ2[e] = ϕ1[e] - Δt*signψ₀[e]*(Hg(signψ₀[e],ϕxM,ϕxm,ϕyM,ϕym) - 1.0)
            end
            # step 3 in RK2
            ϕ = (ϕ1 + ϕ2)/2.0

            if show_errormsg
                error_ = norm(ϕ - ϕold)/norm(ϕold)
                println("error[$(iter)] = ",error_)
            end
        end

    end # if scheme ...

    return ϕ
end

end

max_stars_repo_name: joaquinmura/lsAdvect

max_issues_repo_name: joaquinmura/lsAdvect

max_forks_repo_name: joaquinmura/lsAdvect


text:
using TableReader
using Dates
using Test

@testset "readtsv" begin
    @testset "simple" begin
        # no data
        buffer = IOBuffer("""
        col1
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1]
        @test size(df) == (0, 1)

        buffer = IOBuffer("""
        col1\tcol2\tcol3
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test size(df) == (0, 3)

        # integers
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        1\t23\t456
        -10\t-99\t0
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1, -10]
        @test df[:col2] == [23, -99]
        @test df[:col3] == [456, 0]

        # floats
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        1.0\t1.1\t12.34
        -1.2\t0.0\t-9.
        .000\t.123\t100.000
        1e3\t1.E+123\t-8.2e-00
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1.0, -1.2, 0.000, 1e3]
        @test df[:col2] == [1.1, 0.0, 0.123, 1e123]
        @test df[:col3] == [12.34, -9.0, 100.000, -8.2]

        # bools
        buffer = IOBuffer("""
        col1\tcol2
        true\tfalse
        True\tFalse
        TrUe\tFaLsE
        TRUE\tFALSE
        t\tf
        T\tF
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2]
        @test df[:col1] == [true, true, true, true, true, true]
        @test df[:col2] == [false, false, false, false, false, false]

        # strings
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        a\tb\tc
        foo\tbar\tbaz
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == ["a", "foo"]
        @test df[:col2] == ["b", "bar"]
        @test df[:col3] == ["c", "baz"]

        # dates
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        2015-12-21\t2019-01-01\t1999-09-11
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [Date(2015, 12, 21)]
        @test df[:col2] == [Date(2019, 1, 1)]
        @test df[:col3] == [Date(1999, 9, 11)]

        # datetimes
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        2015-12-21T00:00:00\t2015-12-21T11:22:33\t2015-12-21T11:22:33.444
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [DateTime(2015, 12, 21, 0, 0, 0)]
        @test df[:col2] == [DateTime(2015, 12, 21, 11, 22, 33)]
        @test df[:col3] == [DateTime(2015, 12, 21, 11, 22, 33, 444)]

        # datetimes (delimited by space)
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        2015-12-21 00:00:00\t2015-12-21 11:22:33\t2015-12-21 11:22:33.444
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [DateTime(2015, 12, 21, 0, 0, 0)]
        @test df[:col2] == [DateTime(2015, 12, 21, 11, 22, 33)]
        @test df[:col3] == [DateTime(2015, 12, 21, 11, 22, 33, 444)]
    end

    @testset "tricky float" begin
        # NaN
        buffer = IOBuffer("""
        col1\tcol2\tcol3\tcol4
        nan\tNan\tNaN\tNAN
        -nan\t+Nan\t-NaN\t+NAN
        """)
        df = readtsv(buffer)
        @test all(isnan.(df[:col1]))
        @test all(isnan.(df[:col2]))
        @test all(isnan.(df[:col3]))
        @test all(isnan.(df[:col4]))

        # Inf
        buffer = IOBuffer("""
        col1\tcol2\tcol3\tcol4
        inf\tInf\tInF\tINF
        -inf\t+Inf\t-InF\t+INF
        infinity\tInfinity\tInFiNiTy\tINFINITY
        -infinity\t+Infinity\t-InFiNiTy\t+INFINITY
        """)
        df = readtsv(buffer)
        @test all(isinf.(df[:col1]))
        @test all(isinf.(df[:col2]))
        @test all(isinf.(df[:col3]))
        @test all(isinf.(df[:col4]))
    end

    @testset "quotation" begin
        data = """
        "col1"\t"col2"\t"col3"
        "1"\t"23"\t"456"
        "-10"\t"-99"\t"0"
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, -10]
        @test df[:col2] == [23, -99]
        @test df[:col3] == [456, 0]

        # floats
        buffer = IOBuffer("""
        "col1"\t"col2"\t"col3"
        "1.0"\t"1.1"\t"12.34"
        "-1.2"\t"0.0"\t"-9."
        ".000"\t".123"\t"100.000"
        "1e3"\t"1.E+123"\t"-8.2e-00"
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1.0, -1.2, 0.000, 1e3]
        @test df[:col2] == [1.1, 0.0, 0.123, 1e123]
        @test df[:col3] == [12.34, -9.0, 100.000, -8.2]

        # strings
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        "a"\t"b"\t"c"
        "foo"\t"bar"\t"baz"
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == ["a", "foo"]
        @test df[:col2] == ["b", "bar"]
        @test df[:col3] == ["c", "baz"]

        buffer = IOBuffer("""
        col1\tcol2\tcol3
        foo\tbar\tbaz
        foo\tbar\tbaz
        """)
        df = readtsv(buffer)
        @test df[:col1] == ["foo", "foo"]
        @test df[:col2] == ["bar", "bar"]
        @test df[:col3] == ["baz", "baz"]

        # quotation marks in a quoted string
        buffer = IOBuffer("""
        col1\tcol2
        " ""OK"" "\t"\""OK"\""
        """)
        df = readtsv(buffer)
        @test df[:col1] == [" \"OK\" "]
        @test df[:col2] == ["\"OK\""]
    end

    @testset "trimming" begin
        # trimming space
        data = """
        col1\tcol2\tcol3
        1   \t   2\t   3
           4\t   5\t   6
          7 \t 8  \t 9
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, 4, 7]
        @test df[:col2] == [2, 5, 8]
        @test df[:col3] == [3, 6, 9]

        df = readtsv(IOBuffer(data); trim = true)
        @test df[:col1] == [1, 4, 7]
        @test df[:col2] == [2, 5, 8]
        @test df[:col3] == [3, 6, 9]

        df = readtsv(IOBuffer(data); trim = false)
        @test df[:col1] == ["1   ", "   4", "  7 "]
        @test df[:col2] == ["   2", "   5", " 8  "]
        @test df[:col3] == ["   3", "   6", " 9  "]

        data = """
         col1  \t col2 \t col3
         foo   \t  b  \t baz
        """
        df = readtsv(IOBuffer(data); trim = true)
        @test df[:col1] == ["foo"]
        @test df[:col2] == ["b"]
        @test df[:col3] == ["baz"]

        df = readtsv(IOBuffer(data); trim = false)
        @test df[Symbol(" col1  ")] == [" foo   "]
        @test df[Symbol(" col2 ")] == ["  b  "]
        @test df[Symbol(" col3  ")] == [" baz"]
    end

    @testset "missing" begin
        data = """
        col1\tcol2\tcol3
        1\t2\t3
        \t5\t
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] isa Vector{Union{Int,Missing}}
        @test df[:col2] isa Vector{Int}
        @test df[:col3] isa Vector{Union{Int,Missing}}
        @test df[1,:col1] == 1
        @test ismissing(df[2,:col1])
        @test df[1,:col2] == 2
        @test df[2,:col2] == 5
        @test df[1,:col3] == 3
        @test ismissing(df[2,:col3])

        data = """
        col1\tcol2\tcol3
        1.0\t2.2\t-9.8
        \t10\t
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] isa Vector{Union{Float64,Missing}}
        @test df[:col2] isa Vector{Float64}
        @test df[:col3] isa Vector{Union{Float64,Missing}}
        @test df[1,:col1] == 1.0
        @test ismissing(df[2,:col1])
        @test df[1,:col2] == 2.2
        @test df[2,:col2] == 10.0
        @test df[1,:col3] == -9.8
        @test ismissing(df[2,:col3])

        data = """
        col1\tcol2\tcol3
        foo\t\tbar
        baz\tqux\t
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] isa Vector{String}
        @test df[:col2] isa Vector{Union{String,Missing}}
        @test df[:col3] isa Vector{Union{String,Missing}}
        @test df[1,:col1] == "foo"
        @test df[2,:col1] == "baz"
        @test ismissing(df[1,:col2])
        @test df[2,:col2] == "qux"
        @test df[1,:col3] == "bar"
        @test ismissing(df[2,:col3])

        # NA
        data = """
        col1\tcol2
        1\t2
        NA\tNA
        "NA"\t"NA"
        """
        df = readtsv(IOBuffer(data))
        @test df[1,:col1] == 1
        @test ismissing(df[2,:col1])
        @test ismissing(df[3,:col1])
        @test df[1,:col2] == 2
        @test ismissing(df[2,:col2])
        @test ismissing(df[3,:col2])
    end

    @testset "NA" begin
        # NA is case-sensitive.
        data = """
        col1\tcol2
        na\tNa
        NA\tNA
        "Na"\t"nA"
        "NA"\t"NA"
        """
        df = readtsv(IOBuffer(data))
        @test df[1,:col1] == "na"
        @test ismissing(df[2,:col1])
        @test df[3,:col1] == "Na"
        @test ismissing(df[4,:col1])
        @test df[1,:col2] == "Na"
        @test ismissing(df[2,:col2])
        @test df[3,:col2] == "nA"
        @test ismissing(df[4,:col2])
    end

    @testset "large integer" begin
        data = """
        col1\tcol2
        $(typemax(Int))\t$(typemin(Int))
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [typemax(Int)]
        @test df[:col2] == [typemin(Int)]

        # not supported
        data = """
        col1
        99999999999999999999999
        """
        @test_throws OverflowError readtsv(IOBuffer(data))
    end

    @testset "malformed data" begin
        # empty
        @test_throws TableReader.ReadError("found no column names in the header at line 1") readtsv(IOBuffer(""))

        # less columns than expected
        data = """
        col1\tcol2\tcol3
        foo\tbar
        """
        @test_throws TableReader.ReadError("unexpected number of columns at line 2") readtsv(IOBuffer(data))

        # more columns than expected
        data = """
        col1\tcol2\tcol3
        foo\tbar\tbaz\tqux\tquux
        """
        @test_throws TableReader.ReadError("unexpected number of columns at line 2") readtsv(IOBuffer(data))

        # too large floating-point number
        data = """
        col1\tcol2
        1.0\t1e500
        """
        @test_throws TableReader.ReadError("failed to parse a floating-point number") readtsv(IOBuffer(data))
    end

    @testset "UTF-8 strings" begin
        data = """
        col1\tcol2\tcol3
        甲\t乙\t丙
        👌\t😀😀😀\t🐸🐓
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == ["甲", "👌"]
        @test df[:col2] == ["乙", "😀😀😀"]
        @test df[:col3] == ["丙", "🐸🐓"]
    end

    @testset "CR+LF" begin
        data = """
        col1\tcol2\tcol3\r
        1\t2\t3\r
        4\t5\t6 \r
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, 4]
        @test df[:col2] == [2, 5]
        @test df[:col3] == [3, 6]

        data = """
        col1\tcol2\tcol3\r
        1.0\t2.0\t3.0\r
        4.0\t5.0\t6.0 \r
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1.0, 4.0]
        @test df[:col2] == [2.0, 5.0]
        @test df[:col3] == [3.0, 6.0]

        data = """
        col1\tcol2\tcol3\r
        foo\tbar\tbaz\r
        hoge\tfuga\tpiyo \r
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == ["foo", "hoge"]
        @test df[:col2] == ["bar", "fuga"]
        @test df[:col3] == ["baz", "piyo"]
    end

    @testset "CR" begin
        data = """col1\tcol2\r123\t456\r"""
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [123]
        @test df[:col2] == [456]
    end

    @testset "blank lines" begin
        # before header
        data = """

        col1\tcol2
        1\t2
        3\t4
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, 3]
        @test df[:col2] == [2, 4]

        # after header
        data = """
        col1\tcol2

        1\t2
        3\t4
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, 3]
        @test df[:col2] == [2, 4]

        # among data
        data = """
        col1\tcol2
        1\t2

        3\t4
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, 3]
        @test df[:col2] == [2, 4]

        # end of a file
        data = """
        col1\tcol2
        1\t2
        3\t4

        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, 3]
        @test df[:col2] == [2, 4]
    end

    @testset "invalid argument" begin
        @test_throws ArgumentError readtsv(IOBuffer(""), chunkbits = -1)
    end

    @testset "large data" begin
        buf = IOBuffer()
        m = 10000
        n = 2000
        println(buf, "name", '\t', join(("col$(j)" for j in 1:n), '\t'))
        for i in 1:m
            print(buf, "row$(i)")
            for j in 1:n
                print(buf, '\t', i)
            end
            println(buf)
        end
        data = take!(buf)

        # with chunking
        df = readtsv(IOBuffer(data))
        @test size(df) == (m, n + 1)
        @test df[:name] == ["row$(i)" for i in 1:m]
        @test df[:col1] == 1:m
        @test df[:col2] == 1:m

        # without chunking
        df = readtsv(IOBuffer(data), chunkbits = 0)
        @test size(df) == (m, n + 1)
        @test df[:name] == ["row$(i)" for i in 1:m]
        @test df[:col1] == 1:m
        @test df[:col2] == 1:m
    end

    @testset "from file" begin
        df = readtsv(joinpath(@__DIR__, "test.tsv"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readtsv(joinpath(@__DIR__, "test.tsv.gz"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readtsv(joinpath(@__DIR__, "test.tsv.zst"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readtsv(joinpath(@__DIR__, "test.tsv.xz"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]
    end
end

@testset "readcsv" begin
    @testset "simple" begin
        # integers
        buffer = IOBuffer("""
        col1,col2,col3
        1,23,456
        -10,-99,0
        """)
        df = readcsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1, -10]
        @test df[:col2] == [23, -99]
        @test df[:col3] == [456, 0]

        # floats
        buffer = IOBuffer("""
        col1,col2,col3
        1.0,1.1,12.34
        -1.2,0.0,-9.
        .000,.123,100.000
        1e3,1.E+123,-8.2e-00
        """)
        df = readcsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1.0, -1.2, 0.000, 1e3]
        @test df[:col2] == [1.1, 0.0, 0.123, 1e123]
        @test df[:col3] == [12.34, -9.0, 100.000, -8.2]

        # strings
        buffer = IOBuffer("""
        col1,col2,col3
        a,b,c
        foo,bar,baz
        """)
        df = readcsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == ["a", "foo"]
        @test df[:col2] == ["b", "bar"]
        @test df[:col3] == ["c", "baz"]
    end

    @testset "no header" begin
        data = """
        A,B,C
        a,b,c
        """
        df = readcsv(IOBuffer(data))
        @test df[:A] == ["a"]
        @test df[:B] == ["b"]
        @test df[:C] == ["c"]

        df = readcsv(IOBuffer(data), hasheader = false)
        @test df[:X1] == ["A", "a"]
        @test df[:X2] == ["B", "b"]
        @test df[:X3] == ["C", "c"]

        df = readcsv(IOBuffer(data), colnames = [:C1, :C2, :C3])
        @test df[:C1] == ["A", "a"]
        @test df[:C2] == ["B", "b"]
        @test df[:C3] == ["C", "c"]

        df = readcsv(IOBuffer(data), colnames = [:C1, :C2, :C3], hasheader = true)
        @test df[:C1] == ["a"]
        @test df[:C2] == ["b"]
        @test df[:C3] == ["c"]
    end

    @testset "quotation" begin
        data = """
        "col1","col2"
        "hi, there",","
        "1,2,3,4", ",,,"
        """
        df = readcsv(IOBuffer(data))
        @test df[:col1] == ["hi, there", "1,2,3,4"]
        @test df[:col2] == [",", ",,,"]
    end

    @testset "missing" begin
        data = """
        "col1","col2"
        "",""
        "1","2"
        """
        df = readcsv(IOBuffer(data))
        @test ismissing(df[1,:col1])
        @test df[2,:col1] == 1
        @test ismissing(df[1,:col2])
        @test df[2,:col2] == 2
    end

    @testset "multiline field" begin
        data = """
        col1,col2
        "oh,
        there
        there","
        multi
        line
        field
        "
        """
        df = readcsv(IOBuffer(data))
        @test df[:col1] == ["oh,\nthere\nthere"]
        @test df[:col2] == ["\nmulti\nline\nfield\n"]
    end

    @testset "leading-zero string" begin
        data = """
        col1,col2
        0000,0
        0001,1
        0002,2
        0123,123
        """
        df = readcsv(IOBuffer(data))  # lzstring = true by default
        @test df[:col1] == ["0000", "0001", "0002", "0123"]
        @test df[:col2] == [0, 1, 2, 123]

        df = readcsv(IOBuffer(data), lzstring = true)
        @test df[:col1] == ["0000", "0001", "0002", "0123"]
        @test df[:col2] == [0, 1, 2, 123]

        df = readcsv(IOBuffer(data), lzstring = false)
        @test df[:col1] == [0, 1, 2, 123]
        @test df[:col2] == [0, 1, 2, 123]

        data = """
        col1
        0.3
        00.4
        """
        df = readcsv(IOBuffer(data))  # lzstring = true by default
        @test df[:col1] == ["0.3", "00.4"]

        df = readcsv(IOBuffer(data), lzstring = true)
        @test df[:col1] == ["0.3", "00.4"]

        df = readcsv(IOBuffer(data), lzstring = false)
        @test df[:col1] == [0.3, 0.4]
    end

    @testset "EOF without newline" begin
        data = """
        col1,col2
        1,2"""
        df = readcsv(IOBuffer(data))
        @test df[:col1] == [1]
        @test df[:col2] == [2]
    end

    @testset "unnamed column" begin
        data = """
        col1,,col3,
        1,foo,3,bar
        """
        df = readcsv(IOBuffer(data))
        @test df[:col1] == [1]
        @test df[:UNNAMED_2] == ["foo"]
        @test df[:col3] == [3]
        @test df[:UNNAMED_4] == ["bar"]
    end

    @testset "implicit column" begin
        data = """
        col1,col2,col3
        foo,1,2,3
        """
        df = readcsv(IOBuffer(data))
        @test df[:UNNAMED_0] == ["foo"]
        @test df[:col1] == [1]
        @test df[:col2] == [2]
        @test df[:col3] == [3]
    end

    @testset "skip lines" begin
        data = """
        foobarbaz
        hogehogehoge
        col1,col2,col3
        1,2,3
        """
        df = readcsv(IOBuffer(data), skip = 2)
        @test df[:col1] == [1]
        @test df[:col2] == [2]
        @test df[:col3] == [3]
    end

    @testset "comment" begin
        data = """
        # comment before header
        col1,col2
        # comment after header
        1,2
        # comment among data
        3,4
        # comment after data
        """
        df = readcsv(IOBuffer(data), comment = "#")
        @test df[:col1] == [1, 3]
        @test df[:col2] == [2, 4]

        data = """
        ## comment before header
        col1,col2
        ## comment after header
        foo,bar
        ## comment among data
        #foo,#bar
        ## comment after data
        """
        df = readcsv(IOBuffer(data), comment = "##")
        @test df[:col1] == ["foo", "#foo"]
        @test df[:col2] == ["bar", "#bar"]

        @test_throws ArgumentError("comment cannot contain newline characters") readcsv(IOBuffer(data), comment = "\n")
        @test_throws ArgumentError("comment cannot contain newline characters") readcsv(IOBuffer(data), comment = "\r")
    end

    @testset "malformed data" begin
        # empty
        @test_throws TableReader.ReadError("found no column names in the header at line 1") readcsv(IOBuffer(""))

        # tab
        data = """
        col1,col2
        foo,ba\tr
        """
        @test_throws TableReader.ReadError("invalid file format at line 2, column 2 (found '\\t')") readcsv(IOBuffer(data))
    end

    @testset "from file" begin
        df = readcsv(joinpath(@__DIR__, "test.csv"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readcsv(joinpath(@__DIR__, "test.csv.gz"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readcsv(joinpath(@__DIR__, "test.csv.zst"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readcsv(joinpath(@__DIR__, "test.csv.xz"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]
    end

    @testset "from command" begin
        if Sys.which("echo") === nothing
            @info "skip tests: echo command is not found"
        else
            df = readcsv(`echo $("col1,col2\n1,2")`)
            @test df[:col1] == [1]
            @test df[:col2] == [2]
        end

        if Sys.which("cat") === nothing || Sys.which("gzip") === nothing
            @info "skip tests: cat/gzip commands are not found"
        else
            testfile = joinpath(@__DIR__, "test.csv")
            df = readcsv(pipeline(`cat $(testfile)`, `gzip`))
            @test df[:col1] == [1, 2]
            @test df[:col2] == [1.0, 2.0]
            @test df[:col3] == ["one", "two"]
        end
    end

    @testset "switching from strings to numbers" begin
        buf = IOBuffer()
        println(buf, "col1,col2")
        for i in 1:1_000_000
            if i % 1024 == 0
                println(buf, "foo,")  # missing
            else
                println(buf, "foo,bar")
            end
        end
        for i in 1:1_000_000
            println(buf, "1,2")
        end
        for i in 1:1_000_000
            println(buf, "1.1,2.2")
        end
        seekstart(buf)
        df = readcsv(buf)
        @test df[:col1] isa Vector{String}
        @test df[:col2] isa Vector{Union{String,Missing}}
    end

    @testset "repeated variable names" begin
        data = """
            col1,col1
            1,2
            """
            df = readcsv(IOBuffer(data))
            @test df[:col1_1] == [2]
    end

    @testset "no quotation" begin
        data = """
        col1,col2,col3
        "foo","bar","baz"
        """
        df = readcsv(IOBuffer(data), quot = nothing)
        @test df[:col1] == ["\"foo\""]
        @test df[:col2] == ["\"bar\""]
        @test df[:col3] == ["\"baz\""]

        data = """
        col1,\xffcol2,col3
        "foo","bar","baz"
        """
        @test_throws TableReader.ReadError("invalid file header format") readtsv(IOBuffer(data), quot = nothing)

        data = """
        col1,col2,col3
        "foo",\xff"bar","baz"
        """
        @test_throws TableReader.ReadError("invalid file format at line 2, column 1 (found 0xff)") readtsv(IOBuffer(data), quot = nothing)
    end

    @testset "datetimes" begin
        buf = IOBuffer("""
        col1,col2
        1,
        2,2018-04-10T08:19:30.000
        """)
        df = readcsv(buf)
        @test df[:col1] == [1, 2]
        @test ismissing(df[1,:col2])
        @test df[2,:col2] == DateTime(2018, 4, 10, 8, 19, 30)
    end
end

@testset "readdlm" begin
    @testset "custom parser parameters" begin
        # delimtied by |
        data = """
        col1|col2|col3
        1|2|3
        """
        df = readdlm(IOBuffer(data), delim = '|')
        @test df[:col1] == [1]
        @test df[:col2] == [2]
        @test df[:col3] == [3]

        # quoted by `
        data = """
        col1,col2,col3
        `foo`,`bar`,`baz`
        """
        df = readdlm(IOBuffer(data), delim = ',', quot = '`')
        @test df[:col1] == ["foo"]
        @test df[:col2] == ["bar"]
        @test df[:col3] == ["baz"]

        # invalid chunkbits
        @test_throws ArgumentError readdlm(IOBuffer(""), delim = ',', chunkbits = -1)
        @test_throws ArgumentError readdlm(IOBuffer(""), delim = ',', chunkbits = 100)
    end

    @testset "delimiter detection" begin
        # comma delimited
        data = """
        col1,col2,col3
        1,2,3
        """
        df = readdlm(IOBuffer(data))
        @test df[:col1] == [1]
        @test df[:col2] == [2]
        @test df[:col3] == [3]

        # tab delimited
        data = """
        col1\tcol2\tcol3
        1\t2\t3
        """
        df = readdlm(IOBuffer(data))
        @test df[:col1] == [1]
        @test df[:col2] == [2]
        @test df[:col3] == [3]

        # pipe delimited
        data = """
        col1|col2|col3
        1|2|3
        """
        df = readdlm(IOBuffer(data))
        @test df[:col1] == [1]
        @test df[:col2] == [2]
        @test df[:col3] == [3]
    end
end

@testset "colnames" begin
    @testset "normalize" begin
        # test unusual or unsafe (for dataframe column names) strings
        buffer = IOBuffer("""
        col.1,col|2,col\$3, col4, do, α,_col7
        1,2,3,4,5,6,7
        """)
        df = readcsv(buffer,normalizenames=true)
        @test names(df) == [:col_1, :col_2, :col_3,:col4,:_do,:α,:_col7]
    end
    @testset "don't normalize" begin
        buffer = IOBuffer("""
        col.1,col|2,col\$3, col4, do, α,_col7
        1,2,3,4,5,6,7
        """)
        df = readcsv(buffer,normalizenames=false)
        @test names(df) == [Symbol("col.1"), Symbol("col|2"), Symbol("col\$3"), :col4, :do, :α, :_col7]
    end
end

@testset "download" begin
    if Sys.which("curl") === nothing
        @info "Skipped file downloading tests"
    else
        df = readcsv("https://raw.githubusercontent.com/bicycle1885/TableReader.jl/master/test/test.csv")
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]
    end
end

# Field-level test cases
# ----------------------

function testfield(s::String, expected)
    data = """
    col1,col2,col3
    $(s),$(s),$(s)
      $(s)  ,  $(s)  ,  $(s)
    """
    df = readcsv(IOBuffer(data))
    x1, x2 = df[:col1]
    y1, y2 = df[:col2]
    z1, z2 = df[:col3]
    sametype = typeof(x1) == typeof(x2) ==
               typeof(y1) == typeof(y2) ==
               typeof(z1) == typeof(z2) ==
               typeof(expected)
    if expected isa Float64 && isnan(expected)
        @test sametype && isnan(x1) && isnan(x2) &&
                          isnan(y1) && isnan(y2) &&
                          isnan(z1) && isnan(z2)
    elseif expected isa Missing
        @test sametype && ismissing(x1) && ismissing(x2) &&
                          ismissing(y1) && ismissing(y2) &&
                          ismissing(z1) && ismissing(z2)
    else
        @test sametype && x1 == x2 == y1 == y2 == z1 == z2 == expected
    end
end

@testset "checkfield" begin
    @testset "integer" begin
        testfield("0", 0)
        testfield("+0", 0)
        testfield("-0", 0)
        testfield("1", 1)
        testfield("+1", 1)
        testfield("-1", -1)
        testfield("12", 12)
        testfield("+12", 12)
        testfield("-12", -12)
        testfield("1234567890", 1234567890)
        testfield("+1234567890", 1234567890)
        testfield("-1234567890", -1234567890)
        testfield("\"1234567890\"", 1234567890)
        testfield("\"+1234567890\"", 1234567890)
        testfield("\"-1234567890\"", -1234567890)
    end

    @testset "float" begin
        testfield("0.", 0.0)
        testfield("+0.", 0.0)
        testfield("-0.", -0.0)
        testfield("1.", 1.0)
        testfield("+1.", 1.0)
        testfield("-1.", -1.0)
        testfield(".0", 0.0)
        testfield("+.0", 0.0)
        testfield("-.0", -0.0)
        testfield(".1", 0.1)
        testfield("+.1", 0.1)
        testfield("-.1", -0.1)
        testfield("0.0", 0.0)
        testfield("+0.0", 0.0)
        testfield("-0.0", -0.0)
        testfield("1.0", 1.0)
        testfield("+1.0", 1.0)
        testfield("-1.0", -1.0)
        testfield("1.1", 1.1)
        testfield("+1.1", 1.1)
        testfield("-1.1", -1.1)
        testfield("123.4567890", 123.456789)
        testfield("+123.456789", 123.456789)
        testfield("-123.456789", -123.456789)
        testfield("\"-1.23\"", -1.23)
        testfield("1e3", 1e3)
        testfield("1e+3", 1e+3)
        testfield("1e-3", 1e-3)
        testfield("1E3", 1e3)
        testfield("1E+3", 1e+3)
        testfield("1E-3", 1e-3)
        testfield("+1e3", 1e3)
        testfield("-1e3", -1e3)
        testfield("123.45e67", 123.45e67)
        testfield("123.45e+67", 123.45e67)
        testfield("123.45e-67", 123.45e-67)
        testfield("\"1.2e34\"", 1.2e34)
        testfield("inf", Inf)
        testfield("iNf", Inf)
        testfield("INF", Inf)
        testfield("infinity", Inf)
        testfield("iNfInItY", Inf)
        testfield("INFINITY", Inf)
        testfield("\"inf\"", Inf)
        testfield("nan", NaN)
        testfield("nAn", NaN)
        testfield("NAN", NaN)
        testfield("\"nan\"", NaN)
    end

    @testset "bool" begin
        testfield("t", true)
        testfield("T", true)
        testfield("true", true)
        testfield("True", true)
        testfield("TrUe", true)
        testfield("TRUE", true)
        testfield("\"t\"", true)
        testfield("\"true\"", true)
        testfield("f", false)
        testfield("F", false)
        testfield("false", false)
        testfield("False", false)
        testfield("FaLsE", false)
        testfield("FALSE", false)
        testfield("\"f\"", false)
        testfield("\"false\"", false)
    end

    @testset "missing" begin
        testfield("", missing)
        testfield("\"\"", missing)
        testfield("NA", missing)
        testfield("\"NA\"", missing)
    end

    @testset "string" begin
        testfield("!", "!")
        testfield("~", "~")
        testfield(".", ".")
        testfield("..", "..")
        testfield("+", "+")
        testfield("-+", "-+")
        testfield("i", "i")
        testfield("I", "I")
        testfield("n", "n")
        testfield("N", "N")
        testfield("+12!", "+12!")
        testfield("-12!", "-12!")
        testfield("0!", "0!")
        testfield("1!", "1!")
        testfield("0..", "0..")
        testfield("1..", "1..")
        testfield("1.2.", "1.2.")
        testfield("0123", "0123")
        testfield("1e", "1e")
        testfield("1ea", "1ea")
        testfield("1e-", "1e-")
        testfield("1e+", "1e+")
        testfield("1e-!", "1e-!")
        testfield("1e+!", "1e+!")
        testfield("123 456", "123 456")
        testfield("Inf0", "Inf0")
        testfield("Infinity0", "Infinity0")
        testfield("t0", "t0")
        testfield("f0", "f0")
        testfield("true0", "true0")
        testfield("false0", "false0")
        testfield("\"a\"", "a")
        testfield("\"\"\"\"", "\"")
        testfield("\"\"\"\"\"\"", "\"\"")
    end
end

# Deprecated feature
@testset "deprecated" begin
    @info "Testing deprecated features"

    @testset "chunksize" begin
        data = """
        col1,col2
        1,2
        """
        df = readcsv(IOBuffer(data), chunksize = 0)
        @test df[:col1] == [1]
        @test df[:col2] == [2]
    end
end

ext: jl


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Seaborn

# +
using Seaborn
const sns = Seaborn
using Pandas
using PyCall: pyimport, @py_str

# Refer these issues
# https://stackoverflow.com/questions/27835619/urllib-and-ssl-certificate-verify-failed-error/42334357#42334357
# https://github.com/ytdl-org/youtube-dl/issues/4816#issuecomment-455840010
certifi = pyimport("certifi")
ENV["SSL_CERT_FILE"]=certifi.where();
# -

# ## Scatter plot
#
# - https://seaborn.pydata.org/examples/scatterplot_matrix.html
#
# ```python
# # Python
# import seaborn as sns
# sns.set(style="ticks")
#
# df = sns.load_dataset("iris")
# sns.pairplot(df, hue="species")
# ```

# ### Julia implementation

df = sns.load_dataset("iris")
head(df, 10)

pairplot(df, hue="species")

# # KDE
#
# https://seaborn.pydata.org/examples/multiple_joint_kde.html
#
# ```python
# import seaborn as sns
# import matplotlib.pyplot as plt
#
# sns.set(style="darkgrid")
# iris = sns.load_dataset("iris")
#
# # Subset the iris dataset by species
# setosa = iris.query("species == 'setosa'")
# virginica = iris.query("species == 'virginica'")
#
# # Set up the figure
# f, ax = plt.subplots(figsize=(8, 8))
# ax.set_aspect("equal")
#
# # Draw the two density plots
# ax = sns.kdeplot(setosa.sepal_width, setosa.sepal_length,
#                  cmap="Reds", shade=True, shade_lowest=False)
# ax = sns.kdeplot(virginica.sepal_width, virginica.sepal_length,
#                  cmap="Blues", shade=True, shade_lowest=False)
#
# # Add labels to the plot
# red = sns.color_palette("Reds")[-2]
# blue = sns.color_palette("Blues")[-2]
# ax.text(2.5, 8.2, "virginica", size=16, color=blue)
# ax.text(3.8, 4.5, "setosa", size=16, color=red)
# ```

# ### Julia implementation

sns.set(style="darkgrid")
iris = sns.load_dataset("iris")
setosa = query(iris, :(species=="setosa"))
virginica = query(iris, :(species=="virginica"))
fig, ax = sns.subplots(figsize=(8,8))
ax = sns.kdeplot(setosa.sepal_width, setosa.sepal_length,
            cmap="Reds", shade=true, shade_lowest=false)
ax = sns.kdeplot(virginica.sepal_width, virginica.sepal_length,
                 cmap="Blues", shade=true, shade_lowest=false)
# Add labels to the plot
red = sns.seaborn.color_palette("Reds")[end-2]
blue = sns.seaborn.color_palette("Blues")[end-2]
ax.text(2.5, 8.2, "virginica", size=16, color=blue)
ax.text(3.8, 4.5, "setosa", size=16, color=red)


text:
using Dates
using TimeZones

struct DateWrapper
    days::Integer
end

Base.:(==)(x::DateWrapper, y::DateWrapper) = x.days == y.days

function DateWrapper(y, m, d)
    days = (Date(y,m,d)-Date(1970,1,1)).value
    return DateWrapper(days)
end

# timezone_offset here is in seconds
struct TimeWrapper
    nanoseconds::Integer
    timezone_offset::Union{Integer, Nothing}
end

Base.:(==)(x::TimeWrapper, y::TimeWrapper) = x.nanoseconds == y.nanoseconds && x.timezone_offset == y.timezone_offset

# if o is an integer then it is timezone_offset in mins
function TimeWrapper(h, m, s, n=0; timezone=nothing)
    nanoseconds = n + s*1e+9 + m*6e+10 + h*3.6e+12
    return TimeWrapper(nanoseconds, if isa(timezone, Number) Integer(floor(timezone*60)) else timezone end)
end

TimeWrapper(nanoseconds) = TimeWrapper(nanoseconds, nothing)

struct DateTimeWrapper
    seconds::Integer
    nanoseconds::Integer
    timezone_offset::Union{Integer, Nothing}
    timezone_name::Union{String, Nothing}
end

Base.:(==)(x::DateTimeWrapper, y::DateTimeWrapper) = x.seconds == y.seconds && x.nanoseconds == y.nanoseconds && x.timezone_offset == y.timezone_offset && x.timezone_name == y.timezone_name

DateTimeWrapper(seconds, nanoseconds) = DateTimeWrapper(seconds, nanoseconds, nothing, nothing)
DateTimeWrapper(seconds, nanoseconds, timezone::Integer) = DateTimeWrapper(seconds, nanoseconds, timezone, nothing)
DateTimeWrapper(seconds, nanoseconds, timezone::String) = DateTimeWrapper(seconds, nanoseconds, nothing, timezone)

function DateTimeWrapper(y,m,d,h,i,s=0,n=0;timezone=nothing)
    millisecond = (DateTime(y,m,d,h,i,s)-DateTime(1970,1,1)).value
    seconds = Integer(floor(millisecond/1000))
    nanoseconds = n

    if timezone == nothing
        return DateTimeWrapper(seconds, nanoseconds, nothing, nothing)
    elseif isa(timezone, String)
        return DateTimeWrapper(seconds, nanoseconds, nothing, timezone)
    elseif isa(timezone, Number)
        return DateTimeWrapper(seconds, nanoseconds, Integer(floor(timezone*60)), nothing)
    else
        throw(ErrorException("Not a valid DateTime"))
    end
end

struct DurationWrapper
    months::Integer
    days::Integer
    seconds::Integer
    nanoseconds::Integer
end

Base.:(==)(x::DurationWrapper, y::DurationWrapper) = x.months == y.months && x.days == y.days && x.seconds == y.seconds && x.nanoseconds == y.nanoseconds

function DurationWrapper(years, months, days, hours, minutes, seconds=0, nanoseconds=0)
    m = months + years * 12
    d = days
    s = seconds + hours * 3600 + minutes * 60
    n = nanoseconds
    return DurationWrapper(m, d, s, n)
end

DurationWrapper(;years=0, months=0, days=0, hours=0, minutes=0, seconds=0, nanoseconds=0) = DurationWrapper(years, months, days, hours, minutes, seconds, nanoseconds)

function hydrate_date(days; use_julia_dates=true)
    if use_julia_dates
        return Date(1970,1,1) + Dates.Day(days)
    else
        return DateWrapper(days)
    end
end

function dehydrate_date(value::Date)
    return Structure(UInt8('D'), Any[(value-Date(1970,1,1)).value])
end

function dehydrate_date(value::DateWrapper)
    return Structure(UInt8('D'), Any[value.days])
end

# Julia Times do not support time zones
function hydrate_time(nanoseconds::Integer, timezone::Union{Integer, Nothing}=nothing; use_julia_dates=true)
    if use_julia_dates
        return Time(Dates.Nanosecond(nanoseconds))
    else
        return TimeWrapper(nanoseconds, timezone)
    end
end

function dehydrate_time(value::Time)
    return Structure(UInt8('t'), Any[value.instant.value])
end

function dehydrate_time(value::TimeWrapper)
    if value.timezone_offset != nothing
        return Structure(UInt8('T'), Any[value.nanoseconds, value.timezone_offset])
    else
        return Structure(UInt8('t'), Any[value.nanoseconds])
    end
end

# Julia DateTimes do not support nanosecond accuracy time
function hydrate_datetime(seconds::Integer, nanoseconds::Integer, timezone::Union{String, Integer, Nothing}=nothing; use_julia_dates=true)
    if use_julia_dates
        date_time =  DateTime(1970,1,1) + Dates.Second(seconds) + Dates.Nanosecond(nanoseconds)
        if timezone == nothing
            return date_time
        elseif isa(timezone, String)
            return ZonedDateTime(date_time, TimeZone(timezone))
        else
            return ZonedDateTime(date_time, FixedTimeZone("Time Zone", timezone))
        end
    else
        if timezone == nothing
            return DateTimeWrapper(seconds, nanoseconds, nothing, nothing)
        elseif isa(timezone, String)
            return DateTimeWrapper(seconds, nanoseconds, nothing, timezone)
        else
            return DateTimeWrapper(seconds, nanoseconds, timezone, nothing)
        end
    end
end

function dehydrate_datetime(value::DateTime)
    millisecond = (value-DateTime(1970,1,1)).value # miliseconds
    seconds = Integer(floor(millisecond/1000.0))
    nanoseconds = (millisecond % 1000) * 1000000
    return Structure(UInt8('d'), Any[seconds, nanoseconds])
end

function dehydrate_datetime(value::ZonedDateTime)
    d = TimeZones.localtime(value)
    millisecond = (d-DateTime(1970,1,1)).value # miliseconds
    seconds = Integer(floor(millisecond/1000.0))
    nanoseconds = (millisecond % 1000) * 1000000
    if isa(value.timezone, FixedTimeZone)
        return Structure(UInt8('F'), Any[seconds, nanoseconds, value.timezone.offset.std.value])
    else
        return Structure(UInt8('f'), Any[seconds, nanoseconds, string(value.timezone.name)])
    end
end

function dehydrate_datetime(value::DateTimeWrapper)
    if value.timezone_offset != nothing
        return Structure(UInt8('F'), Any[value.seconds, value.nanoseconds, value.timezone_offset])
    elseif value.timezone_name != nothing
        return Structure(UInt8('f'), Any[value.seconds, value.nanoseconds, value.timezone_name])
    else
        return Structure(UInt8('d'), Any[value.seconds, value.nanoseconds])
    end

end

function hydrate_duration(months, days, seconds, nanoseconds; use_julia_dates=true)
    if use_julia_dates
        return Month(months) + Day(days) + Second(seconds) + Nanosecond(nanoseconds)
    else
        return DurationWrapper(months, days, seconds, nanoseconds)
    end
end

function tons(c::Dates.CompoundPeriod)
    s::Float64 = 0.0
    for p in c.periods
        v = Dates.tons(p)
        s += v
    end
    return s
end

function dehydrate_duration(value::Dates.CompoundPeriod)
    value = canonicalize(value)

    months = 0
    periods = Dates.Period[]

    for p in value.periods
        if isa(p, Dates.Month)
            months += p.value
        elseif isa(p, Dates.Year)
            months += p.value * 12
        else
            push!(periods, p)
        end
    end

    value = Dates.CompoundPeriod(periods)
    v = tons(value)

    nanoseconds = v % 1000000000
    v = (v - nanoseconds) / 1000000000
    seconds = v % 86400
    v = (v - seconds) / 86400
    days = v

    return Structure(UInt8('E'), Any[Integer(months), Integer(days), Integer(seconds), Integer(nanoseconds)])
end

dehydrate_duration(value::Period) = dehydrate_duration(Dates.CompoundPeriod(value))

function dehydrate_duration(value::DurationWrapper)
    return Structure(UInt8('E'), Any[value.months, value.days, value.seconds, value.nanoseconds])
end

dehydrate(value::DateWrapper) = dehydrate_date(value)
dehydrate(value::TimeWrapper) = dehydrate_time(value)
dehydrate(value::DateTimeWrapper) = dehydrate_datetime(value)
dehydrate(value::DurationWrapper) = dehydrate_duration(value)
dehydrate(value::Dates.Date) = dehydrate_date(value)
dehydrate(value::Dates.Time) = dehydrate_time(value)
dehydrate(value::Dates.DateTime) = dehydrate_datetime(value)
dehydrate(value::ZonedDateTime) = dehydrate_datetime(value)
dehydrate(value::Dates.CompoundPeriod) = dehydrate_duration(value)
dehydrate(value::Period) = dehydrate_date(value)

function temporal_hydration_functions()
    return Dict(
        UInt8('D') => (psh, values)->hydrate_date(values..., use_julia_dates=psh.use_julia_dates),
        UInt8('T') => (psh, values)->hydrate_time(values..., use_julia_dates=psh.use_julia_dates),
        UInt8('t') => (psh, values)->hydrate_time(values..., use_julia_dates=psh.use_julia_dates),
        UInt8('F') => (psh, values)->hydrate_datetime(values..., use_julia_dates=psh.use_julia_dates),
        UInt8('f') => (psh, values)->hydrate_datetime(values..., use_julia_dates=psh.use_julia_dates),
        UInt8('d') => (psh, values)->hydrate_datetime(values..., use_julia_dates=psh.use_julia_dates),
        UInt8('E') => (psh, values)->hydrate_duration(values..., use_julia_dates=psh.use_julia_dates)
    )
end


text:
# Implementation of the Next Subvolume Method

############################ NSM ###################################
#NOTE state vector u is a matrix. u[i,j] is species i, site j
#NOTE hopping_constants is a matrix. hopping_constants[i,j] is species i, site j
mutable struct NSMJumpAggregation{J,T,RX,HOP,RNG,DEPGR,VJMAP,JVMAP,PQ,SS} <: AbstractSSAJumpAggregator
    next_jump::SpatialJump{J} #some structure to identify the next event: reaction or hop
    prev_jump::SpatialJump{J} #some structure to identify the previous event: reaction or hop
    next_jump_time::T
    end_time::T
    rx_rates::RX
    hop_rates::HOP
    # rates::F1 #rates for constant-rate jumps
    # affects!::F2 #affects! function determines the effect of constant-rate jumps
    save_positions::Tuple{Bool,Bool}
    rng::RNG
    dep_gr::DEPGR #dep graph is same for each site
    vartojumps_map::VJMAP #vartojumps_map is same for each site
    jumptovars_map::JVMAP #jumptovars_map is same for each site
    pq::PQ
    spatial_system::SS
    numspecies::Int #number of species
end

function NSMJumpAggregation(nj::SpatialJump{J}, njt::T, et::T, rx_rates::RX, hop_rates::HOP, sps::Tuple{Bool,Bool},
                                      rng::RNG, spatial_system::SS; num_specs, vartojumps_map=nothing, jumptovars_map=nothing, dep_graph=nothing, kwargs...) where {J,T,RX,HOP,RNG,SS}

    # a dependency graph is needed
    if dep_graph === nothing
        dg = make_dependency_graph(num_specs, rx_rates.ma_jumps)
    else
        dg = dep_graph
        # make sure each jump depends on itself
        add_self_dependencies!(dg)
    end

    # a species-to-reactions graph is needed
    if vartojumps_map === nothing
        vtoj_map = var_to_jumps_map(num_specs, rx_rates.ma_jumps)
    else
        vtoj_map = vartojumps_map
    end

    if jumptovars_map === nothing
        jtov_map = jump_to_vars_map(rx_rates.ma_jumps)
    else
        jtov_map = jumptovars_map
    end

    pq = MutableBinaryMinHeap{T}()

    NSMJumpAggregation{J,T,RX,HOP,RNG,typeof(dg),typeof(vtoj_map),typeof(jtov_map),typeof(pq),SS}(nj, nj, njt, et, rx_rates, hop_rates, sps, rng, dg, vtoj_map, jtov_map, pq, spatial_system, num_specs)
end

############################# Required Functions ##############################
# creating the JumpAggregation structure (function wrapper-based constant jumps)
function aggregate(aggregator::NSM, starting_state, p, t, end_time, constant_jumps, ma_jumps, save_positions, rng; hopping_constants, spatial_system, kwargs...)
    num_species = size(starting_state,1)
    majumps = ma_jumps
    if majumps === nothing
        majumps = MassActionJump(Vector{typeof(end_time)}(), Vector{Vector{Pair{Int,Int}}}(), Vector{Vector{Pair{Int,Int}}}())
    end

    next_jump = SpatialJump{Int}(typemax(Int),typemax(Int),typemax(Int)) #a placeholder
    next_jump_time = typemax(typeof(end_time))
    rx_rates = RxRates(num_sites(spatial_system), majumps)
    hop_rates = HopRates(hopping_constants, spatial_system)

    NSMJumpAggregation(next_jump, next_jump_time, end_time, rx_rates, hop_rates, save_positions, rng, spatial_system; num_specs = num_species, kwargs...)
end

# set up a new simulation and calculate the first jump / jump time
function initialize!(p::NSMJumpAggregation, integrator, u, params, t)
    p.end_time = integrator.sol.prob.tspan[2]
    fill_rates_and_get_times!(p, u, t)
    generate_jumps!(p, integrator, params, u, t)
    nothing
end

# calculate the next jump / jump time
function generate_jumps!(p::NSMJumpAggregation, integrator, params, u, t)
    p.next_jump_time, site = top_with_handle(p.pq)
    p.next_jump_time >= p.end_time && return nothing
    p.next_jump = sample_jump_direct(p, site)
    nothing
end

# execute one jump, changing the system state
function execute_jumps!(p::NSMJumpAggregation, integrator, u, params, t)
    # execute jump
    update_state!(p, integrator)

    # update current jump rates and times
    update_dependent_rates_and_firing_times!(p, integrator.u, t)
    nothing
end

######################## SSA specific helper routines ########################
"""
    fill_rates_and_get_times!(aggregation::NSMJumpAggregation, u, t)

reset all stucts, reevaluate all rates, recalculate tentative site firing times, and reinit the priority queue
"""
function fill_rates_and_get_times!(aggregation::NSMJumpAggregation, u, t)
    @unpack spatial_system, rx_rates, hop_rates, rng = aggregation

    reset!(rx_rates)
    reset!(hop_rates)

    rxs = 1:num_rxs(rx_rates)
    species = 1:aggregation.numspecies
    num_nodes = num_sites(spatial_system)

    pqdata = Vector{typeof(t)}(undef, num_nodes)
    @inbounds for site in 1:num_nodes
        update_rx_rates!(rx_rates, rxs, u, site)
        update_hop_rates!(hop_rates, species, u, site, spatial_system)
        pqdata[site] = t + randexp(rng) / (total_site_rate(rx_rates, hop_rates, site))
    end

    aggregation.pq = MutableBinaryMinHeap(pqdata)
    nothing
end

"""
    update_dependent_rates_and_firing_times!(p, u, t)

recalculate jump rates for jumps that depend on the just executed jump (p.prev_jump)
"""
function update_dependent_rates_and_firing_times!(p::NSMJumpAggregation, u, t)
    jump = p.prev_jump
    if is_hop(p, jump)
        source_site = jump.src
        target_site = jump.dst
        update_rates_after_hop!(p, u, source_site, target_site, jump.jidx)
        update_site_time!(p, source_site, t)
        update_site_time!(p, target_site, t)
    else
        site = jump.src
        update_rates_after_reaction!(p, u, site, reaction_id_from_jump(p,jump))
        update_site_time!(p, site, t)
    end
    nothing
end

"""
    update_site_time!(p::NSMJumpAggregation, site, t)

update the time of site in the priority queue
"""
function update_site_time!(p::NSMJumpAggregation, site, t)
    site_rate = (total_site_rate(p.rx_rates, p.hop_rates, site))
    if site_rate > zero(typeof(site_rate))
        update!(p.pq, site, t + randexp(p.rng) / site_rate)
    else
        update!(p.pq, site, typemax(t))
    end
    nothing
end

"""
    num_constant_rate_jumps(aggregator::NSMJumpAggregation)

number of constant rate jumps
"""
num_constant_rate_jumps(aggregator::NSMJumpAggregation) = 0


text:
"""
`optimalPLD_solveSAO_LP_Benders( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, x_initial, A_feas, b_feas)`

Solve the integer version of the Bertsimas and Mišić first-choice PLD formulation using Benders decomposition.

This function will formulate the Benders reformulation in Gurobi by specifying a lazy constraint callback to generate the Benders cuts. Cuts are generated using the closed-form expressions for the Benders cuts when the x_i's are integer (see Section 4.1 of Bertsimas and Mišić).

`n` = number of candidate products.

`revenues` = array of n numbers representing marginal revenue/profit of each of the n products. (These are the pi_i's in the Bertsimas and Mišić paper.)

`numPermutations` = number of rankings / customer types.

`lambda` = array specifying the probability of each ranking / customer type; should be nonnegative and sum to 1.

`orderings` = 2D array of integers encoding each ranking. orderings[k,j] contains the index of ranking k's jth top product; (n+1) encodes the no-purchase option. For example, if n = 5 and

	orderings[k,:] = [3,5,2,6,1,4]

then the customer prefers 3 < 5 < 2 < 0 (no-purchase option).

`MIP_timeLimit` = maximum time limit for Gurobi to solve the integer problem. If set to 0, no time limit is imposed.

`x_initial` = array with with n values (between 0 and 1) indicating the initial values of the x_i variables to use to warm-start the Benders decomposition method. If empty, no warm-start solution will

`A_feas`, `b_feas` = 2D array and 1D array specifying the system of linear inequalities Ax <= b for the PLD problem. This can be used to impose side constraints (this is the Cx <= d constraint system in the formulation in Bertsimas and Mišić 2018). If b_feas is empty, then no constraints will be added (i.e., it will solve the unconstrained problem).

`rankingid_list`, `alpha_list`, `beta_list`, `gamma_list` = 1D, 2D, 2D and 1D arrays containing information on cuts generated in Benders LP solution process. These cuts are added to the Benders integer formulation so as to warm-start it.

Returns:

`x_val` = array of optimal x_i variables for the integer problem.

`expRevenue` = expected per-customer revenue/profit of the integer solution.

`upperBound` = best upper bound on the optimal per-customer revenue/profit.

`solutiontime` = time required to solve the problem.

`MIPgap` = optimality gap of the best integer solution found.

`total_num_lazy_checked` = number of lazy constraints that are checked (incremented by 1 for each customer type and each invocation of the lazy constraint callback)

`total_num_lazy_generated` = number of lazy constraints added to the formulation (incremented by 1 every time a lazy constraint for any customer type is found to be violated and is added to the formulation)

"""
function optimalPLD_solveSAO_Benders( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial,
								A_feas, b_feas,
								rankingid_list, alpha_list, beta_list, gamma_list )

	maxRev = maximum(revenues);

	if (MIP_timeLimit > 0)
		m = Model( solver = GurobiSolver(TimeLimit = MIP_timeLimit))
	else
		m = Model( solver = GurobiSolver())
	end

	if (!isRelaxation)
		@variable(m, x[1:n], Bin)
	else
		@variable(m, 0 <= x[1:n] <= 1)
	end

	if (!isRelaxation)
		if (!isempty(x_initial))
			for i=1:n
				setvalue(x[i], x_initial[i]);
			end
		end
	end

	@variable(m, 0 <= subobj[1:numPermutations] <= maxRev)

	for t=1:length(b_feas)
		@constraint(m, sum( A_feas[t,i]*x[i] for i=1:n) <= b_feas[t])
	end

	for t=1:length(rankingid_list)
		k = rankingid_list[t];
		@constraint(m, subobj[k] <= gamma_list[t] + sum( alpha_list[t,i]*x[i] for i=1:n) + sum(beta_list[t,i]*(1-x[i]) for i=1:n) )
	end

	@objective(m, Max, sum( lambda[k] * subobj[k] for k=1:numPermutations))

	const total_num_lazy_generated = 0;
	const total_num_lazy_checked = 0;

	orderings = convert( Array{Int64}, orderings);
	p2r_mat = zeros(Int64, numPermutations, n+1);
	temp_rev = [revenues;0];
	revenues_ordered_mat = zeros(numPermutations, n+1);
	sorted_rev_inds_mat = zeros(Int64, numPermutations, n+1);

	for k in 1:numPermutations
		single_ordering = convert(Array{Int64}, orderings[k,:][:]);
		revenues_ordered =  temp_rev[single_ordering[:]]
		revenues_ordered_mat[k,:] = revenues_ordered
		sorted_rev_inds_mat[k,:] = sortperm(revenues_ordered[:], rev= true);
		p2r_mat[k,:] = sortperm(single_ordering[:]); # single_ordering is r2p
	end

	tol_integrality = 1e-4;
	# Compute current lower bound.
	function benders_cb(cb)
		x_val = getvalue(x)[:]

		sol_is_integral = all( i->( (i>1-tol_integrality) &  (i<tol_integrality)), x_val )

		# println("candSet = $candSet")
		subobj_val = getvalue(subobj)[:]

		candSet_ponly = find(x_val .> 0.5)
		candSet = [candSet_ponly; n+1];
		for k=1:numPermutations
			total_num_lazy_checked = total_num_lazy_checked + 1;
			#p2r = sortperm(orderings[k,:][:])
			gamma_val = 0;
			alpha_val = zeros(n);
			beta_val = zeros(n+1);
			#temp = indmin(p2r_mat[k,candSet[:]])
			temp = indmin(p2r_mat[k,candSet])
			if (candSet[temp] <= n)
				gamma_val = revenues[candSet[temp]];
			else
				gamma_val = 0.0
			end
			# println("\tgamma_val is ")
			# println(gamma_val)

			if (!isempty(candSet_ponly))
				maxrev = maximum(revenues[candSet_ponly[:]]);
				beta_val[candSet[temp]] = maxrev - gamma_val
			end
			# println("\tbeta_val[temp] is (temp = $temp) ")
			# println(beta_val[candSet[temp]])

			#alpha_val = zeros(size(alpha_val))
			for i2=1:(n+1)
				if (orderings[k,i2] <= n)
					alpha_val[orderings[k,i2]] = max(revenues[orderings[k,i2]] - gamma_val - sum( beta_val[orderings[k,1:(i2-1)]]), 0 )
				end
			end

			if (gamma_val < subobj_val[k])
				# println("\tperm $k: gamma_val = $gamma_val, subobj_val[k] = $(subobj_val[k])")
				@lazyconstraint(cb, gamma_val + sum( alpha_val[i] * x[i] for i=1:n) + sum( beta_val[i] * (1 - x[i]) for i=1:n) >= subobj[k])
				total_num_lazy_generated = total_num_lazy_generated+1;
			end
		end
	end

	addlazycallback(m, benders_cb)

	solve(m)

	println(getvalue(x))

	x_val = getvalue(x)[:]

	solutiontime = Gurobi.get_runtime(m.internalModel.inner);
	if (!isRelaxation)
		MIPgap = Gurobi.get_dblattr(m.internalModel.inner, "MIPGap");
	else
		MIPgap = -1;
	end

	expRevenue = getobjectivevalue(m)

	upperBound = Gurobi.get_dblattr(m.internalModel.inner, "ObjBound");

	return x_val, expRevenue, upperBound, solutiontime, MIPgap, total_num_lazy_checked, total_num_lazy_generated
end

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
# using StaticArrays

# Define Rosenbrock function: https://en.wikipedia.org/wiki/Rosenbrock_function
#
# This function can be interpreted as a NLS pb:
# r_1 = (1-θ₁)
# r_2 = 10(θ₂-θ₁²)
#
# Jacobian is: [[-1, 0],[-20, 10]]
#
# Minimum for θ=(1,1) where f=1/2 r^2 = 0
#
# TODO: implement even n (here 2) generalization of the wiki page.
#
struct Rosenbrock <: AbstractNLS
end

parameter_size(::Rosenbrock) = 2
residue_size(::Rosenbrock) = 2

function eval_r(nls::Rosenbrock,θ::AbstractVector{T}) where T
    @assert length(θ)==parameter_size(nls)

    T[ 1-θ[1], 10*(θ[2]-θ[1]^2) ]
end

function eval_r_J(nls::Rosenbrock,θ::AbstractVector{T}) where T
    @assert length(θ)==parameter_size(nls)

    r = T[ 1-θ[1], 10*(θ[2]-θ[1]^2) ]
    J = T[ -1 0; -20*θ[1] 10]

    (r,J)
end

# ================================================================

const Rosenbrock_Static = Rosenbrock

# # Define the same Rosenbrock function but using Statics Arrays... This
# # cause some problems: https://github.com/JuliaArrays/StaticArrays.jl/issues/971
# #
# struct Rosenbrock_Static <: AbstractNLS
# end

# parameter_size(::Rosenbrock_Static) = 2
# residue_size(::Rosenbrock_Static) = 2

# function eval_r(nls::Rosenbrock_Static,θ::AbstractVector{T}) where T
#     @assert length(θ)==parameter_size(nls)

#     r = @SVector T[1-θ[1],
#                    10*(θ[2]-θ[1]^2)]
# end

# function eval_r_J(nls::Rosenbrock_Static,θ::AbstractVector{T}) where T
#     @assert length(θ)==parameter_size(nls)

#     r = eval_r(nls,θ)

#     # CAVEAT: SMatrix are filled column by column
#     #
#     # The Jacobian
#     #
#     # | ∂1r1, ∂2r1 |
#     # | ∂1r2, ∂2r2 |
#     #
#     # must be created by
#     #
#     # SMatrix(∂1r1,∂1r2, # col 1
#     #         ∂2r1,∂2r2) # col 2
#     #
#     # which is not intuitive and error prone.
#     #
#     # The alternative to is to use the @SMatrix macro
#     # that follows the "natural" order (row by row)
#     #
#     J = @SMatrix T[      -1   +0;     # ∂1r1, ∂2r1
#                          -20*θ[1]  +10]     # ∂1r2, ∂2r2

#     (r, J)
# end


text:
type GaussianMixture
    K::Integer
    max_iters::Integer
    method::String
    tolerance::Float64
    weights::Vector
    means::Matrix
    cov_::Dict{Integer,Matrix}
    responsibilities::Matrix
    likelihood::Vector
end

function GaussianMixture(;
                         K::Integer = 4,
                         max_iters::Integer = 500,
                         tolerance::Float64 = 1e-3,
                         method::String = "random",
                         weights::Vector = zeros(4),
                         means::Matrix = zeros(4,4),
                         cov_::Dict{Integer,Matrix} = Dict{Integer,Matrix}(),
                         likelihood::Vector = zeros(4),
                         responsibilities::Matrix = zeros(4,4))
    return GaussianMixture(K, max_iters, method,
     tolerance, weights, means, cov_, responsibilities, likelihood)
end

function train!(model::GaussianMixture, X::Matrix)
    n_sample = size(X,1)
    initialize_(model, X)
    for i = 1:model.max_iters
        E_step!(model,X)
        M_step!(model,X)
        if isconverge_(model)
            break
        end
    end
end

function initialize_(model, X)
    model.weights = ones(model.K) * (1/model.K)
    model.means = X[randperm(size(X,1))[1:model.K],:]'
    rand_cov = cov(X)
    for i = 1:model.K
        model.cov_[i] = rand_cov
    end
    model.likelihood = Vector{Real}()
end

function E_step!(model, X)
    likelihood = zeros(size(X,1),model.K)
    for i = 1:model.K
        model.cov_[i] = Symmetric(model.cov_[i],:L)
        dis = MvNormal(model.means[:,i], model.cov_[i])
        likeli = pdf(dis, X')
        likelihood[:, i] = likeli
    end
    weighted_likelihood = repmat(model.weights',size(X,1),1) .* likelihood
    temp1 = sum(log(sum(weighted_likelihood,2)))
    push!(model.likelihood, temp1)
    for i = 1:size(weighted_likelihood, 1)
        weighted_likelihood[i,:] = weighted_likelihood[i,:] / sum(weighted_likelihood[i,:])
    end
    model.responsibilities = weighted_likelihood
end

function M_step!(model, X)
    n_sample = size(X,1)
    n_feature = size(X,2)
    for i = 1:model.K
        resp = model.responsibilities[:,i]
        resp = resp'
        model.weights[i] = sum(resp)/n_sample
        for j = 1:n_feature
            temp =  resp * X[:,j]/sum(resp)
            model.means[j, i] =  temp[1]
        end
        cov_ = 0
        for j = 1:n_sample
            temp = (X[j,:] - model.means[:,i])
            cov_ += resp[j] .* temp * temp'
        end
        model.cov_[i] = cov_ ./ sum(resp)
    end
end

function isconverge_(model::GaussianMixture)
    if length(model.likelihood) > 1 && model.likelihood[end] - model.likelihood[end-1] <= model.tolerance
        return true
    end
    return false
end
function predict(model::GaussianMixture,
                 x::Matrix)
    n = size(x,1)
    res = zeros(n)
    for i = 1:n
        res[i] = predict(model, x[i,:])
    end
    return res
end

function predict(model::GaussianMixture,
                 x::Vector)
    likelihood = zeros(model.K)
    for i = 1:model.K
        model.cov_[i] = Symmetric(model.cov_[i],:L)
        dis = MvNormal(model.means[:,i], model.cov_[i])
        likeli = pdf(dis, x)
        likelihood[i] = likeli
    end
    weighted_likelihood = model.weights .* likelihood
    return indmax(weighted_likelihood)
end

# maybe some problem

function test_GaussianMixture()
    X_train, y_test= make_blo()
    clu = size(X_train, 2)
    model = GaussianMixture(K = clu)
    train!(model,X_train)
    predictions = predict(model,X_train)
end


text:
export lexicase_evaluate

function mse_valid(h::AbstractArray, y::AbstractArray; margin::Float64=0.1)
    sum((h - y).^2) < margin
end

function classify_valid(h::AbstractArray, y::AbstractArray)
    argmax(h) == argmax(y)
end

"""
lexicase selection for data-based GP problems where each genome G in
e.population corresponds to a data-processing function F, where F =
interpret(G). X and Y are 2D arrays, and individuals are evaluated based on
valid(F(X[:, i]), Y[:, i]), where the valid function determines if individuals
should continue to be evaluated.
"""
function lexicase_evaluate!(e::AbstractEvolution, X::AbstractArray, Y::AbstractArray,
                            interpret::Function; valid::Function=classify_valid,
                            verify_best::Bool=true, seed::Int=e.gen)
    Random.seed!(seed)
    npop = length(e.population)
    ndata = size(Y, 2)
    functions = Array{Function}(undef, npop)
    for i in 1:npop
        functions[i] = interpret(e.population[i])
    end
    data_inds = shuffle!(collect(1:ndata))
    is_valid = trues(npop)
    fits = zeros(npop)
    # if verify_best
    #     fits = [e.population[i].fitness[1] for i in 1:npop]
    #     is_valid = fits .== -Inf # do not re-evaluate expert
    # end
    for i in 1:ndata
        x = X[:, data_inds[i]]
        y = Y[:, data_inds[i]]
        for j in 1:npop
            if is_valid[j]
                h = functions[j](x)
                if valid(h, y)
                    fits[j] = i
                else
                    is_valid[j] = false
                end
            end
        end
        if sum(is_valid) <= 1
            break
        end
    end
    for i in 1:npop
        e.population[i].fitness[1] = fits[i]
    end
end

#     if verify_best
#         best = argmax(fits)
#         if e.population[best] == -Inf
#             # new best individual, re-evaluate
#             fits[best] = 0.0
#             for i in 1:ndata
#                 for j in 1:npop
#                     if valid(functions[j](X[:, i]), Y[:, i])
#                         fits[best] += 1.0
#                     end
#                 end
#             end
#         end
#     end
# end


text:
export build_momentrelaxation

"""
    momentrelaxation = build_momentrelaxation(relax_ctx, problem, moment_param::Dict{String, Tuple{Set{String}, Int}}, max_cliques::Dict{String, Set{Variable}})

    Compute the `momentrelaxation` of `problem` corresponding to the clique decomposition `max_cliques` and parameters `moment_param`.
"""
function build_momentrelaxation(relax_ctx::RelaxationContext,
                                problem::Problem,
                                momentmat_param::Dict{String, Int},
                                localizingmat_param::Dict{String, Tuple{Set{String}, Int}},
                                max_cliques::Dict{String, Set{Variable}};
                                T=Float64)

    var_to_cliques = Dict{Variable, Set{String}}()
    for (clique, vars) in max_cliques
        for var in vars
            haskey(var_to_cliques, var) || (var_to_cliques[var] = Set{String}())
            push!(var_to_cliques[var], clique)
        end
    end

    ## Building linear-in-moments objective
    objective = Dict{Moment, T}()
    for (expo, val) in problem.objective
        clique = get_exponentclique(expo, var_to_cliques)
        objective[Moment(expo, clique)] = val
    end

    ## Building linear matrix inequalities
    momentmatrices = Dict{Tuple{String, String}, MomentMatrix{T}}()

    ## Build moment matrix
    for (cliquename, vars) in max_cliques
        dcl = momentmat_param[cliquename]
        momentmatrices[(get_momentcstrname(), cliquename)] = MomentMatrix{T}(relax_ctx, vars, dcl, relax_ctx.symmetries,
                                                                                relax_ctx.cstrtypes[get_momentcstrname()],
                                                                                default_clique = cliquename)
    end

    ## Build localizing matrices
    for (cstrname, cstr) in problem.constraints

        cstrtype = get_cstrtype(cstr)
        if cstrtype == :ineqdouble
            cstrname_lo, cstrname_up = get_cstrname_lower(cstrname), get_cstrname_upper(cstrname)

            # Deal with lower inequality
            clique_keys, order = localizingmat_param[cstrname_lo]
            vars, cliquename = collect_cliquesvars(clique_keys, max_cliques)
            # length(clique_keys) == 1 || error(LOGGER, "build_momentrelaxation(): constraint $cstrname spans several cliques ($clique_keys).\nNot supported yet.")

            mmt = MomentMatrix{T}(relax_ctx, vars, order, relax_ctx.symmetries,
                                                          relax_ctx.cstrtypes[cstrname_lo],
                                                          var_to_cliques = var_to_cliques)
            # print_with_color(:green, "$cstrname, :Lo\n") ##NOTE: find better logging system.
            product!(mmt, cstr.p - cstr.lb, var_to_cliques)
            momentmatrices[(cstrname_lo, cliquename)] = mmt

            # Deal with upper inequality
            clique_keys, order = localizingmat_param[cstrname_up]
            vars, cliquename = collect_cliquesvars(clique_keys, max_cliques)
            # length(clique_keys) == 1 || error(LOGGER, "build_momentrelaxation(): constraint $cstrname spans several cliques ($clique_keys).\nNot supported yet.")

            mmt = MomentMatrix{T}(relax_ctx, vars, order, relax_ctx.symmetries,
                                                          relax_ctx.cstrtypes[cstrname_up],
                                                          var_to_cliques = var_to_cliques)
            # print_with_color(:green, "$cstrname, :Up\n") ##NOTE: find better logging system.
            product!(mmt, cstr.ub - cstr.p, var_to_cliques)
            momentmatrices[(cstrname_up, cliquename)] = mmt

            # # Deal with upper inequality, no recomputing of variables or moment matrix if possible
            # clique_keys_up, order_up = localizingmat_param[cstrname_up]
            # if collect(clique_keys) != collect(clique_keys_up)
            #     warn(LOGGER, "clique keys different from lower and upper side of double constraint")
            #     length(clique_keys_up) == 1 || error(LOGGER, "build_momentrelaxation(): constraint $cstrname spans several cliques ($clique_keys).\nNot supported yet.")
            #     vars, cliquename = collect_cliquesvars(clique_keys_up, max_cliques)

            #     mmt = MomentMatrix(relax_ctx, vars, order_up, relax_ctx.symmetries,
            #                                                   relax_ctx.cstrtypes[cstrname_hi],
            #                                                   var_to_cliques = var_to_cliques)
            # elseif order_up != order
            #     warn(LOGGER, "order different from lower and upper side of double constraint")
            #     mmt = MomentMatrix(relax_ctx, vars, order_up, relax_ctx.symmetries,
            #                                                   relax_ctx.cstrtypes[cstrname_hi],
            #                                                   var_to_cliques = var_to_cliques)
            # end

        else
            # either cstrtype == :ineqlo, :ineqhi, :eq
            clique_keys, order = localizingmat_param[get_cstrname(cstrname, cstrtype)]
            vars, cliquename = collect_cliquesvars(clique_keys, max_cliques)
            # length(clique_keys) == 1 || error(LOGGER, "build_momentrelaxation(): constraint $cstrname spans several cliques ($clique_keys).\nNot supported yet.")

            mmt = MomentMatrix{T}(relax_ctx, vars, order, relax_ctx.symmetries,
                                                          relax_ctx.cstrtypes[get_cstrname(cstrname, cstrtype)],
                                                          var_to_cliques = var_to_cliques)

            # print_with_color(:green, "$cstrname, :Up\n") ##NOTE: find better logging system.
            product!(mmt, get_normalizedpoly(cstr, cstrtype), var_to_cliques)

            momentmatrices[(get_cstrname(cstrname, cstrtype), cliquename)] = mmt
        end
    end

    ## Locate clique overlapping moments
    expo_to_cliques = Dict{Exponent, Set{String}}()

    # Collect Exponents per clique (moment matrix)
    for ((ctrobj, clique), mmtmat) in momentmatrices
        if ctrobj == get_momentcstrname()

            for (key, momentpoly) in mmtmat.mm
                for (moment, coeff) in momentpoly
                    expo = product(moment.conj_part, moment.expl_part)
                    haskey(expo_to_cliques, expo) || (expo_to_cliques[expo] = Set{String}())
                    push!(expo_to_cliques[expo], moment.clique)
                end
            end
        end
    end

    nb_expos = length(expo_to_cliques)
    for (expo, cliques) in expo_to_cliques
        length(cliques) > 1 || delete!(expo_to_cliques, expo)
    end

    momentrelaxation = SDPDual{T}(objective, momentmatrices, expo_to_cliques)

    print_build_momentrelax(relax_ctx, momentrelaxation, nb_expos)

    return momentrelaxation
end


text:
using OMOPCommonDataModel
using Documenter

makedocs(;
    modules=[OMOPCommonDataModel],
    authors="Dilum Aluthge, Brown Center for Biomedical Informatics, JuliaHealth, and contributors",
    repo="https://github.com/JuliaHealth/OMOPCommonDataModel.jl/blob/{commit}{path}#L{line}",
    sitename="OMOPCommonDataModel.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaHealth.github.io/OMOPCommonDataModel.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "CDM Version" => "version.md",
        "Examples" => "examples.md",
        "API" => "api.md",
        "Auto-Generation" => "generate.md",
    ],
    strict=true,
)

deploydocs(;
    repo="github.com/JuliaHealth/OMOPCommonDataModel.jl",
)