text:
include("flippablelayout.jl")

function initialize!(p::GridVisualizer,::Type{MakieType})

    Makie=p.context[:Plotter]

    # Check for version compatibility
    version_min=v"0.15"
    version_max=v"0.15.99"

    version_installed=PkgVersion.Version(Makie.Makie)

    if version_installed<version_min
        error("Outdated version $(version_installed) of Makie. Please upgrade to at least $(version_min)")
    end

    if version_installed>version_max
        @warn("Possibly breaking version $(version_installed) of Makie.")
    end

    # Prepare flippable layout
    FlippableLayout.setmakie!(Makie)

    layout=p.context[:layout]

    parent,flayout=FlippableLayout.flayoutscene(resolution=p.context[:resolution])

    p.context[:figure]=parent
    p.context[:flayout]=flayout

    # copy arguments to sublayout
    for I in CartesianIndices(layout)
        ctx=p.subplots[I]
        ctx[:figure]=parent
        ctx[:flayout]=flayout
    end

    # Don't call display on pluto
    if !isdefined(Main, :PlutoRunner)
        Makie.display(parent)
    end

    parent
end

# Adding a scene to the layout just adds to the
# flippable layout.
add_scene!(ctx,ax)=ctx[:flayout][ctx[:subplot]...]=ax

# Revealing the  visualizer just returns the figure
reveal(p::GridVisualizer,::Type{MakieType}) = p.context[:figure]

function reveal(ctx::SubVisualizer,TP::Type{MakieType})
    FlippableLayout.yieldwait(ctx[:flayout])
    if ctx[:show]||ctx[:reveal]
        reveal(ctx[:GridVisualizer],TP)
    end
end

save(fname,p::GridVisualizer,::Type{MakieType})= p.context[:Plotter].save(fname, p.context[:figure])
save(fname,scene,Makie,::Type{MakieType})=Makie.save(fname, scene)

"""

     scene_interaction(update_scene,view,switchkeys::Vector{Symbol}=[:nothing])

Control multiple scene elements via keyboard up/down keys.
Each switchkey is assumed to correspond to one of these elements.
Pressing a switch key transfers control to its associated element.

Control of values of the current associated element is performed
by triggering change values via up/down (± 1)  resp. page_up/page_down (±10) keys

The update_scene callback gets passed the change value and the symbol.
"""
function scene_interaction(update_scene,scene,Makie,switchkeys::Vector{Symbol}=[:nothing])

    # Check if pixel position pos sits within the scene
    function _inscene(scene,pos)
        area=scene.px_area[]
        pos[1]>area.origin[1] &&
            pos[1] < area.origin[1]+area.widths[1] &&
            pos[2]>area.origin[2] &&
            pos[2] < area.origin[2]+area.widths[2]
    end

    # Initial active switch key is the first in the vector passed
    activeswitch=Makie.Node(switchkeys[1])

    # Handle mouse position within scene
    mouseposition=Makie.Node((0,0))

    Makie.on(scene.events.mouseposition) do m
        mouseposition[]=m
        false
    end

    # Set keyboard event callback
    Makie.on(scene.events.keyboardbutton) do buttons
        if _inscene(scene,mouseposition[])
            # On pressing a switch key, pass control
            for i=1:length(switchkeys)
                if switchkeys[i]!=:nothing && Makie.ispressed(scene,getproperty(Makie.Keyboard,switchkeys[i]))
                    activeswitch[]=switchkeys[i]
                    update_scene(0,switchkeys[i])
                    return true
                end
            end

            # Handle change values via up/down control
            if Makie.ispressed(scene, Makie.Keyboard.up)
                update_scene(1,activeswitch[])
                return true
            elseif Makie.ispressed(scene, Makie.Keyboard.down)
                update_scene(-1,activeswitch[])
                return true
            elseif Makie.ispressed(scene, Makie.Keyboard.page_up)
                update_scene(10,activeswitch[])
                return true
            elseif Makie.ispressed(scene, Makie.Keyboard.page_down)
                update_scene(-10,activeswitch[])
                return true
            end
        end
        return false
    end
end

# Standard kwargs for Makie scenes
scenekwargs(ctx)=Dict(:xticklabelsize => 0.5*ctx[:fontsize],
                      :yticklabelsize => 0.5*ctx[:fontsize],
                      :zticklabelsize => 0.5*ctx[:fontsize],
                      :xlabelsize => 0.5*ctx[:fontsize],
                      :ylabelsize => 0.5*ctx[:fontsize],
                      :zlabelsize => 0.5*ctx[:fontsize],
                      :xlabeloffset => 20,
                      :ylabeloffset => 20,
                      :zlabeloffset => 20,
                      :titlesize => ctx[:fontsize])

############################################################################################################
#1D grid

# Point list for node markers
function basemesh1d(grid)
    coord=vec(grid[Coordinates])
    ncoord=length(coord)
    points=Vector{Point2f0}(undef,0)
    (xmin,xmax)=extrema(coord)
    h=(xmax-xmin)/40.0
    for i=1:ncoord
        push!(points,Point2f0(coord[i],h))
        push!(points,Point2f0(coord[i],-h))
    end
    points
end

# Point list for intervals
function regionmesh1d(grid,iregion)
    coord=vec(grid[Coordinates])
    points=Vector{Point2f0}(undef,0)
    cn=grid[CellNodes]
    cr=grid[CellRegions]
    ncells=length(cr)
    for i=1:ncells
        if cr[i]==iregion
            push!(points,Point2f0(coord[cn[1,i]],0))
            push!(points,Point2f0(coord[cn[2,i]],0))
        end
    end
    points
end

# Point list for boundary nodes
function bregionmesh1d(grid,ibreg)
    nbfaces=num_bfaces(grid)
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    coord=vec(grid[Coordinates])
    points=Vector{Point2f0}(undef,0)
    (xmin,xmax)=extrema(coord)
    h=(xmax-xmin)/20.0
    for ibface=1:nbfaces
        if bfaceregions[ibface]==ibreg
            push!(points,Point2f0(coord[bfacenodes[1,ibface]],h))
            push!(points,Point2f0(coord[bfacenodes[1,ibface]],-h))
        end
    end
    points
end

# Point list for scene size
function scenecorners1d(grid)
    coord=vec(grid[Coordinates])
    (xmin,xmax)=extrema(coord)
    h=(xmax-xmin)/40.0
    [Point2f0(xmin,-5*h),Point2f0(xmax,5*h)]
end

function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{1}}, grid)

    Makie=ctx[:Plotter]
    nregions=num_cellregions(grid)
    nbregions=num_bfaceregions(grid)

    if !haskey(ctx,:scene)
        ctx[:scene]=Makie.Axis(ctx[:figure];
                               yticklabelsvisible=false,
                               yticksvisible=false,
                               title=ctx[:title],
                               scenekwargs(ctx)...)

        ctx[:grid]=Makie.Node(grid)
        cmap=region_cmap(nregions)
        bcmap=bregion_cmap(nbregions)

        # Set scene size with invisible markers
        Makie.scatter!(ctx[:scene],
                       Makie.lift(g->scenecorners1d(grid),ctx[:grid]),
                       color=:white,
                       markersize=0.0,
                       strokewidth=0)

        # Draw node markers
        Makie.linesegments!(ctx[:scene],
                            Makie.lift(g->basemesh1d(g),ctx[:grid]),
                            color=:black)

        # Colored cell regions
        for i=1:nregions
            Makie.linesegments!(ctx[:scene],
                                Makie.lift(g->regionmesh1d(g,i), ctx[:grid]),
                                color=cmap[i],
                                linewidth=4,
                                label="c $(i)")
        end

        # Colored boundary grid
        for i=1:nbregions
            Makie.linesegments!(ctx[:scene],
                                Makie.lift(g->bregionmesh1d(g,i),ctx[:grid]),
                                color=bcmap[i],
                                linewidth=4,
                                label="b$(i)")
        end

        # Legende
        if ctx[:legend]!=:none
            pos=ctx[:legend]==:best ? :rt : ctx[:legend]
            Makie.axislegend(ctx[:scene],
                             position=pos,
                             labelsize=0.5*ctx[:fontsize],
                             nbanks=5)
        end

        add_scene!(ctx, ctx[:scene])
    else
        ctx[:grid][]=grid
    end
    reveal(ctx,TP)
end

########################################################################
# 1D function

function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{1}}, grid,func)
    Makie=ctx[:Plotter]

    if ctx[:title]==""
        ctx[:title]=" "
    end

    # ... keep this for the case we are unsorted
    function polysegs(grid,func)
        points=Vector{Point2f0}(undef,0)
        cellnodes=grid[CellNodes]
        coord=grid[Coordinates]
        for icell=1:num_cells(grid)
            i1=cellnodes[1,icell]
            i2=cellnodes[2,icell]
            x1=coord[1,i1]
            x2=coord[1,i2]
            push!(points,Point2f0(x1,func[i1]))
            push!(points,Point2f0(x2,func[i2]))
        end
        points
    end

    function polyline(grid,func)
        coord=grid[Coordinates]
        points=[Point2f0(coord[1,i],func[i]) for i=1:length(func)]
    end

    coord=grid[Coordinates]
    xlimits=ctx[:xlimits]
    ylimits=ctx[:flimits]
    xmin=coord[1,1]
    xmax=coord[1,end]
    if xlimits[1]<xlimits[2]
        xmin=xlimits[1]
        xmax=xlimits[2]
    end

    (ymin,ymax)=extrema(func)
    if ylimits[1]<ylimits[2]
        ymin=ylimits[1]
        ymax=ylimits[2]
    end

    function update_lines(ctx)
        if ctx[:markershape]==:none
            #line without marker
            Makie.lines!(ctx[:scene],Makie.lift(a->a, ctx[:lines][end]),
                         linestyle=ctx[:linestyle],
                         linewidth=ctx[:linewidth],
                         color=RGB(ctx[:color]),
                         label=ctx[:label])
        else
            # line with markers separated by markevery

            # draw plain line without the label
            Makie.lines!(ctx[:scene],Makie.lift(a->a, ctx[:lines][end]),
                         linestyle=ctx[:linestyle],
                         color=RGB(ctx[:color]),
                         linewidth=ctx[:linewidth])

            # draw markers without label
            Makie.scatter!(ctx[:scene],
                           Makie.lift(a->a[1:ctx[:markevery]:end],ctx[:lines][end]),
                           color=RGB(ctx[:color]),
                           marker=ctx[:markershape],
                           markercolor=RGB(ctx[:color]),
                           markersize=ctx[:markersize])

            # Draw  dummy line with marker on top ot the first
            # marker position already drawn in order to
            # get the proper legend entry
            if ctx[:label]!=""
                Makie.scatterlines!(ctx[:scene],
                                    Makie.lift(a->a[1:1], ctx[:lines][end]),
                                    linestyle=ctx[:linestyle],
                                    linewidth=ctx[:linewidth],
                                    marker=ctx[:markershape],
                                    markersize=ctx[:markersize],
                                    markercolor=RGB(ctx[:color]),
                                    color=RGB(ctx[:color]),
                                    label=ctx[:label])
            end
        end

        if ctx[:legend]!=:none
            pos=ctx[:legend]==:best ? :rt : ctx[:legend]
            Makie.axislegend(ctx[:scene],position=pos,labelsize=0.5*ctx[:fontsize],backgroundcolor=:transparent)
        end

    end

    if !haskey(ctx,:scene)

        ctx[:xtitle]=Makie.Node(ctx[:title])

        # Axis
        ctx[:scene]=Makie.Axis(ctx[:figure];
                               title=Makie.lift(a->a,ctx[:xtitle]),
                               scenekwargs(ctx)...)
        # Plot size
        Makie.scatter!(ctx[:scene],
                       [Point2f0(xmin,ymin),Point2f0(xmax,ymax)],
                       color=:white,
                       markersize=0.0,
                       strokewidth=0)

        # ctx[:lines]  is an array of lines to draw
        # Here, we start just with the first one.
        ctx[:lines]=[Makie.Node(polyline(grid,func))]

        update_lines(ctx)

        Makie.reset_limits!(ctx[:scene])

        ctx[:nlines]=1

        add_scene!(ctx,ctx[:scene])

    else

        if ctx[:clear]
            ctx[:nlines]=1
        else
            ctx[:nlines]+=1
        end

        p=polyline(grid,func)
        # Either update existing line, or
        # create new one. This works with repeating sequences of
        # updating lines.
        if ctx[:nlines]<=length(ctx[:lines])
            ctx[:lines][ctx[:nlines]][]=p
        else
            push!(ctx[:lines],Makie.Node(p))
            update_lines(ctx)
        end

        Makie.reset_limits!(ctx[:scene])

        ctx[:xtitle][]=ctx[:title]

    end
    reveal(ctx,TP)
end

#######################################################################################
# 2D grid

"""
   makescene2d(ctx)

Complete scene with title and status line showing interaction state.
This uses a gridlayout and its  protrusion capabilities.
"""
function makescene2d(ctx)
    Makie=ctx[:Plotter]
    GL=Makie.GridLayout(ctx[:figure])
    GL[1,1]=ctx[:scene]
    if ctx[:colorbar]==:vertical
        GL[1,2]=Makie.Colorbar(ctx[:figure],ctx[:poly],width=15, textsize=0.5*ctx[:fontsize],ticklabelsize=0.5*ctx[:fontsize])
    elseif ctx[:colorbar]==:horizontal
        GL[2,1]=Makie.Colorbar(ctx[:figure],ctx[:poly],height=15, textsize=0.5*ctx[:fontsize],ticklabelsize=0.5*ctx[:fontsize],vertical=false)
    end
    GL
end

function makescene2d_grid(ctx)
    Makie=ctx[:Plotter]
    GL=Makie.GridLayout(ctx[:figure])
    GL[1,1]=ctx[:scene]
    ncol=length(ctx[:cmap])
    if ctx[:colorbar]==:vertical
        GL[1,2]=Makie.Colorbar(ctx[:figure],
                               colormap=Makie.cgrad(ctx[:cmap],categorical=true),
                               limits=(1,ncol),
                               width=15, textsize=0.5*ctx[:fontsize],ticklabelsize=0.5*ctx[:fontsize])
    elseif ctx[:colorbar]==:horizontal
        GL[2,1]=Makie.Colorbar(ctx[:figure],
                               colormap=Makie.cgrad(ctx[:cmap],categorical=true),
                               limits=(1,ncol),
                               heigth=15, textsize=0.5*ctx[:fontsize],ticklabelsize=0.5*ctx[:fontsize],
                               vertical=false)
    end
    GL
end

function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{2}},grid)

    Makie=ctx[:Plotter]

    nregions=num_cellregions(grid)

    nbregions=num_bfaceregions(grid)

    if !haskey(ctx,:scene)
        ctx[:scene]=Makie.Axis(ctx[:figure];
                               title=ctx[:title],
                               aspect=Makie.DataAspect(),
                               scenekwargs(ctx)...)

        ctx[:grid]=Makie.Node(grid)

        # Draw cells with region mark
        cmap=region_cmap(nregions)
        ctx[:cmap]=cmap
        for i=1:nregions
            Makie.poly!(ctx[:scene],Makie.lift(g->regionmesh(g,i), ctx[:grid]) ,
                        color=cmap[i],
                        strokecolor=:black,
                        strokewidth=ctx[:linewidth])
        end

        # Draw boundary lines
        bcmap=bregion_cmap(nbregions)
        for i=1:nbregions
            Makie.linesegments!(ctx[:scene],
                                Makie.lift(g->bfacesegments(g,i),ctx[:grid]),
                                label="$(i)",
                                color=bcmap[i],
                                linewidth=4)
        end
        # Describe legend
        if ctx[:legend]!=:none
            pos=ctx[:legend]==:best ? :rt : ctx[:legend]
            Makie.axislegend(ctx[:scene],
                             position=pos,
                             labelsize=0.5*ctx[:fontsize],
                             backgroundcolor=:transparent)
        end

        add_scene!(ctx, makescene2d_grid(ctx))
    else
        ctx[:grid][]=grid
    end
    reveal(ctx,TP)
end

# 2D function
function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{2}},grid, func)
    Makie=ctx[:Plotter]

    # Create GeometryBasics.mesh from grid data.
    function make_mesh(grid::ExtendableGrid,func,elevation)
        coord=grid[Coordinates]
        npoints=num_nodes(grid)
        cellnodes=grid[CellNodes]
        if elevation ≈ 0.0
            points=[Point2f0(coord[1,i],coord[2,i]) for i=1:npoints]
        else
            points=[Point3f0(coord[1,i],coord[2,i],func[i]*elevation) for i=1:npoints]
        end
        faces=[TriangleFace(cellnodes[1,i],cellnodes[2,i],cellnodes[3,i]) for i=1:size(cellnodes,2)]
        Mesh(points,faces)
    end

    # Calculate isolevel values
    function isolevels(ctx,func)
        flimits=ctx[:flimits]
        if flimits[1]<flimits[2]
            collect(LinRange(flimits[1],flimits[2],ctx[:isolines]))
        else
            limits=extrema(func)
            collect(LinRange(limits[1],limits[2],ctx[:isolines]))
        end
    end

    flimits=ctx[:flimits]
    if flimits[1]<flimits[2]
        crange=flimits
    else
        crange=extrema(func)
    end

    if !haskey(ctx,:scene)
        # Here, we work with data tuples in the nodes,
        # these  can be easily extended.
        ctx[:data]=Makie.Node((g=grid,f=func,e=ctx[:elevation],t=ctx[:title],l=isolevels(ctx,func),c=crange))

        # would need to switch to Axis3 for supporting elevtion
        ctx[:scene]=Makie.Axis(ctx[:figure];
                               title=Makie.lift(data->data.t,ctx[:data]),
                               aspect=Makie.DataAspect(),
                               scenekwargs(ctx)...)

        # Draw the mesh for the cells
        ctx[:poly]=Makie.poly!(ctx[:scene],
                               Makie.lift(data->make_mesh(data.g,data.f,data.e),ctx[:data]),
                               color=Makie.lift(data->data.f,ctx[:data]),
                               colorrange=Makie.lift(data->data.c,ctx[:data]),
                               colormap=ctx[:colormap])

        # draw the isolines via marching triangles
        Makie.linesegments!(ctx[:scene],
                            Makie.lift(data->marching_triangles(data.g,data.f,data.l),ctx[:data]),
                            color=:black,
                            linewidth=ctx[:linewidth])

        add_scene!(ctx,makescene2d(ctx))

    else
        # Just refresh the data.
        ctx[:data][]=(g=grid,f=func,e=ctx[:elevation],t=ctx[:title],l=isolevels(ctx,func),c=crange)
    end
    reveal(ctx,TP)
end

#######################################################################################
#######################################################################################
# 3D Grid

function xyzminmax(grid::ExtendableGrid)
    coord=grid[Coordinates]
    ndim=size(coord,1)
    xyzmin=zeros(ndim)
    xyzmax=ones(ndim)
    for idim=1:ndim
        @views mn,mx=extrema(coord[idim,:])
        xyzmin[idim]=mn
        xyzmax[idim]=mx
    end
    xyzmin,xyzmax
end

"""
   makeaxis3d(ctx)

Dispatch between LScene and new Axis3. Axis3 does not allow zoom, so we
support LScene in addition.
"""
function makeaxis3d(ctx)
    Makie=ctx[:Plotter]
    if ctx[:scene3d]=="LScene"
        # "Old" LScene with zoom-in functionality
        Makie.LScene(ctx[:figure])
    else
        # "New" Axis3 with prospective new stuff by Julius.
        Makie.Axis3(ctx[:figure];
                    aspect=:data,
                    viewmode=:fitzoom,
                    elevation=ctx[:elev]*π/180,
                    azimuth=ctx[:azim]*π/180,
                    perspectiveness=ctx[:perspectiveness],
                    title=Makie.lift(data->data.t,ctx[:data]),
                    scenekwargs(ctx)...)
    end
end

"""
       makescene3d(ctx)

Complete scene with title and status line showing interaction state.
This uses a gridlayout and its  protrusion capabilities.
"""

function makescene3d(ctx)
    Makie=ctx[:Plotter]
    GL=Makie.GridLayout(ctx[:figure];default_rowgap=0)
    if ctx[:scene3d]=="LScene"
        # LScene has no title, put the title into protrusion space on top  of the scene
        GL[1,1,Makie.Top()]=Makie.Label(ctx[:figure],
                                        " $(Makie.lift(data->data.t,ctx[:data])) ",
                                        tellwidth=false,
                                        height=30,
                                        textsize=ctx[:fontsize])
    end
    GL[1,1]=ctx[:scene]
    # Horizontal or vertical colorbar
    if haskey(ctx,:mesh)
        if ctx[:colorbar]==:vertical
            GL[1,2]=Makie.Colorbar(ctx[:figure],
                                   ctx[:mesh],
                                   width=15,
                                   textsize=0.5*ctx[:fontsize],
                                   ticklabelsize=0.5*ctx[:fontsize])
        elseif ctx[:colorbar]==:horizontal
            GL[2,1]=Makie.Colorbar(ctx[:figure],
                                   ctx[:mesh],
                                   height=15,
                                   textsize=0.5*ctx[:fontsize],
                                   ticklabelsize=0.5*ctx[:fontsize],
                                   vertical=false)
        end
    end
    # Put the status label into protrusion space on the bottom of the scene
    GL[1,1,Makie.Bottom()]=Makie.Label(ctx[:figure],
                                       ctx[:status],
                                       tellwidth=false,
                                       height=30,
                                       textsize=0.5*ctx[:fontsize])
    GL
end

const keyboardhelp=
"""
    Keyboard interactions:
          x: control xplane
          y: control yplane
          z: control zplane
          l: control isolevel
    up/down: fine control control value
pgup/pgdown: coarse control control value
          h: print this message
"""

function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid)

    make_mesh(pts,fcs)=Mesh(meta(pts,normals=normals(pts, fcs)),fcs)

    nregions=num_cellregions(grid)
    nbregions=num_bfaceregions(grid)

    Makie=ctx[:Plotter]
    xyzmin,xyzmax=xyzminmax(grid)
    xyzstep=(xyzmax-xyzmin)/100

    function adjust_planes()
        ctx[:xplane]=max(xyzmin[1],min(xyzmax[1],ctx[:xplane]) )
        ctx[:yplane]=max(xyzmin[2],min(xyzmax[2],ctx[:yplane]) )
        ctx[:zplane]=max(xyzmin[3],min(xyzmax[3],ctx[:zplane]) )
    end

    adjust_planes()

    if !haskey(ctx,:scene)

        ctx[:data]=Makie.Node((g=grid,x=ctx[:xplane],y=ctx[:yplane],z=ctx[:zplane],t=ctx[:title]))

        ctx[:scene]=makeaxis3d(ctx)

        ############# Interior cuts
        # We draw a mesh for each color.
        if ctx[:interior]
            cmap=region_cmap(nregions)
            ctx[:celldata]=Makie.lift(
                d->extract_visible_cells3D(d.g,
                                           (d.x,d.y,d.z),
                                           primepoints=hcat(xyzmin,xyzmax),
                                           Tp=Point3f0,
                                           Tf=GLTriangleFace),
                ctx[:data])

            ctx[:cellmeshes]=Makie.lift(d->[make_mesh(d[1][i],d[2][i]) for i=1:nregions], ctx[:celldata])

            for i=1:nregions
                Makie.mesh!(ctx[:scene],Makie.lift(d->d[i], ctx[:cellmeshes]),
                            color=cmap[i],
                            backlight=1f0
                            )

                if ctx[:linewidth]>0
                    Makie.wireframe!(ctx[:scene],Makie.lift(d->d[i], ctx[:cellmeshes]),
                                     strokecolor=:black,
                                     strokewidth=ctx[:linewidth],
                                     linewidth=ctx[:linewidth],
                                     )
                end
            end
        end

        ############# Visible boundary faces
        bcmap=bregion_cmap(nbregions)

        ctx[:facedata]=Makie.lift(
            d->extract_visible_bfaces3D(d.g,
                                        (d.x,d.y,d.z),
                                        primepoints=hcat(xyzmin,xyzmax),
                                        Tp=Point3f0,
                                        Tf=GLTriangleFace),
            ctx[:data])

        ctx[:facemeshes]=Makie.lift(d->[make_mesh(d[1][i],d[2][i]) for i=1:nbregions], ctx[:facedata])

        for i=1:nbregions
            Makie.mesh!(ctx[:scene],Makie.lift(d->d[i], ctx[:facemeshes]),
                        color=bcmap[i],
                        backlight=1f0
                        )
            if ctx[:linewidth]>0
                Makie.wireframe!(ctx[:scene],Makie.lift(d->d[i], ctx[:facemeshes]),
                                 strokecolor=:black,
                                 linewidth=ctx[:linewidth])
            end
        end

        ############# Transparent outline

        if ctx[:outline]

            ctx[:outlinedata]=Makie.lift(d->extract_visible_bfaces3D(d.g,
                                                                     xyzmax,
                                                                     primepoints=hcat(xyzmin,xyzmax),
                                                                     Tp=Point3f0,
                                                                     Tf=GLTriangleFace),
                                         ctx[:data])
            ctx[:outlinemeshes]=Makie.lift(d->[make_mesh(d[1][i],d[2][i]) for i=1:nbregions], ctx[:outlinedata])

            for i=1:nbregions
                Makie.mesh!(ctx[:scene],Makie.lift(d->d[i], ctx[:outlinemeshes]),
                            color=(bcmap[i],ctx[:alpha]),
                            transparency=true,
                            backlight=1f0
                            )
            end
        end

        ##### Interaction
        scene_interaction(ctx[:scene].scene,Makie,[:z,:y,:x,:q]) do delta,key
            if key==:x
                ctx[:xplane]+=delta*xyzstep[1]
                ctx[:status][]=@sprintf("x=%.3g",ctx[:xplane])
            elseif key==:y
                ctx[:yplane]+=delta*xyzstep[2]
                ctx[:status][]=@sprintf("y=%.3g",ctx[:yplane])
            elseif key==:z
                ctx[:zplane]+=delta*xyzstep[3]
                ctx[:status][]=@sprintf("z=%.3g",ctx[:zplane])
            elseif key==:q
                ctx[:status][]=" "
            end
            adjust_planes()

            ctx[:data][]=(g=grid,x=ctx[:xplane],y=ctx[:yplane],z=ctx[:zplane],t=ctx[:title])
        end

        ctx[:status]=Makie.Node(" ")

        add_scene!(ctx,makescene3d(ctx))

    else
        ctx[:data][]=(g=grid,x=ctx[:xplane],y=ctx[:yplane],z=ctx[:zplane],t=ctx[:title])
    end

    reveal(ctx,TP)
end

# 3d function
function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid , func)

    make_mesh(pts,fcs)=Mesh(pts,fcs)

    function make_mesh(pts,fcs,vals)
        colors = Makie.Makie.interpolated_getindex.((cmap,), vals, (fminmax,))
        GeometryBasics.Mesh(meta(pts, color=colors,normals=normals(pts, fcs)), fcs)
    end

    nregions=num_cellregions(grid)
    nbregions=num_bfaceregions(grid)

    Makie=ctx[:Plotter]
    cmap = Makie.to_colormap(ctx[:colormap])
    xyzmin,xyzmax=xyzminmax(grid)
    xyzstep=(xyzmax-xyzmin)/100

    fminmax=extrema(func)

    flimits=ctx[:flimits]
    if flimits[1]<flimits[2]
        fminmax=(flimits[1],flimits[2])
    end

    fstep=(fminmax[2]-fminmax[1])/100
    if fstep≈0
        fstep=0.1
    end

    function adjust_planes()
        ctx[:xplane]=max(xyzmin[1],min(xyzmax[1],ctx[:xplane]) )
        ctx[:yplane]=max(xyzmin[2],min(xyzmax[2],ctx[:yplane]) )
        ctx[:zplane]=max(xyzmin[3],min(xyzmax[3],ctx[:zplane]) )
        ctx[:flevel]=max(fminmax[1],min(fminmax[2],ctx[:flevel]))
    end

    adjust_planes()

    makeplanes(x,y,z)=[[1,0,0,-x],
                       [0,1,0,-y],
                       [0,0,1,-z]]

    if !haskey(ctx,:scene)

        ctx[:data]=Makie.Node((g=grid,f=func,x=ctx[:xplane],y=ctx[:yplane],z=ctx[:zplane],l=ctx[:flevel],t=ctx[:title]))

        ctx[:scene]=makeaxis3d(ctx)

        #### Transparent outlne
        if ctx[:outline]
            ctx[:outlinedata]=Makie.lift(d->extract_visible_bfaces3D(d.g,
                                                                     xyzmax,
                                                                     primepoints=hcat(xyzmin,xyzmax),
                                                                     Tp=Point3f0,
                                                                     Tf=GLTriangleFace),
                                         ctx[:data])
            ctx[:facemeshes]=Makie.lift(d->[make_mesh(d[1][i],d[2][i]) for i=1:nbregions], ctx[:outlinedata])
            bcmap=bregion_cmap(nbregions)
            for i=1:nbregions
                Makie.mesh!(ctx[:scene],Makie.lift(d->d[i], ctx[:facemeshes]),
                            color=(bcmap[i],ctx[:alpha]),
                            transparency=true,
                            backlight=1f0
                            )
            end
        end

        #### Plane sections and isosurfaces
        Makie.mesh!(ctx[:scene],
                    Makie.lift(d->make_mesh(marching_tetrahedra(d.g,
                                                                d.f,
                                                                makeplanes(d.x,d.y,d.z),
                                                                [d.l],
                                                                primepoints=hcat(xyzmin,xyzmax),
                                                                primevalues=fminmax,
                                                                Tp=Point3f0,
                                                                Tf=GLTriangleFace,
                                                                Tv=Float32)...),ctx[:data]),
                    backlight=1f0)

        #### Interactions
        scene_interaction(ctx[:scene].scene,Makie,[:z,:y,:x,:l,:q]) do delta,key
            if key==:x
                ctx[:xplane]+=delta*xyzstep[1]
                ctx[:status][]=@sprintf("x=%.3g",ctx[:xplane])
            elseif key==:y
                ctx[:yplane]+=delta*xyzstep[2]
                ctx[:status][]=@sprintf("y=%.3g",ctx[:yplane])
            elseif key==:z
                ctx[:zplane]+=delta*xyzstep[3]
                ctx[:status][]=@sprintf("z=%.3g",ctx[:zplane])
            elseif key==:l
                ctx[:flevel]+=delta*fstep
                ctx[:status][]=@sprintf("l=%.3g",ctx[:flevel])
            elseif key==:q
                ctx[:status][]=" "
            end
            adjust_planes()
            ctx[:data][]=(g=grid,f=func,x=ctx[:xplane],y=ctx[:yplane],z=ctx[:zplane],l=ctx[:flevel],t=ctx[:title])
        end

        ctx[:status]=Makie.Node(" ")
        add_scene!(ctx,makescene3d(ctx))
    else
        ctx[:data][]=(g=grid,f=func,x=ctx[:xplane],y=ctx[:yplane],z=ctx[:zplane],l=ctx[:flevel],t=ctx[:title])
    end
    reveal(ctx,TP)
end

        # TODO: allow aspect scaling
        # if ctx[:aspect]>1.0
        #     Makie.scale!(ctx[:scene],ctx[:aspect],1.0)
        # else
        #     Makie.scale!(ctx[:scene],1.0,1.0/ctx[:aspect])
        # end

        # TODO: use distinguishable colors
        # http://juliagraphics.github.io/Colors.jl/stable/colormapsandcolorscales/#Generating-distinguishable-colors-1

# TODO: a priori angles aka pyplot3D
# rect = ctx[:scene]
# azim=ctx[:azim]
# elev=ctx[:elev]
# arr = normalize([cosd(azim/2), 0, sind(azim/2), -sind(azim/2)])
# Makie.rotate!(rect, Makie.Quaternionf0(arr...))

# 3 replies
# Julius Krumbiegel  5 hours ago
# try lines(x, y, axis = (limits = lims,)), the other keyword arguments go to the plot
# Julius Krumbiegel  5 hours ago
# although I think it should be targetlimits because the limits are computed from them usually, considering that there might be aspect constraints that should be met or linked axes (edited)
# Christophe Meyer  4 hours ago
# Thanks!  lines(x, y, axis = (targetlimits = lims,))  indeed makes the limits update.^
# I found that autolimits!(axis) gave good results, even better than me manually computing limits!


text:
# SLAM simulation using extended Kalman filter. Inspired by Matlab implementation
# by T. Bailey et al (https://openslam.org/bailey-slam.html).

# Here are three options to make the SLAM module (or any other) available to your code:
#
# 1. Include it directly: `include("path/to/SLAM.jl/src/SLAM.jl")`
# 2. Add it to Julia's LOAD_PATH environment variable: if `module_dir = "/path/to/SLAM.jl"`,
#    add a line like `module_dir in LOAD_PATH || push!(LOAD_PATH, module_dir)` to your code,
#    your REPL session, or to `~/.juliarc.jl`.
# 3. Use SLAM.jl as a package: Pkg.clone("url/to/SLAM.jl")
include("../src/SLAM.jl")

using SLAM

include("sim-utils.jl")

function ekfsim_setup(n_landmarks::Integer, waypoints_file::AbstractString)

    # Scene boundaries: xmin, xmax, ymin, ymax
    boundaries = [0.; 100; 0; 100]
    scene = Scene(boundaries,
                  get_waypoints(waypoints_file),
                  make_landmarks(n_landmarks, boundaries, 0.05),
                  Array{Float64}(3, 10000),
                  Array{Float64}(3, 10000),
                  0)

    # Vehicle
    vehicle = Vehicle()
    vehicle.wheelbase = 4.0            # [m]
    vehicle.max_gamma = 60*pi/180      # [rad] max steering angle (-max < g < max)
    vehicle.steer_rate = 60*pi/180     # [rad/s] max rate of change in steer angle
    vehicle.sensor_range = 30          # [m] Landmark detection radius
    vehicle.shape = [1. -1 -1; 0 1 -1] # A little triangle for visualization
    vehicle.pose = initial_pose(scene)
    vehicle.target_speed = 8           # [m/s]
    vehicle.waypoint_id = 1            # Initialize to index of first waypoint

    # SLAM state vector
    # First three elements comprise the SLAM vehicle pose; state is augmented as
    # new landmarks are observed. Covariance matrix is also augmented during simulation.
    state = EKFSlamState(vehicle.pose, zeros(3, 3))

    SimData(scene, vehicle, state, zeros(2, n_landmarks), 0, false, false, 2)
end

# This method runs the simulation without monitoring
function sim!(simdata::SimData)
    sim!(simdata::SimData, (x...) -> (), [])
end

function sim!(simdata::SimData, monitor::Function, args::AbstractVector)

    # Use "typedefs" for simdata fields to shorten names.
    # Since these are references, the simdata object is modified in place.
    scene, vehicle, state = simdata.scene, simdata.vehicle, simdata.state

    n_landmarks = size(scene.landmarks, 2)

    d_min = 1.0                        # [m] Once inside d_min, head for next waypoint

    # Simulated control noise covariance matrix
    sigma_speed = 0.5                  # [m/s] Uncertainty about target speed
    sigma_steer = (3.0*pi/180)         # [rad] Uncertainty about target gamma
    Q = [sigma_speed^2 0; 0 sigma_steer^2]

    # Simulated observation noise covariance matrix
    sigmaR = 0.1                       # [m] Range uncertainty
    sigmaB = (1.0*pi/180)              # [rad] Bearing angle uncertainty
    R = [sigmaR^2 0; 0 sigmaB^2]       # Obs. covariance matrix

    # Timing variables
    dt = 0.025                         # [s] Interval between control updates
    dt_obs = 8*dt                      # [s] Interval between state updates
    dtsum = 0                          # Time since last observation
    marker = time()

    while vehicle.waypoint_id != 0

        start = marker

        # Update heading and target waypoint
        steer!(vehicle, scene.waypoints, d_min, dt)

        # If lap through waypoints is finished with more to go, start a new lap
        if vehicle.waypoint_id == 0 && simdata.nlaps > 1
            vehicle.waypoint_id = 1
            simdata.nlaps -= 1
        end

        # Advance vehicle's equations of motion by one timestep
        step_vehicle!(vehicle, dt)

        # Simulate imprecision about control target speed and steering angle
        add_control_noise!(vehicle, Q)

        # Prediction update for state vector and covariance
        state.x, state.cov = predict(state, vehicle, Q, dt)

        dtsum += dt

        # EKF update step - takes place every (dt_obs/dt) timesteps
        if dtsum > dt_obs
            dtsum = 0
            simdata.state_updated = true
            simdata.z, tags = get_observations(vehicle, scene, R)
            simdata.nz = size(simdata.z, 2)

            # Last two parameters are thresholds (Mahalanobis distances)
            # 4.0  [m] max distance for association
            # 25.0 [m] min distance for creation of new feature
            zf, idf, zn = associate(state, simdata.z, R, 4.0, 25.0)

            # Update SLAM state
            state.x, state.cov = update(state, zf, R, idf)

            # Add z measurements to the SLAM state vector and its covariance
            state.x, state.cov = add_features(state, zn, R)
        else
            simdata.state_updated = false
        end

        # Update vehicle tracks in scene
        scene.nsteps += 1
        scene.true_track[:, scene.nsteps] = vehicle.pose
        scene.slam_track[:, scene.nsteps] = state.x[1:3]

        monitor(args...)

        # Limit frame rate for realistic timing
        marker = time()
        if start + dt > marker
            sleep(start + dt - marker)
        end

        while simdata.paused
            wait()
        end

    end
end

ext: jl


text:
using Agents
using BenchmarkTools
using Random
using StatsBase: mean

mutable struct Agent1 <: AbstractAgent
  id::Int
  weight::Float64
  type::String
end
function init1()
    Random.seed!(267)
    model = ABM(Agent1;
                properties = Dict(:year => 0, :tick => 0, :flag => false))
    for a in 1:1000
        add_agent!(model, rand(), randstring(4))
    end
    return model
end
function as1!(agent, model)
    nothing
end
function ms1!(model)
    nothing
end


text:
using ConstrainedSystems
using Test
using CartesianGrids
using LinearAlgebra

import ConstrainedSystems: recursivecopy!, needs_iteration, ArrayPartition

const GROUP = get(ENV, "GROUP", "All")

if GROUP == "All" || GROUP == "Utils"
  include("utils.jl")
end

if GROUP == "All" || GROUP == "Types"
  include("types.jl")
end

if GROUP == "All" || GROUP == "Saddle"
  include("saddle.jl")
end

if GROUP == "All" || GROUP == "Convergence"
  include("algconvergence.jl")
end


text:
"""
Provides Genie with encryption and decryption capabilities.
"""
module Encryption

import Genie, Nettle

const ENCRYPTION_METHOD = "AES256"

"""
    encrypt{T}(s::T) :: String

Encrypts `s`.
"""
function encrypt(s::T)::String where T
  (key32, iv16) = encryption_sauce()
  encryptor = Nettle.Encryptor(ENCRYPTION_METHOD, key32)

  Nettle.encrypt(encryptor, :CBC, iv16, Nettle.add_padding_PKCS5(Vector{UInt8}(s), 16)) |> bytes2hex
end

"""
    decrypt(s::String) :: String

Decrypts `s` (a `string` previously encrypted by Genie).
"""
function decrypt(s::String) :: String
  (key32, iv16) = encryption_sauce()
  decryptor = Nettle.Decryptor(ENCRYPTION_METHOD, key32)
  deciphertext = Nettle.decrypt(decryptor, :CBC, iv16, s |> hex2bytes)

  String(Nettle.trim_padding_PKCS5(deciphertext))
end

function encryption_sauce() :: Tuple{Vector{UInt8},Vector{UInt8}}
  if ! isdefined(Genie, :SECRET_TOKEN)
    @warn "Encryption error"

    if ! Genie.Configuration.isprod()
      @warn "Generating temporary secret token"
      Core.eval(Genie, :(const SECRET_TOKEN = $(Genie.Generator.secret_token())))
    else
      error("Can't encrypt - please make sure you run a Genie project and SECRET_TOKEN is defined in config/secrets.jl")
    end
  end

  passwd = Genie.SECRET_TOKEN[1:32]
  salt = hex2bytes(Genie.SECRET_TOKEN[33:64])

  Nettle.gen_key32_iv16(Vector{UInt8}(passwd), salt)
end

end


text:
@testset "Transition process tests" begin

@testset "Basic transition test" begin
    mpSim = ManpowerSimulation()

    attribute = Attribute( "Level" )
    setPossibleAttributeValues!( attribute, ["A", "B"] )
    addSimulationAttribute!( mpSim, attribute )

    node = BaseNode( "A" )
    setNodeRequirements!( node, ["Level"], ["A"] )
    addSimulationBaseNode!( mpSim, node )

    node = BaseNode( "B" )
    setNodeRequirements!( node, ["Level"], ["B"] )
    setNodeTarget!( node, -1 )
    addSimulationBaseNode!( mpSim, node )

    recruitment = Recruitment( "EW" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "A" )
    setRecruitmentFixed!( recruitment, 10 )
    setRecruitmentAgeFixed!( recruitment, 18 * 12 )
    addSimulationRecruitment!( mpSim, recruitment )

    transition = Transition( "trans", "A", "B" )
    setTransitionSchedule!( transition, 60 )
    addSimulationTransition!( mpSim, transition )

    retirement = Retirement()
    setRetirementCareerLength!( retirement, 120 )
    setRetirementSchedule!( retirement, 12 )
    setSimulationRetirement!( mpSim, retirement )

    setSimulationLength!( mpSim, 300.0 )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )

    report = nodeEvolutionReport( mpSim, 12, "A", "B" )
    popReport = report[1]
    fluxReportA = report[2]["A"]
    fluxReportB = report[2]["B"]

    @test all( popReport[:, :A] .== vcat( repeat( 10:10:50, 5 ), 10 ) ) &&
        all( popReport[:, :B] .==
        vcat( fill( 0, 5 ), fill( 50, 5 ), repeat( 90:-10:50, 3 ), 90 ) )
    @test all( fluxReportA[1][:, Symbol( "EW: external => A" )] .== 10 ) &&
        all( fluxReportA[2][:, Symbol( "trans: A => B" )] .==
        vcat( 0, repeat( [0, 0, 0, 0, 50], 5 ) ) )
    @test all( fluxReportB[1][:, Symbol( "trans: A => B" )] .==
        vcat( 0, repeat( [0, 0, 0, 0, 50], 5 ) ) ) &&
        all( fluxReportB[2][:, Symbol( "retirement: B => external" )] .==
        vcat( fill( 0, 10 ), fill( 10, 16 ) ) )

    clearSimulationTransitions!( mpSim )
    setTransitionSchedule!( transition, 48 )
    addSimulationTransition!( mpSim, transition )

    run( mpSim, saveConfig = false )
    report = nodePopReport( mpSim, 12, "A", "B" )
    @test all( report[:, :A] .== vcat( repeat( 10:10:40, 6 ), 10, 20 ) ) &&
        all( report[:, :B] .==
        vcat( fill( 0, 4 ), fill( 40, 4 ), 80, repeat( [80, 70, 60, 90], 4 ),
        80 ) )

    clearSimulationTransitions!( mpSim )
    setTransitionSchedule!( transition, 60, 24 )
    addSimulationTransition!( mpSim, transition )

    run( mpSim, saveConfig = false )
    report = nodePopReport( mpSim, 12, "A", "B" )
    @test all( report[:, :A] .==
        vcat( 10, 20, repeat( 10:10:50, 4 ), 10:10:40 ) ) &&
        all( report[:, :B] .== vcat( 0, 0, fill( 20, 5 ), 70, 70,
        repeat( [70, 60, 50, 90, 80], 3 ), 70, 60 ) )
end  # @testset "Basic transition test"
#=
@testset "OUT transition test" begin
    mpSim = ManpowerSimulation()

    attribute = Attribute( "Level" )
    setPossibleAttributeValues!( attribute, ["A"] )
    addSimulationAttribute!( mpSim, attribute )

    node = BaseNode( "A" )
    setNodeRequirements!( node, ["Level"], ["A"] )
    addSimulationBaseNode!( mpSim, node )

    recruitment = Recruitment( "EW" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "A" )
    setRecruitmentFixed!( recruitment, 10 )
    setRecruitmentAgeFixed!( recruitment, 18 * 12 )
    addSimulationRecruitment!( mpSim, recruitment )

    transition = Transition( "PE", "A" )
    setTransitionSchedule!( transition, 60 )
    addSimulationTransition!( mpSim, transition )

    setSimulationLength!( mpSim, 300.0 )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )

    report = nodeEvolutionReport( mpSim, 12, "A" )
    popReport = report[1]
    fluxReport = report[2]["A"]
    @test all( popReport[:, :A] .== vcat( repeat( 10:10:50, 5 ), 10 ) )
    @test all( fluxReport[1][:, Symbol( "EW: external => A" )] .== 10 ) &&
        all( fluxReport[2][:, Symbol( "PE: A => external" )] .==
        vcat( 0, repeat( [0, 0, 0, 0, 50], 5 ) ) )

    clearSimulationTransitions!( mpSim )
    setTransitionSchedule!( transition, 48 )
    addSimulationTransition!( mpSim, transition )

    run( mpSim, saveConfig = false )
    report = nodePopReport( mpSim, 12, "A" )
    @test all( report[:, :A] .== vcat( repeat( 10:10:40, 6 ), 10, 20 ) )

    clearSimulationTransitions!( mpSim )
    setTransitionSchedule!( transition, 60, 24 )
    addSimulationTransition!( mpSim, transition )

    run( mpSim, saveConfig = false )
    report = nodePopReport( mpSim, 12, "A" )
    @test all( report[:, :A] .== vcat( 10, 20, repeat( 10:10:50, 4 ),
        10:10:40 ) )
end  # @testset "OUT transition test"

@testset "Time conditions test" begin
    mpSim = ManpowerSimulation()

    attribute = Attribute( "Level" )
    setPossibleAttributeValues!( attribute, ["A", "B", "C"] )
    addSimulationAttribute!( mpSim, attribute )

    node = BaseNode( "A" )
    setNodeRequirements!( node, ["Level"], ["A"] )
    addSimulationBaseNode!( mpSim, node )

    node = BaseNode( "B" )
    setNodeRequirements!( node, ["Level"], ["B"] )
    setNodeTarget!( node, -1 )
    addSimulationBaseNode!( mpSim, node )

    node = BaseNode( "C" )
    setNodeRequirements!( node, ["Level"], ["C"] )
    setNodeTarget!( node, -1 )
    addSimulationBaseNode!( mpSim, node )

    recruitment = Recruitment( "EW" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "A" )
    setRecruitmentFixed!( recruitment, 10 )
    setRecruitmentAgeFixed!( recruitment, 18 * 12 )
    addSimulationRecruitment!( mpSim, recruitment )

    transition = Transition( "trans", "A", "B" )
    setTransitionSchedule!( transition, 24 )
    addTransitionCondition!( transition, MPcondition( "tenure", >=, 36 ) )
    addSimulationTransition!( mpSim, transition )

    retirement = Retirement()
    setRetirementCareerLength!( retirement, 120 )
    setRetirementSchedule!( retirement, 12 )
    setSimulationRetirement!( mpSim, retirement )

    setSimulationLength!( mpSim, 300.0 )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )

    report = nodeEvolutionReport( mpSim, 12, "A", "B" )
    popReport = report[1]
    fluxReportA = report[2]["A"]
    fluxReportB = report[2]["B"]
    @test all( popReport[:, :A] .==
        vcat( 10, 20, repeat( [30, 40], 12 ) ) ) &&
        all( popReport[:, :B] .==
        vcat( fill( 0, 4 ), 20, 20, 40, 40, 60, 60, repeat( [70, 60], 8 ) ) )
    @test all( fluxReportA[1][:, Symbol( "EW: external => A" )] .== 10 ) &&
        all( fluxReportA[2][:, Symbol( "trans: A => B" )] .==
        vcat( fill( 0, 4 ), repeat( [20, 0], 11 ) ) )
    @test all( fluxReportB[1][:, Symbol( "trans: A => B" )] .==
        vcat( fill( 0, 4 ), repeat( [20, 0], 11 ) ) ) &&
        all( fluxReportB[2][:, Symbol( "retirement: B => external" )] .==
        vcat( fill( 0, 10 ), fill( 10, 16 ) ) )

    clearSimulationTransitions!( mpSim )
    clearTransitionConditions!( transition )
    addTransitionCondition!( transition, MPcondition( "age", ==, 21 * 12 ) )
    addSimulationTransition!( mpSim, transition )

    run( mpSim, saveConfig = false )
    report = nodePopReport( mpSim, 12, "A", "B" )
    @test all( report[:, :A] .==
        vcat( 10:10:40, 40, 50, 50, 60, repeat( [60, 70], 9 ) ) ) &&
        all( report[:, :B] .==
        vcat( fill( 0, 4 ), 10, 10, 20, 20, 30, 30, repeat( [40, 30], 8 ) ) )

    clearSimulationTransitions!( mpSim )
    clearTransitionConditions!( transition )
    setTransitionSchedule!( transition, 12 )
    addTransitionCondition!( transition, MPcondition( "tenure", >=, 36 ) )
    addSimulationTransition!( mpSim, transition )

    transition = Transition( "trans", "B", "C" )
    setTransitionSchedule!( transition, 24 )
    addTransitionCondition!( transition, MPcondition( "tenure", >=, 60 ),
        MPcondition( "time in node", <=, 24 ) )
    addSimulationTransition!( mpSim, transition )

    run( mpSim, saveConfig = false )
    report = nodePopReport( mpSim, 12, "A", "B", "C" )
    @test all( report[:, :A] .== vcat( 10, 20, fill( 30, 24 ) ) ) &&
        all( report[:, :B] .==
        vcat( fill( 0, 3 ), 10:10:30, 30, 40, repeat( [40, 50], 9 ) ) ) &&
        all( report[:, :C] .==
        vcat( fill( 0, 6 ), 10, 10, 20, 20, repeat( [30, 20], 8 ) ) )
end  # @testset "Time conditions test"

@testset "Extra conditions test" begin
    mpSim = ManpowerSimulation()

    attribute = Attribute( "Level" )
    setPossibleAttributeValues!( attribute, ["A", "B", "C"] )
    addSimulationAttribute!( mpSim, attribute )

    attribute = Attribute( "Branch" )
    setPossibleAttributeValues!( attribute, ["A", "B"] )
    addSimulationAttribute!( mpSim, attribute )

    node = BaseNode( "A" )
    setNodeRequirements!( node, ["Level", "Branch"], ["A", "A"] )
    addSimulationBaseNode!( mpSim, node )

    node = BaseNode( "B" )
    setNodeRequirements!( node, ["Level", "Branch"], ["A", "B"] )
    addSimulationBaseNode!( mpSim, node )

    node = BaseNode( "C" )
    setNodeRequirements!( node, ["Level"], ["B"] )
    setNodeTarget!( node, -1 )
    addSimulationBaseNode!( mpSim, node )

    node = BaseNode( "D" )
    setNodeRequirements!( node, ["Level"], ["C"] )
    setNodeTarget!( node, -1 )
    addSimulationBaseNode!( mpSim, node )

    recruitment = Recruitment( "EW" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "A" )
    setRecruitmentFixed!( recruitment, 10 )
    setRecruitmentAgeFixed!( recruitment, 18 * 12 )
    addSimulationRecruitment!( mpSim, recruitment )

    recruitment = Recruitment( "EW" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "B" )
    setRecruitmentFixed!( recruitment, 10 )
    setRecruitmentAgeFixed!( recruitment, 18 * 12 )
    addSimulationRecruitment!( mpSim, recruitment )

    transition = Transition( "trans", "A", "C" )
    setTransitionSchedule!( transition, 12 )
    addTransitionCondition!( transition, MPcondition( "time in node", >=, 24 ) )
    addSimulationTransition!( mpSim, transition )

    transition = Transition( "trans", "B", "C" )
    setTransitionSchedule!( transition, 12 )
    addTransitionCondition!( transition, MPcondition( "time in node", >=, 24 ) )
    addSimulationTransition!( mpSim, transition )

    transition = Transition( "trans", "C", "D" )
    setTransitionSchedule!( transition, 12 )
    addTransitionCondition!( transition,
        MPcondition( "time in node", >=, 36 ),
        MPcondition( "Branch", !=, "B" ) )
    addSimulationTransition!( mpSim, transition )

    retirement = Retirement()
    setRetirementCareerLength!( retirement, 120 )
    setRetirementSchedule!( retirement, 12 )
    setSimulationRetirement!( mpSim, retirement )

    setSimulationLength!( mpSim, 300.0 )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )

    subpop1 = Subpopulation( "Subpop1" )
    addSubpopulationCondition!( subpop1, MPcondition( "Branch", ==, "A" ) )
    subpop2 = Subpopulation( "Subpop2" )
    addSubpopulationCondition!( subpop2, MPcondition( "Branch", ==, "B" ) )
    subpop3 = Subpopulation( "Subpop3" )
    addSubpopulationCondition!( subpop3, MPcondition( "Branch", ==, "A" ) )
    subpop4 = Subpopulation( "Subpop4" )
    addSubpopulationCondition!( subpop4, MPcondition( "Branch", ==, "B" ) )

    setSubpopulationSourceNode!.( [subpop1, subpop2, subpop3, subpop4],
        ["C", "C", "D", "D"] )

    report = nodePopReport( mpSim, 12, "A", "B", "C", "D" )
    @test all( report[:, :A] .== vcat( 10, fill( 20, 25 ) ) ) &&
        all( report[:, :B] .== vcat( 10, fill( 20, 25 ) ) ) &&
        all( report[:, :C] .==
        vcat( 0, 0, 20, 40, 60:10:100, fill( 110, 17 ) ) ) &&
        all( report[:, :D] .==
        vcat( fill( 0, 5 ), 10:10:40, fill( 50, 17 ) ) )

    report = subpopulationPopReport( mpSim, 12, subpop1, subpop2, subpop3,
        subpop4 )
    @test all( report[:, :Subpop1] .==
        vcat( 0, 0, 10, 20, fill( 30, 22 ) ) ) &&
        all( report[:, :Subpop2] .==
        vcat( 0, 0, 10:10:70, fill( 80, 17 ) ) ) &&
        all( report[:, :Subpop3] .==
        vcat( fill( 0, 5 ), 10:10:40, fill( 50, 17 ) ) ) &&
        all( report[:, :Subpop4] .== 0 )
end  # @testset "Extra conditions test"

@testset "Min/max flux test" begin
    mpSim = ManpowerSimulation()

    attribute = Attribute( "Level" )
    setPossibleAttributeValues!( attribute, ["A", "B"] )
    addSimulationAttribute!( mpSim, attribute )

    node = BaseNode( "A" )
    setNodeRequirements!( node, ["Level"], ["A"] )
    addSimulationBaseNode!( mpSim, node )

    node = BaseNode( "B" )
    setNodeRequirements!( node, ["Level"], ["B"] )
    setNodeTarget!( node, 40 )
    addSimulationBaseNode!( mpSim, node )

    recruitment = Recruitment( "EW" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "A" )
    setRecruitmentFixed!( recruitment, 20 )
    setRecruitmentAgeFixed!( recruitment, 18 * 12 )
    addSimulationRecruitment!( mpSim, recruitment )

    transition = Transition( "trans", "A", "B" )
    setTransitionSchedule!( transition, 12 )
    setTransitionFluxLimits!( transition, 5, 10 )
    addTransitionCondition!( transition, MPcondition( "time in node", >=, 36 ) )
    addSimulationTransition!( mpSim, transition )

    transition = Transition( "trans", "B" )
    setTransitionSchedule!( transition, 12 )
    addTransitionCondition!( transition, MPcondition( "tenure", >=, 120 ) )
    addSimulationTransition!( mpSim, transition )

    retirement = Retirement()
    setRetirementCareerLength!( retirement, 120 )
    setRetirementSchedule!( retirement, 12 )
    setSimulationRetirement!( mpSim, retirement )

    setSimulationLength!( mpSim, 300.0 )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )

    report = nodeEvolutionReport( mpSim, 12, "A", "B" )
    popReport = report[1]
    fluxReportB = report[2]["B"]
    @test all( popReport[:, :B] .==
        vcat( 0, 0, 0, 10:10:40, 45:5:55, 50, 45, fill( 40, 14 ) ) )
    @test all( fluxReportB[1][:, Symbol( "other => B" )] .==
        vcat( 0, 0, 0, fill( 10, 4 ), fill( 5, 6 ), 10, fill( 5, 6 ), 10,
        fill( 5, 5 ) ) )
end  # @testset "Min/max flux test"

@testset "# of attempts test" begin
    mpSim = ManpowerSimulation()

    attribute = Attribute( "Level" )
    setPossibleAttributeValues!( attribute, ["A", "B"] )
    addSimulationAttribute!( mpSim, attribute )

    node = BaseNode( "A" )
    setNodeRequirements!( node, ["Level"], ["A"] )
    addSimulationBaseNode!( mpSim, node )

    node = BaseNode( "B" )
    setNodeRequirements!( node, ["Level"], ["B"] )
    addSimulationBaseNode!( mpSim, node )

    recruitment = Recruitment( "EW" )
    setRecruitmentSchedule!( recruitment, 360 )
    setRecruitmentTarget!( recruitment, "A" )
    setRecruitmentFixed!( recruitment, 100 )
    setRecruitmentAgeFixed!( recruitment, 18 * 12 )
    addSimulationRecruitment!( mpSim, recruitment )

    transition = Transition( "trans", "A", "B" )
    setTransitionSchedule!( transition, 12 )
    setTransitionFluxLimits!( transition, 10, 10 )
    setTransitionMaxAttempts!( transition, 1 )
    addTransitionCondition!( transition, MPcondition( "time in node", >=, 60 ) )
    addSimulationTransition!( mpSim, transition )

    setSimulationLength!( mpSim, 300.0 )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )

    report = nodePopReport( mpSim, 12, "A", "B" )
    @test all( report[:, :B] .== vcat( fill( 0, 5 ), fill( 10, 21 ) ) )

    clearSimulationTransitions!( mpSim )
    setTransitionMaxAttempts!( transition, 5 )
    addSimulationTransition!( mpSim, transition )
    run( mpSim, saveConfig = false )

    report = nodePopReport( mpSim, 12, "A", "B" )
    @test all( report[:, :B] .==
        vcat( fill( 0, 5 ), 10:10:40, fill( 50, 17 ) ) )

    clearSimulationTransitions!( mpSim )
    setTransitionMaxAttempts!( transition, -1 )
    addSimulationTransition!( mpSim, transition )
    run( mpSim, saveConfig = false )

    report = nodePopReport( mpSim, 12, "A", "B" )
    @test all( report[:, :B] .==
        vcat( fill( 0, 5 ), 10:10:90, fill( 100, 12 ) ) )
end  # @testset "# of attempts test"

@testset "Extra changes test" begin
    mpSim = ManpowerSimulation()

    attribute = Attribute( "Level" )
    setPossibleAttributeValues!( attribute, ["A", "B"] )
    addSimulationAttribute!( mpSim, attribute )

    attribute = Attribute( "Extra" )
    addInitialAttributeValue!( attribute, ("A", 1.0) )
    addPossibleAttributeValue!( attribute, "B" )
    addSimulationAttribute!( mpSim, attribute )

    node = BaseNode( "A" )
    setNodeRequirements!( node, ["Level"], ["A"] )
    addSimulationBaseNode!( mpSim, node )

    node = BaseNode( "B" )
    setNodeRequirements!( node, ["Level"], ["B"] )
    addSimulationBaseNode!( mpSim, node )

    recruitment = Recruitment( "EW" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "A" )
    setRecruitmentFixed!( recruitment, 10 )
    setRecruitmentAgeFixed!( recruitment, 18 * 12 )
    addSimulationRecruitment!( mpSim, recruitment )

    transition = Transition( "trans", "A", "B" )
    setTransitionSchedule!( transition, 12 )
    setTransitionMaxAttempts!( transition, 1 )
    addTransitionCondition!( transition, MPcondition( "time in node", >=, 36 ) )
    addTransitionAttributeChange!( transition, ("Extra", "B") )
    addSimulationTransition!( mpSim, transition )

    retirement = Retirement()
    setRetirementCareerLength!( retirement, 120 )
    setRetirementSchedule!( retirement, 12 )
    setSimulationRetirement!( mpSim, retirement )

    setSimulationLength!( mpSim, 300.0 )

    subpop = Subpopulation( "B" )
    addSubpopulationCondition!( subpop, MPcondition( "Extra", ==, "B" ) )
    setSubpopulationSourceNode!( subpop, "B" )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig = false )

    popReport = nodePopReport( mpSim, 12, "A", "B" )
    subpopReport = subpopulationPopReport( mpSim, 12, subpop )
    @test all( popReport[:, :B] .== subpopReport[:, :B] )

    clearSimulationTransitions!( mpSim )
    addTransitionAttributeChange!( transition, ("Level", "A") )
    addSimulationTransition!( mpSim, transition )

    subpop = Subpopulation( "A" )
    addSubpopulationCondition!( subpop, MPcondition( "Level", ==, "A" ) )
    setSubpopulationSourceNode!( subpop, "B" )

    run( mpSim, saveConfig = false )
    popReport = nodePopReport( mpSim, 12, "A", "B" )
    subpopReport = subpopulationPopReport( mpSim, 12, subpop )
    @test all( popReport[:, :B] .== subpopReport[:, :A] )

end  # @testset "Extra changes test"
=#
end  # @testset "Transition process tests"


text:
using BinaryBuilder

# Collection of sources required to build Binaryen
sources = [
    "https://github.com/WebAssembly/binaryen/archive/1.37.36.tar.gz" =>
    "e77520b31e81befe7767f7179243983d678cda6d22cb07110a92aa2ff38831f7",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
ls -l
cd binaryen-1.37.36/
if [ $target = "x86_64-apple-darwin14" ]; then
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_FLAGS="--target=x86_64-apple-darwin14" -DCMAKE_C_FLAGS="--target=x86_64-apple-darwin14"
elif [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_FLAGS="-static -no-pie" -DCMAKE_C_FLAGS="-static -no-pie"
else
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_BUILD_TYPE=Release
fi
make -j8
make install

if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
mv $WORKSPACE/destdir/lib/*.dll $WORKSPACE/destdir/bin
fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi),
    BinaryProvider.MacOS(),
    BinaryProvider.Windows(:x86_64),
    BinaryProvider.Windows(:i686),
    BinaryProvider.Linux(:i686, :glibc, :blank_abi),
    BinaryProvider.Linux(:aarch64, :glibc, :blank_abi),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf),
    BinaryProvider.Linux(:powerpc64le, :glibc, :blank_abi),
    BinaryProvider.Linux(:i686, :musl, :blank_abi),
    BinaryProvider.Linux(:x86_64, :musl, :blank_abi),
    BinaryProvider.Linux(:aarch64, :musl, :blank_abi),
    BinaryProvider.Linux(:armv7l, :musl, :eabihf)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "", :s2wasm),
    ExecutableProduct(prefix, "", :wasm_reduce),
    ExecutableProduct(prefix, "", :wasm_emscripten_finalize),
    ExecutableProduct(prefix, "", :wasm_metadce),
    ExecutableProduct(prefix, "", :wasm_opt),
    ExecutableProduct(prefix, "", :wasm_merge),
    ExecutableProduct(prefix, "", :wasm_dis),
    ExecutableProduct(prefix, "", :wasm_as),
    ExecutableProduct(prefix, "", :wasm2asm),
    ExecutableProduct(prefix, "", :wasm_ctor_eval),
    LibraryProduct(prefix, "libbinaryen", :libbinaryen),
    ExecutableProduct(prefix, "", :wasm_shell),
    ExecutableProduct(prefix, "", :asm2wasm)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "Binaryen", sources, script, platforms, products, dependencies)

max_stars_repo_name: tshort/BinaryenBuilder

max_issues_repo_name: tshort/BinaryenBuilder

max_forks_repo_name: tshort/BinaryenBuilder


text:
module Heatmap

using Distances, Clustering
using Compose, Gadfly, StatsBase

PositionedGuide = Gadfly.Guide.PositionedGuide
top_guide_position = Gadfly.Guide.top_guide_position
right_guide_position = Gadfly.Guide.right_guide_position

"""
    HeatMapStatistic

A Gadfly `StatisticElement` type for heatmaps.
"""
struct HeatMapStatistic <: Gadfly.StatisticElement
   matrix::String
   metric::Metric
end

"""
    heatmap(; matrix="raw", metric=Euclidean())

Construct a `heatmap` type. `matrix` can be "raw" or "distance". `metric` as in package `Distances`.
"""
heatmap(; matrix::String="raw", metric::Metric=Euclidean()) = HeatMapStatistic(matrix, metric)

Gadfly.Stat.input_aesthetics(stat::HeatMapStatistic) =  [:z, :x, :y]
Gadfly.Stat.output_aesthetics(stat::HeatMapStatistic) = [:xmin, :xmax, :ymin, :ymax]
Gadfly.Stat.default_scales(stat::HeatMapStatistic) = [Scale.z_func(), Scale.x_discrete(), Scale.y_discrete(), Scale.color_continuous()]

function Gadfly.Stat.apply_statistic(stat::HeatMapStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    n, p = size(aes.z)
    if stat.matrix=="distance"
        aes.y = repeat(aes.x, inner=n)
        aes.x = repeat(aes.x, outer=n)
        dist = pairwise(stat.metric, aes.z, dims=1)
        aes.color_key_title = replace(string(typeof(stat.metric), "\n","distance"), "Distances."=>"")
    else
        aes.x = repeat(aes.x, outer=p)
        aes.y = repeat(aes.y, inner=n)
        dist = aes.z
        aes.color_key_title = ""
    end
    Stat.apply_statistic(Stat.rectbin(), scales, coord, aes)
    color_scale = get(scales, :color, Scale.color_continuous)
    Scale.apply_scale(color_scale, [aes], Gadfly.Data(color=vec(dist)))
end

function branches(hc::Hclust, useheight::Bool=true)
    order = StatsBase.indexmap(hc.order)
    nodepos = Dict(-i => (float(order[i]), 0.0) for i in hc.order)

    branches1 = Vector{NTuple{2, Float64}}[]
    branches2 = Vector{NTuple{2, Float64}}[]
    ypos = 0.0

    for i in 1:size(hc.merges, 1)
        x1, y1 = nodepos[hc.merges[i, 1]]
        x2, y2 = nodepos[hc.merges[i, 2]]

        xpos = (x1 + x2) / 2
        h = useheight ? hc.heights[i] : 1.0
        ypos = max(y1, y2) + h

        nodepos[i] = (xpos, ypos)
        push!(branches1, [(x1,y1), (x1,ypos), (x2,ypos), (x2,y2)])
        push!(branches2, [(y1,x1), (ypos,x1), (ypos,x2), (y2,x2)])
    end

    return (branches1, branches2, ypos)
end

"""
    Dendrogram

A Gadfly `GuideElement` type for marginal dendrograms.
"""
struct Dendrogram <: Gadfly.GuideElement
    dendrogram::String
    metric::Metric
end

"""
    dendrogram(; dendrogram="both", metric=Euclidean())

Construct a  dendrogram guide. `dendrogram` can be "row", "column" or "both". `metric` as in package `Distances`.
Note that currently `metric` is independent of `heatmap(matrix="distance", metric= )` hence both need to be specified if plotting a distance matrix heatmap.
"""
dendrogram(; dendrogram::String="both", metric::Metric=Euclidean()) = Dendrogram(dendrogram, metric)

function Gadfly.Guide.render(guide::Dendrogram, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    adict = Dict{String, Vector{Int}}("row"=>[1], "column"=>[2], "both"=>[1,2])
    dims = adict[guide.dendrogram]
    n, p = size(aes.z)

    pgv = PositionedGuide[]
    if in(1, dims)
        hc = hclust(pairwise(guide.metric, aes.z, dims=1))
        branches1, _, ypos = branches(hc)
        ctx = context(units=UnitBox(0.5, ypos, n, -ypos, bottompad=4px), minheight=25)
        compose!(ctx, line(branches1), stroke("black"), linewidth(0.3pt))
        push!(pgv, PositionedGuide([ctx], 0, top_guide_position))
    end
    if in(2, dims)
        d = 2
        if occursin("distance", aes.color_key_title)
            p, d = n, 1
        end
        hc = hclust(pairwise(guide.metric, aes.z, dims=d))
        _, branches2, ypos = branches(hc)
        ctx2 = context(units=UnitBox(0, p+0.5,  ypos, -p, leftpad=4px), minwidth=25)
        compose!(ctx2, line(branches2), stroke("black"), linewidth(0.3pt))
        push!(pgv, PositionedGuide([ctx2], 0, right_guide_position))
    end
    return pgv
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ b4949654-ff30-11ea-225c-2ffc47a35a31
using Pkg, DrWatson

# ╔═╡ b494da42-ff30-11ea-25d2-37857c757c9f
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
	using PlutoUI
end

# ╔═╡ b0e1bc4a-ff2f-11ea-1dca-f78fefb58b8f
md"## Clip-06-16s.jl"

# ╔═╡ b495756a-ff30-11ea-378a-236707f95773
begin
	N = 100
	df = DataFrame(
		:h0 => rand(Normal(10,2 ), N),
		:treatment => vcat(zeros(Int, Int(N/2)), ones(Int, Int(N/2)))
	);
	df.fungus = [rand(Binomial(1, 0.5 - 0.4 * df[i, :treatment]), 1)[1] for i in 1:N]
	df.h1 = [df[i, :h0] + rand(Normal(5 - 3 * df[i, :fungus]), 1)[1] for i in 1:N]
end;

# ╔═╡ 221f1e8a-00f6-11eb-3504-a91f4659a8a6
with_terminal() do
	precis(df; io=String)
end

# ╔═╡ b4a2fe56-ff30-11ea-255b-bbe592530dac
stan6_7 = "
data {
  int <lower=1> N;
  vector[N] h0;
  vector[N] h1;
  vector[N] treatment;
  vector[N] fungus;
}
parameters{
  real a;
  real bt;
  real bf;
  real<lower=0> sigma;
}
model {
  vector[N] mu;
  vector[N] p;
  a ~ lognormal(0, 0.2);
  bt ~ normal(0, 0.5);
  bf ~ normal(0, 0.5);
  sigma ~ exponential(1);
  for ( i in 1:N ) {
    p[i] = a + bt*treatment[i] + bf*fungus[i];
    mu[i] = h0[i] * p[i];
  }
  h1 ~ normal(mu, sigma);
}
";

# ╔═╡ b4a44e3c-ff30-11ea-0868-93935a280415
begin
	m6_7s = SampleModel("m6.7s", stan6_7)
	m6_7_data = Dict(
	  :N => nrow(df),
	  :h0 => df[:, :h0],
	  :h1 => df[:, :h1],
	  :fungus => df[:, :fungus],
	  :treatment => df[:, :treatment]
	)
	rc6_7s = stan_sample(m6_7s; data=m6_7_data)
	success(rc6_7s) && (post6_7s_df = read_samples(m6_7s, :dataframe))
end;

# ╔═╡ b4a98408-ff30-11ea-2918-9b5d39692695
success(rc6_7s) && Particles(post6_7s_df)

# ╔═╡ b4b7702a-ff30-11ea-169a-159df3080e8b
success(rc6_7s) && (Text(precis(post6_7s_df; io=String)))

# ╔═╡ b4b8301e-ff30-11ea-2a3c-a11e3745616e
if success(rc6_7s)
	(part6_7s, fig6_7s) = plot_model_coef([m6_7s], [:a, :bt, :bf])
	fig6_7s
end

# ╔═╡ b4c852fa-ff30-11ea-18d0-47dfe321ba8a
md"## End of clip-06-16s.jl"

# ╔═╡ Cell order:
# ╟─b0e1bc4a-ff2f-11ea-1dca-f78fefb58b8f
# ╠═b4949654-ff30-11ea-225c-2ffc47a35a31
# ╠═b494da42-ff30-11ea-25d2-37857c757c9f
# ╠═b495756a-ff30-11ea-378a-236707f95773
# ╠═221f1e8a-00f6-11eb-3504-a91f4659a8a6
# ╠═b4a2fe56-ff30-11ea-255b-bbe592530dac
# ╠═b4a44e3c-ff30-11ea-0868-93935a280415
# ╠═b4a98408-ff30-11ea-2918-9b5d39692695
# ╠═b4b7702a-ff30-11ea-169a-159df3080e8b
# ╠═b4b8301e-ff30-11ea-2a3c-a11e3745616e
# ╟─b4c852fa-ff30-11ea-18d0-47dfe321ba8a


text:
"""
    compute_bin_averages(O::Array{Float64, 1}, n::Int)

Computes the bin averages, used in computing the autocorrelation time.

Arguments
* `O`       : Array representing an observable, e.g. E or M
* `n`       : Number of observations per bin
Returns
* `bin_avgs`: Bin averages
"""
function compute_bin_averages(O::Array{Float64,1}, n::Int)
    N = length(O)
    N_B = Int(floor(N / n))
    bin_avgs = Array{Float64,1}(undef, N_B)

    for k = 1:N_B
        a = (k - 1) * n + 1
        b = k * n
        bin_avgs[k] = mean(O[a:b])
    end

    return bin_avgs
end

"""
    compute_autocorrelation_time(model::Union{Ising, Potts}, ns::Array{Int, 1})

Computes the integrated autocorrelation time using the binning method. The autocorrelation
time is the value at which the output sequence (either energy `:e` or magnetization `:m` in
the return dictionary) asymptotes at.

Arguments
* `model`: A Potts or Ising type model
* `ns`   : Array of bin sizes
Returns
* Dictionary with keys `[:e, :m]` with bin method values to be plotted against `ns`
"""
function compute_autocorrelation_time(model::Union{Ising,Potts}, ns::Array{Int,1})
    V = length(model.σ)
    e = model.observables.E ./ V
    m = abs.(model.observables.M) ./ V
    N = length(e)
    σ²_e = mean(e .^ 2) - mean(e)^2
    σ²_m = mean(m .^ 2) - mean(m)^2
    τ_int_e = []
    τ_int_m = []

    for n in ns
        N_B = floor(Int, N / n)
        e_bin_avgs = compute_bin_averages(e, n)
        m_bin_avgs = compute_bin_averages(m, n)
        ε²_e = sum((e_bin_avgs .- mean(e_bin_avgs)) .^ 2) / (N_B * (N_B - 1))
        ε²_m = sum((m_bin_avgs .- mean(m_bin_avgs)) .^ 2) / (N_B * (N_B - 1))
        push!(τ_int_e, ε²_e * N / 2σ²_e)
        push!(τ_int_m, ε²_m * N / 2σ²_m)
    end

    return Dict(:e => τ_int_e, :m => τ_int_m)
end

"""
    compute_autocorrelation_time(model::XY, ns::Array{Int, 1})

Computes the integrated autocorrelation time using the binning method. The autocorrelation
time is the value at which the output sequence (either energy `:e` or magnetization `:m` in
the return dictionary) asymptotes at.

Arguments
* `model`: An XY type model
* `ns`   : Array of bin sizes
Returns
* Dictionary with keys `[:e, :mx, my]` with bin method values to be plotted against `ns`
"""
function compute_autocorrelation_time(model::XY, ns::Array{Int,1})
    V = length(model.σ)
    e = model.observables.E ./ V
    mx = abs.(model.observables.Mx) ./ V
    my = abs.(model.observables.My) ./ V
    N = length(e)
    σ²_e = mean(e .^ 2) - mean(e)^2
    σ²_mx = mean(mx .^ 2) - mean(mx)^2
    σ²_my = mean(my .^ 2) - mean(my)^2
    τ_int_e = []
    τ_int_mx = []
    τ_int_my = []

    for n in ns
        N_B = floor(Int, N / n)
        e_bin_avgs = compute_bin_averages(e, n)
        mx_bin_avgs = compute_bin_averages(mx, n)
        my_bin_avgs = compute_bin_averages(my, n)
        ε²_e = sum((e_bin_avgs .- mean(e_bin_avgs)) .^ 2) / (N_B * (N_B - 1))
        ε²_mx = sum((mx_bin_avgs .- mean(mx_bin_avgs)) .^ 2) / (N_B * (N_B - 1))
        ε²_my = sum((my_bin_avgs .- mean(my_bin_avgs)) .^ 2) / (N_B * (N_B - 1))
        push!(τ_int_e, ε²_e * N / 2σ²_e)
        push!(τ_int_mx, ε²_mx * N / 2σ²_mx)
        push!(τ_int_my, ε²_my * N / 2σ²_my)
    end

    return Dict(:e => τ_int_e, :mx => τ_int_mx, :my => τ_int_my)
end

"""
    sim_dict_to_df(sims::Dict)

Converts a dictionary where the values are models into a DataFrame object.

Arguments
* `sims`: Dictionary where the values are simulations
Returns
* `df`  : DataFrame containing the simulation data
"""
function sim_dict_to_df(sims::Dict; include_magnetization::Bool = true)
    if include_magnetization
        df = DataFrame(
            L = Int[],
            β = Float64[],
            e = Float64[],
            c = Float64[],
            m = Float64[],
            χ = Float64[],
            U = Float64[],
        )

        for sim in values(sims)
            push!(
                df,
                [
                    sim.params.L,
                    sim.params.β,
                    sim.observables.statistics[:e],
                    sim.observables.statistics[:c],
                    sim.observables.statistics[:m],
                    sim.observables.statistics[:χ],
                    sim.observables.statistics[:U],
                ],
            )
        end
    else
        df = DataFrame(L = Int[], β = Float64[], e = Float64[], c = Float64[])

        for sim in values(sims)
            push!(
                df,
                [
                    sim.params.L,
                    sim.params.β,
                    sim.observables.statistics[:e],
                    sim.observables.statistics[:c],
                ],
            )
        end
    end

    sort!(df, [:L, :β])

    return df
end


text:
module Vcov

using StatsBase
using DataFrames
using LinearAlgebra
using FixedEffects
using Combinatorics

crossmodelmatrix(x::RegressionModel) = cholesky!(modelmatrix(x)' * modelmatrix(x))
##############################################################################
##
## Any type used for standard errors must define the following methods:
##
##############################################################################
materialize(df::AbstractDataFrame, v::CovarianceEstimator) = v
completecases(df::AbstractDataFrame, ::CovarianceEstimator) = trues(size(df, 1))
S_hat(x::RegressionModel, ::CovarianceEstimator) = error("S_hat not defined for this type")
df_FStat(x::RegressionModel, ::CovarianceEstimator, hasintercept::Bool) = dof_residual(x) - hasintercept

include("utils.jl")
include("vcovsimple.jl")
include("vcovrobust.jl")
include("vcovcluster.jl")

end


text:
## Exercise 2-2
## Repeating my advice from the previous chapter, whenever you learn a new feature, you should try it out in interactive mode and make errors on purpose to see what goes wrong.

## 1. We’ve seen that n = 42 is legal. What about 42 = n?
println("Ans 1:")
# 42 = n  # Error: syntax: invalid assignment location "42"

## 2. How about x = y = 1?
println("Ans 2:")
x = y = 1  # works fine x and y get assigned 1

## 3. In some languages every statement ends with a semi-colon, ;. What happens if you put a semi-colon at the end of a Julia statement?
println("Ans 3:")
x = 17;  # when in interactive environment it does not print the value, in script no change.

## 4. What if you put a period at the end of a statement?
println("Ans 4:")
x + 6.  # working fine
# x + y.  # Error syntax incomplete OR Load Error type Int64 has no field println

## In math notation you can multiply x and y like this: x y. What happens if you try that in Julia? What about 5x?
println("Ans 5:")
x = 3
y = 4
println(5x)  # works, 'yey
# println(xy)  # error UndefVarError, as xy will be different from x * y
# println(x y)  # error syntax, missing `,` or `)`
println(x * y)

println("End.")


text:
include("aerodynamic_torque.jl")
include("magnetic_torque.jl")
include("sensors.jl")

using LinearAlgebra
using SatelliteToolbox
#=引数候補
i_m = [1, 1, 1]
B = [1, 1, 1]
rho = 1
v = 1
e_r = [1, 0, 0]
sun_vecs = [-0.5 -0.5 1.]
v_ecef_pu = [1. 1. 1.; 1. 1. 1.]
attqua1 = SatelliteToolbox.Quaternion(1.0/sqrt(2.0), 1.0/sqrt(2.0), 0.0, 0.0)
attqua2 = SatelliteToolbox.Quaternion(1.0/sqrt(2.0), 0.0, 1.0/sqrt(2.0), 0.0)
dt = 5
=#
"""
torque, ss_out, ms_out, acce_out, gyro_out = internal_model(i_m, B, rho, v, e_r, sun_vecs, v_ecef_pu, attqua1, attqua2, dt)

# Argments
- `i_m`：各磁気トルカに流れる電流
- `B`：磁束密度ベクトル＠SCSF
- `rho`：大気密度
- `v`：大気に対する衛星の速度
- `e_r`：衛星の進行方向の単位ベクトル@SCSF（body）
- `sun_vecs`：太陽方向ベクトル@SCSF (body)
- `v_ecef_pu`：今回ループと1回前ループにおける衛星速度ベクトル@ECEF (3*2 1行目が1回前、2行目が今回)
- `attqua1`：2回前ループにおける衛星姿勢クォータニオン
- `attqua2`：1回前ループにおける衛星姿勢クォータニオン
- `dt`：ループの時間間隔

# Returns
- `torqe`：衛星にかかるトルク@SCSF
- `ss_out`：太陽センサが出力する電流値
- `ms_out`：磁気センサが出力する磁束密度データ（I2Cシリアル出力想定）
- `acce_out`：ジャイロセンサが出力する加速度データ@ECEF（I2Cシリアル出力想定）
- `gyro_out`：ジャイロセンサが出力する角速度データ@SCSF（I2Cシリアル出力想定）
"""
function internal_model(i_m, B, rho, v, e_r, sun_vecs, v_ecef_pu, attqua1, attqua2, dt)
	#=
	設定パラメータ
	=#
	C_d = 1 #衛星のサーフェースの抵抗係数
	sur_num = 6 #衛星の面数
	p_c = [150 0 0; -150 0 0; 0 150 0; 0 -150 0; 0 0 150; 0 0 -150] #各面の圧力中心
	e_a = [1 0 0; -1 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 -1] #各面の法線ベクトル
	a = [90000, 90000, 90000, 90000, 90000, 90000] #衛星の各面積
	torque = zeros(3)

	#トルク計算
	aer_tor = aerodynamic_torque(sur_num, C_d, rho, v, a, e_r, e_a, p_c)
	mag_tor = magnetic_torque(i_m, B)
	torque = aer_tor + mag_tor

	#センサ出力
	ss_out = sun_sensor(sun_vecs, "z-")
	ms_out = mag_sensor(B)
	acce_out, gyro_out = gyro_sensor(v_ecef_pu, attqua1, attqua2, dt)

	return torque, ss_out, ms_out, acce_out, gyro_out
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
export TypedBool
"""
    abstract TypedBool

typed bools, [`True`](@ref) and [`False`](@ref), can guarantee type stability in
cases where constant propogation is not working for Bools.

```jldoctest
julia> using Keys

julia> Bool(False())
false

julia> Bool(True())
true

julia> TypedBool(false)
False()

julia> TypedBool(true)
True()

julia> True() & True() & False()
False()

julia> False() & False() & True()
False()

julia> True() | True() | False()
True()

julia> False() | False() | True()
True()

```
"""
abstract type TypedBool end
@inline TypedBool(b::Bool) =
    if b
        True()
    else
        False()
    end

export True
"a [`TypedBool`](@ref)"
struct True <: TypedBool end

export False
"a [`TypedBool`](@ref)"
struct False <: TypedBool end

Bool(::True) = true
Bool(::False) = false

(&)(::False, ::False) = False()
(&)(::False, ::True) = False()
(&)(::True, ::False) = False()
(&)(::True, ::True) = True()

(|)(::False, ::False) = False()
(|)(::False, ::True) = True()
(|)(::True, ::False) = True()
(|)(::True, ::True) = True()

export not
"""
    not(x)

[`TypedBool`](@ref) aware version of `!`.

```jldoctest
julia> using Keys

julia> not(True())
False()

julia> not(False())
True()
```
"""
not(::False) = True()
not(::True) = False()

"""
    if_else(switch, new, old)

[`TypedBool`](@ref) aware version of `ifelse`.

```julia
julia> using Keys

julia> if_else(true, 1, 0)
1

julia> if_else(True(), 1, 0)
1

julia> if_else(False(), 1, 0)
0
```
"""
if_else(b::Bool, new, old) = ifelse(b, new, old)  # generic fallback
if_else(::True, new, old) = new
if_else(::False, new, old) = old

export same_type
"""
    same_type(a, b)

Check whether `a` and `b` are the same type; return a [`TypedBool`](@ref).

```jldoctest
julia> using Keys

julia> same_type(1, 2)
True()

julia> same_type(1, 2.0)
False()
```
"""
same_type(a::T, b::T) where T = True()
same_type(a, b) = False()