text:
module Fermionfields
    using Random
    using LinearAlgebra
    using SparseArrays

    import ..Actions:FermiActionParam,FermiActionParam_Wilson,
                FermiActionParam_WilsonClover,FermiActionParam_Staggered
    import ..Gaugefields:GaugeFields,GaugeFields_1d,SU3GaugeFields,SU2GaugeFields,SU3GaugeFields_1d,SU2GaugeFields_1d,
                            staggered_phase,SUn,SU2,SU3,SUNGaugeFields,SUNGaugeFields_1d,SU
    #import ..Parallel:get_looprange

    #include("cg.jl")

    import ..AbstractFermion:FermionFields,
        Wx!,Wdagx!,clear!,substitute_fermion!,Dx!,Dxplus!,
        fermion_shift!,fermion_shiftB!,add!,set_wing_fermi!,WdagWx!,apply_periodicity,
        gauss_distribution_fermi!,gauss_distribution_fermi_Z2!,Z4_distribution_fermi!,Ddagx!

    import ..WilsonFermion_module:WilsonFermion,mul_γ5x!

    import ..StaggeredFermion_module:StaggeredFermion

    import ..DomainwallFermion_module:DomainwallFermion,D5DWx!,D5DWdagx!

#    abstract type FermionFields end

    struct RGammamatrix{ν,pm}
        #indices::NTuple{}
        #values::
    end

    function FermionFields(NC,NX,NY,NZ,NT,fparam::FermiActionParam,BoundaryCondition)
        if findfirst("Wilson",fparam.Dirac_operator) != nothing #If string fparam.Dirac_operator includes "Wilson"
            #fparam.Dirac_operator == "Wilson"
            return WilsonFermion(NC,NX,NY,NZ,NT,fparam,BoundaryCondition)
        elseif fparam.Dirac_operator == "Staggered"
            return StaggeredFermion(NC,NX,NY,NZ,NT,fparam,BoundaryCondition)
        elseif fparam.Dirac_operator == "Domainwall"
            return DomainwallFermion(NC,NX,NY,NZ,NT,fparam,BoundaryCondition)

        else
            error(fparam.Dirac_operator,"is not supported!")
        end
    end

    #include("boundary_fermi.jl")

    function Base.setindex!(x::FermionFields,v,i1,i2,i3,i4,i5,i6)
        x.f[i1,i2 + 1,i3 + 1,i4 + 1,i5 + 1,i6] = v
    end

    function Base.getindex(x::FermionFields,i1,i2,i3,i4,i5,i6)
        return x.f[i1,i2 .+ 1,i3 .+ 1,i4 .+ 1,i5 .+ 1,i6]
    end

    function LinearAlgebra.dot(a::FermionFields,b::FermionFields)
        return a*b
    end

    function Base.display(x::FermionFields)

        for it=1:x.NT
            for iz=1:x.NZ
                for iy=1:x.NY
                    for ix=1:x.NX
                        for ic=1:x.NC
                            println("$ic $ix $iy $iz $it")
                            display(x[ic,ix,iy,iz,it,:])
                            println("\t")
                        end
                    end
                end

            end
        end

    end

    function LinearAlgebra.axpby!(a::Number,X::FermionFields,Y::FermionFields)
        add!(Y,a,X)
        return
    end

    function LinearAlgebra.axpby!(a::Number,X::FermionFields,b,Y::FermionFields)
        add!(b,Y,a,X)
        return
    end

    function Wdagx_old!(xout::WilsonFermion,U::Array{G,1},
        x::WilsonFermion,temps::Array{T,1}) where {T <: FermionFields,G <: GaugeFields}
        temp = temps[4]
        temp1 = temps[1]
        temp2 = temps[2]

        clear!(temp)
        x5 = temps[3]

        #D^dag = \gamma_5 D \gamma_5

        mul_γ5x!(x5,x)
        set_wing_fermi!(x5)

        for ν=1:4
            fermion_shift!(temp1,U,ν,x5) #nu shift

            #... Dirac multiplication
            #mul!(temp1,view(x.rminusγ,:,:,ν),temp1) #(1 - \gamma)
            mul!(temp1,view(x.rplusγ,:,:,ν),temp1) #(1 + \gamma)

            #
            fermion_shift!(temp2,U,-ν,x5) #-nu shift

            #mul!(temp2,view(x.rplusγ,:,:,ν),temp2) #(1 + \gamma)
            mul!(temp2,view(x.rminusγ,:,:,ν),temp2) #(1 + \gamma)

            add!(temp,x.hopp[ν],temp1,x.hopm[ν],temp2)
        end
        clear!(temp1)
        add!(temp1,1,x5,-1,temp)

        mul_γ5x!(xout,temp1)
        return
    end

    #=
    function WdagWx!(xout::StaggeredFermion,U::Array{G,1},
        x::StaggeredFermion,temps::Array{T,1},fparam::FermiActionParam_Staggered) where {T <: FermionFields,G <: GaugeFields}
        temp = temps[5]
        evensite = true
        Dx!(temp,evensite,U,x,temps,fparam)
        #Dx!(temp,U,x,temps,fparam)
        temp3 = temps[3]
        evensite = false
        Dx!(temp3,evensite,U,temp,temps,fparam)
        #Dx!(temp3,U,temp,temps,fparam)

        clear!(xout)
        add!(xout,fparam.mass^2,x,-1,temp3)

        return
    end
    =#

    function LinearAlgebra.mul!(xout::FermionFields,A::AbstractMatrix,x::FermionFields)
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        NC = x.NC

        n6 = size(x.f)[6]
        f = zeros(ComplexF64,n6)
        e = zeros(ComplexF64,n6)

        for ic=1:NC
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            for k1=1:n6
                                e[k1] = x[ic,ix,iy,iz,it,k1]
                            end

                            for k1=1:n6
                                f[k1] = 0
                                @simd for k2=1:n6
                                    f[k1] += A[k1,k2]*e[k2]
                                end
                                xout[ic,ix,iy,iz,it,k1] = f[k1]
                            end

                        end
                    end
                end
            end
        end
    end

    function LinearAlgebra.mul!(xout::FermionFields,x::FermionFields,A::AbstractMatrix)
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        NC = x.NC

        n6 = size(x.f)[6]
        f = zeros(ComplexF64,n6)
        e = zeros(ComplexF64,n6)

        for ic=1:NC
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX

                            @simd for k1=1:n6
                                e[k1] = x[ic,ix,iy,iz,it,k1]
                            end

                            for k1=1:n6
                                f[k1] = 0
                                @simd for k2=1:n6
                                    f[k1] += e[k2]*A[k2,k1]
                                end
                                xout[ic,ix,iy,iz,it,k1] = f[k1]
                            end

                        end
                    end
                end
            end
        end
    end

    function LinearAlgebra.mul!(a::FermionFields,α::Number,b::FermionFields)
        n1,n2,n3,n4,n5,n6 = size(a.f)

        for i6=1:n6
            for i5=1:n5
                for i4=1:n4
                    for i3=1:n3
                        for i2=1:n2
                            @simd for i1=1:n1
                                a.f[i1,i2,i3,i4,i5,i6] = α*b.f[i1,i2,i3,i4,i5,i6]
                            end
                        end
                    end
                end
            end
        end
        return
    end
"""
    c-----------------------------------------------------c
c     iflag = 1
c       vv = |v1><v2| i.e., vv(a,b) = v1(a)*v2_adj(b)
c     iflag = 2
c                           vv(a,b) = v1(a)*v2(b)
c-----------------------------------------------------c
    """
    function vvmat!(vv::GaugeFields_1d{T},v1::FermionFields,v2::FermionFields,iflag) where T <: SUn
        NX = v1.NX
        NY = v1.NY
        NZ = v1.NZ
        NT = v1.NT

        if T == SU3
            NC = 3
        elseif T == SU2
            NC = 2
        else
            NC = v1.NC
            #error("NC >3 is not supported")
        end
        #else
        #    error("NC >3 is not supported")
        #end

        if iflag == 1

            is = 0
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            is += 1
                            for ia=1:NC
                                for ib=1:NC
                                    vv[ia,ib,is] = 0
                                    @simd for k=1:4
                                        vv[ia,ib,is] += v1[ia,ix,iy,iz,it,k]*conj(v2[ib,ix,iy,iz,it,k])
                                    end

                                end
                            end
                        end
                    end
                end
            end
        end
        if iflag == 2

                    is = 0

            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            is += 1
                            for ia=1:NC
                                for ib=1:NC
                                    vv[ia,ib,is] = 0
                                    @simd for k=1:4
                                        vv[ia,ib,is] += v1[ia,ix,iy,iz,it,k]*v2[ib,ix,iy,iz,it,k]
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    function vvmat!(vv::GaugeFields_1d{T},v1::StaggeredFermion,v2::StaggeredFermion,iflag) where T <: SUn
        NX = v1.NX
        NY = v1.NY
        NZ = v1.NZ
        NT = v1.NT

        if T == SU3
            NC = 3
        elseif T == SU2
            NC = 2
        else
            NC = v1.NC
            #error("NC >3 is not supported")
        end
        #else
        #    error("NC >3 is not supported")
        #end

        if iflag == 1

            is = 0
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            is += 1
                                for ia=1:NC
                                    for ib=1:NC
                                    vv[ia,ib,is] = 0
                                    @simd for k=1:1
                                        vv[ia,ib,is] += v1[ia,ix,iy,iz,it,k]*conj(v2[ib,ix,iy,iz,it,k])
                                    end

                                end
                            end
                        end
                    end
                end
            end
        end
        if iflag == 2

                    is = 0

            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            is += 1
                            for ia=1:NC
                                for ib=1:NC
                                    vv[ia,ib,is] = 0
                                    @simd for k=1:1
                                        vv[ia,ib,is] += v1[ia,ix,iy,iz,it,k]*v2[ib,ix,iy,iz,it,k]
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

end


text:
using DashBootstrapComponents

color_selector = html_div(
    [
        html_div("Select a colour theme:"),
        dbc_select(
            id = "change-table-color",
            options = [
                Dict("label" => "primary", "value" => "primary"),
                Dict("label" => "secondary", "value" => "secondary"),
                Dict("label" => "success", "value" => "success"),
                Dict("label" => "danger", "value" => "danger"),
                Dict("label" => "warning", "value" => "warning"),
                Dict("label" => "info", "value" => "info"),
                Dict("label" => "light", "value" => "light"),
                Dict("label" => "dark", "value" => "dark"),
            ],
            value = "primary",
        ),
    ],
    className = "p-3 m-2 border",
)

table = html_div([
    color_selector,
    dbc_table(
        # using the same table as in the above example
        table_header + table_body,
        id = "table-color",
        color = "primary",
    ),
])

callback!(
    app,
    Output("table-color", "color"),
    Input("change-table-color", "value"),
) do color
    return color
end;

ext: jl

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
# Julia wrapper for header: cimplot.h
# Automatically generated using Clang.jl

function ImPlotPoint_ImPlotPointNil()
    ccall((:ImPlotPoint_ImPlotPointNil, libcimplot), Ptr{ImPlotPoint}, ())
end

function ImPlotPoint_destroy(self)
    ccall((:ImPlotPoint_destroy, libcimplot), Cvoid, (Ptr{ImPlotPoint},), self)
end

function ImPlotPoint_ImPlotPointdouble(_x, _y)
    ccall((:ImPlotPoint_ImPlotPointdouble, libcimplot), Ptr{ImPlotPoint}, (Cdouble, Cdouble), _x, _y)
end

function ImPlotRange_ImPlotRange()
    ccall((:ImPlotRange_ImPlotRange, libcimplot), Ptr{ImPlotRange}, ())
end

function ImPlotRange_destroy(self)
    ccall((:ImPlotRange_destroy, libcimplot), Cvoid, (Ptr{ImPlotRange},), self)
end

function ImPlotRange_Contains(self, value)
    ccall((:ImPlotRange_Contains, libcimplot), Bool, (Ptr{ImPlotRange}, Cdouble), self, value)
end

function ImPlotRange_Size(self)
    ccall((:ImPlotRange_Size, libcimplot), Cdouble, (Ptr{ImPlotRange},), self)
end

function ImPlotLimits_ImPlotLimits()
    ccall((:ImPlotLimits_ImPlotLimits, libcimplot), Ptr{ImPlotLimits}, ())
end

function ImPlotLimits_destroy(self)
    ccall((:ImPlotLimits_destroy, libcimplot), Cvoid, (Ptr{ImPlotLimits},), self)
end

function ImPlotLimits_ContainsPlotPoInt(self, p)
    ccall((:ImPlotLimits_ContainsPlotPoInt, libcimplot), Bool,
          (Ptr{ImPlotLimits}, ImPlotPoint), self, p)
end

function ImPlotLimits_Containsdouble(self, x, y)
    ccall((:ImPlotLimits_Containsdouble, libcimplot), Bool,
          (Ptr{ImPlotLimits}, Cdouble, Cdouble), self, x, y)
end

function ImPlotStyle_ImPlotStyle()
    ccall((:ImPlotStyle_ImPlotStyle, libcimplot), Ptr{ImPlotStyle}, ())
end

function ImPlotStyle_destroy(self)
    ccall((:ImPlotStyle_destroy, libcimplot), Cvoid, (Ptr{ImPlotStyle},), self)
end

function BeginPlot(title_id, x_label, y_label, size,
                   flags, x_flags, y_flags, y2_flags, y3_flags)::Bool
    ccall((:ImPlot_BeginPlot, libcimplot), Bool,
          (Cstring, Cstring, Cstring, ImVec2,
           ImPlotFlags, ImPlotAxisFlags, ImPlotAxisFlags, ImPlotAxisFlags, ImPlotAxisFlags),
          title_id, x_label, y_label, size, flags, x_flags, y_flags, y2_flags, y3_flags)
end

function EndPlot()
    ccall((:ImPlot_EndPlot, libcimplot), Cvoid, ())
end

function PlotLineFloatPtrInt(label_id, values, count, offset, stride)
    ccall((:ImPlot_PlotLineFloatPtrInt, libcimplot), Cvoid,
          (Cstring, Ptr{Cfloat}, Cint, Cint, Cint),
          label_id, values, count, offset, stride)
end

function PlotLinedoublePtrInt(label_id, values, count, offset, stride)
    ccall((:ImPlot_PlotLinedoublePtrInt, libcimplot), Cvoid,
          (Cstring, Ptr{Cdouble}, Cint, Cint, Cint),
          label_id, values, count, offset, stride)
end

function PlotLineFloatPtrFloatPtr(label_id, xs, ys, count, offset, stride)
    ccall((:ImPlot_PlotLineFloatPtrFloatPtr, libcimplot), Cvoid,
          (Cstring, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint, Cint),
          label_id, xs, ys, count, offset, stride)
end

function PlotLinedoublePtrdoublePtr(label_id, xs, ys, count, offset, stride)
    ccall((:ImPlot_PlotLinedoublePtrdoublePtr, libcimplot), Cvoid,
          (Cstring, Ref{Cdouble}, Ref{Cdouble}, Cint, Cint, Cint),
          label_id, xs, ys, count, offset, stride)
end

function PlotLineVec2Ptr(label_id, data, count, offset)
    ccall((:ImPlot_PlotLineVec2Ptr, libcimplot), Cvoid,
          (Cstring, Ptr{Cint}, Cint, Cint),
          label_id, data, count, offset)
end

function PlotLinePlotPoIntPtr(label_id, data, count, offset)
    ccall((:ImPlot_PlotLinePlotPoIntPtr, libcimplot), Cvoid,
          (Cstring, Ptr{ImPlotPoint}, Cint, Cint),
          label_id, data, count, offset)
end

function PlotLineFnPlotPoIntPtr(label_id, getter, data, count, offset)
    ccall((:ImPlot_PlotLineFnPlotPoIntPtr, libcimplot), Cvoid,
          (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint),
          label_id, getter, data, count, offset)
end

function PlotScatterFloatPtrInt(label_id, values, count, offset, stride)
    ccall((:ImPlot_PlotScatterFloatPtrInt, libcimplot), Cvoid,
          (Cstring, Ptr{Cfloat}, Cint, Cint, Cint),
          label_id, values, count, offset, stride)
end

function PlotScatterdoublePtrInt(label_id, values, count, offset, stride)
    ccall((:ImPlot_PlotScatterdoublePtrInt, libcimplot), Cvoid,
          (Cstring, Ptr{Cdouble}, Cint, Cint, Cint),
          label_id, values, count, offset, stride)
end

function PlotScatterFloatPtrFloatPtr(label_id, xs, ys, count, offset, stride)
    ccall((:ImPlot_PlotScatterFloatPtrFloatPtr, libcimplot), Cvoid,
          (Cstring, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint, Cint),
          label_id, xs, ys, count, offset, stride)
end

function PlotScatterdoublePtrdoublePtr(label_id, xs, ys, count, offset, stride)
    ccall((:ImPlot_PlotScatterdoublePtrdoublePtr, libcimplot), Cvoid,
          (Cstring, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Cint),
          label_id, xs, ys, count, offset, stride)
end

function PlotScatterVec2Ptr(label_id, data, count, offset)
    ccall((:ImPlot_PlotScatterVec2Ptr, libcimplot), Cvoid,
          (Cstring, Ptr{Cint}, Cint, Cint),
          label_id, data, count, offset)
end

function PlotScatterPlotPoIntPtr(label_id, data, count, offset)
    ccall((:ImPlot_PlotScatterPlotPoIntPtr, libcimplot), Cvoid,
          (Cstring, Ptr{ImPlotPoint}, Cint, Cint), label_id, data, count, offset)
end

function PlotScatterFnPlotPoIntPtr(label_id, getter, data, count, offset)
    ccall((:ImPlot_PlotScatterFnPlotPoIntPtr, libcimplot), Cvoid, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint), label_id, getter, data, count, offset)
end

function PlotBarsFloatPtrIntFloat(label_id, values, count, width, shift, offset, stride)
    ccall((:ImPlot_PlotBarsFloatPtrIntFloat, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Cint, Cfloat, Cfloat, Cint, Cint), label_id, values, count, width, shift, offset, stride)
end

function PlotBarsdoublePtrIntdouble(label_id, values, count, width, shift, offset, stride)
    ccall((:ImPlot_PlotBarsdoublePtrIntdouble, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Cint, Cdouble, Cdouble, Cint, Cint), label_id, values, count, width, shift, offset, stride)
end

function PlotBarsFloatPtrFloatPtr(label_id, xs, ys, count, width, offset, stride)
    ccall((:ImPlot_PlotBarsFloatPtrFloatPtr, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cfloat, Cint, Cint), label_id, xs, ys, count, width, offset, stride)
end

function PlotBarsdoublePtrdoublePtr(label_id, xs, ys, count, width, offset, stride)
    ccall((:ImPlot_PlotBarsdoublePtrdoublePtr, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cdouble, Cint, Cint), label_id, xs, ys, count, width, offset, stride)
end

function PlotBarsFnPlotPoIntPtr(label_id, getter, data, count, width, offset)
    ccall((:ImPlot_PlotBarsFnPlotPoIntPtr, libcimplot), Cvoid, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cdouble, Cint), label_id, getter, data, count, width, offset)
end

function PlotBarsHFloatPtrIntFloat(label_id, values, count, height, shift, offset, stride)
    ccall((:ImPlot_PlotBarsHFloatPtrIntFloat, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Cint, Cfloat, Cfloat, Cint, Cint), label_id, values, count, height, shift, offset, stride)
end

function PlotBarsHdoublePtrIntdouble(label_id, values, count, height, shift, offset, stride)
    ccall((:ImPlot_PlotBarsHdoublePtrIntdouble, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Cint, Cdouble, Cdouble, Cint, Cint), label_id, values, count, height, shift, offset, stride)
end

function PlotBarsHFloatPtrFloatPtr(label_id, xs, ys, count, height, offset, stride)
    ccall((:ImPlot_PlotBarsHFloatPtrFloatPtr, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cfloat, Cint, Cint), label_id, xs, ys, count, height, offset, stride)
end

function PlotBarsHdoublePtrdoublePtr(label_id, xs, ys, count, height, offset, stride)
    ccall((:ImPlot_PlotBarsHdoublePtrdoublePtr, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cdouble, Cint, Cint), label_id, xs, ys, count, height, offset, stride)
end

function PlotBarsHFnPlotPoIntPtr(label_id, getter, data, count, height, offset)
    ccall((:ImPlot_PlotBarsHFnPlotPoIntPtr, libcimplot), Cvoid, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cdouble, Cint), label_id, getter, data, count, height, offset)
end

function PlotErrorBarsFloatPtrFloatPtrFloatPtrInt(label_id, xs, ys, err, count, offset, stride)
    ccall((:ImPlot_PlotErrorBarsFloatPtrFloatPtrFloatPtrInt, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint, Cint), label_id, xs, ys, err, count, offset, stride)
end

function PlotErrorBarsdoublePtrdoublePtrdoublePtrInt(label_id, xs, ys, err, count, offset, stride)
    ccall((:ImPlot_PlotErrorBarsdoublePtrdoublePtrdoublePtrInt, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Cint), label_id, xs, ys, err, count, offset, stride)
end

function PlotErrorBarsFloatPtrFloatPtrFloatPtrFloatPtr(label_id, xs, ys, neg, pos, count, offset, stride)
    ccall((:ImPlot_PlotErrorBarsFloatPtrFloatPtrFloatPtrFloatPtr, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint, Cint), label_id, xs, ys, neg, pos, count, offset, stride)
end

function PlotErrorBarsdoublePtrdoublePtrdoublePtrdoublePtr(label_id, xs, ys, neg, pos, count, offset, stride)
    ccall((:ImPlot_PlotErrorBarsdoublePtrdoublePtrdoublePtrdoublePtr, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Cint), label_id, xs, ys, neg, pos, count, offset, stride)
end

function PlotErrorBarsHFloatPtrFloatPtrFloatPtrInt(label_id, xs, ys, err, count, offset, stride)
    ccall((:ImPlot_PlotErrorBarsHFloatPtrFloatPtrFloatPtrInt, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint, Cint), label_id, xs, ys, err, count, offset, stride)
end

function PlotErrorBarsHdoublePtrdoublePtrdoublePtrInt(label_id, xs, ys, err, count, offset, stride)
    ccall((:ImPlot_PlotErrorBarsHdoublePtrdoublePtrdoublePtrInt, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Cint), label_id, xs, ys, err, count, offset, stride)
end

function PlotErrorBarsHFloatPtrFloatPtrFloatPtrFloatPtr(label_id, xs, ys, neg, pos, count, offset, stride)
    ccall((:ImPlot_PlotErrorBarsHFloatPtrFloatPtrFloatPtrFloatPtr, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint, Cint), label_id, xs, ys, neg, pos, count, offset, stride)
end

function PlotErrorBarsHdoublePtrdoublePtrdoublePtrdoublePtr(label_id, xs, ys, neg, pos, count, offset, stride)
    ccall((:ImPlot_PlotErrorBarsHdoublePtrdoublePtrdoublePtrdoublePtr, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Cint), label_id, xs, ys, neg, pos, count, offset, stride)
end

function PlotPieChartFloatPtr(label_ids, values, count, x, y, radius, normalize, label_fmt, angle0)
    ccall((:ImPlot_PlotPieChartFloatPtr, libcimplot), Cvoid, (Ptr{Cstring}, Ptr{Cfloat}, Cint, Cfloat, Cfloat, Cfloat, Bool, Cstring, Cfloat), label_ids, values, count, x, y, radius, normalize, label_fmt, angle0)
end

function PlotPieChartdoublePtr(label_ids, values, count, x, y, radius, normalize, label_fmt, angle0)
    ccall((:ImPlot_PlotPieChartdoublePtr, libcimplot), Cvoid, (Ptr{Cstring}, Ptr{Cdouble}, Cint, Cdouble, Cdouble, Cdouble, Bool, Cstring, Cdouble), label_ids, values, count, x, y, radius, normalize, label_fmt, angle0)
end

function PlotHeatmapFloatPtr(label_id, values, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
    ccall((:ImPlot_PlotHeatmapFloatPtr, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Cint, Cint, Cfloat, Cfloat, Cstring, ImPlotPoint, ImPlotPoint), label_id, values, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end

function PlotHeatmapdoublePtr(label_id, values, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
    ccall((:ImPlot_PlotHeatmapdoublePtr, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Cint, Cint, Cdouble, Cdouble, Cstring, ImPlotPoint, ImPlotPoint), label_id, values, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end

function PlotDigitalFloatPtr(label_id, xs, ys, count, offset, stride)
    ccall((:ImPlot_PlotDigitalFloatPtr, libcimplot), Cvoid, (Cstring, Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint, Cint), label_id, xs, ys, count, offset, stride)
end

function PlotDigitaldoublePtr(label_id, xs, ys, count, offset, stride)
    ccall((:ImPlot_PlotDigitaldoublePtr, libcimplot), Cvoid, (Cstring, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Cint), label_id, xs, ys, count, offset, stride)
end

function PlotDigitalFnPlotPoIntPtr(label_id, getter, data, count, offset)
    ccall((:ImPlot_PlotDigitalFnPlotPoIntPtr, libcimplot), Cvoid, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint), label_id, getter, data, count, offset)
end

function PlotTextFloat(text, x, y, vertical, pixel_offset)
    ccall((:ImPlot_PlotTextFloat, libcimplot), Cvoid, (Cstring, Cfloat, Cfloat, Bool, Cint), text, x, y, vertical, pixel_offset)
end

function PlotTextdouble(text, x, y, vertical, pixel_offset)
    ccall((:ImPlot_PlotTextdouble, libcimplot), Cvoid, (Cstring, Cdouble, Cdouble, Bool, Cint), text, x, y, vertical, pixel_offset)
end

function IsPlotHovered()
    ccall((:ImPlot_IsPlotHovered, libcimplot), Bool, ())
end

function GetPlotMousePos(y_axis)
    ccall((:ImPlot_GetPlotMousePos, libcimplot), ImPlotPoint, (Cint,), y_axis)
end

function GetPlotLimits(y_axis)
    ccall((:ImPlot_GetPlotLimits, libcimplot), ImPlotLimits, (Cint,), y_axis)
end

function IsPlotQueried()
    ccall((:ImPlot_IsPlotQueried, libcimplot), Bool, ())
end

function GetPlotQuery(y_axis)
    ccall((:ImPlot_GetPlotQuery, libcimplot), ImPlotLimits, (Cint,), y_axis)
end

function GetStyle()
    ccall((:ImPlot_GetStyle, libcimplot), Ptr{ImPlotStyle}, ())
end

function PushStyleColorU32(idx, col)
    ccall((:ImPlot_PushStyleColorU32, libcimplot), Cvoid, (ImPlotCol, Cint), idx, col)
end

function PushStyleColorVec4(idx, col)
    ccall((:ImPlot_PushStyleColorVec4, libcimplot), Cvoid, (ImPlotCol, Cint), idx, col)
end

function PopStyleColor(count)
    ccall((:ImPlot_PopStyleColor, libcimplot), Cvoid, (Cint,), count)
end

function PushStyleVarFloat(idx, val)
    ccall((:ImPlot_PushStyleVarFloat, libcimplot), Cvoid, (ImPlotStyleVar, Cfloat), idx, val)
end

function PushStyleVarInt(idx, val)
    ccall((:ImPlot_PushStyleVarInt, libcimplot), Cvoid, (ImPlotStyleVar, Cint), idx, val)
end

function PopStyleVar(count)
    ccall((:ImPlot_PopStyleVar, libcimplot), Cvoid, (Cint,), count)
end

function SetColormapPlotColormap(colormap, samples)
    ccall((:ImPlot_SetColormapPlotColormap, libcimplot), Cvoid, (ImPlotColormap, Cint), colormap, samples)
end

function SetColormapVec4Ptr(colors, num_colors)
    ccall((:ImPlot_SetColormapVec4Ptr, libcimplot), Cvoid, (Ptr{Cint}, Cint), colors, num_colors)
end

function GetColormapSize()
    ccall((:ImPlot_GetColormapSize, libcimplot), Cint, ())
end

function GetColormapColor(pOut, index)
    ccall((:ImPlot_GetColormapColor, libcimplot), Cvoid, (Ptr{Cint}, Cint), pOut, index)
end

function LerpColormap(pOut, t)
    ccall((:ImPlot_LerpColormap, libcimplot), Cvoid, (Ptr{Cint}, Cfloat), pOut, t)
end

function SetNextPlotLimits(x_min, x_max, y_min, y_max, cond)
    ccall((:ImPlot_SetNextPlotLimits, libcimplot), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cint), x_min, x_max, y_min, y_max, cond)
end

function SetNextPlotLimitsX(x_min, x_max, cond)
    ccall((:ImPlot_SetNextPlotLimitsX, libcimplot), Cvoid, (Cdouble, Cdouble, Cint), x_min, x_max, cond)
end

function SetNextPlotLimitsY(y_min, y_max, cond, y_axis)
    ccall((:ImPlot_SetNextPlotLimitsY, libcimplot), Cvoid, (Cdouble, Cdouble, Cint, Cint), y_min, y_max, cond, y_axis)
end

function SetNextPlotTicksXdoublePtr(values, n_ticks, labels, show_default)
    ccall((:ImPlot_SetNextPlotTicksXdoublePtr, libcimplot), Cvoid, (Ptr{Cdouble}, Cint, Ptr{Cstring}, Bool), values, n_ticks, labels, show_default)
end

function SetNextPlotTicksXdouble(x_min, x_max, n_ticks, labels, show_default)
    ccall((:ImPlot_SetNextPlotTicksXdouble, libcimplot), Cvoid, (Cdouble, Cdouble, Cint, Ptr{Cstring}, Bool), x_min, x_max, n_ticks, labels, show_default)
end

function SetNextPlotTicksYdoublePtr(values, n_ticks, labels, show_default, y_axis)
    ccall((:ImPlot_SetNextPlotTicksYdoublePtr, libcimplot), Cvoid, (Ptr{Cdouble}, Cint, Ptr{Cstring}, Bool, Cint), values, n_ticks, labels, show_default, y_axis)
end

function SetNextPlotTicksYdouble(y_min, y_max, n_ticks, labels, show_default, y_axis)
    ccall((:ImPlot_SetNextPlotTicksYdouble, libcimplot), Cvoid, (Cdouble, Cdouble, Cint, Ptr{Cstring}, Bool, Cint), y_min, y_max, n_ticks, labels, show_default, y_axis)
end

function SetPlotYAxis(y_axis)
    ccall((:ImPlot_SetPlotYAxis, libcimplot), Cvoid, (Cint,), y_axis)
end

function GetPlotPos(pOut)
    ccall((:ImPlot_GetPlotPos, libcimplot), Cvoid, (Ptr{Cint},), pOut)
end

function GetPlotSize(pOut)
    ccall((:ImPlot_GetPlotSize, libcimplot), Cvoid, (Ptr{Cint},), pOut)
end

function PixelsToPlot(pix, y_axis)
    ccall((:ImPlot_PixelsToPlot, libcimplot), ImPlotPoint, (Cint, Cint), pix, y_axis)
end

function PlotToPixels(pOut, plt, y_axis)
    ccall((:ImPlot_PlotToPixels, libcimplot), Cvoid, (Ptr{Cint}, ImPlotPoint, Cint), pOut, plt, y_axis)
end

function ShowColormapScale(scale_min, scale_max, height)
    ccall((:ImPlot_ShowColormapScale, libcimplot), Cvoid, (Cdouble, Cdouble, Cfloat), scale_min, scale_max, height)
end

function PushPlotClipRect()
    ccall((:ImPlot_PushPlotClipRect, libcimplot), Cvoid, ())
end

function PopPlotClipRect()
    ccall((:ImPlot_PopPlotClipRect, libcimplot), Cvoid, ())
end

function ShowDemoWindow(p_open)
    ccall((:ImPlot_ShowDemoWindow, libcimplot), Cvoid, (Ptr{Bool},), p_open)
end


text:
module BirdNPC

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Bird NPC" => Ahorn.EntityPlacement(
        Maple.Bird
    )
)

# Values might be wrong
# -1 = left, 1 = right
const modeFacingScale = Dict{String, Integer}(
    "climbingtutorial" => -1,
    "dashingtutorial" => 1,
    "dreamjumptutorial" => 1,
    "superwalljumptutorial" => -1,
    "hyperjumptutorial" => -1,
    "movetonodes" => -1,
    "waitforlightningoff" => -1,
    "flyaway" => -1,
    "sleeping" => 1,
    "none" => -1
)

Ahorn.nodeLimits(entity::Maple.Bird) = 0, -1

Ahorn.editingOptions(entity::Maple.Bird) = return Dict{String, Any}(
    "mode" => Maple.bird_npc_modes
)

sprite = "characters/bird/crow00"

function Ahorn.selection(entity::Maple.Bird)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)
    key = lowercase(get(entity.data, "mode", "Sleeping"))
    scaleX = get(modeFacingScale, key, -1)

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y, sx=scaleX, jx=0.5, jy=1.0)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.getSpriteRectangle(sprite, nx, ny, sx=scaleX, jx=0.5, jy=1.0))
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Bird)
    px, py = Ahorn.position(entity)
    key = lowercase(get(entity.data, "mode", "Sleeping"))
    scaleX = get(modeFacingScale, key, -1)

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)

        theta = atan(py - ny, px - nx)
        Ahorn.drawArrow(ctx, px, py - 8, nx + cos(theta) * 8, ny + sin(theta) * 8 - 8, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny, sx=scaleX, jx=0.5, jy=1.0)

        px, py = nx, ny
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Bird, room::Maple.Room)
    key = lowercase(get(entity.data, "mode", "Sleeping"))
    scaleX = get(modeFacingScale, key, -1)

    Ahorn.drawSprite(ctx, sprite, 0, 0, sx=scaleX, jx=0.5, jy=1.0)
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
using ForwardKinematics
using Test

@testset "Simple FK" begin
  point = forward_kinematics([pi], [1.0])
  @test point.x ≈ -1 atol=0.00001
  @test point.y ≈ 0 atol=0.00001
end

@testset "Exception due to unequal argument sizes" begin
  try
    forward_kinematics([3.1415], [1.0, 2.0])
    @test false
  catch err
  end
end

@testset "Complex FK" begin
  point = forward_kinematics([0.4, 1.0, 2.0, 0.8, 0.0], [1.0, 2.0, 0.5, 2.0, 1.0])
  @test point.x ≈ -0.693186 atol=0.00001
  @test point.y ≈ -0.382180 atol=0.00001
end


text:
using ImageMagick
using Images, TestImages, StaticArrays, ImageTracking, ImageView, LinearAlgebra, CoordinateTransformations, Gtk.ShortNames
#=Image Credit:  C. Liu. Beyond Pixels: Exploring New Representations and
#Applications for Motion Analysis. Doctoral Thesis. Massachusetts Institute of
#Technology. May 2009. =#
img1 = load("./error_evaluation_demo_data/frame1.png")
img2 = load("./error_evaluation_demo_data/frame2.png")

algorithm = Farneback(50, estimation_window = 11,
                         σ_estimation_window = 9.0,
                         expansion_window = 6,
                         σ_expansion_window = 5.0)

estimated_flow = optical_flow(Gray{Float32}.(img1), Gray{Float32}.(img2), algorithm)

# read the ground truth flow
file_path =  "./error_evaluation_demo_data/example.flo"
ground_truth_flow = read_flow_file(file_path)

# Endpoint error statistics
endpoint_error = evaluate_flow_error(ground_truth_flow, estimated_flow, EndpointError())
ep_mean, ep_sd, ep_RX, ep_AX = calculate_statistics(endpoint_error)
println("Endpoint error statistics :")
println("Mean = ", ep_mean, ", standard deviation = ", ep_sd, ", RX Dict = ", ep_RX, ", AX Dict = ", ep_AX)

# Angular error statistics
angular_error = evaluate_flow_error(ground_truth_flow, estimated_flow, AngularError())
ang_mean, ang_sd, ang_RX, ang_AX = calculate_statistics(angular_error)
println("Angular error statistics :")
println("Mean = ", ang_mean, ", standard deviation = ", ang_sd, ", RX Dict = ", ang_RX, ", AX Dict = ", ang_AX)


text:
module StocasExamples

using Stocas

# An example from Aguirregabiria Mira 2010
function am2010(;Nmarket = 5, β = 0.95)
    ### Aguirregabiria & Mira
    ## State 1: Endogenous incumbancy state
    X_entry = EntryState(;dense=false)

    ## State 2: Exogenous market state
    # Construct a market state with typical probability matrix as below
    # (shown for Nmarket = 5)
    # M = [0.8 0.2 0.0 0.0 0.0;
    #      0.2 0.6 0.2 0.0 0.0;
    #      0.0 0.2 0.6 0.2 0.0;
    #      0.0 0.0 0.2 0.6 0.2;
    #      0.0 0.0 0.0 0.2 0.8;]

    X2 = 1:Nmarket
    M = zeros(Nmarket, Nmarket)
    p = [0.2, 0.6, 0.2]
    # First corner
    M[1, 1:2] = [sum(p[1:2]), p[3]]
    # Last corner
    M[end, end-1:end] = [p[1], sum(p[2:3])]
    # Second to second last rows
    M[2:end-1, :] = vcat([vcat(zeros(i-2), p, zeros(Nmarket-1-i))' for i = 2:Nmarket-1]...)

    # Combine states
    S = States(X_entry, CommonState(X2, M))

    ## Utility specification
    # Z matrices for
    # exit choice:
    Z1 = zeros(Nmarket*2, 3)
    # entry choice:
    Z2 = [-ones(Nmarket) log.(X2) -ones(Nmarket);
          -ones(Nmarket) log.(X2)  zeros(Nmarket)]
    Z = (Z1, Z2)

    truepar = [-1.9;1.;2.] # fixed costs, ..., entry costs
    U = LinearUtility(Z, β, truepar)

    return U, S
end

# As AM but with a Tauchen discretized Gaussian market state
function am_tauchen(;β=0.95, # discount factor,
                     N = 500)
    # Aguirregabiria & Magesan 2016
    truepar = [-1.9;1.;2.] # fixed costs, ..., entry costs

    nX2 = N
    X2 = 1:nX2
    F2 = Stocas.tauchen(0.1,0.01,2, nX2)

    S = States(EntryState(),
               CommonState(X2, F2))

    Z = (zeros(nX2*2, 3), # exit
        [-ones(nX2) log.(X2) -ones(nX2); # entry
         -ones(nX2) log.(X2)  zeros(nX2)])

    U = LinearUtility(Z, β, truepar)

    return U, S
end

# A simple market entry model
function dixit(;N = 5, β = 0.99)
    # State 1
    X1 = 0:1
    F1 = [[1. 0.; 1. 0.], [0. 1.; 0. 1.]]
    # State 2
    nX2 = N
    X2 = 1:nX2
    F2 = 1./(1+abs.(ones(length(X2),1)*X2'-X2*ones(1, length(X2))))
    F2 = F2./sum(F2,1)'
    # States
    S = States(State("incumbancy", X1, F1), CommonState("market", X2, F2))

    # Utility variables
    Z1 = zeros(nX2*2, 3)
    Z2 = [ones(nX2) X2 -ones(nX2);
          ones(nX2) X2 zeros(nX2)]
    U = LinearUtility((Z1, Z2), β, [-.1;.2;1])

    return U, S
end

# An example from a presentation by Günter Hitsch
function hitsch(;N = 20)
    n = [0, 2, 5]
    next_i(i, k) = max.(0, min.(I, i+n[k]-1))
    F_X = ["low", "high"]
    F_P = [0.16 0.84; 0.16 0.84]
    SP = CommonState(F_X, F_P)

    I = N
    K = length(n)
    F_i = [spzeros(1+I, 1+I) for k = 1:K]
    for i = 0:I, j = 1:K
            F_i[j][i+1, next_i(i, j)+1]=1
    end
    Si = State(0:I, F_i);

    S = States(SP, Si)

    delta, alpha = 4., 4.

    Z1 = [[zeros(1); ones(I)] zeros(I+1) -next_i(0:I, 1);
          [zeros(1); ones(I)] zeros(I+1) -next_i(0:I, 1)]
    Z2 = [ones(2*(I+1)) -[1.2*ones(I+1);2*ones(I+1)] -kron(ones(2), next_i(0:I, 2))]
    Z3 = [ones(2*(I+1)) -[3.0*ones(I+1);5*ones(I+1)] -kron(ones(2), next_i(0:I, 3))]
    U = LinearUtility(("Buy 0", "Buy 2", "Buy 5"), (Z1, Z2, Z3), 0.998, [delta; alpha; 0.05])
    return U, S
end

# Rust 1987
function rust(;N = 175, β = 0.95, sparse = false)
    # State space
    nX = N
    n = nX
    o_maximum = 450
    X = linspace(0, o_maximum, nX)
    p = [0.0937; 0.4475; 0.4459; 0.0127; 0.0002];

    F1 = zeros(nX, nX)
    offset = -1
    for i = 1:nX, j = 1:length(p)
        i+offset+j > nX && continue
        F1[i,i+offset+j] = i+offset+j == nX ? sum(p[j:end]) : p[j]
    end

    # We can handle sparse matrices
    F = [F1, F1[ones(Int64, nX), :]]
    sparse == true && map!(sparse, F)
    S = State(:mileage, X, F)

    Z1 = [zeros(nX) -0.001*X]
    Z2 = [-ones(nX) zeros(nX)]

    U = LinearUtility(("replace", "repair"),(Z1, Z2), β, [11.;2.5])

    return U, S
end
# A somewhat contrived day effect model
function seven(;N = 5, β = 0.95)
	# Aguirregabiria & Mira
	# 2 Model
	truepar = [-1.9;1.;2.] # fixed costs, ..., entry costs
	# State 1: Exogenous state
	nX1 = 2
	X1 = [0.;1.]
	F1 = [sparse([1. 0.; 1. 0.]), sparse([0. 1.; 0. 1.])]

	nX2 = N
	X2 = 1:nX2
	M = [zeros(nX2-1) eye(nX2-1, nX2-1); [1.0 zeros(1,nX2-1)]]

	S = States(State(X1, F1),
	           CommonState(X2, M))

	Z = (zeros(nX2*2, 3), # exit
	             [-ones(nX2) log.(X2) -ones(nX2); # entry
				  -ones(nX2) log.(X2)  zeros(nX2)])

	U = LinearUtility(Z, β, truepar)

    return U, S
end
end # module


text:
module FANN

# from sundials code
if isfile(joinpath(dirname(dirname(@__FILE__)),"deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("Fann not properly installed. Please run Pkg.build(\"FANN\")")
end

# C FILE pointer
immutable FILE
    ptr::Ptr{Void}
end

# load library wrapper
include("libfann_common.jl")

# load files
include("data.jl")
include("ann.jl")
include("constants.jl")

# exports
export DataSet, savedset, loaddset
export ANN, train!, predict, mse, savenet, loadnet, weights, weights!
export setup_qprop!, setup_rprop!, setup_bprop!

end


text:
using SphericalLearning
using Gadfly
##

##
spherical_plot((x->harmonic(hologram(4,[0,0,.9]),x)),50)
##
spherical_plot((x->helmholtzG([0,0,.9],x;k=5)),50)

##

x = randn(3)
y = randn(3)

## these should show that the functions
Gadfly.plot([
    (z->imag(helmholtzG(z*x,y,k=1+0.2im))  ),
    (z->imag(helmholtzGapprox(z*x,y,k=1+0.2im))),
    (z->real(helmholtzG(z*x,y,k=1+0.2im))  ),
    (z->real(helmholtzGapprox(z*x,y,k=1+0.2im)) )
    ],0.2,7)
##
Gadfly.plot([
    (x->imag(sum(harmonic(hologram(l,[0,0,.9]),[0,1,x]) for l in 0:12))),
    (x->imag(helmholtzG([0,0,.9],normalize([0,1,x])))+0.01),
    (x->real(sum(harmonic(hologram(l,[0,0,.9]),[0,1,x]) for l in 0:12))),
    (x->real(helmholtzG([0,0,.9],normalize([0,1,x]))))
    ], -1,1)


text:
module MyNewPackage

export testPackage

# Write your package code here.
"""
    println("Hello World")

    Comment of function goes here

        # Parameter
        - None

        # Return
        - None

"""
function testPackage()
    println("Hello World")

end

end


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ 5eebac68-891b-11eb-37ba-c1d9481c6134
begin
  using Pkg; Pkg.activate("MLJ_env", shared=true)
	using Test
	using Random
	using PlutoUI
	using Printf
end

# ╔═╡ 87e5e2ec-8915-11eb-362b-6bf13a36b8e4
md"""
## Naive Bayes

ref. from book **"Data Science from Scratch"**, Chap 13

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ cb25501a-8915-11eb-3626-a1ae6d233f92
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ cb0838ae-8915-11eb-082f-3991192a101f
html"""

"""

# ╔═╡ caf03ada-8915-11eb-2935-01e110accba6
html"""

"""

# ╔═╡ cad66356-8915-11eb-0173-ada200678f63
md"""
### The model

First we will create a function to tokenize messages into tokens.
"""

# ╔═╡ dc3b7280-8915-11eb-258a-1789476897b9
function tokenize(text::String)::Set{String}
	lowercase(text) |>
		split |>
	 	a -> filter(s -> occursin(r"\A[\w']+.?\z", s), a) |>
		a -> map(s -> replace(s, r"\A([\w']+)[^w]?\z" => s"\1"), a) |>
		Set
end

# ╔═╡ dc209b9a-8915-11eb-0f63-ab3b7c33793c
begin
	test_msg₀ = "Machine Leaning isn't that fun, if all one does is Data Munging!"
	test_set₀ = tokenize(test_msg₀)

	@test typeof(test_set₀) == Set{String}
	@test length(test_set₀) == 12
	@test test_set₀ == Set{String}(["is", "data", "munging", "one",
		"that", "leaning", "if", "machine", "all", "fun", "isn't", "does"
	])
end

# ╔═╡ 5a207b90-89d3-11eb-2a3e-655804b5bc8b
md"""
Second, let us define a structure for our messages:
"""

# ╔═╡ 6a8e15fc-89d3-11eb-3c28-43b3d0ad0e23
struct Message
	text::String
	is_spam::Bool

	Message(text::String; is_spam::Bool=false) = new(text, is_spam)
end

# ╔═╡ d7aade44-891b-11eb-00d3-abe7406c4b09
md"""
As our classifier needs to keep track of tokens, counts, and labels from the training
data, we will create a (mutable) struct for this and a collection of related functions.
"""

# ╔═╡ 6dc9648a-89d6-11eb-35d0-e30bf7e3d3fc
const TF = Float64

# ╔═╡ 9ae99efe-891c-11eb-3d1f-89182e3f9ff8
mutable struct NaiveBayes
	k::TF
	tokens::Set{String}
	token_spam_cnt::Dict{String, Integer}
	token_ham_cnt::Dict{String, Integer}
	n_spam_msg::Integer
	n_ham_msg::Integer

	function NaiveBayes(;k::TF=0.5)
		@assert k > zero(TF)
		new(k, Set{String}(), Dict{String, Integer}(), Dict{String, Integer}(),
			0, 0)
	end
end

# ╔═╡ d7931368-891b-11eb-2899-f5e5a3f4bd35
md"""
Next, let us define a function to train it on a collection of messages.
"""

# ╔═╡ ea1f4a2c-891d-11eb-3027-bd447dd9d1e1
function train(nb::NaiveBayes, messages::Vector{Message})
	for msg ∈ messages
		if msg.is_spam
			nb.n_spam_msg += 1
		else
			nb.n_ham_msg += 1
		end
		for token ∈ tokenize(msg.text)
			push!(nb.tokens, token)
			if msg.is_spam
				nb.token_spam_cnt[token] = get(nb.token_spam_cnt, token, 0) + 1
			else
				nb.token_ham_cnt[token] = get(nb.token_ham_cnt, token, 0) + 1
			end
		end
	end
end

# ╔═╡ 50fe5114-891e-11eb-18ee-35d98ad32b7b
md"""
We will want to predict $P(spam | token)$. In order to apply Bayes’s theorem we need to know $P(token | spam)$ and $P(token | ham)$ for each token in the vocabulary.
"""

# ╔═╡ 4ffc3b48-895a-11eb-1ed2-f78b522d0ae9
function probabilities(nb::NaiveBayes, token::String)::Tuple{TF, TF}
	spam = get(nb.token_spam_cnt, token, 0)
	ham = get(nb.token_ham_cnt, token, 0)
	#
	p_token_spam = (spam + nb.k) / (nb.n_spam_msg + 2. * nb.k)
	p_token_ham = (ham + nb.k) / (nb.n_ham_msg + 2. * nb.k)
	#
	(p_token_spam, p_token_ham)
end

# ╔═╡ 90f58d04-8958-11eb-170e-0f17904c9f2c
md"""
Finally our predict function, where to avoid numerical underflow, we will sum up the log of probabilites instead of multiplying probabilities, before converting them back using exponentiaL.
"""

# ╔═╡ 123bcdb0-89d6-11eb-18f2-7f1acdb57a51
function predict(nb::NaiveBayes, text::String)::TF
	tokens = tokenize(text)
	log_prob_spam = log_prob_ham = zero(TF)
	#
	for token ∈ nb.tokens
		prob_spam, prob_ham = probabilities(nb, token)
		if token ∈ tokens
			log_prob_spam += log(prob_spam)
			log_prob_ham += log(prob_ham)
		else
			# otherwise add the log probability of _not_ seeing token in message
			log_prob_spam += log(one(TF) - prob_spam)
			log_prob_ham += log(one(TF) - prob_ham)
		end
	end
	prob_spam, prob_ham = exp(log_prob_spam), exp(log_prob_ham)
	prob_spam / (prob_spam + prob_ham)
end

# ╔═╡ dc0c5a5e-8915-11eb-2c62-5328b42add69
html"""

"""

# ╔═╡ dbefe77a-8915-11eb-39f1-5d0c84729739
md"""
### Testing Our Model
"""

# ╔═╡ dbd86924-8915-11eb-351f-8362f09ba984
begin
	test_msgs = [Message("spam rules", is_spam=true),
		Message("ham rules", is_spam=false),
		Message("hello ham", is_spam=false)
               ]
	test_model = NaiveBayes(;k=0.5)
	train(test_model, test_msgs)
	#
	@test test_model.tokens == Set{String}(["spam", "ham", "rules", "hello"])
	@test test_model.n_spam_msg == 1
	@test test_model.n_ham_msg == 2
	@test test_model.token_spam_cnt == Dict("spam" => 1, "rules" => 1)
	@test test_model.token_ham_cnt == Dict("ham" => 2, "rules" => 1, "hello" => 1)
end

# ╔═╡ c7410266-89d6-11eb-36b2-71e5df086e52
md"""
Ok, now let’s make a prediction.
"""

# ╔═╡ d32b86e0-8928-11eb-193c-e3d7c85d23dd
begin
	test_text = "hello spam"
	d_s = (1. + 2. * .5)
	probs_spam = [
		(1. + .5) / d_s,       # 'spam' present
		1. - (0. + .5) / d_s,  # 'ham' not present
		1. - (1. + .5) / d_s,  # 'rules' not present
		(0. + .5) / d_s,       # 'hello' present
	]
	d_h = (2. + 2. * .5)
	probs_ham = [
		(0. + .5) / d_h,       # 'spam' present
		1. - (2. + .5) / d_h,  # 'ham' not present
		1. - (1. + .5) / d_h,  # 'rules' not present
		(1. + .5) / d_h,       # 'hello' present
	]
	p_spam = exp(sum(log.(probs_spam)))
	p_ham = exp(sum(log.(probs_ham)))
	#
	@test predict(test_model, test_text) ≈ p_spam / (p_spam + p_ham)
end

# ╔═╡ ab9ac56e-8969-11eb-3340-337957fd81b7
md"""
And next let us try on some real data.
"""

# ╔═╡ dbbf22a2-8915-11eb-00eb-4b0278c0283d
html"""

"""

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

# ╔═╡ 5f5e6d02-892a-11eb-32d8-dfc1a7f259fa
md"""
### Using our Model

A popular dataset is the SpamAssassin public corpus. We will look at the files prefixed with 20021010.

First let us download and untar the files
"""

# ╔═╡ 073f6442-892e-11eb-0fc2-bbeb1d12619e
begin
	const BASE_URL = "https://spamassassin.apache.org/old/publiccorpus"
	const FILE_SDIRS = [
		("20021010_easy_ham.tar.bz2", "easy_ham"),
		("20021010_hard_ham.tar.bz2", "hard_ham"),
		("20021010_spam.tar.bz2", "spam")
	]
	const OUTPUT_DIR = "./spam_data/"
end

# ╔═╡ 6a454512-89e8-11eb-303a-ad5920ec581d
begin
	prepdir() = !isdir(OUTPUT_DIR) && mkdir(OUTPUT_DIR)

	function download_files()
		for (file, _sdir) ∈ FILE_SDIRS
			fp = string(OUTPUT_DIR, "/", file)
			isfile(fp) && continue
			## Note that this function relies on the availability of external tools
			## such as curl, wget or fetch to download the file and is provided
			## for convenience.
			res = download(string(BASE_URL, "/", file), fp)
		end
	end

	function untar_files()
		"""
		Create 3 subdirs: spam, easy_ham, hard_ham unless already created...
		"""
		for (file, sdir) ∈ FILE_SDIRS
			isdir(sdir) && continue
			cmd = Cmd(`tar xjf $(file)`, ignorestatus=false, detach=false,
				dir=OUTPUT_DIR)
			run(cmd, wait=true)
		end
	end

	function doit()
		prepdir()
		download_files()
		untar_files()
	end
end

# ╔═╡ aa4dae40-89ed-11eb-2cb6-9734911ab226
doit()

# ╔═╡ 2be56d9e-89ee-11eb-0e39-bd2174d78266
md"""
Next we need to load the files and to keep things simple (for now) we will just load the subject of the email...
"""

# ╔═╡ 21668a04-8a02-11eb-07b7-0f83275b0ad6
function cleanup_entry(line::String)::String
	try
	line |>
		s -> replace(s, r"[^\w^\s^']+" => "") |>
		strip |>
		s -> filter(c -> isvalid(c), collect(s)) |>   ## pass invalid char
		join |>                                       ## re-create string
		split |>                                      ## drop word with less...
		a -> filter(s -> 2 < length(s) < 11, a) |>    ## than 2 letters and nore
		#                                             ## than 11
		a -> filter(s -> !occursin(r"\d+", s), a) |>  ## drop word with digit
		a -> join(a, " ") |>
		strip

	catch
		@warn "problem with line: <$(line)>"
		return ""
	end
end

# ╔═╡ aa3808b2-8a0e-11eb-32e8-e70ec61b55c7
function load_data()
	"""
	Select only the Subject line from the message
	"""
	data = Vector{Message}()
	for (root, dirs, files) ∈ walkdir(OUTPUT_DIR;
		topdown=true, follow_symlinks=false)

		for file in files
			fp = joinpath(root, file)
			is_spam = occursin(r"ham", fp)

			open(fp) do fh
				for line ∈ readlines(fh)
					!occursin(r"\ASubject:", line) && continue
					subject = replace(line, "Subject: " => "") |>
						cleanup_entry
					push!(data, Message(string(subject); is_spam))
					break  # we are done, next
				end
			end
		end
	end
	data
end

# ╔═╡ 2bcd8a62-89ee-11eb-0575-dfa083fc40d6
function load_full_data()
	"""
	Select Subject and Body from the message
	Body appears after Subject + 1 empty line
	"""
	data = Vector{Message}()
	for (root, dirs, files) ∈ walkdir(OUTPUT_DIR;
		topdown=true, follow_symlinks=false)

		for file in files
			fp = joinpath(root, file)
			is_spam = occursin(r"ham", fp)
			found_subject = false
			body, subject = "", ""

			open(fp) do fh
				for line ∈ readlines(fh)
					# 78line = strip(line)
					if occursin(r"\ASubject:", line)
						subject = replace(line, "Subject: " => "") |>
							cleanup_entry
						found_subject = true
					elseif length(line) == 0
						continue
					elseif found_subject
						# cumulate cleaned up line found in body
						body = string(body, " ", cleanup_entry(line))
					end
				end
			end
			# finally make up full message as subject + body
			msg = string(subject, " ", body)
			push!(data, Message(msg; is_spam))
		end
	end
	data
end

# ╔═╡ 2bb0e580-89ee-11eb-24ae-a3947d35f204
begin
	data = load_full_data()
	length(data), data
end

# ╔═╡ 2b949b4e-89ee-11eb-0efa-9bcde1769743
function train_test_split(ds::Vector{Message};
                split=0.8, seed=42, shuffled=true)
	Random.seed!(seed)
	nr = length(ds)
	row_ixes = shuffled ? shuffle(1:nr) : collect(1:nr)
	nrp = round(Int, length(row_ixes) * split)
	(ds[row_ixes[1:nrp]], ds[row_ixes[nrp+1:nr]])
end

# ╔═╡ 347e75f2-89f5-11eb-2a08-d9b2bdda1511
begin
	train_messages, test_messages = train_test_split(data; split=0.75)
	model = NaiveBayes()
	train(model, train_messages)
	#
	length(train_messages), length(test_messages)
end

# ╔═╡ 5bcb82ee-89fa-11eb-36d7-c5668ac807d5
ŷ = [(msg, predict(model, msg.text)) for msg in test_messages]

# ╔═╡ 82175c86-89fa-11eb-1703-b5665587eedb
begin
	function confusion_matrix(ŷ)
		# Assume that spam_probability > 0.5 corresponds to spam prediction
		# and count the combinations of (actual is_spam, predicted is_spam)
		conf_matrix = Dict{Tuple{Bool, Bool}, Integer}()
		for (msg, spam_prob) ∈ ŷ
			#          real label   pred
			keypair = (msg.is_spam, spam_prob > 0.5)
			conf_matrix[keypair] = get(conf_matrix, keypair, 0) + 1
		end
		conf_matrix
	end

	tp(cm) = cm[(true, true)]     # True Positive
	tn(cm) = cm[(false, false)]   # True Negative
	fp(cm) = cm[(false, true)]    # False Positive
	fn(cm) = cm[(true, false)]    # False Negative

	function precision(cm)
		tp_, fp_ = tp(cm), fp(cm)
		tp_ / (tp_ + fp_)
	end

	function recall(cm)
		tp_, fn_ = tp(cm), fn(cm)
		tp_ / (tp_ + fn_)
	end

	function accuracy(cm)
		"""correct predictions / total predictions"""
		tp_, tn_ = tp(cm), tn(cm)
		fp_, fn_ = fp(cm), fn(cm)
		(tp_ + tn_) / (tp_ + fp_ + tn_ + fn_)
	end

	function error_rate(cm)
		1. - accuracy(cm)
	end

	function f₁_score(cm)
		tp_ = tp(cm)
		fp_, fn_ = fp(cm), fn(cm)
		2. * tp_ / (2. * tp_ + fp_ + fn_)
	end
end

# ╔═╡ f8501806-89fe-11eb-1703-b5665587eedb
begin
	cm = confusion_matrix(ŷ)
	@test sum(values(cm)) == length(test_messages)
end

# ╔═╡ 686d2c82-89ff-11eb-322c-e5ef7394e0f4
with_terminal() do
	@printf("tp: %4d / fp: %4d\n", tp(cm), fp(cm))
	@printf("tn: %4d / fn: %4d\n", tn(cm), fn(cm))
end

# ╔═╡ 2760a7e8-89fd-11eb-36d7-c5668ac807d5
(precision=precision(cm), recall=recall(cm),
	accuracy=accuracy(cm), f₁_score=f₁_score(cm))

# ╔═╡ 68c746b2-8a00-11eb-2762-ef6e77a8f8c3
function p_spam_given_token(model::NaiveBayes, token::String)::Float64
	prob_spam, prob_ham = probabilities(model, token)
	prob_spam / (prob_spam + prob_ham)
end

# ╔═╡ 9ec18d88-8a00-11eb-322c-e5ef7394e0f4
begin
	words = sort(collect(model.tokens),
		by=t -> p_spam_given_token(model, t), rev=false)

	with_terminal() do
		println("spammiest_words:")
		for w ∈ words[1:10]
			@printf("\t%15s\n", w)
		end

		println("\n\nhammiest_words:")
		for w ∈ words[end-10:end]
			@printf("%15s\n", w)
		end
	end
end

# ╔═╡ 661fff40-8a13-11eb-38c7-e99a5dda7c7a
model.tokens

# ╔═╡ Cell order:
# ╟─87e5e2ec-8915-11eb-362b-6bf13a36b8e4
# ╠═5eebac68-891b-11eb-37ba-c1d9481c6134
# ╟─cb25501a-8915-11eb-3626-a1ae6d233f92
# ╟─cb0838ae-8915-11eb-082f-3991192a101f
# ╟─caf03ada-8915-11eb-2935-01e110accba6
# ╟─cad66356-8915-11eb-0173-ada200678f63
# ╠═dc3b7280-8915-11eb-258a-1789476897b9
# ╠═dc209b9a-8915-11eb-0f63-ab3b7c33793c
# ╟─5a207b90-89d3-11eb-2a3e-655804b5bc8b
# ╠═6a8e15fc-89d3-11eb-3c28-43b3d0ad0e23
# ╟─d7aade44-891b-11eb-00d3-abe7406c4b09
# ╠═6dc9648a-89d6-11eb-35d0-e30bf7e3d3fc
# ╠═9ae99efe-891c-11eb-3d1f-89182e3f9ff8
# ╟─d7931368-891b-11eb-2899-f5e5a3f4bd35
# ╠═ea1f4a2c-891d-11eb-3027-bd447dd9d1e1
# ╟─50fe5114-891e-11eb-18ee-35d98ad32b7b
# ╠═4ffc3b48-895a-11eb-1ed2-f78b522d0ae9
# ╟─90f58d04-8958-11eb-170e-0f17904c9f2c
# ╠═123bcdb0-89d6-11eb-18f2-7f1acdb57a51
# ╟─dc0c5a5e-8915-11eb-2c62-5328b42add69
# ╟─dbefe77a-8915-11eb-39f1-5d0c84729739
# ╠═dbd86924-8915-11eb-351f-8362f09ba984
# ╟─c7410266-89d6-11eb-36b2-71e5df086e52
# ╠═d32b86e0-8928-11eb-193c-e3d7c85d23dd
# ╟─ab9ac56e-8969-11eb-3340-337957fd81b7
# ╟─dbbf22a2-8915-11eb-00eb-4b0278c0283d
# ╟─5f5e6d02-892a-11eb-32d8-dfc1a7f259fa
# ╠═073f6442-892e-11eb-0fc2-bbeb1d12619e
# ╠═6a454512-89e8-11eb-303a-ad5920ec581d
# ╠═aa4dae40-89ed-11eb-2cb6-9734911ab226
# ╟─2be56d9e-89ee-11eb-0e39-bd2174d78266
# ╠═21668a04-8a02-11eb-07b7-0f83275b0ad6
# ╠═aa3808b2-8a0e-11eb-32e8-e70ec61b55c7
# ╠═2bcd8a62-89ee-11eb-0575-dfa083fc40d6
# ╠═2bb0e580-89ee-11eb-24ae-a3947d35f204
# ╠═2b949b4e-89ee-11eb-0efa-9bcde1769743
# ╠═347e75f2-89f5-11eb-2a08-d9b2bdda1511
# ╠═5bcb82ee-89fa-11eb-36d7-c5668ac807d5
# ╠═82175c86-89fa-11eb-1703-b5665587eedb
# ╠═f8501806-89fe-11eb-1703-b5665587eedb
# ╠═686d2c82-89ff-11eb-322c-e5ef7394e0f4
# ╠═2760a7e8-89fd-11eb-36d7-c5668ac807d5
# ╠═68c746b2-8a00-11eb-2762-ef6e77a8f8c3
# ╠═9ec18d88-8a00-11eb-322c-e5ef7394e0f4
# ╠═661fff40-8a13-11eb-38c7-e99a5dda7c7a

max_stars_repo_path: Data Science From Scratch/13_NaiveBayes.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
function asap4(style = Rectangular(2,1); blacklist = ())
    arch = TopLevel{2}("asap4")
    arch.metadata["style"] = style

    # Extra parameters
    num_fifos = 2

    ####################
    # Normal Processor #
    ####################
    processor = build_processor_tile(style)
    for r in 0:23, c in 0:26
        in((r,c), blacklist) && continue
        add_child(arch, processor, CartesianIndex(r,c))
    end
    for r in 0:19, c in 27
        in((r,c), blacklist) && continue
        add_child(arch, processor, CartesianIndex(r,c))
    end

    ####################
    # Memory Processor #
    ####################
    memory_processor = build_processor_tile(style, include_memory = true)
    for r = 24, c = 0:26
        add_child(arch, memory_processor, CartesianIndex(r,c))
    end

    #################
    # 2 Port Memory #
    #################
    memory_2port = build_memory(2)
    for r = 25, c in (0, 2, 4, 6, 9, 11, 13, 15, 18, 20, 22, 24)
        add_child(arch, memory_2port, CartesianIndex(r,c))
    end

    #################
    # 1 Port Memory #
    #################
    memory_1port = build_memory(1)
    for r = 25, c in (8, 17)
        add_child(arch, memory_1port, CartesianIndex(r,c))
    end

    #################
    # Input Handler #
    #################
    input_handler = build_input_handler(1)
    for (r,c) ∈ zip((0,11,13,17), (-1, 28, -1, 28))
        add_child(arch, input_handler, CartesianIndex(r,c))
    end
    # for (r,c) ∈ zip((5,6,5,6), (-1,-1,28,28))
    #     add_child(arch, input_handler, CartesianIndex(r,c))
    # end
    ##################
    # Output Handler #
    ##################
    output_handler = build_output_handler(1)
    for (r,c) ∈ zip((0,11,13,17), (28, -1, 28, -1))
        add_child(arch, output_handler, CartesianIndex(r,c))
    end

    #######################
    # Global Interconnect #
    #######################
    connect_processors(arch, style)
    connect_io(arch, style)
    connect_memories(arch, style)
    return arch
end


text:
import CausalityToolsBase:
    BinningScheme,
    RectangularBinningScheme,
    RectangularBinning,
    TriangulationBinningScheme,
    TriangulationBinning,
    SimplexIntersectionType,
    ExactIntersection,
    ApproximateIntersection

import StaticArrays: SVector, MVector
import DelayEmbeddings: Dataset
import CausalityToolsBase: get_minima_and_edgelengths, encode

"""
    transferoperator(points, binning_scheme::RectangularBinning; kwargs...)

Discretize `points` using the provided `binning_scheme` and compute the transfer operator
over the partition elements.

# Example

Assume we have enough points that a rectangular partition yields a good estimate of the
invariant measure. Then the transfer operator over the partition can be computed as follows.

```julia
pts = [rand(3) for i = 1:2000]

# Use rectangular boxes constructed by subdividing each coordinate
# axis into 10 subintervals of equal length.
transferoperator(pts, RectangularBinning(10))
```
"""
function transferoperator(points, binning_scheme::RectangularBinning; kwargs...) end

########################
# Triangulation binnings
########################
"""
    transferoperator(pts, ϵ::TriangulationBinning,
        simplex_intersection_type::ApproximateIntersection;
        n::Int = 200, sample_randomly::Bool = false) -> TransferOperatorTriangulationApprox

Estimate the invariant measure over the state space defined by `pts` using a triangulation
of the phase space as the partition, using approximate simplex intersections to compute transition
probabilities between the states (simplices).

`n` is the number of points that each simplex is
sampled with, and `sample_randomly` indicates whether points used be sampled randomly
within each simpled (`sample_randomly = true`) or by a regular simplex splitting routine
(`sample_randomly = false`, which is default).

# Example

Assume we have sufficiently few points that a triangulation approach involving simplex
intersections is computationally feasible to compute the invariant measure. Then
a transfer operator can be computed as follows.

```julia
pts = [rand(3) for i = 1:30]
transferoperator(pts, TriangulationBinning(), ApproximateIntersection())
```
"""
function transferoperator(pts, ϵ::TriangulationBinning,
        simplex_intersection_type::ApproximateIntersection;
        n::Int = 200, sample_randomly::Bool = false)

    transferoperator_triangulation_approx(pts,
        n_sample_pts = n,
        sample_randomly = sample_randomly)
end

"""
    transferoperator(pts, ϵ::TriangulationBinning,
        simplex_intersection_type::ExactIntersection) -> TransferOperatorTriangulationExact

Estimate the invariant measure over the state space defined by `pts` using a triangulation
of the phase space as the partition, using exact simplex intersections to compute transition
probabilities between the states (simplices).

# Example

Assume we have sufficiently few points that a triangulation approach involving simplex
intersections is computationally feasible to compute the invariant measure. Then
a transfer operator can be computed as follows.

```julia
pts = [rand(3) for i = 1:30]
transferoperator(pts, TriangulationBinning(), ExactIntersection())
```
"""
function transferoperator(pts, ϵ::TriangulationBinning,
        simplex_intersection_type::ExactIntersection)

    transferoperator_triangulation_exact(pts)
end

########################
# Rectangular binnings
########################

function transferoperator(points::Vector{T}, binning_scheme::RectangularBinning;
    allocate_frac::Float64 = 1.0, boundary_condition = :none) where {T <:Union{SVector, MVector, Vector}}

    # Identify which bins of the partition resulting from using ϵ each
    # point of the embedding visits.
    mini, edgelengths = get_minima_and_edgelengths(points, binning_scheme)
    encoded_pts = encode(points, mini, edgelengths)

    #visited_bins = assign_bin_labels(points, binning_scheme.ϵ)

    # Which are the visited bins, which points
    # visits which bin, repetitions, etc...
    binvisits = get_binvisits(encoded_pts)

    # Use that information to estimate transfer operator
    estimate_transferoperator_from_binvisits(binvisits,
                        allocate_frac = allocate_frac,
                        boundary_condition = boundary_condition)
end

function transferoperator(points::AbstractArray{T, 2}, binning_scheme::RectangularBinning;
        allocate_frac::Float64 = 1.0, boundary_condition = :none) where {T}

    if size(points, 1) > size(points, 2)
        points = transpose(points)
    end

    pts = [points[:, i] for i = 1:maximum(size(points))]

    # Identify which bins of the partition resulting from using ϵ each
    # point of the embedding visits.
    mini, edgelengths = get_minima_and_edgelengths(pts, binning_scheme)
    encoded_pts = encode(pts, mini, edgelengths)

    # Which are the visited bins, which points
    # visits which bin, repetitions, etc...
    binvisits = get_binvisits(encoded_pts)

    # Use that information to estimate transfer operator
    estimate_transferoperator_from_binvisits(binvisits,
                        allocate_frac = allocate_frac,
                        boundary_condition = boundary_condition)
end

function transferoperator(points::Vector{Vector{T}}, binning_scheme::RectangularBinning;
    allocate_frac::Float64 = 1.0, boundary_condition = :none) where {T}

    transferoperator(hcat(points...,), binning_scheme, allocate_frac = allocate_frac,
        boundary_condition = boundary_condition)
end

function transferoperator(points::Vector{SVector{T}}, binning_scheme::RectangularBinning;
    allocate_frac::Float64 = 1.0, boundary_condition = :none) where {T}

    transferoperator(Array(hcat(points...,)), binning_scheme, allocate_frac = allocate_frac,
        boundary_condition = boundary_condition)
end

function transferoperator(points::Vector{MVector{T}}, binning_scheme::RectangularBinning;
    allocate_frac::Float64 = 1.0, boundary_condition = :none) where {T}

    transferoperator(Array(hcat(points...,)), binning_scheme, allocate_frac = allocate_frac,
        boundary_condition = boundary_condition)
end

function transferoperator(points::Dataset, binning_scheme::RectangularBinning;
    allocate_frac::Float64 = 1.0, boundary_condition = :none) where {T}

    transferoperator(transpose(Matrix(points)), binning_scheme, allocate_frac = allocate_frac,
        boundary_condition = boundary_condition)
end

########################
# Triangulation binnings
########################

"""
    transferoperator(points, binning_scheme::TriangulationBinningScheme,
        simplex_intersection::SimplexIntersectionType; kwargs...)

Discretize `points` using the provided `binning_scheme` and compute the transfer operator
over the partition elements using the type of simplex intersections indicated by
`simplex_intersection_type`.
"""
transferoperator(invariant_pts, binning_scheme::TriangulationBinning,
    simplex_intersection_type::SimplexIntersectionType; kwargs...)

export transferoperator


text:
"Logger configuration"
type LoggerConfig <: Filterable
    name::AbstractString
    level::LEVEL
    additive::Bool

    appenders::AppenderReferences
    parent::Nullable{LoggerConfig}
    event::FACTORY

    # Filterable
    filter::FILTER
    state::LifeCycle.State

    #TODO: properties::Dict{Property, Bool}
end

typealias LOGCONFIGS Dict{AbstractString, LoggerConfig}

# Constructors
function LoggerConfig(name::AbstractString;
                      level::LEVEL = LEVEL(),
                      filter::FILTER=FILTER(),
                      additive::Bool = false)
    return LoggerConfig(name, level, additive,
                        AppenderReferences(),
                        Nullable{LoggerConfig}(),
                        FACTORY(LOG4JL_LOG_EVENT),
                        filter, LifeCycle.INITIALIZED)
end
LoggerConfig(lvl::Level.EventLevel) = LoggerConfig(ROOT_LOGGER_NAME, level=LEVEL(lvl))
LoggerConfig() = LoggerConfig(LOG4JL_DEFAULT_STATUS_LEVEL)

"Start the logging configuration."
function start(lc::LoggerConfig)
    trace(LOGGER, "Starting $lc.")
    state!(lc, LifeCycle.STARTING)
    start(filter(lc))
    state!(lc, LifeCycle.STARTED)
    trace(LOGGER, "Started $lc OK.")
end

"Disable the logging configuration."
function stop(lc::LoggerConfig)
    trace(LOGGER, "Stopping $lc.")
    # change state
    state!(lc, LifeCycle.STOPPING)
    # stop and remove appenders
    for (apnm, ref) in appenders
        delete!(appenderRefs, apnm)
        stop(ref)
    end
    # stop filter
    stop(filter(lc))
    # change state
    state!(lc, LifeCycle.STOPPED)
    trace(LOGGER, "Stoped $lc OK.")
end

"Returns the logger name"
name(lc::LoggerConfig) = lc.name

"Returns the logging level"
level(lc::Nullable{LoggerConfig}) = isnull(lc) ? LOG4JL_DEFAULT_STATUS_LEVEL : level(get(lc))
level(lc::LoggerConfig) = get(lc.level, level(lc.parent))
level!(lc::LoggerConfig, lvl::Level.EventLevel) = (lc.level = lvl)

"Returns the value of the additive flag"
isadditive(lc::LoggerConfig) = lc.additive

"Logs an event"
function log(lc::LoggerConfig, evnt::Event)
    isfiltered(lc.filter, evnt) && return # Logger filters are execuded
    for ref in values(lc.appenders)
        append!(ref, evnt)
    end
    lc.additive && !isnull(lc.parent) && log(get(lc.parent), evnt)
end
function log(lc::LoggerConfig, logger, fqmn, level, marker, msg)
    log(lc, LOG4JL_LOG_EVENT(logger, fqmn, marker, level, msg)) #TODO: properties
end

show(io::IO, lc::LoggerConfig) = print(io, "LoggerConfig(", isempty(lc.name) ? "root" : lc.name, ":", level(lc) , ")")

"Returns all appender references"
references(lc::LoggerConfig) = values(lc.appenders)

"Adds an appender reference to configuration"
function reference!(lc::LoggerConfig, apndr::Appender, lvl::LEVEL=LEVEL(), flt::FILTER=FILTER())
    apn = name(apndr)
    lvl = get(lvl, Level.ALL)
    apndref = Appenders.Reference(apndr, level=lvl, filter=flt)
    lc.appenders[apn] = apndref
    return apndref
end

"Return parent logger configuration"
parent(lc::LoggerConfig) = lc.parent

"Set parent logger configuration"
parent!(lc::LoggerConfig, parent::Nullable{LoggerConfig}) = (lc.parent = parent)
parent!(lc::LoggerConfig, parent::LoggerConfig) = parent!(lc, Nullable(parent))


text:
include("expected_score_struct.jl")
mutable struct Constraint
    length_min::Int64
    length_max::Int64
    expected_score::ExpectedScore
    mean_vars::Vector{Symbol} # constrain the mean to be
    mean_vars_min::Vector{Float64}
    mean_vars_max::Vector{Float64}
    constr_A::Matrix{Float64}
    constr_b::Vector{Float64}
    # ol_max::Matrix{Int64}
    Constraint(
        length_min,
        length_max,
        expected_score,
        mean_vars,
        mean_vars_min,
        mean_vars_max,
        constr_A,
        constr_b,
    ) = new(
        length_min,
        length_max,
        expected_score,
        mean_vars,
        mean_vars_min,
        mean_vars_max,
        constr_A,
        constr_b,
    )
    Constraint() = new(
        zero(Int64),
        10000000,
        ExpectedScore(),
        Vector{Symbol}(undef, 0),
        Vector{Float64}(undef, 0),
        Vector{Float64}(undef, 0),
        Matrix{Float64}(undef, 0, 0),
        Vector{Float64}(undef, 0),
    )
end


text:
module TestModelMatching

using MLJ
using Test

X = (a = rand(5), b = categorical(1:5))
y = rand(5)
w = rand(5)

@test matching(X) == MLJ.ModelChecker{false,false,scitype(X),missing}()
@test matching(X, y) == MLJ.ModelChecker{true,false,scitype(X),scitype(y)}()
@test matching(X, y, w) == MLJ.ModelChecker{true,true,scitype(X),scitype(y)}()

@test !matching("RidgeRegressor", pkg="MultivariateStats", X)
@test matching("FeatureSelector", X)

m1 = models(matching(X))
@test issubset([info("FeatureSelector"),
               info("OneHotEncoder"),
                info("Standardizer")], m1)

@test !("PCA" in m1)
@test !(info("PCA") in m1)

m2 = models(matching(X, y))
@test info("ConstantRegressor") in m2
@test !(info("DecisionTreeRegressor") in m2)

end
true


text:
function substep{T}(a::Matrix{T}, beta::T, b::Matrix{T}, result::Matrix{T})
    for j=1:size(a,2)
        @simd for i=1:size(b,1)
            @inbounds result[i,j] = a[i,j] + beta*b[i,j]
        end
    end
end

function scale{T}(alpha::T, a::Matrix{T}, result::Matrix{T})
    for j=1:size(a,2)
        for i=1:size(a,1)
            @inbounds result[i,j] = alpha*a[i,j]
        end
    end
end

function oderk_step{T}(F::Function, t::T, h::T, x::Matrix{T}, a::Matrix{T},
                bs::Vector{T}, bp::Vector{T}, c::Vector{T},
                xs::Matrix{T}, xp::Matrix{T}, dx::Matrix{T}, tmp::Matrix{T},
                k::Vector{Matrix{T}})
    substep(x, h*bs[1], k[1], xs) #xs = x + h*bs[1]*k[1]
    substep(x, h*bp[1], k[1], xp) #xp = x + h*bp[1]*k[1]
    for j = 2:length(c)
        scale(a[j,1], k[1], dx) #dx = a[j,1]*k[1]
        for i = 2:j-1
            substep(dx, a[j,i], k[i], dx) #dx += a[j,i]*k[i]
        end
        #k[j] = F(t + h*c[j], x + h*dx)
        substep(x, h, dx, tmp)
        F(t + h*c[j], tmp, k[j])
        # compute the (p-1)th order estimate
        substep(xs, h*bs[j], k[j], xs) #xs = xs + h*bs[j]*k[j]
        # compute the pth order estimate
        substep(xp, h*bp[j], k[j], xp) #xp = xp + h*bp[j]*k[j]
    end
    # Estimate error
    substep(xs, complex(-1.), xp, tmp) #gamma1 = xs - xp
    return norm(reshape(tmp, length(tmp)), Inf)
end

const dp_coefficients = (5,
                         [    0           0          0         0         0        0
                              1/5         0          0         0         0        0
                              3/40        9/40       0         0         0        0
                             44/45      -56/15      32/9       0         0        0
                          19372/6561 -25360/2187 64448/6561 -212/729     0        0
                           9017/3168   -355/33   46732/5247   49/176 -5103/18656  0
                             35/384       0        500/1113  125/192 -2187/6784  11/84],
                         # 4th order b-coefficients
                         [5179/57600 0 7571/16695 393/640 -92097/339200 187/2100 1/40],
                         # 5th order b-coefficients
                         [35/384 0 500/1113 125/192 -2187/6784 11/84 0],
                         )
const c_compl = vec(complex(sum(dp_coefficients[2], 2)))   # consistency condition
const a_compl = complex(dp_coefficients[2])
const bs_compl = vec(complex(dp_coefficients[3]))
const bp_compl = vec(complex(dp_coefficients[4]))

const a2 = Float64[1/5]
const a3 = Float64[3/40, 9/40]
const a4 = Float64[44/45, -56/15, 32/9]
const a5 = Float64[19372/6561, -25360/2187, 64448/6561, -212/729]
const a6 = Float64[9017/3168, -355/33, 46732/5247, 49/176, -5103/18656]
const a7 = Float64[35/384, 0., 500/1113, 125/192, -2187/6784, 11/84]
const a = {Float64[] a2 a3 a4 a5 a6 a7}
const bp = a7
const bs = Float64[5179/57600, 0., 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]
const c = Float64[0., 1/5, 3/10, 4/5, 8/9, 1., 1.]

function substep2{T}(x::Vector{T}, x0::Vector{T}, h::Float64, coeffs::Vector{Float64}, k::Vector{Vector{T}})
    @inbounds for m=1:length(x0)
        dx::T = 0.
        @inbounds for i=1:length(coeffs)
            dx += coeffs[i]::Float64*k[i][m]
        end
        x[m] = x0[m] + h*dx
    end
    return nothing
end

function dopri_step3{T}(F::Function, t::Float64, h::Float64,
                x0::Vector{T}, xp::Vector{T}, xs::Vector{T}, k::Vector{Vector{T}})
    for i=2:length(c)
        substep2(xp, x0, h, a[i], k)
        F(t + h*c[i], xp, k[i])
    end
    substep2(xs, x0, h, bs, k)
    return nothing
end

function step1{T}(xp::Vector{T}, x::Vector{T}, h::Float64, k1::Vector{T})
    @inbounds for m=1:length(x)
        xp[m] = x[m] + h*a2[1]*k1[m]
    end
end

function step2{T}(xp::Vector{T}, x::Vector{T}, h::Float64, k1::Vector{T}, k2::Vector{T})
    @inbounds for m=1:length(x)
        dx = a3[1]*k1[m] + a3[2]*k2[m]
        xp[m] = x[m] + h*dx
    end
end

function step3{T}(xp::Vector{T}, x::Vector{T}, h::Float64, k1::Vector{T}, k2::Vector{T}, k3::Vector{T})
    @inbounds for m=1:length(x)
        dx = a4[1]*k1[m] + a4[2]*k2[m] + a4[3]*k3[m]
        xp[m] = x[m] + h*dx
    end
end

function step4{T}(xp::Vector{T}, x::Vector{T}, h::Float64, k1::Vector{T}, k2::Vector{T}, k3::Vector{T}, k4::Vector{T})
    @inbounds for m=1:length(x)
        dx = a5[1]*k1[m] + a5[2]*k2[m] + a5[3]*k3[m] + a5[4]*k4[m]
        xp[m] = x[m] + h*dx
    end
end

function step5{T}(xp::Vector{T}, x::Vector{T}, h::Float64, k1::Vector{T}, k2::Vector{T}, k3::Vector{T}, k4::Vector{T}, k5::Vector{T})
    @inbounds for m=1:length(x)
        dx = a6[1]*k1[m] + a6[2]*k2[m] + a6[3]*k3[m] + a6[4]*k4[m]
        dx += a6[5]*k5[m]
        xp[m] = x[m] + h*dx
    end
end

function step6{T}(xp::Vector{T}, x::Vector{T}, h::Float64, k1::Vector{T}, k2::Vector{T}, k3::Vector{T}, k4::Vector{T}, k5::Vector{T}, k6::Vector{T})
    @inbounds for m=1:length(x)
        dx = a7[1]*k1[m] + a7[3]*k3[m] + a7[4]*k4[m]
        dx += a7[5]*k5[m] + a7[6]*k6[m]
        xp[m] = x[m] + h*dx
    end
end

function error_estimates{T}(xp::Vector{T}, x::Vector{T}, h::Float64, k1::Vector{T}, k2::Vector{T}, k3::Vector{T}, k4::Vector{T}, k5::Vector{T}, k6::Vector{T})
    max_difference::Float64 = 0.
    max_value::Float64 = 0.
    @inbounds for m=1:length(x)
        dx = bs[1]*k1[m] + bs[3]*k3[m] + bs[4]*k4[m]
        dx += bs[5]*k5[m] + bs[6]*k6[m] + bs[7]*k7[m]
        v = abs2(x[m] + h*dx - xp[m])
        #v = abs2(xs[m]-xp[m])
        if v > max_difference
            max_difference = v
        end
        v = abs2(xp[m])
        if v > max_value
            max_value = v
        end
    end
    # println("DiffValue: ", sqrt(max_difference))
    # println("MValue: ", sqrt(max_value))
    # println("|xp|", norm(xp,Inf))
    # println("|x|", norm(x,Inf))
    return sqrt(max_value), sqrt(max_difference)
end

function dopri_step{T}(F::Function, t::Float64, h::Float64,
                x::Vector{T}, xp::Vector{T},
                k1::Vector{T}, k2::Vector{T}, k3::Vector{T}, k4::Vector{T}, k5::Vector{T}, k6::Vector{T}, k7::Vector{T})
    step1(xp, x, h, k1)
    F(t + h*c[2], xp, k2)
    step2(xp, x, h, k1, k2)
    F(t + h*c[3], xp, k3)
    step3(xp, x, h, k1, k2, k3)
    F(t + h*c[4], xp, k4)
    step4(xp, x, h, k1, k2, k3, k4)
    F(t + h*c[5], xp, k5)
    step5(xp, x, h, k1, k2, k3, k4, k5)
    F(t + h*c[6], xp, k6)
    step6(xp, x, h, k1, k2, k3, k4, k5, k6)
    F(t + h*c[7], xp, k7)
    return error_estimates(xp, x, h, k1, k2, k3, k4, k5, k6)
end

function dopri_step2{T}(F::Function, t::Float64, h::Float64,
                x::Vector{T}, xp::Vector{T},
                k1::Vector{T}, k2::Vector{T}, k3::Vector{T}, k4::Vector{T}, k5::Vector{T}, k6::Vector{T}, k7::Vector{T})
    # step 2
    #dim = length(x)
    @inbounds for m=1:length(x)
        xp[m] = x[m] + h*a2[1]*k1[m]
    end
    F(t + h*c[2], xp, k2)
    # step 3
    @inbounds for m=1:length(x)
        dx = a3[1]*k1[m] + a3[2]*k2[m]
        xp[m] = x[m] + h*dx
    end
    F(t + h*c[3], xp, k3)
    # step 4
    @inbounds for m=1:length(x)
        dx = a4[1]*k1[m] + a4[2]*k2[m] + a4[3]*k3[m]
        xp[m] = x[m] + h*dx
    end
    F(t + h*c[4], xp, k4)
    # step 5
    @inbounds for m=1:length(x)
        dx = a5[1]*k1[m] + a5[2]*k2[m] + a5[3]*k3[m] + a5[4]*k4[m]
        xp[m] = x[m] + h*dx
    end
    F(t + h*c[5], xp, k5)
    # step 6
    @inbounds for m=1:length(x)
        dx = a6[1]*k1[m] + a6[2]*k2[m] + a6[3]*k3[m] + a6[4]*k4[m]
        dx += a6[5]*k5[m]
        xp[m] = x[m] + h*dx
    end
    F(t + h*c[6], xp, k6)
    # step 7
    @inbounds for m=1:length(x)
        dx = a7[1]*k1[m] + a7[3]*k3[m] + a7[4]*k4[m]
        dx += a7[5]*k5[m] + a7[6]*k6[m]
        xp[m] = x[m] + h*dx
    end
    F(t + h*c[7], xp, k7)
    # Final
    max_difference::Float64 = 0.
    max_value::Float64 = 0.
    @inbounds for m=1:length(x)
        dx = bs[1]*k1[m] + bs[3]*k3[m] + bs[4]*k4[m]
        dx += bs[5]*k5[m] + bs[6]*k6[m] + bs[7]*k7[m]
        v = abs2(x[m] + h*dx - xp[m])
        if v > max_difference
            max_difference = v
        end
        v = abs2(xp[m])
        if v > max_value
            max_value = v
        end
    end
    # println("DiffValue: ", sqrt(max_difference))
    # println("MValue: ", sqrt(max_value))
    # println("|xp|", norm(xp,Inf))
    # println("|x|", norm(x,Inf))
    return sqrt(max_value), sqrt(max_difference)
end

function f(t::Float64, x::Vector{Complex128}, result::Vector{Complex128})
    for i=1:length(x)
        result[i] = cos(x[i]*t)
    end
end

const N = 800
srand(0)
const x0 = complex(rand(N^2))#Complex128[1., 0., 0. ,0.]
const x = zeros(Complex128, N^2)
const xs = zeros(Complex128, N^2)
const xp = zeros(Complex128, N^2)
const k1 = zeros(Complex128, N^2)
const k2 = zeros(Complex128, N^2)
const k3 = zeros(Complex128, N^2)
const k4 = zeros(Complex128, N^2)
const k5 = zeros(Complex128, N^2)
const k6 = zeros(Complex128, N^2)
const k7 = zeros(Complex128, N^2)

function g(t::Complex128, x::Matrix{Complex128}, result::Matrix{Complex128})
    for j=1:size(x,2), i=1:size(x,1)
        result[i,j] = cos(x[i,j]*t)
    end
end

const x0_ = reshape(deepcopy(x0), N, N)
const x_ = zeros(Complex128, N, N)
const dx_ = zeros(Complex128, N, N)
const tmp_ = zeros(Complex128, N, N)
const xs_ = zeros(Complex128, N, N)
const xp_ = zeros(Complex128, N, N)
const k1_ = zeros(Complex128, N, N)
const k2_ = zeros(Complex128, N, N)
const k3_ = zeros(Complex128, N, N)
const k4_ = zeros(Complex128, N, N)
const k5_ = zeros(Complex128, N, N)
const k6_ = zeros(Complex128, N, N)
const k7_ = zeros(Complex128, N, N)

const k_ = Matrix{Complex128}[k1_, k2_, k3_, k4_, k5_, k6_, k7_]

function profile_optimized(runs)
    t = 1.
    h = 0.1
    f(t,x0,k1)
    for i=1:runs
        dopri_step(f, t, h,
                    x0, xp,
                    k1, k2, k3, k4, k5, k6, k7)
        #println("regular |xp|: ", norm(xp,Inf))
        t += h
    end
end

k = Vector{Complex128}[]
push!(k, k1)
push!(k, k2)
push!(k, k3)
push!(k, k4)
push!(k, k5)
push!(k, k6)
push!(k, k7)

function profile_optimized2(runs)
    t = 1.
    h = 0.1
    f(t,x0,k1)
    for i=1:runs
        dopri_step3(f, t, h,
                    x0, xp, xs,
                    k)
        #println("regular |xp|: ", norm(xp,Inf))
        t += h
    end
end

function profile(runs)
    t = 1.+0.im
    h = 0.1+0.im
    g(t,x0_,k1_)
    for i=1:runs
        oderk_step(g, t, h,
                    x0_, a_compl, bs_compl, bp_compl, c_compl,
                    xs_, xp_, dx_, tmp_,
                    k_)
        #println("regular |xp_|: ", norm(reshape(xp_,N^2),Inf))
        t += h
    end
    return nothing
end

const runs = 1

println("=======Regular========")
@time profile(1)
@time profile(runs)

println("=======Optimized========")
@time profile_optimized2(1)
@time profile_optimized2(runs)

# println("=======Optimized2========")
# @time profile_optimized(1)
# @time profile_optimized(runs)

println(norm(xp-reshape(xp_,N^2)))


text:
using DataStructures
import DataStructures: dequeue!, enqueue!
import Base: length, next, start, done

# struct HalfPath
#     q::Queue{Int}
#     HalfPath() = new(Queue(Int))
# end
HalfPath = Deque{Int}

enqueue!(p::HalfPath, x::Int) = push!(p, x)
dequeue!(p::HalfPath) = shift!(p)

struct Path
    head::HalfPath
    tail::HalfPath
end
empty_path() = Path(HalfPath(), HalfPath())

Base.length(P::Path) = length(P.head) + length(P.tail)
Base.start(P::Path) = (false, start(reverse_iter(P.head)))
Base.done(P::Path, state) = state[1] && done(P.tail, state[2])
function Base.next(P::Path, state)
    in_tail, s = state
    if in_tail
        item, ns = next(P.tail, s)
    elseif done(reverse_iter(P.head), s)
        ns = start(P.tail)
        item, ns = next(P.tail, ns)
        in_tail = true
    else
        item, ns = next(reverse_iter(P.head), s)
    end
    return (item, (in_tail, ns))
end
Base.eltype(::Type{Path}) = Int


text:
"""
Wrapper around the C++ libtreelas
"""
module Cxx

import Libdl

const BUILD_DIR = joinpath(@__DIR__, "..", "deps", "build")

const lib = let BUILD_CONFIG = @static Sys.iswindows() ? "Release" : ""
    joinpath(
        BUILD_DIR,
        BUILD_CONFIG,
        (@static Sys.isunix() ? "lib" : "") * "treelas.$(Libdl.dlext)",
    )
end

function _find_python_extension()::String
    let regex = Regex("_treelas.+\\.$(Libdl.dlext)")
        _libs = filter(s -> match(regex, s) != nothing, readdir(BUILD_DIR))
    end
    @assert length(_libs) == 1
    return joinpath(BUILD_DIR, _libs[1])
end

function _find_all_tree_dp()
    ids = readlines(pipeline(`nm -D $(Cxx.lib)`, `grep tree_dp`, `awk '{print $3}'`))
    ids_cxx = map(s->readchomp(`c++filt $s`), ids)
    return ids, ids_cxx
end

function _find_tree_dp(; verbose=false)
    header = """
    const double*
    tree_dp<false, true>(
    const size_t n,
    double *x,
    const double *y,
    const int *parent,
    const double lam,
    const double mu,
    const int root
    );
    """

    header_lines = split(header, "\n")
    sig = findfirst(s -> occursin("tree_dp", s), header_lines)
    sig = header_lines[sig]
    ids, cxx = _find_all_tree_dp()
    i = findfirst(s -> occursin(sig, s), cxx)
    cxx = cxx[i]
    if verbose
        param_names = map(s -> replace(s, r".+[ *](\w+),?" => s"\1"),
                          header_lines[3:end-2])
        param_types = String.(split(replace(cxx, r".+\((.+)\)" => s"\1"), ", "))
        map!(param_types, param_types) do s::String
            for (a, b) in ["const" => "", " *" => "*", "unsigned long" => "size_t"]
                s = replace(s, a => b)
            end
            s
        end
        param_jltypes = map(param_types) do s
            if s[end] == '*'
                Ptr{eval(Symbol("C" * s[1:end-1]))}
            else
                eval(Symbol("C" * s))
            end
        end

        println("ccall(")
        println("    (:$(ids[i]), lib),")
        println("    Ref{Cdouble},")
        println("    (")
        for s in  param_types
            t = if s[end] == '*'
                "Ref{C" * s[1:end-1] * "}"
            else
                "C" * s
            end
            println("        ", t, ",")
        end
        println("    ),")
        for n in param_names
            println("    ", n, ",")
        end
        println(")")
    end
    return ids[i], cxx
end

function cxx_tree_dp(
    y::Array{Float64},
    parent::Vector{I},
    root::I,
    lam::Float64,
    mu::Float64 = 1.0,
) where {I}
    x = similar(y)
    @assert parent[root] == root
    n = length(y)
    @assert length(parent) == n
    parent = Vector{Cint}(copy(parent))
    parent .-= 1
    root -= 1
    ccall((:timer_disable, lib), Cvoid, ())
    ccall(
        (:tree_dp_f64_i32, lib),
        Ref{Cdouble},
        (
            Csize_t,
            Ref{Cdouble},
            Ref{Cdouble},
            Ref{Cint},
            Cdouble,
            Cdouble,
            Cint,
        ),
        n,
        x,
        y,
        parent,
        lam,
        mu,
        root,
    )
    return x
end

function cxx_tree_dual!(
    alpha::Array{Float64},
    parent::Vector{I},
    root::I,
) where {I}
    @assert parent[root] == root
    n = length(parent)
    @assert length(alpha) == n
    parent = Vector{Cint}(copy(parent))
    parent .-= 1
    root -= 1
    x = copy(alpha)
    post_order = Ptr{Cint}(0)
    ccall((:timer_disable, lib), Cvoid, ())
    ccall((:tree_dual_f64_i32, lib),
        Ref{Cdouble},
        (
            Csize_t,
            Ref{Cdouble},
            Ref{Cint},
            Ptr{Cint},
            Ref{Cdouble},
            Cint,
            Cuchar,
        ),
        n,
        x,
        parent,
        post_order,
        alpha,
        root,
        1,  # tree_orientation
    )
    return vec(alpha)
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
using Lynx
using Luxor

Lynx.init("Basic Example", 800, 600)

function setup()
    @info "Canvas width is $(@width) and height is $(@height)"
end

time = 0.

function update(dt)
    background("#111")
    sethue("yellowgreen")
    origin()

    radius = 50 * (cos(time) + 2)
    circle(O, radius, :fill)

    global time += dt
end

run!(update, setup)


text:
function propertynames(output::T) where {T<:Union{HelmholtzOutput,LaplaceOutput}}
    pg  = output.pg
    pgt = output.pgt

    fieldstarg = [:pottarg,:gradtarg,:hesstarg]
    fieldstarg = fieldstarg[1:pgt]

    if pg === nothing
        return tuple(fieldstarg...,)
    end

    fields = [:pot,:grad,:hess]
    fields = fields[1:pg]

    return tuple(cat(fields,fieldstarg,dims=1)...)

end

function propertynames(output::StokesOutput)
    ifppreg  = output.ifppreg
    ifppregtarg = output.ifppregtarg

    fieldstarg = [:pottarg,:pretarg,:gradtarg]
    fieldstarg = fieldstarg[1:ifppregtarg]

    if ifppreg === nothing
        return tuple(fieldstarg...,)
    end

    fields = [:pot,:pre,:grad]
    fields = fields[1:ifppreg]

    return tuple(cat(fields,fieldstarg,dims=1)...)

end

function propertynames(output::MaxwellOutput)
    ifeh  = output.ifeh
    ifehtarg = output.ifehtarg

    fieldstarg = [:etarg,:gradetarg,:htarg,:gradhtarg]
    if ifehtarg <= 2
        fieldstarg = fieldstarg[1:ifehtarg]
    elseif ifehtarg <= 4
        fieldstarg = fieldstarg[3:ifehtarg]
    elseif ifehtarg == 5
        fieldstarg = [:etarg,:htarg]
    end

    fields = [:e,:grade,:h,:gradh]
    if ifeh === nothing
        return tuple(fieldstarg...,)
    elseif ifeh <= 2
        fields = fields[1:ifeh]
    elseif ifeh <= 4
        fields = fields[3:ifeh]
    elseif ifeh == 5
        fields = [:e,:h]
    end

    return tuple(cat(fields,fieldstarg,dims=1)...)

end

# Probably not the prettiest print. Better than nothing
function show(io::IO, ::MIME"text/plain", output::T) where {T<:FMMVals}
    for field in propertynames(output)
        strfield = string(field)
        println("."*strfield*" "^(max(9-length(strfield),0)),
                '\t', "$(typeof(getfield(output,field)))")
    end
end


text:
# Generate observed Born data for Marmousi examples
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using ClusterManagers, Distributed, ArgParse; addprocs_slurm(2)

### Process command line args
function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--time_order", "-t"
            help = "Time order"
            arg_type = Int
            default = 2
        "--velmod_dir"
            help = "Velocity model directory"
            default = ""
        "--dobs_dir"
            help = "Outpur data directory"
            default = ""
    end
    return parse_args(s)
end

@everywhere using JUDI, SegyIO, HDF5, JLD, LinearAlgebra
parsed_args = parse_commandline()

time_order = parsed_args["time_order"]
velmod_dir = parsed_args["velmod_dir"]
dobs_dir   = parsed_args["dobs_dir"]

# Load migration velocity model
n1, d1, o1, m, m0, rho, rho0, qp, qp0, dm = read(h5open(velmod_dir, "r"), "n", "d", "o", "m", "m0", "rho", "rho0", "qp", "qp0", "dm")

# Set up model structure
n=(n1[1], n1[2]);
d=(d1[1], d1[2]);
o=(o1[1], o1[2]);

model = Model(n, d, o, m0, rho0, qp0);
dm = vec(dm)

nsrc = 2
xsrc = convertToCell(range(model.o[1], stop=model.n[1]*model.d[1], length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(30f0, stop=30f0, length=nsrc))
# receiver sampling and recording time
timeR = 4000f0   # receiver recording time [ms]
dtR = 4f0    # receiver sampling interval
# Source sampling and number of time steps
timeS = 4000f0
dtS = 4f0
# nxrec = 369
nxrec = model.n[1] - 2
# xrec = collect(range(model.o[1], stop=model.n[1]*model.d[1], length=nxrec))
xrec = collect(range(model.d[1], stop=(model.n[1]-2)*model.d[1], length=nxrec))
yrec = collect(range(0f0, stop=0f0, length=nxrec))
zrec = collect(range(30f0, stop=30f0, length=nxrec))
# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)
# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = 0.015f0     # kHz
wavelet = ricker_wavelet(timeS, dtS, f0)

q = judiVector(srcGeometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

opt = Options(f0=f0)

Pr = judiProjection(recGeometry)
F = judiModeling(model)
Ps = judiProjection(srcGeometry)
J = judiJacobian(Pr*F*adjoint(Ps), q)

# Born modeling
dobs = J*dm
block_out = judiVector_to_SeisBlock(dobs, q; source_depth_key="SourceDepth")
segy_write(dobs_dir, block_out)


text:
module Functional

using ..NeuralAttentionlib: @imexport

@imexport import ..NeuralAttentionlib:
    generic_qkv_attention, generic_multihead_qkv_attention,
    naive_qkv_attention, multihead_qkv_attention,
    mixing, weighted_sum_mixing, attention_score,
    scaled_dot_product_score, dot_product_score,
    masked_score, normalized_score

import ..NeuralAttentionlib:
    split_head, move_head_dim_out_perm, move_head_dim_out,
    move_head_dim_in_perm, move_head_dim_in, merge_head

using ..NeuralAttentionlib: SymLengthMask, BiLengthMask

"""
    generic_qkv_attention(mixingf, scoref, q, k, v, args...)

Generic version of [`naive_qkv_attention`](@ref). Need to specify mixing and score function.
"""
generic_qkv_attention

"""
    generic_multihead_qkv_attention(mixingf, scoref, head, q, k, v, args...)

Generic version of [`multihead_qkv_attention`](@ref). Need to specify mixing and score function.
"""
generic_multihead_qkv_attention

"""
    multihead_qkv_attention(head, q, k, v, mask=nothing)

Multihead version of [`naive_qkv_attention`](@ref). The core operation for implement a regular transformer layer.
"""
multihead_qkv_attention

@doc raw"""
    naive_qkv_attention(q, k, v, mask=nothing)

The scaled dot-product attention of a regular transformer layer.

``Attention(Q, K, V) = softmax(\frac{QK^T}{\sqrt{d_k}})V``

It's equivalent to `generic_qkv_attention(weighted_sum_mixing, normalized_score(NNlib.softmax) $ masked_score(mask) $ scaled_dot_product_score, q, k, v)`.

#Example

```julia
julia> fdim, ldim, bdim = 32, 10, 4;

julia> x = randn(fdim, ldim, bdim);

julia> y = naive_qkv_attention(x, x, x); # simple self attention

# no mask here
julia> z = generic_qkv_attention(weighted_sum_mixing, normalized_score(NNlib.softmax) $ scaled_dot_product_score, x, x, x);

julia> y ≈ z
true

```

See also: [`generic_qkv_attention`](@ref)
"""
naive_qkv_attention

"""
    mixing(f, v, g, args...) = f(attention_score(g, args...), v)

`Mixing` function api. Can be overload for doing custom implementation with [`generic_qkv_attention`](@ref).
 `f` is the mixing function and `g` is score function.

See also: [`generic_qkv_attention`](@ref), [`generic_multihead_qkv_attention`](@ref), [`attention_score`](@ref)
"""
mixing

"""
    weighted_sum_mixing(s, v)

The mixing function of a regular transformer layer. `s` is the attention score and `v` is the value of QKV attention.
"""
weighted_sum_mixing

"""
    attention_score(f, args...) = f(args...)

Attention score api. Can be overload for doing custom implementation with [`generic_qkv_attention`](@ref).
 `f` is the score function.

See also: [`generic_qkv_attention`](@ref), [`generic_multihead_qkv_attention`](@ref), [`mixing`](@ref)
"""
attention_score

@doc raw"""
    normalized_score(norm) = normalized_score $ norm
    normalized_score(norm, score, args...)

Normalized attenion score api. `norm` is the normalize function (like `softmax`) and `score` is the function
 that compute attention score from `args...`.

See also: [`naive_qkv_attention`](@ref)
"""
normalized_score

@doc raw"""
    masked_score(mask) = masked_score $ mask
    masked_score(maskop::AbstractAttenMaskOp, mask::AbstractAttenMask, score, args...)

Masked attention score api. Applying the `mask` according to `maskop` on the attention score
 compute from `score(args...)`.

See also: [`naive_qkv_attention`](@ref), [`SymLengthMask`](@ref), [`BiLengthMask`](@ref)
"""
masked_score

@doc raw"""
     scaled_dot_product_score(q, k, s = sqrt(inv(size(k, 1))))

The scaled dot-product attention score function of a regular transformer layer.

``Score(Q, K) = \frac{QK^T}{\sqrt{d_k}}``

See also: [`naive_qkv_attention`](@ref)
"""
scaled_dot_product_score

"""
    dot_product_score(q, k)

Dot-product attention score function. Equivalent to `scaled_dot_product_score(q, k, 1)`.

See also: [`scaled_dot_product_score`](@ref)
"""
dot_product_score

"""
    split_head(head::Int, x)

Split the first dimension into `head` piece of small vector. Equivalent to
 `reshape(x, :, head, tail(size(x))...)`.
"""
split_head

"""
    move_head_dim_out_perm(x::AbstractArray{T, N}, nobatch=false)
    move_head_dim_out_perm(N::Int, nobatch=false)

Dimension order for `permutedims` to move the `head` dimension (created by [`split_head`](@ref)) to batch dimension.
 Return a tuple of integer of length `n`. `nobatch` specify where `x` is a batch of data.

# Example

```julia
julia> Functional.move_head_dim_out_perm(5, false)
(1, 3, 4, 2, 5)

julia> Functional.move_head_dim_out_perm(5, true)
(1, 3, 4, 5, 2)

```

See also: [`split_head`](@ref), [`move_head_dim_out`](@ref)
"""
move_head_dim_out_perm

"""
    move_head_dim_out(x::AbstractArray, nobatch=false)

Equivanlent to `permutedims(x, move_head_dim_out_perm(x, nobatch)))`

See also: [`split_head`](@ref), [`move_head_dim_out_perm`](@ref)
"""
move_head_dim_out

"""
    move_head_dim_in_perm(x::AbstractArray{T, N}, nobatch=false)
    move_head_dim_in_perm(N::Int, nobatch=false)

Dimension order for `permutedims` to move the `head` dimension (created by [`split_head`](@ref)) from batch dimension
 to feature dimension (for [`merge_head`](@ref)). Return a tuple of integer of length `n`.
 `nobatch` specify where `x` is a batch of data.

# Example

```julia
julia> Functional.move_head_dim_in_perm(5, false)
(1, 4, 2, 3, 5)

julia> Functional.move_head_dim_in_perm(5, true)
(1, 5, 2, 3, 4)

```

See also: [`merge_head`](@ref), [`move_head_dim_in`](@ref)
"""
move_head_dim_in_perm

"""
    move_head_dim_in(x::AbstractArray, nobatch=false)

Equivanlent to `permutedims(x, move_head_dim_in_perm(x, nobatch)))`

See also: [`merge_head`](@ref), [`move_head_dim_in_perm`](@ref)
"""
move_head_dim_in

"""
    merge_head(x)

merge the `head` dimension split by [`split_head`](@ref).
"""
merge_head

end


text:
# defines the functions which actually run simulations

export rebound_time, # outputs a single rebound time
    viral_rebound_times, # rebound vector as a funciton of a single viral population
    trial_rebound_times, # matrix of rebound times
    trial_rebound_prob, # function of a theta vector
	virus_time_trace # full trace for drawing rebound trajectories

function wt_esc_counts(vp)
	# note that this name is exactly backward
	# returns a tuple
	esc::Int = 0
	for virus in vp.pop
		if virus.escape
			esc += 1
		end
	end
	(esc, length(vp.pop) - esc)
end

function virus_time_trace(vp::ViralPopulation, time_points; breakpoint = 5.0)
	# records the viral state (wt,esc) counts as a tuple
	# breakpoint sets a threshold at which to end the simulation and save on time
	# because escape has already been acieved
	# Default traces are set way above capacity, assuring that the trace has
	# the same length as time_points.
	out = Tuple{Int64, Int64}[]
	for time in time_points
		run_viruses_stop(vp, time)
		(esc, wt) = wt_esc_counts(vp)
		push!(out, (esc,wt))
		if esc/vp.capacity > breakpoint
			break
		end
	end
	return out
end

function rebound_time(vp::ViralPopulation; breakpoint = 0.8)
	# function to get a single rebound time and trace
	# breakpoint is a keyword that determines when a simulation should be treated deterministically
	# it is set very conservatively here. The main purpose is to prevent running simulations
	# at full capacity all the way to the endpoint time (56 days)
	# returns a tuple of rebound time and trace if required.
	(esc, _) = wt_esc_counts(vp)
	trace = virus_time_trace(vp, 0:56; breakpoint)
	if esc/vp.capacity > breakpoint
		t_rebound = 0.0
	elseif trace[end][1]/vp.capacity < breakpoint
		t_rebound = 57
	else
		ind = findfirst(col->col[1]/vp.capacity > breakpoint, trace)
		t_cross = (0:56)[ind]
		q_cross = min(trace[ind][1]/vp.capacity,1)
		t_rebound = 3*log(1 + exp(t_cross / 3) * (1-q_cross)/q_cross)
	end
	return (t_rebound,trace)
end

function viral_rebound_times(vp::ViralPopulation;
    breakpoint = 0.8,
    n_samples = 10)
	# Simulation with reused viral population structs
	# Single threaded loop used to rerun a population multiple times
    rebounds = Float64[]
    for _ in 1:n_samples
        (t_rebound, _) = rebound_time(vp::ViralPopulation; breakpoint)
        push!(rebounds, t_rebound)
        restart_viral_population!(vp)
    end
    return rebounds
end

function trial_rebound_times(theta_vec, ab_profile_list;
	# ab_profile_list is a list of lists of antibody mutation and selection properties
	# the lengths of the sublists determines the escape function in the dynamics.
    n_samples = 100, # number of samples per bayesian realization * n_runs
	n_runs = 10, # parameter for better thread efficiency
    breakpoint = 0.8, pars...)
    times = zeros(length(theta_vec)*n_runs, n_samples) # float matrix of rebound times
	Threads.@threads for (ii, θ) in collect(enumerate(repeat(theta_vec, n_runs)))
        vp = initialize_viral_population(θ, ab_profile_list; pars...)
        times[ii,:] .= viral_rebound_times(vp; n_samples,  breakpoint)
    end
    return reshape(times,(length(theta_vec), n_runs*n_samples))
end

function trial_rebound_prob(theta_vec, ab_profile_list;
    breakpoint = 0.1, kwds...)
    mean(vec(trial_rebound_times(theta_vec, ab_profile_list; breakpoint, kwds...)) .< 56.0)
end

max_stars_repo_name: StatPhysBio/EscapeSimulator

max_issues_repo_name: StatPhysBio/EscapeSimulator

max_forks_repo_name: StatPhysBio/EscapeSimulator