text:
Flux.@nograd Base.similar
Base.copyto!(A::UpperTriangular{<:Any, <:Flux.CuArrays.GPUArrays.AbstractGPUArray}, B::UpperTriangular{<:Any, <:Flux.CuArrays.GPUArrays.AbstractGPUArray}) = LinearAlgebra.triu!(copyto!(parent(A), parent(B)))

module GPU

using LinearAlgebra
using Flux.CuArrays: CuMatrix
using Flux.CuArrays.CUBLAS: CublasFloat
import Base: +

function LinearAlgebra.diag(x::Adjoint{T,<:CuMatrix{T}}, i::Int) where T<:CublasFloat
  diag(parent(x), -i)
end

function +(x::UnitUpperTriangular{T,<:Adjoint{T,<:CuMatrix{T}}}, y::UniformScaling) where T<:CublasFloat
  LowerTriangular(parent(parent(x)) + y)'
end

end # module


text:
export ISSCorner

mutable struct ISSCorner{T<:AbstractFloat} <: Environment
  worldAABBmin::Vector{T}
  worldAABBmax::Vector{T}
  keepin_zones::Vector
  keepout_zones::Vector
  obstacle_set::Vector
end

function ISSCorner{T}() where T
  keepin_zones = Vector{HyperRectangle}()
  vars = matread(joinpath(abspath(joinpath(dirname(Base.find_package("GuSTO")), "..")), "src", "environment","iss_corner.mat"))
  for zone in vars["keepin_zones"]
    push!(keepin_zones,
      HyperRectangle(Vec3f0(zone["corner1"][:]),Vec3f0(zone["corner2"][:]-zone["corner1"][:])))
  end

  keepout_zones = Vector{GeometryTypes.GeometryPrimitive}()
  for zone in vars["keepout_zones"]
    push!(keepout_zones,
      HyperRectangle(Vec3f0(zone["corner1"][:]),Vec3f0(zone["corner2"][:]-zone["corner1"][:])))
  end

  obstacle_set = Vector{GeometryTypes.GeometryPrimitive}()

  worldAABBmin = Inf*ones(T,3)
  worldAABBmax = -Inf*ones(T,3)
  for zone in (keepin_zones..., keepout_zones...)
    corner1, corner2 = zone.origin, zone.origin + zone.widths
    zone_min = [min(corner1[i], corner2[i]) for i in 1:3]
    worldAABBmin[worldAABBmin .> zone_min] = zone_min[worldAABBmin .> zone_min]
    zone_max = [max(corner1[i], corner2[i]) for i in 1:3]
    worldAABBmax[worldAABBmax .< zone_max] = zone_max[worldAABBmax .< zone_max]
  end

  return ISSCorner{T}(worldAABBmin, worldAABBmax, keepin_zones, keepout_zones, obstacle_set)
end
ISSCorner(::Type{T} = Float64; kwargs...) where {T} = ISSCorner{T}(; kwargs...)

function update_aabb!(env::ISSCorner)
  for zone in env.keepin_zones
    corner1, corner2 = zone.origin, zone.origin + zone.widths
    zone_min = [min(corner1[i], corner2[i]) for i in 1:3]
    env.worldAABBmin[env.worldAABBmin .> zone_min] = zone_min[env.worldAABBmin .> zone_min]
    zone_max = [max(corner1[i], corner2[i]) for i in 1:3]
    env.worldAABBmax[env.worldAABBmax .< zone_max] = zone_max[env.worldAABBmax .< zone_max]
  end
  warn("Overriding collision world in env()")
end

function add_obstacles!(env::ISSCorner)
  vars = matread(joinpath(abspath(joinpath(dirname(Base.find_package("GuSTO")), "..")), "src", "environment","iss_corner.mat"))
  for zone in vars["rectangles"]
    push!(env.obstacle_set,
      HyperRectangle(Vec3f0(zone["corner1"][:]),Vec3f0(zone["corner2"][:]-zone["corner1"][:])))
  end

  for zone in vars["spheres"]
    push!(env.obstacle_set,
      HyperSphere(Point3f0(zone["center"][:]), Float32(zone["radius"])))
  end
end

ext: jl


text:
using DocumentationGenerator, DataStructures, Pkg

packages = DocumentationGenerator.installable_on_version(VERSION)

max_packages = if isempty(ARGS)
    length(packages)
else
    parse(Int, ARGS[1])
end
# max_packages = 60

processes = parse(Int, get(ENV, "NUM_PKG_PROCESSES", "8"))

docspath = normpath(joinpath(@__DIR__, ".."))

tstart = time()
DocumentationGenerator.build_documentations(packages[1:max_packages], processes = processes, basepath = docspath)
@info("DONE in $(round(time() - tstart, digits=1)) seconds")

##
@info("Building Website...")

distpath = joinpath(@__DIR__, "..", "dist")
distdocpath = joinpath(@__DIR__, "..", "dist", "packages")

rm(distpath, recursive=true, force=true)
isdir(distpath) || mkpath(distpath)
isdir(distdocpath) || mkdir(distdocpath)

printheader(io::IO, name; level = 1) = print(io, "#"^level, " ", name)

function printlink(io::IO, name, link)
    print(io, "[", name, "](", link, ")")
end

## get packages and versions
built_packages = OrderedDict{String, Vector{String}}()
for pkg in sort!(readdir(joinpath(docspath, "build")), by=x->lowercase(x))
    pkgpath = joinpath(docspath, "build", pkg)
    isdir(pkgpath) || continue

    cp(pkgpath, joinpath(distdocpath, pkg); force=true)
    open(joinpath(distdocpath, pkg, "versions.js"), "w") do io
        println(io, "var DOC_VERSIONS = [")
        vers = sort!(readdir(pkgpath), by=x->VersionNumber(x))
        for ver in vers
            isdir(joinpath(pkgpath, ver)) || continue
            push!(get!(built_packages, pkg, String[]), ver)
            println(io, "\"", ver, "\",")
            open(joinpath(distdocpath, pkg, ver, "siteinfo.js"), "w") do io
                println(io, """
                var DOCUMENTER_CURRENT_VERSION = "$(ver)";
                """)
            end
        end
        print(io, "];")
    end
    open(joinpath(distdocpath, pkg, "index.html"), "w") do io
        println(io, """
        <head>
        <meta http-equiv="refresh" content="0; url=$(last(built_packages[pkg]))" />
        </head>
        """)
    end
end

include("generate_json.jl")


text:
using FunctionOperators, LinearAlgebra, Test

@testset "Helpers.jl" begin
    Op₁ = FunctionOperator{Float64}(x -> x, (1,), (1,))
    Op₂ = FunctionOperator{Int64}(x -> x, (1,), (1,))
    Op₃ = FunctionOperator{Float64}(x -> x, (1,), (1,))
    Op₄ = FunctionOperator{Float64}(x -> x, (1,), (2,))
    buffer = Array{Float64}(undef, 1)
    @testset "eltype" begin
        @test eltype(Op₁) == Float64
        @test eltype(Op₁ * Op₃) == Float64
    end
    @testset "assertions" begin
        @testset "TypeError" begin
            @test_throws TypeError Op₁ * Op₂
            @test_throws TypeError Op₁ + Op₂
            @test_throws TypeError Op₁ - Op₂
            @test_throws TypeError (Op₁ + Op₁) * Op₂
            @test_throws TypeError (Op₁ + Op₁) + Op₂
            @test_throws TypeError (Op₁ + Op₁) - Op₂
            @test_throws TypeError Op₁ * (Op₂ + Op₂)
            @test_throws TypeError Op₁ + (Op₂ + Op₂)
            @test_throws TypeError Op₁ - (Op₂ + Op₂)
            @test_throws TypeError (Op₁ + Op₁) * (Op₂ + Op₂)
            @test_throws TypeError (Op₁ + Op₁) + (Op₂ + Op₂)
            @test_throws TypeError (Op₁ + Op₁) - (Op₂ + Op₂)
            @test_throws TypeError mul!(Array{Int64}(undef, 1), Op₁, [1.])
            @test_throws TypeError mul!(Array{Int64}(undef, 1), Op₁ * Op₂, [1.])
            @test_throws TypeError mul!(buffer, Op₁, [1])
            @test_throws TypeError mul!(buffer, Op₁ * Op₂, [1])
        end
        @testset "DimensionError" begin
            @test_throws DimensionMismatch Op₁ * [1.,2.]
            @test_throws DimensionMismatch Op₁ * Op₃ * [1.,2.]
            @test_throws DimensionMismatch (Op₁ + Op₃) * [1.,2.]
            @test_throws DimensionMismatch (Op₁ - Op₃) * [1.,2.]
            @test_throws DimensionMismatch Op₃ * Op₄
            @test_throws DimensionMismatch (Op₄ + Op₃)
            @test_throws DimensionMismatch (Op₄ - Op₃)
            @test_throws DimensionMismatch (Op₄ + I)
            @test_throws DimensionMismatch (Op₄ - I)
            @test_throws DimensionMismatch (Op₄ + 3I)
            @test_throws DimensionMismatch (Op₄ - 3I)
            @test_throws DimensionMismatch mul!(buffer, Op₁, [1., 2.])
            @test_throws DimensionMismatch mul!(buffer, Op₁ * Op₃, [1.,2.])
            @test_throws DimensionMismatch mul!(Array{Float64}(undef, 2), Op₁, [1.])
            @test_throws DimensionMismatch mul!(Array{Float64}(undef, 2), Op₁ * Op₃, [1.])
        end
    end
end;


text:
using Test
using GraphPlot # For plotting tests
using Neo4j
push!(ENV,"DFG_USE_CGDFG"=>"true")
using DistributedFactorGraphs
using Pkg
using Dates
using TimeZones
using SHA
using UUIDs

## To run the IIF tests, you need a local Neo4j with user/pass neo4j:test
# To run a Docker image
# NOTE: that Neo4j.jl doesn't currently support > Neo4j 3.x
# Install: docker pull neo4j:3.5.6
# Run: docker run -d --publish=7474:7474 --publish=7687:7687 --env NEO4J_AUTH=neo4j/test neo4j:3.5.6
##

# If you want to enable debugging logging (very verbose!)
# using Logging
# logger = SimpleLogger(stdout, Logging.Debug)
# global_logger(logger)

include("test_defVariable.jl")

include("testBlocks.jl")

@testset "Test generated ==" begin
    include("compareTests.jl")
end

@testset "Testing LightDFG.FactorGraphs functions" begin
    include("LightFactorGraphsTests.jl")
end

if get(ENV, "DO_CGDFG_TESTS", "") == "true"
    apis = [
        LightDFG,
        CloudGraphsDFG,
        ]
    @warn "TEST: Removing all data for user 'testUserId'!"
    clearUser!!(CloudGraphsDFG(userId="testUserId"))

else
    apis = [
        LightDFG,
        ]
end

for api in apis
    @testset "Testing Driver: $(api)" begin
        @info "Testing Driver: $(api)"
        global testDFGAPI = api
        include("interfaceTests.jl")
    end
end

# Test special cases
@testset "Plotting Tests" begin
    include("plottingTest.jl")
end

@testset "Data Store Tests" begin
    include("consol_DataEntryBlobTests.jl")
end

@testset "LightDFG subtype tests" begin
    for type in [(var=DFGVariableSummary, fac=DFGFactorSummary), (var=SkeletonDFGVariable,fac=SkeletonDFGFactor)]
        @testset "$(type.var) and $(type.fac) tests" begin
            @info "Testing $(type.var) and $(type.fac)"
            global VARTYPE = type.var
            global FACTYPE = type.fac
            include("LightDFGSummaryTypes.jl")
        end
    end
end

if get(ENV, "IIF_TEST", "") == "true"

    # Switch to our upstream test branch.
    Pkg.add(PackageSpec(name="IncrementalInference", rev="upstream/dfg_integration_test"))
    @info "------------------------------------------------------------------------"
    @info "These tests are using IncrementalInference to do additional driver tests"
    @info "------------------------------------------------------------------------"

    using IncrementalInference

    apis = Vector{AbstractDFG}()
    push!(apis, LightDFG(solverParams=SolverParams(), userId="testUserId"))
    get(ENV, "DO_CGDFG_TESTS", "false") == "true" && push!(apis, CloudGraphsDFG(solverParams=SolverParams(), userId="testUserId"))

    for api in apis
        @testset "Testing Driver: $(typeof(api))" begin
            @info "Testing Driver: $(api)"
            global dfg = deepcopy(api)
            include("iifInterfaceTests.jl")
        end

        @testset "FileDFG Testing Driver: $(typeof(api))" begin
            @info "FileDFG Testing Driver: $(typeof(api))"
            global dfg = deepcopy(api)
            include("fileDFGTests.jl")
        end
    end

    @testset "IIF Compare Tests" begin
        #run a copy of compare tests from IIF
        include("iifCompareTests.jl")
    end

    @testset "CGStructure Tests for CGDFG" begin
        # Run the CGStructure tests
        include("CGStructureTests.jl")
    end

    # Simple graph solving test
    @testset "Simple graph solving test" begin
        # This is just to validate we're not going to blow up downstream.
        apis = [
            # GraphsDFG{SolverParams}(),
            LightDFG(solverParams=SolverParams(), userId="testUserId"),
            CloudGraphsDFG(solverParams=SolverParams(), userId="testUserId")
            ]
        for api in apis
            @info "Running simple solver test: $(typeof(api))"
            global dfg = deepcopy(api)
            include("solveTest.jl")
        end
    end
else
    @warn "Skipping IncrementalInference driver tests"
end

struct NotImplementedDFG{T} <: AbstractDFG{T} end

@testset "No Interface tests" begin
    dfg = NotImplementedDFG{NoSolverParams}()
    v1 = SkeletonDFGVariable(:v1)
    f1 = SkeletonDFGFactor(:f1)

    @test_throws ErrorException exists(dfg, v1)
    @test_throws ErrorException exists(dfg, f1)

    @test_throws ErrorException exists(dfg, :s)
    @test_throws ErrorException addVariable!(dfg, v1)

    @test_throws ErrorException getVariable(dfg, :a)
    @test_throws ErrorException getFactor(dfg, :a)
    @test_throws ErrorException updateVariable!(dfg, v1)
    @test_throws ErrorException updateFactor!(dfg, f1)

    @test_throws ErrorException deleteVariable!(dfg, :a)
    @test_throws ErrorException deleteFactor!(dfg, :a)
    @test_throws ErrorException getVariables(dfg)
    @test_throws ErrorException getFactors(dfg)
    @test_throws ErrorException isConnected(dfg)
    @test_throws ErrorException getNeighbors(dfg, v1)
    @test_throws ErrorException getNeighbors(dfg, :a)

    @test_throws ErrorException _getDuplicatedEmptyDFG(dfg)

    @test_throws ErrorException isVariable(dfg, :a)
    @test_throws ErrorException isFactor(dfg, :a)

end


text:
function getDivGrad(n1,n2,n3)
	# D = getDivGrad(n1,n2,n3)
	# builds 3D divergence operator
	D1 = kron(sparse(1.0I,n3*n2,n3*n2),ddx(n1))
	D2 = kron(sparse(1.0I,n3,n3),kron(ddx(n2),sparse(1.0I,n1,n1)))
	D3 = kron(ddx(n3),sparse(1.0I,n2*n1,n2*n1))

	Div = [D1 D2 D3]
	return Div*Div'
end

function ddx(n)
# generate 1D finite difference on staggered grid
I,J,V = SparseArrays.spdiagm_internal(0=>-ones(n),1=>ones(n) )
return sparse(I,J,V,n,n+1)
	# return d = spdiags(ones(n)*[-1 1],[0,1],n,n+1)
end

function spdiags(B,d,m,n)
# A = spdiags(B,d,m,n)
# creates a sparse matrix from its diagonals
	d = d[:]
	p = length(d)
	len = zeros(p+1,1)
	for k = 1:p
		len[k+1] = round.(Int,len[k]+length(max(1,1-d[k]):min(m,n-d[k])))
	end
	a = zeros(round.(Int,len[p+1]),3)
	for k = 1:p
		# Append new d[k]-th diagonal to compact form
		i = max(1,1-d[k]):min(m,n-d[k])
		a[(round.(Int,len[k])+1):round.(Int,len[k+1]),:] = [i i+d[k] B[i+(m>=n)*d[k],k]]
	end

	A = sparse(round.(Int,a[:,1]),round.(Int,a[:,2]),a[:,3],m,n)
	return A
end


text:
#=
Popis
  Funkcia vypocita radius grafu.

Parametre
  g : datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozena vlastnost vrchola :: String
  components : parameter pre prepinanie graf/max_komponent :: Bool

Navratova hodnota
  : radius maximalneho komponentu :: Float64
  g_radius : radius grafu :: Float64

Poznamky
  Funkcia podla hodnoty parametra components vypocita hodnotu
  pre cely graf alebo pre najvacsi komponent.
=#

function graph_radius(g::GenericGraph; dist_key::AbstractString="weight", components::Bool=false)

  if components
    comp_radius = Dict{Array{ExVertex,1}, Float64}()
    comp_length = Int64[]

    for comp in connected_components(g)
      comp_radius[comp] = minimum([vertex_eccentricity(v.index, g) for v in comp])
      push!(comp_length, length(comp))
    end

    max_comp = maximum(comp_length)

    return [r[2] for r=filter!(x -> length(x[1]) == max_comp, collect(zip(keys(comp_radius), values(comp_radius))))]
  else
    g_radius = [vertex_eccentricity(v.index, g) for v in vertices(g)]

    return minimum(g_radius)
  end

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
abstract type Bound{T} end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, bound::Bound{T}) where {T}
    summary(io, bound)
end

@inline getT(bound::Bound{T}) where T = T

@inline g(bound::Bound{T}) where T = zero(T)

@inline ∂g∂pos(bound::Bound{T}) where T = szeros(T, 6)
@inline ∂g∂vel(bound::Bound{T}) where T = szeros(T, 6)
@inline schurf(bound::Bound{T}) where T = szeros(T, 6)
@inline schurD(bound::Bound{T}) where T = szeros(T, 6, 6)

@inline additionalforce(bound::Bound{T}) where T = szeros(T, 6)


text:
#=
    - Julia Version : 1.4.0
    - Author        : Alexander Bills and Abhinav Misalkar
    - Date          : 14th May 2020
    - Title         : Unit Converter
=#

struct unit_conversion

 ################################################### Design Input Unit Conversion ###################################################

    fc_pos_sp_cap           # mAh/g --> Ah/g
    fc_neg_sp_cap           # mAh/g --> Ah/g
    rev_pos_sp_cap          # mAh/g --> Ah/g
    rev_neg_sp_cap          # mAh/g --> Ah/g

    pos_th                  # microns --> cm
    neg_th                  # microns --> cm

    pos_extra_l             # mm --> cm
    neg_extra_l             # mm --> cm

    pos_CC_th               # microns --> cm
    neg_CC_th               # microns --> cm

    pos_tab_t               # microsn --> cm
    neg_tab_t               # microsn --> cm

    pos_tab_width           # mm --> cm
    neg_tab_width           # mm --> cm

    sep_th                  # microns --> cm

 ################################################### Design Output Unit Conversion ###################################################

    energy_cell                  #  Wh   -->   kWh

    mass_pos_AM                  #  g    -->   kg
    mass_neg_AM                  #  g    -->   kg

    pos_CC_area                  #  cm2  -->   m2
    neg_CC_area                  #  cm2  -->   m2

    sep_area                     #  cm2  -->   m2

    electrolyte_volm             #  mL   -->   L

    total_coated_area_cathode    #  cm2  -->   m2

 ################################################### Cost Input Unit Conversion ###################################################

    units_mfg               # unit to million units

end

mult = unit_conversion(
                            1000.0,
                            1000.0,
                            1000.0,
                            1000.0,

                            10000.0,
                            10000.0,

                            10.0,
                            10.0,

                            10000.0,
                            10000.0,

                            10000.0,
                            10000.0,

                            10.0,
                            10.0,

                            10000.0,

                            #---------------------------------------------------------#

                            1000.0,
                            1000.0,
                            1000.0,
                            10000.0,
                            10000.0,
                            10000.0,
                            1000.0,
                            10000.0,

                            #---------------------------------------------------------#

                            1/1000000.0
                        )

function converter(initial_value,mult)
   return initial_value ./ mult
end

function convert_all(cell_general, cost, mult)
   cell_general.cathode.fc_sp_cap =
      converter(cell_general.cathode.fc_sp_cap, mult.fc_pos_sp_cap)
   cell_general.anode.fc_sp_cap =
      converter(cell_general.anode.fc_sp_cap, mult.fc_neg_sp_cap)
   cell_general.cathode.rev_sp_cap =
      converter(cell_general.cathode.rev_sp_cap, mult.rev_pos_sp_cap)
   cell_general.anode.rev_sp_cap =
      converter(cell_general.anode.rev_sp_cap, mult.rev_neg_sp_cap)

   cell_general.cathode.th = converter(cell_general.cathode.th, mult.pos_th)
   cell_general.anode.th = converter(cell_general.anode.th, mult.neg_th)

   cell_general.cathode.extra_length =
      converter(cell_general.cathode.extra_length, mult.pos_extra_l)
   cell_general.anode.extra_length =
      converter(cell_general.anode.extra_length, mult.neg_extra_l)

   cell_general.cathode.CC_th =
      converter(cell_general.cathode.CC_th, mult.pos_CC_th)
   cell_general.anode.CC_th =
      converter(cell_general.anode.CC_th, mult.neg_CC_th)

   cell_general.cathode.tab_th =
      converter(cell_general.cathode.tab_th, mult.pos_tab_t)
   cell_general.anode.tab_th =
      converter(cell_general.anode.tab_th, mult.neg_tab_t)

   cell_general.cathode.tab_width =
      converter(cell_general.cathode.tab_width, mult.pos_tab_width)
   cell_general.anode.tab_width =
      converter(cell_general.anode.tab_width, mult.neg_tab_width)

   cell_general.sep_th = converter(cell_general.sep_th, mult.sep_th)

   cost.general_costs.no_units_mfg =
      converter(cost.general_costs.no_units_mfg, mult.units_mfg)
   return cell_general, cost
end

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
abstract type SceneOverlay end

@with_kw mutable struct TextParams
    size::Int = 12
    color::Colorant = colorant"white"
    x::Int = 20
    y_start::Int = 20
    y_jump::Int = round(Int, size*1.2)
end
function drawtext(text::AbstractString, y::Int, rendermodel::RenderModel, t::TextParams; coordinate_system::Symbol=:camera_pixels)
    add_instruction!(rendermodel, render_text, (text, t.x, y, t.size, t.color), coordinate_system=coordinate_system)
    y + t.y_jump
end

"""
    TextOverlay

Displays some text at the desired location.
The coordinates and size units are in pixels by default.
The option `coordinate_system` allows to use different units.

# Fields
- `text::Vector{String}`
- `color::Colorant = colorant"white"`
- `font_size::Int = 10`
- `pos::VecE2 = VecE2(10, font_size)`
- `line_spacing::Float64 = 1.5` multiple of `font_size`
- `coordinate_system::Symbol=:camera_pixels`
"""
@with_kw mutable struct TextOverlay <: SceneOverlay
    text::Vector{String}
    color::Colorant = colorant"white"
    font_size::Int = 10 # [pix]
    pos::VecE2 = VecE2(10, font_size)
    line_spacing::Float64 = 1.5 # multiple of font_size
    coordinate_system::Symbol=:camera_pixels
end
function add_renderable!(rendermodel::RenderModel, overlay::TextOverlay)
    x = overlay.pos.x
    y = overlay.pos.y
    y_jump = overlay.line_spacing*overlay.font_size
    for line in overlay.text
        add_instruction!(rendermodel, render_text, (line, x, y, overlay.font_size, overlay.color), coordinate_system=overlay.coordinate_system)
        y += y_jump
    end
    rendermodel
end

# method for new interface
function add_renderable!(rendermodel::RenderModel, overlay::TextOverlay, scene)
    x = overlay.pos.x
    y = overlay.pos.y
    y_jump = overlay.line_spacing*overlay.font_size
    for line in overlay.text
        add_instruction!(rendermodel, render_text, (line, x, y, overlay.font_size, overlay.color), coordinate_system=overlay.coordinate_system)
        y += y_jump
    end
    rendermodel
end

"""
    Overwash
overlay that renders a plain color on the whole canvas.
"""
mutable struct Overwash <: SceneOverlay
    color::Colorant
end
function add_renderable!(rendermodel::RenderModel, overlay::Overwash)
    add_instruction!(rendermodel, render_paint, (overlay.color,))
    rendermodel
end

"""
    HistogramOverlay

Display a bar at the specified position `pos`, the bar is of size `width`, `height` and is filled up to a given proportion of its height.
The fill proportion is set using `val`, it should be a number between 0 and 1. If it is 0, the bar is not filled, if it is 1 it is filled to the top.

# Fields
- `pos::VecE2{Float64} = VecE2(0.,0.)`
- `coordinate_system::Symbol = :scene`
- `label::String = "histogram"`
- `val::Float64 = 0.5` should be between 0 and 1
- `width::Float64 = 2.`
- `height::Float64 = 5.`
- `fill_color::Colorant = colorant"blue"`
- `line_color::Colorant = colorant"white"`
- `font_size::Int64 = 15 # [pix]`
- `label_pos::VecSE2{Float64} = pos + VecSE2(0., -height/2)` position of the label

"""
@with_kw mutable struct HistogramOverlay <: SceneOverlay
    pos::VecE2{Float64} = VecE2(10.,10.)
    coordinate_system::Symbol = :camera_pixels
    label::String = "histogram"
    val::Float64 = 0.5 # should be between 0 and 1
    width::Float64 = 2.
    height::Float64 = 5.
    fill_color::Colorant = colorant"blue"
    line_color::Colorant = colorant"white"
    font_size::Int64 = 15 # [pix]
    label_pos::VecSE2{Float64} = pos + VecSE2(0., -height/2)
end

function AutoViz.add_renderable!(rendermodel::RenderModel, overlay::HistogramOverlay)
    # render value
    add_instruction!(rendermodel, render_rect, (overlay.pos.x, overlay.pos.y, overlay.width, overlay.val*overlay.height,overlay.fill_color, true, false), coordinate_system=overlay.coordinate_system)
    # render histogram outline
    add_instruction!(rendermodel, render_rect, (overlay.pos.x, overlay.pos.y, overlay.width, overlay.height, overlay.line_color), coordinate_system=overlay.coordinate_system)
     # label
    add_instruction!(rendermodel, render_text, (overlay.label, overlay.label_pos.x, overlay.label_pos.y, overlay.font_size, overlay.line_color), coordinate_system=overlay.coordinate_system)
end

"""
    IDOverlay

Display the ID on top of each entity in a scene.
The text can be customized with the `color::Colorant` (default=white) and `font_size::Int64` (default=15) keywords.
The position of the ID can be adjusted using `x_off::Float64` and `y_off::Float64` (in camera coordinates).

# Fields
- `color::Colorant = colorant"white"`
- `font_size::Int = 15`
- `x_off::Float64 = 0.`
- `y_off::Float64 = 0.`
"""
@with_kw mutable struct IDOverlay <: SceneOverlay
    color::Colorant = colorant"white"
    font_size::Int = 15
    x_off::Float64 = 0.
    y_off::Float64 = 0.
end

function AutoViz.add_renderable!(rendermodel::RenderModel, overlay::IDOverlay, scene::Frame{Entity{S,D,I}}, env::E) where {S,D,I,E}
    font_size = overlay.font_size
    for veh in scene
        add_instruction!(rendermodel, render_text, ("$(veh.id)", veh.state.posG.x + overlay.x_off, veh.state.posG.y + overlay.y_off, font_size, overlay.color), coordinate_system=:scene)
    end
    return rendermodel
end

@with_kw mutable struct LineToCenterlineOverlay <: SceneOverlay
    target_id::Int # if -1 does it for all
    line_width::Float64 = 0.5
    color::Colorant = colorant"blue"
end
function add_renderable!(rendermodel::RenderModel, overlay::LineToCenterlineOverlay, scene::Frame{Entity{S,D,I}}, roadway::Any) where {S,D,I}

    if overlay.target_id < 0
        target_inds = 1:length(scene)
    else
        target_inds = overlay.target_id:overlay.target_id
    end

    for ind in target_inds
        veh = scene[ind]
        footpoint = get_footpoint(veh)
        add_instruction!(rendermodel, render_line_segment,
            (veh.state.posG.x, veh.state.posG.y, footpoint.x, footpoint.y, overlay.color, overlay.line_width))
    end

    rendermodel
end

@with_kw mutable struct LineToFrontOverlay <: SceneOverlay
    target_id::Int # if -1 does it for all
    line_width::Float64 = 0.5
    color::Colorant = colorant"blue"
end
function add_renderable!(rendermodel::RenderModel, overlay::LineToFrontOverlay, scene::Frame{Entity{S,D,I}}, roadway::Roadway) where {S,D,I}

    if overlay.target_id < 0
        target_inds = 1:length(scene)
    else
        target_inds = overlay.target_id:overlay.target_id
    end

    for ind in target_inds
        veh = scene[ind]
        veh_ind_front = get_neighbor_fore_along_lane(scene, ind, roadway).ind
        if veh_ind_front != nothing
            v2 = scene[veh_ind_front]
            add_instruction!(rendermodel, render_line_segment,
                (veh.state.posG.x, veh.state.posG.y, v2.state.posG.x, v2.state.posG.y, overlay.color, overlay.line_width))
        end
    end

    rendermodel
end

"""
    BlinkerOverlay
Displays a circle on one of the top corner of a vehicle to symbolize a blinker.
fields:
- on: turn the blinker on
- right: blinker on the top right corner, if false, blinker on the left
- veh: the vehicle for which to display the blinker
- color: the color of the blinker
- size: the size of the blinker
"""
@with_kw struct BlinkerOverlay <: SceneOverlay
    on::Bool = false
    right::Bool = true
    veh::Vehicle = Vehicle(VehicleState(), VehicleDef(), 0)
    color::Colorant = colorant"0xFFEF00" # yellow
    size::Float64 = 0.3
end

function add_renderable!(rendermodel::RenderModel, overlay::BlinkerOverlay)
    if !overlay.on
        return nothing
    end
    if overlay.right
        pos = get_front(overlay.veh) + polar(overlay.veh.def.width/2, overlay.veh.state.posG.θ - pi/2)
    else
        pos = get_front(overlay.veh) + polar(overlay.veh.def.width/2, overlay.veh.state.posG.θ + pi/2)
    end
    add_instruction!(rendermodel, render_circle, (pos.x, pos.y, overlay.size, overlay.color))
end

"""
    CarFollowingStatsOverlay

Displays statistics about the front neighbor of the car of id `target_id`.

# Constructor

`CarFollowingStatsOverlay(;target_id, verbosity=1, color=colorant"white", font_size=10)`
"""
@with_kw mutable struct CarFollowingStatsOverlay <: SceneOverlay
    target_id::Int
    verbosity::Int = 1
    color::Colorant = colorant"white"
    font_size::Int = 10
end
function add_renderable!(rendermodel::RenderModel, overlay::CarFollowingStatsOverlay, scene::Frame{Entity{S,D,I}}, roadway::Roadway) where {S,D,I}

    font_size = overlay.font_size
    text_y = font_size
    text_y_jump = round(Int, font_size*1.2)
    fmt_txt = @sprintf("id %d", overlay.target_id)
    add_instruction!( rendermodel, render_text, (fmt_txt, 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
        text_y += text_y_jump

    veh_index = findfirst(overlay.target_id, scene)
    if veh_index != nothing
        veh = scene[veh_index]

        if overlay.verbosity ≥ 2
            add_instruction!( rendermodel, render_text, ("posG: " * string(veh.state.posG), 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
            text_y += text_y_jump
            add_instruction!( rendermodel, render_text, ("posF: " * string(veh.state.posF), 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
            text_y += text_y_jump
        end
        fmt_txt = @sprintf("speed: %0.3f", veh.state.v)
        add_instruction!( rendermodel, render_text, (fmt_txt, 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
        text_y += text_y_jump

        foreinfo = get_neighbor_fore_along_lane(scene, veh_index, roadway; max_distance_fore=Inf)
        if foreinfo.ind != nothing
            v2 = scene[foreinfo.ind]
            rel_speed = v2.state.v - veh.state.v
            fmt_txt = @sprintf("Δv = %10.3f m/s", rel_speed)
            add_instruction!( rendermodel, render_text, (fmt_txt, 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
            text_y += text_y_jump
            fmt_txt = @sprintf("Δs = %10.3f m/s", foreinfo.Δs)
            add_instruction!( rendermodel, render_text, (fmt_txt, 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
            text_y += text_y_jump

            if overlay.verbosity ≥ 2
                add_instruction!( rendermodel, render_text, ("posG: " * string(v2.state.posG), 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
                text_y += text_y_jump
                add_instruction!( rendermodel, render_text, ("posF: " * string(v2.state.posF), 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
                text_y += text_y_jump
                fmt_txt = @sprintf("speed: %.3f", v2.state.v)
                add_instruction!( rendermodel, render_text, (fmt_txt, 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
                text_y += text_y_jump
            end
        else
            fmt_txt = @sprintf("no front vehicle")
            add_instruction!( rendermodel, render_text, (fmt_txt, 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
        end
    else
        fmt_txt = @sprintf("vehicle %d not found", overlay.target_id)
        add_instruction!( rendermodel, render_text, (fmt_txt, 10, text_y, font_size, overlay.color), coordinate_system=:camera_pixels)
    end

    rendermodel
end
"""
    NeighborsOverlay

Draws a line between a vehicle and its neighbors. The neighbors are linked with different colors depending on their lanes.
This overlay needs to be wrapped as a RenderableOverlay which needs a scene and a roadway to perform the calculation of the neighbors.

# Fields

-  `target_id::Int`
-  `color_L::Colorant = colorant"blue"`
-  `color_M::Colorant = colorant"green`
-  `color_R::Colorant = colorant"red"`
-  `line_width::Float64 = 0.5`
-  `textparams::TextParams = TextParams()`
"""
@with_kw mutable struct NeighborsOverlay <: SceneOverlay
    target_id::Int
    color_L::Colorant = colorant"blue"
    color_M::Colorant = colorant"green"
    color_R::Colorant = colorant"red"
    line_width::Float64 = 0.5 # [m]
    textparams::TextParams = TextParams()
end
function add_renderable!(rendermodel::RenderModel, overlay::NeighborsOverlay, scene::Frame{Entity{S,D,I}}, roadway::Roadway) where {S,D,I}

    textparams = overlay.textparams
    yₒ = textparams.y_start
    Δy = textparams.y_jump

    vehicle_index = findfirst(overlay.target_id, scene)
    if vehicle_index != nothing

        veh_ego = scene[vehicle_index]
        t = veh_ego.state.posF.t
        ϕ = veh_ego.state.posF.ϕ
        v = veh_ego.state.v
        len_ego = veh_ego.def.length

        fore_L = get_neighbor_fore_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
        if fore_L.ind != nothing
            veh_oth = scene[fore_L.ind]
            A = get_front(veh_ego)
            B = get_rear(veh_oth)
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_L, overlay.line_width))
            drawtext(@sprintf("d fore left:   %10.3f", fore_L.Δs), yₒ + 0*Δy, rendermodel, textparams)
        end

        fore_M = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
        if fore_M.ind != nothing
            veh_oth = scene[fore_M.ind]
            A = get_front(veh_ego)
            B = get_rear(veh_oth)
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_M, overlay.line_width))
            drawtext(@sprintf("d fore middle: %10.3f", fore_M.Δs), yₒ + 1*Δy, rendermodel, textparams)
        end

        fore_R = get_neighbor_fore_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
        if fore_R.ind != nothing
            veh_oth = scene[fore_R.ind]
            A = get_front(veh_ego)
            B = get_rear(veh_oth)
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_R, overlay.line_width))
            drawtext(@sprintf("d fore right:  %10.3f", fore_R.Δs), yₒ + 2*Δy, rendermodel, textparams)
        end

        rear_L = get_neighbor_rear_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointFront(), VehicleTargetPointRear())
        if rear_L.ind != nothing
            veh_oth = scene[rear_L.ind]
            A = get_rear(veh_ego)
            B = get_front(veh_oth)
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_L, overlay.line_width))
            drawtext(@sprintf("d rear left:   %10.3f", rear_L.Δs), yₒ + 3*Δy, rendermodel, textparams)
        end

        rear_M = get_neighbor_rear_along_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointFront(), VehicleTargetPointRear())
        if rear_M.ind != nothing
            veh_oth = scene[rear_M.ind]
            A = get_rear(veh_ego)
            B = get_front(veh_oth)
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_M, overlay.line_width))
            drawtext(@sprintf("d rear middle: %10.3f", rear_M.Δs), yₒ + 4*Δy, rendermodel, textparams)
        end

        rear_R = get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointFront(), VehicleTargetPointRear())
        if rear_R.ind != nothing
            veh_oth = scene[rear_R.ind]
            A = get_rear(veh_ego)
            B = get_front(veh_oth)
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_R, overlay.line_width))
            drawtext(@sprintf("d rear right:  %10.3f", rear_R.Δs), yₒ + 5*Δy, rendermodel, textparams)
        end
    end

    rendermodel
end

mutable struct MarkerDistOverlay <: SceneOverlay
    target_id::Int
    textparams::TextParams
    rec::SceneRecord
    function MarkerDistOverlay(target_id::Int;
        textparams::TextParams=TextParams(),
        )

        new(target_id, textparams, SceneRecord(1, 0.1))
    end
end

function add_renderable!(rendermodel::RenderModel, overlay::MarkerDistOverlay, scene::Frame{Entity{S,D,I}}, roadway::Roadway) where {S,D,I}

    textparams = overlay.textparams
    yₒ = textparams.y_start
    Δy = textparams.y_jump

    update!(overlay.rec, scene)

    vehicle_index = findfirst(scene, overlay.target_id)
    if vehicle_index != nothing
        veh_ego = scene[vehicle_index]
        drawtext(@sprintf("lane offset:       %10.3f", veh_ego.state.posF.t), yₒ + 0*Δy, rendermodel, textparams)
        drawtext(@sprintf("markerdist left:   %10.3f", convert(Float64, get(MARKERDIST_LEFT, overlay.rec, roadway, vehicle_index))), yₒ + 1*Δy, rendermodel, textparams)
        drawtext(@sprintf("markerdist right:  %10.3f", convert(Float64, get(MARKERDIST_RIGHT, overlay.rec, roadway, vehicle_index))), yₒ + 2*Δy, rendermodel, textparams)
    end

    rendermodel
end

"""
    RenderableOverlay

Decorator which allows to use `SceneOverlay` objects together with the method
    render([Renderables])

This is required primarily for allowing backward compatibility with overlays
that use the old rendering interface.

usage:  `RenderableOverlay(o::Overlay, scene::Frame, roadway::Roadway)`
"""
struct RenderableOverlay{O,S,D,I} <: Renderable where {O<:SceneOverlay,S,D,I}
    overlay::O
    scene::Union{Nothing, Frame{Entity{S,D,I}}}
    roadway::Union{Nothing, Roadway}
end
RenderableOverlay(overlay::O) where {O<:SceneOverlay} = RenderableOverlay(overlay, nothing, nothing)

function add_renderable!(rm::RenderModel, ro::RenderableOverlay{O,S,D,I}) where {O<:SceneOverlay,S,D,I}
    if (ro.scene === nothing) && (ro.roadway === nothing)
        # These are the overlays that could also be used without scene and roadway,
        # they do not require the RenderableOverlay wrapper
        add_renderable!(rm, ro.overlay)
    else
        # These overlays rely on the `RenderableOverlay` wrapper to provide scene and roadway
        @assert (ro.scene !== nothing) && (ro.roadway !== nothing)
        add_renderable!(rm, ro.overlay, ro.scene, ro.roadway)
    end
end

function add_renderable!(rm::RenderModel, iterable::Union{Array{<:RenderableOverlay},Tuple{<:RenderableOverlay}})
    for ro in iterable add_renderable!(rm, ro) end
end
isrenderable(::Type{Array{<:RenderableOverlay}}) = true
isrenderable(::Type{Tuple{<:RenderableOverlay}}) = true


text:
# FASTA Reader
# ============

struct Reader{S <: TranscodingStream} <: BioGenerics.IO.AbstractReader
    state::State{S}
    index::Union{Index, Nothing}
end

"""
    FASTA.Reader(input::IO; index = nothing)

Create a data reader of the FASTA file format.

# Arguments
* `input`: data source
* `index=nothing`: filepath to a random access index (currently *fai* is supported)
"""
function Reader(input::IO; index = nothing)
    if isa(index, AbstractString)
        index = Index(index)
    elseif index != nothing
        throw(ArgumentError("index must be a filepath or nothing"))
    end
    if !(input isa TranscodingStream)
        stream = TranscodingStreams.NoopStream(input)
        return Reader(State(stream, 1, 1, false), index)
    else
        return Reader(State(input, 1, 1, false), index)
    end
end

function Base.eltype(::Type{<:Reader})
    return Record
end

function BioGenerics.IO.stream(reader::Reader)
    return reader.state.stream
end

function Base.read!(rdr::Reader, rec::Record)
    cs, ln, f = readrecord!(rdr.state.stream, rec, (rdr.state.state, rdr.state.linenum))
    rdr.state.state = cs
    rdr.state.linenum = ln
    rdr.state.filled = f
    if !f
        cs == 0 && throw(EOFError())
        throw(ArgumentError("malformed FASTA file"))
    end
    return rec
end

function Base.close(reader::Reader)
    if reader.state.stream isa IO
        close(reader.state.stream)
    end
    return nothing
end

function Base.getindex(reader::Reader, name::AbstractString)
    if reader.index == nothing
        throw(ArgumentError("no index attached"))
    end
    #seekrecord(reader.state.stream, reader.index, name)
    seekrecord(reader.state.stream, reader.index, name) # TODO: reader.index may not be type stable.
    #reader.state.cs = file_machine.start_state
    #reader.state.finished = false
    #return read(reader)
    record = Record()
    cs, linenum, found = readrecord!(TranscodingStreams.NoopStream(reader.state.stream), record, (1, 1))
    @assert cs ≥ 0 && found
    return record
end

function index!(record::Record)
    stream = TranscodingStreams.NoopStream(IOBuffer(record.data))
    cs, linenum, found = readrecord!(stream, record, (1, 1))
    if !found || !allspace(stream)
        throw(ArgumentError("invalid FASTA record"))
    end
    return record
end

function allspace(stream)
    while !eof(stream)
        if !isspace(read(stream, Char))
            return false
        end
    end
    return true
end


text:
module Lorenzo

# package code goes here
include("scraper.jl")
export test
end # module


text:
#Unscented Kalman Inversion: specific structures and function definitions

"""
    Unscented{FT<:AbstractFloat, IT<:Int} <: Process

An unscented Kalman Inversion process
"""
mutable struct Unscented{FT <: AbstractFloat, IT <: Int} <: Process
    "a vector of arrays of size `N_parameters` containing the mean of the parameters (in each `uki` iteration a new array of mean is added)"
    u_mean::Vector{Array{FT, 1}}
    "a vector of arrays of size (`N_parameters x N_parameters`) containing the covariance of the parameters (in each `uki` iteration a new array of `cov` is added)"
    uu_cov::Vector{Array{FT, 2}}
    "a vector of arrays of size `N_y` containing the predicted observation (in each `uki` iteration a new array of predicted observation is added)"
    obs_pred::Vector{Array{FT, 1}}
    "weights in UKI"
    c_weights::Array{FT, 1}
    mean_weights::Array{FT, 1}
    cov_weights::Array{FT, 1}
    "covariance of the artificial evolution error"
    Σ_ω::Array{FT, 2}
    "covariance of the artificial observation error"
    Σ_ν_scale::FT
    "regularization parameter"
    α_reg::FT
    "regularization vector"
    r::Array{FT, 1}
    "update frequency"
    update_freq::IT
    "current iteration number"
    iter::IT
end

# outer constructors
function EnsembleKalmanProcess(
    obs_mean::Array{FT, 1},
    obs_noise_cov::Array{FT, 2},
    process::Unscented{FT, IT};
    Δt = FT(1),
) where {FT <: AbstractFloat, IT <: Int}

    #initial parameters stored as columns
    init_params = [DataContainer(update_ensemble_prediction!(process), data_are_columns = true)]
    # Number of parameter
    N_u = length(process.u_mean[1])
    # ensemble size
    N_ens = 2N_u + 1 #stored with data as columns
    #store for model evaluations
    g = []
    # error store
    err = FT[]
    # timestep store
    Δt = Array([Δt])

    EnsembleKalmanProcess{FT, IT, Unscented}(init_params, obs_mean, obs_noise_cov, N_ens, g, err, Δt, process)
end

"""
EnsembleKalmanProcess Constructor
u0_mean::Array{FT} : prior mean
uu0_cov::Array{FT, 2} : prior covariance
obs_mean::Array{FT,1} : observation
obs_noise_cov::Array{FT, 2} : observation error covariance
α_reg::FT : regularization parameter toward `u0` (0 < `α_reg` ≤ 1), default should be 1, without regulariazion
update_freq::IT : set to 0 when the inverse problem is not identifiable,
                  namely the inverse problem has multiple solutions,
                  the covariance matrix will represent only the sensitivity of the parameters,
                  instead of posterior covariance information;
                  set to 1 (or anything > 0) when the inverse problem is identifiable, and
                  the covariance matrix will converge to a good approximation of the
                  posterior covariance with an uninformative prior.

"""
function Unscented(
    u0_mean::Array{FT, 1},
    uu0_cov::Array{FT, 2},
    α_reg::FT,
    update_freq::IT;
    modified_uscented_transform::Bool = true,
    κ::FT = 0.0,
    β::FT = 2.0,
) where {FT <: AbstractFloat, IT <: Int}

    N_u = size(u0_mean, 1)
    # ensemble size
    N_ens = 2 * N_u + 1

    c_weights = zeros(FT, N_u)
    mean_weights = zeros(FT, N_ens)
    cov_weights = zeros(FT, N_ens)

    # todo parameters λ, α, β

    α = min(sqrt(4 / (N_u + κ)), 1.0)
    λ = α^2 * (N_u + κ) - N_u

    c_weights[1:N_u] .= sqrt(N_u + λ)
    mean_weights[1] = λ / (N_u + λ)
    mean_weights[2:N_ens] .= 1 / (2 * (N_u + λ))
    cov_weights[1] = λ / (N_u + λ) + 1 - α^2 + β
    cov_weights[2:N_ens] .= 1 / (2 * (N_u + λ))

    if modified_uscented_transform
        mean_weights[1] = 1.0
        mean_weights[2:N_ens] .= 0.0
    end

    u_mean = Array{FT, 1}[]  # array of Array{FT, 2}'s
    push!(u_mean, u0_mean) # insert parameters at end of array (in this case just 1st entry)
    uu_cov = Array{FT, 2}[] # array of Array{FT, 2}'s
    push!(uu_cov, uu0_cov) # insert parameters at end of array (in this case just 1st entry)

    obs_pred = Array{FT, 1}[]  # array of Array{FT, 2}'s
    err = FT[]

    Σ_ω = (2 - α_reg^2) * uu0_cov
    Σ_ν_scale = 2.0

    r = u0_mean
    iter = 0

    Unscented(
        u_mean,
        uu_cov,
        obs_pred,
        c_weights,
        mean_weights,
        cov_weights,
        Σ_ω,
        Σ_ν_scale,
        α_reg,
        r,
        update_freq,
        iter,
    )
end

"""
    construct_sigma_ensemble(process::Unscented, x_mean::Array{FT}, x_cov::Array{FT, 2}) where {FT <: AbstractFloat, IT <: Int}

Construct the sigma ensemble based on the mean `x_mean` and covariance `x_cov`.
"""
function construct_sigma_ensemble(
    process::Unscented,
    x_mean::Array{FT},
    x_cov::Array{FT, 2},
) where {FT <: AbstractFloat, IT <: Int}

    N_x = size(x_mean, 1)

    c_weights = process.c_weights

    chol_xx_cov = cholesky(Hermitian(x_cov)).L

    x = zeros(FT, N_x, 2 * N_x + 1)
    x[:, 1] = x_mean
    for i in 1:N_x
        x[:, i + 1] = x_mean + c_weights[i] * chol_xx_cov[:, i]
        x[:, i + 1 + N_x] = x_mean - c_weights[i] * chol_xx_cov[:, i]
    end

    return x
end

"""
construct_mean `x_mean` from ensemble `x`.
"""
function construct_mean(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    x::Array{FT, 2},
) where {FT <: AbstractFloat, IT <: Int}
    N_x, N_ens = size(x)

    @assert(uki.N_ens == N_ens)

    x_mean = zeros(FT, N_x)

    mean_weights = uki.process.mean_weights

    for i in 1:N_ens
        x_mean += mean_weights[i] * x[:, i]
    end

    return x_mean
end

"""
construct_cov `xx_cov` from ensemble `x` and mean `x_mean`.
"""
function construct_cov(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    x::Array{FT, 2},
    x_mean::Array{FT},
) where {FT <: AbstractFloat, IT <: Int}
    N_ens, N_x = uki.N_ens, size(x_mean, 1)

    cov_weights = uki.process.cov_weights

    xx_cov = zeros(FT, N_x, N_x)

    for i in 1:N_ens
        xx_cov .+= cov_weights[i] * (x[:, i] - x_mean) * (x[:, i] - x_mean)'
    end

    return xx_cov
end

"""
construct_cov `xy_cov` from ensemble x and mean `x_mean`, ensemble `obs_mean` and mean `y_mean`.
"""
function construct_cov(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    x::Array{FT, 2},
    x_mean::Array{FT},
    obs_mean::Array{FT, 2},
    y_mean::Array{FT},
) where {FT <: AbstractFloat, IT <: Int, P <: Process}
    N_ens, N_x, N_y = uki.N_ens, size(x_mean, 1), size(y_mean, 1)

    cov_weights = uki.process.cov_weights

    xy_cov = zeros(FT, N_x, N_y)

    for i in 1:N_ens
        xy_cov .+= cov_weights[i] * (x[:, i] - x_mean) * (obs_mean[:, i] - y_mean)'
    end

    return xy_cov
end

"""
uki prediction step : generate sigma points
"""
function update_ensemble_prediction!(process::Unscented) where {FT <: AbstractFloat, IT <: Int}

    process.iter += 1
    # update evolution covariance matrix
    if process.update_freq > 0 && process.iter % process.update_freq == 0
        process.Σ_ω = (2 - process.α_reg^2) * process.uu_cov[end]
    end

    u_mean = process.u_mean[end]
    uu_cov = process.uu_cov[end]

    α_reg = process.α_reg
    r = process.r
    Σ_ω = process.Σ_ω

    N_u = length(process.u_mean[1])
    ############# Prediction step:

    u_p_mean = α_reg * u_mean + (1 - α_reg) * r
    uu_p_cov = α_reg^2 * uu_cov + Σ_ω

    ############ Generate sigma points
    u_p = construct_sigma_ensemble(process, u_p_mean, uu_p_cov)
    return u_p
end

"""
uki analysis step
g is the predicted observations  `Ny x N_ens` matrix
"""
function update_ensemble_analysis!(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    u_p::Array{FT, 2},
    g::Array{FT, 2},
) where {FT <: AbstractFloat, IT <: Int}

    obs_mean = uki.obs_mean
    Σ_ν = uki.process.Σ_ν_scale * uki.obs_noise_cov

    ############# Prediction step:

    u_p_mean = construct_mean(uki, u_p)
    uu_p_cov = construct_cov(uki, u_p, u_p_mean)

    ###########  Analysis step

    g_mean = construct_mean(uki, g)
    gg_cov = construct_cov(uki, g, g_mean) + Σ_ν
    ug_cov = construct_cov(uki, u_p, u_p_mean, g, g_mean)

    tmp = ug_cov / gg_cov

    u_mean = u_p_mean + tmp * (obs_mean - g_mean)
    uu_cov = uu_p_cov - tmp * ug_cov'

    ########### Save results
    push!(uki.process.obs_pred, g_mean) # N_ens x N_data
    push!(uki.process.u_mean, u_mean) # N_ens x N_params
    push!(uki.process.uu_cov, uu_cov) # N_ens x N_data

    push!(uki.g, DataContainer(g, data_are_columns = true))

    compute_error!(uki)
end

function update_ensemble!(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    g_in::Array{FT, 2},
) where {FT <: AbstractFloat, IT <: Int}
    #catch works when g_in non-square
    if !(size(g_in)[2] == uki.N_ens)
        throw(DimensionMismatch("ensemble size in EnsembleKalmanProcess and g_in do not match, try transposing or check ensemble size"))
    end

    u_p_old = get_u_final(uki)

    #perform analysis on the model runs
    update_ensemble_analysis!(uki, u_p_old, g_in)
    #perform new prediction output to model parameters u_p
    u_p = update_ensemble_prediction!(uki.process)

    push!(uki.u, DataContainer(u_p, data_are_columns = true))

    return u_p

end

function get_u_mean_final(uki::EnsembleKalmanProcess{FT, IT, Unscented}) where {FT <: AbstractFloat, IT <: Int}
    return uki.process.u_mean[end]
end

function compute_error!(uki::EnsembleKalmanProcess{FT, IT, Unscented}) where {FT <: AbstractFloat, IT <: Int}
    mean_g = uki.process.obs_pred[end]
    diff = uki.obs_mean - mean_g
    X = uki.obs_noise_cov \ diff # diff: column vector
    newerr = dot(diff, X)
    push!(uki.err, newerr)
end

function get_error(uki::EnsembleKalmanProcess{FT, IT, Unscented}) where {FT <: AbstractFloat, IT <: Int}
    return uki.err
end

function Gaussian_2d(
    u_mean::Array{FT, 1},
    uu_cov::Array{FT, 2},
    Nx::IT,
    Ny::IT;
    xx = nothing,
    yy = nothing,
) where {FT <: AbstractFloat, IT <: Int}
    # 2d Gaussian plot
    u_range = [min(5 * sqrt(uu_cov[1, 1]), 5); min(5 * sqrt(uu_cov[2, 2]), 5)]

    if xx === nothing
        xx = Array(LinRange(u_mean[1] - u_range[1], u_mean[1] + u_range[1], Nx))
    end
    if yy == nothing
        yy = Array(LinRange(u_mean[2] - u_range[2], u_mean[2] + u_range[2], Ny))
    end
    X, Y = repeat(xx, 1, Ny), Array(repeat(yy, 1, Nx)')
    Z = zeros(FT, Nx, Ny)

    det_uu_cov = det(uu_cov)

    for ix in 1:Nx
        for iy in 1:Ny
            Δxy = [xx[ix] - u_mean[1]; yy[iy] - u_mean[2]]
            Z[ix, iy] = exp(-0.5 * (Δxy' / uu_cov * Δxy)) / (2 * pi * sqrt(det_uu_cov))
        end
    end

    return xx, yy, Z

end


text:
const simplegaussstanmodel = "
data {
  int<lower=0>  N;
  real x[N];
}
parameters {
  real<lower=0> s;
  real m;
}
model {
  s ~ inv_gamma(2,3);
  m ~ normal(0,sqrt(s));
    x ~ normal(m,sqrt(s));
}
"


text:
using Test
include("MotifCommunities.jl")
include("Utilities.jl")

karate_edges = [[2,1],[3,1],[3,2],[4,1],[4,2],[4,3],[5,1],[6,1],[7,1],[7,5],[7,6],[8,1],[8,2],[8,3],[8,4],[9,1],
[9,3],[10,3],[11,1],[11,5],[1,6],[12,1],[13,1],[13,4],[14,1],[14,2],[14,3],[14,4],[17,6],[17,7],[18,1],[18,2],
[20,1],[20,2],[22,1],[22,2],[26,24],[26,25],[28,3],[28,24],[28,25],[29,3],[30,24],[30,27],[31,2],[31,9],[32,1],
[32,25],[32,26],[32,29],[33,3],[33,9],[33,15],[33,16],[33,19],[33,21],[33,23],[33,24],[33,30],[33,31],[33,32],[34,9],
[34,10],[34,14],[34,15],[34,16],[34,19],[34,20],[34,21],[34,23],[34,24],[34,27],[34,28],[34,29],[34,30],[34,31],[34,32],[34,33]]
karate_edges = reduce(vcat, [[i, reverse(i)] for i in karate_edges])
K = MatrixNetwork(map(x->x[1],karate_edges),map(x->x[2],karate_edges)) # Zachary's Karate Club

A = MatrixNetwork([1,1,2,2,3,3,3,4,4,4,5,5,5,6,6,7,7,8,8,8,9,9,10,10],
                  [2,3,1,3,1,2,4,3,5,8,4,6,7,5,7,5,6,4,9,10,8,10,8,9])

@testset "Modularity" begin
    mod(A::MatrixNetwork,
        C::Vector{Vector{Int64}};
        M::MatrixNetwork=MatrixNetwork([1],[2])) = motif_modularity(A,C;M=M)[1]

    @test mod(A, [[i for i = 1:10]]) == 0.
    @test mod(A, [[1,2,3,4], [5,6,7], [8,9,10]]) == 0.48958333333333337
    @test mod(A, [[1,2,3,4], [5,6,7], [8,9,10]]) ==
          mod(A, [[1,2,3], [4,5,6,7], [8,9,10]]) ==
          mod(A, [[1,2,3], [5,6,7], [4,8,9,10]])
    @test mod(A, [[1,2,3,4], [5,6,7], [8,9,10]]; M=clique(3)) == 0.8624989471979267
    @test mod(K, [[1,2,3,4,5,6,7,8,9,11,12,13,14,17,18,20,22],
                  [10,15,16,19,21,23,24,25,26,27,28,29,30,31,32,33,34]]) == 0.3568055321302074
    @test mod(K, [[1,2,3,4,5,6,7,8,9,11,12,13,14,17,18,20,22],
                  [10,15,16,19,21,23,24,25,26,27,28,29,30,31,32,33,34]]; M=path(2)) == 0.4814330690954068
end

@testset "Community Detection" begin
    @test_throws AssertionError louvain_motif(MatrixNetwork(Vector{Tuple{Int64,Int64}}(), 0)) == []
    @test louvain_motif(MatrixNetwork(Vector{Tuple{Int64,Int64}}(), 1)) == [[1]]
    @test louvain_motif(A) == [[1,2,3], [4,8,9,10], [5,6,7]]
    @test louvain_motif(A; M=clique(3)) == [[1,2,3], [4], [5,6,7], [8,9,10]]
    @test louvain_motif(K; M=path(2)) == [[1,2,3,4,5,6,7,8,9,11,12,13,14,17,18,20,22],
                                          [10,15,16,19,21,23,24,25,26,27,28,29,30,31,32,33,34]]

    D = MatrixNetwork([(1,2),(1,3),(1,4),(2,3),(2,4),(3,5),(3,10),(5,6),(5,8),(6,7),(6,8),(6,9),
                       (7,9),(8,9),(8,12),(9,11),(10,11),(10,13),(10,14),(11,12),(11,13),(12,14),(13,14)], 14)
    D = make_undirected(D)
    @test louvain_motif(D) == [[1,2,3,4],[5,6,7,8,9],[10,11,12,13,14]]
    @test louvain_motif(D; M=clique(3)) == [[1,2,3,4],[5,6,7,8,9],[10,11,13,14],[12]]
end


text:
function MORFE_integrate_rdyn_frc(analysis_name,
  zero_amplitude,harmonics_init,param_init,cont_param,time_integration_length,
  forward=true,MaxNumPoints=100,minstep=1e-8,maxstep=20.0,initstep=0.1,ncol=4.0,ntst=40.0,analysis_number=1)
    #
    # check for continuation direction
    if (forward)
      var = 0
    else
      var = 1
    end
    #
    # output folder
    analysis_folder = "./"*analysis_name*"/matcont_automatic"
    analysis_output = "./"*analysis_name*"/frc_"*string(analysis_number)*".mat"
    #
    nΩ = Int(size(harmonics_init)[1]/2)
    control_parameters = "mu"
    for i = 1:nΩ
      control_parameters *= ",beta"*string(i)
    end
    for i = 1:nΩ
      control_parameters *= ",w"*string(i)
    end
    #
    # initialise normal coordinates and harmonics amplitudes
    X0 = zeros(Float64,size(zero_amplitude)[1]+size(harmonics_init)[1])
    X0[1:size(zero_amplitude)[1]] = zero_amplitude
    X0[size(zero_amplitude)[1]+1:end] = harmonics_init
    #
    # inizialise control parameters
    param = ""
    param *= "mu = "*string(param_init[1])
    for i = 1:nΩ
      param *= ";beta"*string(i) * " = " * string(param_init[i+1])
    end
    for i = 1:nΩ
      param *= ";w"*string(i) * " = " * string(param_init[i+1+nΩ])
    end
    param *= ";"
    #
    # Matlab session
    ms = MSession();
    @mput analysis_folder
    @mput analysis_output
    @mput X0
    @mput time_integration_length
    @mput cont_param
    @mput ncol
    @mput ntst
    @mput MaxNumPoints
    @mput maxstep
    @mput minstep
    @mput initstep
    @mput var
    #
    eval_string(ms,"
    warning off;
    addpath(genpath('./MatCont7p3'))
    addpath(analysis_folder);
    ndofs = size(X0);
    ndofs = ndofs(1);
    init
    tfin = time_integration_length;
    ")
    #
    eval_string(ms,param)
    #
    eval_string(ms,"
    %%           Find periodic limit cycle for first continuation point
    ndofs = size(X0);
    ndofs = ndofs(1);
    hls = feval(@MORFEsystem);
    options=odeset('RelTol',1e-8);
    [t,y] = ode45(hls{2},[0,tfin],X0,options,"*control_parameters*");
    x1 = y(end,:);
    [vl,pk] = findpeaks(y(:,1));
    rel_err = (vl(end)-vl(end-1))/vl(end);
    period = t(pk(end))-t(pk(end-1));
    %
    %
    disp(\"-----------------------------------\")
    disp(\"| Time integration to find periodic limit cycle at first continuation point\")
    disp(\"| Relative error between last two consecutive peaks:\")
    disp(strcat(\"| \",num2str(rel_err)))
    disp(\"| Frequency (rad,s) found at supposed steady state:\")
    disp(strcat(\"| \",num2str(2*pi/period)))
    disp(\"| Consider enlarging the time_integration_length if the transient doesn't appear decayed\")
    disp(\"| Figure 100 shows the behaviour of the first normal coordinates in time\")
    figure(100);plot(t,y(:,1));
    xlabel('\$t\$','Interpreter','latex');
    ylabel('\$a_1 \$','Interpreter','latex');
    disp(\"-----------------------------------\")
    %
    %              integrate for a single period
    [t,y] = ode45(hls{2},[0 period],x1,options,"*control_parameters*");
    tolerance=1e-4;
    %
    %              use limit cycle as initial orbit for continuation
    [x0,v0]=initOrbLC(@MORFEsystem,t,y,["*control_parameters*"],[cont_param],ntst,ncol,tolerance);
    %
    %%           Continuation of Periodic orbits
    opt=contset;
    opt=contset(opt,'MaxNumPoints',MaxNumPoints);
    opt=contset(opt, 'InitStepsize' , initstep);
    opt=contset(opt,'MaxStepsize'  , maxstep);
    opt=contset(opt,'MinStepsize'  , minstep);
    opt=contset(opt,'Backward',var);
    opt=contset(opt,'FunTolerance', 1e-6);
    opt=contset(opt,'VarTolerance', 1e-6);
    opt=contset(opt,'ActiveParams', cont_param);
    %
    %              launch continuation
    [xlcc,vlcc,slcc,hlcc,flcc]=cont(@limitcycle,x0,v0,opt);
    %
    %              save output
    save(analysis_output,'xlcc','ncol','ntst')
    ")
    #
    return ms
    #
  end

function MORFE_integrate_rdyn_backbone(analysis_name,
  X0,time_integration_length,forward=true,MaxNumPoints=100,
  minstep=1e-8,maxstep=20.0,initstep=0.1,ncol=4.0,ntst=40.0,analysis_number=1)

  if (forward)
    var = 0
  else
    var = 1
  end

  analysis_folder = "./"*analysis_name*"/matcont_automatic"
  analysis_output = "./"*analysis_name*"/backbone_"*string(analysis_number)*".mat"

  ms = MSession();
  @mput analysis_folder
  @mput analysis_output
  @mput X0
  @mput time_integration_length
  @mput ncol
  @mput ntst
  @mput MaxNumPoints
  @mput maxstep
  @mput minstep
  @mput initstep
  @mput var
  eval_string(ms,"
  warning off
  addpath(genpath('./MatCont7p3'))
  addpath(analysis_folder)
  init
  mu = 0.0;
  ndofs = size(X0);
  ndofs = ndofs(1);
  tfin=time_integration_length;
  %
  hls=feval(@MORFEsystem);
  options=odeset('RelTol',1e-8);
  %
  [t,y]=ode45(hls{2},[0,tfin],X0,options,mu);
  %
  x1 = y(end,:);
  %
  [vl,pk]=findpeaks(y(:,1));
  rel_err = (vl(end)-vl(end-1))/vl(end);
  period = t(pk(end))-t(pk(end-1));
  %
  %
  disp(\"-----------------------------------\")
  disp(\"| Time integration to find periodic limit cycle at first continuation point\")
  disp(\"| Relative error between last two consecutive peaks:\")
  disp(strcat(\"| \",num2str(rel_err)))
  disp(\"| Frequency (rad,s) found at supposed steady state:\")
  disp(strcat(\"| \",num2str(2*pi/period)))
  disp(\"| Consider enlarging the time_integration_length if the transient doesn't appear decayed\")
  disp(\"| Figure 100 shows the behaviour of the first normal coordinates in time\")
  figure(100);plot(t,y(:,1));
  xlabel('\$t\$','Interpreter','latex');
  ylabel('\$a_1 \$','Interpreter','latex');
  disp(\"-----------------------------------\")
  %
  %
  [t,y] = ode45(hls{2},[0 period],x1,options,mu);
  %
  %
  active_pars=[1];
  tolerance=1e-4;
  [x0,v0]=initOrbLC(@MORFEsystem,...
  t,y,...
  [mu],active_pars,ntst,ncol,...
  tolerance);
  %
  opt=contset;
  opt=contset(opt,'MaxNumPoints',MaxNumPoints);
  opt=contset(opt, 'InitStepsize' , initstep);
  opt=contset(opt,'MaxStepsize'  , maxstep);
  opt=contset(opt,'MinStepsize'  , minstep);
  opt=contset(opt,'Backward',var);
  opt=contset(opt,'FunTolerance', 1e-6);
  opt=contset(opt,'VarTolerance', 1e-6);
  [xlcc,vlcc,slcc,hlcc,flcc]=cont(@limitcycle,x0,v0,opt);
  %
  save(analysis_output,'xlcc','ncol','ntst')
  ")
  #
  return ms
  #

end

function MORFE_compute_backbone_modal(analysis_name,analysis_number=1)
  #
  comb = readdlm(analysis_name*"/comb_a_0/monomials.txt")
  maps = readdlm(analysis_name*"/manifold_a_0/psi.txt")
  #
  po_dir = "./"*analysis_name*"/backbone_"*string(analysis_number)*".mat"
  #
  nc = size(comb)[1]
  ndofs = size(comb)[2]
  neig = size(maps)[2]
  #
  data = matread(po_dir)
  #
  xlcc = data["xlcc"]
  ncol = data["ncol"]
  ntst = data["ntst"]
  #
  npoints = Int(ncol*ntst)
  npo = size(xlcc)[2]
  vars = zeros(Float64,ndofs)
  field = zeros(Float64,neig)
  frf = zeros(npo,neig+1)
  #
  for po = 1:npo
  #
    ω = 2*π/xlcc[end-1,po]
    frf[po,1] = ω
    for p = 1:npoints
      # extract variables
      fill!(field,0.0)
      #
      for d = 1:ndofs
        vars[d] = xlcc[d+(p-1)*ndofs,po]
      end
      #
      for c = 1:nc
        #
        Πi = 1.0
        for d = 1:ndofs
          Πi *= vars[d]^comb[c,d]
        end
        #
        for d = 1:neig
          field[d] += maps[c,d]*Πi
        end
        #
      end
      #
      for d = 1:neig
        if (abs(field[d])>frf[po,1+d])
          frf[po,1+d] = abs(field[d])
        end
      end
      #
    end
    #
  end
  #
  return frf
  #
end

# reconstruction with only autonomous vector till the paper on the forcing is
# becomes public
function MORFE_compute_frc_modal(analysis_name,Ω_list,analysis_number=1)
  #
  comb = readdlm(analysis_name*"/comb_a_0/monomials.txt")
  maps = readdlm(analysis_name*"/manifold_a_0/psi.txt")
  #
  po_dir = "./"*analysis_name*"/frc_"*string(analysis_number)*".mat"
  #
  nc = size(comb)[1]
  ndofs = size(comb)[2]
  nΩ = 2*size(Ω_list)[1]
  neig = size(maps)[2]
  #
  data = matread(po_dir)
  #
  xlcc = data["xlcc"]
  ncol = data["ncol"]
  ntst = data["ntst"]
  #
  npoints = Int(ncol*ntst)
  npo = size(xlcc)[2]
  vars = zeros(Float64,ndofs)
  field = zeros(Float64,neig)
  frf = zeros(npo,neig+1)
  #
  for po = 1:npo
  #
    ω = xlcc[end,po]
    frf[po,1] = ω
    for p = 1:npoints
      # extract variables
      fill!(field,0.0)
      #
      for d = 1:ndofs
        vars[d] = xlcc[d+(p-1)*(ndofs+nΩ),po]
      end
      #
      for c = 1:nc
        #
        Πi = 1.0
        for d = 1:ndofs
          Πi *= vars[d]^comb[c,d]
        end
        #
        for d = 1:neig
          field[d] += maps[c,d]*Πi
        end
        #
      end
      #
      for d = 1:neig
        if (abs(field[d])>frf[po,1+d])
          frf[po,1+d] = abs(field[d])
        end
      end
      #
    end
    #
  end
  #
  return frf
  #
end


text:
module test_triangulate_catch_exit
using Triangulate

function test()
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Float64}([1.0 0.0 ; 0.0 1.0 ;])
    try
        (triout, vorout)=triangulate("Q", triin)
    catch err
        if typeof(err)==TriangulateError
            println("Catched TriangulateError")
            return true
        end
    end
    return false
end
end


text:
## Run operations
function numruns(re::AbstractVector)
  n = 1
  current = re[1]
  for i in 2:length(re)
    if re[i] != current
      n += 1
      current = re[i]
    end
  end
  return(n)
end

# Run End Encode: Like RLE, but return (runvalues,runends) rather than (runvalues,runlengths)
function ree{T}(x::AbstractVector{T})
  xlen = length(x)
  xlen < 2 && return( (x,[xlen]) )
  nrun = numruns(x)
  runvalues = similar(x,nrun)
  runends = Vector{Int}(nrun)
  run = 1
  current = x[1]
  for i in 2:xlen
    if x[i] != current
      runvalues[run] = current
      runends[run] = i-1
      current = x[i]
      run = run + 1
    end
  end
  runvalues[nrun] = current
  runends[nrun] = xlen
  return( (runvalues,runends) )
end

# Recompress runvalues and runends for an RLEVector
function numruns(runvalues, runends)
  len = length(runends)
  length(runends) != len && throw(ArgumentError("runvalues and runends must be the same length."))
  len < 2 && return(len)
  n = 1
  current_val = runvalues[1]
  current_end = runends[1]
    @inbounds for i in 2:len
        rv = runvalues[i]
        re = runends[i]
    if rv != current_val && re != current_end
      n = n + 1
      re < current_end && throw(ArgumentError("The provided runends were not sorted, please use cumsum(runlengths) to get the right values."))
      current_val = rv
      current_end = re
    end
  end
  return(n)
end

# Tidy up an existing (mostly) Run End Encoded vector pair, dropping zero length runs and fixing any runvalue runs
function ree(runvalues, runends)
  ree(runvalues, runends, numruns(runvalues, runends))
end

function ree(runvalues, runends, nrun)
  newv = similar(runvalues,nrun)
  newe = similar(runends,nrun)
  current_val = runvalues[1]
  current_end = runends[1]
  n = 1
    @inbounds for i in 2:length(runvalues)
        rv = runvalues[i]
        re = runends[i]
    if runends[i] != current_end
      if rv != current_val
        newv[n] = current_val
        newe[n] = current_end
        n = n + 1
        current_val = rv
      end
      current_end = re
    end
  end
  @inbounds newv[n] = current_val
  @inbounds newe[n] = current_end
  return( (newv,newe))
end

function ree(x)
  return( ([x],[1]) )
end

function inverse_ree(runvalues,runends)
  len = length(runvalues)
  len != length(runends) && throw(ArgumentError("runvalues and runends must be of the same length."))
  len == 0 && return(similar(runvalues,0))
  n = runends[end]
  rval = Array(eltype(runvalues),n)
  j=1
  @inbounds for i in 1:n
    rval[i] = runvalues[j]
    if runends[j] == i
      j = j + 1
    end
  end
  return(rval)
end


text:
@testset "Dynamic Oracles" begin

    error_cb(args...) = nothing

    oracle = DynamicOracle(ArcEager(), arc=typed)

    model(cfg) = nothing

    graph = test_sentence("economicnews.conll")

    model(cfg) = static_oracle(cfg, graph, typed)
    function error_cb(x, ŷ, y)
        @assert false
    end

    cfg = initconfig(oracle, graph)
    nocost(t, cfg) = DependencyTrees.cost(t, cfg, graph) == 0
    while !isfinal(cfg)
        pred = model(cfg)
        G = gold_transitions(oracle, cfg, graph)
        @test pred in G
        @test all(t -> nocost(t, cfg), G)
        cfg = pred(cfg)
    end

    # make sure this works the same for untyped oracles too
    oracle_ut = DynamicOracle(ArcEager(), arc=untyped)
    DependencyTrees.xys(oracle_ut, graph)
    model(cfg) = static_oracle(cfg, graph, untyped)
    cfg = initconfig(oracle_ut.system, graph)
    while !isfinal(cfg)
        pred = model(cfg)
        G = gold_transitions(oracle_ut, cfg, graph)
        @test pred in G
        @test any(t -> !nocost(t, cfg), [Shift(), Reduce(), LeftArc(), RightArc()])
        @test all(t -> nocost(t, cfg), G)
        cfg = pred(cfg)
    end

    tbfile = joinpath(@__DIR__, "data", "wsj_0001.dp")
    treebank = Treebank{TypedDependency}(tbfile, add_id=true)

    for (cfg, gold) in xys(oracle, graph)
        @test length(gold) >= 1
    end

    @testset "Projectivity" begin
        oracle = DynamicOracle(ArcHybrid())
        tb = Treebank{CoNLLU}(joinpath(@__DIR__, "data", "nonprojective1.conll"))
        @test length(collect(tb)) == 1
        @test length(collect(DependencyTrees.xys(oracle, tb))) == 0

        for tree in tb
            if !isprojective(tree)
                @test isempty(oracle(tree))
            else
                @test !isempty(oracle(tree))
            end
        end
    end

    @testset "Dynamic Iteration" begin
        oracle = DynamicOracle(ArcHybrid(), arc=untyped)
        policy = NeverExplore()

        tb = Treebank{CoNLLU}(joinpath(@__DIR__, "data", "hybridtests.conll"))

        for tree in treebank
            for (cfg, G) in xys(oracle, tree)
                @test cfg isa DT.ArcHybridConfig
                for t in G
                    T = typeof(t)
                    @test T <: LeftArc || T <: RightArc || T <: Shift
                end
            end

            for state in oracle(tree)
                @test state.G ⊆ state.A
                t = policy(state)
                @test t in state.A
            end
        end
    end
end

@testset "Transition Output Spaces" begin
    @test transition_space(ArcEager()) == [LeftArc(), RightArc(), Reduce(), Shift()]
    @test transition_space(ArcEager(), ["a","b"]) == [LeftArc("a"),  LeftArc("b"),
                                                      RightArc("a"), RightArc("b"),
                                                      Reduce(), Shift()]

    @test transition_space(ArcHybrid()) == [LeftArc(), RightArc(), Shift()]
    @test transition_space(ArcHybrid(), ["a","b"]) == [LeftArc("a"),  LeftArc("b"),
                                                       RightArc("a"), RightArc("b"),
                                                       Shift()]

    @test transition_space(ArcStandard()) == [LeftArc(), RightArc(), Shift()]
    @test transition_space(ArcStandard(), ["a","b"]) == [LeftArc("a"),  LeftArc("b"),
                                                         RightArc("a"), RightArc("b"),
                                                         Shift()]

    @test transition_space(ArcSwift(); max_k=2) == [LeftArc(1),  LeftArc(2),
                                                    RightArc(1), RightArc(2),
                                                    Shift()]
    ts1 = Set(transition_space(ArcSwift(), ["a","b"]; max_k=2))
    ts2 = Set([LeftArc(1, "a"),  LeftArc(2, "a"),
               LeftArc(1, "b"),  LeftArc(2, "b"),
               RightArc(1, "a"), RightArc(2, "a"),
               RightArc(1, "b"), RightArc(2, "b"),
               Shift()])
    @test  ts1 == ts2

    @test transition_space(ListBasedNonProjective()) == [LeftArc(), RightArc(), NoArc(), Shift()]
    @test transition_space(ListBasedNonProjective(), ["a","b"]) == [LeftArc("a"), LeftArc("b"),
                                                                    RightArc("a"), RightArc("b"),
                                                                    NoArc(), Shift()]
end

@testset "Exploration Policies" begin
    always1, never1 = AlwaysExplore(), NeverExplore()
    always2, never2 = ExplorationPolicy(1), ExplorationPolicy(0)
    @test showstr(always1) == "AlwaysExplore"
    @test showstr(never1) == "NeverExplore"
    for i = 1:10
        @test always1() == always2() == true
        @test never1()  == never2()  == false
    end
end

@testset "Feature Extraction" begin

    function check(oracle, gold_tree)
        for state in oracle(gold_tree)
            cfg = state.cfg
            stk, buf = stack(cfg), buffer(cfg)
            ts = [buffertoken(cfg, b) for b in buf]
            @test id.(buffertoken(cfg, i) for i in 1:length(buf)) == buf
            @test id.(stacktoken(cfg, i) for i in length(stk):-1:1) == stk
            @test id(stacktoken(cfg, 10000)) == -1
        end
        return true
    end

    treebank = Treebank{CoNLLU}(joinpath(@__DIR__, "data", "hybridtests.conll"))
    for system in (ArcEager(), ArcHybrid())
        for O in (StaticOracle, DynamicOracle)
            oracle = O(system, arc=untyped)
            @test all(check(oracle, tree) for tree in treebank)
        end
    end
end


text:
module DickeModel
  export ClassicalDicke
  export ClassicalSystems
  export DickeBCE
  export UPOs
  export TWA
  export FOTOCTWA
  export EnergyShellProjections
  export ClassicalLMG
  export PhaseSpaces
      include("PhaseSpaces.jl")
      include("ClassicalSystems.jl")
      include("TWA.jl")
      include("ClassicalDicke.jl")
      include("DickeBCE.jl")
      include("EnergyShellProjections.jl")
      include("ClassicalLMG.jl")
      include("UPOs.jl")
end


text:
import Pkg

Pkg.update()

root_directory = dirname(@__DIR__)

bpfnative = Pkg.PackageSpec(path = root_directory)
Pkg.develop(bpfnative)
Pkg.build()
Pkg.precompile()
if Sys.islinux()
Pkg.test("BPFnative")
else
import BPFnative
end


text:
abstract type FluidFluidOperator{N} <: AbstractSparseArray{Float64,Int,N} end

# show(io::IO, mime::MIME"text/plain", B::FluidFluidOperator) = show(io, mime, B.array)

SparseArrays.nonzeros(B::FluidFluidOperator) = nonzeros(B.array)
SparseArrays.nonzeroinds(R::FluidFluidOperator) = SparseArrays.nonzeroinds(R.array)

struct FluidFluidGenerator <: FluidFluidOperator{2}
    array::SparseMatrixCSC{Float64,Int}
    ffq::FluidFluidQueue
end
FluidFluidGenerator(ffq::FluidFluidQueue) = FluidFluidGenerator(build_full_generator(ffq.dq).B,ffq)

struct RatesOperator <: FluidFluidOperator{1}
    array::SparseVector{Float64,Int}
    ffq::FluidFluidQueue
end

function RatesOperator(ffq::FluidFluidQueue)
    r = abs.([ffq.rates.lwr[:];
        repeat(ffq.rates.interior[:],inner=n_bases_per_cell(ffq.dq));
        ffq.rates.upr[:]])
    return RatesOperator(SparseVector{Float64,Int}([r[i]!=0.0 ? 1.0./r[i] : 0.0 for i in 1:length(r)]), ffq)
end

struct InOutGenerator <: FluidFluidOperator{2}
    array::SparseMatrixCSC{Float64,Int}
    ffq::FluidFluidQueue
end

function InOutGenerator(B::FluidFluidGenerator,s::Float64)
    R = RatesOperator(B.ffq)
    num_nz = size(B[NotZero,:],1)
    num_z = size(B,1) - num_nz
    Ds = (SparseMatrixCSC{Float64,Int}(diagm(R.array))*B.array) - s*I(num_nz+num_z) +
        (SparseMatrixCSC{Float64,Int}(diagm(R.array))*B[:,Zero])*
        SparseMatrixCSC{Float64,Int}(inv(s*I(num_z) - Matrix(B[Zero,Zero])))*B[Zero,:]
    return InOutGenerator(Ds,B.ffq)
end

function InOutGenerator(ffq::FluidFluidQueue,s::Float64)
    B = FluidFluidGenerator(ffq)
    return InOutGenerator(B,s)
end

function build_psi(
    D::InOutGenerator;
    MaxIters::Int = 1000,
    err::Float64 = 1e-9,
    v::Bool = false,
)

    exitflag = ""

    Psi = zeros(Float64, size(D[Plus,Minus]))
    A = Matrix(D[Plus,Plus])
    B = Matrix(D[Minus,Minus])
    C = Matrix(D[Plus,Minus])
    # RA, QA = LinearAlgebra.schur(Matrix(A)) # uncomment for algorithm 1
    # RB, QB = LinearAlgebra.schur(Matrix(B)) # uncomment for algorithm 1
    OldPsi = Psi
    flag = 1
    for n = 1:MaxIters
        ## line below goes with algorithm 4, is quadratically convergent
       @time Psi = sylvester(A, B, C)
        ## line below goes with algorithm 1
        # Psi = mysyl(RA,QA,RB,QB,C)
        @show maximum(abs.(OldPsi - Psi))
        @show n
        @time if maximum(abs.(OldPsi - Psi)) < err
            flag = 0
            exitflag = string(
                "Reached err tolerance in ",
                n,
                " iterations with error ",
                string(maximum(abs.(OldPsi - Psi))),
            )
            break
        elseif any(isnan.(Psi))
            flag = 0
            exitflag = string("Produced NaNs at iteration ", n)
            break
        end
        @time OldPsi = copy(Psi)
        ## Algorithm 4 of Bean, O'Reilly, Taylor, 2008,
        ## Algorithms for the Laplace–Stieltjes Transforms of First Return Times for Stochastic Fluid Flows,
        ## Methodol Comput Appl Probab (2008) 10:381–408, DOI 10.1007/s11009-008-9077-3
        @time A = Psi * D[Minus,Plus]
        @time C = D[Plus,Minus] - A * Psi
        @time A += D[Plus,Plus]
        @time B = D[Minus,Minus] + D[Minus,Plus] * Psi

        ## Algorithm 1 of the above citation
        # A, B dont change, need only comput schur decomp once
        # C = D["+-"] + Psi * D["-+"] * Psi
    end
    if flag == 1
        exitflag = string(
            "Reached Max Iters ",
            MaxIters,
            " with error ",
            string(maximum(abs.(OldPsi - Psi))),
        )
    end
    v && println("UPDATE: Iterations for Ψ exited with flag: ", exitflag)
    return Psi
end

function build_psi(ffq::FluidFluidQueue, s::Union{Real,Complex}; kwargs...)
    return build_psi(InOutGenerator(ffq,s; kwargs...))
end

function build_xi(B::FluidFluidGenerator,Ψ::Array{Float64,2})
    # BBullet = [B["--"] B["-0"]; B["0-"] B["00"]]
    # invB = inv(Matrix(Bbullet))

    # BBulletPlus = [B["-+"]; B["0+"]]

    # solve the linear system -[ξ 0] [Bmm Bm0; B₀₋ B00]^-1 [B₋₊; B₀₊]Ψ = ξ
    # writing out the system it turns out we only need the index -- and -0
    # blocks of the inverse. Wikipedia tells us that these are

    # the next two lines are a fancy way to do an inverse of a block matrix
    tempMat = inv(Matrix(B[Zero,Zero]))
    invBmm = inv(B[Minus,Minus] - B[Minus,Zero]*tempMat*B[Zero,Minus])
    invBm0 = -invBmm*B[Minus,Zero]*tempMat

    A = -(invBmm*B[Minus,Plus]*Ψ + invBm0*B[Zero,Plus]*Ψ + I)
    b = zeros(1,size(invBmm,1))

    A[:,1] .= 1.0 # normalisation conditions
    b[1] = 1.0 # normalisation conditions

    ξ = b/A
    coeffs = zeros(size(B,1))
    coeffs[index(B.ffq,Minus)] = ξ
    return SFMDistribution(coeffs,B.ffq.dq)
end

function build_limit_dist_operators(B::FluidFluidGenerator,D::InOutGenerator,
    R::RatesOperator,Ψ::Array{<:Real,2},ξ::SFMDistribution)

    B00inv = inv(Matrix(B[Zero,Zero]))
    invBmm = inv(Matrix(D[Minus,Minus]))# inv(B[Minus,Minus] - B[Minus,Zero]*B00inv*B[Zero,Minus])
    invBm0 = -invBmm*B[Minus,Zero]*B00inv

    # B00inv = inv(Matrix(B["00"]))
    # invBmm = inv(B["--"] - B["-0"]*B00inv*B["0-"])
    # invBm0 = -invBmm*B["-0"]*B00inv

    αp = Array(transpose(ξ.coeffs[index(B.ffq,Minus)])) * -[invBmm invBm0]

    K = D[Plus,Plus] + Ψ * D[Minus,Plus]

    n₊ = size(K,1)
    n₀ = size(invBm0,2)
    n₋ = size(invBm0,1)

    # BBulletPlus = [B[("-","+"),(:,:)]; B[("0","+"),(:,:)]]

    αintegralPibullet = ((αp * [B[Minus,Plus]; B[Zero,Plus]]) / -K) * [diagm(R[Plus]) Ψ*diagm(R[Minus])]
    αintegralPi0 = -αintegralPibullet * [B[Plus,Zero];B[Minus,Zero]] * B00inv

    α = sum(αintegralPibullet) + sum(αintegralPi0) + sum(αp)

    p = αp ./ α
    integralPibullet = αintegralPibullet ./ α
    integralPi0 = αintegralPi0 ./ α

    marginalX = zeros(Float64, n₊ + n₋ + n₀)
    marginalX[index(B.ffq,Plus)] = integralPibullet[1:n₊]
    marginalX[index(B.ffq,Minus)] = integralPibullet[(n₊+1):end] + p[1:n₋]
    marginalX[index(B.ffq,Zero)] = integralPi0[:] + p[(n₋+1):end]

    Y_point_mass = zeros(Float64, n₊ + n₋ + n₀)
    Y_point_mass[index(B.ffq,Minus)] = p[1:n₋]
    Y_point_mass[index(B.ffq,Zero)] = p[n₋+1:end]
    return SFMDistribution(marginalX,B.ffq.dq), SFMDistribution(Y_point_mass,B.ffq.dq), K
end


text:
using Distributed
addprocs(2)

using Test
using MLJTuning
using MLJBase

# Display Number of processes and if necessary number
# of Threads
@info "nworkers: $(nworkers())"
@static if VERSION >= v"1.3.0-DEV.573"
    @info "nthreads: $(Threads.nthreads())"
else
    @info "Running julia $(VERSION). Multithreading tests excluded. "
end

include("test_utilities.jl")

print("Loading some models for testing...")
# load `Models` module containing models implementations for testing:
include("models.jl") #precompile first in master node
print("\r                                           \r")
@everywhere include("models.jl")

@testset "utilities" begin
    @test include("utilities.jl")
end

@testset "selection heuristics" begin
    include("selection_heuristics.jl")
end

@testset "tuned_models.jl" begin
    @test include("tuned_models.jl")
end

@testset "range_methods" begin
    @test include("range_methods.jl")
end

@testset "grid" begin
    @test include("strategies/grid.jl")
end

@testset "random search" begin
    @test include("strategies/random_search.jl")
end

@testset "learning curves" begin
        @test include("learning_curves.jl")
end

# @testset "julia bug" begin
#     @test include("julia_bug.jl")
# end