text:
"""
$(TYPEDEF)

Structure containig plot information for vtk.js.
"""
mutable struct PlutoVTKPlot  <: AbstractPlutoVistaBackend
    # command list passed to javascript
    jsdict::Dict{String,Any}

    # size in screen coordinates
    w::Float64
    h::Float64

    # update of a already created plot ?

    update::Bool

    args

    # uuid for identifying html element
    uuid::UUID
    PlutoVTKPlot(::Nothing)=new()
end

"""
$(SIGNATURES)

Create a vtk plot with given resolution in the notebook.
"""
function PlutoVTKPlot(;resolution=(300,300), kwargs...)
    p=PlutoVTKPlot(nothing)
    p.uuid=uuid1()
    p.jsdict=Dict{String,Any}("cmdcount" => 0,"cbar" => 0)
    p.w=resolution[1]
    p.h=resolution[2]

    default_args=(title="",
                  titlefontsize=12,
                  axisfontsize=10,
                  tickfontsize=10,
                  xlabel="x",
                  ylabel="y",
                  zlabel="z",
                  aspect=1.0,
                  zoom=1.0,
                  legendfontsize=10,
                  colorbarticks=:default,
                  clear=false)

    p.args=merge(default_args,kwargs)

    p.update=false
    p
end

const canvascolorbar = read(joinpath(@__DIR__, "..", "assets", "canvascolorbar.js"), String)
const plutovtkplot = read(joinpath(@__DIR__, "..", "assets", "plutovtkplot.js"), String)

"""
    Base.show(io::IO,::MIME"text/html",p::PlutoVTKPlot)

Show plot in html. This creates a vtk.js based renderer along with a canvas
for handling the colorbar.
"""
function Base.show(io::IO, ::MIME"text/html", p::PlutoVTKPlot)
    uuidcbar="$(p.uuid)"*"cbar"
    div=""
    if !p.update
    div="""

        <canvas id="$(uuidcbar)" width=60, height="$(p.h-25)"  style="display: inline-block; "></canvas>

    """
    end
    result="""

        """
     p.update=true
     write(io,result*div)
end

"""
    axis3d!(vtkplot)
Add 3D coordinate system axes to the plot.
Sets camera handling to 3D mode.
"""
function axis3d!(p::PlutoVTKPlot; kwargs...)
    args=merge(p.args,kwargs)

    command!(p,"axis")
    parameter!(p,"axisfontsize"   ,kwargs[:axisfontsize])
    parameter!(p,"tickfontsize"   ,kwargs[:tickfontsize])
    parameter!(p,"zoom"   ,kwargs[:zoom])
    parameter!(p,"cam","3D")
    parameter!(p,"xlabel",args[:xlabel])
    parameter!(p,"ylabel",args[:ylabel])
    parameter!(p,"zlabel",args[:zlabel])

end

"""
    axis2d!(vtkplot)
Add 2D coordinate system axes to the plot.
Sets camera handling to 2D mode.
"""
function axis2d!(p::PlutoVTKPlot; kwargs...)
    args=merge(p.args,kwargs)

    command!(p,"axis")
    parameter!(p,"axisfontsize"   ,kwargs[:axisfontsize])
    parameter!(p,"tickfontsize"   ,kwargs[:tickfontsize])
    parameter!(p,"zoom"   ,kwargs[:zoom])
    parameter!(p,"xlabel",args[:xlabel])
    parameter!(p,"ylabel",args[:aspect]==1  ?  args[:ylabel] : "$(args[:ylabel])*$(args[:aspect])" )
    parameter!(p,"zlabel",args[:zlabel])
    parameter!(p,"cam","2D")
end

"""
       vtkpolys(tris; offset=0)
Set up  polygon (triangle) data for vtk.
Coding is   [3, i11, i12, i13,   3 , i21, i22 ,i23, ...]
Careful: js indexing counts from zero.
"""
function vtkpolys(tris; offset=0)
   ipoly=1
    ntri=size(tris,2)
    off=offset-1
    polys=Vector{Int32}(undef,4*ntri)
    for itri=1:ntri
        polys[ipoly] = 3
        polys[ipoly+1] = tris[1,itri]+off
        polys[ipoly+2] = tris[2,itri]+off
        polys[ipoly+3] = tris[3,itri]+off
        ipoly+=4
    end
    polys
end

"""
    outline!(p::PlutoVTKPlot,pts,faces,facemarkers,facecolormap,nbregions,xyzmin,xyzmax;alpha=0.1)

Plot transparent outline of grid boundaries.
"""
function outline!(p::PlutoVTKPlot,pts,faces,facemarkers,facecolormap,nbregions,xyzmin,xyzmax;alpha=0.1)
    bregpoints0,bregfacets0=GridVisualize.extract_visible_bfaces3D(pts,faces,facemarkers,nbregions,
                                                                   xyzmax,
                                                                   primepoints=hcat(xyzmin,xyzmax)
                                                                   )
    bregpoints=hcat([reshape(reinterpret(Float32,bregpoints0[i]),(3,length(bregpoints0[i]))) for i=1:nbregions]...)
    bregfacets=vcat([vtkpolys(reshape(reinterpret(Int32,bregfacets0[i]),(3,length(bregfacets0[i]))),
                              offset= ( i==1 ? 0 : sum(k->length(bregpoints0[k]),1:i-1) ) )
                     for i=1:nbregions]...)
    bfacemarkers=vcat([fill(i,length(bregfacets0[i])) for i=1:nbregions]...)

    if typeof(facecolormap)==Symbol
        facecmap=colorschemes[facecolormap]
    else
        facecmap=ColorScheme(facecolormap)
    end
    facergb=reinterpret(Float64,get(facecmap,bfacemarkers,(1,size(facecmap))))
    nfaces=length(facergb)÷3
    facergba=zeros(UInt8,nfaces*4)
    irgb=0
    irgba=0
    for i=1:nfaces
        facergba[irgba+1]=UInt8(floor(facergb[irgb+1]*255))
        facergba[irgba+2]=UInt8(floor(facergb[irgb+2]*255))
        facergba[irgba+3]=UInt8(floor(facergb[irgb+3]*255))
        facergba[irgba+4]=UInt8(floor(alpha*255))
        irgb+=3
        irgba+=4
    end
    parameter!(p,"opolys",bregfacets)
    parameter!(p,"opoints",vec(bregpoints))
    parameter!(p,"ocolors",facergba)

end

"""
$(SIGNATURES)

Plot piecewise linear function on  triangular grid given as "heatmap".
Isolines can be given as a number or as a range.
"""
function tricontour!(p::PlutoVTKPlot, pts, tris,f;kwargs...)

    default_args=(colormap=:viridis, levels=0, limits=:auto,gridscale=1.0)
    args=merge(p.args,default_args)
    args=merge(args,kwargs)

    colormap=args[:colormap]

    p.jsdict=Dict{String,Any}("cmdcount" => 0)

    command!(p,"tricontour")

    levels,crange,colorbarticks=GridVisualize.isolevels(args,f)
    crange=Float64.(crange)

    parameter!(p,"points",vec(vcat(pts,zeros(eltype(pts),length(f))')))
    parameter!(p,"polys",vtkpolys(tris))
    parameter!(p,"gridscale",args[:gridscale])

    rgb=reinterpret(Float64,Base.get(colorschemes[colormap],f,crange))
    parameter!(p,"colors",UInt8.(floor.(rgb*255)))

    parameter!(p,"isopoints","none")
    parameter!(p,"isolines","none")
    parameter!(p,"aspect",args[:aspect])

    iso_pts=GridVisualize.marching_triangles(pts,tris,f,collect(levels))
    niso_pts=length(iso_pts)
    iso_pts=vcat(reshape(reinterpret(Float32,iso_pts),(2,niso_pts)),zeros(niso_pts)')
    iso_lines=Vector{UInt32}(undef,niso_pts+Int32(niso_pts//2))
    iline=0
    ipt=0
    for i=1:niso_pts//2
        iso_lines[iline+1]=2
        iso_lines[iline+2]=ipt
        iso_lines[iline+3]=ipt+1
        iline=iline+3
        ipt=ipt+2
    end
    parameter!(p,"isopoints",vec(iso_pts))
    parameter!(p,"isolines",iso_lines)

    # It seems a colorbar is best drawn via canvas...
    # https://github.com/Kitware/vtk-js/issues/1621
    bar_stops=collect(0:0.01:1)
    bar_rgb=reinterpret(Float64,get(colorschemes[colormap],bar_stops,(0,1)))
    bar_rgb=UInt8.(floor.(bar_rgb*255))
    p.jsdict["cbar"]=1
    p.jsdict["cbar_stops"]=bar_stops
    p.jsdict["cbar_colors"]=bar_rgb
    p.jsdict["cbar_levels"]=collect(colorbarticks)
    p.jsdict["cbar_fontsize"]=args[:legendfontsize]

    axis2d!(p; args...)
    p
end

"""
$(SIGNATURES)

Plot piecewise linear function on  triangular grid created from the tensor product of X and Y arrays as "heatmap".
Levels can be given as a number or as a range.
"""
contour!(p::PlutoVTKPlot,X,Y,f; kwargs...)=tricontour!(p,triang(X,Y)...,vec(f);kwargs...)

"""
$(SIGNATURES)

Plot isosurfaces given by `levels` and contour maps on planes given by the `*planes` parameters
for piecewise linear function on  tetrahedral mesh.
"""
function tetcontour!(p::PlutoVTKPlot, pts, tets,func; kwargs...)

    default_args=(colormap=:viridis,
                  levels=5,
                  limits=:auto,
                  faces=nothing,
                  facemarkers=nothing,
                  facecolormap=nothing,
                  xplanes=[prevfloat(Inf)],
                  yplanes=[prevfloat(Inf)],
                  zplanes=[prevfloat(Inf)],
                  levelalpha=0.25,
                  outlinealpha=0.1)
    args=merge(p.args,default_args)
    args=merge(args,kwargs)

    levels,crange=GridVisualize.isolevels(args,func)

    colormap=args[:colormap]
    faces=args[:faces]
    facemarkers=args[:facemarkers]
    facecolormap=args[:facecolormap]

    p.jsdict=Dict{String,Any}("cmdcount" => 0)
    command!(p,"tetcontour")
    xyzmin=zeros(3)
    xyzmax=ones(3)

    nbregions= facemarkers==nothing ? 0 :  maximum(facemarkers)

    if faces!=nothing && nbregions==0
        nbregions=1
        facemarkers=fill(Int32(1),size(faces,2))
    end

    if facecolormap==nothing
        facecolormap=GridVisualize.bregion_cmap(nbregions)
    end

    @views for idim=1:3
        xyzmin[idim]=minimum(pts[idim,:])
        xyzmax[idim]=maximum(pts[idim,:])
    end

    xplanes=args[:xplanes]
    yplanes=args[:yplanes]
    zplanes=args[:zplanes]

    cpts0,faces0,values=GridVisualize.marching_tetrahedra(pts,tets,func,
                                                          primepoints=hcat(xyzmin,xyzmax),
                                                          primevalues=crange,
                                                          GridVisualize.makeplanes(xyzmin,xyzmax,xplanes,yplanes,zplanes),
                                                          levels;
                                                          tol=0.0
                                                          )

    cfaces=reshape(reinterpret(Int32,faces0),(3,length(faces0)))
    cpts=copy(reinterpret(Float32,cpts0))
    parameter!(p,"points",cpts)
    parameter!(p,"polys",vtkpolys(cfaces))
    nan_replacement=0.5*(crange[1]+crange[2])
    for i=1:length(values)
        if isnan(values[i]) || isinf(values[i])
            values[i]=nan_replacement
        end
    end
    # nan_replacement=0.0
    # for i=1:length(cpts)
    #     if isnan(cpts[i]) || isinf(cpts[i])
    #         cpts[i]=nan_replacement
    #     end
    # end

    crange=Float64.(crange)
    rgb=reinterpret(Float64,get(colorschemes[colormap],values,crange))

    if args[:levelalpha]>0
        nfaces=length(rgb)÷3
        rgba=zeros(UInt8,nfaces*4)
        irgb=0
        irgba=0
        for i=1:nfaces
            rgba[irgba+1]=UInt8(floor(rgb[irgb+1]*255))
            rgba[irgba+2]=UInt8(floor(rgb[irgb+2]*255))
            rgba[irgba+3]=UInt8(floor(rgb[irgb+3]*255))
            rgba[irgba+4]=UInt8(floor(args[:levelalpha]*255))
            irgb+=3
            irgba+=4
        end
        parameter!(p,"transparent",1)
        parameter!(p,"colors",rgba)
    else
        parameter!(p,"transparent",0)
        parameter!(p,"colors",UInt8.(floor.(rgb*255)))
    end

    # It seems a colorbar is best drawn via canvas...
    # https://github.com/Kitware/vtk-js/issues/1621
    bar_stops=collect(0:0.01:1)
    bar_rgb=reinterpret(Float64,get(colorschemes[colormap],bar_stops,(0,1)))
    bar_rgb=UInt8.(floor.(bar_rgb*255))
    p.jsdict["cbar"]=1
    p.jsdict["cbar_stops"]=bar_stops
    p.jsdict["cbar_colors"]=bar_rgb
    p.jsdict["cbar_levels"]=vcat([crange[1]],levels,[crange[2]])
    p.jsdict["cbar_fontsize"]=args[:legendfontsize]

    if args[:outlinealpha]>0 && faces!=nothing
        parameter!(p,"outline",1)
        outline!(p,pts,faces,facemarkers,facecolormap,nbregions,xyzmin,xyzmax;alpha=args[:outlinealpha])
    else
        parameter!(p,"outline",0)
    end

    axis3d!(p; args...)
    p

end

"""
$(SIGNATURES)

Plot  triangular grid with optional region and boundary markers.
"""
function trimesh!(p::PlutoVTKPlot,pts, tris; kwargs...)

    default_args=(markers=nothing,
                  colormap=nothing,
                  edges=nothing,
                  gridscale=1.0,
                  edgemarkers=nothing,
                  edgecolormap=nothing)

    args=merge(p.args,default_args)
    args=merge(args,kwargs)
    parameter!(p,"aspect",args[:aspect])

    colormap=args[:colormap]
    markers=args[:markers]
    edgemarkers=args[:edgemarkers]
    edgecolormap=args[:edgecolormap]
    edges=args[:edges]

    ntri=size(tris,2)
    command!(p,"trimesh")

    cminmax=extrema(pts, dims=(2,))
    extent=maximum([cminmax[i][2]-cminmax[i][1] for i=1:2])
    zcoord=fill(extent/100,size(pts,2))

    parameter!(p,"points",vec(vcat(pts,zcoord')))
    parameter!(p,"polys",vtkpolys(tris))
    parameter!(p,"gridscale",args[:gridscale])
    parameter!(p,"aspect",args[:aspect])

    if markers!=nothing
        nregions=maximum(markers)
        if colormap==nothing
            colormap=GridVisualize.region_cmap(nregions)
        end
        if typeof(colormap)==Symbol
            cmap=colorschemes[colormap]
        else
            cmap=ColorScheme(colormap)
        end
        rgb=reinterpret(Float64,get(cmap,markers,(1,size(cmap))))
        parameter!(p,"colors",UInt8.(floor.(rgb*255)))

        bar_stops=collect(1:size(cmap))
        bar_rgb=reinterpret(Float64,get(cmap,bar_stops,(1,size(cmap))))
        bar_rgb=UInt8.(floor.(bar_rgb*255))
        p.jsdict["cbar"]=2
        p.jsdict["cbar_stops"]=bar_stops
        p.jsdict["cbar_colors"]=bar_rgb
        p.jsdict["cbar_levels"]=collect(1:size(cmap))
        p.jsdict["cbar_fontsize"]=args[:legendfontsize]

    else
        parameter!(p,"colors","none")
    end

    if edges!=nothing
        nedges=size(edges,2)
        lines=Vector{UInt32}(undef,3*nedges)
        iline=0
        for i=1:nedges
            lines[iline+1]=2
            lines[iline+2]=edges[1,i]-1  #  0-1 discrepancy between jl and js...
            lines[iline+3]=edges[2,i]-1
            iline=iline+3
        end
        parameter!(p,"lines",lines)

        if edgemarkers!=nothing
            (fmin,nbregions)=Int64.(extrema(edgemarkers))
            if edgecolormap==nothing
                edgecolormap=GridVisualize.bregion_cmap(nbregions)
            end
            if typeof(edgecolormap)==Symbol
                ecmap=colorschemes[edgecolormap]
            else
                ecmap=ColorScheme(edgecolormap)
            end
            edgergb=reinterpret(Float64,get(ecmap,edgemarkers,(1,size(ecmap))))
            parameter!(p,"linecolors",UInt8.(floor.(edgergb*255)))

            ebar_stops=collect(1:size(ecmap))
            ebar_rgb=reinterpret(Float64,get(ecmap,ebar_stops,(1,size(ecmap))))
            ebar_rgb=UInt8.(floor.(ebar_rgb*255))
            p.jsdict["ecbar_stops"]=ebar_stops
            p.jsdict["ecbar_colors"]=ebar_rgb
            p.jsdict["ecbar_levels"]=collect(1:size(ecmap))
        else
            parameter!(p,"linecolors","none")
        end
    else
        parameter!(p,"lines","none")
        parameter!(p,"linecolors","none")
    end

    axis2d!(p; args...)
    p
end

"""
$(SIGNATURES)

Plot parts of tetrahedral mesh below the planes given by the `*plane` parameters.
"""
function tetmesh!(p::PlutoVTKPlot, pts, tets;kwargs...)

    default_args=(markers=nothing,
                  colormap=nothing,
                  faces=nothing,
                  facemarkers=nothing,
                  facecolormap=nothing,
                  xplanes=[prevfloat(Inf)],
                  yplanes=[prevfloat(Inf)],
                  zplanes=[prevfloat(Inf)],
                  outlinealpha=0.1)

    args=merge(p.args,default_args)
    args=merge(args,kwargs)

    markers=args[:markers]
    colormap=args[:colormap]
    faces=args[:faces]
    facemarkers=args[:facemarkers]
    facecolormap=args[:facecolormap]

    xplane=args[:xplanes][1]
    yplane=args[:yplanes][1]
    zplane=args[:zplanes][1]

    ntet=size(tets,2)
    command!(p,"tetmesh")
    nregions=  markers==nothing  ? 0 : maximum(markers)
    nbregions= facemarkers==nothing ? 0 :  maximum(facemarkers)

    if nregions==0
        nregions=1
        markers=fill(Int32(1),ntet)
    end

    if faces!=nothing && nbregions==0
        nbregions=1
        facemarkers=fill(Int32(1),size(faces,2))
    end

    if colormap==nothing
        colormap=GridVisualize.region_cmap(nregions)
    end

    if facecolormap==nothing
        facecolormap=GridVisualize.bregion_cmap(nbregions)
    end

    xyzmin=zeros(3)
    xyzmax=ones(3)

    @views for idim=1:3
        xyzmin[idim]=minimum(pts[idim,:])
        xyzmax[idim]=maximum(pts[idim,:])
    end

    xyzcut=[xplane,yplane,zplane]

    regpoints0,regfacets0=GridVisualize.extract_visible_cells3D(pts,tets,markers,nregions,
                                                                xyzcut,
                                                                primepoints=hcat(xyzmin,xyzmax)
                                                                )

    points=hcat( [          reshape(reinterpret(Float32,regpoints0[i]),(3,length(regpoints0[i])))   for i=1:nregions]... )
    facets=vcat( [vtkpolys( reshape(reinterpret(Int32,  regfacets0[i]),(3,length(regfacets0[i]))),
                            offset= ( i==1 ? 0 : sum(k->length(regpoints0[k]),1:i-1) ) )
                            for i=1:nregions]... )

    regmarkers=vcat([fill(i,length(regfacets0[i])) for i=1:nregions]...)

    if typeof(colormap)==Symbol
        cmap=colorschemes[colormap]
    else
        cmap=ColorScheme(colormap)
    end
    rgb=reinterpret(Float64,get(cmap,regmarkers,(1,size(cmap))))
    nfaces=length(rgb)÷3

    bar_stops=collect(1:size(cmap))
    bar_rgb=reinterpret(Float64,get(cmap,bar_stops,(1,size(cmap))))
    bar_rgb=UInt8.(floor.(bar_rgb*255))
    p.jsdict["cbar"]=2
    p.jsdict["cbar_stops"]=bar_stops
    p.jsdict["cbar_colors"]=bar_rgb
    p.jsdict["cbar_levels"]=collect(1:size(cmap))
    p.jsdict["cbar_fontsize"]=args[:legendfontsize]

    if faces!=nothing
        bregpoints0,bregfacets0=GridVisualize.extract_visible_bfaces3D(pts,faces,facemarkers,nbregions,
                                                                       xyzcut,
                                                                       primepoints=hcat(xyzmin,xyzmax)
                                                                       )
        bregpoints=hcat([reshape(reinterpret(Float32,bregpoints0[i]),(3,length(bregpoints0[i]))) for i=1:nbregions]...)
        bregfacets=vcat([vtkpolys(reshape(reinterpret(Int32,bregfacets0[i]),(3,length(bregfacets0[i]))),
                                  offset= size(points,2) + ( i==1 ? 0 : sum(k->length(bregpoints0[k]),1:i-1) ) )
                                  for i=1:nbregions]...)
        bfacemarkers=vcat([fill(i,length(bregfacets0[i])) for i=1:nbregions]...)

        if typeof(facecolormap)==Symbol
            facecmap=colorschemes[facecolormap]
        else
            facecmap=ColorScheme(facecolormap)
        end
        facergb=reinterpret(Float64,get(facecmap,bfacemarkers,(1,size(facecmap))))
        facets=vcat(facets,bregfacets)
        points=hcat(points,bregpoints)
        rgb=vcat(rgb,facergb)

        ecmap=facecmap
        ebar_stops=collect(1:size(ecmap))
        ebar_rgb=reinterpret(Float64,get(ecmap,ebar_stops,(1,size(ecmap))))
        ebar_rgb=UInt8.(floor.(ebar_rgb*255))
        p.jsdict["ecbar_stops"]=ebar_stops
        p.jsdict["ecbar_colors"]=ebar_rgb
        p.jsdict["ecbar_levels"]=collect(1:size(ecmap))
    end

    if args[:outlinealpha]>0 && faces!=nothing
        parameter!(p,"outline",1)
        outline!(p,pts,faces,facemarkers,facecolormap,nbregions,xyzmin,xyzmax;alpha=args[:outlinealpha])
    else
        parameter!(p,"outline",0)
    end

    parameter!(p,"polys",facets)
    parameter!(p,"points",vec(points))
    parameter!(p,"colors",UInt8.(floor.(rgb*255)))

    axis3d!(p; args...)
    p
end

"""
$(SIGNATURES)

2D quiver.
"""
function quiver2d!(p::PlutoVTKPlot, pts, qvec; kwargs...)

    args=merge((kwargs...),p.args)
    command!(p,"quiver")
    zcoord=zeros(size(pts,2))

    cminmax=extrema(pts, dims=(2,))

    # length scale for arrowhead drawing
    extent=maximum([cminmax[i][2]-cminmax[i][1] for i=1:2])
    l0=extent/150

    pts3=vcat(pts,zcoord')
    vec3=vcat(qvec,zcoord')

    nvec=size(qvec,2)
    # Calculate points for arrowheads
    apts=copy(pts3)
    bpts=copy(pts3)

    # possibly make them parameters to satisfy user's taste
    arrowtip_width_factor=1.5
    arrowtip_length_factor=4.0

    for i=1:nvec

        #normal to vectors
        #     vnorm=sqrt(qvec[1,i]^2+qvec[2,i]^2)
        # We don't normalise it as everything should be scaled with vector length
        normal1= arrowtip_width_factor*l0*qvec[2,i]
        normal2=-arrowtip_width_factor*l0*qvec[1,i]

        # distance from tip
        dist1=arrowtip_length_factor*l0*qvec[1,i]
        dist2=arrowtip_length_factor*l0*qvec[2,i]

        # side points for arrowhead
        apts[1,i]=pts[1,i]+qvec[1,i]-dist1+normal1
        apts[2,i]=pts[2,i]+qvec[2,i]-dist2+normal2

        bpts[1,i]=pts[1,i]+qvec[1,i]-dist1-normal1
        bpts[2,i]=pts[2,i]+qvec[2,i]-dist2-normal2
    end

    qpts=hcat(pts3,pts3+vec3,apts,bpts)

    npts=size(pts,2)
    lines=Vector{UInt32}(undef,3*npts+4*npts)
    iline=0
    ipt=0
    for i=1:npts

        # vector line
        lines[iline+1]=2
        lines[iline+2]=ipt      # pts3
        lines[iline+3]=npts+ipt # pts3+vec3
        iline=iline+3

        # arrowhead
        lines[iline+1]=3
        lines[iline+2]=2*npts+ipt  # apts
        lines[iline+3]=npts+ipt    # pts3+vec3
        lines[iline+4]=3*npts+ipt  # bpts

        iline=iline+4

        ipt=ipt+1
    end

    parameter!(p,"points",vec(qpts))
    parameter!(p,"lines",lines)
    axis2d!(p; args...)
    p
end

#####################################
# Experimental part
"""
$(SIGNATURES)

Experimental: Plot piecewise linear function on  triangular grid given by points and triangles
as matrices
"""
function triplot!(p::PlutoVTKPlot,pts, tris,f; kwargs...)
    args=merge((kwargs...),p.args)
    p.jsdict=Dict{String,Any}("cmdcount" => 0)
    command!(p,"triplot")
    # make 3D points from 2D points by adding function value as
    # z coordinate
    p.jsdict["cbar"]=0
    parameter!(p,"points",vec(vcat(pts,f')))
    parameter!(p,"polys",vtkpolys(tris))
    axis3d!(p;args...)
    p
end

function plot!(p::PlutoVTKPlot,x,y; kwargs...)
    command!(p,"plot")
    n=length(x)
    points=vec(vcat(x',y',zeros(n)'))
    lines=collect(UInt16,0:n)
    lines[1]=n
    parameter!(p,"points",points)
    parameter!(p,"lines",lines)
    parameter!(p,"cam","2D")
    p
end


text:
struct MeshEntity{dim}
  index::Int
end
const MeshVertex = MeshEntity{0}
const MeshEdge = MeshEntity{1}
const MeshFace = MeshEntity{2}

#M = dim + 1
struct MeshConectivity
  indices::NTuple
  offsets::NTuple
end

#dim1 = dim+1 (better solutions?)
struct PolytopalMesh2{dim,T,dim1} <: AbstractPolytopalMesh{dim,T}
  entities::NTuple{dim1, Int}
  vertices::Vector{Tensors.Vec{dim, T}} #   NTuple{dim, T}}
  geometry::Dict{NTuple{2,Int},MeshConectivity}
  # Sets
  entitysets::Dict{Int,Dict{String,Set{NTuple{2,Int}}}}
end

#Generic Interface
getfacet(mesh::PolytopalMesh2{dim}, facet) where {dim} = getcellsubentities(mesh,facet.cellidx,dim-1)[facet.idx]
getnfacetvertices(mesh::PolytopalMesh2{dim}, cell_idx, facet_lidx) where {dim} = getnvertices(mesh,getcellsubentities(mesh,cell_idx,dim-1)[facet_lidx])
reference_facet_vertices(mesh::PolytopalMesh2{dim},cell::MeshEntity{dim})  where {dim} = reference_facet_vertices(getCellType(mesh,cell.index))

"Return coordinates keeping orientation"
function getverticescoords(mesh::PolytopalMesh2{d}, edge::EdgeIndex) where {d}
  return [mesh.vertices[i] for i in _local_edgesD(mesh,MeshEntity{d}(edge.cellidx))[edge.idx]]
end

function getverticescoords(mesh::PolytopalMesh2{d}, cell::Int) where {d}
  return getverticescoords(mesh, MeshEntity{d}(cell))
end

getdim(mesh::PolytopalMesh2{dim}) where {dim} = dim

function get_cell_connectivity_list(mesh::PolytopalMesh2{dim}) where {dim}
  _unpack_connectivity(get_connectivity!(mesh,dim,0))
end

getnvertices(mesh::PolytopalMesh2{dim}, entity::MeshEntity{dim}) where {dim} = getnsubentities(mesh,entity,0)
getncellvertices(mesh::PolytopalMesh2{dim}, cell_idx::Int) where {dim} = getnvertices(mesh,MeshEntity{dim}(cell_idx))
getnfacets(mesh::PolytopalMesh2{dim}, cell_idx) where {dim} = getncellsubentities(mesh,cell_idx,dim-1)

function gettopology(mesh::PolytopalMesh2{dim}, cell::MeshEntity{dim}) where {dim}
  topology = Dict{Int,Int}()
  for i in 0:(dim-1)
    push!(topology, i => getncellsubentities(mesh,cell.index,i))
  end
  topology[dim] = 1
  return topology
end

getnedges(mesh::PolytopalMesh2{dim}, cell::MeshEntity{dim}) where {dim} = getncellsubentities(mesh,cell.index,1)
getnfaces(mesh::PolytopalMesh2{dim}, cell::MeshEntity{dim}) where {dim} = getncellsubentities(mesh,cell.index,2)

function reference_edge_vertices(mesh::PolytopalMesh2{dim}, cell::MeshEntity{dim}) where {dim}
  N = getnvertices(mesh,cell)
  return Tuple((i,mod1(i+1,N)) for i in 1:N)
end

getverticesindices(mesh::PolytopalMesh2{dim}, cell::MeshEntity{dim}) where {dim} = getsubentitiesindexes(mesh,cell,0)

#Common API
@inline getnentities(mesh::PolytopalMesh2,d) = mesh.entities[d+1]
function getentities(mesh::PolytopalMesh2, d)
  return [MeshEntity{d}(idx) for idx in 1:getnentities(mesh,d)]
end
getentityset(mesh::PolytopalMesh2, entity::Int, set::String) = mesh.entitysets[entity][set]

getncells(mesh::PolytopalMesh2{dim}) where {dim} = getnentities(mesh,dim)
getnvertices(mesh::PolytopalMesh2) = getnentities(mesh,0)
getvertexset(mesh::PolytopalMesh2, set::String) = getentityset(mesh,0,set)
getedgeset(mesh::PolytopalMesh2, set::String) = getentityset(mesh,1,set)
getfacetset(mesh::PolytopalMesh2{dim}, set::String) where {dim} = getentityset(mesh,dim-1,set)
getcells(mesh::PolytopalMesh2{dim}) where {dim} = getentities(mesh,dim)
getcell(mesh::PolytopalMesh2{dim}, idx::Int) where {dim} = (1<= idx <= getncells(mesh) ? MeshEntity{dim}(idx) : error("cell index out of bounds"))
function getCellType(mesh::PolytopalMesh2{dim}, idx::Int) where {dim}
  return Cell{dim, getncellsubentities(mesh,idx,0), getncellsubentities(mesh,idx,1), getncellsubentities(mesh,idx,2)}
end
getvertices(mesh::PolytopalMesh2) = getentities(mesh,0)
getedges(mesh::PolytopalMesh2) = getentities(mesh,1)
getfacets(mesh::PolytopalMesh2{dim}) where {dim} = getentities(mesh,dim-1)

function getcellsubentities(mesh::PolytopalMesh2{dim},cellidx,entity::Int) where {dim}
  connectivity = get_connectivity!(mesh,dim,entity)
  return [MeshEntity{entity}(idx) for idx in get_entity_indices(connectivity,cellidx)]
end

function getsubentitiesindexes(mesh::PolytopalMesh2,entity::MeshEntity{d1},d2::Int) where {d1}
  connectivity = get_connectivity!(mesh,d1,d2)
  return get_entity_indices(connectivity,entity.index)
end

function getncellsubentities(mesh::PolytopalMesh2{dim},cellidx,entity::Int) where {dim}
  entity == dim && return 1
  connectivity = get_connectivity!(mesh,dim,entity)
  return connectivity.offsets[cellidx+1] - connectivity.offsets[cellidx]
end

function getnsubentities(mesh::PolytopalMesh2,entity::MeshEntity{d1},d2::Int) where {d1}
  connectivity = get_connectivity!(mesh,d1,d2)
  return connectivity.offsets[entity.index+1] - connectivity.offsets[entity.index]
end

getvertexcoords(mesh::PolytopalMesh2, vertex_idx::Int) = mesh.vertices[vertex_idx]

function getverticescoords(mesh::PolytopalMesh2, entity::MeshEntity)
  return [mesh.vertices[i] for i in getverticesidx(mesh, entity)]
end

function getvertexcoords(mesh::PolytopalMesh2, entity::MeshEntity, vidx::Int)
  return mesh.vertices[getverticesidx(mesh, entity)[vidx]]
end
function getnvertices(mesh::PolytopalMesh2,entity::MeshEntity{d}) where{d}
  connectivity = get_connectivity!(mesh,d,0)
  idx = entity.index
  return connectivity.offsets[idx+1] - connectivity.offsets[idx]
end

entityeltype(mesh::PolytopalMesh2, dim) = MeshEntity{dim}

#Internal
function get_entity_name(entity,dim)
  if entity == 0
    return "vertex"
  elseif entity > 0 && entity < dim
    if entity == 1
      return "edge"
    elseif entity == 2
      return "face"
    end
  elseif entity == dim
    return "cell"
  else
    return "dim $entity"
  end
end
@inline get_entity_indices(connectivity, idx) = connectivity.indices[connectivity.offsets[idx]:connectivity.offsets[idx+1]-1]
function getverticesidx(mesh::PolytopalMesh2{dim}, entity::MeshEntity{d}) where {dim,d}
  connectivity = get_connectivity!(mesh,d,0)
  return get_entity_indices(connectivity,entity.index)
end

function getentities(mesh::PolytopalMesh2, d2, d1, entity)
  connectivity = get_connectivity!(mesh,d2,d1)
  return [MeshEntity{d1}(i) for i in get_entity_indices(connectivity,entity.index)]
end

function _pack_connectivity(indices_mat)
  offsets = [1]
  k = 1
  indices = []
  for x in indices_mat
    k = k + size(x,1)
    push!(offsets,k)
    push!(indices,x...)
  end
  MeshConectivity(Tuple(indices), Tuple(offsets))
end

function _unpack_connectivity(connectivity::MeshConectivity)
  return [get_entity_indices(connectivity,i) for i in 1:(size(connectivity.offsets,1)-1)]
end

# Compute d1 -> d2 from d2 -> d1
function _transpose(mesh, d2, d1)
  indices = [Int[] for _ in 1:getnentities(mesh,d1)]
  for entity_j in getentities(mesh,d2)
    for entity_i in getentities(mesh,d2,d1,entity_j)
      push!(indices[entity_i.index], entity_j.index)
    end
  end
  _pack_connectivity(indices)
end

# Compute d1 -> d2 from d1 -> d3 and d3 -> d2
function _intersection(mesh,d1,d2,d3)
  indices = [Int[] for _ in 1:getnentities(mesh,d1)]
  for entity_i in getentities(mesh,d1)
    for entity_k in getentities(mesh,d1,d3,entity_i)
      for entity_j in getentities(mesh,d3,d2,entity_k)
        if (d1 == d2 && entity_i.index != entity_j.index) ||
          (d1 > d2 && all(x in getentities(mesh,d1,0,entity_i) for x in getentities(mesh,d2,0,entity_j)))
          if !(entity_j.index in indices[entity_i.index])
            push!(indices[entity_i.index], entity_j.index)
          end
        end
      end
    end
  end
  _pack_connectivity(indices)
end

function _local_edgesD(mesh::PolytopalMesh2{2},cell)
  indices = getverticesidx(mesh,cell)
  N = size(indices,1)
  idx = Tuple((i,mod1(i+1,N)) for i in 1:N)
  return [[indices[x[1]],indices[x[2]]] for x in idx]
end

#Compute D -> d and d -> 0 from D -> 0 and D -> D for 0 < d < D
function _build!(mesh::PolytopalMesh2{D},d) where {D}
  indices1 = [Int[] for _ in 1:getnentities(mesh,D)]
  indices2 = [Int[] for _ in 1:getnentities(mesh,d)]
  if d == 1
    edgesDict = Dict{Set{Int},Int}()
    nextidx = 1 #first edge index
    for cell in getentities(mesh,D)
      Vi = _local_edgesD(mesh, cell)
      for vi in Vi
        token = Base.ht_keyindex2!(edgesDict, Set(vi))
        if token > 0 # reuse edge index
          reuse_idx = edgesDict.vals[token]
          push!(indices1[cell.index], reuse_idx)    #sign is used for orientation
        else
          push!(indices1[cell.index], nextidx)
          indices2[nextidx] = vi
          Base._setindex!(edgesDict, nextidx, Set(vi), -token)
          nextidx += 1
        end
      end
    end
  #TODO if d == 2
  else
    error("Automatic construction of entities of dim $d failed")
  end
  V = _pack_connectivity(indices1)
  push!(mesh.geometry,(D,d)=>V)
  V = _pack_connectivity(indices2)
  push!(mesh.geometry,(d,0)=>V)
end

function get_connectivity!(mesh::PolytopalMesh2{D},d1::Int,d2::Int) where {D}
  if (0 < d1 < D) && !haskey(mesh.geometry,(d1,0))
      _build!(mesh,d1)
  elseif (d1 != d2) && (0 < d2 < D) && !haskey(mesh.geometry,(d2,0))
      _build!(mesh,d2)
  end
  if haskey(mesh.geometry,(d1,d2))
    get(mesh.geometry,(d1,d2),nothing)
  else
    if d1 < d2
      V = get_connectivity!(mesh,d2,d1)
      V = _transpose(mesh,d2,d1)
    else
      d3 = (d1 == 0 && d2 == 0) ? D : 0
      V = get_connectivity!(mesh,d1,d3)
      V = get_connectivity!(mesh,d3,d2)
      V = _intersection(mesh,d1,d2,d3)
    end
    push!(mesh.geometry,(d1,d2)=>V)
    return V
  end
end

ext: jl

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
import Base: iterate
import Base.Iterators.take

"""
Structure to hold the iterates of the lifted coordinate descent method.

LiftedCCIterable(M, b, A, λ; optTol=√eps)

Construct a lifted coordinate descent iterator for the problem

    minimize  ½‖Mx-b‖² + λ⋅gauge(A|x)

The absolute tolerance `optTol` is used to determine the stopping condition.

"""
mutable struct LiftedCCIterable{T1<:Float64, T2<:Vector{T1}, T3<:AbstractArray{T1}, T4<:AbstractLinearOp, T5<:AbstractAtomicSet}
    M::T4           # linear operator
    b::T2           # observation
    A::T5           # atomic set
    r::T2           # residual
    z::T3           # negative gradient
    Ma::T3          # image of the atom
    Mx::T3          # image of primal variable
    λ::T1           # current scale of atomic set
    optTol::T1      # tolerance
    function LiftedCCIterable(M::AbstractLinearOp, b::Vector{Float64}, A::AbstractAtomicSet, λ::Float64, optTol = sqrt(eps(1.0)))
        r = copy(b)
        z = M'*r
        Ma = M*expose(A, z)
        Mx = copy(Ma); fill!(Mx, 0.0)
        new{Float64, Vector{Float64}, AbstractArray{Float64}, AbstractLinearOp, AbstractAtomicSet}(M, b, A, r, z, Ma, Mx, λ, optTol)
    end
end

"""
Defines one complete iteration of the lifted coordinate descent method.
"""
function iterate(p::LiftedCCIterable, k::Int=0)

    # Compute  M*a  where  a = expose(A, z)
    s = support(p.A, p.z)
    a = expose(p.A, p.z)
    p.Ma .= p.M*a

    # Compute g = λ - support(A|z) and possible exit
    g = p.λ - s

    # Compute step size
    δ = -g / (norm(p.Ma)^2)

    # update iterates
    p.Mx .+= δ*p.Ma
    p.r .= p.b - p.Mx
    p.z .= p.M' * p.r

    return (g, p.r), k+1
end

"""
Return an exit flag for the lifted coordinate descent algorithm.
"""
function checkExitCd(g::Float64, u::Float64, k::Int64, α::Float64, feaTol::Float64, maxIts::Int64)
    u < α + feaTol && return :feasible                  # feasible
    g > -sqrt(eps(1.0))/2 && return :optimal                    # optimal
    k ≥ maxIts && return :iterations                    # out of iterations
    return :noerror                                     # otherwise
end

"""
Lifted coordinate descent method for the problem

    min ½‖Mx-b‖² + λ⋅gauge(A|x)
"""
function coordinate_descent(M::AbstractLinearOp, b::Vector{Float64}, A::AbstractAtomicSet, λ::Float64;
                   α::Float64 = 0.0,
                   tol::Float64 = 1e-12,
                   maxIts::Int = size(M,2),
                   logger::Bool = true
                   )

    feaTol = tol*(1+norm(b))
    # Structure to hold the iterates
    lcc = LiftedCCIterable(M, b, A, λ)
    # Structure to hold the solution
    sol = Solution(M, b, A)

    # head logging
    m, n = size(M)
    logger && logger_head_cd(m, n, feaTol, α, maxIts)

    # Check whether b is zero
    if norm(b) < tol
        return sol
    end

    # ----------------------------------------------------------------
    # min loop
    # ----------------------------------------------------------------
    k = 0
    exitFlag = :noerror
    for (g, r) in lcc

        exitFlag = primalrecover!(lcc, sol, α, feaTol, exitFlag)
        exitFlag == :noerror || break
        k = k + 1
        u = norm(r)^2/2
        exitFlag = checkExitCd(g, u, k, α, feaTol, maxIts)

        logger && logger_level_cd(α, k, u, g, exitFlag, sol.feas)
    end

    primalrecover!(lcc, sol, α, feaTol, exitFlag)

    # foot logging
    logger&& logger_foot_cd(b, sol.feas, k)

    return sol
end

function logger_head_cd(m, n, feaTol, α, maxIts)
    println()
    println("  -------------------------------------------------------------------------")
    println("  Lifted Coordinate Descent Method")
    println("  -------------------------------------------------------------------------")
    @printf("  %-22s %7d %7s %-20s %7d\n","number of variables",n,"","number of constraints",m)
    @printf("  %-22s %7.2e %7s %-20s %7.2e\n","feasibility tolerance",feaTol,"","α",α)
    @printf("  %-22s %7d \n","max iterations",maxIts)
    println("  -------------------------------------------------------------------------")
    @printf("  %8s   %8s   %8s   %8s   %10s\n",
            "iteration","u-α","g","infeas-α","exitFlag")
end

function logger_foot_cd(b, feas, totIts)
    scale = 1+norm(b)
    normr = sqrt(2*feas)
    println("  -------------------------------------------------------------------------")
    @printf("  %-22s %8.1e \n","residual (abs)",normr)
    @printf("  %-22s %8.1e \n","residual (rel)",normr/scale)
    @printf("  %-22s %8d \n","total iterations",totIts)
    println("  -------------------------------------------------------------------------")
end

function logger_level_cd(α, k, u, g, exitFlag, feas)
    @printf("  %8d   %8.2e   %8.2e   %8.2e   %s\n",
            k, u-α, g, feas-α, exitFlag)
end

"""
    primal recover
"""
function primalrecover!(p::LiftedCCIterable, sol::Solution, α::Float64, feaTol::Float64, exitFlag)
    flag = exitFlag
    s = support(p.A, p.z)
    F = face(p.A, p.z/s)
    ϵ = copy(p.r); ϵ .*= sqrt(2*α)/norm(ϵ)
    c, r = face_project(p.M, F, p.b - ϵ)
    feas = norm(r + ϵ)^2/2
    # @show feas
    if feas ≤ sol.feas
        sol.F = F
        sol.c = c
        sol.feas = feas
    end
    if feas - α ≤ feaTol
        flag = :feasible
    elseif flag == :optimal
        flag = :noerror
        replace_iterate!(p, sol)
    end
    return flag
end

"""
    replace iterates
"""
function replace_iterate!(p::LiftedCCIterable, sol::Solution)
    MF = p.M * sol.F
    p.Mx .= MF * sol.c
    p.r .= p.b - p.Mx
    p.z .= p.M' * p.r
end


text:
using Makie

scene = Scene();


text:
JopReshape(dom::JetAbstractSpace, rng::JetAbstractSpace) = JopLn(dom = dom, rng = rng, df! = JopReshape_df!)

export JopReshape

function JopReshape_df!(d, m; kwargs...)
    d[:] .= m[:]
    d
end


text:
module DifferentialEvolutionMCMC
    using Random, StatsBase, ProgressMeter, Parameters, Distributions
    using LinearAlgebra, ConcreteStructs
    import AbstractMCMC: AbstractSampler, AbstractModel
    import AbstractMCMC: bundle_samples, sample, MCMCThreads
    import MCMCChains: Chains
    import LinearAlgebra: norm
    export DE, Particle, DEModel, sample, MCMCThreads, fixed_gamma
    export variable_gamma, random_gamma, evaluate_fun!, compute_posterior!
    export greedy_update!, optimize, get_optimal, resample, as_union
    export DEMCMC

    include("structs.jl")
    include("main.jl")
    include("optimize.jl")
    include("migration.jl")
    include("crossover.jl")
    include("mutation.jl")
    include("utilities.jl")
end


text:
module ReadableRegex

export RegexString
export @rs_str

export WORD
export NON_WORD
export DIGIT
export NON_DIGIT
export WHITESPACE
export NON_WHITESPACE
export NON_LINEBREAK
export ANY
export BEGIN
export END
export WORDBOUND
export NON_WORDBOUND

export one_or_more
export at_least
export between
export exactly
export maybe
export zero_or_more
export look_for
export either
export capture
export reference
export char_in
export char_not_in
export @compile
export lazy_zero_or_more
export lazy_one_or_more
export lazy_between
export lazy_maybe
export lazy_at_least
export possessive_zero_or_more
export possessive_one_or_more
export possessive_between
# export possessive_maybe
export possessive_at_least

export LETTER
export UPPERCASE
export LOWERCASE
export TITLECASE
export MODIFIER_LETTER
export OTHER_LETTER
export MARK
export NONSPACING_MARK
export SPACING_COMBINING_MARK
export ENCLOSING_MARK
export NUMBER
export DEC_DIGIT_NUMBER
export LETTER_NUMBER
export OTHER_NUMBER
export SYMBOL
export MATH_SYMBOL
export CURRENCY
export MODIFIER_SYMBOL
export OTHER_SYMBOL
export PUNCTUATION
export CONNECTOR_PUNCT
export DASH_PUNCT
export OPEN_PUNCT
export CLOSE_PUNCT
export INITIAL_PUNCT
export FINAL_PUNCT
export OTHER_PUNCT
export SEPARATOR
export SPACE_SEP
export LINE_SEP
export PARAGRAPH_SEP
export OTHER
export CONTROL
export FORMAT
export SURROGATE
export PRIVATE_USE
export UNASSIGNED
export NON_LETTER
export NON_UPPERCASE
export NON_LOWERCASE
export NON_TITLECASE
export NON_MODIFIER_LETTER
export NON_OTHER_LETTER
export NON_MARK
export NON_NONSPACING_MARK
export NON_SPACING_COMBINING_MARK
export NON_ENCLOSING_MARK
export NON_NUMBER
export NON_DEC_DIGIT_NUMBER
export NON_LETTER_NUMBER
export NON_OTHER_NUMBER
export NON_SYMBOL
export NON_MATH_SYMBOL
export NON_CURRENCY
export NON_MODIFIER_SYMBOL
export NON_OTHER_SYMBOL
export NON_PUNCTUATION
export NON_CONNECTOR_PUNCT
export NON_DASH_PUNCT
export NON_OPEN_PUNCT
export NON_CLOSE_PUNCT
export NON_INITIAL_PUNCT
export NON_FINAL_PUNCT
export NON_OTHER_PUNCT
export NON_SEPARATOR
export NON_SPACE_SEP
export NON_LINE_SEP
export NON_PARAGRAPH_SEP
export NON_OTHER
export NON_CONTROL
export NON_FORMAT
export NON_SURROGATE
export NON_PRIVATE_USE
export NON_UNASSIGNED

"""
    RegexString(s::String)

Holds a string meant as the raw format of a regular expression.
This should always be convertible to a valid regex, so it's best
not to construct one yourself. Use the builder functions instead.
"""
struct RegexString
    s::String
end

macro rs_str(exp)
    :(RegexString($exp))
end

function noncapturing_group_or_append(inner::String, outer::String)
    if isempty(inner)
        error("Empty inner string")
    else
        noncapturing_group_or_token(inner) * outer
    end
end

# This helps keeping the nesting code simple,
# because there can be no interactions beyond the group barrier.
function noncapturing_group_or_token(s::String)
    if length(s) == 1
        s
    else
        "(?:" * s * ")"
    end
end

# Define helper constants.
const WORD = rs"\w"
const NON_WORD = rs"\W"
const DIGIT = rs"\d"
const NON_DIGIT = rs"\D"
const WHITESPACE = rs"\s"
const NON_WHITESPACE = rs"\S"
const NON_LINEBREAK = rs"."
const ANY = rs"[\s\S]"
const BEGIN = rs"^"
const END = rs"$"
const WORDBOUND = rs"\b"
const NON_WORDBOUND = rs"\B"

const LETTER = rs"\p{L}"
const UPPERCASE = rs"\p{Lu}"
const LOWERCASE = rs"\p{Ll}"
const TITLECASE = rs"\p{Lt}"
const MODIFIER_LETTER = rs"\p{Lm}"
const OTHER_LETTER = rs"\p{Lo}"
const MARK = rs"\p{M}"
const NONSPACING_MARK = rs"\p{Mn}"
const SPACING_COMBINING_MARK = rs"\p{Mc}"
const ENCLOSING_MARK = rs"\p{Me}"
const NUMBER = rs"\p{N}"
const DEC_DIGIT_NUMBER = rs"\p{Nd}"
const LETTER_NUMBER = rs"\p{Nl}"
const OTHER_NUMBER = rs"\p{No}"
const SYMBOL = rs"\p{S}"
const MATH_SYMBOL = rs"\p{Sm}"
const CURRENCY = rs"\p{Sc}"
const MODIFIER_SYMBOL = rs"\p{Sk}"
const OTHER_SYMBOL = rs"\p{So}"
const PUNCTUATION = rs"\p{P}"
const CONNECTOR_PUNCT = rs"\p{Pc}"
const DASH_PUNCT = rs"\p{Pd}"
const OPEN_PUNCT = rs"\p{Ps}"
const CLOSE_PUNCT = rs"\p{Pe}"
const INITIAL_PUNCT = rs"\p{Pi}"
const FINAL_PUNCT = rs"\p{Pf}"
const OTHER_PUNCT = rs"\p{Po}"
const SEPARATOR = rs"\p{Z}"
const SPACE_SEP = rs"\p{Zs}"
const LINE_SEP = rs"\p{Zl}"
const PARAGRAPH_SEP = rs"\p{Zp}"
const OTHER = rs"\p{C}"
const CONTROL = rs"\p{Cc}"
const FORMAT = rs"\p{Cf}"
const SURROGATE = rs"\p{Cs}"
const PRIVATE_USE = rs"\p{Co}"
const UNASSIGNED = rs"\p{Cn}"

const NON_LETTER = rs"\P{L}"
const NON_UPPERCASE = rs"\P{Lu}"
const NON_LOWERCASE = rs"\P{Ll}"
const NON_TITLECASE = rs"\P{Lt}"
const NON_MODIFIER_LETTER = rs"\P{Lm}"
const NON_OTHER_LETTER = rs"\P{Lo}"
const NON_MARK = rs"\P{M}"
const NON_NONSPACING_MARK = rs"\P{Mn}"
const NON_SPACING_COMBINING_MARK = rs"\P{Mc}"
const NON_ENCLOSING_MARK = rs"\P{Me}"
const NON_NUMBER = rs"\P{N}"
const NON_DEC_DIGIT_NUMBER = rs"\P{Nd}"
const NON_LETTER_NUMBER = rs"\P{Nl}"
const NON_OTHER_NUMBER = rs"\P{No}"
const NON_SYMBOL = rs"\P{S}"
const NON_MATH_SYMBOL = rs"\P{Sm}"
const NON_CURRENCY = rs"\P{Sc}"
const NON_MODIFIER_SYMBOL = rs"\P{Sk}"
const NON_OTHER_SYMBOL = rs"\P{So}"
const NON_PUNCTUATION = rs"\P{P}"
const NON_CONNECTOR_PUNCT = rs"\P{Pc}"
const NON_DASH_PUNCT = rs"\P{Pd}"
const NON_OPEN_PUNCT = rs"\P{Ps}"
const NON_CLOSE_PUNCT = rs"\P{Pe}"
const NON_INITIAL_PUNCT = rs"\P{Pi}"
const NON_FINAL_PUNCT = rs"\P{Pf}"
const NON_OTHER_PUNCT = rs"\P{Po}"
const NON_SEPARATOR = rs"\P{Z}"
const NON_SPACE_SEP = rs"\P{Zs}"
const NON_LINE_SEP = rs"\P{Zl}"
const NON_PARAGRAPH_SEP = rs"\P{Zp}"
const NON_OTHER = rs"\P{C}"
const NON_CONTROL = rs"\P{Cc}"
const NON_FORMAT = rs"\P{Cf}"
const NON_SURROGATE = rs"\P{Cs}"
const NON_PRIVATE_USE = rs"\P{Co}"
const NON_UNASSIGNED = rs"\P{Cn}"

# Convert a string with special regex char_in to one where they are all escaped with backslashes.
escaped(s::String) = replace(s, r"([\\\.\+\-\^\$\\|\?\*[\]\(\)\{\}])" => s"\\\1")

# Define handy conversion routines for RegexStrings.
Base.convert(::Type{RegexString}, s::String) = RegexString(escaped(s))
Base.convert(::Type{RegexString}, rs::RegexString) = rs
Base.convert(::Type{RegexString}, c::Char) = RegexString(escaped(string(c)))
Base.convert(::Type{RegexString}, sr::StepRange{Char, Int}) = RegexString("[$(escaped(string(sr.start)))-$(escaped(string(sr.stop)))]")
Base.convert(::Type{RegexString}, list::Union{AbstractVector, Tuple}) = either(list...)

_c(x) = convert(RegexString, x)

# These are functions that give the typical regex logic building blocks.
one_or_more(r) = RegexString(noncapturing_group_or_append(_c(r).s, "+"))
at_least(n::Int, r) = RegexString(noncapturing_group_or_append(_c(r).s, "{$n,}"))
between(low::Int, high::Int, r) = RegexString(noncapturing_group_or_append(_c(r).s, "{$low,$high}"))
exactly(n::Int, r) = RegexString(noncapturing_group_or_append(_c(r).s, "{$n}"))

maybe(r) = RegexString(noncapturing_group_or_append(_c(r).s, "?"))
zero_or_more(r) = RegexString(noncapturing_group_or_append(_c(r).s, "*"))
followed_by(r, by) = RegexString(noncapturing_group_or_token(_c(r).s) * "(?=$(_c(by).s))")
not_followed_by(r, by) = RegexString(noncapturing_group_or_token(_c(r).s) * "(?!$(_c(by).s))")
preceded_by(r, by) = RegexString("(?<=$(_c(by).s))" * noncapturing_group_or_token(_c(r).s))
not_preceded_by(r, by) = RegexString("(?<!$(_c(by).s))" * noncapturing_group_or_token(_c(r).s))

lazy_zero_or_more(r) = zero_or_more(r) * rs"?"
lazy_one_or_more(r) = one_or_more(r) * rs"?"
lazy_between(low, high, r) = between(low, high, r) * rs"?"
lazy_maybe(r) = maybe(r) * rs"?"
lazy_at_least(n, r) = at_least(n, r) * rs"?"

possessive_zero_or_more(r) = zero_or_more(r) * rs"+"
possessive_one_or_more(r) = one_or_more(r) * rs"+"
possessive_between(low, high, r) = between(low, high, r) * rs"+"
# possessive_maybe(r) = maybe(r) * rs"+"
possessive_at_least(n, r) = at_least(n, r) * rs"+"

function look_for(r;
        after = nothing,
        before = nothing,
        not_after = nothing,
        not_before = nothing)

    current = r
    if !isnothing(after)
        current = preceded_by(current, after)
    end
    if !isnothing(before)
        current = followed_by(current, before)
    end
    if !isnothing(not_after)
        current = not_preceded_by(current, not_after)
    end
    if !isnothing(not_before)
        current = not_followed_by(current, not_before)
    end
    current
end

either(args...) = RegexString(noncapturing_group_or_token(join([noncapturing_group_or_token(_c(r).s) for r in args], "|")))

function capture(r; as = nothing)
    if isnothing(as)
        RegexString("($(_c(r).s))")
    else
        RegexString("(?<$as>$(_c(r).s))")
    end
end

reference(i::Int) = RegexString("\\$i")
reference(name) = RegexString("\\k<$name>")

char_in(args...) =     RegexString("""[$(join((_charset(arg) for arg in args), ""))]""")
char_not_in(args...) = RegexString("""[^$(join((_charset(arg) for arg in args), ""))]""")

# allowed characters are different within a set, only \ - ^ and ] have to be escaped
_escape_for_set(s::String) = replace(s, r"([\\\-\^\]])" => s"\\\1")

_charset(s::String) = _escape_for_set(s)
_charset(c::Char) = c
_charset(sr::StepRange{Char, Int}) = "$(_escape_for_set(string(sr.start)))-$(_escape_for_set(string(sr.stop)))"
_charset(rs::RegexString) = rs.s

# Define the multiplication operator on RegexStrings as concatenation like Strings.
# Anything can be concatenated if it can be converted to a RegexString.
Base.:*(r1::RegexString, r2::RegexString) = RegexString(r1.s * r2.s)
Base.:*(anything, r::RegexString) = convert(RegexString, anything) * r
Base.:*(r::RegexString, anything) = r * convert(RegexString, anything)

Base.Regex(r::RegexString) = Regex(r.s)
Base.Regex(r::RegexString, flags) = Regex(r.s, flags)
Base.match(rs::RegexString, x) = Base.match(Regex(rs), x)
Base.eachmatch(rs::RegexString, x) = Base.eachmatch(Regex(rs), x)

macro compile(exp)
    regexstring = @eval $exp
    Regex(regexstring.s)
end

end # module

max_stars_repo_name: jkrumbiegel/ReadableRegex

max_issues_repo_name: jkrumbiegel/ReadableRegex

max_forks_repo_name: jkrumbiegel/ReadableRegex


text:
export
    isvalidtree, fulltree, transformtree

function isvalidtree(signal, tree::BitVector)
    maxlevels = levelmax(signal)
    treesize = length(tree)
    treesize == (1 << maxlevels) - 1 || return false
    maxleaf = ((1 << (maxlevels - 1)) - 1)
    @assert maxleaf << 1 + 1 <= treesize
    leaves = firstindex(tree):maxleaf
    return all(leaves) do leaf
        return @inbounds tree[leaf] || !(tree[leaf << 1] || tree[leaf << 1 + 1])
    end
end

function maketree(signal, level)
    maxlevels = levelmax(signal)
    treesize = (1 << maxlevels) - 1
    maxleaf = ((1 << (maxlevels - 1)) - 1)
    @assert 0 <= level <= maxlevels
    @assert maxleaf << 1 + 1 <= treesize
    tree = falses(treesize)
    return tree
end

function fulltree(signal, level)
    tree = maketree(signal, level)
    return filltree!(tree, level)
end

function filltree!(tree, level)
    @inbounds for leaf in 1:((1 << level) - 1)
        tree[leaf] = true
    end
    return tree
end

function transformtree(signal, level)
    tree = maketree(signal, level)
    return transformtree!(tree, level)
end

function transformtree!(tree, level)
    @inbounds for eachlevel in 1:level
        tree[(1 << eachlevel) - 1] = true
    end
    return tree
end


text:
import Base: string, print, show, isless, hash, ==, read

if VERSION < v"0.7-"

    """
        SHA1(x::Vector{UInt8})
        SHA1(str::AbstractString)

    A 20-byte sha1 hash.
    """
    struct SHA1
        bytes::Vector{UInt8}
        function SHA1(bytes::Vector{UInt8})
            length(bytes) == 20 ||
                throw(ArgumentError("wrong number of bytes for SHA1 hash: $(length(bytes))"))
            return new(bytes)
        end
    end
    SHA1(s::AbstractString) = SHA1(hex2bytes(s))

    string(hash::SHA1) = bytes2hex(hash.bytes)
    print(io::IO, hash::SHA1) = print(io, string(hash))

    show(io::IO, hash::SHA1) = print(io, "SHA1(\"", string(hash), "\")")
    isless(a::SHA1, b::SHA1) = lexless(a.bytes, b.bytes)
    hash(a::SHA1, h::UInt) = hash((SHA1, a.bytes), h)
    ==(a::SHA1, b::SHA1) = a.bytes == b.bytes

end

macro sha1_str(str)
    SHA1(str)
end

function Base.read(io::IO, ::Type{SHA1})
    data = Array{UInt8}(undef, 20)
    SHA1(read!(io, data))
end

max_stars_repo_name: simonbyrne/GitX

max_issues_repo_name: simonbyrne/GitX

max_forks_repo_name: simonbyrne/GitX


text:
import Base.Iterators: product
####
#### Output Structure
####

"""
    struct RNNOutput
        y
        hidden
        memory
        indices
    end

RNN layer outputs are always `RNNOutput`
`hidden`,`memory` and `indices` may be nothing depending on the keyword arguments you used in forward.

`y` is last hidden states of each layer. `size(y)=(H/2H,[B,T])`.
If you use unequal length instances in a batch input, `y` becomes 2D array `size(y)=(H/2H,sum_of_sequence_lengths)`.
See `indices` below to get correct time outputs for a specific instance or
see  `PadRNNOutput` to pad whole output.
`h` is the hidden states in each timesstep. `size(h) = h,B,L/2L`

`c` is the hidden states in each timesstep. `size(h) = h,B,L/2L`

`indices` is corresponding instace indices for your `RNNOutput.y`. You may call `yi = y[:,indices[i]]`.
"""
struct RNNOutput{T,V,Z,S<:VVecOrNothing}
    y::T
    hidden::V
    memory::Z
    indices::S
end

function Base.show(io::IO,R::RNNOutput{T,V,Z,S}) where {T,V,Z,S}
    print(io,"RNNOutput(y: $T$(size(R.y))")
    V !== Nothing && print(io,", hidden : $V$(size(R.hidden))")
    Z !== Nothing && print(io,", memory : $Z$(size(R.memory))")
    S !== Nothing && print(io,", indices: $S$(size(R.indices))")
    print(io,")")
end

####
#### RNN Types
####
"""
    SRNN(;input=inputSize, hidden=hiddenSize, activation=:relu, options...)
    LSTM(;input=inputSize, hidden=hiddenSize, options...)
    GRU(;input=inputSize, hidden=hiddenSize, options...)

    (1) (l::T)(x; kwargs...) where T<:AbstractRNN
    (2) (l::T)(x::Array{Int}; batchSizes=nothing, kwargs...) where T<:AbstractRNN
    (3) (l::T)(x::Vector{Vector{Int}}; sorted=false, kwargs...) where T<:AbstractRNN

All RNN layers has above forward run(1,2,3) functionalities.

(1) `x` is an AbstractArray{T,N} with size equals d,[B,T]

(2) For this, an RNN with embedding layer is needed. see `embed` option below.
`x` is an integer array and inputs corresponds one-hot indices.
One can give 2D array for minibatching as rows corresponds to one instance, or
One can give 1D array with minibatching by specifying batch batchSizes argument.
see `Knet.rnnforw` for the details of batchSizes.

(3) For this, RNN with embedding layer is needed.
`x` is an vector of integer vectors. Every integer vector corresponds to an
instance. It automatically batches inputs for the user. It is better to give inputs as sorted.
If your inputs sorted you can make `sorted` argument true to increase performance.

see `RNNOutput` for output details

# options

* `embed=nothing`: embedding size or and embedding layer
* `numLayers=1`: Number of RNN layers.
* `bidirectional=false`: Create a bidirectional RNN if `true`.
* `dropout=0`: Dropout probability. Ignored if `numLayers==1`.
* `skipInput=false`: Do not multiply the input with a matrix if `true`.
* `algo=0`: Algorithm to use, see CUDNN docs for details.
* `seed=0`: Random number seed for dropout. Uses `time()` if 0.
* `winit=xavier`: Weight initialization method for matrices.
* `binit=zeros`: Weight initialization method for bias vectors.
* `atype=KnetLayers.arrtype`: array type for model weights.

# Keywords

* hx=nothing : initial hidden states
* cx=nothing : initial memory cells
* hy=false   : if true returns h
* cy=false   : if true returns c

"""
AbstractRNN

for layer in (:SRNN, :LSTM, :GRU)
    layername=string(layer)
    @eval begin

        mutable struct $layer{P,E<:LayerOrNothing} <: AbstractRNN{P,E}
            embedding::E
            params::P
            specs::RNN
            gatesview::DictOrNothing
        end

        @inline (m::$layer)(x,h...;o...) = RNNOutput(_forw(m,x,h...;o...)...)
        # FIXME: activation input is not compatible with rest of the package
        function $layer(;input::Integer, hidden::Integer, embed=nothing, activation=:relu, atype=arrtype, o...)
            embedding,inputSize = _getEmbed(input,embed)
            rnnType = $layer==SRNN ? activation : Symbol(lowercase($layername))
            r = RNN(inputSize, hidden; rnnType, atype, o...)
            $layer(embedding,r.w,r,gatesview($layer,r))
        end
    end
end

###
### Internal Utils
###
function Base.show(io::IO, m::AbstractRNN{P,E}) where {P,E}
    r = m.specs
    embedSize,inputSize = m.embedding !== nothing ? size(m.embedding.weight) : (nothing,r.inputSize)
    print(io, ("SRNN(ReLU){","SRNN(Tanh){","LSTM{","GRU{")[r.mode+1], "$P, $E}(input=", inputSize, ",hidden=", r.hiddenSize)
    if embedSize!=nothing; print(io,",embed=",embedSize); end
    if r.direction == 1; print(io, ",bidirectional"); end
    if r.numLayers > 1; print(io, ",layers=", r.numLayers); end
    if r.dropout > 0; print(io, ",dropout=", r.dropout); end
    if r.inputMode == 1; print(io, ",skipinput"); end
    if r.dataType != Float32; print(io, ',', r.dataType); end
    print(io, ')')
end

@inline _getEmbed(input::Integer,embed::Nothing) = (nothing,input)
@inline _getEmbed(input::Integer,embed::Embed) =
    size(embed.weight,2) == input ? (embed,size(embed.weight,1)) : error("dimension mismatch in embeddings")
@inline _getEmbed(input::Integer,embed::Integer) = (Embed(input=input,output=embed),embed)

gate_mappings(::Type{<:SRNN}) = Dict(:h=>(1,2))
gate_mappings(::Type{<:GRU})  = Dict(:r=>(1,4),:u=>(2,5),:n=>(3,6))
gate_mappings(::Type{<:LSTM}) = Dict(:i=>(1,5),:f=>(2,6),:n=>(3,7),:o=>(4,8))
const input_mappings = Dict(:i=>1,:h=>2)
const param_mappings = Dict(:w=>1,:b=>2)

gatesview(T::Type{<:AbstractRNN},r::RNN) =
    Dict((Symbol(ty,ih,g,l,d),rnnparam(r, (r.direction+1)*(l-1)+d+1, id[ihid], param; useview=true))
          for (l,d,(g,id),(ih,ihid),(ty,param)) in
              product(1:r.numLayers,0:r.direction,gate_mappings(T),
                      input_mappings,param_mappings))

# Saves from unnecessary memory taken by gatesview
function _ser(r::T, s::IdDict, m::typeof(JLDMODE)) where T <: AbstractRNN
    if !haskey(s,r)
        if r.gatesview !== nothing
            s[r] = T(_ser(r.embedding,s,m), _ser(r.params,s,m), _ser(r.specs,s,m), nothing)
        else
            s[r] = T(_ser(r.embedding,s,m), _ser(r.params,s,m), _ser(r.specs,s,m), gatesview(T,s[r.specs]))
        end
    end
    return s[r]
end

####
#### User Utils
####
"""

    PadSequenceArray(batch::Vector{Vector{T}}; direction=:Right, pad=0) where T<:Integer

Pads a batch of integer arrays with zeros

```julia
julia> PadSequenceArray([[1,2,3],[1,2],[1]])
3×3 Array{Int64,2}:
 1  2  3
 1  2  0
 1  0  0

 julia> PadSequenceArray([[1,2,3],[1,2],[1]], direction=:Left)
 3×3 Array{Int64,2}:
  1  2  3
  0  1  2
  0  0  1
```

"""
function PadSequenceArray(batch::Vector{Vector{T}}; direction=:Right, pad=0) where T<:Integer
    B      = length(batch)
    lngths = length.(batch)
    Tmax   = maximum(lngths)
    padded = Array{T}(undef,B,Tmax)
    @inbounds for n = 1:B
        if direction == :Right
            padded[n,1:lngths[n]] = batch[n]
            padded[n,lngths[n]+1:end] .= pad
        else direction == :Left
            padded[n,1:end-lngths[n]] .= pad
            padded[n,end-lngths[n]+1:end] = batch[n]
        end
    end
    return padded
end

# FIXME: long
"""
    PadRNNOutput(s::RNNOutput)
Pads a RNNOutput if it is produced by unequal length batches
`size(s.y)=(H/2H,sum_of_sequence_lengths)` becomes `(H/2H,B,Tmax)`
"""
@inline PadRNNOutput(s::RNNOutput{<:Any,<:Any,<:Any,Nothing}) = s,nothing

function PadRNNOutput(s::RNNOutput)
    d = size(s.y,1)
    B = length(s.indices)
    lngths = length.(s.indices)
    Tmax = maximum(lngths)
    mask = trues(d,B,Tmax)
    cw = []
    @inbounds for i=1:B
        y1 = s.y[:,s.indices[i]]
        df = Tmax-lngths[i]
        if df > 0
            mask[:,:,end-df+1:end] .= false
            pad = fill!(arrtype(undef,d*df),zero(eltype(arrtype)))
            ypad = reshape(cat1d(y1,pad),d,Tmax) # hcat(y1,kpad)
            push!(cw,ypad)
        else
            push!(cw,y1)
        end
    end
    RNNOutput(reshape(vcat(cw...),d,B,Tmax),s.hidden,s.memory,nothing),mask
end

"""
    _pack_sequence(batch::Vector{Vector{T}}) where T<:Integer

Packs unequal length of sequence batches for cuDNN format.
It return a tuple of tokens and batchSizes.
`tokens` consists of all input tokens in the order of cuDNN format. So, length(tokens) == sum(length,batch)
`batchSizes` keeps the information of how many instance available for a time sequence. So, lengh(bathSizes) = maximum(length,batch)

#Example
```julia
julia> _pack_sequence([[1,2,3,4],[7,6],[8,5],[9,10]])
(tokens = [1, 7, 8, 9, 2, 6, 5, 10, 3, 4], batchSizes = [4, 4, 1, 1])
```
"""
function _pack_sequence(batch::Vector{Vector{T}}) where T<:Integer
    B      = length(batch)
    Lmax   = length(first(batch))
    tokens = zeros(Int,sum(length,batch))
    bsizes = zeros(Int,Lmax)
    i = 1
    @inbounds for t = 1:Lmax
        bs = 0
        @inbounds for k = 1:B
            if t<=length(batch[k])
                tokens[i] = batch[k][t]
                bs += 1
                i  += 1
            end
        end
        bsizes[t] = bs
    end
    return (tokens=tokens,batchSizes=bsizes)
end

"""
    _batchSizes2indices(batchSizes::Vector{<:Integer})

    Finds the indices of the `tokens` belongs to for each instance.
    see `_pack_sequence` for details
"""
@inline _batchSizes2indices(batchSizes) =
    map(1:batchSizes[1]) do i
        @inbounds [i;i.+cumsum(filter(x->(x>=i),batchSizes)[1:end-1])]
    end

# FIXME: long
function _forw(rnn::AbstractRNN{<:Any,<:Embed},batch::Vector{Vector{T}},h...;
               sorted=issorted(batch,by=length),o...) where T<:Integer

    if all(length.(batch).==length(batch[1]))
        return _forw(rnn,permutedims(cat(batch...;dims=2),(2,1)),h...;o...)
    end

    if !sorted
        v     = sortperm(batch; by=length, rev=true, alg=MergeSort)
        batch = batch[v]
    end

    tokens, bsizes = _pack_sequence(batch)
    inds = _batchSizes2indices(bsizes)
    y,hidden,memory,_ = _forw(rnn,tokens,h...;batchSizes=bsizes,o...)

    if !sorted
        rev  = invperm(v)
        return y,_sort3D(hidden,rev),_sort3D(memory,rev),inds[rev]
    end
    return y,hidden,memory,inds
end

@inline function _forw(specs::RNN,params,x,h...;o...)
    y,hidden,memory,_ = rnnforw(specs,params,x,h...;o...)
    return y,hidden,memory,nothing
end

@inline _forw(rnn::AbstractRNN{<:Any,Nothing},x,h...;o...) =
    _forw(rnn.specs,rnn.params,x,h...;o...)

@inline _forw(rnn::AbstractRNN,x,h...;o...) =
    _forw(rnn.specs, rnn.params, rnn.embedding(x),h...;o...)

@inline function _forw(rnn::AbstractRNN{<:Any,<:Embed}, seq::Vector{T}, h...; batchSizes=nothing, o...) where T<:Integer
     inp = batchSizes===nothing ? reshape(seq,1,length(seq)) : seq
    _forw(rnn.specs, rnn.params, rnn.embedding(inp), h...; batchSizes=batchSizes,o...)
end

@inline _forw(rnn::AbstractRNN{<:Any,<:Nothing}, seq::Array{T},h...;o...) where T<:Integer =
    error("Integer inputs can only be used with RNNs that has embeddings.")

@inline _sort3D(hidden::Nothing,inds) = nothing
@inline _sort3D(hidden::AbstractArray,inds) = hidden[:,inds,:]
@inline _sort3D(h::KnetArray,inds)  =
    reshape(cat1d(map(i->h[:,:,i][:,inds],1:size(h,3))...),size(h))


text:
using DataFrames
using CSV
using Plots
#%%
data_k1 = open("wave1.dat")
df_k1=CSV.read(data_k1, header=false, delim=" ")
r_k1=df_k1[1]
phi_k1=df_k1[2]

data_k2 = open("wave2.dat")
df_k2=CSV.read(data_k2, header=false, delim=" ")
r_k2=df_k2[1]
phi_k2=df_k2[2]

data_k3 = open("wave3.dat")
df_k3=CSV.read(data_k3, header=false, delim=" ")
r_k3=df_k3[1]
phi_k3=df_k3[2]

data_k4 = open("wave4.dat")
df_k4=CSV.read(data_k4, header=false, delim=" ")
r_k4=df_k4[1]
phi_k4=df_k4[2]

#%%
plot(df_k1[1], df_k1[2], label="k 1", dpi=600)
plot!(df_k2[1], df_k2[2], label="k 2", dpi=600)
plot!(df_k3[1], df_k3[2], label="k 3", dpi=600)
plot!(df_k4[1], df_k4[2], label="k 4", dpi=600)
plot!(grid=true, legend=:topleft)
xlabel!("r")
ylabel!("phi")
savefig("wave.png")
# xaxis!(:log10)
# yaxis!(:log10)

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
# An explicitly unrolled LSTM with fixed sequence length.
using MXNet

#--LSTMState
immutable LSTMState
  c :: mx.SymbolicNode
  h :: mx.SymbolicNode
end
#--/LSTMState

#--LSTMParam
immutable LSTMParam
  i2h_W :: mx.SymbolicNode
  h2h_W :: mx.SymbolicNode
  i2h_b :: mx.SymbolicNode
  h2h_b :: mx.SymbolicNode
end
#--/LSTMParam

#--lstm_cell
function lstm_cell(data::mx.SymbolicNode, prev_state::LSTMState, param::LSTMParam;
                   num_hidden::Int=512, dropout::Real=0, name::Symbol=gensym())

  if dropout > 0
    data = mx.Dropout(data, p=dropout)
  end

  i2h = mx.FullyConnected(data=data, weight=param.i2h_W, bias=param.i2h_b,
                          num_hidden=4num_hidden, name=symbol(name, "_i2h"))
  h2h = mx.FullyConnected(data=prev_state.h, weight=param.h2h_W, bias=param.h2h_b,
                          num_hidden=4num_hidden, name=symbol(name, "_h2h"))

  gates = mx.SliceChannel(i2h + h2h, num_outputs=4, name=symbol(name, "_gates"))

  in_gate     = mx.Activation(gates[1], act_type=:sigmoid)
  in_trans    = mx.Activation(gates[2], act_type=:tanh)
  forget_gate = mx.Activation(gates[3], act_type=:sigmoid)
  out_gate    = mx.Activation(gates[4], act_type=:sigmoid)

  next_c = (forget_gate .* prev_state.c) + (in_gate .* in_trans)
  next_h = out_gate .* mx.Activation(next_c, act_type=:tanh)

  return LSTMState(next_c, next_h)
end
#--/lstm_cell

#--LSTM-part1
function LSTM(n_layer::Int, seq_len::Int, dim_hidden::Int, dim_embed::Int, n_class::Int;
              dropout::Real=0, name::Symbol=gensym(), output_states::Bool=false)

  # placeholder nodes for all parameters
  embed_W = mx.Variable(symbol(name, "_embed_weight"))
  pred_W  = mx.Variable(symbol(name, "_pred_weight"))
  pred_b  = mx.Variable(symbol(name, "_pred_bias"))

  layer_param_states = map(1:n_layer) do i
    param = LSTMParam(mx.Variable(symbol(name, "_l$(i)_i2h_weight")),
                      mx.Variable(symbol(name, "_l$(i)_h2h_weight")),
                      mx.Variable(symbol(name, "_l$(i)_i2h_bias")),
                      mx.Variable(symbol(name, "_l$(i)_h2h_bias")))
    state = LSTMState(mx.Variable(symbol(name, "_l$(i)_init_c")),
                      mx.Variable(symbol(name, "_l$(i)_init_h")))
    (param, state)
  end
  #...
  #--/LSTM-part1

  #--LSTM-part2
  # now unroll over time
  outputs = mx.SymbolicNode[]
  for t = 1:seq_len
    data   = mx.Variable(symbol(name, "_data_$t"))
    label  = mx.Variable(symbol(name, "_label_$t"))
    hidden = mx.FullyConnected(data=data, weight=embed_W, num_hidden=dim_embed,
                               no_bias=true, name=symbol(name, "_embed_$t"))

    # stack LSTM cells
    for i = 1:n_layer
      l_param, l_state = layer_param_states[i]
      dp = i == 1 ? 0 : dropout # don't do dropout for data
      next_state = lstm_cell(hidden, l_state, l_param, num_hidden=dim_hidden, dropout=dp,
                             name=symbol(name, "_lstm_$t"))
      hidden = next_state.h
      layer_param_states[i] = (l_param, next_state)
    end

    # prediction / decoder
    if dropout > 0
      hidden = mx.Dropout(hidden, p=dropout)
    end
    pred = mx.FullyConnected(data=hidden, weight=pred_W, bias=pred_b, num_hidden=n_class,
                             name=symbol(name, "_pred_$t"))
    smax = mx.SoftmaxOutput(pred, label, name=symbol(name, "_softmax_$t"))
    push!(outputs, smax)
  end
  #...
  #--/LSTM-part2

  #--LSTM-part3
  # append block-gradient nodes to the final states
  for i = 1:n_layer
    l_param, l_state = layer_param_states[i]
    final_state = LSTMState(mx.BlockGrad(l_state.c, name=symbol(name, "_l$(i)_last_c")),
                            mx.BlockGrad(l_state.h, name=symbol(name, "_l$(i)_last_h")))
    layer_param_states[i] = (l_param, final_state)
  end

  # now group all outputs together
  if output_states
    outputs = outputs ∪ [x[2].c for x in layer_param_states] ∪
                        [x[2].h for x in layer_param_states]
  end
  return mx.Group(outputs...)
end
#--/LSTM-part3

# Negative Log-likelihood
type NLL <: mx.AbstractEvalMetric
  nll_sum  :: Float64
  n_sample :: Int

  NLL() = new(0.0, 0)
end

function mx.update!(metric :: NLL, labels :: Vector{mx.NDArray}, preds :: Vector{mx.NDArray})
  @assert length(labels) == length(preds)
  nll = 0.0
  for (label, pred) in zip(labels, preds)
    @mx.nd_as_jl ro=(label, pred) begin
      nll -= sum(log(max(broadcast_getindex(pred, round(Int,label+1), 1:length(label)), 1e-20)))
    end
  end

  nll = nll / length(labels)
  metric.nll_sum += nll
  metric.n_sample += length(labels[1])
end

function mx.get(metric :: NLL)
  nll  = metric.nll_sum / metric.n_sample
  perp = exp(nll)
  return [(:NLL, nll), (:perplexity, perp)]
end

function mx.reset!(metric :: NLL)
  metric.nll_sum  = 0.0
  metric.n_sample = 0
end


text:
# ---
# title: 739. Daily Temperatures
# id: problem739
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Stack
# link: <https://leetcode.com/problems/daily-temperatures/description/>
# hidden: true
# ---
#
# Given a list of daily temperatures `T`, return a list such that, for each day
# in the input, tells you how many days you would have to wait until a warmer
# temperature. If there is no future day for which this is possible, put `0`
# instead.
#
# For example, given the list of temperatures `T = [73, 74, 75, 71, 69, 72, 76,
# 73]`, your output should be `[1, 1, 4, 2, 1, 1, 0, 0]`.
#
# **Note:** The length of `temperatures` will be in the range `[1, 30000]`. Each
# temperature will be an integer in the range `[30, 100]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function _precompile()
    precompile(ParseState, (String,))
    precompile(ParseState, (String, Int))
    precompile(next, (ParseState,))
    precompile(parse_doc, (ParseState,))
    precompile(Expr, (EXPR,))
    precompile(parse_string_or_cmd, (ParseState, Bool))

    # precompile(setbinding!, (EXPR,))
    # precompile(setbinding!, (EXPR, EXPR))
    # precompile(Binding, (EXPR,))
    # precompile(Binding, (String, EXPR, Nothing, Vector{Any}, Nothing))
    # precompile(Scope, (Nothing, Dict{String,Binding}, Nothing, Bool))
    # precompile(Scope, (Scope, Dict{String,Binding}, Nothing, Bool))
    # precompile(Scope, (Nothing, Dict{String,Binding}, Dict{String,Any}, Bool))
    # precompile(Scope, (Scope, Dict{String,Binding}, Dict{String,Any}, Bool))
    # precompile(Scope, ())
    precompile(EXPR, (Head, Nothing, Int, Int, Nothing, Tokens.Kind, Bool, Nothing, Nothing, Nothing, Nothing))
    precompile(EXPR, (Head, Vector{EXPR}, Int, Int, Nothing, Tokens.Kind, Bool, Nothing, Nothing, Nothing, Nothing))
    precompile(EXPR, (Head, Vector{EXPR}, Int, Int, String, Tokens.Kind, Bool, Nothing, Nothing, Nothing, Nothing))
    precompile(EXPR, (Head, Vector{EXPR}, Int, Int))
    precompile(EXPR, (Head, Vector{EXPR}))

    precompile(mIDENTIFIER, (ParseState,))
    precompile(mPUNCTUATION, (ParseState,))
    precompile(mPUNCTUATION, (Tokens.Kind, Int, Int))
    precompile(mKEYWORD, (ParseState,))
    precompile(mKEYWORD, (Tokens.Kind, Int, Int))
    precompile(mOPERATOR, (ParseState,))
    precompile(mOPERATOR, (Tokens.Kind, Int, Int, Bool))
    precompile(mLITERAL, (ParseState,))
    precompile(mLITERAL, (Int, Int, String, Tokens.Kind))
    precompile(INSTANCE, (ParseState,))

    precompile(tostr, (IOBuffer,))
    precompile(str_value, (EXPR,))

    precompile(CSTParser.parse_expression, (ParseState,))
    precompile(CSTParser.parse_paren, (ParseState,))
    precompile(CSTParser._continue_doc_parse, (ParseState, EXPR))
    precompile(CSTParser.parse_compound, (ParseState, EXPR))
    precompile(CSTParser.parse_operator, (ParseState, EXPR, EXPR))
    precompile(parse_unary, (ParseState, EXPR))
    precompile(parse_unary_colon, (ParseState, EXPR))
    precompile(parse_operator_anon_func, (ParseState, EXPR, EXPR))
    precompile(parse_operator_colon, (ParseState, EXPR, EXPR))
    precompile(parse_operator_cond, (ParseState, EXPR, EXPR))
    precompile(parse_operator_dot, (ParseState, EXPR, EXPR))
    precompile(parse_operator_eq, (ParseState, EXPR, EXPR))
    precompile(parse_operator_power, (ParseState, EXPR, EXPR))
    precompile(parse_operator_where, (ParseState, EXPR, EXPR))

    precompile(CSTParser.parse_block, (ParseState, Vector{EXPR}))
    precompile(CSTParser.parse_block, (ParseState, Vector{EXPR}, Tuple{Tokens.Kind}))
    precompile(CSTParser.parse_block, (ParseState, Vector{EXPR}, Tuple{Tokens.Kind,Bool}))
    precompile(CSTParser.parse_iterator, (ParseState,))
    precompile(CSTParser.parse_iterators, (ParseState, Bool))
    precompile(CSTParser.parse_call, (ParseState, EXPR, Bool))
    precompile(CSTParser.parse_comma_sep, (ParseState, Vector{EXPR}, Bool, Bool, Bool))
    precompile(CSTParser.parse_parameters, (ParseState, Vector{EXPR}, Vector{EXPR}))
    precompile(CSTParser.parse_macrocall, (ParseState,))
    precompile(CSTParser.parse_generator, (ParseState, EXPR))
    precompile(CSTParser.parse_dot_mod, (ParseState, Bool))

    precompile(CSTParser.parse_tuple, (ParseState, EXPR))
    precompile(CSTParser.parse_array, (ParseState, Bool))
    precompile(CSTParser.parse_ref, (ParseState, EXPR))
    precompile(CSTParser.parse_curly, (ParseState, EXPR))
    precompile(CSTParser.parse_braces, (ParseState,))

    precompile(CSTParser.parse_kw, (ParseState,))
    precompile(CSTParser.parse_const, (ParseState,))
    precompile(CSTParser.parse_global, (ParseState,))
    precompile(CSTParser.parse_local, (ParseState,))
    precompile(CSTParser.parse_return, (ParseState,))
    precompile(CSTParser.parse_abstract, (ParseState,))
    precompile(CSTParser.parse_primitive, (ParseState,))
    precompile(CSTParser.parse_imports, (ParseState,))
    precompile(CSTParser.parse_export, (ParseState,))
    precompile(CSTParser.parse_blockexpr, (ParseState, CSTParser.Head))
    precompile(CSTParser.parse_if, (ParseState, Bool))
    precompile(CSTParser.parse_try, (ParseState,))
    precompile(CSTParser.parse_do, (ParseState, EXPR))
    precompile(CSTParser.parse_mutable, (ParseState,))

    precompile(CSTParser.closer, (ParseState,))
    precompile(CSTParser.create_tmp, (Closer,))
    precompile(CSTParser.update_from_tmp!, (Closer, Closer_TMP))
    precompile(CSTParser.update_to_default!, (Closer,))

    precompile(Tokenize.Lexers.Lexer, (IOBuffer, Type{Tokens.RawToken}))
    precompile(Tokenize.Lexers.next_token, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(get, (Dict{Char,Tokenize.Tokens.Kind}, Char, Tokenize.Tokens.Kind))
    precompile(getindex, (Dict{Tokenize.Tokens.Kind,Symbol}, Tokenize.Tokens.Kind))

    precompile(Tokenize.Lexers.next_token, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_amper, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_backslash, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_bar, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_circumflex, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_cmd, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_colon, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_comment, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_digit, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken}, Tokens.Kind))
    precompile(Tokenize.Lexers.lex_division, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_dollar, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_dot, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_equal, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_exclaim, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_forwardslash, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_greater, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_identifier, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken}, Char))
    precompile(Tokenize.Lexers.lex_less, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_minus, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_percent, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_plus, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_prime, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_quote, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_star, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_whitespace, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
    precompile(Tokenize.Lexers.lex_xor, (Tokenize.Lexers.Lexer{IOBuffer,Tokens.RawToken},))
end


text:
using CuArrays
using Flux
using Flux: @functor, Recur, LSTMCell, dropout, mse, logitbinarycrossentropy
using Zygote
using Zygote: Buffer
using OMEinsum
using NamedTupleTools
using Statistics

CuArrays.allowscalar(false)

include("utils.jl")
include("dataprepLJSpeech/dataprep.jl")
include("model.jl")

###
# datadir = "/Users/aza/Projects/TTS/data/LJSpeech-1.1"
datadir = "/home/azamat/Projects/TTS/LJSpeech-1.1"
metadatapath = joinpath(datadir, "metadata.csv")
melspectrogramspath = joinpath(datadir, "melspectrograms.jld2")

batchsize = 11
batches_trn, alphabet = build_batches(metadatapath, melspectrogramspath, batchsize)

batch = batches_trn[argmin(map(x -> size(last(x), 2), batches_trn))]
textindices, meltarget, stoptarget = batch
meltarget, stoptarget = meltarget |> gpu, stoptarget |> gpu
time_out = size(stoptarget, 2)

m = Tacotron₂(alphabet)

Flux.reset!(m, textindices)

##
m(textindices, time_out)

##
@time loss(m, textindices, meltarget, stoptarget)

θ = Flux.params(m)
# #=check=# length(θ) == 71 == 1 + 3(2 + 2 + 0) + 10 + ((2 + 4) + 4 + 10 + 2 + 4 + 0) + 2 + 5(2 + 2 + 0)

@time gs = gradient(θ) do
   loss(m, textindices, meltarget, stoptarget)
end

Juno.@profiler gradient(θ) do
   loss(m, batch)
end


text:
# Basic operations
const OP1 = 1 # >
const OP2 = 2 # <
const OP3 = 3 # +
const OP4 = 4 # -
const OP5 = 5 # ,
const OP6 = 6 # .
const OP7 = 7 # [
const OP8 = 8 # ]

function char2op(chr::Char)
    const chars2ops = Dict(['>', '<', '+', '-', ',', '.', '[', ']'], 1:8)
    return chars2ops[chr]
end

function lexandparse(s::String)
    ops = Array(Int, 0)

    i = nextind(s, 0)

    line, character = 1, 0

    while i <= endof(s)
        chr = s[i]
        character += 1

        while isspace(chr)
            if chr == '\n'
                line += 1
                character = 0
            end
            i = nextind(s, i)
            if i <= endof(s)
                chr = s[i]
                character += 1
            else
                return ops
            end
        end

        if !(chr in ['>', '<', '+', '-', ',', '.', '[', ']'])
            msg = @sprintf(
                "Invalid token found at line %d, character %d",
                line,
                character
            )
            throw(ArgumentError(msg))
        else
            push!(ops, char2op(chr))
        end

        i = nextind(s, i)
    end

    return ops
end


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

Clip the range of grid values. will set values < low to below and/or values > high to above.
You can also specify one or more intervals where all values should be set to ``between``,
or replace individual values.

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

Parameters
----------

- **cmd0** :: [Type => Str]

    Either the input file name or the full monolitic options string. Do not use this
    when the grid (a GMTgrid type) is passed via the ``arg1`` argument.
- **G** | **save** | **outgrid** | **outfile** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdclip(....) form.
    ($(GMTdoc)grdclip.html#g)
- $(GMT.opt_R)
- **above** | **high** :: [Type => Array | Str]

    Two elements array with ``high`` and ``above`` or a string with "high/above".
    It sets all data[i] > ``high`` to ``above``.
- **below** | **low** :: [Type => Array | Str]

    Two elements array with ``low`` and ``below`` or a string with "low/below".
    It sets all data[i] < ``low`` to ``below``.
- **between** :: [Type => Array | Str]

    Three elements array with ``low, high`` and ``between`` or a string with "low/high/between".
    It sets all data[i] >= ``low`` and <= ``high`` to ``between``.
- **old** | **new** :: [Type => Array | Str]

    Two elements array with ``old`` and ``new`` or a string with "old/new".
    It sets all data[i] == ``old`` to ``new``.
- **S** :: [Type => Str]

    Condense all replacement options above in a single string.
    ($(GMTdoc)grdclip.html#s)
- $(GMT.opt_V)

-   Examples:

        G=gmt("grdmath", "-R0/10/0/10 -I1 X");
        G2=grdclip(G, above=[5 6], low=[2 2], between="3/4/3.5")

    or (note the use of -S for second on options because we can't repeat a kwarg name)

        G2=grdclip(G, S="a5/6 -Sb2/2 -Si3/4/3.5")
"""
function grdclip(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("grdclip", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd, = parse_common_opts(d, "", [:G :R :V_params])
	cmd  = add_opt(d, cmd, "S", [:S])
	cmd  = opt_S(d, cmd, [:high :above], 'a')
	cmd  = opt_S(d, cmd, [:low :below], 'b')
	cmd  = opt_S(d, cmd, [:old :new], 'r')
	cmd  = opt_S(d, cmd, [:between], 'i')

	common_grd(d, cmd0, cmd, "grdclip ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdclip(arg1, cmd0::String=""; kw...) = grdclip(cmd0, arg1; kw...)

# ---------------------------------------------------------------------------------------------------
function opt_S(d::Dict, cmd::String, symbs, flag::Char)
	# This is common to the 4 cases
	val, symb = find_in_dict(d, symbs)
	if (val !== nothing)
		cmd = string(cmd, " -S", flag)
		if (isa(val, String))
			cmd *= val
		elseif (isa(val, Array) || isa(val, Tuple))
			if (symb == :between && length(val) == 3)
				cmd = @sprintf("%s%.16g/%.16g/%.16g", cmd, val[1], val[2], val[3])
			elseif (length(val) == 2)
				cmd = @sprintf("%s%.16g/%.16g", cmd, val[1], val[2])
			else
				error("Wrong number of elements in S option")
			end
		else
			error("OPT_S: argument must be a string or a two elements array.")
		end
	end
	return cmd
end


text:
using Base: String

struct Style{T} end

Style(s::Symbol) = Style{s}()

"""
    polish_names!(table::CleanTable; style::Symbol=:snake_case)

Return a `CleanTable` where column names have been replaced to be unique and formated using
the `style` selected.

# Styles
- snake_case
- camelCase
"""
function polish_names!(table::CleanTable; style::Symbol=:snake_case)
    rename!(table, generate_polished_names(names(table); style=style))

    return table
end

"""
    polish_names(table; style=:snake_case)

Create and return a `CleanTable` with copied columns having column names replaced to be unique and formated
using the `style` selected.

# Styles
- snake_case
- camelCase
"""
function polish_names(table; style::Symbol=:snake_case)
    return polish_names!(CleanTable(table); style=style)
end

"""
    generate_polished_names(names; style::Symbol=:snake_case)

Return a vector of symbols containing new names that are unique and formated using the `style` selected.
"""
function generate_polished_names(names; style::Symbol=:snake_case)
    return generate_polished_names(names, Style(style))
end

function generate_polished_names(names, ::Style{:snake_case})
    new_names = Vector{String}()

    for name in names
        new_name = _sanitize_snake_case(
            join(split(_replace_uppers(String(name)), r"[\s\-.]"; keepempty=false), "_")
        )
        push!(new_names, new_name)
    end

    return _sanitize_dupes(new_names)
end

function generate_polished_names(names, ::Style{:camelCase})
    new_names = Vector{String}()

    for name in names
        new_name = lowercasefirst(
            join(uppercasefirst.(split(String(name), r"[\s\-._]"; keepempty=false)), "")
        )
        push!(new_names, new_name)
    end

    return _sanitize_dupes(new_names)
end

function generate_polished_names(names, ::Style)
    return error("Invalid style selected. Options are :snake_case, :camelCase")
end

function _replace_uppers(word)
    fixed_word = ""

    for letter in word
        if isuppercase(letter)
            fixed_word = fixed_word * "_" * lowercase(letter)
        else
            fixed_word = fixed_word * letter
        end
    end

    return fixed_word
end

function _sanitize_snake_case(dirty_snake)
    new_name = ""
    under_score_behind = false

    for letter in strip(dirty_snake, '_')
        if letter != '_'
            new_name = new_name * letter
            under_score_behind = false
        elseif !under_score_behind
            new_name = new_name * letter
            under_score_behind = true
        end
    end

    return new_name
end

function _sanitize_dupes(names)
    new_names = Vector{Symbol}()
    dupes = Dict{String,Int64}()

    for name in names
        if !(Symbol(name) in new_names)
            push!(new_names, Symbol(name))
        else
            dupes["$name"] = get(dupes, name, 0) + 1
            push!(new_names, Symbol(name * "_" * string(get(dupes, name, 0))))
        end
    end

    return new_names
end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
using HTTPClient.HTTPC
using JSON
using Compat
using Base.Test

import Compat.String

const upload_file = joinpath(dirname(@__FILE__), "uploadfile.txt")

# get a request bin name
r=HTTPC.post("http://requestb.in/api/v1/bins", "")
@test r.http_code == 200

jdict = JSON.parse(String(r.body))
@test haskey(jdict, "name")

RB = "http://requestb.in/" * jdict["name"]
HB = "http://httpbin.org"

println("URL : $RB")

r=HTTPC.get(RB * "?test=1")
@test r.http_code == 200
println("Test 1 passed, http_code : " * string(r.http_code))

r=HTTPC.get(RB * "?test=1.hdrs", RequestOptions(headers=[("Foo", "Bar"),("Baz", "Wux")]))
@test r.http_code == 200
println("Test 1.hdrs passed, http_code : " * string(r.http_code))

ostream=open("/tmp/google.txt", "w+")
r=HTTPC.get("http://www.google.com/", RequestOptions(ostream=ostream))
@test r.http_code == 200
println("Test 1.ostream passed, http_code : " * string(r.http_code))
close(ostream)

r=HTTPC.get("http://www.google.com/", RequestOptions(ostream="/tmp/google2.txt"))
@test r.http_code == 200
println("Test 1.ostream_file passed, http_code : " * string(r.http_code))

rr=HTTPC.get(RB * "?test=1nb", RequestOptions(blocking=false))
r = fetch(rr)
@test r.http_code == 200
println("Test 1nb passed, http_code : " * string(r.http_code))

r=HTTPC.get(HB * "/cookies/set?k1=v1&k2=v2")
@test r.http_code == 200
@test Set(r.headers["Set-Cookie"]) == Set(["k1=v1; Path=/", "k2=v2; Path=/"])
println("Test.headers passed")

r=HTTPC.get(RB, RequestOptions(query_params=[(:test, 1.1), (:Hello, "\"World\""), ("_rt", "&!@#%")]))
@test r.http_code == 200
println("Test 1.1 passed, http_code : " * string(r.http_code))

rr=HTTPC.get(RB, RequestOptions(query_params=[(:test, 1.1), (:Hello, "\"World\""), ("_rt", "&!@#%")], blocking=false))
r = fetch(rr)
@test r.http_code == 200
println("Test 1.2 passed, http_code : " * string(r.http_code))

r=HTTPC.get(RB * "?test=2", RequestOptions(request_timeout=60.0))
@test r.http_code == 200
println("Test 2 passed, http_code : " * string(r.http_code))

@test_throws ErrorException HTTPC.get(RB * "?test=low_rto", RequestOptions(request_timeout=0.001))

r=HTTPC.put(RB * "?test=3", "some text")
@test r.http_code == 200
println("Test 3 passed, http_code : " * string(r.http_code))

rr=HTTPC.put(RB * "?test=3nb", "some text", RequestOptions(blocking=false))
r = fetch(rr)
@test r.http_code == 200
println("Test 3nb passed, http_code : " * string(r.http_code))

r=HTTPC.put(RB, "some text", RequestOptions(query_params=[(:test, 3.1), (:Hello, "\"World\"")]))
@test r.http_code == 200
println("Test 3.1 passed, http_code : " * string(r.http_code))

rr=HTTPC.put(RB, "some text", RequestOptions(query_params=[(:test, 3.2), (:Hello, "\"World\"")], blocking=false))
r = fetch(rr)
@test r.http_code == 200
println("Test 3.2 passed, http_code : " * string(r.http_code))

r=HTTPC.put(RB * "?test=4", (:file, upload_file))
@test r.http_code == 200
println("Test 4 passed, http_code : " * string(r.http_code))

r=HTTPC.post(RB * "?test=5", "[1,2]")
@test r.http_code == 200
println("Test 5 passed, http_code : " * string(r.http_code))

rr=HTTPC.post(RB * "?test=5nb", "[1,2]", RequestOptions(blocking=false))
r = fetch(rr)
@test r.http_code == 200
println("Test 5nb passed, http_code : " * string(r.http_code))

r=HTTPC.post(RB, "[1,2]", RequestOptions(query_params=[(:test, 5.1), (:Hello, "\"World\"")]))
@test r.http_code == 200
println("Test 5.1 passed, http_code : " * string(r.http_code))

r=HTTPC.post(RB * "?test=6", "[1,2]", RequestOptions(content_type="application/json"))
@test r.http_code == 200
println("Test 6 passed, http_code : " * string(r.http_code))

rr=HTTPC.post(RB * "?test=6.1", "[1,2]", RequestOptions(content_type="application/json", blocking=false))
r = fetch(rr)
@test r.http_code == 200
println("Test 6.1 passed, http_code : " * string(r.http_code))

r=HTTPC.post(RB * "?test=6.1.0", @compat Dict("a" => 1, "b" => 2))
@test r.http_code == 200
println("Test 6.1.0 passed, http_code : " * string(r.http_code))

r=HTTPC.post(RB * "?test=6.1.1", [("a",1), ("b",2)])
@test r.http_code == 200
println("Test 6.1.1 passed, http_code : " * string(r.http_code))

ios = IOBuffer()
write(ios, "Hello World!")
seekstart(ios)
r=HTTPC.post(RB * "?test=6.2", ios)
@test r.http_code == 200
println("Test 6.2 passed, http_code : " * string(r.http_code))

r=HTTPC.post(RB * "?test=7", (:file, upload_file), RequestOptions(content_type="text/plain"))
@test r.http_code == 200
println("Test 7 passed, http_code : " * string(r.http_code))

r=HTTPC.post(RB * "?test=8", (:file, upload_file))
@test r.http_code == 200
println("Test 8 passed, http_code : " * string(r.http_code))

rr=HTTPC.post(RB * "?test=8.1", (:file, upload_file), RequestOptions(blocking=false))
r = fetch(rr)
@test r.http_code == 200
println("Test 8.1 passed, http_code : " * string(r.http_code))

r=HTTPC.head(RB * "?test=9")
@test r.http_code == 200
println("Test 9 passed, http_code : " * string(r.http_code))

rr=HTTPC.head(RB * "?test=9nb", RequestOptions(blocking=false))
r = fetch(rr)
@test r.http_code == 200
println("Test 9nb passed, http_code : " * string(r.http_code))

r=HTTPC.trace(RB * "?test=10")
println("Test 10 passed, http_code : " * string(r.http_code))
#@test r.http_code == 200

rr=HTTPC.trace(RB * "?test=10nb", RequestOptions(blocking=false))
r = fetch(rr)
println("Test 10nb passed, http_code : " * string(r.http_code))

r=HTTPC.delete(RB * "?test=11")
println("Test 11 passed, http_code : " * string(r.http_code))
@test r.http_code == 200

rr=HTTPC.delete(RB * "?test=11nb", RequestOptions(blocking=false))
r = fetch(rr)
println("Test 11nb passed, http_code : " * string(r.http_code))
@test r.http_code == 200

r=HTTPC.options(RB * "?test=12")
println("Test 12 passed, http_code : " * string(r.http_code))
@test r.http_code == 200

rr=HTTPC.options(RB * "?test=12nb", RequestOptions(blocking=false))
r = fetch(rr)
println("Test 12nb passed, http_code : " * string(r.http_code))
@test r.http_code == 200

function waitnexec(id)
    tname = "async" * string(id)
    global trigger
    while (trigger != :go)
        sleep(0.001)
    end
    r=HTTPC.get(RB * "?test=" * tname)
    println("Test " * tname * " passed, http_code : " * string(r.http_code))
    r
end

# Run 100 requests in parallel asynchronously
trigger = :wait
rrefs = [remotecall(waitnexec, 1, i) for i in 1:100]

trigger = :go
# wait for all of them to finish
for ref in rrefs
    wait(ref)
end


text:
# # Testing the stability of CKSVAR models

# Here, we test the stability of the models estimated in:
# - Ikeda, Li, Mavroeidis and Zanetti (2021), "[Testing the effectiveness of monetary policy in Japan and the United States](https://arxiv.org/abs/2012.15158)." Working paper.
# This paper is hereafter referred to as [ILMZ21].
##
# # Japan
# In [ILMZ21], the effectiveness of unconventional monetary policy in Japan is investigated via estimation of a CKSVAR(2) model in inflation, output gap and the Call Rate for Japan, over the period 1985q3-2019q1, with shadow call rate as the latent variable.
#
# Converting this model into TAR form, we determine upper bounds on the state-constrained joint spectral radius (SCJSR)
##
using DataFrames, DataFramesMeta, CSV, HTTP
using ThresholdStability

extract_est(parname, df) = @subset(df, in([parname]).(:Parameter)).:Estimate[1]

JP = CSV.read(HTTP.get("https://raw.githubusercontent.com/samwycherley/ThresholdStability.jl/master/examples/src/estimates/estimates_JP.csv").body, DataFrame, copycols=true)  # Japan data from CSV

JPβtilde = zeros(2)
for i in 1:2
    JPβtilde[i] = extract_est("tildebeta_$i", JP)
end

JPCstar = zeros(3, 2)
for i in 1:3
    for j in 1:2
        JPCstar[i, j] = extract_est("Eq.$i lSR_$j", JP)
    end
end

Cbars = []
for j in 1:2
    Cbar_j = zeros(3, 3)
    for i in 1:3
        Cbar_j[i, 1] = extract_est("Eq.$i INFL_$j", JP)
        Cbar_j[i, 2] = extract_est("Eq.$i YGAP_BOJ_$j", JP)
        Cbar_j[i, 3] = extract_est("Eq.$i SR_$j", JP)
    end
    push!(Cbars, Cbar_j)
end
JPCbar = hcat(Cbars[1], Cbars[2])

JPC = copy(JPCbar); k = size(JPC, 1)  # convert C̄ to C
for i in 1:2
    JPC[:, i*k] -= JPCstar[:, i]
end

Σ_JP, X_JP = CKSVAR_to_TAR(JPC, JPCstar, JPβtilde, 2)

G_JP = automaton_constructor(Σ_JP)
s_JP = discreteswitchedsystem(Σ_JP, G_JP, X_JP)
##
@show sosbound_γ(s_JP, 2)
##
# # United States
# For the United States, [ILMZ21] estimate a CKSVAR(3) model in inflation, output gap, and the Federal Funds Rate in first differences. The shadow Federal Funds Rate is the latent variable, and also enters in first differences.
#
# Because the Federal Funds Rate enters in first differences, the SCJSR is no longer appropriate, since `Σ` is now a continuum of matrices. Instead we rely on the joint spectral radius.
##
USfd = CSV.read(HTTP.get("https://raw.githubusercontent.com/samwycherley/ThresholdStability.jl/master/examples/src/estimates/estimates_US-fd.csv").body, DataFrame, copycols=true)  # US data from CSV

USfdβtilde = zeros(2)
for i in 1:2
    USfdβtilde[i] = extract_est("tildebeta_$i", USfd)
end

Cstar = zeros(3, 3)
for i in 1:3
    for j in 1:3
        Cstar[i, j] = extract_est("Eq.$i lFEDFUNDS_$j", USfd)
    end
end
USfdFstar = zeros(3, 2)
for i in 1:2
    for j in 1:i
        USfdFstar[1:2, i] += Cstar[1:2, j]
        USfdFstar[3, i] += Cstar[3, j]
    end
    USfdFstar[1:2, i] -= USfdβtilde
    USfdFstar[3, i] -= 1.
end

Cbars = []
for j in 1:3
    Cbar_j = zeros(3, 3)
    for i in 1:3
        Cbar_j[i, 1] = extract_est("Eq.$i INFL_$j", USfd)
        Cbar_j[i, 2] = extract_est("Eq.$i YGAP_$j", USfd)
        Cbar_j[i, 3] = extract_est("Eq.$i FEDFUNDS_$j", USfd)
    end
    push!(Cbars, Cbar_j)
end
Cs = []
for j in 1:3
    C_j = copy(Cbars[j])
    C_j[:, 3] -= Cstar[:, j]
    push!(Cs, C_j)
end
Fs = []
for i in 1:3
    F_i = zeros(3, 3)
    F_i[1:3, 1:2] = Cs[i][1:3, 1:2]
    for j in 1:i
        F_i[1:2, 3] += Cs[j][1:2, 3]
    end
    F_i[1:2, 3] += USfdβtilde
    for j in 1:i
        F_i[3, 3] += Cs[j][3, 3]
    end
    push!(Fs, F_i)
end
USfdF = hcat(Fs[1], Fs[2], Fs[3])

Σ_US = CKSVAR_to_companionFD(USfdF, USfdFstar, USfdβtilde, 3, diff=true)
##
@show jsr(Σ_US)
##
# Both models are stable.