text:
function Session(connections = WebSocket[])
    Session(
        Ref(false),
        connections,
        Dict{String, Tuple{Bool, Observable}}(),
        Dict{Symbol, Any}[],
        Set{Asset}(),
        JSCode[]
    )
end

function Base.close(session::Session)
    foreach(close, session.connections)
    empty!(session.connections)
    empty!(session.observables)
    empty!(session.on_document_load)
    empty!(session.message_queue)
    empty!(session.dependencies)
end

function Base.copy(session::Session)
    obs = Dict((k => (true, map(identity, obs)) for (k, (regs, obs)) in session.observables))
    return Session(
        WebSocket[],
        session.observables,
        Dict{Symbol, Any}[],
        copy(session.dependencies),
        JSCode[]
    )
end

function Base.push!(session::Session, x::Observable)
    session.observables[x.id] = (false, x)
end

function Base.push!(session::Session, dependency::Dependency)
    for asset in dependency.assets
        push!(session, asset)
    end
    return dependency
end

function Base.push!(session::Session, asset::Asset)
    push!(session.dependencies, asset)
    if asset.onload !== nothing
        on_document_load(session, asset.onload)
    end
    return asset
end

function Base.push!(session::Session, websocket::WebSocket)
    push!(session.connections, websocket)
    filter!(isopen, session.connections)
    return session
end

"""
    send_queued(session::Session)

Sends all queued operations to the frontend
"""
function send_queued(session::Session)
    if isopen(session)
        # send all queued messages
        for message in session.message_queue
            send(session, message)
        end
        empty!(session.message_queue)
    else
        error("To send queued messages make sure that session is open.")
    end
end

"""
    queued_as_script(session::Session)

Returns all queued messages as a script that can be included into html
"""
function queued_as_script(io::IO, session::Session)
    # send all queued messages
    # first register observables
    for (id, (registered, observable)) in session.observables
        if !registered
            # Register on the JS side by sending the current value
            updater = JSUpdateObservable(session, id)
            # Make sure we update the Javascript values!
            on(updater, observable)
            session.observables[id] = (true, observable)
            serialize_string(io, js"    registered_observables[$(observable)] = $(observable[]);")
            println(io)
        end
    end
    for message in session.message_queue
        serialize_string(
            io,
            js"    process_message(deserialize_js($(AsJSON(message))));"
        )
        println(io)
    end
    empty!(session.message_queue)
end
queued_as_script(session::Session) = sprint(io-> queued_as_script(io, session))

"""
    send(session::Session; attributes...)

Send values to the frontend via JSON for now
"""
Sockets.send(session::Session; kw...) = send(session, Dict{Symbol, Any}(kw))

function Sockets.send(session::Session, message::Dict{Symbol, Any})
    if isopen(session) && !session.fusing[]
        # send all queued messages
        # send_queued(session)
        # sent the actual message
        for connection in session.connections
            serialize_websocket(connection, message)
        end
    else
        push!(session.message_queue, message)
    end
end

function fuse(f, session::Session)
    session.fusing[] = true
    f()
    session.fusing[] = false
    evaljs(session, JSCode([JSString(queued_as_script(session))]))
end

function Base.isopen(session::Session)
    return !isempty(session.connections) && isopen(session.connections[1])
end

"""
    onjs(session::Session, obs::Observable, func::JSCode)

Register a javascript function with `session`, that get's called when `obs` gets a new value.
If the observable gets updated from the JS side, the calling of `func` will be triggered
entirely in javascript, without any communication with the Julia `session`.
"""
function onjs(session::Session, obs::Observable, func::JSCode)
    # register the callback with the JS session
    register_resource!(session, (obs, func))

    send(
        session,
        type = OnjsCallback,
        id = obs.id,
        # eval requires functions to be wrapped in ()
        payload = js"($func)"
    )
end

"""
    onload(session::Session, node::Node, func::JSCode)

calls javascript `func` with node, once node has been displayed.
"""
function onload(session::Session, node::Node, func::JSCode)
    on_document_load(session, js"""
        // on document load, call func with the node
        ($(func))($node)
    """)
end

"""
    on_document_load(session::Session, js::JSCode)

executes javascript after document is loaded
"""
function on_document_load(session::Session, js::JSCode)
    register_resource!(session, js)
    push!(session.on_document_load, js)
end

"""
    linkjs(session::Session, a::Observable, b::Observable)

for an open session, link a and b on the javascript side. This will also
Link the observables in Julia, but only as long as the session is active.
"""
function linkjs(session::Session, a::Observable, b::Observable)
    # register the callback with the JS session
    onjs(
        session,
        a,
        js"""
        function (value){
            // update_obs will return false once b is gone,
            // so this will automatically deregister the link!
            return update_obs($b, value)
        }
        """
    )
end

"""
    evaljs(session::Session, jss::JSCode)

Evaluate a javascript script in `session`.
"""
function evaljs(session::Session, jss::JSCode)
    register_resource!(session, jss)
    send(session, type = EvalJavascript, payload = jss)
end

"""
    active_sessions(app::Application)

Returns all active sessions of an Application
"""
function active_sessions(app::Application)
    collect(filter(app.sessions) do (k, v)
        any(x-> isopen(x[2]), v) # leave not yet started connections
    end)
end

"""
    register_resource!(session::Session, domlike)

Walks dom like structures and registers all resources (Observables, Assets Depencies)
with the session.
"""
register_resource!(session::Session, @nospecialize(jss)) = nothing # do nothing for unknown type

function register_resource!(session::Session, list::Union{Tuple, AbstractVector, Pair})
    for elem in list
        register_resource!(session, elem)
    end
end

function register_resource!(session::Session, jss::JSCode)
    register_resource!(session, jss.source)
end

function register_resource!(session::Session, asset::Union{Asset, Dependency, Observable})
    push!(session, asset)
end

function register_resource!(session::Session, node::Node)
    walk_dom(session, node) do x
        register_resource!(session, x)
    end
end

function update_dom!(session::Session, dom)
    # empty!(session.on_document_load)
    dom = jsrender(session, dom)
    register_resource!(session, dom)
    innerhtml = repr(MIME"text/html"(), dom)
    new_deps = session.dependencies
    new_jss = JSCode(Any[])
    for jss in session.on_document_load
        append_source!(new_jss, jss)
    end
    register_obs!(session)
    script_urls = url.(new_deps)
    update_script = js"""
        var dom = document.getElementById('application-dom')
        dom.innerHTML = $(innerhtml)
        var urls = $(script_urls)
        for (var i = 0; i < urls.length; i++) {
            var s = document.createElement("script");
            s.type = "text/javascript";
            s.async = false
            s.src = urls[i];
            document.head.appendChild(s);
        }
        var s = document.createElement("script");
        s.type = "text/javascript";
        s.async = false
        s.text = $(serialize_string(new_jss));
        document.head.appendChild(s);
    """
    println(serialize_string(update_script))
    evaljs(session, update_script)
end


text:
# This is an example that shows you how to create VTK data types natively in Julia

using VTKDataTypes
import GeometryTypes

function create_image()
    origin = (0.0, 0.0, 0.0)
    spacing = (1.0, 2.0, 1.0)
    _extents = (2, 3, 5) #Must be Vector{Int}

    image = VTKImageData(origin, spacing, _extents)
    @test dim(image) == 3
    @test extents(image) == image.extents
    @test cell_extents(image) == image.extents .- 1

    image.point_data["Point scalar data"] = zeros(extents(image)...)
    image.cell_data["Cell scalar data"] = zeros(cell_extents(image)...)
    image.point_data["Point vector data"] = zeros(3, extents(image)...)
    image.cell_data["Cell vector data"] = zeros(3, cell_extents(image)...)

    @test cell_connectivity(image, (1,1,1)) == [1, 2, 3, 4, 7, 8, 9, 10]
        #Used to obtain the points which make up the cell in VTK convention
    @test cell_type(image, (1,1,1)) == 11
        #Used to obtain the VTK cell type

    @test image == deepcopy(image)
    return image
end

function create_rectilinear_data1()
    x = y = z = collect(range(-2, stop=2, length=5))

    rectilinear = VTKRectilinearData((x, y, z))
    @test dim(rectilinear) == 3
    @test extents(rectilinear) == (5,5,5)
    @test cell_extents(rectilinear) == (4,4,4)

    rectilinear.point_data["Point scalar data"] = zeros(extents(rectilinear)...)
    rectilinear.cell_data["Cell scalar data"] = zeros(cell_extents(rectilinear)...)
    rectilinear.point_data["Point vector data"] = zeros(3, extents(rectilinear)...)
    rectilinear.cell_data["Cell vector data"] = zeros(3, cell_extents(rectilinear)...)

    @test cell_connectivity(rectilinear, (1,1,1)) == [1, 2, 7, 6, 26, 27, 32, 31]
        #Used to obtain the points which make up the cell in VTK convention
    @test cell_type(rectilinear, (1,1,1)) == 12
        #Used to obtain the VTK cell type

    @test rectilinear == deepcopy(rectilinear)
    return rectilinear
end

function create_rectilinear_data2()
    image = create_image()
    rectilinear = VTKRectilinearData(image)

    @test dim(rectilinear) == dim(image) == 3
    @test extents(rectilinear) == extents(image) == (2,3,5)
    @test cell_extents(rectilinear) == cell_extents(image) == (1,2,4)

    return rectilinear
end

function create_structured_data()
    rectilinear = create_rectilinear_data1()
    structured = VTKStructuredData(rectilinear)
    @test size(structured.point_coords) == (3,extents(rectilinear)...) == (3,5,5,5)

    @test extents(structured) == extents(rectilinear) == (5,5,5)
    @test cell_extents(structured) == cell_extents(rectilinear) == (4,4,4)
    @test dim(structured) == dim(rectilinear) == 3

    structured.point_data["Point scalar data"] = zeros(extents(structured)...)
    structured.cell_data["Cell scalar data"] = zeros(cell_extents(structured)...)
    structured.point_data["Point vector data"] = zeros(3, extents(structured)...)
    structured.cell_data["Cell vector data"] = zeros(3, cell_extents(structured)...)

    @test cell_connectivity(structured, (1,1,1)) == [1, 2, 7, 6, 26, 27, 32, 31]
        #Used to obtain the points which make up the cell in VTK convention
    @test cell_type(structured, (1,1,1)) == 12
        #Used to obtain the VTK cell type

    @test structured == deepcopy(structured)
    return structured
end

function create_unstructured_data()
    rectilinear = create_rectilinear_data1()
    structured = create_structured_data()
    unstruct1 = VTKUnstructuredData(rectilinear)
    unstruct2 = VTKUnstructuredData(structured)

    @test length(unstruct1.cell_connectivity) == length(unstruct1.cell_types)
    @test size(unstruct1.point_coords) == (dim(unstruct1), num_of_points(unstruct1))
    @test size(unstruct1.point_data["Point scalar data"], 1) == size(unstruct1.point_coords, 2)
    @test size(unstruct1.cell_data["Cell vector data"]) == (3, num_of_cells(unstruct1))

    @test unstruct1 == unstruct2

    @test cell_connectivity(unstruct1, 1) == unstruct1.cell_connectivity[1]
    @test cell_type(unstruct1, 1) == unstruct1.cell_types[1]
    @test is_homogeneous(unstruct1)
    return unstruct1
end

function create_poly_data()
    rectilinear = create_rectilinear_data1()
    polydata = VTKPolyData(rectilinear)
        # automatically decomposes linear volume cells to faces and averages cell data

    @test length(polydata.cell_connectivity) == length(polydata.cell_types)
    @test size(polydata.point_coords) == (dim(polydata), num_of_points(polydata))
    @test size(polydata.point_data["Point scalar data"], 1) == size(polydata.point_coords, 2)
    @test size(polydata.cell_data["Cell vector data"]) == (3, num_of_cells(polydata))

    @test polydata == polydata
    return polydata
end

function create_multiblock_data()
    mb = VTKMultiblockData((create_image(), create_rectilinear_data1(),
        create_structured_data(), create_unstructured_data(), create_poly_data()))
    @test length(mb) == 5
    for b in mb
        @test isa(b, AbstractVTKSimpleData)
    end
    mb2 = VTKMultiblockData((mb, create_image()))
    for b in simple_block_generator(mb2) #Recursively unwraps any VTKMultiblockData
        @test isa(b, AbstractVTKSimpleData)
    end

    return mb2[1]
end

function create_timeseries_data()
    timemarkers = [0, 0.5, 1, 1.5]

    #2D polydata
    x = y = [0., 1., 2.];
    point_coords = (x, y);
    point_data = Dict{String, Array{Float64}}();
    cell_data = Dict{String, Array{Float64}}();
    a = VTKRectilinearData(point_coords, point_data, cell_data);
    a.cell_data["Color"] = reshape(rand(num_of_cells(a)), cell_extents(a));
    b = VTKPolyData(a);

    x = y = [3., 4., 5.];
    point_coords = (x, y);
    c = VTKPolyData(VTKRectilinearData(point_coords, point_data, cell_data));
    c.cell_data["Color"] = rand(num_of_cells(c));

    m = VTKMultiblockData((b, c));
    timeseries = VTKTimeSeriesData([0.], [m]);

    no_of_timesteps = 5;
    timestep = 0.5;
    speed = 1;
    i = 0
    while i < no_of_timesteps
        i += 1;
        new_data = deepcopy(timeseries.data[i]);

        #Random walk of first dataset
        new_data[1].point_coords = new_data[1].point_coords +
            vcat(fill(speed*rand()-speed/2, (1,num_of_points(new_data[1]))),
            fill(speed*rand()-speed/2, (1,num_of_points(new_data[1]))));

        #Random walk of second dataset
        new_data[2].point_coords = new_data[2].point_coords +
            vcat(fill(speed*rand()-speed/2, (1,num_of_points(new_data[2]))),
            fill(speed*rand()-speed/2, (1,num_of_points(new_data[2]))));

        insert_timed_data!(timeseries, i*timestep, new_data);
    end

    for ts in timeseries
        @test isa(ts, VTKMultiblockData)
    end

    @test isa(timeseries[1.1], VTKMultiblockData) #Indexing by time, does linear interpolation
    @test isa(timeseries[1], VTKMultiblockData) #Indexing by frame index

    return timeseries
end

function triangulation_and_decompose()
    rectilinear = create_rectilinear_data1()
    tridata = triangulate(rectilinear)
    @test triangular(tridata)
    @test isa(GLMesh(rectilinear), GeometryTypes.HomogenousMesh)
    @test isa(decompose(rectilinear, "Faces"), VTKPolyData)
    @test triangulate(decompose(rectilinear, "Faces")) == tridata

    return tridata
end

create_image()
create_rectilinear_data1()
create_rectilinear_data2()
create_structured_data()
create_unstructured_data()
create_poly_data()
create_multiblock_data()
create_timeseries_data()
triangulation_and_decompose()

ext: jl


text:
# Arvind Prasadan
# 2019 January
# DMF Package
# Finds the squared error between two sets of eigenvectors, accounting for permutations

"""
	eigenvector_error(Q, Q_hat)

Finds the squared error
``\\|q_i - p_i \\widehat{q}_i\\|_2^2 = q_i^T q_i + \\widehat{q}_i^T \\widehat{q}_i - 2 p_i q_i^T\\widehat{q}_i = 2 [1 - |q_i^T \\widehat{q}_i|]`` for each q_i

# Arguments
- `Q`: The true vectors (columns)
- `Q_hat`: The estimated vectors (columns)

# Output
- Error
"""
function eigenvector_error(Q, Q_hat)
	@assert(size(Q, 1) == size(Q_hat, 1), "Q and Q_hat have different dimensions")
	@assert(size(Q, 2) == size(Q_hat, 2), "Q and Q_hat have different dimensions")

	return 2.0 * sum(1.0 .- abs.(permute_IP(Q, Q_hat)))
end

# Old function
"""
function eigenvector_error(Q, Q_hat)

    if size(Q, 1) != size(Q_hat, 1)
       return Inf
    end

    Q = mapslices(normalize, Q; dims = 1)
    Q_hat = mapslices(normalize, Q_hat; dims = 1)
    k = size(Q, 2)

    return abs(2.0 * (k - norm(mapslices(maximum, abs.(Q_hat' * Q); dims = 1), 2)^2.0))
end
"""


text:
using RecipesBase: apply_recipe

@testset "BLSPeriodogram plotting" begin
    t, y, dy, params = make_data()
    result = BLS(t, y, dy; params.duration)
    recipes = apply_recipe(Dict{Symbol,Any}(), result)
    for rec in recipes
        @test getfield(rec, 1) == Dict{Symbol,Any}(
            :seriestype => :line,
            :xguide => "period",
            :yguide => "log-likeilhood")

            # test to make sure our position is correct (should be +0.5 the given (x,y))
        period = rec.args[1]
        power = rec.args[2]

        @test period == result.periods
        @test power == result.power
    end
end


text:
plot_panel(panel :: Panel3D) = let ps = panel_coordinates(panel); Tuple.([ ps; [ps[1]] ]) end

"""
    plot_panels(panels :: Array{Panel3D})

Get the coordinates of an array of `Panel3D` for plotting.
"""
plot_panels(panels) = plot_panel.(vec(panels))

# foil_coords = [ [ [coord[1]; 0; coord[2]] .* chord .+ loc for coord in foil.coordinates ] for (chord, foil, loc) in zip(wing.right.chords[end:-1:1], wing.right.foils[end:-1:1], wing_coords) ]

function plot_planform(mesh :: Matrix{SVector{3,T}}) where T <: Real
    coords  =   [
                    mesh[1,1:end-1];
                    mesh[1:end-1,end];
                    mesh[end,end:-1:2];
                    mesh[end:-1:1,1]
                ]

    permutedims(combinedimsview(coords))
end

"""
    plot_planform(wing :: AbstractWing)

Get the planform coordinates of an `AbstractWing` for plotting.
"""
plot_planform(wing :: AbstractWing) = plot_planform(coordinates(wing))

plot_streamlines(system :: VortexLatticeSystem, points, length, num_steps) = Tuple.(streamlines(system, points, length, num_steps))

"""
    plot_surface(wing :: AbstractWing)
    plot_surface(wing :: AbstractWing,
                 span_num :: Vector{Integer},
                 chord_num :: Integer)

Get the surface coordinates of an `AbstractWing` for plotting, optionally with specified spanwise and chordwise discretizations.
"""
plot_surface(wing :: AbstractWing) = plot_panels(mesh_wing(wing, [length(spans(wing))], minimum(length, foil.x for foil in foils(wing))))

plot_surface(wing :: AbstractWing, span_num, chord_num) = plot_panels(mesh_wing(wing, span_num, chord_num))


text:
# This file inludes the test set for coding

@testset "Coding" begin
    # Fields
    @test hasfield(GrayCoding, :pairs)

    # Construction
    @test_throws Exception GrayCoding(3)
    @test_throws Exception GrayCoding(7)
    coding = GrayCoding(4)
    coding = GrayCoding(8)

    # Check bit to symbol mapping
    coding = GrayCoding(4)
    @test coding.pairs[[0, 0]] == 1
    @test coding.pairs[[0, 1]] == 2
    @test coding.pairs[[1, 1]] == 3
    @test coding.pairs[[1, 0]] == 4

    # Check bit stream mapping
    bits = [0, 1, 1, 1, 1, 1, 0, 0, 1, 0]
    symbols = coding(bits)
    @test symbols == [2, 3, 3, 1, 4]
end


text:
############
# BuildRef #
############

mutable struct BuildRef
    repo::String  # the build repo
    sha::String   # the build + status SHA
    vinfo::String # versioninfo() taken during the build
end

BuildRef(repo, sha) = BuildRef(repo, sha, "retrieving versioninfo() failed")

function Base.:(==)(a::BuildRef, b::BuildRef)
    return (a.repo == b.repo &&
            a.sha == b.sha &&
            a.vinfo == b.vinfo)
end

Base.summary(build::BuildRef) = string(build.repo, SHA_SEPARATOR, snipsha(build.sha))

# if a PR number is included, attempt to build from the PR's merge commit
# FIXME: re-use PkgEval's BinaryBuilder-based build
function build_julia!(config::Config, build::BuildRef, logpath, prnumber::Union{Int,Nothing}=nothing)
    # make a temporary workdir for our build
    builddir = mktempdir(workdir(config))
    cd(workdir(config))

    # clone/fetch the appropriate Julia version
    if prnumber !== nothing
        # clone from `trackrepo`, not `build.repo`, since that's where the merge commit is
        gitclone!(config.trackrepo, builddir)
        cd(builddir)
        try
            run(`git fetch --quiet origin +refs/pull/$(prnumber)/merge:`)
        catch
            # if there's not a merge commit on the remote (likely due to
            # merge conflicts) then fetch the head commit instead.
            run(`git fetch --quiet origin +refs/pull/$(prnumber)/head:`)
        end
        run(`git checkout --quiet --force FETCH_HEAD`)
        build.sha = readchomp(`git rev-parse HEAD`)
    else
        gitclone!(build.repo, builddir)
        cd(builddir)
        run(`git checkout --quiet $(build.sha)`)
    end

    # set up logs for STDOUT and STDERR
    logname = string(build.sha, "_build")
    outfile = joinpath(logpath, string(logname, ".out"))
    errfile = joinpath(logpath, string(logname, ".err"))

    # run the build
    cpus = mycpus(config)
    run(pipeline(`make -j$(length(cpus)) USECCACHE=1`,
                 stdout=outfile, stderr=errfile))
    cd(workdir(config))
    return builddir
end


text:
# tmux attach
# source activate julia
# julia -p 3 experiments/exp_surface_extraction_run.jl
# Ctrl-b d

# paper note
# the deinit with nothing, automatic end of while true
# using Revise
using Test
using Logging
using Distributed
@info "distributed loaded"
if nprocs() == 1
    addprocs(3)
end

using Dates
using ExSup

@everywhere using LarSurf

@info "LarSurf everywhere"
# global_logger(SimpleLogger(stdout, Logging.Debug))
# # set logger on all workers
# for wid in workers()
#     @spawnat wid global_logger(SimpleLogger(stdout, Logging.Debug))
# end

fn = "exp_surface_extraction6_test.csv"

# Number of logical CPU cores available in the system.
# println("CPU cores: ", Sys.CPU_CORES)
Sys.cpu_summary()
# Number of available processes.
number_procs = nprocs()
number_workers = nworkers()
println("nprocs: ", nprocs())
println("nworkers: ", nworkers())

include("prepare_data.jl")

general_prepare_data = LarSurf.generate_almost_cube
general_prepare_data = LarSurf.generate_cube

segmentation = LarSurf.generate_cube(10)

fcns = [
    # TODO prepare
    # LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_parallel,
    # LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_fixed_block_size,
    # LarSurf.get_surface_grid_per_block_Vreduced_FVreduced_old,
    # LarSurf.get_surface_grid_per_block_Vreduced_FVreduced,
    # LarSurf.get_surface_grid_per_block_FVreduced,
    # LarSurf.get_surface_grid_per_block_full,
    # LarSurf.get_surface_grid_old,
    # LarSurf.get_surface_grid,
    LarSurf.lsp_get_surface,
    # LarSurf.get_surface_grid_per_voxel,
]
# set last two are with one parameter
nargs = 1 * ones(Int64, length(fcns))
# nargs[end-1:end] .= 1
# inargs = [segmentation, block_size]
fcns_nargs = collect(zip(fcns, nargs))

# it is list with just fast functions
# fcns_fast = fcns_nargs[1:end-4]
fcns_fast = fcns_nargs
fcns_all  = fcns_nargs

@info "Nuber of functions per experiment: $(length(fcns_all))"

# = 17.6.
# # set last two are with one parameter
# nargs = 2 * ones(Int64, length(fcns))
# nargs[end-1:end] .= 1
# # inargs = [segmentation, block_size]
# fcns_nargs = collect(zip(fcns, nargs))
#
# # it is list with just fast functions
# fcns_fast = fcns_nargs[1:end-4]
# fcns_all  = fcns_nargs

# fcns_nargs_local = fcns_nargs[1:end]
# end
# fcns_nargs = [
#     (fcns[i], nargs[])
#     # (fcns[i], inargs[1:nargs[i]])
#     for i=1:length(fcns)
# ]

# Warm up on small data
@info "First Warm Up..."
block_size = [2,2,2]
LarSurf.lsp_setup(block_size)
for (fcni, nargs) in fcns_nargs
    @debug "running on small segmentation with size $(size(segmentation))"
    argsi = [segmentation, block_size]
    fcni(argsi[1:nargs]...)
end
@info "...done"

function save_data(experiment, timed, segmentation, b3_size, append_dct)
    # @info "save data "
    tm = timed[2]
    al = timed[3]
    println(" time=$tm, alloc=$al ", append_dct["fcn"])
    dct = Dict()
    # dct = ExSup.timed_and_more_to_dict!(dct, timed;experiment=experiment)
    dct = ExSup.timed_to_dict!(dct, timed)
    dct = ExSup.datetime_to_dict!(dct)
    dct["experiment"] = experiment
    dct = ExSup.segmentation_description_to_dict!(dct, segmentation)
    dct = ExSup.size_to_dict!(dct, b3_size, "boundary3_")
    # @info "save data done before append" append_dct dct
    if append_dct != nothing
        merge!(dct, append_dct)
    end
    # @info "save data done append done"

    # @info "save data done append done" dct
    # @info "save data before save" keys(dct) dct
    ExSup.add_to_csv(dct, fn)
    # @info "save data done append done"
end

"""
fcns_nargs_local
"""
function run_measurement(
    fcns_nargs_local, prepare_data_parameter, block_size,
    experiment=nothing; append_dct=nothing, skip_slow=false, data_fcn=nothing)

    if experiment == nothing
        experiment = "time measurement"
    end

    # LasRurf.set_block_size(block_size) # this is done by fallowing function too
    LarSurf.lsp_setup(block_size)
    # TODO set time
    LarSurf.set_time_data(Dict())

    println(experiment)
    if data_fcn == nothing
        prepare_data = general_prepare_data
    else
        prepare_data = data_fcn
    end

    segmentation = prepare_data(prepare_data_parameter)
    for (fcni, nargs) in fcns_nargs_local
        @info "==== Running $(String(Symbol(fcni)))"

        argsi = [segmentation, block_size]
        tmd = @timed(fcni(argsi[1:nargs]...))
        if append_dct == nothing
            append_dct = Dict()
        end
        append_dct["fcn"] = String(Symbol(fcni))
        append_dct["nprocs"] = number_procs
        append_dct["nworkers"] = number_workers
        append_dct["jlfile"] = @__FILE__
        append_dct["hostname"] = gethostname()
        append_dct["ncores"] = length(Sys.cpu_info())
        append_dct["data parameter"] = prepare_data_parameter
        tme_data = LarSurf.get_time_data()
        # @info "merging append dct" append_dct tme_data
        merge!(append_dct, tme_data)

        # @info "time data done"
        save_data(experiment, tmd, segmentation, block_size, append_dct)
    end
    println("============================")
end

# Warming
@info "Warming..."
run_measurement(fcns_all , 10, [1,1,1] .*  8, "warming")
run_measurement(fcns_fast, 20, [1,1,1] .* 16, "warming"; skip_slow=true)
run_measurement(fcns_fast, 40, [1,1,1] .* 32, "warming"; skip_slow=true)
run_measurement(fcns_fast, 70, [1,1,1] .* 64, "warming"; skip_slow=true)
@info "...done"

for i=1:0
    @info "small scale experiments"
    ## Small
    block_size = [1,1,1] .* 16
    # LarSurf.lsp_setup(block_size)
    run_measurement(fcns_fast,  40, block_size, "small b3")
    run_measurement(fcns_fast,  40, block_size, "small b3")
    run_measurement(fcns_fast,  40, block_size, "small b3")

    ## Datasize, constant b3

    block_size = [1,1,1] .* 64
    # LarSurf.lsp_setup(block_size)
    run_measurement(fcns_fast,  20, block_size, "data size")
    run_measurement(fcns_fast,  40, block_size, "data size")
    run_measurement(fcns_fast,  80, block_size, "data size")
    run_measurement(fcns_fast, 160, block_size, "data size")
#     block_size = [1,1,1] .* 32
#     # LarSurf.lsp_setup(block_size)
#     run_measurement(fcns_fast, 512, [1,1,1] .*  8, "boundary size big 32")
#     run_measurement(fcns_fast, 512, [1,1,1] .* 16, "boundary size big 32")
#     run_measurement(fcns_fast, 512, [1,1,1] .* 32, "boundary size big 32")
#     run_measurement(fcns_fast, 512, [1,1,1] .* 64, "boundary size big 32")
end

# Experiments
for i=1:0
    @info "big data experiments"
    block_size = [1,1,1] .* 64
    run_measurement(fcns_fast, 320, block_size, "data size")
    run_measurement(fcns_fast, 512, block_size, "data size")
end

for i=1:0
    ## Boundary matrix size
    @info "boundary matrix size experiments"

    run_measurement(fcns_fast, 100, [1,1,1] .*  8, "boundary size")
    run_measurement(fcns_fast, 100, [1,1,1] .* 16, "boundary size")
    run_measurement(fcns_fast, 100, [1,1,1] .* 32, "boundary size")
    run_measurement(fcns_fast, 100, [1,1,1] .* 64, "boundary size")
    # run_measurement(fcns_fast,100, [1,1,1] .* 64, "data size"; skip_slow=true)
end

for i=1:1
    @info "slow experiments"
    # block_size = [1,1,1] .* 32
    # LarSurf.lsp_setup(block_size)
    # run_measurement(fcns_fast, 512, [1,1,1] .*  8, "boundary size big 32")
    # run_measurement(fcns_fast, 512, [1,1,1] .* 16, "boundary size big 32")
    # run_measurement(fcns_fast, 512, [1,1,1] .* 32, "boundary size big 32")
    # run_measurement(fcns_fast, 512, [1,1,1] .* 64, "boundary size big 32")

    run_measurement(fcns_fast, 512, [1,1,1] .* 64, "boundary size big")
    run_measurement(fcns_fast, 512, [1,1,1] .* 32, "boundary size big")
    run_measurement(fcns_fast, 512, [1,1,1] .* 16, "boundary size big")
    run_measurement(fcns_fast, 512, [1,1,1] .*  8, "boundary size big")
    run_measurement(fcns_fast, 512, [1,1,1] .* 128, "boundary size big")
end

# CT data
for i=1:0
    for id=1:0
        @info "CT data"
        # block_size = [1,1,1] .* 32
        # LarSurf.lsp_setup(block_size)
        # run_measurement(fcns_fast, 512, [1,1,1] .*  8, "boundary size big 32")
        # run_measurement(fcns_fast, 512, [1,1,1] .* 16, "boundary size big 32")
        # run_measurement(fcns_fast, 512, [1,1,1] .* 32, "boundary size big 32")
        # run_measurement(fcns_fast, 512, [1,1,1] .* 64, "boundary size big 32")

        run_measurement(fcns_fast, id, [1,1,1] .*  64, "Ircadb1"; data_fcn=prepare_ircad)
    end
end


text:
function RiccatiEquationSolverImpl(arg0::RealMatrix, arg1::RealMatrix, arg2::RealMatrix, arg3::RealMatrix)
    return RiccatiEquationSolverImpl((RealMatrix, RealMatrix, RealMatrix, RealMatrix), arg0, arg1, arg2, arg3)
end

function get_k(obj::RiccatiEquationSolverImpl)
    return jcall(obj, "getK", RealMatrix, ())
end

function get_p(obj::RiccatiEquationSolverImpl)
    return jcall(obj, "getP", RealMatrix, ())
end


text:
module PoleExploration

using AbstractPlotting.MakieLayout
using Makie
using ControlSystems
using UnicodeFun

export start

include("roots.jl")
include("plotting.jl")

function scenesetup()
    # Variables
    roots = Node(Root[Root(Point2f0(-1, 0), true, false)])

    zeros = lift(get_zeros, roots)
    poles = lift(get_poles, roots)
    gain = Node(1.0)

    selected_data = lift(get_selected, roots)
    selected_idx = lift(x -> x[1], selected_data)
    selected_pos = lift(x -> x[2], selected_data)
    selected_token = lift(x -> x[3], selected_data)

    sys = lift((z, p, k) -> begin
        zpk([a + b*im for (a, b) in z], [a + b*im for (a, b) in p], k)
    end, zeros, poles, gain)

    # Parent scene
    outer_padding = 30
    scene, layout = layoutscene(outer_padding, resolution = (1200, 1000), backgroundcolor = RGBf0(0.98, 0.98, 0.98))

    # General settings
    kwargs = Dict(:xtickalign => 1, :ytickalign => 1, :yticklabelspace => 25)

    # Root locus
    root_ax = layout[1:2, 1] = LAxis(scene; title = "Root locus", kwargs...)
    scatter!(root_ax, poles, color=:red, marker='+', markersize=10)
    scatter!(root_ax, zeros, color=:red, marker='o', markersize=10)
    scatter!(root_ax, selected_pos, color=:green, marker=selected_token, markersize=10)

    # Step plot
    step_ax = layout[3, 1] = LAxis(scene; title = "Step", kwargs...)
    step_points = lift(sys -> begin
        y, t, x = step(sys)
        limits!(step_ax, find_limits(t), find_limits(y))
        convert.(Point2f0, zip(t, vec(y)))
    end, sys)
    lines!(step_ax, step_points)

    # Impulse plot
    impulse_ax = layout[4, 1] = LAxis(scene; title = "Impulse", kwargs...)
    impulse_points = lift(sys -> begin
        y, t, x = impulse(sys)
        limits!(impulse_ax, find_limits(t), find_limits(y))
        #axes[4].xticks = 0:pi:2pi
        convert.(Point2f0, zip(t, vec(y)))
    end, sys)
    lines!(impulse_ax, impulse_points)

    # Bode plot
    bodemag_ax = layout[1, 2] = LAxis(scene; title = "Magnitude", kwargs...)
    bodephase_ax = layout[2, 2] = LAxis(scene; title = "Phase", kwargs...)
    linkxaxes!(bodemag_ax, bodephase_ax)
    bodevars = lift(sys -> begin
        mag, phase, w = bode(sys)
        logmag = log10.(mag)
        logw = log10.(w)
        wlims = find_limits(logw)
        maglims = find_limits(logmag)
        phaselims = find_limits(phase)
        wticks = log_ticks(wlims)
        magticks = log_ticks(maglims)
        bodemag_ax.xticks = wticks
        bodemag_ax.yticks = magticks
        bodephase_ax.xticks = wticks
        limits!(bodemag_ax, wlims, maglims)
        limits!(bodephase_ax, wlims, phaselims)
        logw, logmag, phase
    end, sys) # TODO use nyquist instead to calculate this?
    bodemag_points = lift(x -> convert.(Point2f0, zip(x[1], x[2])), bodevars)
    bodephase_points = lift(x -> convert.(Point2f0, zip(x[1], x[3])), bodevars)
    lines!(bodemag_ax, bodemag_points)
    lines!(bodephase_ax, bodephase_points)

    # Nyquist plot
    nyquist_ax = layout[3:4, 2] = LAxis(scene; title = "Nyquist", kwargs...)
    nyquist_points = lift(sys -> begin
        a, b, _ = nyquistv(sys)
        limits!(nyquist_ax, find_limits(a, start=[-1, 1]), find_limits(b, start=[-1, 1]))
        convert.(Point2f0, zip(a, b))
    end, sys)
    lines!(nyquist_ax, nyquist_points)
    lines!(nyquist_ax, cos.(0:0.01:2pi), sin.(0:0.01:2pi), color=:red)

    gain_slider = LSlider(scene, range=-10:0.01:10, startvalue=gain[])
    gain_label = LText(scene, lift(x -> "K=$(x)", gain_slider.value))
    on(gain_slider.value) do value
        gain[] = value
    end
    layout[0, 1:2] = hbox!(gain_slider, gain_label)

    # Other
    tf_text = lift(print_tf, roots, gain)
    #tf_label = layout[0, 2] = MakieTeX.LTeX(scene, raw"\int \mathbf E \cdot d\mathbf a = \frac{Q_{encl}}{4\pi\epsilon_0}", tellwidth=false)
    tf_label = layout[0, 1:2] = LText(scene, text=tf_text, tellwidth=false, font="Courier")

    mousestate = addmousestate!(root_ax.scene)
    onmouseleftclick(mousestate) do state
        # Find closest point, if point is within reasonable distance given scale select it
        root = find_close(state.pos, roots[], root_ax.limits[].widths ./ 50)
        # TODO set selected
        unselect_all!(roots) # Unselects without sending out update
        if !isnothing(root)
            root.selected = true
            roots[] = roots[] # Update
            #@show "selected", root
        end
    end
    onmouseleftdoubleclick(mousestate) do state
        x = state.pos[1]
        y = state.pos[2]
        if abs(y) < root_ax.limits[].widths[2] / 100
            y = 0
        end
        unselect_all!(roots) # Unselects without sending out update
        newroot = Root(Point2f0(x, abs(y)), true, true) # Add new root and send update
        roots[] = push!(roots[], newroot)
    end
    onmouserightclick(mousestate) do state
        if selected_idx[] != 0
            x = state.pos[1]
            y = state.pos[2]
            if roots.val[selected_idx[]].pos[2] == 0
                y = 0
            else
                y = abs(y)
            end
            roots.val[selected_idx[]].pos = Point2f0(x, y)
            roots[] = roots[]
        end
    end

    on(scene.events.keyboardbuttons) do button
        if ispressed(button, Keyboard.r)
            roots[] = Root[Root(Point2f0(-1, 0), true, false)]
            set_close_to!(gain_slider, 1.0)
            autolimits!(root_ax)
        elseif ispressed(button, Keyboard.space)
            if selected_idx[] != 0
                if roots.val[selected_idx[]].pole && length(poles[]) - length(zeros[]) < 2 * (1 + (roots.val[selected_idx[]].pos[2] != 0))
                    println("Can't change that since it would create an system with more poles than zeros.")
                else
                    roots.val[selected_idx[]].pole = !roots.val[selected_idx[]].pole
                    roots[] = roots[]
                end
            end
        elseif ispressed(button, Keyboard.delete)
            if selected_idx[] != 0
                if roots.val[selected_idx[]].pole && length(poles[]) - length(zeros[]) < (1 + (roots.val[selected_idx[]].pos[2] != 0))
                    println("Can't remove that since it would create an system with more poles than zeros.")
                else
                    roots[] = [roots[][1:selected_idx[]-1]; roots[][selected_idx[]+1:end]]
                    selected_idx[] = 0
                end
            end
        elseif ispressed(button, Keyboard.z)
            sys[] = sys[]
            autolimits!(root_ax)
        end
    end

    scene
end

function start()
    println("""
    This is a tool for pole/zero exploration.

    Double left click will add poles in the root locus.
    Select roots with left click to modify them.
    * Space will switch between pole and zero for the selected root.
    * Right click moves the selected root.
    * Delete will remove the selected root.

    Panning is done with left click. Zooming can be done with scroll or right click drag.
    Zooming or panning while holding x or y will constrain the zoom/pan to that axis.
    Pressing z will reset all zoom/pan to the automatic value.

    Pressing r will reset everything to the start configuration.
    """)
    scenesetup()
end

end

max_stars_repo_name: albheim/PoleExploration

max_issues_repo_name: albheim/PoleExploration

max_forks_repo_name: albheim/PoleExploration


text:
# SPDX-License-Identifier: X11
# 2020-09-23
# : (Colon) (300pt)

lawofcosines(a, b, θ) = sqrt(a * a + b * b - 2a * b * cos(θ))
angleofhour(h,m) = π * (60h + m) / (6 * 60)
angleofminute(m) = 2π * m / 60

A, B, H, M = parse.(Int, split(readline()))

theta = abs(angleofminute(M) - angleofhour(H,M))
println(lawofcosines(A, B, theta))

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using UnivariateFunctions
using Dates
using Random
using GLM

tol = 10*eps()

Random.seed!(1)
obs = 1000
X = rand(obs)
y = X .+ rand(Normal(),obs) .+ 7

# Basic use case with 2 degrees
lm1 = fit(LinearModel, hcat(ones(obs), X, X .^ 2), y)
glm_preds = predict(lm1, hcat(ones(obs), X, X .^ 2))
package_approximation = UnivariateFunctions.create_ols_approximation(y, X, 0.0, 2, true)
package_predictions = UnivariateFunctions.evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

# Degree of 1 with no intercept
lm1 = fit(LinearModel, hcat(X), y)
glm_preds = predict(lm1, hcat(X))
package_approximation = UnivariateFunctions.create_ols_approximation(y, X, 0.0, 1, false)
package_predictions = UnivariateFunctions.evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

# Degree of 1 with no intercept
lm1 = fit(LinearModel, hcat(ones(obs)), y)
glm_preds = predict(lm1, hcat(ones(obs)))
package_approximation = UnivariateFunctions.create_ols_approximation(y, X, 0.0, 0, true)
package_predictions = UnivariateFunctions.evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

# With applying a base
Xmin = X .- 5
lm1 = fit(LinearModel, hcat(ones(obs), Xmin, Xmin .^ 2), y)
glm_preds = predict(lm1, hcat(ones(obs), Xmin, Xmin .^ 2))
package_approximation = UnivariateFunctions.create_ols_approximation(y, X, 5.0, 2, true)
package_predictions = UnivariateFunctions.evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

# With x as dates and a base.
X = Array{Date}(undef, obs)
baseDate =  Date(2016, 7, 21)
StartDate = Date(2018, 7, 21)
for i in 1:obs
    X[i] = StartDate +Dates.Day(2* (i-1))
end
XConverted = years_between.(X, baseDate)
lm1 = fit(LinearModel, hcat(ones(obs), XConverted, XConverted .^ 2), y)
glm_preds = predict(lm1, hcat(ones(obs), XConverted, XConverted .^ 2))
package_approximation = UnivariateFunctions.create_ols_approximation(y, X, baseDate, 2, true)
package_predictions = UnivariateFunctions.evaluate.(package_approximation, X)
sum(abs.(glm_preds .- package_predictions)) < 1e-10

# Chebyshev approximation
func = sin
nodes  =  12
degree =  8
left   = -2.0
right  =  5.0
approxim = UnivariateFunctions.create_chebyshev_approximation(func, nodes, degree, left, right)
X = convert(Array{Float64,1}, left:0.01:right)
y = func.(X)
y_approx = evaluate.(Ref(approxim), X)
maximum(abs.(y .- y_approx)) < 0.01

func = exp
nodes  =  12
degree =  8
left   =  1.0
right  =  5.0
approxim = UnivariateFunctions.create_chebyshev_approximation(func, nodes, degree, left, right)
X = convert(Array{Float64,1}, left:0.01:right)
y = func.(X)
y_approx = evaluate.(Ref(approxim), X)
maximum(abs.(y .- y_approx)) < 0.01


text:
export coordinates, PhysicalGrid, limits, origin, cellsize

"""
    coordinates(w::GridData;[dx=1.0],[I0=(1,1)])

Return a tuple of the ranges of the physical coordinates in each direction for grid
data `w`. If `w` is of `Nodes` type, then it returns a tuple of the form
`xg,yg`. If `w` is of `Edges` or `NodePair` type, then it returns a tuple of
the form `xgu,ygu,xgv,ygv`.

The optional keyword argument `dx` sets the grid spacing; its default is `1.0`. The
optional keyword `I0` accepts a tuple of integers to set the index pair of the
primal nodes that coincide with the origin. The default is `(1,1)`.

# Example

```jldoctest
julia> w = Nodes(Dual,(12,22));

julia> xg, yg = coordinates(w,dx=0.1)
(-0.05:0.1:1.05, -0.05:0.1:2.0500000000000003)
```
"""
function coordinates end

for (gridtype,ctype,dnx,dny,shiftx,shifty) in @generate_scalarlist(SCALARLIST)
   @eval coordinates(w::$gridtype{$ctype,NX,NY,T};dx::Float64=1.0,I0::Tuple{Int,Int}=(1,1)) where {NX,NY,T} =
    dx.*((1-I0[1]-$shiftx):(NX-$dnx-I0[1]-$shiftx),
         (1-I0[2]-$shifty):(NY-$dny-I0[2]-$shifty))

end

for (gridtype,ctype,dunx,duny,dvnx,dvny,shiftux,shiftuy,shiftvx,shiftvy) in @generate_collectionlist(VECTORLIST)
   @eval coordinates(w::$gridtype{$(ctype...),NX,NY,T};dx::Float64=1.0,I0::Tuple{Int,Int}=(1,1)) where {NX,NY,T,DDT} =
    dx.*((1-I0[1]-$shiftux):(NX-$dunx-I0[1]-$shiftux),
         (1-I0[2]-$shiftuy):(NY-$duny-I0[2]-$shiftuy),
         (1-I0[1]-$shiftvx):(NX-$dvnx-I0[1]-$shiftvx),
         (1-I0[2]-$shiftvy):(NY-$dvny-I0[2]-$shiftvy))

end

struct PhysicalGrid{ND}
  N :: NTuple{ND,Int}
  I0 :: NTuple{ND,Int}
  Δx :: Float64
  xlim :: NTuple{ND,Tuple{Real,Real}}
end

"""
    PhysicalGrid(xlim::Tuple{Real,Real},ylim::Tuple{Real,Real},Δx::Float64)

Constructor to set up a grid connected to physical space. The region to be
discretized by the grid is defined by the limits `xlim` and `ylim`, and the
cell spacing (uniform and indentical in each direction) is specified by `Δx`.
The constructor uses this information to determine the number of
cells in each direction, expanding the given range if necessary to accommodate
an integer number. It also pads each side with a ghost cell.
It also determines the indices corresponding to the corner
of the cell to which the physical origin corresponds. Note that the corner
corresponding to the lowest limit in each direction has indices (1,1).
"""
function PhysicalGrid(xlim::Tuple{Real,Real},
                      ylim::Tuple{Real,Real},Δx::Float64)

  #= set grid spacing and the grid position of the origin
  In case the physical limits are not consistent with an integer number of dual cells, based on
  the given Δx, we adjust them outward a bit in all directions. We also seek to place the
  origin on the corner of a cell.
  =#
  xmin, xmax = xlim
  ymin, ymax = ylim
  @assert xmax >= xmin && ymax >= ymin "Maximum limits must exceed minimum limits"
  Lx = xmax-xmin
  Ly = ymax-ymin

  NX, i0, xlimnew = _find_efficient_1d_grid(xmin,xmax,Δx)
  NY, j0, ylimnew = _find_efficient_1d_grid(ymin,ymax,Δx)

  PhysicalGrid((NX,NY),(i0,j0),Δx,(xlimnew,ylimnew))
end

function _set_1d_grid(xmin::Real,xmax::Real,Δx::Float64)
  NL, NR = floor(Int,xmin/Δx), ceil(Int,xmax/Δx)
  return NR-NL+2, 1-NL, (Δx*NL, Δx*NR)
end

function _factor_prime(N::Integer)
    # Determine if N can be factorized into small primes
    # Returns the list of powers and the remaining factor
    # after factorization.
    Nr = N
    blist = [2,3,5,7,11,13]
    pow = zero(blist)
    for (i,b) in enumerate(blist)
        while mod(Nr,b) == 0
            pow[i] += 1
            Nr = Int(Nr/b)
        end
    end
    return pow, Nr
end

function _find_efficient_1d_grid(xmin::Real,xmax::Real,Δx::Float64)
    # Based on the provided grid dimensions and grid spacing,
    # find the minimally larger grid that has a number of cells
    # that can be factorized into small primes for efficient FFT calculations.
    N, i0, xlimnew = _set_1d_grid(xmin,xmax,Δx)
    pow, Nr = _factor_prime(N)
    while Nr > 1 || sum(pow[end-1:end]) > 1
        xminnew, xmaxnew = xlimnew
        xminnew -= Δx
        xmaxnew += Δx
        N, i0, xlimnew = _set_1d_grid(xminnew,xmaxnew,Δx)
        pow, Nr = _factor_prime(N)
    end
    return N, i0, xlimnew
end

"""
    size(g::PhysicalGrid,d::Int) -> Int

Return the number of cells in direction `d` in grid `g`.
"""
Base.size(g::PhysicalGrid,d::Int) = g.N[d]

"""
    size(g::PhysicalGrid) -> Tuple

Return a tuple of the number of cells in all directions in grid `g`.
"""
Base.size(g::PhysicalGrid) = g.N

"""
    length(g::PhysicalGrid,d::Int) -> Int

Return the total number of cells in grid `g`.
"""
Base.length(g::PhysicalGrid) = prod(size(g))

"""
    limits(g::PhysicalGrid,d::Int) -> Tuple

Return the minimum and maximum physical dimensions in direction `d` for grid `g`.
"""
limits(g::PhysicalGrid,d::Int) = g.xlim[d]

"""
    coordinates(w::Nodes/Edges,g::PhysicalGrid) -> Range

Return coordinate data range for type of `w`.
"""
coordinates(w,g::PhysicalGrid) = coordinates(w,dx=g.Δx,I0=g.I0)

"""
    origin(g::PhysicalGrid) -> Tuple{Int,Int}

Return a tuple of the indices of the primal node that corresponds to the
physical origin of the coordinate system used by `g`. Note that these
indices need not lie inside the range of indices occupied by the grid.
For example, if the range of physical coordinates occupied by the grid
is (1.0,3.0) x (2.0,4.0), then the origin is not inside the grid.
"""
origin(g::PhysicalGrid) = g.I0

"""
    cellsize(g::PhysicalGrid) -> Float64

Return the grid cell size of system `sys`
"""
cellsize(g::PhysicalGrid) = g.Δx


text:
module volume_computation3D_obj

using Modia3D
using Modia3D.ModiaMath
using Modia3D.Unitful

# filename = joinpath(Modia3D.path, "objects", "pascal", "pascal.obj")
# filename = joinpath(Modia3D.path, "objects", "pascal", "convexSolids_pascal.obj", "pascal-0.obj")
# filename = joinpath(Modia3D.path, "objects", "test2", "test2.obj")
 filename = joinpath(Modia3D.path, "objects", "primitiveFiles", "cube.obj")
# filename = joinpath(Modia3D.path, "objects", "test1", "test1.obj")

# Properties
vmat1 = Modia3D.Material(color="LightBlue", transparency=0.5, wireframe=false)
vmat2 = deepcopy(vmat1)
vmat2.transparency = 0.7
smat1 = Modia3D.MassProperties(m=1.0)
font  = Modia3D.Font(fontFamily="TimesNewRoman", color="Black", bold=false, charSize=0.2)
rtext = [0.0, 0.0, 1.0]
ltext = [0.0, 0.0, -0.5]
Rmesh = ModiaMath.rot123(90u"°", 180u"°", 0.0)

LxBox = 10.0; LyBox = 0.01; LzBox = 0.01
DxCap = 0.4; DyCap = 0.3; LzCap = 0.45

DxCyl = 2.0; DyCyl = 2.0; LzCyl = 2.0

# Beam: just a box
#LxBeamBox = 10.0; DyBeamBox = 0.01; LzBeamBox = 0.01

LxBeam = 0.001; DyBeam = 2.0; LzBeam = 2.0

# Solid objects
@assembly Solids begin
    world     = Modia3D.Object3D(visualizeFrame=false)
    #box       = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(LxBox,LyBox,LzBox)      , smat1, vmat1); r=[ 4.5,0.0,0.0])
    #capsule   = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidCapsule(DxCap,LzCap; Dy=DyCap)     , smat1      , vmat1); r=[-3.0,0.0,0.0])
    beam      = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBeam(LxBeam,DyBeam,LzBeam)      , smat1, vmat1); r=[ 3.0,0.0,0.0])
    cylinder  = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidCylinder(DxCyl,LzCyl; Dy=DyCyl) , smat1, vmat1); r=[ 0.0,0.0,0.0])

    #fileMesh1 = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidFileMesh(filename,1.0), smat1      , vmat1);  R=Rmesh, r=[ 0.0,0.0,0.0])
    # fileMesh2 = Modia3D.Object3D(world, Modia3D.SolidWithConvexDecomposition(Modia3D.SolidFileMesh(filename,0.2), smat1, vmat1, vmat2); r=[-1.5,0.0,-2.5], R=Rmesh)
end
solid = Solids(sceneOptions=Modia3D.SceneOptions(visualizeFrames=false, defaultFrameLength=0.7) )
Modia3D.visualizeAssembly!( solid )

#=
println("inertia ", Modia3D.inertiaMatrix(solid.fileMesh1.data.geo,1) )
println("volume ", Modia3D.volume(solid.fileMesh1.data.geo) )
println("centroid ", Modia3D.centroid(solid.fileMesh1.data.geo) )

println("solid.fileMesh1.data.geo.centroidAlgo =  ", solid.fileMesh1.data.geo.centroidAlgo)
println("solid.fileMesh1.data.geo.centroid =  ", solid.fileMesh1.data.geo.centroid)
=#

# println("Modia3D.inertiaMatrix(capsule) = ",   Modia3D.inertiaMatrix(solid.capsule.data.geo,1))

#println("Modia3D.inertiaMatrix(box) = ",      Modia3D.inertiaMatrix(solid.box.data.geo,1))
#println("Modia3D.inertiaMatrix(beam) = ",      Modia3D.inertiaMatrix(solid.beam.data.geo,1))
#println("Modia3D.inertiaMatrix(cylinder) = ",      Modia3D.inertiaMatrix(solid.cylinder.data.geo,1))

println("... success of volume_computation3D_obj.jl!")
end


text:
using LinearAlgebra
function Newton(f, f_prime; x_0, tolerance=1.0E-10, maxiter=10000)
    x_old = x_0
    normdiff = Inf
    iter = 1
    while normdiff > tolerance && iter <= maxiter
        x_new = x_old - (f(x_old) / f_prime(x_old))
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end
    return (value=x_old, normdiff=normdiff, iter=iter)
end

f(x) = (x - 1)^3
f_prime(x) = 3 * (x - 1)^2
g(x) = x^2 + 3 * x - 10
g_prime(x) = 2 * x + 3
sol_f = Newton(f, f_prime, x_0 = 0)
sol_g = Newton(g, g_prime, x_0 = 0)
sol_f, sol_g

max_stars_repo_path: Problem Sets/Newton.jl

max_issues_repo_path: Problem Sets/Newton.jl

max_forks_repo_path: Problem Sets/Newton.jl


text:
#=
VecE3: a 3d euclidean vector
=#

immutable VecE3 <: VecE
    x :: Float64
    y :: Float64
    z :: Float64

    VecE3() = new(0.0,0.0,0.0)
    VecE3(x::Real, y::Real, z::Real) = new(x,y,z)
end

Base.length(::VecE3) = 3
Base.copy(a::VecE3) = VecE3(a.x, a.y, a.z)
Base.convert(::Type{Vector{Float64}}, a::VecE3) = [a.x, a.y, a.z]
function Base.convert{R<:Real}(::Type{VecE3}, a::AbstractArray{R})
    @assert(length(a) == 3)
    VecE3(a[1], a[2], a[3])
end
Base.show(io::IO, a::VecE3) = @printf(io, "VecE3(%.3f, %.3f, %.3f)", a.x, a.y, a.z)

@compat Base.:+(a::VecE3, b::Real) = VecE3(a.x+b, a.y+b, a.z+b)
@compat Base.:+(a::VecE3, b::VecE3) = VecE3(a.x+b.x, a.y+b.y, a.z+b.z)

@compat Base.:-(a::VecE3) = VecE3(-a.x, -a.y, -a.z)
@compat Base.:-(a::VecE3, b::Real) = VecE3(a.x-b, a.y-b, a.z-b)
@compat Base.:-(a::VecE3, b::VecE3) = VecE3(a.x-b.x, a.y-b.y, a.z-b.z)

@compat Base.:*(a::VecE3, b::Real) = VecE3(a.x*b, a.y*b, a.z*b)

@compat Base.:/(a::VecE3, b::Real) = VecE3(a.x/b, a.y/b, a.z/b)

@compat Base.:^(a::VecE3, b::Integer) = VecE3(a.x^b, a.y^b, a.z^b)
@compat Base.:^(a::VecE3, b::AbstractFloat) = VecE3(a.x^b, a.y^b, a.z^b)

# %(a::VecE3, b::Real) = VecE3(a.x%b, a.y%b, a.z%b)

@compat Base.:(==)(a::VecE3, b::VecE3) = isequal(a.x, b.x) && isequal(a.y, b.y) && isequal(a.z, b.z)
Base.isequal(a::VecE3, b::VecE3) = isequal(a.x, b.x) && isequal(a.y, b.y) && isequal(a.z, b.z)

Base.isfinite(a::VecE3) = isfinite(a.x) && isfinite(a.y) && isfinite(a.z)
Base.isinf(a::VecE3) = isinf(a.x) || isinf(a.y) || isinf(a.z)
Base.isnan(a::VecE3) = isnan(a.x) || isnan(a.y) || isnan(a.z)

Base.round(a::VecE3) = VecE3(round(a.x), round(a.y), round(a.z))
Base.floor(a::VecE3) = VecE3(floor(a.x), floor(a.y), floor(a.z))
Base.ceil(a::VecE3) = VecE3(ceil(a.x), ceil(a.y), ceil(a.z))
Base.trunc(a::VecE3) = VecE3(trunc(a.x), trunc(a.y), trunc(a.z))
Base.clamp(a::VecE3, lo::Real, hi::Real) = VecE3(clamp(a.x, lo, hi), clamp(a.y, lo, hi), clamp(a.z, lo, hi))

Base.abs(a::VecE3) = sqrt(a.x*a.x + a.y*a.y + a.z*a.z)
Base.abs2(a::VecE3) = a.x*a.x + a.y*a.y + a.z*a.z
function Base.norm(a::VecE3)
    m = abs(a)
    VecE3(a.x/m, a.y/m, a.z/m)
end

function dist(a::VecE3, b::VecE3)
    Δx = a.x-b.x
    Δy = a.y-b.y
    Δz = a.z-b.z
    sqrt(Δx*Δx + Δy*Δy + Δz*Δz)
end
function dist2(a::VecE3, b::VecE3)
    Δx = a.x-b.x
    Δy = a.y-b.y
    Δz = a.z-b.z
    Δx*Δx + Δy*Δy + Δz*Δz
end

Base.dot(a::VecE3, b::VecE3) = a.x*b.x + a.y*b.y + a.z*b.z
proj(a::VecE3, b::VecE3, ::Type{Float64}) = (a.x*b.x + a.y*b.y + a.z*b.z) / sqrt(b.x*b.x + b.y*b.y + b.z*b.z) # dot(a,b) / |b|
function proj(a::VecE3, b::VecE3, ::Type{VecE3})
    # dot(a,b) / dot(b,b) ⋅ b
    s = (a.x*b.x + a.y*b.y + a.z*b.z) / (b.x*b.x + b.y*b.y + b.z*b.z)
    VecE3(s*b.x, s*b.y, s*b.z)
end

function Base.cross(a::VecE3, b::VecE3)
    x = a.y*b.z - a.z*b.y
    y = a.x*b.z - a.z*b.x
    z = a.x*b.y - a.y*b.x
    VecE3(x,y,z)
end

lerp(a::VecE3, b::VecE3, t::Real) = VecE3(a.x + (b.x-a.x)*t, a.y + (b.y-a.y)*t, a.z + (b.z-a.z)*t)

function rot(a::VecE3, axis::VecE3, θ::Real)

    #=
    rotate the point a: (x,y,z) about the axis: (u,v,w) by the angle θ following the
    right hand rule
    =#

    x,y,z = a.x, a.y, a.z
    u,v,w = axis.x, axis.y, axis.z

    u² = u*u
    v² = v*v
    w² = w*w

    m = u² + v² + w²

    c = cos(θ)
    s = sin(θ)
    mc = 1.0 - c

    rtms = sqrt(m)*s

    magc = m*c
    ux_vy_wz_mc = (u*x + v*y + w*z)*mc

    new_x = (u*ux_vy_wz_mc + x*magc + rtms*(-w*y + v*z)) / m
    new_y = (v*ux_vy_wz_mc + y*magc + rtms*( w*x - u*z)) / m
    new_z = (w*ux_vy_wz_mc + z*magc + rtms*(-v*x + u*y)) / m

    VecE3(new_x, new_y, new_z)
end
function rot_normalized(a::VecE3, axis::VecE3, θ::Real)

    #=
    rotate the point a: (x,y,z) about the axis: (u,v,w) by the angle θ following the
    right hand rule

    Here we assume axis is normalized
    =#

    x,y,z = a.x, a.y, a.z
    u,v,w = axis.x, axis.y, axis.z

    u² = u*u
    v² = v*v
    w² = w*w

    m = u² + v² + w²

    c = cos(θ)
    s = sin(θ)
    mc = 1.0 - c

    ux_vy_wz_mc = (u*x + v*y + w*z)*mc

    new_x = u*ux_vy_wz_mc + x*c + (-w*y + v*z)*s
    new_y = v*ux_vy_wz_mc + y*c + ( w*x - u*z)*s
    new_z = w*ux_vy_wz_mc + z*c + (-v*x + u*y)*s

    VecE3(new_x, new_y, new_z)
end


text:
include("LinearAlgebra/FakeFmpqMat.jl")
include("LinearAlgebra/LatEnum.jl")
include("LinearAlgebra/Howell.jl")
include("LinearAlgebra/Solve.jl")
#include("LinearAlgebra/PseudoMatrix.jl")


text:
module HRV

end


text:
using Documenter, JuliaPetra

makedocs(
    modules = [JuliaPetra],
    sitename="JuliaPetra",
    pages = [
        "Home" => "index.md",
        "Communcation Layer" => "CommunicationLayer.md",
        "Problem Distribution Layer" => "ProblemDistributionLayer.md",
        "Linear Algebra Layer" => "LinearAlgebraLayer.md"
    ])

deploydocs(
    repo = "github.com/Collegeville/JuliaPetra.jl.git"
)


text:
abstract type TauLeapFormula end # <: Function

"""
Gillespie 2001, Eq 26(a)
"""
struct DG2001Eq26a{T1,T2} <: TauLeapFormula
  ξ::T1
  b::T2
  ϵ::Float64

  DG2001Eq26a(ξ::T1, b::T2, ϵ::Real) where {T1,T2} = new{T1,T2}(ξ, b, ϵ)
end

function DG2001Eq26a(number_species::Integer, number_jumps::Integer, ϵ::Real)
  ξ = zeros(number_species)
  b = zeros(number_species, number_jumps)

  return DG2001Eq26a(ξ, b, ϵ)
end

function (formula::DG2001Eq26a)(rates, total_rate)
  ξ = formula.ξ
  b = formula.b
  ϵ = formula.ϵ

  τ = typemax(eltype(rates))

  for j in eachindex(rates)
    proposal = zero(τ)

    for i in eachindex(ξ)
      proposal = proposal + ξ[i] * b[i, j]
    end

    proposal = ϵ * total_rate / abs(proposal)

    if proposal < τ
      τ = proposal
    end
  end

  return τ
end

function update!(formula::DG2001Eq26a, state, model, stoichiometry, rates, total_rate)
  ξ = formula.ξ
  b = formula.b

  mul!(ξ, stoichiometry, rates)

  n, m = size(b)

  for j in 1:m
    for i in 1:n
      b[i,j] = rate_derivative(model, state, i, j)
    end
  end

  return nothing
end

"""
  Gillespie and Petzold 2003, Eq 6
"""
struct DGLP2003Eq6{U,B} <: TauLeapFormula
  V::U
  b::B
  ϵ::Float64
end

function DGLP2003Eq6(number_species, number_jumps, V, ϵ)
  b = zeros(number_species, number_jumps)

  return DGLP2003Eq6(V, b, ϵ)
end

function (formula::DGLP2003Eq6)(rates, total_rate)
  T = eltype(rates)
  V = formula.V
  b = formula.b
  ϵ = formula.ϵ
  τ = typemax(T)

  A1 = ϵ * total_rate
  A2 = A1^2

  for j in eachindex(rates)
    μ, σ² = zero(T), zero(T)
    f = f_calculation(V, b, j)
    for k in eachindex(rates)
      μ  = μ + f * rates[k]
      σ² = σ² + f^2 * rates[k]
    end
    τ′ = min(A1 / abs(μ), A2 / σ²)
    τ  = min(τ, τ′)
  end

  return τ
end

function f_calculation(V::DenseMatrix{Int}, b, j)
  f = zero(eltype(b))

  @simd for i in 1:size(b, 1)
    @inbounds @fastmath f = f + b[i, j] * V[i, j]
  end

  return f
end

function f_calculation(V::SparseMatrixCSC{Int,Int}, b, j)
  f = zero(eltype(b))
  rv = rowvals(V)
  nz = nonzeros(V)

  for i in nzrange(V, j)
    @inbounds @fastmath f = f + b[rv[i], j] * nz[i]
  end

  return f
end

function update!(formula::DGLP2003Eq6, state, model, stoichiometry, rates, total_rate)
  b = formula.b

  n, m = size(b)

  for j in 1:m
    for i in 1:n
      b[i,j] = rate_derivative(model, state, i, j)
    end
  end

  return nothing
end

"""
Unknown formula?
"""
struct GenericLeapFormula
  k::Vector{Float64}
  U::SparseMatrixCSC{Int,Int}
  V::SparseMatrixCSC{Int,Int}
  dxdt::Vector{Float64}
  drdt::Vector{Float64}
  ϵ::Float64
end

# function GenericLeapFormula(number_species::Integer, number_jumps::Integer, U, V, ϵ)
#   dxdt = zeros(number_species)
#   drdt = zeros(number_jumps)

#   return GenericLeapFormula(U, V, dxdt, drdt, ϵ)
# end

function (formula::GenericLeapFormula)(rates, total_rate)
  drdt = formula.drdt
  k = formula.k
  τ = typemax(total_rate)
  ϵ = formula.ϵ

  for j in eachindex(rates)
    A = ϵ * max(rates[j], k[j])
    B = abs(drdt[j])

    τ = min(τ, A / B)
  end

  return τ
end

function update!(formula::GenericLeapFormula, state, model, V, rates, total_rate)
  dxdt = formula.dxdt
  drdt = formula.drdt
  # V = formula.V
  U = formula.U

  update_mean_derivatives!(dxdt, V, rates)
  update_time_derivatives!(drdt, state, model, dxdt)
end

update_mean_derivatives!(dxdt, V, rates) = (dxdt .= V * rates)

function update_time_derivatives!(drdt, state, model, dxdt)
  for j in eachindex(drdt)
    drdt[j]  = 0.0
    for i in eachindex(state)
      drdx = rate_derivative(model, state, i, j)
      drdt[j] = drdt[j] + drdx * dxdt[i]
    end
  end

  return drdt
end

##### post-leap check for valid step #####

struct RejectionThinning{T1,T2,F1,F2,F3,F4}
  threshold::T1
  proposal::T2
  is_accepted::F1
  is_invalid::F2
  execute_leap!::F3
  reverse_leap!::F4
end

function RejectionThinning(threshold::Real, proposal, execute_leap!, reverse_leap!)
  is_accepted = Base.Fix2(<, threshold)     # x -> x < threshold
  is_negative = Base.Fix2(<, 0)             # x -> x < 0
  is_invalid  = Base.Fix1(any, is_negative) # x -> any(is_negative, x)

  RejectionThinning(threshold, proposal, is_accepted, is_invalid, execute_leap!, reverse_leap!)
end

function (f::RejectionThinning)(v::AbstractVector, s::Real)
  # unpack
  threshold = f.threshold
  is_invalid = f.is_invalid
  is_accepted = f.is_accepted
  execute_leap! = f.execute_leap!
  reverse_leap! = f.reverse_leap!
  proposal = f.proposal

  execute_leap!(proposal, v)

  while is_invalid(proposal)
    reverse_leap!(proposal, v)

    # reject events in the proposed leap
    for j in eachindex(v)
      # generate a stream of iid uniform deviates for the rejection step
      event = (rand() for j in 1:v[j])

      # count the number of accepted events
      v[j] = count(is_accepted, event)
    end

    # contract the leap by the threshold parameter
    s = s * threshold

    # apply new update
    execute_leap!(proposal, v)
  end

  return v, s
end

##### extracting stoichiometry matrix #####

function extract_net_stoichiometry(model::ReactionSystem)
  reaction = model.reactions

  number_entries = sum(length(reaction[j].net_change) for j in eachindex(reaction))

  I = zeros(Int, number_entries) # species index
  J = zeros(Int, number_entries) # reaction index
  V = zeros(Int, number_entries) # net stoichiometry

  idx = 1

  for j in eachindex(reaction)
    for (i, v) in reaction[j].net_change

      I[idx] = i
      J[idx] = j
      V[idx] = v

      idx += 1
    end
  end

  return sparse(I, J, V)
end

function extract_coefficients(model::ReactionSystem)
  reaction = model.reactions

  number_entries = sum(length(reaction[j].reactants) for j in eachindex(reaction))

  I = zeros(Int, number_entries) # species index
  J = zeros(Int, number_entries) # reaction index
  V = zeros(Int, number_entries) # net stoichiometry

  idx = 1

  for j in eachindex(reaction)
    for (i, v) in reaction[j].reactants

      I[idx] = i
      J[idx] = j
      V[idx] = v

      idx += 1
    end
  end

  return sparse(I, J, V)
end

##### applying leap updates #####

"""
`ApplyLeapUpdate`

A type-stable closure used to update system state.
Specifically, this type is used to cache the net stoichiometry matrix.

##### Usage

```
execute_leap! = ApplyLeapUpdate(mul!, V)
execute_leap!(x, n) # mul!(x, V, n)
```
"""
struct ApplyLeapUpdate{F,T} <: Function
  f::F
  V::T

  ApplyLeapUpdate(f::F, V::T) where {F,T} = new{F,T}(f, V)
  ApplyLeapUpdate(f::Type{F}, x::T) where {F,T} = new{Type{F},T}(f, V)
end

(g::ApplyLeapUpdate)(x, n) = g.f(x, g.V, n)

function forward_leap!(x, V, n)
  x .+= V*n
end

function backward_leap!(x, V, n)
  x .-= V*n
end


text:
#=
code to handle parsing and manipulating of expressions defining nodes,
factors, naturals, and value functions
=#

"""
Recursively wrap variables in an expression so that in the new expression,
they are projected down to their index tuples in a factor with variables
given in vars.
"""
macro wrapvars(vars, ex, wrapper, varargs...)
    esc(_wrapvars(vars, ex, wrapper, varargs...))
end

_wrapvars(vars, ex, wrapper, varargs...) = ex

function _wrapvars(vars::Vector{Symbol}, ex::Expr, wrapper, varargs...)
    # copy AST
    new_ex = copy(ex)

    # recursively wrap variables
    for i in eachindex(ex.args)
        new_ex.args[i] = _wrapvars(vars, ex.args[i], wrapper, varargs...)
    end

    # return new expression
    new_ex
end

function _wrapvars(vars::Vector{Symbol}, s::Symbol, wrapper, varargs...)
    # if s is a variable in the approved list of vars
    if s in vars
        sym = Expr(:quote, s)
        # :(project(f, $sym, $indtup))
        :($wrapper($sym, $(varargs...)))
    else
        s
    end
end

"""
Recursively make a list of all symbols used as arguments in a given expression.
"""
get_all_syms(x) = Set(Symbol[])
get_all_syms(s::Symbol) = Set(Symbol[s])
function get_all_syms(ex::Expr)
    symset = Set{Symbol}()
    _get_all_syms(ex, symset)
    symset
end
function _get_all_syms(ex::Expr, symset)
    if ex.head in [:call, :quote]
        arglist = ex.args[2:end]
    else
        arglist = ex.args
    end
    for arg in arglist
        if isa(arg, Symbol)
            push!(symset, arg)
        elseif isa(arg, Expr)
            _get_all_syms(arg, symset)
        end
    end
end

"""
Recursively get a list of all index symbols used in a given expression.
"""
get_all_inds(x) = Set(Symbol[])
get_all_inds(s::Symbol) = Set(Symbol[])
function get_all_inds(ex::Expr)
    indset = Set{Symbol}()
    _get_all_inds(ex, indset)
    indset
end
function _get_all_inds(ex::Expr, indset)
    if ex.head == :ref
        for arg in ex.args[2:end]
            if isa(arg, Symbol)
                push!(indset, arg)
            end
        end
    else
        for arg in ex.args
            if isa(arg, Expr)
                _get_all_inds(arg, indset)
            end
        end
    end
end

"""
Remove all indexing from an expression.
"""
strip_inds(x) = x
strip_inds(s::Symbol) = s
function strip_inds(ex::Expr)
    _strip_inds(ex)
end
_strip_inds(x) = x
function _strip_inds(ex::Expr)
    if ex.head == :ref
        out = ex.args[1]
    else
        out = copy(ex)
        out.args = map(_strip_inds, ex.args)
    end
    out
end


text:
# Note that the cartpole scenario can actually be efficiently
# hardcoded as a proper MDP. This implementation is more of a proof
# of concept on how systems can be specified.

mutable struct CartpoleSimulator{M<:Mechanism,S<:MechanismState,I} <: RigidBodySimulator
    mechansim::M
    state::S
    integrator::I
    torques::Vector{Float64}
    time::Float64
    cart_mass::Float64
    pole_mass::Float64
    pole_length::Float64
    pole_com::Float64 # center of mass
    simstep::Float64 # timestep for advancing the physics simulation
    dt::Float64 # timestep for observing next state (also duration for applying an action)
    xlimits::Tuple{Float64,Float64}
end

function CartpoleSimulator(;
        M=1., m=.1, l=1., com=l, xlim::Tuple=(-2.5,2.5),
        dt::Float64=2e-2, simstep::Float64=2e-2, g=-9.8)
    cart_mass = Float64(M)
    pole_mass = Float64(m)
    pole_length = Float64(l)
    pole_com = Float64(com)
    xlimits = map(Float64, xlim)
    gravity = Float64(g)
    CartpoleSimulator(cart_mass, pole_mass, pole_length, pole_com, xlimits, dt, simstep, gravity)
end

function CartpoleSimulator(
        cart_mass::Float64,
        pole_mass::Float64,
        pole_length::Float64,
        pole_com::Float64,
        xlimits::Tuple{Float64,Float64} = (-2.5,2.5),
        dt::Float64 = 2e-2,
        simstep::Float64 = 2e-2,
        gravity::Float64 = -9.8)
    @assert simstep <= dt
    world = RigidBody{Float64}("world")
    mechanism = Mechanism(world; gravity = SVector(0., 0., gravity))

    # create massless slidebar that will restrict cart to x-axis
    inertia1 = SpatialInertia(CartesianFrame3D("slidebar"), zeros(SMatrix{3,3}), zeros(SVector{3}), 0.)
    slidebar = RigidBody(inertia1)
    slidebar_to_world = Joint("slidebar_to_world", Fixed{Float64}())
    before_slidebar_to_world = eye(Transform3D, frame_before(slidebar_to_world), default_frame(world))
    attach!(mechanism, world, slidebar, slidebar_to_world, joint_pose = before_slidebar_to_world)

    # create cart of mass M and connect to slidebar
    inertia2 = SpatialInertia(CartesianFrame3D("cart"), zeros(SMatrix{3,3}), zeros(SVector{3}), cart_mass)
    cart = RigidBody(inertia2)
    cart_to_slidebar = Joint("cart_to_slidebar", Prismatic(SVector(1.,0.,0.)))
    before_cart_to_slidebar = Transform3D(frame_before(cart_to_slidebar), frame_after(slidebar_to_world), zeros(SVector{3}))
    attach!(mechanism, slidebar, cart, cart_to_slidebar, joint_pose = before_cart_to_slidebar)
    # set some joint limits (currently not enforced I think)
    cart_to_slidebar.effort_bounds[1] = RigidBodyDynamics.Bounds{Float64}(-1e3,1e3)
    cart_to_slidebar.position_bounds[1] = RigidBodyDynamics.Bounds{Float64}(xlimits...)
    cart_to_slidebar.velocity_bounds[1] = RigidBodyDynamics.Bounds{Float64}(-5,5)

    # create pole of length l with center of mass in com
    inertia3 = SpatialInertia(CartesianFrame3D("pole"),
                              @SMatrix([0. 0. 0.; 0. pole_mass*pole_length^2 0.; 0. 0. 0.]),
                              pole_mass * SVector(0., 0., pole_com),
                              pole_mass)
    pole = RigidBody(inertia3)
    pole_to_cart = Joint("pole_to_cart", Revolute(SVector(0.,1.,0.)))
    before_pole_to_cart = eye(Transform3D, frame_before(pole_to_cart), frame_after(cart_to_slidebar))
    attach!(mechanism, cart, pole, pole_to_cart, joint_pose = before_pole_to_cart)

    # Set up simulator
    state = MechanismState(mechanism)
    TE = RigidBodyDynamics.cache_eltype(state)
    result = RigidBodyDynamics.DynamicsResult{TE}(state.mechanism)
    torques = zeros(velocity(state))
    closed_loop_dynamics! = (v̇::AbstractArray, ṡ::AbstractArray, t, state) -> begin
        RigidBodyDynamics.dynamics!(result, state, torques)
        #RigidBodyDynamics.zero_torque!(torques, t, state) # reset after applied
        copy!(v̇, result.v̇)
        copy!(ṡ, result.ṡ)
        nothing
    end
    tableau = RigidBodyDynamics.runge_kutta_4(Float64)
    integrator = RigidBodyDynamics.MuntheKaasIntegrator(closed_loop_dynamics!, tableau, ThrowAwaySink())
    OdeIntegrators.set_num_positions!(integrator.stages, length(configuration(state)))
    OdeIntegrators.set_num_velocities!(integrator.stages, length(velocity(state)))
    OdeIntegrators.set_num_additional_states!(integrator.stages, length(additional_state(state)))
    sim = CartpoleSimulator{typeof(mechanism),typeof(state),typeof(integrator)}(
        mechanism,
        state,
        integrator,
        torques,
        0.,
        cart_mass,
        pole_mass,
        pole_length,
        pole_com,
        simstep,
        dt,
        xlimits
    )
    reset!(sim)
    sim
end

function reset!(sim::CartpoleSimulator, s0 = SVector(0.,0.,0.,0.))
    configuration(sim.state)[1] = Float64(s0[1])
    configuration(sim.state)[2] = Float64(s0[3])
    velocity(sim.state)[1] = Float64(s0[2])
    velocity(sim.state)[2] = Float64(s0[4])
    setdirty!(sim.state)
    sim.time = 0.
    sim.torques .= 0.
    sim.state
end

function step!(sim::CartpoleSimulator, force::Number = 0.)
    steps = floor(Int, sim.dt / sim.simstep)
    for i in 1:steps
        @inbounds sim.torques[1] = Float64(force)
        OdeIntegrators.step(sim.integrator, sim.time, sim.state, sim.simstep)
        sim.time += sim.simstep
    end
    sim.state
end


text:
function molecule(::Molecule"benzene")
  return """
C	0.0000000	1.3916730	0.0000000
C	1.2052240	0.6958360	0.0000000
C	1.2052240	-0.6958360	0.0000000
C	0.0000000	-1.3916730	0.0000000
C	-1.2052240	-0.6958360	0.0000000
C	-1.2052240	0.6958360	0.0000000
H	0.0000000	2.4695880	0.0000000
H	2.1387260	1.2347940	0.0000000
H	2.1387260	-1.2347940	0.0000000
H	0.0000000	-2.4695880	0.0000000
H	-2.1387260	-1.2347940	0.0000000
H	-2.1387260	1.2347940	0.0000000
  """
end


text:
module SDPSymmetryReduction

using SparseArrays
using LinearAlgebra

export Partition, admPartSubspace, blockDiagonalize

# Stores a partition of [m]×[m] in a single matrix
# Entries of P should always be 1,…,n
"""
    Partition

A partition subspace. `P.n` is the number of parts, and `P.P` an integer matrix defining the basis elements.
"""
struct Partition
    n::Int64 # Number of parts
    P::Matrix{Int64} # Matrix with entries 1,...,n
end

"""
    roundToZero(f::Number)

Round numbers near zero to zero.
"""
function roundToZero(f::Number)
    if abs(f) < 1e-8
        return 0.0
    end
    return f
end

"""
    part(M::Matrix{T}) where T

Create a partition from the unique entries of `M`.
"""
function part(M::Matrix{T}) where T
    u = unique(M)
    filter!(e -> e ≠ 0, u)
    d = Dict([(u[i], i) for i in eachindex(u)])
    d[0] = 0

    return Partition(size(u, 1), [d[i] for i in M])
end

"""
    coarsestPart(P1::Partition, P2::Partition)

Find the coarsest partition refining `P1` and `P2`.
"""
function coarsestPart(P1::Partition, P2::Partition)
    return part(P1.P * (P2.n + 1) + P2.P)
end

"""
    rndPart(P::Partition)

Returns a random linear combination in the partition space `P`.
"""
function rndPart(P::Partition)
    r = [rand() for i = 1:P.n+1]
    r[1] = 0
    return [r[i+1] for i in P.P]
end

"""
    roundMat(M)

Rounds the matrix, sets entries near 0.0 to 0.0.
"""
function roundMat(M)
    tmp = round.(M, sigdigits = 5)
    tmp = roundToZero.(tmp)
    return tmp
end

"""
    orthProject(A::AbstractMatrix{T}, v::AbstractVector{T}) where T

Project v orthogonally to the span of columns of A
"""
function orthProject(A::AbstractMatrix{T}, v::AbstractVector{T}) where T
    return A * ((A' * A) \ Vector(A' * v))
end

"""
    projectAndRound(M::AbstractMatrix{T}, A::AbstractMatrix{T}; round = true) where T

Projects M to the nullspace of A and rounds the matrix afterwards.
"""
function projectAndRound(M::AbstractMatrix{T}, A::AbstractMatrix{T}; round = true) where T
    n = Int64(sqrt(length(M)))
    tmp = vec(Matrix(M))
    tmp = tmp - orthProject(A, tmp)
    if round
        tmp = roundMat(tmp)
    end
    tmp = reshape(tmp, n, n)
    return Float64.(tmp)
end

"""
    admPartSubspace(C::Vector{T}, A::Matrix{T}, b::Vector{T}, verbose::Bool = false)

Returns the optimal admissible partition subspace for the SDP

``\\inf\\{C^Tx, Ax = b, \\mathrm{Mat}(x) \\succcurlyeq 0\\}.``

This is done using a randomized Jordan-reduction algorithm, and it returns a Jordan algebra (closed under linear combinations and squaring). SDPs can be restricted to such a subspace without changing their optimal value.

## Output
A `Partition` subspace `P`.
"""
function admPartSubspace(C::AbstractVector{T}, A::AbstractMatrix{T}, b::AbstractVector{T}, verbose::Bool = false) where T<:AbstractFloat

    n = Int(sqrt(Float64(length(C))))

    verbose && print("\nStarting the reduction. Original dimension: $(Int64((n^2+n)/2))\n",)

    C = roundMat(C)
    CL = C - orthProject(A', C)
    CL = reshape(roundMat(CL), n, n)
    CL = Matrix(Symmetric(0.5 * (CL + CL')))

    X0 = sparse(qr(sparse(A)) \ b)

    X0Sym = vec(0.5 * (reshape(X0, n, n) + reshape(X0, n, n)'))
    X0Lp = orthProject(A', X0Sym)
    X0Lp = roundMat(X0Lp)
    X0Lp = Matrix(Symmetric(reshape(X0Lp, n, n)))

    # Initialize P as the coarsest partition refining the partitions of two matrices
    P = coarsestPart(part(X0Lp), part(CL))
    dim = P.n

    it = 0

    # Iterate until converged
    while true
        it += 1

        verbose && print("Iteration $it, Current dimension: $dim\n",)

        # Add a random projection to L
        randomizedP = rndPart(P)
        BL = projectAndRound(randomizedP, A')
        P = coarsestPart(P, part(BL))

        # If the last projection didn't change the partition, reuse the random linear combination
        if (dim != P.n)
            randomizedP = rndPart(P)
        end

        # Add random square
        P2 = randomizedP^2
        P2 = roundMat(P2)

        P = coarsestPart(P, part(P2))

        # Check if converged (not changed or reached full partition)
        if (dim == P.n) | (P.n == Int64((n^2 + n) / 2))
            dim = P.n
            break
        end

        dim = P.n
    end

    verbose &&
        print("$it Total iterations, Final dimension: $dim, Old dimension: $(Int64((n^2+n)/2))\n")
    return P
end

"""
    unSymmetrize(P::Partition)

WL algorithm to "un-symmetrize" the Jordan algebra corresponding to `P`.
"""
function unSymmetrize(P::Partition)
    P = deepcopy(P)
    dim = P.n

    it = 0

    # Iterate until converged
    while true
        it += 1

        randomizedP1 = rndPart(P)
        randomizedP2 = rndPart(P)

        P2 = randomizedP1 * randomizedP2
        P2 = roundMat(P2)

        P = coarsestPart(P, part(P2))

        # Check if converged
        if dim == P.n
            break
        end

        dim = P.n
    end
    return P
end

"""
    blockDiagonalize(P::Partition, verbose = true; epsilon = 1e-8, complex = false)

Determines a block-diagonalization of a (Jordan)-algebra given by a partition `P` using a randomized algorithm. `blockDiagonalize(P)` returns a real block-diagonalization `blkd`, if it exists, otherwise `nothing`.

`blockDiagonalize(P; complex = true)` returns the same, but with complex valued matrices, and should be used if no real block-diagonalization was found. To use the complex matrices practically, remember that a Hermitian matrix `A` is positive semidefinite iff `[real(A) -imag(A); imag(A) real(A)]` is positive semidefinite.

## Output

* `blkd.blkSizes` is an integer array of the sizes of the blocks.
* `blkd.blks` is an array of length `P.n` containing arrays of (real/complex) matrices of sizes `blkd.blkSizes`. I.e. `blkd.blks[i]` is the image of the basis element `P.P .== i`.

"""
function blockDiagonalize(P::Partition, verbose = true; epsilon = 1e-8, complex = false)
    P2 = P
    if complex
        P2 = unSymmetrize(P)
    end

    function getRandomMatrix()
        if !complex
            return rndPart(P)
        else
            tmp = rndPart(P2) + im * rndPart(P2)
            return tmp + tmp'
        end
    end

    verbose && println("Determining block sizes...")

    A = getRandomMatrix()
    F = eigen(A)
    Q = F.vectors

    # split by eigenvalues
    roundedEV = round.(F.values, digits = 10)
    uniqueEV = unique(roundedEV)
    countEV = [count(roundedEV .== u) for u in uniqueEV]

    QSplit = [Q[:, [i for i = 1:length(roundedEV) if roundedEV[i] == u]] for u in uniqueEV]

    PSplit = [P.P .== i for i = 1:P.n]

    K = collect(1:length(uniqueEV))
    tmp = getRandomMatrix()
    for i = 1:length(uniqueEV)
        for j = i:length(uniqueEV)
            if K[i] != K[j] && countEV[i] == countEV[j]
                if any(x -> x >= epsilon, abs.(QSplit[i]' * tmp * QSplit[j]))
                    K[K.==K[j]] .= K[i]
                end
            end
        end
    end

    verbose && println("Block sizes are $(sort!([count(K.==Ki) for Ki in unique(K)]))")

    if !complex &&
       sum([count(K .== Ki) * (count(K .== Ki) + 1) / 2 for Ki in unique(K)]) != P.n
        if verbose
            @error("Dimensions do not match up. Rounding error (try different epsilons and/or try again) or not block-diagonalizable over the reals (try parameter complex = true).")
            @show sum([count(K .== Ki) * (count(K .== Ki) + 1) / 2 for Ki in unique(K)])
            @show P.n
        end
        return nothing
    elseif complex && sum([count(K .== Ki)^2 for Ki in unique(K)]) != P2.n
        if verbose
            @error("Dimensions do not match up. Probably a rounding error (try different epsilons and/or try again).")
            @show sum([count(K .== Ki)^2 for Ki in unique(K)])
            @show P2.n
        end
        return nothing
    end

    verbose && println("Determining the algebra-isomorphism...")

    uniqueKs = unique(K)
    reducedQis = []
    for Ki in uniqueKs
        countKi = count(K .== Ki)

        QKi = hcat([QSplit[i] for i = 1:length(uniqueEV) if K[i] == Ki]...)
        B1 = Symmetric(QKi' * getRandomMatrix() * QKi)
        QKi3 = zeros(complex ? Complex{Float64} : Float64, size(B1))

        mult = countEV[Ki]

        for j = 1:countKi
            if j == 1

                QKi3[
                    (mult*(j-1)+1):(mult*(j-1)+mult),
                    (mult*(j-1)+1):(mult*(j-1)+mult),
                ] = Matrix(I, mult, mult)
            else
                QKi3[
                    (mult*(j-1)+1):(mult*(j-1)+mult),
                    (mult*(j-1)+1):(mult*(j-1)+mult),
                ] =
                    B1[
                        1:mult,
                        (mult*(j-1)+1):(mult*(j-1)+mult),
                    ]^(-1)

                QKi3[
                    (mult*(j-1)+1):(mult*(j-1)+mult),
                    (mult*(j-1)+1):(mult*(j-1)+mult),
                ] ./= norm(QKi3[
                    (mult*(j-1)+1),
                    (mult*(j-1)+1):(mult*(j-1)+mult),
                ])

            end
        end

        Per = zeros(size(B1))
        for i = 1:countKi
            for j = 1:mult
                Per[i+countKi*(j-1), j+mult*(i-1)] = 1
            end
        end

        reducedQi = (QKi*QKi3*Per')[:, 1:countKi]

        push!(reducedQis, reducedQi)

    end

    verbose && println("Calculating image of the basis of the algebra...")
    blockDiagonalization = [
        [
            (
                complex ?
                real(B) .* (abs.(real(B)) .>= epsilon) +
                im * imag(B) .* (abs.(imag(B)) .>= epsilon) :
                B .* (abs.(B) .>= epsilon)
            ) for B in [Qi' * P * Qi for Qi in reducedQis]
        ] for P in PSplit
    ]

    blockSizes = [size(b)[1] for b in blockDiagonalization[1]]
    return (blkSizes = blockSizes, blks = blockDiagonalization)
end

end


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModels(Distribution).jl for Static Power System #
# State Estimation.                                                            #
################################################################################
# Contributor: Sander Claeys (@sanderclaeys)                                   #
################################################################################

## Electricity North West Data Parser
# constants
const ENWL_ROOT_PATH    = joinpath(BASE_DIR,"test/data/enwl")
const ENWL_PROFILE_PATH = "$ENWL_ROOT_PATH/profiles"
const ENWL_NETWORK_PATH = "$ENWL_ROOT_PATH/networks"
const ENWL_PROFILES     = Dict(
    ("summer","load") =>"Summer_Load_Profiles.csv",
    ("summer","pv")   =>"Summer_PV_Profiles.csv",
    ("winter","ehp")  =>"Winter_EHP_Profiles.csv",
    ("winter","ev")   =>"Winter_EV_Profiles.csv",
    ("winter","load") =>"Winter_Load_Profiles.csv",
    ("winter","uchp") =>"Winter_uCHP_Profiles.csv",
)
const ENWL_PROFILE_SIGN = Dict(
    ("summer","load") =>  1,
    ("summer","pv")   => -1,
    ("winter","ehp")  =>  1,
    ("winter","ev")   =>  1,
    ("winter","load") =>  1,
    ("winter","uchp") => -1,
)

# functions
get_profiles(ssn, dev) =
    convert(Matrix, _CSV.read("$ENWL_PROFILE_PATH/$(ENWL_PROFILES[(ssn, dev)])"))
function get_enwl_dss_path(network::Int, feeder::Int)
    return "$ENWL_NETWORK_PATH/network_$network/Feeder_$feeder/Master.dss"
end

"""
    insert_profiles!(data, season, devices, pfs; t=missing, useactual=true)

This function adds the load profiles to the parsed ENWL ENGINEERING data `data`.
# Arguments
- data: ENGINEERING data model of the feeder
- season: `"summer"`, `"winter"`
- devices: `"load"`, `"pv"`, `"ev"`, `"ehp"`, `"uchp"`
- pfs: power factor of the devices
- t: time-step

# Example
```julia-repl
julia> data = parse_file(get_enwl_dss_path(10, 1))
julia> insert_profiles!(data, "summer", ["load", "pv"], [0.95, 0.90], t = 144)
```
"""
function insert_profiles!(data, season, devices, pfs; t=missing, useactual=true)
    N = length(data["load"])
    mult = 1E3/data["settings"]["power_scale_factor"]
    sd_tot = zeros(287, N)
    for i in 1:length(devices)
        device = devices[i]; pf = pfs[i]
        linds  = ((1:N).-1).%100 .+1
        sign   = ENWL_PROFILE_SIGN[(season, device)]
        pd     = mult*sign*get_profiles(season, device)[:,linds]
        sd_tot += (1 + im*sqrt(1-pf^2))*pd
    end

    if ismissing(t)
        t = argmax(abs.(reshape(sum(real.(sd_tot), dims=2), 287)))
    end

    for i in 1:N
        load = data["load"]["load$i"]
        dims = length(load["pd_nom"])
        if useactual
            load["pd_nom"] = fill(real(sd_tot[t,i])/dims, dims)
            load["qd_nom"] = fill(imag(sd_tot[t,i])/dims, dims)
        else
            load["pd_nom"] *= real(sd_tot[t,i])
            load["qd_nom"] *= imag(sd_tot[t,i])
        end
    end
end
"""
    rm_enwl_transformer!(data_eng)

This function removes the transformer from a parsed ENWL `ENGINEERING` data file.
"""
function rm_enwl_transformer!(data_eng)
    if haskey(data_eng, "transformer")
        line1 = data_eng["line"]["line1"]
        trans = data_eng["transformer"]["tr1"]
        vprim_scale = trans["vm_nom"][2]/trans["vm_nom"][1]

        vsource = data_eng["voltage_source"]["source"]

        vsource["vm"] *= vprim_scale
        vsource["rs"] *= vprim_scale^2
        vsource["xs"] *= vprim_scale^2
        vsource["bus"] = "1"

        delete!(data_eng, "transformer")
        delete!(data_eng["bus"], "sourcebus")

        vbases_default = data_eng["settings"]["vbases_default"]
        vbases_default["1"] = vbases_default["sourcebus"]*vprim_scale
        delete!(vbases_default, "sourcebus")
    end
end
"""
    reduce_enwl_lines_eng!(data_eng)

This function removes all trailing lines from a parsed ENWL `ENGINEERING` data
file.
"""
function reduce_enwl_lines_eng!(data_eng)
    rm_trailing_lines_eng!(data_eng)
    join_lines_eng!(data_eng)
end
function rm_trailing_lines_eng!(data_eng)

    buses_exclude = []
    for comp_type in ["load", "shunt", "generator", "voltage_source"]
        if haskey(data_eng, comp_type)
            buses_exclude = union(buses_exclude, [comp["bus"] for (_, comp) in data_eng[comp_type]])
        end
    end
    if haskey(data_eng, "transformer")
        buses_exclude = union(buses_exclude, hcat([tr["bus"] for (_, tr) in data_eng["transformer"]]...))
    end

    line_has_shunt = Dict()
    bus_lines = Dict(k=>[] for k in keys(data_eng["bus"]))
    for (id, line) in data_eng["line"]
        lc = data_eng["linecode"][line["linecode"]]
        line_has_shunt[id] = !all(iszero(lc[k]) for k in ["b_fr", "b_to", "g_fr", "g_to"])
        push!(bus_lines[line["f_bus"]], id)
        push!(bus_lines[line["t_bus"]], id)
    end

    eligible_buses = [bus_id for (bus_id, line_ids) in bus_lines if length(line_ids)==1 && !(bus_id in buses_exclude) && !line_has_shunt[line_ids[1]]]

    while !isempty(eligible_buses)
        for bus_id in eligible_buses
            # this trailing bus has one associated line
            line_id = bus_lines[bus_id][1]
            line = data_eng["line"][line_id]

            delete!(data_eng["line"], line_id)
            delete!(data_eng["bus"],  bus_id)

            other_end_bus = line["f_bus"]==bus_id ? line["t_bus"] : line["f_bus"]
            bus_lines[other_end_bus] = setdiff(bus_lines[other_end_bus], [line_id])
            delete!(bus_lines,  bus_id)
        end

        eligible_buses = [bus_id for (bus_id, line_ids) in bus_lines if length(line_ids)==1 && !(bus_id in buses_exclude) && !line_has_shunt[line_ids[1]]]
    end
end
function _line_reverse_eng!(line)
    prop_pairs = [("f_bus", "t_bus")]

    for (x,y) in prop_pairs
        tmp = line[x]
        line[x] = line[y]
        line[y] = tmp
    end
end
function join_lines_eng!(data_eng)
    # a bus is eligible for reduction if it only appears in exactly two lines
    buses_all = collect(keys(data_eng["bus"]))
    buses_exclude = []

    # start by excluding all buses that appear in components other than lines
    for comp_type in ["load", "shunt", "generator", "voltage_source"]
        if haskey(data_eng, comp_type)
            buses_exclude = union(buses_exclude, [comp["bus"] for (_, comp) in data_eng[comp_type]])
        end
    end

    # per bus, list all inbound or outbound lines
    bus_lines = Dict(bus=>[] for bus in buses_all)
    for (id, line) in data_eng["line"]
        push!(bus_lines[line["f_bus"]], id)
        push!(bus_lines[line["t_bus"]], id)
    end

    # exclude all buses that do not have exactly two lines connected to it
    buses_exclude = union(buses_exclude, [bus for (bus, lines) in bus_lines if length(lines)!=2])

    # now loop over remaining buses
    candidates = setdiff(buses_all, buses_exclude)
    for bus in candidates
        line1_id, line2_id = bus_lines[bus]
        line1 = data_eng["line"][line1_id]
        line2 = data_eng["line"][line2_id]

        # reverse lines if needed to get the order
        # (x)--fr-line1-to--(bus)--to-line2-fr--(x)
        if line1["f_bus"]==bus
            _line_reverse_eng!(line1)
        end
        if line2["f_bus"]==bus
            _line_reverse_eng!(line2)
        end

        reducable = true
        reducable = reducable && line1["linecode"]==line2["linecode"]
        reducable = reducable && all(line1["t_connections"].==line2["t_connections"])
        if reducable

            line1["length"] += line2["length"]
            line1["t_bus"] = line2["f_bus"]
            line1["t_connections"] = line2["f_connections"]

            delete!(data_eng["line"], line2_id)
            delete!(data_eng["bus"], bus)
            for x in candidates
                if line2_id in bus_lines[x]
                    bus_lines[x] = [setdiff(bus_lines[x], [line2_id])..., line1_id]
                end
            end
        end
    end

    return data_eng
end
"""
    reduce_enwl_lines_math!(data_math)

This function removes all trailing lines from a parsed ENWL `MATHEMATICAL` data
file.
"""
function reduce_enwl_lines_math!(data_math)

    # a bus is eligible for reduction if it only appears in exactly two lines
    buses_all = collect(keys(data_math["bus"]))
    buses_exclude = []

    # start by excluding all buses that appear in components other than lines
    for comp_type in ["load", "shunt", "gen"]
        if haskey(data_math, comp_type)
            buses_exclude = union(buses_exclude, [string(comp["$(comp_type)_bus"]) for (_, comp) in data_math[comp_type]])
        end
    end

    # per bus, list all inbound or outbound lines
    bus_lines = Dict(bus=>[] for bus in buses_all)
    for (id, line) in data_math["branch"]
        push!(bus_lines[string(line["f_bus"])], id)
        push!(bus_lines[string(line["t_bus"])], id)
    end

    # exclude all buses that do not have exactly two lines connected to it
    buses_exclude = union(buses_exclude, [bus for (bus, lines) in bus_lines if length(lines)!=2])

    # now loop over remaining buses
    candidates = setdiff(buses_all, buses_exclude)
    for bus in candidates
        line1_id, line2_id = bus_lines[bus]
        line1 = data_math["branch"][line1_id]
        line2 = data_math["branch"][line2_id]

        # reverse lines if needed to get the order
        # (x)--fr-line1-to--(bus)--to-line2-fr--(x)
        if line1["f_bus"]==bus
            _line_reverse!(line1)
        end
        if line2["f_bus"]==bus
            _line_reverse!(line2)
        end

        reducable = true
        reducable = reducable && iszero(line1["g_to"]) && iszero(line1["b_to"])
        reducable = reducable && iszero(line2["g_fr"]) && iszero(line2["b_fr"])
        if reducable
            line1["br_r"] = line1["br_r"] .+ line2["br_r"]
            line1["br_x"] = line1["br_x"] .+ line2["br_x"]
            line1["g_to"] = line2["g_fr"]
            line1["b_to"] = line2["b_fr"]
            line1["t_bus"] = line2["f_bus"]

            delete!(data_math["branch"], line2_id)
            for x in candidates
                if line2_id in bus_lines[x]
                    bus_lines[x] = [setdiff(bus_lines[x], [line2_id])..., line1_id]
                end
            end
        end
    end

    return data_math
end
function _line_reverse!(line)
    prop_pairs = [("f_bus", "t_bus"), ("g_fr", "g_to"), ("b_fr","b_to")]

    for (x,y) in prop_pairs
        tmp = line[x]
        line[x] = line[y]
        line[y] = tmp
    end
end


text:
struct ComputExpDB{AT, PT, PFT, LVT, IT, ET, #=EPT, EIT,=# RT, IPT, IRT, RIRT, RIPT}
    algorithms_table::AT
    problems_table::PT
    problem_features_table::PFT
    language_versions_table::LVT
    implementations_table::IT
    experiments_table::ET
    #experiment_problem_table::EPT
    #experiment_implementation_table::EIT
    runs_table::RT
    implementation_parameters_table::IPT
    implementation_results_table::IRT
    run_implementation_result_table::RIRT
    run_implementation_parameter_table::RIPT
end

import Base: show, display
show(io::IO, db::ComputExpDB) = println("ComputExpDB")

function ComputExpDB()
    algorithms_table = table(@NT(ID = Int[], Name = String[],
        Description = String[]), pkey=:ID)

    problems_table = table(@NT(ID = Int[], Description = String[]),
        pkey=:ID)

    problem_features_table = table(@NT(ID = Int[], Problem_ID = Int[],
        Feature_symbol = Symbol[], Feature_value = Any[]), pkey=:ID)

    language_versions_table = table(@NT(ID = Int[], Name = String[],
        Version = String[], Commit = U(String)[]), pkey=:ID)

    implementations_table = table(@NT(ID = Int[], Algorithm_ID = Int[],
        Language_ID = Int[], Repository = String[], Commit = U(String)[],
        Script = String[], Function = Symbol[]), pkey = :ID)

    experiments_table = table(@NT(ID = Int[], Directory = String[],
        Execute = Bool[]), pkey = :ID)

    #experiment_problem_table = table(@NT(ID = Int[], Experiment_ID = Int[],
    #    Problem_ID = Int[]), pkey = :ID)

    #experiment_implementation_table = table(@NT(ID = Int[],
    #    Experiment_ID = Int[], Implementation_ID = Int[]), pkey = :ID)

    runs_table = table(@NT(ID = Int[], Experiment_ID = Int[],
        Problem_ID = Int[], Implementation_ID = Int[], Executing = Bool[],
        Executed = Bool[]), pkey = :ID)

    implementation_parameters_table = table(@NT(ID = Int[],
        Implementation_ID = Int[], Parameter_symbol = Symbol[],
        Parameter_description = String[]), pkey = :ID)

    implementation_results_table = table(@NT(ID = Int[],
        Implementation_ID = Int[], Result_name = String[]), pkey = :ID)

    run_implementation_result_table = table(@NT(ID = Int[],
        Implementation_result_ID = Int[], Run_ID = [],
        Value = Any[]), pkey = :ID)

    run_implementation_parameter_table = table(@NT(ID = Int[],
        Implementation_parameter_ID = Int[], Run_ID = [],
        Value = Any[]), pkey = :ID)

    return ComputExpDB(algorithms_table, problems_table,
        problem_features_table, language_versions_table,
        implementations_table, experiments_table#=, experiment_problem_table=#,
        #=experiment_implementation_table,=# runs_table,
        implementation_parameters_table, implementation_results_table,
        run_implementation_result_table, run_implementation_parameter_table)
end


text:
using Test
using ConvDiffMIPDECO
using jInv.ForwardShare
using jInv.LinearSolvers
using KrylovMethods
using jInv.Mesh
using SparseArrays
using LinearAlgebra
using MUMPSjInv

domain = [0 1. 0 3. 0 2.]

N = [8;16;32]
v = randn(3,1)
v /= norm(v)
V = X -> ones(size(X,1))*v'
bc = (:neu,:dir,:dir,:neu,:dir,:dir)

u  = X -> cos.(pi*X[:,1]).*cos.(pi*X[:,2]).*cos.(pi*X[:,3])
ux = X -> -pi*sin.(pi*X[:,1]).*cos.(pi*X[:,2]).*cos.(pi*X[:,3])
uy = X -> -pi*cos.(pi*X[:,1]).*sin.(pi*X[:,2]).*cos.(pi*X[:,3])
uz = X -> -pi*cos.(pi*X[:,1]).*cos.(pi*X[:,2]).*sin.(pi*X[:,3])
sig = .05;

f  = X->  3*sig*pi^2*u(X)+v[1]*ux(X)+v[2]*uy(X)+v[3]*uz(X)

err = zeros(length(N))

times = zeros(length(N),3)
for j=1:3
	for k=1:length(N)
		Ainvs = (getJuliaSolver(), getMUMPSsolver(),ConvDiffMIPDECO.getBICGSTB())

		nk      = [N[k];N[k];N[k]]
		Mk  = getRegularMesh(domain,nk)
		xck = getCellCenteredGrid(Mk)

		pFork   = getConvDiffParam(Mk,V,sig=sig,gd=u,bc=bc,Ainv=Ainvs[j])

		dobsk,pFork = getData(vec(f(xck)), pFork)

		utruek = u(xck)
		times[k,j] = @elapsed begin
			err[k] = norm(pFork.Fields - utruek,Inf)/norm(utruek,Inf)
		end
		println("err[$k] = $(err[k])")
	end
	@test all(abs.(diff(log2.(err))).>0.7)
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO