text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    LRNCType

Enum representing the column types for LRNData:

`ignore = 0, data = 1, class = 3, key = 9`
"""
@enum LRNCType begin
    ignore = 0
    data = 1
    class = 3
    key = 9
end

"""
    LRNData

`LRNData` represents the contents of a `*.lrn` file.

# Fields
- `data::Matrix{Float64}`: Matrix of data, cases in rows, variables in columns.
- `column_types::Array{LRNCType, 1}`: Column types, see `LRNCType`.
- `key::Array{Integer, 1}`: Unique key for each line.
- `names::Array{String, 1}`: Column names.
- `key_name::String`: Name for key column.
- `comment::String`: Comments about the data.
"""
struct LRNData
    data::Matrix{Float64}
    column_types::Array{LRNCType,1}
    key::Array{Int64,1}
    names::Array{String,1}
    key_name::String
    comment::String

    function LRNData(data::Matrix{Float64}, column_types, key, names, key_name, comment)
        (nrow, ncol) = size(data)
        # Enforcing invariants
        @assert length(names) == ncol
        @assert length(key) == nrow
        @assert length(column_types) == ncol
        @assert allunique(key)
        @assert all(i -> 1 <= i <= nrow, key)
        new(data, column_types, key, names, key_name, comment)
    end
end

"""
    LRNData(data::AbstractMatrix{Float64})

Convenience constructor for `LRNData`. Uses sensible defaults:
```julia
column_types = [9,1,1...]
key          = [1,2,3...]
names        = ["C1","C2","C3"...]
key_name     = "Key"
comment      = ""
```
"""
function LRNData(
    data::AbstractMatrix{Float64};
    column_types = LRNCType.(fill(1, size(data, 2))),
    key          = collect(1:size(data, 1)),
    names        = fill("C", size(data, 2)) .* string.(1:size(data, 2)),
    key_name     = "Key",
    comment      = "",
)
    LRNData(data, column_types, key, names, key_name, comment)
end

Base.firstindex(D::LRNData) = 1
Base.lastindex(D::LRNData) = size(D, 2)
Base.getindex(D::LRNData, i::Int) = D.data[D.key[i], :]
Base.getindex(D::LRNData, I) = [D[i] for i in I]

"""
    writeLRN(filename::String, lrn::LRNData, directory=pwd())

Write the contents of a `LRNData` struct into a file.
"""
function writeLRN(lrn::LRNData, filename::String, directory = pwd())
    (nrow, ncol) = size(lrn.data)
    filename = prepare_path(filename, "lrn", directory)
    open(filename, "w") do f
        # write comment if it exists
        if !isempty(lrn.comment)
            write(f, "# $(lrn.comment)\n#\n")
        end
        # write number of cases
        write(f, "% $(nrow)\n")
        # write number of variables (+1 for key)
        write(f, "% $(ncol+1)\n")
        # write column types (9 for key)
        write(f, "% $(Int(key))\t$(join(map(Int,lrn.column_types), '\t'))\n")
        # write names
        write(f, "% $(lrn.key_name)\t$(join(lrn.names, '\t'))\n")
        # write data
        for (index, row) in enumerate(eachrow(lrn.data))
            new_row = replace(row, Inf => NaN)
            write(f, "$(lrn.key[index])\t$(join(new_row,'\t'))\n")
        end
    end
end

"""
    readLRN(filename::String, directory=pwd())

Read the contents of a `*.lrn` and return a `LRNData` struct.
"""
function readLRN(filename::String, directory = pwd())
    data = []
    column_types = []
    key = []
    names = []
    key_name = ""
    comment = ""
    key_index = 0

    # There is currently no way to have a native file chooser dialogue
    # without building a whole Gtk/Tk GUI

    filename = prepare_path(filename, "lrn", directory)
    open(filename, "r") do f
        line = readline(f)
        # Comments
        while startswith(line, '#')
            comment *= strip(line, [' ', '\t', '#'])
            line = readline(f)
        end
        strip_header = l -> strip(l, [' ', '\t', '%'])
        # Number of datasets
        nrow = parse(Int, strip_header(line))
        line = readline(f)
        # Number of columns
        ncol = parse(Int, strip_header(line))
        line = readline(f)
        # Column types
        column_types = split(strip_header(line), '\t')
        column_types = map(x -> LRNCType(parse(Int, x)), column_types)
        key_index = findfirst(x -> x == DataIO.key, column_types)
        deleteat!(column_types, key_index)
        line = readline(f)
        # Names
        names = split(strip_header(line), '\t')
        key_name = names[key_index]
        deleteat!(names, key_index)
        # Data
        data = readdlm(f, Float64, skipblanks = true)
        key = map(Int, data[:, key_index])
        data = data[:, deleteat!(collect(1:ncol), key_index)] # remove key column
    end

    LRNData(data, column_types, key, names, key_name, comment)
end

max_stars_repo_name: ckafi/DataIO

max_issues_repo_name: ckafi/DataIO

max_forks_repo_name: ckafi/DataIO


text:
const FOLDER_LIST =
    ["CU",
    "EV_CHARGE",
    "EV_DISCHARGE",
    "EV_INF",
    "EV_L",
    "EV_PHEVFUEL",
    "G",
    "G_INF",
    "H2_G2P",
    "H2_GS_IN",
    "H2_GS_L",
    "H2_GS_OUT",
    "H2_P2G",
    "HEAT_HP_IN",
    "HEAT_INF",
    "HEAT_STO_L",
    "N",
    "N_G2P",
    "N_GS",
    "N_P2G",
    "N_STO_E",
    "N_STO_P",
    "STO_IN",
    "STO_L",
    "STO_OUT"]

"""
Remove subfolders of `rdir`, the results directory.
Note that the default is for individual variable data files to be saved to the directory
".../<save_name>/single_results/<name>"
when the function save_results(dtr::DieterModel,"<save_name>") has been called.
"""
function rm_folders(rdir::AbstractString; singledir::AbstractString="single_results")
    for folder in FOLDER_LIST
        path = joinpath(rdir,singledir,folder)
        isdir(path) && rm(path, recursive=true)
    end
end

function merge_results(rdir::AbstractString; singledir::AbstractString="single_results")

    prog = Progress(length(FOLDER_LIST), 0.01, "Merging results...   ", 40)

    for folder in FOLDER_LIST
        f = joinpath(rdir, singledir, folder)
        if isdir(f)
            dfs = [Feather.read(joinpath(f,file))
                for file in readdir(f) if rsplit(file,".";limit=2)[end] == "feather"]
            merged_df = vcat(dfs...)
            Feather.write(joinpath(rdir,folder*".feather"), merged_df)
        end
        next!(prog)
    end

    return nothing
end

function join_dfs(arr::Vector{DataFrame}; on::Array{Symbol, 1})

    merge_df = join(arr[1], arr[2], on=on)
    for df in arr[3:end]
        merge_df = join(merge_df, df, on=on, kind = :outer)
    end

    return merge_df
end

function need_to_merge(rdir::AbstractString; singledir::AbstractString="single_results")

    for f in FOLDER_LIST
        p = joinpath(rdir, singledir, f)
        pp = readdir(p)
        m = maximum([mtime(joinpath(p,x)) for x in pp])
        merged_file = joinpath(rdir, f*".feather")

        if isfile(merged_file)
            m2 = mtime(merged_file)
        else
            println("No merged file yet.")
            return true
        end

        if m2 < m
            println("Newer version of $f found.")
            return true
        end
    end

    return false
end

const COLS = [:timestep, :scen, :min_res, :ev, :heat, :h2]

function create_generation_summary(rdir::AbstractString)

    g = Feather.read(joinpath(rdir,"G.feather"))
    cu = Feather.read(joinpath(rdir,"CU.feather"))
    sto_in = Feather.read(joinpath(rdir,"STO_IN.feather"))
    sto_out = Feather.read(joinpath(rdir,"STO_OUT.feather"))

    ev_charge = Feather.read(joinpath(rdir,"EV_CHARGE.feather"))
    ev_discharge = Feather.read(joinpath(rdir,"EV_DISCHARGE.feather"))

    heat_demand = Feather.read(joinpath(rdir,"HEAT_HP_IN.feather"))

    p2g = Feather.read(joinpath(rdir,"H2_P2G.feather"))
    g2p = Feather.read(joinpath(rdir,"H2_G2P.feather"))

    gen_summary = DataFrame[]
    push!(gen_summary, by(g, COLS, Generation_total = :Value => sum))

    gen_by_fuel = by(g, vcat(COLS,:fuel), Generation_by_fuel = :Value => sum)
    push!(gen_summary, unstack(gen_by_fuel,  :fuel, :Generation_by_fuel))

    gen_by_res = by(g, vcat(COLS,:renewable), Generation_by_res = :Value => sum)
    gen_by_res= unstack(gen_by_res,  :renewable, :Generation_by_res)
    rename!(gen_by_res, Dict(Symbol(:false) => :conventional, Symbol(:true) => :renewable))
    push!(gen_summary, gen_by_res)

    push!(gen_summary, by(cu, COLS, Curtailment = :Value => sum))

    push!(gen_summary, by(ev_charge, COLS, Charged_by_ev = :Value => sum))
    push!(gen_summary, by(ev_discharge, COLS, Discharged_by_ev = :Value => sum))

    push!(gen_summary, by(heat_demand, COLS, Demand_from_hp = :Value => sum))

    push!(gen_summary, by(sto_in, COLS, Storage_in = :Value => sum))
    push!(gen_summary, by(sto_out, COLS, Storage_out = :Value => sum))

    push!(gen_summary, by(p2g, COLS, Power2Gas = :Value => sum))
    push!(gen_summary, by(g2p, COLS, Gas2power = :Value => sum))

    summary_generation = join_dfs(gen_summary, on=COLS)

    summary_generation[!,:res_share] =
        summary_generation[!,:renewable].*100 ./ summary_generation[!,:Generation_total]

    summary_generation[!,:relative_curtailment] =
        summary_generation[!,:Curtailment].*100 ./ summary_generation[!,:Generation_total]

    Feather.write(joinpath(rdir,"summary_generation.feather"), summary_generation)
    CSV.write(joinpath(rdir,"summary_generation.csv"), summary_generation)

end

function create_investment_summary(rdir::AbstractString)

    inv = Feather.read(joinpath(rdir,"N.feather"))
    inv_power = Feather.read(joinpath(rdir,"N_STO_P.feather"))
    inv_storage = Feather.read(joinpath(rdir,"N_STO_E.feather"))

    n_p2g = Feather.read(joinpath(rdir,"N_P2G.feather"))
    n_g2p = Feather.read(joinpath(rdir,"N_G2P.feather"))
    n_gs = Feather.read(joinpath(rdir,"N_GS.feather"))

    inv_summary = DataFrame[]

    inv = by(inv, vcat(COLS,:fuel), Investment_by_fuel = :Value => sum)
    push!(inv_summary, unstack(inv,  :fuel, :Investment_by_fuel))
    push!(inv_summary, by(inv_power, COLS, Investment_power = :Value => sum))
    push!(inv_summary, by(inv_storage, COLS, Investment_storage = :Value => sum))

    push!(inv_summary, by(n_p2g, COLS, Investment_p2g = :Value => sum))
    push!(inv_summary, by(n_g2p, COLS, Investment_g2p = :Value => sum))
    push!(inv_summary, by(n_gs, COLS, Investment_gs = :Value => sum))

    summary_investment = join_dfs(inv_summary, on=COLS)

    Feather.write(joinpath(rdir,"summary_investment.feather"), summary_investment)
    CSV.write(joinpath(rdir,"summary_investment.csv"), summary_investment)

end

function post_process_results(rdir::AbstractString)

    if need_to_merge(rdir)
        merge_results(rdir)
        create_generation_summary(rdir)
        create_investment_summary(rdir)
    else
        println("Did not find any new files.")
    end

    return nothing
end

ext: jl


text:
include("./../../src/Lar.jl")

include("./randomarrangement2d.jl")
include("./randomshapes.jl")
include("./stresstest2d.jl")
include("./svg2lar.jl")

randomarrangement2d(100)
randomshapes(100)
stresstest2d()

#svgarrangement("./test/svg/curved.svg")
svgarrangement("./test/svg/holes.svg")
#svgarrangement("./test/svg/howto.svg")
#svgarrangement("./test/svg/howto2.svg")
#svgarrangement("./test/svg/Lar.svg")
#svgarrangement("./test/svg/new.svg")
#svgarrangement("./test/svg/paths.svg")
#svgarrangement("./test/svg/twopaths.svg")


text:
# # NICE network
# For the definition of this network and concepts of normalizing flow,
# please refer this nice blog: https://lilianweng.github.io/lil-log/2018/10/13/flow-based-deep-generative-models.html,
# and the pytorch notebook: https://github.com/GiggleLiu/marburg/blob/master/notebooks/nice.ipynb

using NiLang, NiLang.AD
using LinearAlgebra
using DelimitedFiles
using Plots

# `include` the optimizer, you can find it under the `Adam.jl` file in the `examples/` folder.
include(NiLang.project_relative_path("examples", "Adam.jl"))

# ## Model definition
# First, define the single layer transformation and its behavior under `GVar` - the gradient wrapper.
struct NiceLayer{T}
    W1::Matrix{T}
    b1::Vector{T}
    W2::Matrix{T}
    b2::Vector{T}
    y1::Vector{T}
    y1a::Vector{T}
end

"""Apply a single NICE transformation."""
@i function nice_layer!(x::AbstractVector{T}, layer::NiceLayer{T},
                y!::AbstractVector{T}) where T
    @routine @invcheckoff begin
        i_affine!(layer.y1, layer.W1, layer.b1, x)
        @inbounds for i=1:length(layer.y1)
            if (layer.y1[i] > 0, ~)
                layer.y1a[i] += layer.y1[i]
            end
        end
    end
    i_affine!(y!, layer.W2, layer.b2, layer.y1a)
    ~@routine
    ## clean up accumulated rounding error, since this memory is reused.
    @safe layer.y1 .= zero(T)
end

# Here, in each layer, we use the information in `x` to update `y!`.
# During computing, we use the `y1` and `y1a` fields of the network as ancilla space,
# both of them can be uncomputed at the end of the function.
# However, we need to erase small numbers to make sure the rounding error does not accumulate.

# A nice network always transforms inputs reversibly.
# We update one half of `x!` a time, so that input and output memory space do not clash.
const NiceNetwork{T} = Vector{NiceLayer{T}}

"""Apply a the whole NICE network."""
@i function nice_network!(x!::AbstractVector{T}, network::NiceNetwork{T}) where T
    @invcheckoff for i=1:length(network)
        np ← length(x!)
        if (i%2==0, ~)
            @inbounds nice_layer!(x! |> subarray(np÷2+1:np), network[i], x! |> subarray(1:np÷2))
        else
            @inbounds nice_layer!(x! |> subarray(1:np÷2), network[i], x! |> subarray(np÷2+1:np))
        end
        np → length(x!)
    end
end

function random_nice_network(nparams::Int, nhidden::Int, nlayer::Int; scale=0.1)
    random_nice_network(Float64, nparams, nhidden, nlayer; scale=scale)
end

function random_nice_network(::Type{T}, nparams::Int, nhidden::Int, nlayer::Int; scale=0.1) where T
    nin = nparams÷2
    scale = T(scale)
    y1 = zeros(T, nhidden)
    NiceLayer{T}[NiceLayer(randn(T, nhidden, nin)*scale, randn(T, nhidden)*scale,
            randn(T, nin, nhidden)*scale, randn(T, nin)*scale, y1, zero(y1)) for _ = 1:nlayer]
end

# ## Parameter management

nparameters(n::NiceLayer) = length(n.W1) + length(n.b1) + length(n.W2) + length(n.b2)
nparameters(n::NiceNetwork) = sum(nparameters, n)

"""collect parameters in the `layer` into a vector `out`."""
function collect_params!(out, layer::NiceLayer)
    a, b, c, d = length(layer.W1), length(layer.b1), length(layer.W2), length(layer.b2)
    out[1:a] .= vec(layer.W1)
    out[a+1:a+b] .= layer.b1
    out[a+b+1:a+b+c] .= vec(layer.W2)
    out[a+b+c+1:end] .= layer.b2
    return out
end

"""dispatch vectorized parameters `out` into the `layer`."""
function dispatch_params!(layer::NiceLayer, out)
    a, b, c, d = length(layer.W1), length(layer.b1), length(layer.W2), length(layer.b2)
    vec(layer.W1) .= out[1:a]
    layer.b1 .= out[a+1:a+b]
    vec(layer.W2) .= out[a+b+1:a+b+c]
    layer.b2 .= out[a+b+c+1:end]
    return layer
end

function collect_params(n::NiceNetwork{T}) where T
    out = zeros(T, nparameters(n))
    k = 0
    for layer in n
        np = nparameters(layer)
        collect_params!(view(out, k+1:k+np), layer)
        k += np
    end
    return out
end

function dispatch_params!(network::NiceNetwork, out)
    k = 0
    for layer in network
        np = nparameters(layer)
        dispatch_params!(layer, view(out, k+1:k+np))
        k += np
    end
    return network
end

# ## Loss function

# To obtain the log-probability of a data.

@i function logp!(out!::T, x!::AbstractVector{T}, network::NiceNetwork{T}) where T
    (~nice_network!)(x!, network)
    @invcheckoff for i = 1:length(x!)
        @routine begin
            xsq ← zero(T)
            @inbounds xsq += x![i]^2
        end
        out! -= 0.5 * xsq
        ~@routine
    end
end

# The negative-log-likelihood loss function

@i function nice_nll!(out!::T, cum!::T, xs!::Matrix{T}, network::NiceNetwork{T}) where T
    @invcheckoff for i=1:size(xs!, 2)
        @inbounds logp!(cum!, xs! |> subarray(:,i), network)
    end
    out! -= cum!/(@const size(xs!, 2))
end

# ## Training

function train(x_data, model; num_epochs = 800)
    num_vars = size(x_data, 1)
    params = collect_params(model)
    optimizer = Adam(; lr=0.01)
    for epoch = 1:num_epochs
        loss, a, b, c = nice_nll!(0.0, 0.0, copy(x_data), model)
        if epoch % 50 == 1
            println("epoch = $epoch, loss = $loss")
            display(showmodel(x_data, model))
        end
        _, _, _, gmodel = (~nice_nll!)(GVar(loss, 1.0), GVar(a), GVar(b), GVar(c))
        g = grad.(collect_params(gmodel))
        update!(params, grad.(collect_params(gmodel)), optimizer)
        dispatch_params!(model, params)
    end
    return model
end

function showmodel(x_data, model; nsamples=2000)
    scatter(x_data[1,1:nsamples], x_data[2,1:nsamples]; xlims=(-5,5), ylims=(-5,5))
    zs = randn(2, nsamples)
    for i=1:nsamples
        nice_network!(view(zs, :, i), model)
    end
    scatter!(zs[1,:], zs[2,:])
end

# you can find the training data in `examples/` folder
x_data = Matrix(readdlm(NiLang.project_relative_path("examples", "train.dat"))')

import Random; Random.seed!(22)
model = random_nice_network(Float64, size(x_data, 1), 10, 4; scale=0.1)

# Before training, the distribution looks like
# ![before](../../asset/nice_before.png)
model = train(x_data, model; num_epochs=800)

# After training, the distribution looks like
# ![before](../../asset/nice_after.png)


text:
export
    nexthailstone, hailstonelength, allhailstone

nexthailstone(n::Integer) = begin
    n <= 0 && throw(DomainError())
    iseven(n) ? n ÷ 2 : 3*n + 1
end

hailstonelengthcalc(x, cache) = 1 + hailstonelength(nexthailstone(x), cache)
hailstonelength(n::T, cache::Array{T,1} = zeros(T, 0)) where {T<:Integer} = begin
    n <= 0 && throw(DomainError())
    n == 1 && return one(T)

    n > length(cache) && return hailstonelengthcalc(n, cache)
    (ans = cache[n]) > 0 && return ans
    cache[n] = hailstonelengthcalc(n, cache)
end

struct HailstoneIterator{T<:Integer}
    n::T
end

allhailstone(n::Integer) = begin
    n <= 0 && throw(DomainError())
    HailstoneIterator(n)
end

Base.start(it::HailstoneIterator) = it.n
Base.next(::HailstoneIterator{T}, state) where {T<:Integer} =
    (state, state == 1 ? zero(T) : nexthailstone(state))
Base.done(::HailstoneIterator, state) = state <= 0

Base.eltype(it::HailstoneIterator) = Base.eltype(typeof(it))
Base.eltype(::Type{HailstoneIterator{T}}) where {T} = T

Base.iteratorsize(::HailstoneIterator) = Base.SizeUnknown()


text:
"""
    savetaskmodel(path, task, model[; force = false])

Save a trained `model` along with a `task` to `path` for later inference.
Use [`loadtaskmodel`](#) for loading both back into a session. If `path`
already exists, only write to it if `force = true`.

If `model` weights are on a GPU, they will be moved to the CPU before saving
so they can be loaded in a non-GPU environment.

[JLD2.jl](https://github.com/JuliaIO/JLD2.jl) is used for serialization.
"""
function savetaskmodel(path, task::LearningTask, model; force = false)
    if !force && isfile(path)
        error("$path already exists. Use `force = true` to overwrite.")
    end
    jldsave(string(path); model=cpu(model), task=task)
end

"""
    loadtaskmodel(path) -> (task, model)

Load a trained `model` along with a `task` from `path` that were saved
using [`savetaskmodel`](#).

[JLD2.jl](https://github.com/JuliaIO/JLD2.jl) is used for serialization.
"""
function loadtaskmodel(path)
    isfile(path) || error("\"$path\" is not an existing file.")
    task, model = jldopen(string(path), "r") do f
        f["task"], f["model"]
    end
    return task, model
end


text:
module ProximalAlgorithms

using ProximalOperators
import ProximalOperators: prox!, gradient!

const RealOrComplex{R} = Union{R,Complex{R}}
const Maybe{T} = Union{T,Nothing}

"""
    prox!(y, f, x, gamma)

Compute the proximal mapping of `f` at `x`, with stepsize `gamma`, and store the result in `y`.
Return the value of `f` at `y`.
"""
prox!(y, f, x, gamma)

"""
    gradient!(g, f, x)

Compute the gradient of `f` at `x`, and stores it in `y`. Return the value of `f` at `x`.
"""
gradient!(y, f, x)

# TODO move out
ProximalOperators.is_quadratic(::Any) = false

# utilities

include("utilities/ad.jl")
include("utilities/conjugate.jl")
include("utilities/fb_tools.jl")
include("utilities/iteration_tools.jl")

# acceleration operators

include("accel/traits.jl")
include("accel/lbfgs.jl")
include("accel/anderson.jl")
include("accel/nesterov.jl")
include("accel/broyden.jl")
include("accel/noaccel.jl")

# algorithms

include("algorithms/forward_backward.jl")
include("algorithms/fast_forward_backward.jl")
include("algorithms/zerofpr.jl")
include("algorithms/panoc.jl")
include("algorithms/douglas_rachford.jl")
include("algorithms/drls.jl")
include("algorithms/primal_dual.jl")
include("algorithms/davis_yin.jl")
include("algorithms/li_lin.jl")
include("algorithms/fista.jl")
include("algorithms/panocplus.jl")
include("algorithms/pnp_douglas_rachford.jl")
include("algorithms/admm.jl")

end # module


text:
###################################
# Sim statistics aggregator
###################################

"""
node statistics
Creates a dictionary with keys as unique intersections ids and values as NodeStat struct containing coordinates of each intersection, number of agents driving by each intersection and their demographic profiles stored in data frame.

**Arguments**
* `sim_data` : `SimData` object used in sumulation

"""
function node_statistics(sim_data::OpenStreetMapXSim.SimData)::Dict{Int,OpenStreetMapXSim.NodeStat}
    nodes_stats = Dict{Int,OpenStreetMapXSim.NodeStat}()
    for (key,value) in sim_data.map_data.intersections
        coords = OpenStreetMapX.LLA(sim_data.map_data.nodes[key],OpenStreetMapX.center(sim_data.map_data.bounds))
        latitude, longitude = coords.lat, coords.lon
        nodes_stats[key] = OpenStreetMapXSim.NodeStat(0,latitude,longitude,nothing)
    end
    return nodes_stats
end

"""
Statistics aggregator

Aggregates data for each route intersection

**Arguments**
* `nodes_stats` : dictionary with keys as unique nodes_id and values as NodeStat struct
* `agent_profile` : selected agent profile as a DemoProfile object
* `route` : route represented by nodes ids

"""
function stats_aggregator!(nodes_stats::Dict{Int,OpenStreetMapXSim.NodeStat},
                        agent_profile::DataFrames.DataFrame,
                        route::Array{Int,1} )
    for indice in route
        node = nodes_stats[indice]
        node.count += 1
		if isa(node.agents_data,Nothing)
			node.agents_data = deepcopy(agent_profile)
		else
			append!(node.agents_data,agent_profile)
		end
    end
end


text:
## test DoublePendulum after creating a joint between two single pendulums
# include("../src/plotSol.jl")
include("../src/plotSol.jl")
include("../src/simulateDiff.jl")
include("../src/OrientationConversion.jl")
using Revise
using JLD
using BenchmarkTools
##
clearconsole()
j1 = Nothing
j2 = Nothing
## Pendulum 1
m = 1.0; l = 1.0; # Mass and length of bar
# Assuming bar revolves about Y axis
I1 = [1 0 0; 0 m*l^2/12 0; 0 0 1]

# Testing Dynamics with Revolute Joint
Pend1 = RigidBody(m,I1,2)
RbI = InertialFrameAsRB()
x0R1 = [l/2;zeros(2);[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Pend1,x0R1)

axisY = [0.0 1.0 0.0][:] # Axis about which bar is revolving

rj1 = [0.0 0.0 0.0][:] # Joint Location in body frame of first body
rj2 = [-l/2 0.0 0.0][:] # Joint location in body frame of second body

j1_1 = Joint(RbI,Pend1,rj1,rj2,type="Revolute",axis=axisY)
j1 = [j1_1];

# just simulate without any external force or gravity first
g = zeros(3);
tEnd = 1.0; tInt = 0.01;
tSim, solFinal = simulateDiff(tEnd,tInt,j1,g=g);
solPend1 = solFinal[1];

## Pendulum 2
Pend2 = RigidBody(m,I1,2);
x0R2 = [3*l/2;zeros(2);[1;zeros(3)];zeros(3);zeros(4)]
initialiseRigidBody!(Pend2,x0R2)

rj1 = [l 0.0 0.0][:] #Joint location wrt inertial frame
rj2 = [-l/2 0.0 0.0][:] #Joint location in body frame of pend2

j2_1 = Joint(RbI,Pend2,rj1,rj2,type="Revolute",axis=axisY);
j2 = [j2_1];

# just simulate without any external force or gravity first
g = zeros(3);
tEnd = 1.0; tInt = 0.01;
tSim, solFinal = simulateDiff(tEnd,tInt,j1,g = g);
solPend2 = solFinal[1];

## Join the two pendulums
println("########################")
println("Joining the two pendulums")
rj1 = [l/2;zeros(2)];
rj2 = [-l/2;zeros(2)];
joinBodies!(j1,j2,2,2,rj1,rj2, type = "Revolute", axis = axisY);

## Simulation
g = [0.0,0.0,-9.806]; # Gravity Vector.
tEnd = 1.0
tSpan = 0.01
tSim, solFinal = simulateDiff(tEnd,tSpan,j1,g=g)

solPend1  = solFinal[1]
solPend2  = solFinal[2]
##
# Check if joint location has moved
jointLoc1 = Matrix{Float64}(undef,length(tSim),3)
jointLoc2 = Matrix{Float64}(undef,length(tSim),3)
ωSol = Matrix{Float64}(undef,length(tSim),3)
for i=1:length(tSim)
    Pend1.dcm = quat2dcm(solPend1.β[i,:])
    Pend2.dcm = quat2dcm(solPend2.β[i,:])
    # ωSol[i,:] = angVel(solPend1.β[i,:],solPend1.βdot[i,:])
    jointLoc1[i,:] = solPend1.r[i,:] + transpose(Pend1.dcm)*rj1
    jointLoc2[i,:] = solPend2.r[i,:] + transpose(Pend2.dcm)*rj2
end
##
close("all");
# plotErrNorm(tSim,solPend.β)
plotPos(tSim,jointLoc1 .- jointLoc2)
# plotPos(tSim,jointLoc2)
# plotQuat(tSim, solPend1.β)
# plotAngVel(tSim, ωSol)
# plotPos(tSim, solPend.r)
## Compare with MinReal
solMin = load("DoublePendulumMin.jld")
r1Min = solMin["r1Min"]
r2Min = solMin["r2Min"]

# Plotting errors
plotPos(tSim,r1Min-solPend1.r)
plotPos(tSim,r2Min-solPend2.r)
# plotAngVel(tSim,ωMin+ωSol)
# plotPos(tSim,solPend.r)

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
#
# Correctness Tests
#

fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
quiet = length(ARGS) > 0 && ARGS[1] == "-q"
anyerrors = false

using InMemoryDatasets, Test, Random

if VERSION > v"1.3"
    if Threads.nthreads() < 2
        @warn("Running with only one thread: correctness of parallel operations is not tested")
    else
        @show Threads.nthreads()
    end
end

my_tests = ["constructors.jl",
            "characters.jl",
            "utils.jl",
            "subdataset.jl",
            "broadcasting.jl",
            "data.jl",
            "byrow.jl",
            "sort.jl",
            "grouping.jl",
            "join.jl",
            "transpose.jl"
            ]

println("Running tests:")

for my_test in my_tests
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch e
        global anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        if fatalerrors
            rethrow(e)
        elseif !quiet
            showerror(stdout, e, backtrace())
            println()
        end
    end
end

if anyerrors
    throw("Tests failed")
end


text:
export AcousticPropagatorParams, solve

struct AcousticPropagatorParams{DIM, AT<:AbstractArray{Float64,DIM}}
    # number of grids along x,y axis and time steps
    NX::Int
    NY::Int
    NSTEP::Int

    # size of grid cell and time step
    DELTAX::Float64
    DELTAY::Float64
    DELTAT::Float64

    # Auxilliary Data
    Σx::AT
    Σy::AT
end

function AcousticPropagatorParams(; nx::Int, ny::Int, nstep::Int,
        dx::Float64, dy::Float64, dt::Float64,
        Rcoef::Float64=0.001, # Relative reflection coefficient
        vp_ref::Float64=1000.0,
        npoints_PML::Int=12,
        USE_PML_XMAX::Bool = true,
        USE_PML_XMIN::Bool = true,
        USE_PML_YMAX::Bool = true,
        USE_PML_YMIN::Bool = true)
    # computing damping coefficient
    Lx = npoints_PML * dx
    Ly = npoints_PML * dy
    damping_x = vp_ref/Lx*log(1/Rcoef)
    damping_y = vp_ref/Ly*log(1/Rcoef)

    Σx, Σy = zeros(nx+2, ny+2), zeros(nx+2, ny+2)
    for i = 1:nx+2
        for j = 1:ny+2
            Σx[i,j] = pml_helper((i-1)*dx, nx, dx,
                damping_x, npoints_PML,
                USE_PML_XMIN, USE_PML_XMAX)
            Σy[i,j] = pml_helper((j-1)*dy, ny, dy,
                damping_y, npoints_PML,
                USE_PML_YMIN, USE_PML_YMAX)
        end
    end
    return AcousticPropagatorParams(nx, ny, nstep, dx, dy, dt, Σx, Σy)
end

function pml_helper(x::Float64, nx::Int, dx::Float64, ξx::Float64, npoints_PML::Int,
        USE_PML_XMIN, USE_PML_XMAX)
    Lx = npoints_PML * dx
    out = 0.0
    if x<Lx && USE_PML_XMIN
        d = abs(Lx-x)
        out = ξx * (d/Lx - sin(2π*d/Lx)/(2π))
    elseif x>dx*(nx+1)-Lx && USE_PML_XMAX
        d = abs(x-(dx*(nx+1)-Lx))
        out = ξx * (d/Lx - sin(2π*d/Lx)/(2π))
    end
    return out
end

function one_step!(param::AcousticPropagatorParams, u, w, wold, φ, ψ, σ, τ, c)
    Δt = param.DELTAT
    hx, hy = param.DELTAX, param.DELTAY

    @inbounds for j=2:param.NY+1, i=2:param.NX+1
        uij = (2 - σ[i,j]*τ[i,j]*Δt^2 - 2*Δt^2/hx^2 * c[i,j] - 2*Δt^2/hy^2 * c[i,j]) * w[i,j] +
            c[i,j] * (Δt/hx)^2  *  (w[i+1,j]+w[i-1,j]) +
            c[i,j] * (Δt/hy)^2  *  (w[i,j+1]+w[i,j-1]) +
            (Δt^2/(2hx))*(φ[i+1,j]-φ[i-1,j]) +
            (Δt^2/(2hy))*(ψ[i,j+1]-ψ[i,j-1]) -
            (1 - (σ[i,j]+τ[i,j])*Δt/2) * wold[i,j]
        u[i,j] = uij / (1 + (σ[i,j]+τ[i,j])/2*Δt)
    end
    @inbounds for j=2:param.NY+1, i=2:param.NX+1
        φ[i,j] = (1. -Δt*σ[i,j]) * φ[i,j] + Δt * c[i,j] * (τ[i,j] -σ[i,j])/2hx *
            (u[i+1,j]-u[i-1,j])
        ψ[i,j] = (1-Δt*τ[i,j]) * ψ[i,j] + Δt * c[i,j] * (σ[i,j] -τ[i,j])/2hy *
            (u[i,j+1]-u[i,j-1])
    end
end

function solve(param::AcousticPropagatorParams, src,
            srcv::Array{Float64, 1}, c::Array{Float64, 2})

    tu = zeros(param.NX+2, param.NY+2, param.NSTEP+1)
    tφ = zeros(param.NX+2, param.NY+2)
    tψ = zeros(param.NX+2, param.NY+2)

    for i = 3:param.NSTEP+1
        one_step!(param, view(tu,:,:,i), view(tu,:,:,i-1), view(tu,:,:,i-2), tφ, tψ, param.Σx, param.Σy, c)
        tu[src[1], src[2], i] += srcv[i-2]*param.DELTAT^2
    end
    tu
end

function solve_final(param::AcousticPropagatorParams, src,
            srcv::Array{Float64, 1}, c::Array{Float64, 2})

    tupre = zeros(param.NX+2, param.NY+2)
    tu = zeros(param.NX+2, param.NY+2)
    tφ = zeros(param.NX+2, param.NY+2)
    tψ = zeros(param.NX+2, param.NY+2)

    for i = 3:param.NSTEP+1
        tu_ = zeros(param.NX+2, param.NY+2)
        one_step!(param, tu_, tu, tupre, tφ, tψ, param.Σx, param.Σy, c)
        tu, tupre = tu_, tu
        tu[src...] += srcv[i-2]*param.DELTAT^2
    end
    tu
end


text:
import Base: kill
import Sockets: connect, listenany, accept, IPv4, getsockname, getaddrinfo, wait_connected

################################################################################
# MPI Cluster Manager
# Note: The cluster manager object lives only in the manager process,
# except for MPI_TRANSPORT_ALL

# There are three different transport modes:

# MPI_ON_WORKERS: Use MPI between the workers only, not for the manager. This
# allows interactive use from a Julia shell, using the familiar `addprocs`
# interface.

# MPI_TRANSPORT_ALL: Use MPI on all processes; there is no separate manager
# process. This corresponds to the "usual" way in which MPI is used in a
# headless mode, e.g. submitted as a script to a queueing system.

# TCP_TRANSPORT_ALL: Same as MPI_TRANSPORT_ALL, but Julia uses TCP for its
# communication between processes. MPI can still be used by the user.

@enum TransportMode MPI_ON_WORKERS MPI_TRANSPORT_ALL TCP_TRANSPORT_ALL

mutable struct MPIManager <: ClusterManager
    np::Int # number of worker processes (excluding the manager process)
    mpi2j::Dict{Int,Int} # map MPI ranks to Julia processes
    j2mpi::Dict{Int,Int} # map Julia to MPI ranks
    mode::TransportMode

    launched::Bool # Are the MPI processes running?
    mpirun_cmd::Cmd # How to start MPI processes
    launch_timeout::Int # seconds

    initialized::Bool # All workers registered with us
    cond_initialized::Condition # notify this when all workers registered

    # TCP Transport
    port::UInt16
    ip::UInt32
    stdout_ios::Array

    # MPI transport
    rank2streams::Dict{Int,Tuple{IO,IO}} # map MPI ranks to (input,output) streams
    ranks_left::Array{Int,1} # MPI ranks for which there is no stream pair yet

    # MPI_TRANSPORT_ALL
    comm::MPI.Comm
    initiate_shutdown::Channel{Nothing}
    sending_done::Channel{Nothing}
    receiving_done::Channel{Nothing}

    function MPIManager(; np::Integer = Sys.CPU_THREADS,
                          mpirun_cmd::Cmd = `$mpiexec -n $np`,
                          launch_timeout::Real = 60.0,
                          mode::TransportMode = MPI_ON_WORKERS,
                          master_tcp_interface::String="" )
        mgr = new()
        mgr.np = np
        mgr.mpi2j = Dict{Int,Int}()
        mgr.j2mpi = Dict{Int,Int}()
        mgr.mode = mode

        # Only start MPI processes for MPI_ON_WORKERS
        mgr.launched = mode != MPI_ON_WORKERS
        @assert MPI.Initialized() == mgr.launched
        mgr.mpirun_cmd = mpirun_cmd
        mgr.launch_timeout = launch_timeout

        mgr.initialized = false
        mgr.cond_initialized = Condition()
        if np == 0
            # Special case: no workers
            mgr.initialized = true
            if mgr.mode != MPI_ON_WORKERS
                # Set up mapping for the manager
                mgr.j2mpi[1] = 0
                mgr.mpi2j[0] = 1
            end
        end

        # Listen to TCP sockets if necessary
        if mode != MPI_TRANSPORT_ALL
            # Start a listener for capturing stdout from the workers
            if master_tcp_interface != ""
               # Listen on specified server interface
               # This allows direct connection from other hosts on same network as
               # specified interface.
               port, server =
                listenany(getaddrinfo(master_tcp_interface), 11000)
            else
               # Listen on default interface (localhost)
               # This precludes direct connection from other hosts.
               port, server = listenany(11000)
            end
            ip = getsockname(server)[1].host
            @async begin
                while true
                    sock = accept(server)
                    push!(mgr.stdout_ios, sock)
                end
            end
            mgr.port = port
            mgr.ip = ip
            mgr.stdout_ios = IO[]
        else
            mgr.rank2streams = Dict{Int,Tuple{IO,IO}}()
            size = MPI.Comm_size(MPI.COMM_WORLD)
            mgr.ranks_left = collect(1:size-1)
        end

        if mode == MPI_TRANSPORT_ALL
            mgr.sending_done = Channel{Nothing}(np)
            mgr.receiving_done = Channel{Nothing}(1)
        end
        mgr.initiate_shutdown = Channel{Nothing}(1)
        global initiate_shutdown = mgr.initiate_shutdown

        return mgr
    end
end

function Base.show(io::IO, mgr::MPIManager)
    print(io, "MPI.MPIManager(np=$(mgr.np),launched=$(mgr.launched),mode=$(mgr.mode))")
end

################################################################################
# Cluster Manager functionality required by Base, mostly targeting the
# MPI_ON_WORKERS case

# Launch a new worker, called from Base.addprocs
function Distributed.launch(mgr::MPIManager, params::Dict,
                instances::Array, cond::Condition)
    try
        if mgr.mode == MPI_ON_WORKERS
            # Start the workers
            if mgr.launched
                println("Reuse of an MPIManager is not allowed.")
                println("Try again with a different instance of MPIManager.")
                throw(ErrorException("Reuse of MPIManager is not allowed."))
            end
            cookie = string(":cookie_",Distributed.cluster_cookie())
            setup_cmds = `import MPIClusterManagers\;MPIClusterManagers.setup_worker'('$(mgr.ip),$(mgr.port),$cookie')'`
            mpi_cmd = `$(mgr.mpirun_cmd) $(params[:exename]) -e $(Base.shell_escape(setup_cmds))`
            open(detach(mpi_cmd))
            mgr.launched = true
        end

        if mgr.mode != MPI_TRANSPORT_ALL
            # Wait for the workers to connect back to the manager
            t0 = time()
            while (length(mgr.stdout_ios) < mgr.np &&
                   time() - t0 < mgr.launch_timeout)
                sleep(1.0)
            end
            if length(mgr.stdout_ios) != mgr.np
                error("Timeout -- the workers did not connect to the manager")
            end

            # Traverse all worker I/O streams and receive their MPI rank
            configs = Array{WorkerConfig}(undef, mgr.np)
            @sync begin
                for io in mgr.stdout_ios
                    @async let io=io
                        config = WorkerConfig()
                        config.io = io
                        # Add config to the correct slot so that MPI ranks and
                        # Julia pids are in the same order
                        rank = Serialization.deserialize(io)
                        idx = mgr.mode == MPI_ON_WORKERS ? rank+1 : rank
                        configs[idx] = config
                    end
                end
            end

            # Append our configs and notify the caller
            append!(instances, configs)
            notify(cond)
        else
            # This is a pure MPI configuration -- we don't need any bookkeeping
            for cnt in 1:mgr.np
                push!(instances, WorkerConfig())
            end
            notify(cond)
        end

    catch e
        println("Error in MPI launch $e")
        rethrow(e)
    end
end

# Entry point for MPI worker processes for MPI_ON_WORKERS and TCP_TRANSPORT_ALL
setup_worker(host, port) = setup_worker(host, port, nothing)
function setup_worker(host, port, cookie)
    !MPI.Initialized() && MPI.Init()
    # Connect to the manager
    io = connect(IPv4(host), port)
    wait_connected(io)
    redirect_stdout(io)
    redirect_stderr(io)

    # Send our MPI rank to the manager
    rank = MPI.Comm_rank(MPI.COMM_WORLD)
    Serialization.serialize(io, rank)

    # Hand over control to Base
    if cookie == nothing
        Distributed.start_worker(io)
    else
        if isa(cookie, Symbol)
            cookie = string(cookie)[8:end] # strip the leading "cookie_"
        end
        Distributed.start_worker(io, cookie)
    end
end

# Manage a worker (e.g. register / deregister it)
function Distributed.manage(mgr::MPIManager, id::Integer, config::WorkerConfig, op::Symbol)
    if op == :register
        # Retrieve MPI rank from worker
        # TODO: Why is this necessary? The workers already sent their rank.
        rank = remotecall_fetch(()->MPI.Comm_rank(MPI.COMM_WORLD), id)
        mgr.j2mpi[id] = rank
        mgr.mpi2j[rank] = id

        if length(mgr.j2mpi) == mgr.np
            # All workers registered
            mgr.initialized = true
            notify(mgr.cond_initialized)
            if mgr.mode != MPI_ON_WORKERS
                # Set up mapping for the manager
                mgr.j2mpi[1] = 0
                mgr.mpi2j[0] = 1
            end
        end
    elseif op == :deregister
        info("pid=$(getpid()) id=$id op=$op")
        # TODO: Sometimes -- very rarely -- Julia calls this `deregister`
        # function, and then outputs a warning such as """error in running
        # finalizer: ErrorException("no process with id 3 exists")""". These
        # warnings seem harmless; still, we should find out what is going wrong
        # here.
    elseif op == :interrupt
        # TODO: This should never happen if we rmprocs the workers properly
        info("pid=$(getpid()) id=$id op=$op")
        @assert false
    elseif op == :finalize
        # This is called from within a finalizer after deregistering; do nothing
    else
        info("pid=$(getpid()) id=$id op=$op")
        @assert false # Unsupported operation
    end
end

# Kill a worker
function kill(mgr::MPIManager, pid::Int, config::WorkerConfig)
    # Exit the worker to avoid EOF errors on the workers
    @spawnat pid begin
        MPI.Finalize()
        exit()
    end
    Distributed.set_worker_state(Distributed.Worker(pid), Distributed.W_TERMINATED)
end

# Set up a connection to a worker
function connect(mgr::MPIManager, pid::Int, config::WorkerConfig)
    if mgr.mode != MPI_TRANSPORT_ALL
        # Forward the call to the connect function in Base
        return invoke(connect, Tuple{ClusterManager, Int, WorkerConfig},
                      mgr, pid, config)
    end

    rank = MPI.Comm_rank(mgr.comm)
    if rank == 0
        # Choose a rank for this worker
        to_rank = pop!(mgr.ranks_left)
        config.connect_at = to_rank
        return start_send_event_loop(mgr, to_rank)
    else
        return start_send_event_loop(mgr, config.connect_at)
    end
end

# Event loop for sending data to one other process, for the MPI_TRANSPORT_ALL
# case
function start_send_event_loop(mgr::MPIManager, rank::Integer)
    try
        r_s = Base.BufferStream()
        w_s = Base.BufferStream()
        mgr.rank2streams[rank] = (r_s, w_s)

        # TODO: There is one task per communication partner -- this can be
        # quite expensive when there are many workers. Design something better.
        # For example, instead of maintaining two streams per worker, provide
        # only abstract functions to write to / read from these streams.
        @async begin
            rr = MPI.Comm_rank(mgr.comm)
            reqs = MPI.Request[]
            while !isready(mgr.initiate_shutdown)
                # When data are available, send them
                while bytesavailable(w_s) > 0
                    data = take!(w_s.buffer)
                    push!(reqs, MPI.Isend(data, rank, 0, mgr.comm))
                end
                if !isempty(reqs)
                    (indices, stats) = MPI.Testsome!(reqs)
                    filter!(req -> req != MPI.REQUEST_NULL, reqs)
                end
                # TODO: Need a better way to integrate with libuv's event loop
                yield()
            end
            put!(mgr.sending_done, nothing)
        end
        (r_s, w_s)
    catch e
        Base.show_backtrace(stdout, catch_backtrace())
        println(e)
        rethrow(e)
    end
end

################################################################################
# Alternative startup model: All Julia processes are started via an external
# mpirun, and the user does not call addprocs.

# Enter the MPI cluster manager's main loop (does not return on the workers)
function start_main_loop(mode::TransportMode=TCP_TRANSPORT_ALL;
                         comm::MPI.Comm=MPI.COMM_WORLD)
    !MPI.Initialized() && MPI.Init()
    @assert MPI.Initialized() && !MPI.Finalized()
    if mode == TCP_TRANSPORT_ALL
        # Base is handling the workers and their event loop
        # The workers have no manager object where to store the communicator.
        # TODO: Use a global variable?
        comm = MPI.COMM_WORLD
        rank = MPI.Comm_rank(comm)
        size = MPI.Comm_size(comm)
        if rank == 0
            # On the manager: Perform the usual steps
            # Create manager object
            mgr = MPIManager(np=size-1, mode=mode)
            mgr.comm = comm
            # Needed because of Julia commit https://github.com/JuliaLang/julia/commit/299300a409c35153a1fa235a05c3929726716600
            if isdefined(Distributed, :init_multi)
                Distributed.init_multi()
            end
            # Send connection information to all workers
            # TODO: Use Bcast
            for j in 1:size-1
                cookie = Distributed.cluster_cookie()
                MPI.send((mgr.ip, mgr.port, cookie), j, 0, comm)
            end
            # Tell Base about the workers
            addprocs(mgr)
            return mgr
        else
            # On a worker: Receive connection information
            (obj, status) = MPI.recv(0, 0, comm)
            (host, port, cookie) = obj
            # Call the regular worker entry point
            setup_worker(host, port, cookie) # does not return
        end
    elseif mode == MPI_TRANSPORT_ALL
        comm = MPI.Comm_dup(comm)
        rank = MPI.Comm_rank(comm)
        size = MPI.Comm_size(comm)
        # We are handling the workers and their event loops on our own
        if rank == 0
            # On the manager:
            # Create manager object
            mgr = MPIManager(np=size-1, mode=mode)
            mgr.comm = comm

            # Send the cookie over. Introduced in v"0.5.0-dev+4047". Irrelevant under MPI
            # transport, but need it to satisfy the changed protocol.
            Distributed.init_multi()
            MPI.bcast(Distributed.cluster_cookie(), 0, comm)
            # Start event loop for the workers
            @async receive_event_loop(mgr)
            # Tell Base about the workers
            addprocs(mgr)
            return mgr
        else
            # On a worker:
            # Create a "fake" manager object since Base wants one
            mgr = MPIManager(np=size-1, mode=mode)
            mgr.comm = comm
            # Recv the cookie
            cookie = MPI.bcast(nothing, 0, comm)
            Distributed.init_worker(cookie, mgr)
            # Start a worker event loop
            receive_event_loop(mgr)
            if isdefined(MPI, :free) && hasmethod(MPI.free, Tuple{MPI.Comm})
                MPI.free(comm)
            end
            MPI.Finalize()
            exit()
        end
    else
        error("Unknown mode $mode")
    end
end

# Event loop for receiving data, for the MPI_TRANSPORT_ALL case
function receive_event_loop(mgr::MPIManager)
    num_send_loops = 0
    while !isready(mgr.initiate_shutdown)
        (hasdata, stat) = MPI.Iprobe(MPI.MPI_ANY_SOURCE, 0, mgr.comm)
        if hasdata
            count = MPI.Get_count(stat, UInt8)
            buf = Array{UInt8}(undef, count)
            from_rank = MPI.Get_source(stat)
            MPI.Recv!(buf, from_rank, 0, mgr.comm)

            streams = get(mgr.rank2streams, from_rank, nothing)
            if streams == nothing
                # This is the first time we communicate with this rank.
                # Set up a new connection.
                (r_s, w_s) = start_send_event_loop(mgr, from_rank)
                Distributed.process_messages(r_s, w_s)
                num_send_loops += 1
            else
                (r_s, w_s) = streams
            end
            write(r_s, buf)
        else
            # TODO: Need a better way to integrate with libuv's event loop
            yield()
        end
    end

    for i in 1:num_send_loops
        fetch(mgr.sending_done)
    end
    put!(mgr.receiving_done, nothing)
end

# Stop the main loop
# This function should be called by the main process only.
function stop_main_loop(mgr::MPIManager)
    if mgr.mode == TCP_TRANSPORT_ALL
        # Shut down all workers
        rmprocs(workers())
        # Poor man's flush of the send queue
        sleep(1)
        put!(mgr.initiate_shutdown, nothing)
        MPI.Finalize()
    elseif mgr.mode == MPI_TRANSPORT_ALL
        # Shut down all workers, but not ourselves yet
        for i in workers()
            if i != myid()
                @spawnat i begin
                    global initiate_shutdown
                    put!(initiate_shutdown, nothing)
                end
            end
        end
        # Poor man's flush of the send queue
        sleep(1)
        # Shut down ourselves
        put!(mgr.initiate_shutdown, nothing)
        wait(mgr.receiving_done)
        MPI.Finalize()
    else
        @assert false
    end
end

################################################################################
# MPI-specific communication methods

# Execute a command on all MPI ranks
# This uses MPI as communication method even if @everywhere uses TCP
function mpi_do(mgr::MPIManager, expr)
    !mgr.initialized && wait(mgr.cond_initialized)
    jpids = keys(mgr.j2mpi)
    refs = Array{Any}(undef, length(jpids))
    for (i,p) in enumerate(Iterators.filter(x -> x != myid(), jpids))
        refs[i] = remotecall(expr, p)
    end
    # Execution on local process should be last, since it can block the main
    # event loop
    if myid() in jpids
        refs[end] = remotecall(expr, myid())
    end

    # Retrieve remote exceptions if any
    @sync begin
        for r in refs
            @async begin
                resp = remotecall_fetch(r.where, r) do rr
                    wrkr_result = rr[]
                    # Only return result if it is an exception, i.e. don't
                    # return a valid result of a worker computation. This is
                    # a mpi_do and not mpi_callfetch.
                    isa(wrkr_result, Exception) ? wrkr_result : nothing
                end
                isa(resp, Exception) && throw(resp)
            end
        end
    end
    nothing
end

macro mpi_do(mgr, expr)
    quote
        # Evaluate expression in Main module
        thunk = () -> (Core.eval(Main, $(Expr(:quote, expr))); nothing)
        mpi_do($(esc(mgr)), thunk)
    end
end

# All managed Julia processes
Distributed.procs(mgr::MPIManager) = sort(collect(keys(mgr.j2mpi)))

# All managed MPI ranks
mpiprocs(mgr::MPIManager) = sort(collect(keys(mgr.mpi2j)))


text:
module MetidaCu
    using LinearAlgebra, CUDA
    import MetidaNLopt, Metida
    import MetidaNLopt: reml_sweep_β_cuda#, cudata
    import Metida: LMM, AbstractLMMDataBlocks, rmat_base_inc!, zgz_base_inc!, logerror!

    struct LMMDataBlocks{T1, T2} <: AbstractLMMDataBlocks
        # Fixed effect matrix views
        xv::T1
        # Responce vector views
        yv::T2
        function LMMDataBlocks(xv::Matrix{T}, yv::Vector{T}, vcovblock::Vector) where T
            x = Vector{CuArray{T, 2}}(undef, length(vcovblock))
            y = Vector{CuArray{T, 1}}(undef, length(vcovblock))
            for i = 1:length(vcovblock)
                x[i] = CuArray(view(xv, vcovblock[i],:))
                y[i] = CuArray(view(yv, vcovblock[i]))
            end
            new{typeof(x), typeof(y)}(x, y)
        end
        function LMMDataBlocks(lmm)
            return LMMDataBlocks(lmm.data.xv, lmm.data.yv, lmm.covstr.vcovblock)
        end
    end

    function MetidaNLopt.cudata(lmm::LMM)
        LMMDataBlocks(lmm.data.xv, lmm.data.yv, lmm.covstr.vcovblock)
    end

    function MetidaNLopt.reml_sweep_β_cuda(lmm::LMM, θ::Vector{T}) where T
        reml_sweep_β_cuda(lmm,  MetidaNLopt.cudata(lmm), θ)
    end

    function MetidaNLopt.reml_sweep_β_cuda(lmm::LMM, data::AbstractLMMDataBlocks, θ::Vector{T}) where T
        n             = length(lmm.covstr.vcovblock)
        N             = length(lmm.data.yv)
        c             = (N - lmm.rankx)*log(2π)
        #-----------------------------------------------------------------------
        θ₁            = zero(T)
        θ₂            = zeros(T, lmm.rankx, lmm.rankx)
        θ₂tc          = CUDA.zeros(T, lmm.rankx, lmm.rankx)
        θ₃            = zero(T)
        βtc           = CUDA.zeros(T, lmm.rankx)
        β             = Vector{T}(undef, lmm.rankx)
        A             = Vector{CuArray{T, 2}}(undef, n)
        logdetθ₂      = zero(T)
        try
            @inbounds @simd for i = 1:n
                q    = length(lmm.covstr.vcovblock[i])
                qswm = q + lmm.rankx
                V    = zeros(T, q, q)
                Metida.zgz_base_inc!(V, θ, lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
                Metida.rmat_base_inc!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
                A[i] = CuArray(V)
            #-------------------------------------------------------------------
            # Cholesky
                A[i] = LinearAlgebra.LAPACK.potrf!('U', A[i])[1]
                θ₁  += sum(log.(diag(A[i])))*2
                vX   = LinearAlgebra.LAPACK.potrs!('U', A[i], copy(data.xv[i]))
                vy   = LinearAlgebra.LAPACK.potrs!('U', A[i], copy(data.yv[i]))
                CUDA.CUBLAS.gemm!('T', 'N', one(T), data.xv[i], vX, one(T), θ₂tc)
                CUDA.CUBLAS.gemv!('T', one(T), data.xv[i], vy, one(T), βtc)
            #-------------------------------------------------------------------
            end
        # Beta calculation
            copyto!(θ₂, θ₂tc)
            LinearAlgebra.LAPACK.potrf!('U', θ₂tc)
            copyto!(β, LinearAlgebra.LAPACK.potrs!('U', θ₂tc, βtc))
        # θ₃ calculation
            @inbounds @simd for i = 1:n
                r    = CUDA.CUBLAS.gemv!('N', -one(T), data.xv[i], βtc, one(T), copy(data.yv[i]))
                vr   = LinearAlgebra.LAPACK.potrs!('U', A[i], copy(r))
                θ₃  += r'*vr
            end
            logdetθ₂ = logdet(θ₂)
        catch e
            logerror!(e, lmm)
            return (Inf, nothing, nothing, nothing, false)
        end
        return   θ₁ + logdetθ₂ + θ₃ + c, β, θ₂, θ₃, true
    end
end # module


text:
module GLParallel

using GLVisualize
using GLFW
using GLWindow
using GLAbstraction
using ModernGL
using GeometryTypes
using Colors
using Reactive
using Quaternions
using FixedPointNumbers
using FileIO
using Packing
using SignedDistanceFields
using FreeType

# Parallel helper function to launch our glvisualize operations on the process,
# on which we run GLVisualize
# I will move these functions to GLVisualize once this is tested better!
function set_arg!(obj, key, value)
    @spawnat workerid[] begin
        GLAbstraction.set_arg!(fetch(obj), key, value)
        #Reactive.post_empty()
        nothing
    end
end

function _view(main, style=:default, cam=:perspective; kw_args...)
    @spawnat workerid[] begin
        obj = GLVisualize.visualize(main, style; kw_args...)
        GLVisualize._view(obj, camera=cam)
        obj
    end
end

function is_current_open()
    remotecall_fetch(workerid[]) do
        isempty(GLVisualize.get_screens()) && return false
        isopen(GLVisualize.current_screen())
    end
end

function empty_current!()
    remotecall(workerid[]) do
        empty!(GLVisualize.current_screen())
        nothing
    end
end

global const workerid = Ref(0)

function glscreen(
        id = workers()[end];
        color = RGBA(1f0, 1f0, 1f0, 1f0)
    )
    if workerid[] == 0 || id != workers()[end]
        workerid[] = id
    end
    # start GLVisualize worker process
    @spawnat workerid[] begin
        window = GLVisualize.glscreen(color = color)
        @async GLWindow.waiting_renderloop(window)
        nothing
    end
end

end


text:
# parse regex into an AST
function parse_regex(regex::AbstractString, datatype::DataType=String)
	child, pos = parse_regex_expression(regex, datatype)
	ast = ASTNode(:regex, [child,], Dict{Symbol, Any}(:datatype=>datatype))
end

# parse an expression within the regex
function parse_regex_expression(regex::AbstractString, datatype::DataType, pos=1)

	rootexpr = (pos==1)
	ornode = ASTNode(:or)
	andnode = ASTNode(:and)
	push!(ornode.children, andnode)

	escaped = false
	endofexpr = isempty(regex)

	while (!endofexpr)

		if pos > sizeof(regex)
			error("end of regex reached before parsing of expression complete")
		end

		chr = regex[pos]
		pos = nextind(regex, pos)

		if (chr=='\\') && (!escaped)
			escaped = true
		else

			if (chr=='|') && (!escaped)

				finish_and_node(andnode, datatype)
				andnode = ASTNode(:and)
				push!(ornode.children, andnode)

			elseif (chr=='(') && (!escaped)

				exprnode, pos = parse_regex_expression(regex, datatype, pos)
				push!(andnode.children, exprnode)

			elseif (chr==')') && (!escaped)
				if rootexpr
					error("additional right parenthesis encountered")
				end
				endofexpr = true

			elseif (chr in ['?']) && (!escaped)

				if isempty(andnode.children)
					error("no preceding item to apply $(chr) to")
				end

				precedingnode = pop!(andnode.children)

				optionalnode = ASTNode(:optional)
				push!(optionalnode.children, precedingnode)

				push!(andnode.children, optionalnode)

			elseif (chr in ['+', '*', '{']) && (!escaped)

				if isempty(andnode.children)
					error("no preceding item to apply $(chr) to")
				end

				precedingnode = pop!(andnode.children)

				quantifiernode = ASTNode(:quantifier)
				push!(quantifiernode.children, precedingnode)

				if chr=='{'

					quantifiernode.args[:lowerbound], pos = parse_regex_int(regex::AbstractString, pos)
					if (pos<=sizeof(regex)) && (regex[pos]==',')
						pos = nextind(regex,pos)
						if pos<=sizeof(regex)
							if regex[pos]=='}'
								quantifiernode.args[:upperbound] = typemax(Int)
							else
								quantifiernode.args[:upperbound], pos = parse_regex_int(regex::AbstractString, pos)
							end
						end
					else
						quantifiernode.args[:upperbound] = quantifiernode.args[:lowerbound]
					end
					if (pos>sizeof(regex)) || (regex[pos]!='}')
						error("end of regex reached before parsing of quantifier expression complete")
					else
						pos = nextind(regex,pos)
					end

				else

					quantifiernode.args[:lowerbound] = (chr=='+') ? 1 : 0
					quantifiernode.args[:upperbound] = typemax(Int)

				end

				push!(andnode.children, quantifiernode)

			elseif (chr=='[') && (!escaped)

				bracketnode, pos = parse_regex_bracket(regex, pos)
				push!(andnode.children, bracketnode)

			elseif ((chr in ['s', 'S', 'd', 'D', 'w', 'W', 'i', 'I', 'c', 'C']) && (escaped)) ||
							(chr in ['.'] && (!escaped))

				# TODO implicitly here, set of tokenspace is 9:10,13:13,32:126 - tab, new line, carriage return, plus printable ASCII
				# TODO to what extent should ASCII control characters be included? - here have tab (9), new line (10) and carriage return (13)
				# TODO strictly whitespace can include vertical tab (11) and form feed (13), but XML definition excludes them
				# TODO strictly wildcard (.) excludes just new line; but XML also excludes carriage return

				# TODO need to handle for non-ASCII datatypes

				if chr == 's'  # whitespace: space, tab, carriage return, new line
					classspace = (UnitRange)[9:10,13:13,32:32]
	      elseif chr == 'S'  # non-whitespace
	        classspace = (UnitRange)[33:126]
				elseif chr == 'd'  # digits (currently just ASCII digits)
		      classspace = (UnitRange)[48:57]
		    elseif chr == 'D'  # non-digits
					classspace = (UnitRange)[9:10,13:13,32:47,58:126]
				elseif chr == 'w'  # word characters: letters, digits, plus underscore
		      classspace = (UnitRange)[48:57,65:90,95:95,97:122]
		    elseif chr == 'W'  # non-word characters
					classspace = (UnitRange)[9:10,13:13,32:47,58:64,91:94,96:96,123:126]
				elseif chr == 'i'  # XSD extension - initial name characters: letters, plus hyphen
		      classspace = (UnitRange)[45:45,65:90,97:122]
		    elseif chr == 'I'  # XSD extension - not initial name characters
					classspace = (UnitRange)[9:10,13:13,32:44,46:64,91:96,123:126]
				elseif chr == 'c'  # XSD extension - name characters: letters, digits, plus hyphen, period, colon
		      classspace = (UnitRange)[45:46,48:58,65:90,97:122]
		    elseif chr == 'C'  # XSD extension - not name characters
					classspace = (UnitRange)[9:10,13:13,32:44,47:47,59:64,91:96,123:126]
	      else
				if chr != '.'
					warn("cannot process character class $(node.args[:value]) - using wildcard instead")
				end
	        	classspace = (UnitRange)[9:9,32:126] # wildcard (.) - excludes new line and -- for XML -- carriage return
	      end

				bracketnode = ASTNode(:bracket)
	      bracketnode.func = :bracket

	      for range in classspace
	        rangenode = ASTNode(:range)
	        rangenode.args[:value] = range
	        push!(bracketnode.children, rangenode)
				end

				push!(andnode.children, bracketnode)

			elseif (chr in ['^','$'] && (!escaped))

				# ignore: ^ and $ have no meaning in generation context TODO?

			else

				terminalnode = ASTNode(:terminal)
				terminalnode.args[:value] = string(chr)
				push!(andnode.children, terminalnode)

			end

			escaped = false

		end

		if rootexpr && (pos > sizeof(regex))
			endofexpr = true
		end

	end

	finish_and_node(andnode, datatype)

	ornode, pos

end

# at end of an 'and' node in the AST
function finish_and_node(andnode::ASTNode, datatype::DataType)

	if isempty(andnode.children)

		# special case of empty expression: we intepret this as any string of any characters from the tokenspace
		# for XML Schema, this is used for xs:string datatype
		# it differs from the regex ".*" in that all characters (including newlines) may be included

		bracketnode = ASTNode(:bracket)
	    bracketnode.func = :bracket

		# if datatype in [ASCIIString,]
  #       	classspace = (UnitRange)[9:10, 13:13, 32:126]
  #       elseif datatype in [UTF8String, UTF16String, UTF32String,]
			classspace = (UnitRange)[9:10, 13:13, 32:55295, 65536:131071, 131072:196607]
			# 55295 (0xD7FF) is last character before UTF-16 surrogates in plane 0 (BMP)
			# 65536-131071 (0x10000 - 0x1FFFF) is plane 1 (SMP)
			# 131072-196607 (0x20000 - 0x2FFFF) is plane 2 (SIP)
		# else
		# 	error("Do not know how to handle empty regex for string datatype $(datatype)")
		# end

	    for range in classspace # currently just ASCII
	      rangenode = ASTNode(:range)
	      rangenode.args[:value] = range
	      push!(bracketnode.children, rangenode)
		end

		quantifiernode = ASTNode(:quantifier)
		quantifiernode.args[:lowerbound] = 0
		quantifiernode.args[:upperbound] = typemax(Int)
		push!(quantifiernode.children, bracketnode)

		push!(andnode.children, quantifiernode)

	end

end

# parse integer value in the regex starting at pos
function parse_regex_int(regex::AbstractString, pos)

	startpos = pos
	while (pos<=sizeof(regex)) && isdigit(regex[pos])
		pos = nextind(regex, pos)
	end

	if pos > sizeof(regex)
		error("end of regex reached before parsing integer in regex qualifier")
	end

	endpos = prevind(regex, pos)
	if startpos > endpos
			error("empty string when parsing integer in regex qualifier")
	end

	parse(Int,regex[startpos:endpos]), pos

end

# parse a bracket expression in the regex
function parse_regex_bracket(regex::AbstractString, pos)

	bracketnode = ASTNode(:bracket)

	firstchar = true
	endofexpr = false

	while (!endofexpr)

		if pos > sizeof(regex)
			error("end of regex reached before parsing of bracket expression complete")
		end

		chr = regex[pos]
		pos = nextind(regex, pos)

		if (chr=='^') && (firstchar)

       process_error("cannot process negated bracket expressions")

		elseif (chr==']') && (!firstchar)

			endofexpr = true

		else

			if (nextind(regex,pos)<=sizeof(regex)) && (regex[pos]=='-') && (regex[nextind(regex,pos)]!=']')

				rangestart = convert(Int,chr)
				pos = nextind(regex,pos)
				rangeend = convert(Int,regex[pos])
				if rangeend < rangestart
					error("end of range in bracket expression is before the start of the range in $(regex)")
				end
				pos = nextind(regex,pos)

			else

				rangestart = convert(Int,chr)
				rangeend = rangestart

			end

			rangenode = ASTNode(:range)
			rangenode.args[:value] = rangestart:rangeend
			push!(bracketnode.children, rangenode)

		end

		firstchar = false

	end

	bracketnode, pos

end


text:
function get_cartesian_mesh_size(dim, order, prob_size)
    dims = zeros(CeedInt, dim)
    # Use the approximate formula:
    #    prob_size ~ num_elem * order^dim
    num_elem = div(prob_size, order^dim)
    s = 0 # find s: num_elem/2 < 2^s <= num_elem
    while num_elem > 1
        num_elem = div(num_elem, 2)
        s += 1
    end
    r = s%dim
    for d = 1:dim
        sd = div(s, dim)
        if r > 0
            sd += 1
            r -= 1
        end
        dims[d] = 2^sd
    end
    dims
end

struct FormRestrictionMode{T} end
const RestrictionOnly = FormRestrictionMode{:restr}()
const StridedOnly = FormRestrictionMode{:restr_i}()
const RestrictionAndStrided = FormRestrictionMode{:both}()

function build_cartesian_restriction(
    c::Ceed,
    dim,
    nxyz,
    order,
    ncomp,
    num_qpts;
    mode::Mode=RestrictionOnly,
) where {Mode}
    p::CeedInt = order
    pp1::CeedInt = p + 1
    nnodes::CeedInt = pp1^dim # number of scal. nodes per element
    elem_qpts::CeedInt = num_qpts^dim # number of qpts per element

    nd = CeedInt.(p*nxyz .+ 1)
    num_elem::CeedInt = prod(nxyz)
    scalar_size::CeedInt = prod(nd)
    size::CeedInt = scalar_size*ncomp

    form_restr = (Mode() != StridedOnly)
    form_strided = (Mode() != RestrictionOnly)

    # elem:         0             1                 n-1
    #        |---*-...-*---|---*-...-*---|- ... -|--...--|
    # nnodes:   0   1    p-1  p  p+1       2*p             n*p

    if form_restr
        el_nodes = zeros(CeedInt, num_elem*nnodes)
        exyz = zeros(CeedInt, dim)
        @inbounds @simd for e = 0:(num_elem-1)
            re::CeedInt = e
            for d::CeedInt = 1:dim
                exyz[d] = re%nxyz[d]
                re ÷= nxyz[d]
            end
            for lnodes::CeedInt = 0:(nnodes-1)
                gnodes::CeedInt = 0
                gnodes_stride::CeedInt = 1
                rnodes::CeedInt = lnodes
                for d::CeedInt = 1:dim
                    q = rnodes÷pp1
                    r = rnodes - pp1*q
                    gnodes::CeedInt += (exyz[d]*p + r)*gnodes_stride
                    gnodes_stride::CeedInt *= nd[d]
                    rnodes = q
                end
                el_nodes[e*nnodes+lnodes+1] = gnodes
            end
        end
    end

    restr =
        form_restr ?
        create_elem_restriction(
            c,
            num_elem,
            nnodes,
            ncomp,
            scalar_size,
            ncomp*scalar_size,
            el_nodes,
        ) : nothing
    restr_i =
        form_strided ?
        create_elem_restriction_strided(
            c,
            num_elem,
            elem_qpts,
            ncomp,
            ncomp*elem_qpts*num_elem,
            STRIDES_BACKEND,
        ) : nothing

    return size, restr, restr_i
end

function set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords)
    p = mesh_order
    nd = p*nxyz .+ 1
    num_elem::CeedInt = prod(nxyz)
    scalar_size::CeedInt = prod(nd)

    # The H1 basis uses Lobatto quadrature points as nodes.
    nodes::Vector{CeedScalar} = lobatto_quadrature(p + 1) # nodes are in [-1,1]
    nodes = 0.5 .+ 0.5*nodes

    @witharray coords = mesh_coords begin
        @inbounds @simd for gsnodes = 0:scalar_size-1
            rnodes = gsnodes
            for d = 1:dim
                ndd = nd[d]
                q = rnodes÷ndd
                r = rnodes - ndd*q

                d1d = r
                coords[gsnodes+scalar_size*(d-1)+1] = ((d1d÷p) + nodes[d1d%p+1])/nxyz[d]
                rnodes = q
            end
        end
    end
end

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
# The target is to develop an efficient MsFEM for quadrilateral mesh, for research purposes. Here, the quadrilaternal mesh is axis parallel, so it is easy to describe its elements, boundary and neighboring geometry

using SparseArrays
using LinearAlgebra
using Logging
import Base.Threads: nthreads, @threads

# 2d 2scale mesh
# element easily described by x and y coordinates [Nce*Nce elements]
# element boundary easily described by a vector of boundary location (ordering: x increases first y second) [write a function/attribute]

### simple boundary treatment: just run over interior and boundaries separately (only for the boundary part use the if-else to select the correct b.c.). Each time run over it, construct the RHS and boundary data on the fly (not costly, as only once)

# basis function change when it is near to the boudary
# write a MsFEM that can support subsequent solve for the SAME boundary type!
# the key is the get a fine scale solver (surrogate model) for general problem; provide a coarse/fine scale solver for multiscale problems

struct MsFEM_2d2ScaleUnifQuadMesh{Ti,Tf}
    Nce::Ti # number of coarse elements in each dimension
    Nfe::Ti # number of fine elements in each coarse element
    Ne::Ti # number of total fine elements in each dimension
    CGrid_x::Vector{Tf} # coarse uniform grid in x axis, boundary included
    CGrid_y::Vector{Tf} # coarse uniform grid in y axis, boundary included
    ElemNode_loc2glo::Function # local index (node of an element) to global index (global indexed node); one runs over x first then y
    # all but the first attribute are automatically generated
    LocalBdyIndice::Vector{Ti}
    LocalBdyCondition::Matrix{Tf}
end

struct MsFEM_store{Ti,Tf}
    BasisFuns::Array{Tf, 4} # basisfuns in each element
    Fine_localAs::Matrix{SparseMatrixCSC{Tf,Ti}}
    Fine_localMs::Matrix{SparseMatrixCSC{Tf,Ti}}
    A::SparseMatrixCSC{Tf,Ti}
end

# constructor
function MsFEM_2d2ScaleUnifQuadMesh(Nce, Nfe)
    x = collect(LinRange(0, 1, Nce+1))
    y = copy(x)
    function ElemNode_loc2glo(N, i, j, ind_node)
        # N: number of elements in one direction
        # i,j is the location of the element
        # edge marks the index of the edge (from 1 to 4, from bottom to top)
        if ind_node <= 2 # bottom
            global_idx = (N+1)*(j-1) + i + ind_node - 1;
        else # top
            global_idx = (N+1)*j + i + 4 - ind_node;
        end
        return global_idx
    end

    # local indices for local element
    LocalBdyIndice = reduce(vcat,collect.(
        [1:Nfe+1,
        Nfe+2:Nfe+1:(Nfe+1)*(Nfe+1),
        2*(Nfe+1):Nfe+1:(Nfe+1)*(Nfe+1),
        Nfe*Nfe+Nfe+2:(Nfe+1)*(Nfe+1)-1]
        )
    )

    # linear basis functions on the boundary (4 basis functions, equal 1 at one node and 0 at other nodes)
    LocalBdyCondition = reduce(vcat, collect.(
        [LinRange(1, 0, Nfe+1), LinRange(1-1/Nfe, 0, Nfe), zeros(Nfe), zeros(Nfe-1),
        LinRange(0, 1, Nfe+1), zeros(Nfe), LinRange(1-1/Nfe, 0, Nfe), zeros(Nfe-1),
        zeros(Nfe+1), zeros(Nfe), LinRange(1/Nfe, 1, Nfe), LinRange(1/Nfe, 1-1/Nfe, Nfe-1),
        zeros(Nfe+1), LinRange(1/Nfe, 1, Nfe), zeros(Nfe), LinRange(1-1/Nfe, 1/Nfe, Nfe-1)]
        )
    )
    LocalBdyCondition = reshape(LocalBdyCondition, 4*Nfe, 4)

    @info "[Mesh generation] mesh generated, $(Nce+1) coarse nodes a in each dimension"
    @info "[Mesh generation] in each coarse element, $(Nfe+1) fine nodes; in total $(Nce*Nfe+1) nodes in each dimension"

    return MsFEM_2d2ScaleUnifQuadMesh(Nce,Nfe,Nce*Nfe,x,y,ElemNode_loc2glo, LocalBdyIndice, LocalBdyCondition)
end

function MsFEM_StiffnMassAssembly(MsFEMparam, PDEparam)
    # Nce num of coarse elements each dimension
    # Nfe num of fine elements in each coarse element
    Nce = MsFEMparam.Nce

    # sparse assembling
    Irow = zeros(16*(Nce)^2)
    Jcol = copy(Irow)
    Aval = copy(Irow)
    count = 4;
    LocalBasisFuns = zeros((Nfe+1)^2,count,Nce,Nce);
    Fine_localAs = Matrix{SparseMatrixCSC{Float64,Int}}(undef,Nce,Nce)
    Fine_localMs = copy(Fine_localAs)
    # there might be inner patches, patches on the edge, and nodal patches
    # run over coarse patches
    println("[multithreading] using ", Threads.nthreads(), " threads")
    @threads for ci = 1:Nce
        for cj = 1:Nce
            coarse_localbasis, coarse_localA, fine_localA, fine_localM = MsFEM_LocalBasis(MsFEMparam, PDEparam, ci, cj);
            LocalBasisFuns[:,:,ci,cj] = coarse_localbasis; # store basis functions
            Fine_localAs[ci,cj] = fine_localA
            Fine_localMs[ci,cj] = fine_localM

            for p = 1:count
                global_p = MsFEMparam.ElemNode_loc2glo(Nce, ci, cj, p);
                for q = 1:count
                    global_q = MsFEMparam.ElemNode_loc2glo(Nce, ci, cj, q);
                    index = 16*(Nce)*(ci-1)+16*(cj-1)+(4)*(p-1)+q;
                    Irow[index] = global_p;
                    Jcol[index] = global_q;
                    Aval[index] = coarse_localA[p, q];
                end
            end
        end
    end
    A = sparse(Irow,Jcol,Aval,(Nce+1)^2,(Nce+1)^2)
    @info "[Assembly] finish assembly of stiffness matrice"

    return MsFEM_store(LocalBasisFuns,Fine_localAs,Fine_localMs,A)
end

function MsFEM_LocalBasis(MsFEMparam, PDEparam, ci, cj)

    Nfe = MsFEMparam.Nfe

    # boundary location
    b = MsFEMparam.LocalBdyIndice
    # b = reduce(vcat,collect.(
    #     [1:Nfe+1,
    #     Nfe+2:Nfe+1:(Nfe+1)*(Nfe+1),
    #     2*(Nfe+1):Nfe+1:(Nfe+1)*(Nfe+1),
    #     Nfe*Nfe+Nfe+2:(Nfe+1)*(Nfe+1)-1]
    #     )
    # )

    # linear boundary conditions
    # vec_bd_f = reduce(vcat, collect.(
    #     [LinRange(1, 0, Nfe+1), LinRange(1-1/Nfe, 0, Nfe), zeros(Nfe), zeros(Nfe-1),
    #     LinRange(0, 1, Nfe+1), zeros(Nfe), LinRange(1-1/Nfe, 0, Nfe), zeros(Nfe-1),
    #     zeros(Nfe+1), zeros(Nfe), LinRange(1/Nfe, 1, Nfe), LinRange(1/Nfe, 1-1/Nfe, Nfe-1),
    #     zeros(Nfe+1), LinRange(1/Nfe, 1, Nfe), zeros(Nfe), LinRange(1-1/Nfe, 1/Nfe, Nfe-1)]
    #     )
    # )
    # bd_f = reshape(vec_bd_f, 4*Nfe, 4) # 4 edges
    bd_f = MsFEMparam.LocalBdyCondition

    fine_localA, fine_localM = MsFEM_LocalHarmExt(MsFEMparam, PDEparam, ci, cj); # harmonic extension matrix
    Diri_A = copy(fine_localA); # local inner product A
    Diri_A[b,:] .= 0;
    F = zeros((Nfe+1)^2, 4)
    F[b,:] = bd_f
    Diri_A[b,b] .= sparse(I, length(b),length(b));

    coarse_localbasis = Diri_A\F
    coarse_localA = coarse_localbasis'*fine_localA*coarse_localbasis; # energy inner product

    return coarse_localbasis, coarse_localA, fine_localA, fine_localM
end

function MsFEM_LocalHarmExt(MsFEMparam, PDEparam, ci, cj)

    Nfe = MsFEMparam.Nfe

    xlow = MsFEMparam.CGrid_x[ci];
    xhigh = MsFEMparam.CGrid_x[ci+1];
    ylow = MsFEMparam.CGrid_y[cj];
    yhigh = MsFEMparam.CGrid_y[cj+1];
    x = collect(LinRange(xlow, xhigh, Nfe+1))
    y = collect(LinRange(ylow, yhigh, Nfe+1))

    # sparse assembling
    Irow = zeros(16*Nfe^2);
    Jcol = copy(Irow)
    Aval = copy(Irow)
    Mval = copy(Irow)

    for fi = 1:Nfe
        for fj = 1:Nfe
            local_A, local_M = MsFEM_LocalHarmExt_MatrixAssemby(MsFEMparam, PDEparam, x, y, fi, fj);
            for p = 1:4
                global_p = MsFEMparam.ElemNode_loc2glo(Nfe, fi, fj, p);
                for q = 1:4
                    index=16*Nfe*(fi-1)+16*(fj-1)+4*(p-1)+q;
                    global_q = MsFEMparam.ElemNode_loc2glo(Nfe, fi, fj, q);
                    Irow[index] = global_p
                    Jcol[index] = global_q
                    Aval[index] = local_A[p, q]
                    Mval[index] = local_M[p, q]
                end
            end
        end
    end
    A = sparse(Irow,Jcol,Aval,(Nfe+1)^2,(Nfe+1)^2)
    M = sparse(Irow,Jcol,Mval,(Nfe+1)^2,(Nfe+1)^2)
    return A, M
end

function MsFEM_LocalHarmExt_MatrixAssemby(MsFEMparam, PDEparam, x, y, i, j)

    Ne = MsFEMparam.Ne
    local_A = zeros(4,4)
    local_M = copy(local_A)

    xlow = x[i];
    xhigh = x[i+1];
    ylow = y[j];
    yhigh = y[j+1];
    xmid = (xlow+xhigh)/2;
    ymid = (ylow+yhigh)/2;

    for i = 1:4
        for j = 1:4
            if i==j
                local_A[i,j] = 2/3*PDEparam.a(xmid,ymid);
                local_M[i,j] = 1/9/Ne^2 # exact
            elseif i==j+2 || i==j-2
                local_A[i,j] = -1/3*PDEparam.a(xmid,ymid);
                local_M[i,j] = 1/36/Ne^2
            else
                local_A[i,j] = -1/6*PDEparam.a(xmid,ymid);
                local_M[i,j] = 1/18/Ne^2
            end
        end
    end
    return local_A, local_M
end

function MsFEM_BdyRhsAssembly(MsFEMparam, PDEparam, MsFEMstore)
    A = copy(MsFEMstore.A)

    BasisFuns = MsFEMstore.BasisFuns
    Fine_localAs = MsFEMstore.Fine_localAs
    Fine_localMs = MsFEMstore.Fine_localMs

    # bubble part
    Ne = MsFEMparam.Ne
    sol_bubble = zeros(Ne+1,Ne+1)

    ### right hand side
    Nce = MsFEMparam.Nce
    F = zeros((Nce+1)^2)
    count = 4

    # boundary location for bubble part
    b = MsFEMparam.LocalBdyIndice
    # b = reduce(vcat,collect.(
    #             [1:Nfe+1,
    #             Nfe+2:Nfe+1:(Nfe+1)*(Nfe+1),
    #             2*(Nfe+1):Nfe+1:(Nfe+1)*(Nfe+1),
    #             Nfe*Nfe+Nfe+2:(Nfe+1)*(Nfe+1)-1]
    #     )
    # )

    @threads for ci = 1:Nce
        for cj = 1:Nce
            xlow = MsFEMparam.CGrid_x[ci];
            xhigh = MsFEMparam.CGrid_x[ci+1];
            ylow = MsFEMparam.CGrid_y[cj];
            yhigh = MsFEMparam.CGrid_y[cj+1];
            x = collect(LinRange(xlow, xhigh, Nfe+1))
            y = collect(LinRange(ylow, yhigh, Nfe+1))
            f = [PDEparam.rhs(x[i],y[j]) for j in 1:Nfe+1 for i in 1:Nfe+1]

            # assembly global rhs
            @views val = f'*Fine_localMs[ci,cj]*BasisFuns[:,:,ci,cj]
            for p = 1:count
                global_p = MsFEMparam.ElemNode_loc2glo(Nce, ci, cj, p);
                F[global_p] += val[p]
            end

            # assembly local rhs for bubble part
            Diri_A = copy(Fine_localAs[ci, cj]);
            Diri_A[b,:] .= 0;
            Diri_A[b,b] .= sparse(I, length(b),length(b));
            local_F = Fine_localMs[ci,cj]*f
            local_F[b] .= 0

            @views sol_bubble[(ci-1)*Nfe+1:(ci)*Nfe+1, (cj-1)*Nfe+1:(cj)*Nfe+1] = Diri_A\(local_F)
        end
    end

    ### global boundary part
    # location
    b = reduce(vcat,collect.([
        1:Nce+1,Nce+2:Nce+1:(Nce+1)*(Nce+1),2*(Nce+1):Nce+1:(Nce+1)*(Nce+1),
        Nce*Nce+Nce+2:(Nce+1)*(Nce+1)-1,(Nce+1)^2+1:(Nce+1)^2
        ]));

    # assemby
    A[b,:] .= 0;

    # for Dirichlet boundary condition
    F[1:Nce+1] .= PDEparam.bdy_Diri.(MsFEMparam.CGrid_x,0.0)
    F[Nce+2:Nce+1:(Nce+1)*(Nce+1)] .= PDEparam.bdy_Diri.(0.0, MsFEMparam.CGrid_y[2:end])
    F[2*(Nce+1):Nce+1:(Nce+1)*(Nce+1)] .= PDEparam.bdy_Diri.(1.0, MsFEMparam.CGrid_y[2:end])
    F[Nce*Nce+Nce+2:(Nce+1)*(Nce+1)-1] .= PDEparam.bdy_Diri.(MsFEMparam.CGrid_x[2:end-1],1.0)

    A[b,b] .= sparse(I, length(b),length(b));
    return A, F, sol_bubble

end

function MsFEM_CoarseSolver(MsFEMparam,PDEparam)
    MsFEMstore = MsFEM_StiffnMassAssembly(MsFEMparam,PDEparam)
    A, F, sol_bubble = MsFEM_BdyRhsAssembly(MsFEMparam,PDEparam,MsFEMstore)
    coarse_sol = A\F
    @info "[Linear solver] linear system solved, coarse sol obtained"

    return coarse_sol, sol_bubble, MsFEMstore
end

function MsFEM_FineConstruct(coarse_sol, sol_bubble, MsFEMstore)
    Nce = MsFEMparam.Nce
    Nfe = MsFEMparam.Nfe
    BasisFuns = MsFEMstore.BasisFuns
    fine_sol = zeros(Nce*Nfe+1,Nce*Nfe+1)
    count = 4

    @threads for ci in 1:Nce
        for cj in 1:Nce
            val = zeros(Nfe+1,Nfe+1)

            for p = 1:count
                global_p = MsFEMparam.ElemNode_loc2glo(Nce, ci, cj, p);
                val += coarse_sol[global_p]*reshape(BasisFuns[:,p,ci,cj],Nfe+1,Nfe+1)
            end
            # for p = 1:count
            #     global_p = MsFEMparam.ElemNode_loc2glo(Nce, ci, cj, p);
            #     fine_sol[(ci-1)*Nfe+1:(ci)*Nfe+1, (cj-1)*Nfe+1:(cj)*Nfe+1] += coarse_sol[global_p]*reshape(BasisFuns[:,p,ci,cj],Nfe+1,Nfe+1)
            # end
            fine_sol[(ci-1)*Nfe+1:(ci)*Nfe+1, (cj-1)*Nfe+1:(cj)*Nfe+1] = val
        end
    end

    # repeated counting
    # edges
    # fine_sol[1:end,Nfe+1:Nfe:end-Nfe] /= 2
    # fine_sol[Nfe+1:Nfe:end-Nfe,1:end] /= 2

    @info "[Reconstruction] fine scale solution reconstructed"
    return reshape(fine_sol+sol_bubble,(Nce*Nfe+1)^2)
    # return reshape(fine_sol,(Nce*Nfe+1)^2)
end

function MsFEM_SubsequentSolve(MsFEMparam,PDEparam,MsFEMstore)
    A, F, sol_bubble = MsFEM_BdyRhsAssembly(MsFEMparam,PDEparam,MsFEMstore)
    coarse_sol = A\F
    @info "[Linear solver] linear system solved, coarse sol obtained"
    return coarse_sol, sol_bubble, MsFEMstore
end


text:
using QuadGK, HCubature

"""
Create a new directory `dirname` in the current path.
Append a number to the name if the directory already exists.
"""
function make_new_dir(dirname::String)
    if ispath(dirname)
        k=1
        while ispath(dirname*"_$k")
            k+=1
        end
        return mkdir(dirname*"_$k")
    else
        return mkdir(dirname)
    end
end

"""
    nice_number(x,K::Int)

Return a string which is the rounded value of `x` with `K` (decimal) digits and
if necessary, pad `"0"` to the right so that there are exactly `K` decimal places.
"""
function nice_float2string(x,K::Int)
    integer_part = trunc(Int,x)
    y = x - integer_part
    decimal_part = rpad(round(Int,y*10^K),K,"0")
    return "$(integer_part).$decimal_part"
end

"""
Return the coefficients `c1,c2,λ` of the SOH model computed for the `model` `"Fokker-Planck"` (default)
or `"BGK"` with the concentration parameters `κ`.
These coefficients are given by

c₁ = ∫\\_[0,π] cos(θ)exp(κcos(θ))dθ / ∫\\_[0,π] exp(κcos(θ))dθ

c₂ = ∫\\_[0,π] cos(θ)sin(θ)g(θ)exp(κcos(θ))dθ / ∫\\_[0,π] sin(θ)g(θ)exp(κcos(θ))dθ

λ = 1/κ

where

g(θ) = θ/κ - π/κ * ∫\\_[0,θ] exp(-κcos(ψ))dψ / ∫\\_[0,π] exp(-κcos(ψ))dψ in the Fokker-Planck case

and

g(θ) = 1 in the BGK case.
"""
function coefficients_Vicsek(κ;model="Fokker-Planck")
    λ = 1/κ
    I1 = quadgk(θ->cos(θ)*exp(κ*cos(θ)),0,pi,rtol=1e-5)[1]
    Z1 = quadgk(θ->exp(κ*cos(θ)),0,pi,rtol=1e-5)[1]
    c1 = I1/Z1
    if model=="BGK"
        I2 = quadgk(θ->cos(θ)*sin(θ)^2*exp(κ*cos(θ)),0,pi,rtol=1e-5)[1]
        Z2 = quadgk(θ->sin(θ)^2*exp(κ*cos(θ)),0,pi,rtol=1e-5)[1]
        c2 = I2/Z2
        return c1,c2,λ
    elseif model=="Fokker-Planck"
        Zinv = quadgk(θ->exp(-κ*cos(θ)),0,pi,rtol=1e-5)[1]
        I21 = quadgk(θ->θ*cos(θ)*sin(θ)*exp(κ*cos(θ)),0,pi,rtol=1e-5)[1]/κ
        I22 = pi/(κ*Zinv) * hcubature(x->cos(x[1])*sin(x[1])*exp(κ*cos(x[1]))*exp(-κ*cos(x[2]))*(x[2]<x[1]),(0,0),(pi,pi),rtol=1e-3)[1]
        I2 = I21 - I22
        Z21 = quadgk(θ->θ*sin(θ)*exp(κ*cos(θ)),0,pi,rtol=1e-5)[1]/κ
        Z22 = pi/(κ*Zinv) * hcubature(x->sin(x[1])*exp(κ*cos(x[1]))*exp(-κ*cos(x[2]))*(x[2]<x[1]),(0,0),(pi,pi),rtol=1e-3)[1]
        Z2 = Z21 - Z22
        c2 = I2/Z2
        return c1,c2,λ
    else
        ArgumentError("Model not defined!")
    end
end


text:
###
  Name     : c5_14_OLS.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 1/25/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

 use GLM, DataFrames
 data = DataFrame(X=[1,2,3], Y=[2,4,7])
 OLS = glm(@formula(Y ~ X), data, Normal(), IdentityLink())


text:
using SimpleSMC
using RNGPool
import SequentialMonteCarlo.SMCModel
using SMCExamples.FiniteFeynmanKac
import SMCExamples.Particles.Int64Particle
using Test
using Random

function Id(p::Int64Particle)
  return p.x
end

function One(p::Int64Particle)
  return 1.0
end

function IdMinus2(p::Int64Particle)
  return p.x - 2
end

function testEtas(smcio, ffkout, verbose::Bool)
  testapproxequal(
    SimpleSMC.allEtas(smcio, Id, false),
    FiniteFeynmanKac.allEtas(ffkout, Id, false), 0.01, verbose)
  testapproxequal(
    SimpleSMC.allEtas(smcio, Id, true),
    FiniteFeynmanKac.allEtas(ffkout, Id, true), 0.01, verbose)
end

function testZs(smcio, ffkout, verbose::Bool)
  testapproxequal(smcio.logZhats, ffkout.logZhats, 0.1, verbose)
end

function testapproxequal(x::Vector{Tuple{Bool,Float64}},
  y::Vector{Tuple{Bool,Float64}}, tol::Float64, verbose::Bool)
  xNonneg = (p->p[1]).(x)
  xLogVals = (p->p[2]).(x)
  yNonneg = (p->p[1]).(y)
  yLogVals = (p->p[2]).(y)
  @test xNonneg == yNonneg
  testapproxequal(xLogVals, yLogVals, tol, verbose)
end

function testGammas(smcio, ffkout, verbose::Bool)
  testapproxequal(SimpleSMC.allGammas(smcio, IdMinus2, false),
    FiniteFeynmanKac.allGammas(ffkout, IdMinus2, false), 0.1, verbose)
  testapproxequal(SimpleSMC.allGammas(smcio, IdMinus2, true),
    FiniteFeynmanKac.allGammas(ffkout, IdMinus2, true), 0.1, verbose)
end

function testFullOutput(ffk::FiniteFeynmanKac.FiniteFK, verbose::Bool)
  ffkout = FiniteFeynmanKac.calculateEtasZs(ffk)
  model = FiniteFeynmanKac.makeSMCModel(ffk)
  n = model.maxn
  smcio = SMCIO{model.particle, model.pScratch}(2^20, n)
  smc!(model, smcio)

  @testset "Finite FK: eta" begin
    testEtas(smcio, ffkout, verbose)
  end

  @testset "Finite FK: Zhat" begin
    testZs(smcio, ffkout, verbose)
  end

  @testset "Finite FK: slgamma" begin
    testGammas(smcio, ffkout, verbose)
  end

end

function _getFreqs(model::SMCModel, smcio::SMCIO, m::Int64, d::Int64)
  v::Vector{Int64Particle} = FiniteFeynmanKac.Int642Path(1, d, smcio.n)
  counts = zeros(Int64, d^smcio.n)
  result = Vector{Float64}(undef, d^smcio.n)
  for i = 1:m
    csmc!(model, smcio, v, v)
    counts[FiniteFeynmanKac.Path2Int64(v, d)] += 1
  end
  result .= counts ./ m
  return result
end

function _getFreqs(model::SMCModel, smcio::SMCIO, m::Int64, d::Int64, lM)
  v::Vector{Int64Particle} = FiniteFeynmanKac.Int642Path(1, d, smcio.n)
  counts = zeros(Int64, d^smcio.n)
  result = Vector{Float64}(undef, d^smcio.n)
  for i = 1:m
    csmc!(model, lM, smcio, v, v)
    counts[FiniteFeynmanKac.Path2Int64(v, d)] += 1
  end
  result .= counts ./ m
  return result
end

function testcsmc()
  d = 3
  n = 3
  ffk = FiniteFeynmanKac.randomFiniteFK(d, n)

  model = FiniteFeynmanKac.makeSMCModel(ffk)
  densities = Vector{Float64}(undef, d^n)
  for i = 1:length(densities)
    densities[i] = FiniteFeynmanKac.fullDensity(ffk, FiniteFeynmanKac.Int642Path(i, d, n))
  end
  densities ./= sum(densities)

  nsamples = 2^14

  smcio = SMCIO{model.particle, model.pScratch}(4, model.maxn)
  freqs = _getFreqs(model, smcio, nsamples, d)
  testapproxequal(freqs, densities, 0.05, false)

  lM = FiniteFeynmanKac.makelM(ffk)
  freqs = _getFreqs(model, smcio, nsamples, d, lM)
  testapproxequal(freqs, densities, 0.05, false)
end

setRNGs(0)

verbose = false

d = 3
ffk = FiniteFeynmanKac.randomFiniteFK(d, 10)

@time testFullOutput(ffk, verbose)

@time @testset "Finite FK: csmc" begin
  testcsmc()
end


text:
type PI_simple
	PI_simple()=begin
		[
			:(intTerm = 0),
		],
		[
			"",
		],
		new(
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Controller gain",
				:Default=>0.5
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Integral time constant",
				:Unit=>"s"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Previous scaled bias",
				:Default=>0.5
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>-1000
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>1000
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Previous scaled input signal",
				:Default=>0.5,
				:PosX=>0,
				:PosY=>0.5,
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Scaled output signal",
				:Default=>0,
				:PosX=>0.54,
				:PosY=>1,
				:Protected=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Scaled input variable",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Scaled set point",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Integral term",
				:Default=>0,
				:Protected=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Variable outp scaled between 0 and 1",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Scaled setPoint",
				:Lower=>0,
				:Upper=>1,
				:Default=>0.5
			)),
			[
				:(Ki*diff(intTerm) = setPoint - input),
				:(outps = bias + Kp*(setPoint - input) + intTerm),
				:(input*(MaxInput-MinInput) = Input-MinInput),
				:(setPoint*(MaxInput-MinInput) = SetPoint-MinInput),
				:(Output = 1),
				:(Output = 0),
				:(Output = outps),
			],
			[
				"Calculate integral term","Sum of proportional, integral and derivative terms","","","","","",
			],
			[:Kp,:Ki,:bias,:MinInput,:MaxInput,],
			[:Input,:Output,:input,:setPoint,:intTerm,:outps,:SetPoint,]
		)
	end
	Kp::DanaReal
	Ki::DanaReal
	bias::DanaReal
	MinInput::control_signal
	MaxInput::control_signal
	Input::control_signal
	Output::control_signal
	input::DanaReal
	setPoint::DanaReal
	intTerm::DanaReal
	outps::DanaReal
	SetPoint::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export PI_simple
function setEquationFlow(in::PI_simple)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	if outps > 1
		addEquation(5)
	else
		if outps < 0
			addEquation(6)
		else
			addEquation(7)
		end
	end
end
function initial(in::PI_simple)
	addEquation(1)
end
function atributes(in::PI_simple,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/PI"
	drive!(fields,_)
	return fields
end
PI_simple(_::Dict{Symbol,Any})=begin
	newModel=PI_simple()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
"""
$(TYPEDEF)
"""
struct StandardSDEProblem end

"""
$(TYPEDEF)
"""
struct SDEProblem{uType,tType,isinplace,P,NP,F,G,K,ND} <: AbstractSDEProblem{uType,tType,isinplace,ND}
  f::F
  g::G
  u0::uType
  tspan::tType
  p::P
  noise::NP
  kwargs::K
  noise_rate_prototype::ND
  seed::UInt64
  @add_kwonly function SDEProblem{iip}(f::AbstractSDEFunction{iip},g,u0,
          tspan,p=NullParameters();
          noise_rate_prototype = nothing,
          noise= nothing, seed = UInt64(0),
          kwargs...) where {iip}
    _tspan = promote_tspan(tspan)

    new{typeof(u0),typeof(_tspan),
        isinplace(f),typeof(p),
        typeof(noise),typeof(f),typeof(f.g),
        typeof(kwargs),
        typeof(noise_rate_prototype)}(
        f,f.g,u0,_tspan,p,
        noise,kwargs,
        noise_rate_prototype,seed)
  end

  function SDEProblem{iip}(f,g,u0,tspan,p=NullParameters();kwargs...) where {iip}
    SDEProblem(convert(SDEFunction{iip},f,g),g,u0,tspan,p;kwargs...)
  end
end

#=
function SDEProblem(f::AbstractSDEFunction,u0,tspan,p=NullParameters();kwargs...)
  SDEProblem(f,f.g,u0,tspan,p;kwargs...)
end
=#

function SDEProblem(f::AbstractSDEFunction,g,u0,tspan,p=NullParameters();kwargs...)
  SDEProblem{isinplace(f)}(f,g,u0,tspan,p;kwargs...)
end

function SDEProblem(f,g,u0,tspan,p=NullParameters();kwargs...)
  SDEProblem(convert(SDEFunction,f,g),g,u0,tspan,p;kwargs...)
end

"""
$(TYPEDEF)
"""
abstract type AbstractSplitSDEProblem end

"""
$(TYPEDEF)
"""
struct SplitSDEProblem{iip} <: AbstractSplitSDEProblem end
# u' = Au + f
function SplitSDEProblem(f1,f2,g,u0,tspan,p=NullParameters();kwargs...)
  SplitSDEProblem(SplitSDEFunction(f1,f2,g),g,u0,tspan,p;kwargs...)
end

SplitSDEProblem(f::SplitSDEFunction,g,u0,tspan,p=NullParameters();kwargs...) =
  SplitSDEProblem{isinplace(f)}(f,g,u0,tspan,p;kwargs...)

function SplitSDEProblem{iip}(f1,f2,g,u0,tspan,p=NullParameters();kwargs...) where iip
  SplitSDEProblem(SplitSDEFunction(f1,f2,g),g,u0,tspan,p;kwargs...)
end
function SplitSDEProblem{iip}(f::SplitSDEFunction,g,u0,tspan,p=NullParameters();
                                     func_cache=nothing,kwargs...) where iip
  if f.cache === nothing && iip
    cache = similar(u0)
    _f = SplitSDEFunction{iip}(f.f1, f.f2, f.g; mass_matrix=f.mass_matrix,
                              _func_cache=cache, analytic=f.analytic)
  else
    _f = f
  end
  SDEProblem(_f,g,u0,tspan,p;kwargs...)
end

"""
$(TYPEDEF)
"""
abstract type AbstractDynamicalSDEProblem end

"""
$(TYPEDEF)
"""
struct DynamicalSDEProblem{iip} <: AbstractDynamicalSDEProblem end

function DynamicalSDEProblem(f1,f2,g,v0,u0,tspan,p=NullParameters();kwargs...)
  DynamicalSDEProblem(DynamicalSDEFunction(f1,f2,g),g,v0,u0,tspan,p;kwargs...)
end

DynamicalSDEProblem(f::DynamicalSDEFunction,g,v0,u0,tspan,p=NullParameters();kwargs...) =
  DynamicalSDEProblem{isinplace(f)}(f,g,v0,u0,tspan,p;kwargs...)

function DynamicalSDEProblem{iip}(f1,f2,g,v0,u0,tspan,p=NullParameters();kwargs...) where iip
  DynamicalSDEProblem(DynamicalSDEFunction(f1,f2,g),g,v0,u0,tspan,p;kwargs...)
end
function DynamicalSDEProblem{iip}(f::DynamicalSDEFunction,g,v0,u0,tspan,p=NullParameters();
                                     func_cache=nothing,kwargs...) where iip
  if f.cache === nothing && iip
    cache = similar(u0)
    _f = DynamicalSDEFunction{iip}(f.f1, f.f2, f.g; mass_matrix=f.mass_matrix,
                              _func_cache=cache, analytic=f.analytic)
  else
    _f = f
  end
  SDEProblem(_f,g,ArrayPartition(v0,u0),tspan,p;kwargs...)
end


text:
module AlignedArrays
	import Mmap

	export AlignedArray, AlignedVector, AlignedMatrix, AlignedVecOrMat
	export PageAlignedArray, PageAlignedVector, PageAlignedMatrix, PageAlignedVecOrMat

	const PAGESIZE = Mmap.PAGESIZE

	struct AlignedArray{T, N, A} <: DenseArray{T, N}
		parent::Array{T, N}
		addr::Ref{Ptr{Cvoid}}

		function AlignedArray{T, N, A}(::UndefInitializer, dims::NTuple{N, Integer}) where {T, N, A}
			ispow2(A) || error("Alignment must be a power of two")
			isconcretetype(T) || error("Element type must be a concrete type")

			size = isempty(dims) ? 0 : reduce(*, dims)
			@static if Sys.islinux()
				addr = Ref(C_NULL)
				ccall(:posix_memalign, Cint, (Ptr{Ptr{Cvoid}}, Csize_t, Csize_t), addr, A, size) == 0 || error("Failed to allocate aligned memory")
			else
				error("Operating system not yet supported")
			end

			a = new{T, N, A}(unsafe_wrap(Array{T, N}, reinterpret(Ptr{T}, addr[]), dims, own = false), addr)
			finalizer(a.addr) do x
				@static if Sys.islinux()
					ccall(:free, Cvoid, (Ptr{Cvoid},), x[])
				end
			end
			return a
		end
	end

	AlignedArray{T, N, A}(u::UndefInitializer, dims::Integer...) where {T, N, A} = AlignedArray{T, N, A}(u, dims)

	const AlignedVector{T, A} = AlignedArray{T, 1, A}
	const AlignedMatrix{T, A} = AlignedArray{T, 2, A}
	const AlignedVecOrMat{T, A} = Union{AlignedVector{T, A}, AlignedMatrix{T, A}}

	const PageAlignedArray{T, N} = AlignedArray{T, N, PAGESIZE}
	const PageAlignedVector{T} = AlignedVector{T, PAGESIZE}
	const PageAlignedMatrix{T} = AlignedMatrix{T, PAGESIZE}
	const PageAlignedVecOrMat{T} = AlignedVecOrMat{T, PAGESIZE}

	Base.parent(a::AlignedArray) = a.parent

	Base.pointer(a::AlignedArray) = pointer(parent(a))

	Base.size(a::AlignedArray) = size(parent(a))
	Base.length(a::AlignedArray) = length(parent(a))
	Base.axes(a::AlignedArray) = axes(parent(a))

	Base.IndexStyle(::Type{A}) where {T, N, A<:AlignedArray{T, N}} = IndexStyle(Array{T, N})
	Base.getindex(a::AlignedArray, args...) = getindex(parent(a), args...)
	Base.setindex!(a::AlignedArray, args...) = setindex!(parent(a), args...)
	Base.iterate(a::AlignedArray, args...) = iterate(parent(a), args...)

	Base.similar(a::AlignedArray, args...) = similar(parent(a), args...)

	Base.show(io::IO, m::MIME"text/plain", a::AlignedArray) = show(io, m, parent(a))
end


text:
# Command: Get Title
"""
	document_title(session::Session)::String

This command returns the document title of the current top-level browsing context, equivalent to calling document.title.
"""
function document_title(session::Session)::String
    @unpack addr, id = session
    response = HTTP.get("$addr/session/$id/title", [("Content-Type" => "application/json; charset=utf-8")])
    @assert response.status == 200
    JSON3.read(response.body).value
end