text:
#=
Human friendly printing.
=#

"""No more red backgrounds in IJulia output. Also, `@info` will be less verbose."""
function prettify_logging_in_IJulia()
    global_logger(get_interactive_logger())
    redirect_stderr(stdout)
end

# The default logger in IJulia is a `SimpleLogger(stderr)`, which prints source information
# not only for `@warn` but also for `@info`. The default logger in the Julia REPL is a
# `ConsoleLogger(stdout)`, which is smarter.
get_interactive_logger() = ConsoleLogger(stdout)
    # `stdout` is a writeable global (e.g. IJulia modifies it). Hence this is a getter and
    # not a constant.

"""Like `dump` but without types, with colors, and with `:compact` printing by default."""
function dumps(io::IO, @nospecialize(x); depth = 0)
    if fieldcount(typeof(x)) == 0
        printstyled(io, x, "\n", color = :light_blue)
    else
        printstyled(io, nameof(typeof(x)), "\n", color = :light_black)
        depth += 1
        indent = "  "^depth
        for fieldname in fieldnames(typeof(x))
            print(io, indent, fieldname, ": ")
            dumps(io, getproperty(x, fieldname); depth)
                # No infinite recursion guard. CBA.
        end
    end
end
dumps(x) = dumps(IOContext(stdout, :compact => true, :limit => true), x)

"""
    showsome(x; kw...)
    showsome(io::IO, x; kw...)

Print an overview of the vector `x`, by showing example entries. The first and last
entries are shown, and a few randomly sampled entries in between.

Keyword arguments `nfirst`, `nlast`, and `nsample` determine how many samples to print at
the beginning, end, and in between. Each defaults to `2`. `io` is `stdout` if not specified.
"""
showsome(x        ; nfirst = 2, nlast = 2, nsample = 2) = showsome(stdout, x; nfirst, nlast, nsample)
showsome(io::IO, x; nfirst = 2, nlast = 2, nsample = 2) = begin
    if length(x) ≤ nfirst + nlast + nsample
        nfirst = length(x)
        nlast, nsample = 0, 0
    end
    println(io, summary(x), ":")  # eg "640-element Vector{String}:"
    all_i = 1:length(x)  # Always `Int`s
    all_ix = eachindex(x)  # Generally same, but can also be `CartesianIndex`, `Symbol`, …
    shown_i = @views vcat(
        all_i[1:nfirst],
        all_i[nfirst+1:end-nlast] |> (is -> sample(is, nsample)) |> sort,
        all_i[end-nlast+1:end],
    )
    shown_ix = all_ix[shown_i]
    padlen = 1 + maximum(length, string.(shown_ix))
        # The initial extra `1` is the one-space indent of vanilla Vector printing.
    printrow(ix) = println(io, lpad(ix, padlen), ": ", repr(x[ix]))
    printdots() = println(io, lpad("⋮", padlen))
    first(shown_i) == first(all_i) || printdots()
    for (i, inext) in ziplongest(shown_i, shown_i[2:end])
        printrow(all_ix[i])
        isnothing(inext) || inext == i + 1 || printdots()
    end
    last(shown_i) == last(all_i) || printdots()
    return nothing
end


text:
using FractalTools
using Test

@testset "FractalTools.jl" begin
    # Write your own tests here.
end

ext: jl


text:
mutable struct Node
	nodeindex::Int64
	outedges::Array{Int64,1}
	inedges::Array{Int64,1}
end

mutable struct Edge
	edgeindex::Int64
	origin::Int64
	target::Int64
	active::Bool
end

mutable struct Graph
	nodelist::Array{Node,1}
	edgelist::Array{Edge,1}
	edgecounter::Int64
	nodecounter::Int64
end

function addnode!(graph::Graph)
	i = graph.nodecounter
	v = Node(i, Int64[], Int64[])
	push!(graph.nodelist, v)
	graph.nodecounter += 1
end

function addedge!(graph::Graph, start::Int64, fin::Int64)
	i = graph.edgecounter
	e = Edge(i, start, fin, true)
	push!(graph.edgelist, e)
	push!(graph.nodelist[start].outedges, i)
	push!(graph.nodelist[fin].inedges, i)
	graph.edgecounter += 1
end

function numedges(graph::Graph)
	return length(graph.edgelist)
end

function numnodes(graph::Graph)
	return length(graph.nodelist)
end

mutable struct NetStore
	adjmat::Matrix{Float64}
	pvec::Vector{Float64}
	kappa::Float64
end

mutable struct NetStore2
	adjmat::Matrix{Float64}
	pvec::Vector{Float64}
	incidence::Matrix{Float64}
	kappa::Float64
	synctol::Float64
	alpha::Float64
	maxflow::Float64
	sync::Bool
	trip::Bool
end

max_stars_repo_name: o-smith/GridResilience

max_issues_repo_name: o-smith/GridResilience

max_forks_repo_name: o-smith/GridResilience


text:
module UAPIC

using FFTW
using LinearAlgebra

include("ua_type.jl")
include("meshfields.jl")
include("integrate.jl")
include("gnuplot.jl")
include("poisson.jl")
include("particles.jl")
include("ua_steps.jl")
include("landau.jl")
include("compute_rho_cic.jl")

end # module


text:
module Run

using IOCapture: IOCapture
using TimerOutputs
import Term: RenderableText, Panel
import MyterialColors: green_lighter

using jcontrol
import jcontrol: Track
import ..bicycle: State, Bicycle
import ..control: ControlOptions, realistict_control_options, default_control_options

const to = TimerOutput()

export run_mtm

function run_mtm(
    problemtype::Symbol,
    supports_density::Number;
    track::Union{Nothing,Track}=nothing,
    icond::Union{Nothing,State}=nothing,
    fcond::Union{Nothing,State, Symbol}=nothing,
    control_options::Union{ControlOptions,Symbol}=:default,
    showtrials::Union{Nothing,Int64}=nothing,
    n_iter::Int=1000,
    tol::Float64=1e-12,
    verbose::Int=0,
    timed::Bool=false,
    showplots::Bool=true,
    quiet::Bool=false,
    α::Float64=0.0, # cost of Fu
    γ::Float64=0.0,  # cost of δ̇
    bike::Union{Nothing, Bicycle}=nothing,
    waypoint=nothing
)
    problemtype = problemtype == :kinematics ? KinematicsProblem() : DynamicsProblem()
    δt = 0.01 # Δt for forward integration

    # ---------------------------------------------------------------------------- #
    #                                DEFINE PARAMERS                               #
    # ---------------------------------------------------------------------------- #
    # create track object
    track = isnothing(track) ? FULLTRACK : track

    # create bike
    bike = isnothing(bike) ? Bicycle() : bike

    # get control options
    if control_options == :default
        control_options = default_control_options
    elseif control_options == :realistic
        control_options = realistict_control_options
    end
    @assert control_options isa ControlOptions "Control options is not a ControlOptions type: $(typeof(control_options)) $control_options"

    # define initial and final conditions
    icond = isnothing(icond) ? State(; x=track.X[1], y=track.Y[1], u=25, ω=2, ψ=.1) : icond
    # fcond = isnothing(fcond) ? State(; u=40, n=0, ψ=0) : fcond

    # ---------------------------------------------------------------------------- #
    #                                   FIT MODEL                                  #
    # ---------------------------------------------------------------------------- #
    # supports_density = 1 -> 100 supports for the whole track, adjust by track length
    n_supports = (Int ∘ round)(supports_density * 100 * (track.S[end] - track.S[1]) / 261)
    # @info "Running with" n_supports
    control_model = @timeit to "solve control" create_and_solve_control(
        problemtype,
        n_supports,
        track,
        bike,
        control_options,
        icond,
        fcond;
        n_iter=n_iter,
        tollerance=tol,
        verbose=verbose,
        quiet=quiet,
        α=α,
        γ=γ,
        waypoint=waypoint,
    )

    # ---------------------------------------------------------------------------- #
    #                              FORWARD INTEGRATION                             #
    # ---------------------------------------------------------------------------- #
    solution = @timeit to "run forward model" run_forward_model(
        problemtype, track, control_model; δt=δt
    )

    # ---------------------------------------------------------------------------- #
    #                                 VISUALIZATION                                #
    # ---------------------------------------------------------------------------- #
    if showplots
        # visualize results
        @timeit to "plot control" summary_plot(problemtype, control_model, :time)

        # visualize
        trials = if !isnothing(showtrials)
            jcontrol.load_trials(; keep_n=showtrials, method=:efficient)
        else
            nothing
        end
        @timeit to "plot ODEs" summary_plot(
            solution, control_model, track, bike; trials=trials
        )
    end

    # -------------------------------- timing info ------------------------------- #
    # show timing/allocation data
    if timed && !quiet
        c = IOCapture.capture() do
            show(to)
        end
        print(
            "\n\n" / Panel(
                RenderableText(c.output, green_lighter);
                style="green dim",
                title="Timing/Allocations info",
                title_style="green underline bold",
                justify=:center,
            ),
        )
    end

    return track, bike, control_model, solution
end
end

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
JGE.initialize(glfw = true, debug = false)
JGE.test(debug = true)


text:
"""This julia script creates a uniform map with a single initialization spot
in the middle. It takes arguments from the commandline and creates a map in
the directory it is run from.
Arguments in order are map size, simulation length, temperature and above ground carbon
--fn or --filename overrides the autogenerated map name which is plainmap_ whereto
the numeric parameters follow in order

Robin Rölz, 22/2/2021"""

using ArgParse

s = ArgParseSettings()
@add_arg_table! s begin
    "--fn", "--filename"
        help = "sets name for new mapfile" #so far only puts it in wd
        arg_type = String
    "s"
        help = "map size in x and y"
        arg_type = Int
        required = true
    "t"
        help = "Simulation length"
        arg_type = Int
        default = 10
    "temp"
        help = "temperature in Kelvin"
        default = 293
    "agc"
        help = "above ground carbon"
        default = 80 #Not sure about the default value I should choose
end

"""
    getfilename (size, run_length, temperature, agc)

Short utility function to generate a filename based on the input
parameters.
"""
function getfilename(size, run_length, temperature, agc)
    newfilename = string("plainmap_",
                            size, "_",
                            run_length, "_",
                            temperature, "_",
                            agc, ".map")
    return(newfilename)
end
    """
    createplainmap(size, run_length, temperature, agc)

creation function. prints a homogenous map of the specified size
with a single initialization patch
"""
function createplainmap(size::Int, run_length::Int, temperature, agc)
    map_text= string("## ZOSTEROPS EXPERIMENT MAP\n\n",
                    "# Timesteps \n$run_length\n\n",
                    "# Simulation arena - created by 'createplainmap.jl'\n\n",
                    "# <id> <x> <y> <temperature> <agc> [parameters]")
    i = 1:size^2
    x = repeat(1:size, size)
    y = repeat(1:size, inner=size)
    middle = ceil(size/2)
    for eachi in i
        line = string(
            "\n",
            eachi, " ",
            x[eachi], " ",
            y[eachi], " ",
            "temp=", temperature, " ",
            "prec=", agc, " ")
        #TODO: fix map parsing. So far first 'initpop' is ignored. This works but shouldn't
        if ((x[eachi]==middle) && (y[eachi]== middle)) || (i[eachi]==1)
            line=line*"initpop"
        end
        map_text=map_text*line
    end
    return(map_text)
end

println("Creating map")
arg = parse_args(s)
filename = arg["fn"]
if isnothing(filename)
    filename = getfilename(arg["s"], arg["t"], arg["temp"], arg["agc"])
end
path = joinpath(pwd(), filename)
open(path, "w") do map
    write(map, createplainmap( arg["s"],
                arg["t"],
                arg["temp"],
                arg["agc"]
                )
    )
end
println("Finished creating map")

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
include("basisfunctions.jl")
include("fit.jl")
include("designmatrix.jl")
include("splines.jl")
include("predict.jl")
#include("clusterpermutation.jl")


text:
module TestLoggedStatic

using Test
import ThreadPools: LoggedStaticPool
using ThreadPools

include("util.jl")

@testset "LoggedStaticPool" begin

    @testset "tforeach" begin
        @testset "foreground" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            primary = Threads.nthreads() == 1
            fn! = (x) -> begin
                Threads.threadid() == 1 && (primary = true)
                x.data += 1
            end
            pool = LoggedStaticPool()
            tforeach(fn!, pool, objs)
            close(pool)
            @test [x.data for x in objs] == collect(2:N+1)
            @test primary
            @test length(pool.recs) == N*2
            @inferred tforeach(fn!, pool, objs)
        end

        @testset "background" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            fn! = (x) -> begin
                Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                x.data += 1
            end
            pool = LoggedStaticPool(2)
            tforeach(fn!, pool, objs)
            close(pool)
            @test [x.data for x in objs] == collect(2:N+1)
            @test length(pool.recs) == N*2
            @inferred tforeach(fn!, pool, objs)
        end
    end

    @testset "tmap" begin
        @testset "foreground" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            primary = Threads.nthreads() == 1
            fn! = (x) -> begin
                Threads.threadid() == 1 && (primary = true)
                x.data
            end
            pool = LoggedStaticPool()
            @test tmap(fn!, pool, objs) == collect(1:N)
            close(pool)
            @test primary
            @test length(pool.recs) == N*2
            @inferred tmap(fn!, pool, objs)
        end

        @testset "background" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            fn! = (x) -> begin
                Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                x.data
            end
            pool = LoggedStaticPool(2)
            @test tmap(fn!, pool, objs) == collect(1:N)
            close(pool)
            @test length(pool.recs) == N*2
            @inferred tmap(fn!, pool, objs)
        end
    end

    @testset "twith" begin
        @testset "foreground" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            primary = Threads.nthreads() == 1
            pool = twith(LoggedStaticPool()) do pool
                tforeach(pool, objs) do x
                    Threads.threadid() == 1 && (primary = true)
                    x.data += 1
                end
            end
            @test [x.data for x in objs] == collect(2:N+1)
            @test primary
            @test length(pool.recs) == N*2
        end

        @testset "background" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            pool = twith(LoggedStaticPool(2)) do pool
                tforeach(pool, objs) do x
                    Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                    x.data += 1
                end
            end
            @test [x.data for x in objs] == collect(2:N+1)
            @test length(pool.recs) == N*2
        end
    end

    @testset "@tthreads" begin
        @testset "foreground" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            primary = Threads.nthreads() == 1
            pool = LoggedStaticPool()
            @tthreads pool for obj in objs
                Threads.threadid() == 1 && (primary = true)
                obj.data += 1
            end
            close(pool)
            @test primary
            @test [x.data for x in objs] == collect(2:N+1)
            @test length(pool.recs) == N*2
        end

        @testset "background" begin
            N = 2 * Threads.nthreads()
            objs = [TestObj(x) for x in 1:N]
            pool = LoggedStaticPool(2)
            @tthreads pool for obj in objs
                Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
                obj.data += 1
            end
            close(pool)
            @test [x.data for x in objs] == collect(2:N+1)
            @test length(pool.recs) == N*2
        end
    end

end

end # module


text:
# Julia wrapper for header: /usr/include/scip/cons_linear.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrLinear(scip)
    ccall((:SCIPincludeConshdlrLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPincludeLinconsUpgrade(scip, linconsupgd, priority, conshdlrname)
    ccall((:SCIPincludeLinconsUpgrade, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cvoid}, Cint, Cstring), scip, linconsupgd, priority, conshdlrname)
end

function SCIPcreateConsLinear(scip, cons, name, nvars, vars, vals, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nvars, vars, vals, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicLinear(scip, cons, name, nvars, vars, vals, lhs, rhs)
    ccall((:SCIPcreateConsBasicLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cdouble, Cdouble), scip, cons, name, nvars, vars, vals, lhs, rhs)
end

function SCIPcopyConsLinear(scip, cons, sourcescip, name, nvars, sourcevars, sourcecoefs, lhs, rhs, varmap, consmap, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode, _global, valid)
    ccall((:SCIPcopyConsLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Ptr{SCIP_}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cdouble, Cdouble, Ptr{SCIP_HASHMAP}, Ptr{SCIP_HASHMAP}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, Ptr{UInt32}), scip, cons, sourcescip, name, nvars, sourcevars, sourcecoefs, lhs, rhs, varmap, consmap, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode, _global, valid)
end

function SCIPaddCoefLinear(scip, cons, var, val)
    ccall((:SCIPaddCoefLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, val)
end

function SCIPchgCoefLinear(scip, cons, var, val)
    ccall((:SCIPchgCoefLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, val)
end

function SCIPdelCoefLinear(scip, cons, var)
    ccall((:SCIPdelCoefLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}), scip, cons, var)
end

function SCIPgetLhsLinear(scip, cons)
    ccall((:SCIPgetLhsLinear, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetRhsLinear(scip, cons)
    ccall((:SCIPgetRhsLinear, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPchgLhsLinear(scip, cons, lhs)
    ccall((:SCIPchgLhsLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cdouble), scip, cons, lhs)
end

function SCIPchgRhsLinear(scip, cons, rhs)
    ccall((:SCIPchgRhsLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cdouble), scip, cons, rhs)
end

function SCIPgetNVarsLinear(scip, cons)
    ccall((:SCIPgetNVarsLinear, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarsLinear(scip, cons)
    ccall((:SCIPgetVarsLinear, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetValsLinear(scip, cons)
    ccall((:SCIPgetValsLinear, libscip), Ptr{Cdouble}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetActivityLinear(scip, cons, sol)
    ccall((:SCIPgetActivityLinear, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}), scip, cons, sol)
end

function SCIPgetFeasibilityLinear(scip, cons, sol)
    ccall((:SCIPgetFeasibilityLinear, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}), scip, cons, sol)
end

function SCIPgetDualsolLinear(scip, cons)
    ccall((:SCIPgetDualsolLinear, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetDualfarkasLinear(scip, cons)
    ccall((:SCIPgetDualfarkasLinear, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetRowLinear(scip, cons)
    ccall((:SCIPgetRowLinear, libscip), Ptr{SCIP_ROW}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPupgradeConsLinear(scip, cons, upgdcons)
    ccall((:SCIPupgradeConsLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_CONS}}), scip, cons, upgdcons)
end

function SCIPclassifyConstraintTypesLinear(scip, linconsstats)
    ccall((:SCIPclassifyConstraintTypesLinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_LINCONSSTATS}), scip, linconsstats)
end


text:
import Pkg
let
    pkgs = ["Revise", "Parameters", "OhMyREPL"]
    for pkg in pkgs
    if Base.find_package(pkg) === nothing
        Pkg.add(pkg)
    end
    end
end

println("go!")

#Pkg.activate(".")
#using Revise
#using OhMyREPL

max_stars_repo_name: timm/misc

max_issues_repo_name: timm/misc

max_forks_repo_name: timm/misc


text:
using Test
using DFTK: load_psp, list_psp

@testset "Check reading all HGH pseudos" begin
    for record in list_psp()
        psp = load_psp(record.identifier)
        @test psp.identifier == record.identifier
        @test psp.Zion == record.n_elec_valence
    end
end

@testset "Check load_psp routine selectively" begin
    psp_Cu = load_psp(:Cu, functional="lda", family="hgh", core=:semicore)
    @test psp_Cu.identifier == "hgh/lda/cu-q19.hgh"

    psp_Au = load_psp(:Au, functional="pbe", family="hgh", core=:fullcore)
    @test psp_Au.identifier == "hgh/pbe/au-q11.hgh"

    for element in (:Cu, :Ni, :Au)
        sc = load_psp(element, functional="lda", family="hgh", core=:semicore)
        fc = load_psp(element, functional="lda", family="hgh", core=:fullcore)
        @test sc.Zion > fc.Zion
    end
end


text:
"""
    Proof of concept for Caesar-ROS integration
    (check Caesar Docs for details)
    https://juliarobotics.org/Caesar.jl/latest/examples/using_ros/

    Prerequisites:
    - source /opt/ros/noetic/setup.bash
    - cd ~/thecatkin_ws
        - source devel/setup.bash in all 3 terminals
    - run roscore in one terminal
    - Then run this Julia in another terminal/process.

    Input:
    - Make sure the rosbag is in ~/data/Marine/philos_car_far.bag

    Output:
    - Generates output dfg tar and data folder at /tmp/caesar/philos
        containing data from the bagfile, see below for details.

    Future:
    - ROS msg replies
    - periodic export of factor graph object
"""

## Prepare python version
using Distributed
# addprocs(4)

using Pkg
Distributed.@everywhere using Pkg

Distributed.@everywhere begin
  ENV["PYTHON"] = "/usr/bin/python3"
  Pkg.build("PyCall")
end

using PyCall
Distributed.@everywhere using PyCall

## INIT
using RobotOS

# Also rosnode info
# standard types
@rosimport sensor_msgs.msg: PointCloud2
@rosimport sensor_msgs.msg: NavSatFix
@rosimport sensor_msgs.msg: CompressedImage
# @rosimport nmea_msgs.msg: Sentence
# seagrant type

# Application specific ROS message types from catkin workspace
# @rosimport seagrant_msgs.msg: radar

rostypegen()

## Load Caesar with additional tools

using Colors
using Caesar
# using Caesar._ROS

##

# using RoME
# using DistributedFactorGraphs

using DistributedFactorGraphs.DocStringExtensions
using Dates
using JSON2
using BSON
using Serialization
using FixedPointNumbers
using StaticArrays
using ImageMagick, FileIO
using Images

using ImageDraw

##

# /gps/fix              10255 msgs    : sensor_msgs/NavSatFix
# /gps/nmea_sentence    51275 msgs    : nmea_msgs/Sentence
# /radar_pointcloud_0    9104 msgs    : sensor_msgs/PointCloud2
# /velodyne_points      20518 msgs    : sensor_msgs/PointCloud2

# function handleGPS(msg, fg)
# end

"""
    $TYPEDEF
Quick placeholder for the system state - we're going to use timestamps to align all the data.
"""
Base.@kwdef mutable struct SystemState
    curtimestamp::Float64 = -1000
    cur_variable::Union{Nothing, DFGVariable} = nothing
    var_index::Int = 0
    lidar_scan_index::Int = 0
    max_lidar::Int = 3
    radar_scan_queue::Channel{sensor_msgs.msg.PointCloud2} = Channel{sensor_msgs.msg.PointCloud2}(64)
    # SystemState() = new(-1000, nothing, 0, 0, 3)
end

##

"""
    $(SIGNATURES)
Update the system state variable if the timestamp has changed (increment variable)
"""
function updateVariableIfNeeded(fg::AbstractDFG, systemstate::SystemState, newtimestamp::Float64)
    # Make a new variable if so.
    if systemstate.curtimestamp == -1000 || systemstate.cur_variable === nothing || systemstate.curtimestamp < newtimestamp
        systemstate.curtimestamp = newtimestamp
        systemstate.cur_variable = addVariable!(fg, Symbol("x$(systemstate.var_index)"), Pose2, timestamp = unix2datetime(newtimestamp))
        systemstate.var_index += 1
        systemstate.lidar_scan_index = 0
    end
    return nothing
end

"""
    $SIGNATURES

Message callback for Radar pings. Adds a variable to the factor graph and appends the scan as a bigdata element.
"""
function handleRadarPointcloud!(msg::sensor_msgs.msg.PointCloud2, fg::AbstractDFG, systemstate::SystemState)
    @info "handleRadarPointcloud!" maxlog=10

    # assume there is still space (previously cleared)
    # add new piece of radar point cloud to queue for later processing.
    put!(systemstate.radar_scan_queue, msg)

    # check if the queue still has space
    if length(systemstate.radar_scan_queue.data) < systemstate.radar_scan_queue.sz_max
        # nothing more to do
        return nothing
    end

    # Full sweep, lets empty the queue and add a variable
    # type instability
    queueScans = Vector{Any}(undef, systemstate.radar_scan_queue.sz_max)

    # get the first
    md = take!(systemstate.radar_scan_queue)
    pc2 = Caesar._PCL.PCLPointCloud2(md)
    pc_cat = Caesar._PCL.PointCloud(pc2)

    queueScans[1] = pc2

    for i in 1:length(systemstate.radar_scan_queue.data)
        # something minimal, will do util for transforming PointCloud2 next
        md = take!(systemstate.radar_scan_queue)
        # @info typeof(md) fieldnames(typeof(md))
        pc2 = Caesar._PCL.PCLPointCloud2(md)
        pc_ = Caesar._PCL.PointCloud(pc2)
        pc_cat = cat(pc_cat, pc_; reuse=true)

        queueScans[i] = (pc2)
    end

    # add a new variable to the graph
    timestamp = Float64(msg.header.stamp.secs) + Float64(msg.header.stamp.nsecs)/1.0e9
    systemstate.curtimestamp = timestamp
    systemstate.cur_variable = addVariable!(fg, Symbol("x$(systemstate.var_index)"), Pose2, timestamp = unix2datetime(timestamp), tags=[:POSE])
    systemstate.var_index += 1

    io = IOBuffer()
    serialize(io, queueScans)

    # @show datablob = pc # queueScans
    # and add a data blob of all the scans
    # Make a data entry in the graph
    addData!(   fg, :radar, systemstate.cur_variable.label, :RADAR_PC2s,
                take!(io), # get base64 binary
                # Vector{UInt8}(JSON2.write(datablob)),
                mimeType="application/octet-stream/julia.serialize",
                description="queueScans = Serialization.deserialize(PipeBuffer(readBytes))")
    #

    io = IOBuffer()
    serialize(io, pc_cat)

    addData!(   fg, :radar, systemstate.cur_variable.label, :RADAR_SWEEP,
                take!(io),
                mimeType="application/octet-stream/julia.serialize",
                description="queueScans = Serialization.deserialize(PipeBuffer(readBytes))" )
    #

    # also make and add an image of the radar sweep
    img = makeImage!(pc_cat)
    addData!(   fg, :radar, systemstate.cur_variable.label, :RADAR_IMG,
                Caesar.toFormat(format"PNG", img),
                mimeType="image/png",
                description="ImageMagick.readblob(imgBytes)" )
    #

    nothing
end

"""
    $SIGNATURES

Message callback for LIDAR point clouds. Adds a variable to the factor graph and appends the scan as a bigdata element.
Note that we're just appending all the LIDAR scans to the variables because we are keying by RADAR.
"""
function handleLidar!(msg::sensor_msgs.msg.PointCloud2, fg::AbstractDFG, systemstate::SystemState)
    @info "handleLidar" maxlog=10
    # Compare systemstate and add the LIDAR scans if we want to.
    if systemstate.cur_variable === nothing
        return nothing
    end
    timestamp = Float64(msg.header.stamp.secs) + Float64(msg.header.stamp.nsecs)/1.0e9
    @info "[$timestamp] LIDAR pointcloud sample on $(systemstate.cur_variable.label) (sample $(systemstate.lidar_scan_index+1))"

    # Check if we have enough LIDAR's for this variable
    if systemstate.lidar_scan_index >= systemstate.max_lidar
        @warn "Ditching LIDAR sample for this variable, already have enough..."
        return nothing
    end

    # Make a data entry in the graph
    ade,adb = addData!(fg, :lidar, systemstate.cur_variable.label, Symbol("LIDAR$(systemstate.lidar_scan_index)"), Vector{UInt8}(JSON2.write(msg)), mimeType="/velodyne_points;dataformat=Float32*[[X,Y,Z]]*32")

    # NOTE: If JSON, then do this to get to Vector{UInt8} - # byteData = Vector{UInt8}(JSON2.write(xyzLidarF32))

    # Increment LIDAR scan count for this timestamp
    systemstate.lidar_scan_index += 1
end

"""
    $SIGNATURES

Message callback for Radar pings. Adds a variable to the factor graph and appends the scan as a bigdata element.
"""
function handleGPS!(msg::sensor_msgs.msg.NavSatFix, fg::AbstractDFG, systemstate::SystemState)
    @info "handleGPS" maxlog=10
    if systemstate.cur_variable === nothing
        # Keyed by the radar, skip if we don't have a variable yet.
        return nothing
    end
    timestamp = Float64(msg.header.stamp.secs) + Float64(msg.header.stamp.nsecs)/10^9
    # Update the variable if needed
    # updateVariableIfNeeded(fg, systemstate, timestamp)
    @info "[$timestamp] GPS sample on $(systemstate.cur_variable.label)"

    if :GPS in listDataEntries(fg, systemstate.cur_variable.label)
        @warn "GPS sample on $(systemstate.cur_variable.label) already exist, dropping"
        return nothing
    end

    io = IOBuffer()
    JSON2.write(io, msg)
    ade,adb = addData!(fg, :gps_fix, systemstate.cur_variable.label, :GPS, take!(io),  mimeType="application/json", description="JSON2.read(IOBuffer(datablob))")

end

"""
    $SIGNATURES

Message callback for camera images.
"""
function handleCamera_Center!(msg::sensor_msgs.msg.CompressedImage, fg::AbstractDFG, systemstate::SystemState)
  @info "handleCamera_Center!" maxlog=10

  lbls = ls(fg, tags=[:POSE;])
  if length(lbls) == 0
    return nothing
  end

  lb = sortDFG(lbls)[end]

  addData!(   fg, :camera, systemstate.cur_variable.label, Symbol(:IMG_CENTER_, msg.header.seq),
              msg.data,
              mimeType="image/jpeg",
              description="ImageMagick.readblob(imgBytes); # "*msg.format*"; "*string(msg.header.stamp) )

  nothing
end

##

function main(;iters::Integer=50)
    dfg_datafolder = "/tmp/caesar/philos"
    if isdir(dfg_datafolder)
        println("Deleting old contents at: ",dfg_datafolder)
        rm(dfg_datafolder; force=true, recursive=true)
    end
    mkdir(dfg_datafolder)

    @info "Hit CTRL+C to exit and save the graph..."

    init_node("asv_feed")
    # find the bagfile
    bagfile = joinpath(ENV["HOME"],"data","Marine","philos_car_far.bag")
    bagSubscriber = RosbagSubscriber(bagfile)

    # Initialization
    fg = initfg()
    getSolverParams(fg).inflateCycles=1

    ds = FolderStore{Vector{UInt8}}(:radar, "$dfg_datafolder/data/radar")
    addBlobStore!(fg, ds)

    ds = FolderStore{Vector{UInt8}}(:gps_fix, "$dfg_datafolder/data/gps")
    addBlobStore!(fg, ds)

    # add if you want lidar also
    ds = FolderStore{Vector{UInt8}}(:lidar, "$dfg_datafolder/data/lidar")
    addBlobStore!(fg, ds)

    # add if you want lidar also
    ds = FolderStore{Vector{UInt8}}(:camera, "$dfg_datafolder/data/camera")
    addBlobStore!(fg, ds)

    # System state
    systemstate = SystemState()

    # Enable and disable as needed.
    camcen_sub = bagSubscriber("/center_camera/image_color/compressed", sensor_msgs.msg.CompressedImage, handleCamera_Center!, (fg, systemstate) )
    radarpc_sub = bagSubscriber("/broadband_radar/channel_0/pointcloud", sensor_msgs.msg.PointCloud2, handleRadarPointcloud!, (fg, systemstate) )
    # Skipping LIDAR
    # lidar_sub = Subscriber{sensor_msgs.msg.PointCloud2}("/velodyne_points", sensor_msgs.msg.PointCloud2, handleLidar!, (fg,systemstate), queue_size = 10)
    gps_sub = bagSubscriber("/gnss", sensor_msgs.msg.NavSatFix, handleGPS!, (fg, systemstate))

    @info "subscribers have been set up; entering main loop"
    # loop_rate = Rate(20.0)
    while loop!(bagSubscriber)
        iters -= 1
        iters < 0 ? break : nothing
    end

    @info "Exiting"
    # After the graph is built, for now we'll save it to drive to share.
    # Save the DFG graph with the following:
    @info "Saving DFG to $dfg_datafolder/dfg"
    saveDFG(fg, "$dfg_datafolder/dfg")

end

##

# Actually run the program and build
main(iters=100000)

## ===========================================================================================
## after the graph is saved it can be loaded and the datastores retrieved

dfg_datafolder = "/tmp/caesar/philos"

fg = loadDFG("$dfg_datafolder/dfg")

ds = FolderStore{Vector{UInt8}}(:radar, "$dfg_datafolder/data/radar")
addBlobStore!(fg, ds)

ds = FolderStore{Vector{UInt8}}(:gps_fix, "$dfg_datafolder/data/gps")
addBlobStore!(fg, ds)

# add if you want lidar also
ds = FolderStore{Vector{UInt8}}(:lidar, "$dfg_datafolder/data/lidar")
addBlobStore!(fg, ds)

# add if you want lidar also
ds = FolderStore{Vector{UInt8}}(:camera, "$dfg_datafolder/data/camera")
addBlobStore!(fg, ds)

##


text:
get_value(tuple::Tuple{DetailedOrganism, T}) where {T <: Number} = tuple[2]
get_value(organism::DetailedOrganism) = get_fitness(organism)

# mean

function mean(list)
    return sum(list) / length(list)
end

function mean(organisms::Union{Vector{DetailedOrganism}, Vector{Tuple{DetailedOrganism, T}}}) where {T <: Number}
    mean(get_value.(organisms))
end

function mean(organisms::Vector{DetailedOrganism}, conditionals)
    group_function = get_group_function(conditionals)
    mean_per_group = Dict([key => mean(group) for (key, group) in _group_by(organisms, group_function)])
    return [(organism, mean_per_group[group_function(organism)]) for organism in organisms]
end

function mean(organisms::Vector{Tuple{DetailedOrganism, T}}, conditionals) where {T <: Number}
    group_function = x -> get_group_function(conditionals)(x[1])
    mean_per_group = Dict([key => mean(group) for (key, group) in _group_by(organisms, group_function)])
    return [(organism[1], mean_per_group[group_function(organism)]) for organism in organisms]
end

# variance

function var(list)
    avg = mean(list)
    return sum([(item - avg)^2 for item in list]) / (length(list) - 1)
end

function var(organisms::Union{Vector{DetailedOrganism}, Vector{Tuple{DetailedOrganism, T}}}) where {T <: Number}
    var(get_value.(organisms))
end

function var(organisms::Vector{DetailedOrganism}, conditionals)
    group_function = get_group_function(conditionals)
    var_per_group = Dict([key => var(group) for (key, group) in _group_by(organisms, group_function)])
    return [(organism, var_per_group[group_function(organism)]) for organism in organisms]
end

function var(organisms::Vector{Tuple{DetailedOrganism, T}}, conditionals)  where {T <: Number}
    group_function = x -> get_group_function(conditionals)(x[1])
    var_per_group = Dict([key => var(group) for (key, group) in _group_by(organisms, group_function)])
    return [(organism[1], var_per_group[group_function(organism)]) for organism in organisms]
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
function initialize_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, PSY.SingleMass, A, TG, P}},
    inner_vars::AbstractVector,
) where {M <: PSY.Machine, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS} end

function initialize_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, PSY.FiveMassShaft, A, TG, P}},
    inner_vars::AbstractVector,
) where {M <: PSY.Machine, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    shaft_ix = get_local_state_ix(dynamic_device, PSY.FiveMassShaft)
    shaft_states = @view device_states[shaft_ix]
    δ0 = shaft_states[1]
    δ = δ0
    ω = shaft_states[2]
    ω_hp = ω
    ω_ip = ω
    ω_lp = ω
    ω_ex = ω
    ω_sys = ω

    #Obtain parameters
    shaft = PSY.get_shaft(dynamic_device)
    D = PSY.get_D(shaft)
    D_hp = PSY.get_D_hp(shaft)
    D_ip = PSY.get_D_ip(shaft)
    D_lp = PSY.get_D_lp(shaft)
    D_ex = PSY.get_D_ex(shaft)
    D_12 = PSY.get_D_12(shaft)
    D_23 = PSY.get_D_23(shaft)
    D_34 = PSY.get_D_34(shaft)
    D_45 = PSY.get_D_45(shaft)
    K_hp = PSY.get_K_hp(shaft)
    K_ip = PSY.get_K_ip(shaft)
    K_lp = PSY.get_K_lp(shaft)
    K_ex = PSY.get_K_ex(shaft)

    #Obtain inner vars
    τe = inner_vars[τe_var]
    τm0 = inner_vars[τm_var]

    function f!(out, x)
        τm = x[1]
        δ_hp = x[2]
        δ_ip = x[3]
        δ_lp = x[4]
        δ_ex = x[5]

        out[1] = (
            -τe - D * (ω - ω_sys) - D_34 * (ω - ω_lp) - D_45 * (ω - ω_ex) +
            K_lp * (δ_lp - δ) +
            K_ex * (δ_ex - δ)
        )

        out[2] = τm - D_hp * (ω_hp - ω_sys) - D_12 * (ω_hp - ω_ip) + K_hp * (δ_ip - δ_hp)

        out[3] = (
            -D_ip * (ω_ip - ω_sys) - D_12 * (ω_ip - ω_hp) - D_23 * (ω_ip - ω_lp) +
            K_hp * (δ_hp - δ_ip) +
            K_ip * (δ_lp - δ_ip)
        )

        out[4] = (
            -D_lp * (ω_lp - ω_sys) - D_23 * (ω_lp - ω_ip) - D_34 * (ω_lp - ω) +
            K_ip * (δ_ip - δ_lp) +
            K_lp * (δ - δ_lp)
        )
        out[5] = -D_ex * (ω_ex - ω_sys) - D_45 * (ω_ex - ω) + K_ex * (δ - δ_ex)
    end

    x0 = [τm0, δ0, δ0, δ0, δ0]
    sol = NLsolve.nlsolve(f!, x0, ftol = STRICT_NLSOLVE_F_TOLERANCE)
    if !NLsolve.converged(sol)
        @warn("Initialization in Shaft failed")
    else
        sol_x0 = sol.zero
        inner_vars[τm_var] = sol_x0[1] #τm
        shaft_states[3] = sol_x0[2] #δ_hp
        shaft_states[4] = ω #ω_hp
        shaft_states[5] = sol_x0[3] #δ_ip
        shaft_states[6] = ω #ω_ip
        shaft_states[7] = sol_x0[4] #δ_lp
        shaft_states[8] = ω #ω_lp
        shaft_states[9] = sol_x0[5] #δ_ex
        shaft_states[10] = ω #ω_ex
    end
    return
end


text:
using Reexport

@reexport module InvariantDistributions
    include("AbstractInvariantDistribution.jl")
    include("composite_types/InvariantDistribution.jl")
	include("plot_recipes.jl")
end # module

"""
    InvariantDistributions

A module containing abstract and composite types for invariant distributions. These
composite types hold the distributions obtained by the estimators in the parent
`InvariantMeasures` module.
"""
InvariantDistributions


text:
using QCMaterial
using PyCall
using LinearAlgebra
using HDF5
import Random
import PyCall: pyimport

"tau plus"
nsite = 2
n_qubit　= 2*nsite
U = 1.0
V = 1.0
μ = 1.0
ε = 1.0
d_theta = 1e-5
verbose = QCMaterial.MPI_rank == 0

#Hamiltonian
ham_op = generate_impurity_ham_with_1imp_multibath(U, V, μ, ε, nsite)

n_electron_gs = 2
@assert mod(n_electron_gs, 2) == 0
sparse_mat = get_number_preserving_sparse_operator(ham_op, n_qubit, n_electron_gs);
enes_ed = eigvals(sparse_mat.toarray());

#debug
println("Ground energy_ED=",minimum(enes_ed))

#ansatz
state0 = create_hf_state(n_qubit, n_electron_gs)
vc = uccgsd(n_qubit, orbital_rot=true, conserv_Sz_singles=true)
theta_init = rand(num_theta(vc))

#Perform VQE
cost_history, thetas_opt =
QCMaterial.solve_gs(jordan_wigner(ham_op), vc, state0, theta_init=theta_init, verbose=true,
    comm=QCMaterial.MPI_COMM_WORLD
)

#debug
println("Ground energy_VQE=",cost_history[end])

#c^{dag},c
up1 = up_index(1)
down1 = down_index(1)
up2 = up_index(2)
down2 = down_index(2)
#right_op = FermionOperator("$(up1)^ ", 1.0)
#right_op = jordan_wigner(right_op)
#left_op = FermionOperator("$(up1) ", 1.0)
#left_op = jordan_wigner(left_op)

#1down
right_op = FermionOperator("$(down1)^ ", 1.0)
right_op = jordan_wigner(right_op)
left_op = FermionOperator("$(down1) ", 1.0)
left_op = jordan_wigner(left_op)

##Ansatz -> apply_qubit_op & imag_time_evolve
vc_ex = uccgsd(n_qubit, orbital_rot=true)

#state_gs = QulacsQuantumState(n_qubit,0b0000)
state_gs = create_hf_state(n_qubit, n_electron_gs)
update_circuit_param!(vc, thetas_opt)
update_quantum_state!(vc, state_gs)

n_electron_ex = 3
state0_ex = create_hf_state(n_qubit, n_electron_ex)

taus1 = collect(range(0.0, 0.02, length=2))
println("taus1=",taus1)
#beta = taus[end]

#make taus list file
#length=8
function taus_list(file_name)
	x = Float64[]
	for i in 1:10
        tau = 0.005*(i-1) + exp((0.1*(i-1))) -1
		push!(x,tau)
	end
	fid = h5open(file_name,"w")
	fid["/test/data"] = x
	close(fid)
end

#read taus list file
function read_taus_list(file_name)
	fid = h5open(file_name,"r")
	x = fid["/test/data"][:]
	close(fid)
	return x
end

#make taus list file
taus_list("test2.h5")
#generate taus list file
taus = read_taus_list("test2.h5")
println("taus=",taus)

Gfunc_ij_list = compute_gtau(jordan_wigner(ham_op), left_op, right_op, vc_ex,  state_gs, state0_ex, taus, d_theta, verbose=verbose)
println("Gfunc_ij_list=", Gfunc_ij_list)

function write_to_txt(file_name, x, y)
    open(file_name,"w") do fp
        for i in 1:length(x)
            println(fp, x[i], " ", real(y[i]))
        end
    end
end

write_to_txt("g_of_tau_dimer_points_10.txt", taus, Gfunc_ij_list)

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
module DH32ParallelUtils

# package code goes here
export split_int, split_int_lengths

"""
Splits an integer number N into a set of n_procs contiguous, non-overlapping integer intervals [ix_start[i], ix_end[i]].

Input:

- N: number to be split.

- n_procs: number of processes.

Output: ix_start, ix_end

- ix_start: int array of start indexes.

- ix_end: int array of end indexes.
"""
function split_int(N::Int64, n_procs::Int64)
    ix_start = zeros(Int64,n_procs)
    ix_end = zeros(Int64,n_procs)
    float_bounds = range(1.0, stop=Float64(N), length=n_procs+1)
    ix_start[1] = 1
    ix_end[n_procs] = N
    for i in 2:n_procs
        ix_start[i] = round(Int64,float_bounds[i])
        ix_end[i-1] = ix_start[i]-1
    end
    ix_start, ix_end
end

"""
Split a integer N in n_procs integer pieces of approximately the same lengths.

Input: as for split_int.

Output:

- Array of n_procs integer lengths summing up to N

"""
function split_int_lengths(N::Int64, n_procs::Int64)
    ix_start, ix_end = split_int(N, n_procs)
    ix_end - ix_start .+ 1
end

end # module


text:
export IdentityLink, LogLink, LogitLink
export link, unlink, jacobian_link

abstract type Link end

"""
    link(args...)

The link function is a map that brings a parameter ``f`` in a subspace ``\\mathcal{F} \\subset \\mathbb{R}`` to ``\\mathbb{R}``.
"""
function link end

"""
    unlink(args...)

The unlink function is the inverse map of [`link`](@ref link). It brings ``\\tilde{f}`` in ``\\mathbb{R}`` to the subspace ``\\mathcal{F} \\subset \\mathbb{R}``.
"""
function unlink end

"""
    jacobian_link(args...)

Evaluates the derivative of the [`link`](@ref link) with respect to the parameter ``f``.
"""
function jacobian_link end

"""
    IdentityLink <: Link

Define the map ``\\tilde{f} = f`` where ``f \\in \\mathbb{R}`` and ``\\tilde{f} \\in \\mathbb{R}``
"""
struct IdentityLink <: Link end

link(::Type{IdentityLink}, param::T) where T = param
unlink(::Type{IdentityLink}, param_tilde::T) where T = param_tilde
jacobian_link(::Type{IdentityLink}, param::T) where T = one(T)

"""
    LogLink <: Link

Define the map ``\\tilde{f} = \\ln(f - a)`` where ``f \\in [a, \\infty), a \\in \\mathbb{R}`` and ``\\tilde{f} \\in \\mathbb{R}``
"""
struct LogLink <: Link end

link(::Type{LogLink}, param::T, lower_bound::T) where T = log(param - lower_bound)
unlink(::Type{LogLink}, param_tilde::T, lower_bound::T) where T = exp(param_tilde) + lower_bound
jacobian_link(::Type{LogLink}, param::T, lower_bound::T) where T = 1/(param - lower_bound)

"""
    LogitLink <: Link

Define the map ``\\tilde{f} = \\ln(\\frac{f - a}{b - f})`` where ``f \\in [a, b], a, b \\in \\mathbb{R}`` and ``\\tilde{f} \\in \\mathbb{R}``
"""
struct LogitLink <: Link end

function link(::Type{LogitLink}, param::T, lower_bound::T, upper_bound::T) where T
    return log((param - lower_bound)/(upper_bound - param))
end
function unlink(::Type{LogitLink}, param_tilde::T, lower_bound::T, upper_bound::T) where T
    return lower_bound + ((upper_bound - lower_bound)/(1 + exp(-param_tilde)))
end
function jacobian_link(::Type{LogitLink}, param::T, lower_bound::T, upper_bound::T) where T
    return (upper_bound + lower_bound)/((upper_bound - param) * (param - lower_bound))
end

const LINKS = [
    IdentityLink;
    LogLink;
    LogitLink
]


text:
module Roofline

using LinuxPerf

@noinline function escape(x)
    Base.inferencebarrier(nothing)::Nothing
end

struct RooflineBench
    counters::Vector{LinuxPerf.Counter}
    times::Vector{UInt64}
    repeats::Int64
end

function benchmark(func::F, setup::S, args...; repeats=10) where {F, S}
    counters = LinuxPerf.Counter[]
    times = UInt64[]
    val = func(setup(args...)...)
    escape(val)

    for benchPerf in (intel_roofline_bw_use(),
                      intel_roofline_memory(),
                      intel_roofline_double(),
                      intel_roofline_single(),
                      )
        data = setup(args...)
        for i in 1:repeats
            LinuxPerf.enable!(benchPerf)
            start = Base.time_ns()
            val = func(data...)
            stop = Base.time_ns()
            LinuxPerf.disable!(benchPerf)
            escape(val)
            push!(times, stop-start)
        end
        append!(counters, LinuxPerf.counters(benchPerf).counters)
    end
    return RooflineBench(counters, times, repeats)
end

function Base.show(io::IO, bench::RooflineBench)
   t = (sum(bench.times)/length(bench.times))*10e-9
   show(io, typeof(bench))
   println(io)
   show(io, "Mean duration $(t)s")
   println(io)
   results = summarize_intel(bench)
   show(io, "Double GFLOP/s $(results.dflops/t/bench.repeats * 1e-9) ")
   println(io)
   show(io, "Single GFLOP/s $(results.sflops/t/bench.repeats * 1e-9)")
   println(io)
   show(io, "DRAM BW GB/s $(64*results.bw_ops / t / bench.repeats* 1e-9)") # cachelines to GB/s
   println(io)
   show(io, "MemOP per FLOP: $(results.mops / (results.sflops + results.dflops))")
   println(io)
end

function summarize_intel(bench::RooflineBench)
    dflops = 0
    dfinst = 0
    sflops = 0
    sfinst = 0
    mops = 0
    bw_ops = 0
    for counter in bench.counters
        prop = intel_properties[counter.event]
        if counter.enabled == 0 || counter.running == 0
            @info "Not counting" counter
            continue
        end
        if prop.kind === :double
            dflops += prop.mul * Int64(counter.value)
            dfinst += prop.ops * Int64(counter.value)
        elseif prop.kind === :single
            sflops += prop.mul * Int64(counter.value)
            sfinst += prop.ops * Int64(counter.value)
        elseif prop.kind === :mem
            mops += Int64(counter.value)
        elseif prop.kind === :bw
            bw_ops += Int64(counter.value)
        else
            @error "Unkown" counter prop
        end
    end
    return (;dflops=dflops, dfinst=dfinst, sflops=sflops, sinst=sfinst, mops=mops, bw_ops=bw_ops)
end

# Reading:
# - https://crd.lbl.gov/assets/Uploads/CS267-2019-Roofline-SWWilliams.pdf
# - https://www.mjr19.org.uk/sw/mflops/
# - https://software.intel.com/sites/default/files/managed/8b/6e/335279_performance_monitoring_events_guide.pdf

import LinuxPerf: EventType, enable!, disable!, counters

const intel_properties = Dict(
    EventType(:raw, :fp_arith_inst_retired_scalar_double) => (;mul=1, vlen=1, ops=1, kind=:double),
    EventType(:raw, :fp_arith_inst_retired_128B_packed_double) => (;mul=2, vlen=2, ops=1, kind=:double),
    EventType(:raw, :fp_arith_inst_retired_256B_packed_double) => (;mul=4, vlen=4, ops=1, kind=:double),
    # EventType(:raw, :fp_arith_inst_retired_512B_packed_double) => (;mul=8, vlen=8, ops=1, kind=:double),
    EventType(:raw, :fp_arith_inst_retired_scalar_single) => (;mul=1, vlen=1, ops=1, kind=:single),
    EventType(:raw, :fp_arith_inst_retired_128B_packed_single) => (;mul=4, vlen=2, ops=1, kind=:single),
    EventType(:raw, :fp_arith_inst_retired_256B_packed_single) => (;mul=8, vlen=4, ops=1, kind=:single),
    # EventType(:raw, :fp_arith_inst_retired_512B_packed_single) => (;mul=16, vlen=8, ops=1, kind=:single),
    EventType(:raw, :mem_inst_retired_all_loads) => (;mul=1, vlen=1, ops=1, kind=:mem),
    EventType(:raw, :mem_inst_retired_all_stores) => (;mul=1, vlen=1, ops=1, kind=:mem),
    EventType(:raw, :unc_arb_trk_requests_all) => (;mul=1, vlen=1, ops=1, kind=:bw),
    EventType(:raw, :unc_arb_coh_trk_requests_all) =>  (;mul=1, vlen=1, ops=1, kind=:bw),
)

# Can't schedule them all thogether, so seperate double from single from memory.
const intel_roofline_events_double = [
    [
        EventType(:raw, :fp_arith_inst_retired_scalar_double),
        EventType(:raw, :fp_arith_inst_retired_128B_packed_double),
        EventType(:raw, :fp_arith_inst_retired_256B_packed_double),
        # EventType(:raw, :fp_arith_inst_retired_512B_packed_double),
    ],
]

const intel_roofline_events_single = [
    [
        EventType(:raw, :fp_arith_inst_retired_scalar_single),
        EventType(:raw, :fp_arith_inst_retired_128B_packed_single),
        EventType(:raw, :fp_arith_inst_retired_256B_packed_single),
        # EventType(:raw, :fp_arith_inst_retired_512B_packed_single),
    ]
]

const intel_roofline_events_memory = [
    [
        EventType(:raw, :mem_inst_retired_all_loads),
        EventType(:raw, :mem_inst_retired_all_stores)
    ]
]

# https://github.com/andikleen/pmu-tools/blob/bb9e7429903e3a9c9a14c951b1013aa991bf119c/skl_client_ratios.py#L348
# Average external Memory Bandwidth Use for reads and writes [GB / sec]
# def DRAM_BW_Use(self, EV, level):
#    return 64 *(EV("UNC_ARB_TRK_REQUESTS.ALL", level) + EV("UNC_ARB_COH_TRK_REQUESTS.ALL", level)) / OneMillion / EV("interval-s", 0) / 1000

const intel_roofline_events_bw_use = [
    [
        EventType(:raw, :unc_arb_trk_requests_all),
        EventType(:raw, :unc_arb_coh_trk_requests_all)
    ]
]

function intel_roofline_double()
    LinuxPerf.make_bench(intel_roofline_events_double)
end
function intel_roofline_single()
    LinuxPerf.make_bench(intel_roofline_events_single)
end

function summarize_intel(c)
    flops = 0
    finst = 0
    mops = 0
    for counter in c.counters
        prop = intel_properties[c.event]
        if counter.enable == 0 || c.running == 0
            @info "Not counting" c
            continue
        end
        if prop.kind === :double
            flops += prop.mul * Int64(counter.value)
            finst += prop.ops * Int64(counter.value)
        elseif prop.kind === :mem
            mops += Int64(counter.value)
        else
            @error "Unkown" c prop
        end
    end
    return (;flops=flops, finst=finst, mops=mops)
end

end # module


text:
using Plots
using CSV
using Printf
using DataFrames
using LaTeXStrings
using LsqFit            # Required for the statistics on the fitting process
using FrictionModels

# Read in a CSV file containing the data to fit.
#
# The file is structured as follows:
# Column 1: Time
# Column 2: Velocity
# Column 3: Normal force
# Column 4: Friction Force
println(pwd())
data = CSV.read("examples//test_data_1.csv", DataFrame)

# Construct an initial estimate for model parameters
init_guess = HyperbolicModel(
    0.15,
    2.0,
    20.0,
    1.0,
    0.1,
    0.1
)

# Fit the model
results = fit_model(
    init_guess,     # Initial Guess
    data[:,1],      # Time
    data[:,4],      # Friction Force
    data[:,3],      # Normal Force
    data[:,2]       # Velocity
)

# Solve the model
F = friction.(results.model, data[:,3], data[:,2])

# Plot the data
plt = plot(
    data[:,2], F,
    xlabel = L"v(t)",
    ylabel = L"F(t)",
    label = "Fitted Model",
    lw = 2,
    legend = :bottomright
)
plot!(
    plt,
    data[:,2], data[:,4],
    st = :scatter,
    label = "Measured Data"
)
display(plt)

# Print out the coefficients of the fitted model
@printf("Fitted Model Coefficients:\n")

# Print out statistitics related to the fit
sigma = stderror(results.fit)
confidence = confidence_interval(results.fit)

# Print out the results
@printf("friction_coefficient: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.friction_coefficient,
    sigma[1],
    confidence[1][1],
    confidence[1][2]
)
@printf("normalization_coefficient: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.normalization_coefficient,
    sigma[2],
    confidence[2][1],
    confidence[2][2]
)
@printf("dissipation_coefficient: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.dissipation_coefficient,
    sigma[3],
    confidence[3][1],
    confidence[3][2]
)
@printf("hysteresis_coefficient: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.hysteresis_coefficient,
    sigma[4],
    confidence[4][1],
    confidence[4][2]
)
@printf(
    "stribeck_velocity: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.stribeck_velocity,
    sigma[5],
    confidence[5][1],
    confidence[5][2]
)
@printf(
    "viscous_damping: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.viscous_damping,
    sigma[6],
    confidence[6][1],
    confidence[6][2]
)


text:
###
### iterate
###
checkindexlo(r, i::AbstractVector) = checkindexlo(r, minimum(i))
checkindexlo(r, i) = firstindex(r) <= i
checkindexlo(r, i::CartesianIndex{1}) = firstindex(r) <= first(i.I)

checkindexhi(r, i::AbstractVector) = checkindexhi(r, maximum(i))
checkindexhi(r, i) = lastindex(r) >= i
checkindexhi(r, i::CartesianIndex{1}) = firstindex(r) <= first(i.I)

###
### Generic array traits
###
# TODO should this be in Static.jl
static_isempty(x::OrdinalRange) = _static_isempty(static_first(x), static_step(x), static_last(x))
function static_isempty(x)
    len = static_length(x)
    return Static.eq(zero(len), len)
end
function _static_isempty(start::F, step::S, stop::L) where {F,S,L}
    return Static.ne(start, stop) & Static.ne(Static.gt(step, zero(step)), Static.gt(stop, start))
end

# FIXME these absolutely needs to go in ArrayInterface
function ArrayInterface.known_length(::Type{T}) where {T<:AbstractRange}
    if parent_type(T) <: T
        return nothing
    else
        return known_length(parent_type(T))
    end
end

#Base.rem(::Static.StaticFloat64{X}, ::Static.StaticFloat64{Y}) where {X,Y} = static(rem(X, Y))
static_div(x::X, y::Y) where {X,Y} = _div(is_static(X) & is_static(Y), x, y)
_div(::True, x, y) = static(div(known(x), known(y)))
_div(::False, x, y) = div(dynamic(x), dynamic(y))

static_rem(x::X, y::Y) where {X,Y} = _rem(is_static(X) & is_static(Y), x, y)
_rem(::True, x, y) = static(rem(known(x), known(y)))
_rem(::False, x, y) = rem(dynamic(x), dynamic(y))

_sub1(x::T) where {T} = x - oneunit(T)
_add1(x::T) where {T} = x + oneunit(T)
_int(idx) = round(Integer, idx, RoundToZero)::Int
_int(idx::Integer) = Int(idx)::Int
_int(idx::StaticInt{N}) where {N} = idx
_int(idx::TwicePrecision{T}) where {T} = round(Integer, T(idx), RoundToZero)

_drop_unit(x::X) where {X} = div(x, oneunit(x))
_drop_unit(x::Real) = x

### empty
_empty_ur(::Type{T}) where {T} = one(T):zero(T)

_empty(x::X, y::Y) where {X,Y} = Vector{Int}()
@inline function _empty(x::X, y::Y) where {X<:AbstractRange,Y<:AbstractRange}
    if (known_step(X) === nothing) | (known_step(Y) === nothing)
        return 1:1:0
    else
        if known_first(x) === one(eltype(X))  && known_first(y) === one(eltype(Y))
            if known_last(x) isa Nothing || known_last(y) isa Nothing
                return static(1):0
            else
                return static(1):static(0)
            end
        else
            return 1:0
        end
    end
end

const Equal{T} = Union{Fix2{typeof(==),T},Fix2{typeof(isequal),T}}
const NotEqual{T} = Fix2{typeof(!=),T}
const NotIn{T} = (typeof(!in(Any)).name.wrapper){Fix2{typeof(in),T}}
const In{T} = Fix2{typeof(in),T}

_maybe_static(::True, x::Int) = static(x)
_maybe_static(::True, x::StaticInt) = x
_maybe_static(::False, x::Int) = x
_maybe_static(::False, x::StaticInt) = dynamic(x)

ifelseop(::True, t, f, args...) = t(args...)
ifelseop(::False, t, f, args...) = f(args...)
@inline function ifelseop(b::Bool, t, f, args...)
    if b
        return dynamic(t(args...))
    else
        return dynamic(f(args...))
    end
end

return_nothing(args...) = nothing

return_static_one(args...) = static(1)

ifadd1(x, val) = ifelseop(x, _add1, identity, val)
ifsub1(x, val) = ifelseop(x, _sub1, identity, val)

_static_length(_, r) = static_length(r)
_static_length(_, _, r) = static_length(r)


text:
export HierarchicalOperator

##
# Represent a binary hierarchical operator
# [ x x \ /
#   x x / \
#   \ / x x
#   / \ x x ]
# where the diagonal blocks are of one type, and
# the off-diagonal blocks are of a type with low-rank structure.
#
# To have simple recursive construction, we use the ordering:
# [ D₁   L₄
#      ⋯      L₂
#   L₃   D₂
#          D₃   L₆
#      L₁     ⋯
#          L₅   D₄ ]
# The diagonal blocks are ordered along the diagonal,
# and the off-diagonal blocks are ordered from: bottom left,
# top right, then followed recursively by top left and bottom right.
##

mutable struct HierarchicalOperator{S,V,T,HS,HV} <: Operator{T}
    diagonaldata::HS
    offdiagonaldata::HV

    A::Matrix{T} # Cache of matrix for pivot computation
    factorization::PivotLDU{T,Matrix{T}} # Cache of factorization of A for pivot computation
    factored::Bool

    function HierarchicalOperator{S,V,T,HS,HV}(diagonaldata::HS,offdiagonaldata::HV) where {S,V,T,HS,HV}
        H = new{S,V,T,HS,HV}()
        H.diagonaldata = diagonaldata
        H.offdiagonaldata = offdiagonaldata

        r = mapreduce(rank,+,offdiagonaldata)
        H.A = eye(T,r,r)
        H.factored = false

        H
    end
end

HierarchicalOperator(diagonaldata::Tuple{S1,S2},offdiagonaldata::Tuple{V1,V2}) where {S1,S2,V1,V2} = HierarchicalOperator{promote_type(S1,S2),promote_type(V1,V2),promote_type(eltype(S1),eltype(S2),eltype(V1),eltype(V2)),Tuple{S1,S2},Tuple{V1,V2}}(diagonaldata,offdiagonaldata)
HierarchicalOperator(diagonaldata::Tuple{HierarchicalOperator{S1,V1,T1,HS1,HV1},HierarchicalOperator{S2,V2,T2,HS2,HV2}},offdiagonaldata::Tuple{V3,V4}) where {S1,S2,V1,V2,V3,V4,T1,T2,HS1,HS2,HV1,HV2} = HierarchicalOperator{promote_type(S1,S2),promote_type(V1,V2,V3,V4),promote_type(eltype(S1),eltype(S2),eltype(V1),eltype(V2),eltype(V3),eltype(V4),T1,T2),Tuple{HierarchicalOperator{S1,V1,T1,HS1,HV1},HierarchicalOperator{S2,V2,T2,HS2,HV2}},Tuple{V3,V4}}(diagonaldata,offdiagonaldata)
HierarchicalOperator(diagonaldata::Tuple{S1,HierarchicalOperator{S,V,T,HS,HV}},offdiagonaldata::Tuple{V1,V2}) where {S,S1,V,V1,V2,T,HS,HV} = HierarchicalOperator{promote_type(S,S1),promote_type(V,V1,V2),promote_type(eltype(S),eltype(S1),eltype(V),eltype(V1),eltype(V2),T),Tuple{S1,HierarchicalOperator{S,V,T,HS,HV}},Tuple{V1,V2}}(diagonaldata,offdiagonaldata)
HierarchicalOperator(diagonaldata::Tuple{HierarchicalOperator{S,V,T,HS,HV},S1},offdiagonaldata::Tuple{V1,V2}) where {S,S1,V,V1,V2,T,HS,HV} = HierarchicalOperator{promote_type(S,S1),promote_type(V,V1,V2),promote_type(eltype(S),eltype(S1),eltype(V),eltype(V1),eltype(V2),T),Tuple{HierarchicalOperator{S,V,T,HS,HV},S1},Tuple{V1,V2}}(diagonaldata,offdiagonaldata)

HierarchicalOperator(diagonaldata::Vector,offdiagonaldata::Vector)=HierarchicalOperator(diagonaldata,offdiagonaldata,round(Int,log2(length(diagonaldata))))

function HierarchicalOperator(diagonaldata::Vector,offdiagonaldata::Vector,n::Int)
    @assert length(diagonaldata) == 2^n
    @assert length(offdiagonaldata) == 2^(n+1)-2

    if n == 1
        return HierarchicalOperator(tuple(diagonaldata...),tuple(offdiagonaldata...))
    elseif n ≥ 2
        dldp2 = div(length(offdiagonaldata)+2,2)
        return HierarchicalOperator((HierarchicalOperator(diagonaldata[1:2^(n-1)],offdiagonaldata[3:dldp2],n-1),
                                   HierarchicalOperator(diagonaldata[1+2^(n-1):end],offdiagonaldata[dldp2+1:end],n-1)),
                                   (offdiagonaldata[1],offdiagonaldata[2]))
    end
end

Base.similar(H::HierarchicalOperator) = HierarchicalOperator(map(similar,diagonaldata(H)),map(similar,offdiagonaldata(H)))
Base.similar(H::HierarchicalOperator{SS,V,T}, S) where {SS,V,T} = HierarchicalOperator(map(A->similar(A,S),diagonaldata(H)),map(A->similar(A,S),offdiagonaldata(H)))

diagonaldata(H::HierarchicalOperator) = H.diagonaldata
offdiagonaldata(H::HierarchicalOperator) = H.offdiagonaldata

nlevels(H::HierarchicalOperator) = 1+mapreduce(nlevels,max,diagonaldata(H))

partition(H::HierarchicalOperator) = diagonaldata(H),offdiagonaldata(H)

collectoffdiagonaldata(H::HierarchicalOperator{S,V,T,NTuple{2,S}}) where {S,V,T} = collect(offdiagonaldata(H))
function collectoffdiagonaldata(H::HierarchicalOperator{S,V,T,HS}) where {S,V,T,HS}
    data = collect(offdiagonaldata(H))
    append!(data,mapreduce(collectoffdiagonaldata,vcat,diagonaldata(H)))
    data
end

collectdiagonaldata(H::HierarchicalOperator{S,V,T,NTuple{2,S}}) where {S,V,T} = collect(diagonaldata(H))
function collectdiagonaldata(H::HierarchicalOperator{S,V,T,HS}) where {S,V,T,HS}
    data = S[]
    append!(data,mapreduce(collectdiagonaldata,vcat,diagonaldata(H)))
    data
end

isfactored(H::HierarchicalOperator) = H.factored

function condest(H::HierarchicalOperator)
    !isfactored(H) && factorize!(H)
    return cond(H.A)*mapreduce(condest,+,diagonaldata(H))
end

convert(::Type{HierarchicalOperator{S,V,T,HS,HV}},M::HierarchicalOperator) where {S,V,T,HS,HV} = HierarchicalOperator(convert(Vector{S},collectdiagonaldata(M)),convert(Vector{V},collectoffdiagonaldata(M)))
Base.promote_rule(::Type{HierarchicalOperator{S,V,T,HS,HV}},::Type{HierarchicalOperator{SS,VV,TT,HSS,HVV}}) where {S,V,T,HS,HV,SS,VV,TT,HSS,HVV}=HierarchicalOperator{promote_type(S,SS),promote_type(V,VV),promote_type(T,TT),promote_type(HS,HSS),promote_type(HV,HVV)}

Base.transpose(H::HierarchicalOperator) = HierarchicalOperator(map(transpose,diagonaldata(H)),map(transpose,reverse(offdiagonaldata(H))))
Base.ctranspose(H::HierarchicalOperator) = HierarchicalOperator(map(ctranspose,diagonaldata(H)),map(ctranspose,reverse(offdiagonaldata(H))))

Base.copy(H::HierarchicalOperator) = HierarchicalOperator(map(copy,diagonaldata(H)),map(copy,offdiagonaldata(H)))
Base.copy!(H::HierarchicalOperator,J::HierarchicalOperator) = (map(copy!,diagonaldata(H),diagonaldata(J));map(copy!,offdiagonaldata(H),offdiagonaldata(J));H)

rank(A::Operator) = ∞

function blockrank(H::HierarchicalOperator)
    m,n = blocksize(H)
    A = Array{Any}(undef,m,n)
    (m1,n1),(m2,n2) = map(blocksize,diagonaldata(H))
    r1,r2 = map(rank,offdiagonaldata(H))
    for j=1:n1,i=1:m2
        A[i+m1,j] = r1
    end
    for j=1:n2,i=1:m1
        A[i,j+n1] = r2
    end
    A11,A22 = map(blockrank,diagonaldata(H))
    if typeof(A11) <: Infinity{Bool}
        for j=1:n1,i=1:m1
            A[i,j] = ∞
        end
    else
        for j=1:n1,i=1:m1
            A[i,j] = A11[i,j]
        end
    end
    if typeof(A22) <: Infinity{Bool}
        for j=1:n2,i=1:m2
            A[i+m1,j+n1] = ∞
        end
    else
        for j=1:n2,i=1:m2
            A[i+m1,j+n1] = A22[i,j]
        end
    end
    A
end

for OP in (:domainspace,:rangespace)
    @eval $OP(H::HierarchicalOperator{S,V}) where {S<:Operator,V<:AbstractLowRankOperator}=PiecewiseSpace(map($OP,diagonaldata(H))...)
end

blocksize(H::HierarchicalOperator{U,V}) where {U<:Operator,V<:AbstractLowRankOperator} =
    ncomponents(domainspace(H)),ncomponents(rangespace(H)) # TODO: check it's not rangespace...domainspace

for op in (:+,:-)
    @eval begin
        $op(H::HierarchicalOperator) = HierarchicalOperator(map($op,diagonaldata(H)),map($op,offdiagonaldata(H)))
        $op(H::HierarchicalOperator,J::HierarchicalOperator) = HierarchicalOperator(map($op,diagonaldata(H),diagonaldata(J)),map($op,offdiagonaldata(H),offdiagonaldata(J)))
        #$op(H::HierarchicalOperator,L::LowRankMatrix) = $op(promote(H,L)...)
        #$op(L::LowRankMatrix,H::HierarchicalOperator) = $op(promote(L,H)...)
        #$op(H::HierarchicalOperator,A::Operator) = $op(promote(H,A)...)
        #$op(A::Operator,H::HierarchicalOperator) = $op(promote(A,H)...)
    end
end

#=
function *(H::HierarchicalOperator, x::HierarchicalVector)
    T = promote_type(eltype(H),eltype(x))
    mul!(similar(x,T),H,x)
end

function mul!(b::AbstractVector,D::Base.LinAlg.Diagonal,x::AbstractVector)
    d = D.diag
    for i=1:min(length(b),length(x))
        @inbounds b[i] += d[i]*x[i]
    end
    b
end

function mul!(b::HierarchicalVector,H::HierarchicalOperator,h::HierarchicalVector)
    (H11,H22),(H21,H12) = partition(H)
    h1,h2 = partition(h)
    b1,b2 = partition(b)
    mul!(b1,H12,h2)
    mul!(b1,H11,h1)
    mul!(b2,H21,h1)
    mul!(b2,H22,h2)
    b
end

function *(H::HierarchicalOperator,b::Vector)
    (H11,H22),(H21,H12) = partition(H)
    m1,m2 = size(H12,1),size(H21,1)
    (b1,b2) = (b[1:m1],b[1+m1:m1+m2])
    HierarchicalVector((H11*b1⊕H12*b2,H21*b1⊕H22*b2))
end
=#


text:
function noisy_circle_data(N, σ)
    r = 1 .+ σ*randn(1,N)
    θ = rand(Uniform(0, 2π), (1,N))

    rcosθ = r .* cos.(θ)
    rsinθ = r .* sin.(θ)

    data = [rcosθ; rsinθ]

    return data
end


text:
function plot_formation(
    data::TrajectoryData, freq::ℝ, ctrl::FormationControl, formation_at_t
)
    scene, layout = layoutscene(resolution = (1600, 1600))

    ax_traj = layout[1,1] = LAxis(
        scene, title = "Trajectories", aspect = DataAspect(),
        backgroundcolor = RGBf0(1,1,1))

    times = data.times
    time_ls = labelslider!(
        scene, "time", times;
        format = x -> "$(round((x - 1) / freq, digits = 2))s")
    layout[2,1] = time_ls.layout
    t = time_ls.slider.value

    N = size(data.values, 2)
    colors = [get(ColorSchemes.Spectral_10, i/N) for i in 1:N]

    function plot_snapshot(t)
        # draw cars
        for id in 1:N
            car = @lift let
                x, y, θ = map(l -> data[l][$t,id], ("x", "y", "θ"))
                car_triangle(x, y, θ)
            end
            poly!(ax_traj, car; color = colors[id], strokecolor=:black, strokewidth=2)
        end
        # draw ref points
        refpoints = @lift let
            x, y, θ, ψ, v = map(l -> data[l][$t,1], ("x", "y", "θ", "ψ", "v"))
            leader = CarX(;x, y, θ, ψ, v)
            @_ (formation_at_t($t)
                |> map(ref_point(ctrl.K, _),__)
                |> map(to_formation_frame(ctrl, leader), __))
        end
        points = @lift (x -> Point2f0(x[1])).($refpoints)
        directions = @lift (x -> Point2f0(x[2])).($refpoints)
        arrows!(ax_traj, points, directions;
            linecolor=:red, arrowcolor=:red, lengthscale=0.15, arrowsize=0.10, linewidth=1)
        scatter!(ax_traj, points; color=:red)
    end

    # draw trajectories
    for id in 1:N
        xs, ys, θs = map(l -> data[l][:,id], ("x", "y", "θ"))
        lines!(ax_traj, xs, ys; color=colors[id], linewidth=4)
    end
    # plot current snapshot
    plot_snapshot(t)
    # plot additional snapshots (used for paper figures)
    interesting_times = [] # [386, 833, 963, 1254, 1610]
    foreach(plot_snapshot, Observable.(interesting_times))

    scene
end

struct FormationObsDynamics{Fu,Ff} <: ObsDynamics
    "`external_u(x, t) -> u`"
    external_u::Fu
    "`external_formation(x, t) -> formation_idx`"
    external_formation::Ff
end

function OneVision.obs_forward(dy::FormationObsDynamics, x, z, t::𝕋)
    HVec(dy.external_u(x, t), dy.external_formation(x, t))
end

"""
Compute the time-averaged distance between each follower and their expected
formation position using a quadratic norm. Mathematically, this metric is defined as

     √(1/T ∫ ∑_i |x_i - x⋆_i|^2 dt) .
"""
function avg_formation_deviation(
    data::TrajectoryData, ctrl::FormationControl, formation_at_t
)::ℝ
    function deviation_at(t)
        x, y, θ, ψ, v = map(l -> data[l][t, 1], ("x", "y", "θ", "ψ", "v"))
        leader = CarX(;x, y, θ, ψ, v)
        N = size(data.values, 2)
        form = formation_at_t(t)[2:end]
        @asserteq length(form) (N-1)
        ref_pos =
            @_ (form
            |> map(ref_point(ctrl.K, _),__)
            |> map(to_formation_frame(ctrl, leader), __)
            |> map(_.pos, __))
        x_data = data["x"][t, :]
        y_data = data["y"][t, :]
        actual_pos = [@SVector[x_data[i], y_data[i]] for i in 2:N]
        mean(norm.(actual_pos .- ref_pos))
    end

    mean(deviation_at.(eachindex(data.times)))
end

function formation_example(;
        setting::ExampleSetting,
        CF::CFName = onevision_cf,
        switch_formation = true,
        seed = 1234,
        leader_noise = false,
        track_config = false,
        plot_result = true,
    )
    @unpack time_end, freq, noise, sensor_noise, delays, H, model_error = setting
    X, U = CarX{ℝ}, CarU{ℝ}
    Z = HVec{U, ℕ}
    t_end = 𝕋(ceil(time_end * freq))
    delta_t = 1 / freq

    # Car dynamics parameters
    max_ψ = 45°
    dy_model = CarDynamics(l = 0.1 * (1+model_error) ;delta_t)
    rng = MersenneTwister(seed)
    function add_noise(x::X, t)::X where X
        x + CarX(x=0.0, y=0.0, θ=0.0,
                v = randn(rng, ℝ), ψ = randn(rng, ℝ)) * noise
    end
    dy_actual = CarDynamics(; delta_t, max_ψ, add_noise)
    dy_actual_leader = leader_noise ? dy_actual : CarDynamics(; delta_t, max_ψ)
    xs_observer =
        if leader_noise
            (xs, t) -> @_ (_ + randn(rng, X) * sensor_noise).(xs)
        else
            (xs, t) ->
                [i == 1 ? x : x + randn(rng, X) * sensor_noise
                for (i, x) in enumerate(xs)]
        end

    function external_control(_, t)
        (if t ≤ 3 * freq
            U(v̂ = 1.0, ψ̂ = 0.0)
        elseif t ≤ 5 * freq
            U(v̂ = 1.0, ψ̂ = 4°)
        elseif t ≤ 9 * freq
            U(v̂ = 1.0, ψ̂ = 0°)
        elseif t ≤ 15 * freq
            U(v̂ = 2.0, ψ̂ = 2°)
        else
            U(v̂ = 2(1 - (t/freq-15)/5), ψ̂ = 2°)
        end) + U(v̂ = randn(rng, ℝ), ψ̂ = randn(rng, ℝ)) * noise/4
    end

    function external_formation(_, t)
        if t ≥ 8 * freq && switch_formation
            2
        else
            1
        end
    end

    leader_z_dy = FormationObsDynamics(external_control, external_formation)

    # running at 20Hz
    delays_model  = delays
    delays_actual = delays_model
    ΔT = delays_model.ΔT

    N = 4
    triangle_formation = let
        l = 0.8
        Δϕ = 360° / (N-1)
        circle = [X(x = l * cos(Δϕ * i), y = l * sin(Δϕ * i), θ = 0) for i in 1:N-1]
        [[zero(X)]; circle]
    end
    vertical_formation = let
        l = 0.6
        leader_idx = round_ceil(N/2)
        line = [X(x = l * (i - leader_idx), y = 0, θ = 0) for i in 1:N]
        [line[mod1(leader_idx + j - 1, N)] for j in 1:N]
    end

    formations = [triangle_formation, vertical_formation]
    form_from_id(i) = formations[i]

    RefK = if track_config
        ConfigTrackControl(dy_model, 1.0, 1.0, 1.0)
    else
        RefPointTrackControl(;
            dy = dy_model, ref_pos = dy_model.l, ctrl_interval = delta_t * ΔT,
            kp = 1.0, ki = 0.0, kd = 0.0)
    end
    avoidance = CollisionAvoidance(scale=1.0, min_r=dy_model.l, max_r=3*dy_model.l)
    central = FormationControl(
        (_, zs, _) -> form_from_id(zs[1].d),
        RefK, dy_model, avoidance)

    init = let
        formation = rotate_formation(form_from_id(1), 90°)
        map(1:N) do i
            x = formation[i]
            x, HVec(zero(U), 1), zero(U)
        end
    end

    world_model = WorldDynamics(fill((dy_model, StaticObsDynamics()), N))
    loss_model = let
        x_weights = SVector{N}(fill(X(x = 100, y = 100, θ = 1), N))
        u_weights = SVector{N}(fill(U(v̂ = 1, ψ̂ = 10), N))
        RegretLossModel(central, world_model, x_weights, u_weights)
    end
    framework = mk_cf(CF, world_model, central, delays, loss_model; X, Z, H)

    comps = ["x", "y", "θ", "ψ", "v"]
    function record_f(xs, zs, us)
        comps = [:x, :y, :θ, :ψ, :v]
        hcat(((p -> getfield(p,c)).(xs) for c in comps)...)
    end

    formation_at_t(t) = form_from_id(external_formation(missing, t))

    world = WorldDynamics([
        (dy_actual_leader, leader_z_dy)
        fill((dy_actual, StaticObsDynamics()), N-1)])
    result, (logs, loss) = simulate(
        world, delays_actual, framework, init, (comps, record_f), 1:t_end
        ; xs_observer, loss_model)
    if plot_result
        visualize(result; delta_t = 1 / freq, loss) |> display
        plot_formation(result, freq, central, formation_at_t) |> display
    end
    metrics = Dict(
        "avg deviation" => avg_formation_deviation(result, central, formation_at_t))
    loss, metrics
end


text:
using Test, FFAST

@testset "Iron" begin
    @test length(eachedge(26)) == 9

    @test hasedge(26, 1)
    @test hasedge(26, 2)
    @test hasedge(26, 3)
    @test hasedge(26, 4)
    @test hasedge(26, 5)
    @test hasedge(26, 6)
    @test hasedge(26, 7)
    @test hasedge(26, 8)
    @test hasedge(26, 9)
    @test !hasedge(26, 10)

    @test !(10 in eachedge(26))
    @test 2 in eachedge(26)
    @test 7 in eachedge(26)

    @test edgeenergy(26, 1) == 7112.0
    @test edgeenergy(26, 9) == 3.6
    @test edgeenergy(26, 4) ≈ 708.1

    @test isapprox(atomicweight(26), 55.845, atol = 0.003)

    @test density(26) ≈ 7.860

    @test EV(26) ≈ 7.53493E+05

    @test crosssectionfactor(26) ≈ 9.27362E+01 * 1.0e-24

    @test relativisticcorrection(26)[1] ≈ -1.0812E-01
    @test relativisticcorrection(26)[2] ≈ -6.7800E-02

    @test nuclearthompsoncorrection(26) ≈ -6.6403E-03

    @test isapprox(mac(PhotoElectricMAC, 26, 30.0), 106592.1003, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 102.0), 46393.2381, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 304.3), 10243.4145, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 1234.1), 5297.2337, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 2342.0), 1013.5935, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 4523.0), 167.5849, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 12340.0), 95.1269, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 34000.0), 5.2825, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 9630.0), 186.2521, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 201000.0), 0.0233, rtol = 0.001)
end

@testset "Oxygen" begin
    @test length(eachedge(8)) == 4

    @test hasedge(8, 1)
    @test hasedge(8, 2)
    @test hasedge(8, 3)
    @test hasedge(8, 4)
    @test !hasedge(8, 5)
    @test !hasedge(8, 6)
    @test !hasedge(8, 7)
    @test !hasedge(8, 8)
    @test !hasedge(8, 9)
    @test !hasedge(8, 10)

    @test edgeenergy(8, 1) == 5.32000E-01 * 1000.0
    @test edgeenergy(8, 2) == 2.37000E-02 * 1000.0
    @test edgeenergy(8, 4) ≈ 7.10000E-03 * 1000.0

    @test isapprox(atomicweight(8), 15.99940, atol = 0.003)

    @test density(8) ≈ 1.3310E-03

    @test EV(8) ≈ 2.63012E+06

    @test crosssectionfactor(8) ≈ 2.65676e1 * 1.0e-24

    @test relativisticcorrection(8)[1] ≈ -7.7133E-03
    @test relativisticcorrection(8)[2] ≈ -4.2000E-03

    @test nuclearthompsoncorrection(8) ≈ -2.1944E-03

    @test isapprox(mac(PhotoElectricMAC, 8, 30.0), 303584.6463, rtol = 0.001)

    @test isapprox(mac(PhotoElectricMAC, 8, 30.0), 303584.6463, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 102.0), 49996.7695, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 304.3), 4191.9561, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 1234.1), 2491.7539, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 2342.0), 421.9132, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 4523.0), 59.6995, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 12340.0), 2.7359, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 34000.0), 0.1055, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 9630.0), 6.0225, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 201000.0), 0.00030926, rtol = 0.001)
end

@testset "Lead" begin
    @test length(eachedge(82)) == 23

    @test hasedge(82, 1)
    @test hasedge(82, 2)
    @test hasedge(82, 3)
    @test hasedge(82, 4)
    @test hasedge(82, 5)
    @test hasedge(82, 6)
    @test hasedge(82, 7)
    @test hasedge(82, 8)
    @test !hasedge(82, 23)
    @test !hasedge(82, 25)
    @test hasedge(82, 26)
    @test hasedge(82, 27)

    @test isapprox(mac(PhotoElectricMAC, 82, 30.0), 74494.9013, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 102.0), 11406.5505, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 304.3), 12149.8231, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 1234.1), 3192.3245, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 2342.0), 841.2044, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 4523.0), 893.4338, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 12340.0), 69.5927, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 34000.0), 20.0928, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 9630.0), 134.1960, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 201000.0), 0.8207, rtol = 0.001)
end

@testset "MACpe" begin # generated using NIST DTSA-II
    @test isapprox(mac(PhotoElectricMAC, 1, 1825.6371), 0.9074, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 2, 67143.2841), 2.768340e-5, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 3, 52.8894), 6341.1206, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 4, 323.5202), 13139.1498, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 5, 72.4552), 22406.1014, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 6, 131.0993), 12738.7273, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 7, 54.3167), 111770.5028, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 8, 149.9255), 22397.5204, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 9, 34281.3188), 0.1490, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 10, 202.6272), 23577.7535, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 11, 91714.0014), 0.0129137, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 12, 51.6654), 57373.9013, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 13, 93123.8498), 0.0230, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 14, 80174.5908), 0.0503, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 15, 19484.3456), 5.1816, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 16, 51.6164), 12953.7948, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 17, 92.6523), 14039.6832, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 18, 286446.7751), 0.0019365, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 19, 1282.5060), 2023.8581, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 20, 58.1960), 21012.6801, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 21, 2663.9168), 380.4134, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 22, 52.0073), 12148.6626, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 23, 51.4227), 56667.2105, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 24, 22456.4234), 13.8940, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 25, 11995.0019), 91.0401, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 26, 55.5593), 106613.1135, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 27, 22266.9550), 19.8706, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 28, 205.8346), 22273.6361, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 29, 1039.2688), 9356.0732, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 30, 485.3171), 6755.4381, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 31, 56.2430), 75027.5850, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 32, 164.5878), 47424.7131, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 33, 14923.1883), 97.4613, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 34, 442.2493), 12244.0970, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 35, 30402.3441), 16.1595, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 36, 1763.8442), 4211.7645, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 37, 446.8578), 16086.7643, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 38, 96545.4058), 0.7178, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 39, 822.5155), 5704.5987, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 40, 52279.3424), 5.0292, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 41, 127578.1715), 0.4067, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 42, 18169.7086), 15.0481, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 43, 119.6601), 5432.8402, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 44, 3314.3650), 1661.1260, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 45, 613.7420), 15720.1429, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 46, 297.1361), 2936.5826, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 47, 1798.6100), 1682.0306, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 48, 201700.6027), 0.1723, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 49, 68.2226), 74163.3914, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 50, 227.3590), 6691.7553, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 51, 91067.8017), 1.9996, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 52, 7352.2301), 327.2407, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 53, 271.7124), 4697.0153, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 54, 96.1322), 31927.3174, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 55, 452.1765), 4740.8624, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 56, 2334.6617), 1522.5183, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 57, 103.9157), 39300.3286, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 58, 4223.4834), 392.9113, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 59, 140717.3325), 0.8885, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 60, 261513.7888), 0.1654, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 61, 156.0607), 11951.7922, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 62, 106138.7999), 2.2293, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 63, 189663.2918), 0.4635, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 64, 6110.7096), 204.0632, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 65, 44195.9461), 4.8034, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 66, 20052.2771), 44.3789, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 67, 71.5978), 28290.5339, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 68, 3221.8084), 1245.9590, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 69, 58.3034), 31253.1709, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 70, 84.7028), 26314.9599, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 71, 95.3856), 25954.8160, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 72, 10678.1649), 184.5189, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 73, 86298.9086), 5.8990, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 74, 59.5091), 19804.6389, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 75, 1742.5269), 1093.5987, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 76, 874.9993), 4872.1773, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 77, 778.1277), 6112.9332, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 78, 77.9444), 15733.9354, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 79, 53.5598), 19230.5901, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 80, 313.9529), 15487.9170, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 81, 148.5939), 4791.9895, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 82, 169292.0552), 1.2947, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 83, 158.7399), 4193.2797, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 84, 11247.1253), 97.4316, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 85, 64952.8990), 4.0008, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 86, 24619.4642), 53.9372, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 87, 52.2036), 4222.7289, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 88, 51.5550), 5607.8556, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 89, 1035.3394), 5796.1285, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 90, 142.6819), 6556.7729, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 91, 51.2632), 8270.9798, rtol = 0.001)
    @test isapprox(mac(PhotoElectricMAC, 92, 51.8461), 13886.1866, rtol = 0.001)
end


text:
function transform_density!(img::AbstractArray, edges::AbstractArray, cdf::AbstractArray, minval::Union{Real,AbstractGray}, maxval::Union{Real,AbstractGray})
    first_edge = first(edges)
    step_size = step(edges)
    T = eltype(img)
    map!(img,img) do val
        if val >= edges[end]
            newval = cdf[end]
        elseif val < first_edge
            newval = first(cdf)
        else
            index = Int(Base.div(val-first_edge,step_size)) + 1
            newval = cdf[index]
        end
        # Scale the new intensity value to so that it lies in the range [minval, maxval].
        if T <: Integer
            newval = ceil(minval + ((newval - first(cdf)) * (maxval - minval) / (cdf[end] - first(cdf))))
        else
            newval = minval + ((newval - first(cdf)) * (maxval - minval) / (cdf[end] - first(cdf)))
        end
    end
end

function construct_pdfs(img::AbstractArray, targetimg::AbstractArray, edges::AbstractRange)
    _, histogram = build_histogram(img, edges)
    _, target_histogram = build_histogram(targetimg, edges)
    return edges, histogram / sum(histogram), target_histogram / sum(target_histogram)
end

function construct_pdfs(img::AbstractArray, targetimg::AbstractArray, nbins::Integer = 256)
    if eltype(img) <: AbstractGray
        imin, imax = 0, 1
    else
        imin, imax = min(minfinite(img), minfinite(targetimg)), max(maxfinite(img), maxfinite(targetimg))
    end
    edges, histogram = build_histogram(img, nbins, minval = imin, maxval = imax)
    _, target_histogram = build_histogram(targetimg, edges)
    return edges, histogram / sum(histogram), target_histogram / sum(target_histogram)
end

function lookup_icdf(cdf::AbstractArray, targetcdf::AbstractArray)
    lookup_table = zeros(Int, length(cdf))
    i = 1
    for j = 1:length(cdf)
        p = cdf[j]
        while i < length(targetcdf) && targetcdf[i+1] <= p
            i += 1
        end
        lookup_table[j] = i
    end
    lookup_table
end

function cdf2pdf!(pdf::AbstractArray, cdf::AbstractArray)
    pdf[1] = cdf[1]
    for i = length(cdf)-1:-1:2
        pdf[i] = cdf[i] - cdf[i-1]
    end
end

function linear_stretch(x, A, B, a, b)
    return (x-A) * ((b-a)/(B-A)) + a
end


text:
# ---
# title: 838. Push Dominoes
# id: problem838
# author: zhwang
# date: 2022-02-21
# difficulty: Medium
# categories: Two Pointers, Dynamic Programming
# link: <https://leetcode.com/problems/push-dominoes/description/>
# hidden: true
# ---
#
# There are `N` dominoes in a line, and we place each domino vertically upright.
#
# In the beginning, we simultaneously push some of the dominoes either to the
# left or to the right.
#
# ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/05/18/domino.png)
#
# After each second, each domino that is falling to the left pushes the adjacent
# domino on the left.
#
# Similarly, the dominoes falling to the right push their adjacent dominoes
# standing on the right.
#
# When a vertical domino has dominoes falling on it from both sides, it stays
# still due to the balance of the forces.
#
# For the purposes of this question, we will consider that a falling domino
# expends no additional force to a falling or already fallen domino.
#
# Given a string "S" representing the initial state. `S[i] = 'L'`, if the i-th
# domino has been pushed to the left; `S[i] = 'R'`, if the i-th domino has been
# pushed to the right; `S[i] = '.'`, if the `i`-th domino has not been pushed.
#
# Return a string representing the final state.
#
# **Example 1:**
#
#
#
#     Input: ".L.R...LR..L.."
#     Output: "LL.RR.LLRRLL.."
#
#
# **Example 2:**
#
#
#
#     Input: "RR.L"
#     Output: "RR.L"
#     Explanation: The first domino expends no additional force on the second domino.
#
#
# **Note:**
#
#   1. `0 <= N <= 10^5`
#   2. String `dominoes` contains only `'L`', `'R'` and `'.'`
#
#
## @lc code=start
using LeetCode

function push_dominoes(dominoes::String)::String
    ans = ""
    rules = ("R.L" => "T", ".L" => "LL", "R." => "RR", "T" => "R.L")
    while dominoes != ans
        ans = dominoes
        dominoes = foldl(replace, [dominoes, rules...])
    end
    ans
end
## @lc code=end

## @lc test=start
@testset "838.push-dominoes.jl" begin
    @test push_dominoes("RR.L") == "RR.L"
    @test push_dominoes(".L.R...LR..L..") == "LL.RR.LLRRLL.."
    @test push_dominoes(".L.R.") == "LL.RR"
end
## @lc test=end

## @lc info=start
# link: [solution to 838](https://leetcode-cn.com/problems/push-dominoes/solution/pythonjuliazi-dong-ji-shuang-zhi-zhen-by-d4rb/)
## @lc info=end


text:
# K. A. Fujita, Y. Toyoshima, S. Uda, Y. I. Ozaki, H. Kubota, and S. Kuroda
# Decoupling of receptor and downstream signals in the Akt pathway by its low-pass filter characteristics
# https://doi.org/10.1126/scisignal.2000810

using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Debug)
global_logger(logger)

ode = @ODEmodel(
    EGFR'(t) = EGFR_turnover * pro_EGFR(t) + EGF_EGFR(t) * reaction_1_k2 - EGFR(t) * EGFR_turnover - EGF_EGFR(t) * reaction_1_k1,
    pEGFR'(t) = EGF_EGFR(t) * reaction_9_k1 - pEGFR(t) * reaction_4_k1 + pEGFR_Akt(t) * reaction_2_k2 + pEGFR_Akt(t) * reaction_3_k1 - Akt(t) * pEGFR(t) * reaction_2_k1,
    pEGFR_Akt'(t) = Akt(t) * pEGFR(t) * reaction_2_k1 - pEGFR_Akt(t) * reaction_3_k1 - pEGFR_Akt(t) * reaction_2_k2,
    Akt'(t) = pAkt(t) * reaction_7_k1 + pEGFR_Akt(t) * reaction_2_k2 - Akt(t) * pEGFR(t) * reaction_2_k1,
    pAkt'(t) = pAkt_S6(t) * reaction_5_k2 - pAkt(t) * reaction_7_k1 + pAkt_S6(t) * reaction_6_k1 + pEGFR_Akt(t) * reaction_3_k1 - S6(t) * pAkt(t) * reaction_5_k1,
    S6'(t) = pAkt_S6(t) * reaction_5_k2 + pS6(t) * reaction_8_k1 - S6(t) * pAkt(t) * reaction_5_k1,
    pAkt_S6'(t) = S6(t) * pAkt(t) * reaction_5_k1 - pAkt_S6(t) * reaction_6_k1 - pAkt_S6(t) * reaction_5_k2,
    pS6'(t) = pAkt_S6(t) * reaction_6_k1 - pS6(t) * reaction_8_k1,
    EGF_EGFR'(t) = EGF_EGFR(t) * reaction_1_k1 - EGF_EGFR(t) * reaction_9_k1 - EGF_EGFR(t) * reaction_1_k2,
    y1(t) = a1 * (pEGFR(t) + pEGFR_Akt(t)),
    y2(t) = a2 * (pAkt(t) + pAkt_S6(t)),
    y3(t) = a3 * pS6(t)
)

@time println(assess_identifiability(ode))


text:
@testset "type_utils.jl" begin
    @testset "parameters" begin
        # basic case
        @test collect(parameters(AbstractArray{Float32, 3})) == [Float32, 3]
        # Type-alias
        @test collect(parameters(Vector{Float64})) == [Float64, 1]

        # Tuple
        @test collect(parameters(Tuple{Int8, Bool})) == [Int8, Bool]

        # Tuple with fixed count Vararg
        @test collect(parameters(Tuple{Int8, Vararg{Bool, 3}})) == [Int8, Bool, Bool, Bool]

        # Tuple with varadic Vararg
        a, b = collect(parameters(Tuple{Int8, Vararg{Bool}}))
        @test a == Int8
        @test b <: Vararg{Bool}

        # TypeVar
        tvar1 = parameters(Tuple{T} where T<:Number)[1]
        @test tvar1 isa TypeVar
        @test tvar1.name == :T
        @test tvar1.lb == Union{}
        @test tvar1.ub == Number

        # Shared TypeVar
        tvar2, tvar3 = parameters(Tuple{X,X} where X<:Integer)
        @test tvar2 === tvar3
        @test tvar2.name == :X
        @test tvar2.lb == Union{}
        @test tvar2.ub == Integer

        # Shared TypeVar in different parameter
        tvar4, part = parameters(Tuple{Y,Tuple{Y}} where Integer <: Y <: Real)
        @test part <: Tuple
        tvar5 = parameters(part)[1]
        @test tvar4 === tvar5
        @test tvar4.name == :Y
        @test tvar4.lb == Integer
        @test tvar4.ub == Real

        # Union
        @test Set(parameters(Union{Int8, Bool})) == Set([Int8, Bool])
        @test Set(parameters(Union{Int8, Bool, Set})) == Set([Int8, Bool, Set])

        # Partially collapsing Union
        @test Set(parameters(Union{Int8, Real, Set})) == Set([Real, Set])

        # Unions with type-vars
        umem1, umem2 = parameters(Union{Tuple{Z}, Set{Z}} where Z)
        utvar1 = parameters(umem1)[1]
        utvar2 = parameters(umem2)[1]
        @test utvar1 == utvar2
        @test utvar1 isa TypeVar
        @test utvar1.name == :Z
        @test utvar1.lb == Union{}
        @test utvar1.ub == Any

        # Non-parametric type
        @test isempty(parameters(Bool))
    end
end