text:
"""
    ForwardMapTriangulation{DT, TT}

A triangulation containing a triangulation of the

## Fields
- **`data`**. The data
- **`triang`**. A triangulation of points ``p_1, p_2, \\ldots, p_{N-1}``.
- **`forwardmap_triang`**. A triangulation of points ``p_2, p_2, \\ldots, p_{N}``.

"""
abstract type AbstractForwardMapTriangulation{DT, TT} end

ext: jl


text:
function parse_scalar(::Type{T}, e::XMLElement, name::String) where {T}
    parse(T, attribute(e, name))
end

function parse_scalar(::Type{T}, e::XMLElement, name::String, default::String) where {T}
    parse(T, e == nothing ? default : attribute(e, name))
end

function parse_vector(::Type{T}, e::Union{XMLElement, Nothing}, name::String, default::String) where {T}
    usedefault = e == nothing || attribute(e, name) == nothing # TODO: better handling of required attributes
    [parse(T, str) for str in split(usedefault ? default : attribute(e, name))]
end

function parse_inertia(::Type{T}, xml_inertia::XMLElement) where {T}
    ixx = parse_scalar(T, xml_inertia, "ixx", "0")
    ixy = parse_scalar(T, xml_inertia, "ixy", "0")
    ixz = parse_scalar(T, xml_inertia, "ixz", "0")
    iyy = parse_scalar(T, xml_inertia, "iyy", "0")
    iyz = parse_scalar(T, xml_inertia, "iyz", "0")
    izz = parse_scalar(T, xml_inertia, "izz", "0")
    @SMatrix [ixx ixy ixz; ixy iyy iyz; ixz iyz izz]
end

function parse_pose(::Type{T}, xml_pose::Nothing) where {T}
    rot = one(RotMatrix3{T})
    trans = zero(SVector{3, T})
    rot, trans
end

function parse_pose(::Type{T}, xml_pose::XMLElement) where {T}
    rpy = parse_vector(T, xml_pose, "rpy", "0 0 0")
    rot = RotMatrix(RotZYX(rpy[3], rpy[2], rpy[1]))
    trans = SVector{3}(parse_vector(T, xml_pose, "xyz", "0 0 0"))
    rot, trans
end

function parse_joint_type(::Type{T}, xml_joint::XMLElement,
        floating_joint_type::Type{<:JointType{T}}, revolute_joint_type::Type{<:JointType{T}}) where {T}
    joint_type = attribute(xml_joint, "type")
    if joint_type == "revolute" || joint_type == "continuous" # TODO: handle joint limits for revolute
        axis = SVector{3}(parse_vector(T, find_element(xml_joint, "axis"), "xyz", "1 0 0"))
        return revolute_joint_type(axis)
    elseif joint_type == "prismatic"
        axis = SVector{3}(parse_vector(T, find_element(xml_joint, "axis"), "xyz", "1 0 0"))
        return Prismatic(axis)
    elseif joint_type == "floating"
        return floating_joint_type()
    elseif joint_type == "fixed"
        return Fixed{T}()
    elseif joint_type == "planar"
        urdf_axis = SVector{3}(parse_vector(T, find_element(xml_joint, "axis"), "xyz", "1 0 0"))
        # The URDF spec says that a planar joint allows motion in a
        # plane perpendicular to the axis.
        R = Rotations.rotation_between(SVector(0, 0, 1), urdf_axis)
        return Planar{T}(R * SVector(1, 0, 0), R * SVector(0, 1, 0))
    else
        error("joint type $(string(joint_type)) not recognized")
    end
end

function parse_joint_bounds(jtype::JT, xml_joint::XMLElement) where {T, JT <: JointType{T}}
    position_bounds = fill(Bounds{T}(), num_positions(jtype))
    velocity_bounds = fill(Bounds{T}(), num_velocities(jtype))
    effort_bounds = fill(Bounds{T}(), num_velocities(jtype))
    for element in get_elements_by_tagname(xml_joint, "limit")
        if has_attribute(element, "lower")
            position_bounds .= Bounds.(parse_scalar(T, element, "lower"), upper.(position_bounds))
        end
        if has_attribute(element, "upper")
            position_bounds .= Bounds.(lower.(position_bounds), parse_scalar(T, element, "upper"))
        end
        if has_attribute(element, "velocity")
            v = parse_scalar(T, element, "velocity")
            velocity_bounds .= Bounds(-v, v)
        end
        if has_attribute(element, "effort")
            e = parse_scalar(T, element, "effort")
            effort_bounds .= Bounds(-e, e)
        end
    end
    position_bounds, velocity_bounds, effort_bounds
end

function parse_joint(::Type{T}, xml_joint::XMLElement, floating_joint_type, revolute_joint_type) where {T}
    name = attribute(xml_joint, "name")
    joint_type = parse_joint_type(T, xml_joint, floating_joint_type, revolute_joint_type)
    position_bounds, velocity_bounds, effort_bounds = parse_joint_bounds(joint_type, xml_joint)
    return Joint(name, joint_type; position_bounds=position_bounds, velocity_bounds=velocity_bounds, effort_bounds=effort_bounds)
end

function parse_inertia(::Type{T}, xml_inertial::XMLElement, frame::CartesianFrame3D) where {T}
    urdf_frame = CartesianFrame3D("inertia urdf helper")
    moment = parse_inertia(T, find_element(xml_inertial, "inertia"))
    com = zero(SVector{3, T})
    mass = parse_scalar(T, find_element(xml_inertial, "mass"), "value", "0")
    inertia = SpatialInertia(urdf_frame, moment, com, mass)
    pose = parse_pose(T, find_element(xml_inertial, "origin"))
    transform(inertia, Transform3D(urdf_frame, frame, pose...))
end

function parse_body(::Type{T}, xml_link::XMLElement, frame::CartesianFrame3D = CartesianFrame3D(attribute(xml_link, "name"))) where {T}
    xml_inertial = find_element(xml_link, "inertial")
    inertia = xml_inertial == nothing ? zero(SpatialInertia{T}, frame) : parse_inertia(T, xml_inertial, frame)
    linkname = attribute(xml_link, "name") # TODO: make sure link name is unique
    RigidBody(linkname, inertia)
end

function parse_root_link(mechanism::Mechanism{T}, xml_link::XMLElement, root_joint_type::JointType{T}=Fixed{T}()) where {T}
    parent = root_body(mechanism)
    body = parse_body(T, xml_link)
    joint = Joint("$(string(body))_to_world", root_joint_type)
    joint_to_parent = one(Transform3D{T}, frame_before(joint), default_frame(parent))
    attach!(mechanism, parent, body, joint, joint_pose = joint_to_parent)
end

function parse_joint_and_link(mechanism::Mechanism{T}, xml_parent::XMLElement, xml_child::XMLElement, xml_joint::XMLElement,
        floating_joint_type, revolute_joint_type) where {T}
    parentname = attribute(xml_parent, "name")
    candidate_parents = collect(filter(b -> string(b) == parentname, non_root_bodies(mechanism))) # skip root (name not parsed from URDF)
    length(candidate_parents) == 1 || error("Duplicate name: $(parentname)")
    parent = first(candidate_parents)
    joint = parse_joint(T, xml_joint, floating_joint_type, revolute_joint_type)
    pose = parse_pose(T, find_element(xml_joint, "origin"))
    joint_to_parent = Transform3D(frame_before(joint), default_frame(parent), pose...)
    body = parse_body(T, xml_child, frame_after(joint))
    attach!(mechanism, parent, body, joint, joint_pose = joint_to_parent)
end

"""
$(SIGNATURES)

Create a `Mechanism` by parsing a [URDF](https://wiki.ros.org/urdf/XML/model) file.

Keyword arguments:

* `scalar_type`: the scalar type used to store the `Mechanism`'s kinematic and inertial properties. Default: `Float64`.
* `floating`: whether to use a floating joint as the root joint. Default: false.
* `floating_joint_type`: what `JointType` to use for floating joints. Default: `QuaternionFloating{scalar_type}`.
* `revolute_joint_type`: what `JointType` to use for revolute joints. Default: `Revolute{scalar_type}`.
* `root_joint_type`: the joint type used to connect the parsed `Mechanism` to the world. Default: `floating_joint_type()` if `floating`, `Fixed{scalar_type}()` otherwise.
* `remove_fixed_tree_joints`: whether to remove any fixed joints present in the kinematic tree using [`remove_fixed_tree_joints!`](@ref). Default: `true`.
* `gravity`: gravitational acceleration as a 3-vector expressed in the `Mechanism`'s root frame. Default: `$(DEFAULT_GRAVITATIONAL_ACCELERATION)`.
"""
function parse_urdf(filename::AbstractString;
        scalar_type::Type{T}=Float64,
        floating=false,
        floating_joint_type::Type{<:JointType{T}}=QuaternionFloating{scalar_type},
        revolute_joint_type::Type{<:JointType{T}}=Revolute{scalar_type},
        root_joint_type::JointType{T}=floating ? floating_joint_type() : Fixed{scalar_type}(),
        remove_fixed_tree_joints=true,
        gravity::AbstractVector=DEFAULT_GRAVITATIONAL_ACCELERATION) where T

    if floating && !isfloating(root_joint_type)
        error("Ambiguous input arguments: `floating` specified, but `root_joint_type` is not a floating joint type.")
    end

    xdoc = parse_file(filename)
    xroot = LightXML.root(xdoc)
    @assert LightXML.name(xroot) == "robot"

    xml_links = get_elements_by_tagname(xroot, "link")
    xml_joints = get_elements_by_tagname(xroot, "joint")

    # create graph structure of XML elements
    graph = DirectedGraph{Vertex{XMLElement}, Edge{XMLElement}}()
    vertices = Vertex.(xml_links)
    for vertex in vertices
        add_vertex!(graph, vertex)
    end
    name_to_vertex = Dict(attribute(v.data, "name") => v for v in vertices)
    for xml_joint in xml_joints
        parent = name_to_vertex[attribute(find_element(xml_joint, "parent"), "link")]
        child = name_to_vertex[attribute(find_element(xml_joint, "child"), "link")]
        add_edge!(graph, parent, child, Edge(xml_joint))
    end

    # create a spanning tree
    roots = collect(filter(v -> isempty(in_edges(v, graph)), vertices))
    length(roots) != 1 && error("Can only handle a single root")
    tree = SpanningTree(graph, first(roots))

    # create mechanism from spanning tree
    rootbody = RigidBody{T}("world")
    mechanism = Mechanism(rootbody, gravity=gravity)
    parse_root_link(mechanism, Graphs.root(tree).data, root_joint_type)
    for edge in edges(tree)
        parse_joint_and_link(mechanism, source(edge, tree).data, target(edge, tree).data, edge.data, floating_joint_type, revolute_joint_type)
    end

    if remove_fixed_tree_joints
        remove_fixed_tree_joints!(mechanism)
    end

    mechanism
end

@noinline function parse_urdf(scalar_type::Type, filename::AbstractString)
    # TODO: enable deprecation:
    # replacement = if scalar_type == Float64
    #     "parse_urdf(filename, remove_fixed_tree_joints=false)"
    # else
    #     "parse_urdf(filename, scalar_type=$scalar_type, remove_fixed_tree_joints=false)"
    # end
    # msg = """
    # `parse_urdf(scalar_type, filename)` is deprecated, use $replacement instead.
    # This is to reproduce the exact same behavior as before.
    # You may want to consider leaving `remove_fixed_tree_joints` to its default value (`true`).
    # """
    # Base.depwarn(msg, :parse_urdf)
    parse_urdf(filename; scalar_type=scalar_type, remove_fixed_tree_joints=false)
end


text:
# These routines were modified from similar routines in MATLAB.jl
# E.G., @mput, @mget, _mput_multi, _mget_multi, make_getvar_statement

function put_var_multi(vs::Symbol...)
    nv = length(vs)
    if nv == 1
        v = vs[1]
        :(put_var($(v), string($(Meta.quot(v)))))
    else
        stmts = Array(Expr, nv)
        for i = 1:nv
            v = vs[i]
            stmts[i] = :(put_var($(v), string($(Meta.quot(v)))))
        end
        Expr(:block, stmts...)
    end
end

macro put_var(vs...)
    esc(put_var_multi(vs...))
end

function make_getvar_statement(v::Symbol)
    :($(v) = get_var(string($(Meta.quot(v)))))
end

function make_getvar_statement(ex::Expr)
    if !(ex.head == :(=))
        if ex.head == :kw
            error("Must call @get_var without parenthesis if using statements")
        else
            erorr("Invalid expression for @get_var: " * string(ex))
        end
    end
    v::Symbol = ex.args[1]
    k::Symbol = ex.args[2]
    :($(v) = get_var(string($(Meta.quot(k)))))
end

function get_var_multi(vs::Union{Symbol, Expr}...)
    # NOTE: I supressed output by adding :nothing
    nv = length(vs)
    if nv == 1
        stmt = make_getvar_statement(vs[1])
        Expr(:block, stmt, :nothing)
    else
        stmts = Array(Expr, nv)
        for i = 1:nv
            stmts[i] = make_getvar_statement(vs[i])
        end
        Expr(:block, stmts..., :nothing)
    end
end

macro get_var(vs...)
    esc(get_var_multi(vs...))
end

# END: routines modified from similar routines in MATLAB.jl


text:
function setup_allsky_as(
    ds_in,
    idx_gases,
    lkp_lw,
    lkp_sw,
    lkp_lw_cld,
    lkp_sw_cld,
    use_lut,
    ncol,
    FT,
    DA,
    max_threads,
)

    deg2rad = FT(π) / FT(180)
    nlay = Int(ds_in.dim["lay"])
    #ncol = Int(ds_in.dim["col"]) # col#1 repeated 128 times, per RRTMGP example
    nlev = nlay + 1
    ngas = lkp_lw.n_gases
    nbnd_lw = lkp_lw.n_bnd
    nbnd_sw = lkp_sw.n_bnd

    #---no lat / long information for this
    lon = nothing
    lat = nothing
    # The example only reads the first column and
    # replicates it ncol times.

    sfc_emis = DA{FT,2}(undef, nbnd_lw, ncol)
    sfc_alb_direct = DA{FT,2}(undef, nbnd_sw, ncol)
    sfc_alb_diffuse = DA{FT,2}(undef, nbnd_sw, ncol)
    zenith = DA{FT,1}(undef, ncol)
    irrad = DA{FT,1}(undef, ncol)
    # these values are taken from the example
    sfc_emis .= FT(0.98)
    sfc_alb_direct .= FT(0.06)
    sfc_alb_diffuse .= FT(0.06)
    zenith .= acos(FT(0.86))
    irrad .= FT(lkp_sw.solar_src_tot)

    p_lev = Array{FT}(reshape(ds_in["p_lev"][:][1, :], nlev, 1))

    bot_at_1 = p_lev[1, 1] > p_lev[end, 1]
    lev_ind = bot_at_1 ? (1:nlev) : (nlev:-1:1)
    lay_ind = bot_at_1 ? (1:nlay) : (nlay:-1:1)

    p_lev = Array{FT}(reshape(ds_in["p_lev"][:][1, lev_ind], nlev, 1))
    p_lay = Array{FT}(reshape(ds_in["p_lay"][:][1, lay_ind], nlay, 1))
    t_lev = Array{FT}(reshape(ds_in["t_lev"][:][1, lev_ind], nlev, 1))
    t_lay = Array{FT}(reshape(ds_in["t_lay"][:][1, lay_ind], nlay, 1))
    t_sfc = Array{FT}(reshape([t_lev[1, 1]], 1))

    p_lev = repeat(p_lev, 1, ncol)
    p_lay = repeat(p_lay, 1, ncol)
    t_lev = repeat(t_lev, 1, ncol)
    t_lay = repeat(t_lay, 1, ncol)
    t_sfc = repeat(t_sfc, ncol)
    #col_dry = DA{FT,2}(transpose(ds_in["col_dry"][:][:, lay_ind]))
    #col_dry from the dataset not used in the FORTRAN RRTMGP example

    # Reading volume mixing ratios
    vmrat = zeros(FT, nlay, ncol, ngas)

    vmrat[:, 1, idx_gases["h2o"]] .= Array{FT}(ds_in["vmr_h2o"][:][1, lay_ind])
    vmrat[:, 1, idx_gases["o3"]] .= Array{FT}(ds_in["vmr_o3"][:][1, lay_ind])
    vmrat[:, 1, idx_gases["co2"]] .= Array{FT}(ds_in["vmr_co2"][:][1, lay_ind])
    vmrat[:, 1, idx_gases["n2o"]] .= Array{FT}(ds_in["vmr_n2o"][:][1, lay_ind])
    vmrat[:, 1, idx_gases["co"]] .= Array{FT}(ds_in["vmr_co"][:][1, lay_ind])
    vmrat[:, 1, idx_gases["ch4"]] .= Array{FT}(ds_in["vmr_ch4"][:][1, lay_ind])
    vmrat[:, 1, idx_gases["o2"]] .= Array{FT}(ds_in["vmr_o2"][:][1, lay_ind])
    vmrat[:, 1, idx_gases["n2"]] .= Array{FT}(ds_in["vmr_n2"][:][1, lay_ind])

    for icol = 2:ncol
        vmrat[:, icol, :] .= vmrat[:, 1, :]
    end
    vmr = Vmr(DA(vmrat))
    col_dry = DA{FT,2}(undef, nlay, ncol)
    vmr_h2o = view(vmr.vmr, :, :, idx_gases["h2o"])

    cld_mask = zeros(Bool, nlay, ncol)
    cld_r_eff_liq = zeros(FT, nlay, ncol)
    cld_r_eff_ice = zeros(FT, nlay, ncol)
    cld_path_liq = zeros(FT, nlay, ncol)
    cld_path_ice = zeros(FT, nlay, ncol)

    if use_lut
        r_eff_liq = (lkp_lw_cld.radliq_lwr + lkp_lw_cld.radliq_upr) / FT(2)
        r_eff_ice = (lkp_lw_cld.radice_lwr + lkp_lw_cld.radice_upr) / FT(2)
    else
        pade_sizreg_extliq = Array(lkp_lw_cld.pade_sizreg_extliq) # TODO temp fix to avoid scalar indexing on GPU
        pade_sizreg_extice = Array(lkp_lw_cld.pade_sizreg_extice)
        r_eff_liq = sum(extrema(pade_sizreg_extliq)) / FT(2)
        r_eff_ice = sum(extrema(pade_sizreg_extice)) / FT(2)
    end
    # Restrict clouds to troposphere (> 100 hPa = 100*100 Pa)
    # and not very close to the ground (< 900 hPa), and
    # put them in 2/3 of the columns since that's roughly the
    # total cloudiness of earth
    for icol = 1:ncol, ilay = 1:nlay
        if p_lay[ilay, icol] > FT(10000) &&
           p_lay[ilay, icol] < FT(90000) &&
           icol % 3 ≠ 0
            cld_mask[ilay, icol] = true
            if t_lay[ilay, icol] > FT(263)
                cld_path_liq[ilay, icol] = FT(10)
                cld_r_eff_liq[ilay, icol] = r_eff_liq
            end
            if t_lay[ilay, icol] < FT(273)
                cld_path_ice[ilay, icol] = FT(10)
                cld_r_eff_ice[ilay, icol] = r_eff_ice
            end
        end
    end

    p_lay = DA(p_lay)
    p_lev = DA(p_lev)
    t_lay = DA(t_lay)
    t_lev = DA(t_lev)

    compute_col_dry!(p_lev, col_dry, param_set, vmr_h2o, lat) # the example skips lat based gravity calculation

    t_sfc = DA(t_sfc)

    cld_mask = DA(cld_mask)
    cld_r_eff_liq = DA(cld_r_eff_liq)
    cld_r_eff_ice = DA(cld_r_eff_ice)
    cld_path_liq = DA(cld_path_liq)
    cld_path_ice = DA(cld_path_ice)
    ice_rgh = 2 # medium ice roughness

    return (
        AtmosphericState{
            FT,
            typeof(t_sfc),
            typeof(lat),
            typeof(p_lev),
            typeof(cld_r_eff_liq),
            typeof(cld_mask),
            typeof(vmr),
            Int,
        }(
            lon,
            lat,
            p_lay,
            p_lev,
            t_lay,
            t_lev,
            t_sfc,
            col_dry,
            vmr,
            cld_r_eff_liq,
            cld_r_eff_ice,
            cld_path_liq,
            cld_path_ice,
            cld_mask,
            ice_rgh,
            nlay,
            ncol,
            ngas,
        ),
        sfc_emis,
        sfc_alb_direct,
        sfc_alb_diffuse,
        zenith,
        irrad,
        bot_at_1,
    )
end


text:
using DifferentialGeometries
using ForwardDiff, LinearAlgebra, Random, Test

Random.seed!(0)

@testset "DifferentialGeometries tests" begin
  include("Cartesian.jl")
  include("Cylindrical.jl")
  include("Toroidal.jl")
  include("ChangeOfCoordinates.jl")
  include("AutodiffPropagation.jl")
end

max_stars_repo_name: jwscook/DifferentialGeometries

max_issues_repo_name: jwscook/DifferentialGeometries

max_forks_repo_name: jwscook/DifferentialGeometries


text:
function _fatal_err(f::Function, gwt::GWTask)
    try; f()
    catch err
        _printerr(err)
        _print_eotask(gwt)

        _flush()

        if is_worker_mode(gwt)
            _up_task_status!(gwt, _GW_TASK_ERROR_STATUS)
            # TODO: del proc reg
        end

        exit()
    end
end

function _expr_src(ex::Expr)
    buf = IOBuffer()
    Base.show_unquoted(buf, ex)
    src = String(take!(buf))
    src = replace(src, r"#.*#\n" => "")
    return string("# NOTE: This source is not 'original', it was generated from an `Expr`.", "\n", src)
end

_expr_src(gwt::GWTask) = _expr_src(_task_expr(gwt))

## ------------------------------------------------------
const _GW_TASK_T0_KEY = "_T0"
_tic!(gwt::GWTask) = set!(gwt, _GW_TASK_T0_KEY, now())
_toc(gwt::GWTask) = now() - get(gwt, _GW_TASK_T0_KEY, now())

## ------------------------------------------------------
const _GW_WELCOME_TOKEN = rpad("RUNNING TASK ", 60, ">")
const _GW_EOTASK_TOKEN = rpad("END OF TASK ", 60, ">")
const _GW_HEAD_SEP = rpad("", 60, "<")

function _print_welcome(gwt::GWTask)
    _tic!(gwt)
    println(_GW_WELCOME_TOKEN)
    # println()
    println("taskid             ", task_id(gwt))
    println("pid                ", getpid())
    println("pwd                ", pwd())
    println("wroker mode        ", is_worker_mode(gwt))
    println("start time         ", now())
    # println()
    println(_GW_HEAD_SEP)
    println("\n"^2)
    _flush()
end

function _print_eotask(gwt::GWTask)
    _flush()
    println("\n"^2)
    println(_GW_EOTASK_TOKEN)
    println("end time           ", now())
    cantime = Dates.canonicalize(Dates.CompoundPeriod(_toc(gwt)))
    println("tot time           ", cantime)
    println(_GW_HEAD_SEP)
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
__precompile__(true)

module Printff

# the implementations here is copied from
# what is left in base/printf.jl, and uses the utility there

export @printf, @sprintf, @format_str
export format, format_expr, printf, sprintf

import Base.Printf: is_str_expr, fix_dec, DIGITS, print_fixed, decode_dec, decode_hex,
                   ini_hex, ini_HEX, print_exp_a, decode_0ct, decode_HEX, ini_dec, print_exp_e,
                   decode_oct, _limit

using Base.Printf: gen_f, gen_e, gen_a, gen_g, gen_c, gen_s, gen_p, gen_d
using Unicode:  lowercase, textwidth

# copied from Base.Printf - amended by additional features:
# 1. optional argument positions %(\d+[$&])
# 2. special term %\d+% standing for %\d+&s

### generate code from parser output ###

function gen(s::AbstractString)
    args = []
    perm = Int[]
    blk = Expr(:block, :(local neg, pt, len, exp, do_out, args))
    blkst = blk.args
    for x in parse(s)
        if isa(x, AbstractString)
            push!(blkst, :(write(out, $(length(x)==1 ? x[1] : x))))
        else
            c = lowercase(x[end])
            f = c == 'f' ? gen_f :
                c == 'e' ? gen_e :
                c == 'a' ? gen_a :
                c == 'g' ? gen_g :
                c == 'c' ? gen_c :
                c == 's' ? gen_s :
                c == 'p' ? gen_p :
                           gen_d
            arg, ex = f(x[2:end]...)
            push!(args, arg)    # variable name and type ':(var::type)'
            push!(perm, x[1])   # position number
            push!(blkst, ex)    # executable statement write to 'out'
        end
    end
    push!(blkst, :nothing) # end of statement list in blk
    return args, blk, perm
end

### printf format string parsing ###

function parse(s::AbstractString)
    # parse format string into strings and format tuples
    curpos = 0
    list = []
    i = j = start(s)
    j1 = 0 # invariant: j1 == prevind(s, j)
    while !done(s,j)
        c, k = next(s,j)
        if c == '%'
            i > j1 || push!(list, s[i:j1])
            flags, width, precision, conversion, k, pos = parse1(s,k)
            '\'' in flags && throw(ArgumentError("printf format flag ' not yet supported"))
            conversion == 'n'    && throw(ArgumentError("printf feature %n not supported"))
            pos != 0 && conversion == '%' && ( conversion = 's' )
            if conversion == '%'
                push!(list, "%")
            else
                curpos = pos == 0 ? curpos + 1 : pos
                push!(list, (curpos,flags,width,precision,conversion))
            end
            i = k
        end
        j1 = j
        j = k
    end
    i > endof(s) || push!(list, s[i:end])
    # coalesce adjacent strings
    i = 1
    while i < length(list)
        if isa(list[i], AbstractString)
            for outer j = i+1:length(list)
                if !isa(list[j], AbstractString)
                    j -= 1
                    break
                end
                list[i] *= list[j]
            end
            deleteat!(list, i+1:j)
        end
        i += 1
    end
    return list
end

## step forward or bail out at end
function next_or_die(s::AbstractString, k)
    if !done(s,k)
        next(s,k)
    else
        throw(ArgumentError("invalid printf format string: $(repr(s))"))
    end
end

## parse a single printf specifier ##
# printf specifiers:
#   %                       # start
#   (\d+\[$&])?             # arg position
#   [\-\+#0' ]*             # flags
#   (\d+)?                  # width
#   (\.\d*)?                # precision
#   (h|hh|l|ll|L|j|t|z|q)?  # modifier (ignored)
#   [diouxXeEfFgGaAcCsSp%]  # conversion
#
# returns tuple (flags, width, precision, c, k, pos)
function parse1(s::AbstractString, k::Integer)
    j = k
    width = 0
    precision = -1
    pos = 0
    c, k = next_or_die(s,k)
    # handle %%
    if c == '%'
        return "", width, precision, c, k, pos
    end
    # look for optional argument position
    ca, ka = c, k
    # if we didn't allow 0 as first char of position no traceback rqd.
    while '0' <= ca <= '9'
        pos = 10*pos + ca-'0'
        ca, ka = next_or_die(s,ka)
    end
    if ka != k
        if ca == '$' || ca == '&' # ending position number - continue format
            j = ka
            c, k = next_or_die(s,ka)
        elseif ca == '%'          # ending position and default format
            return "s", width, precision, ca, ka, pos
        else
            pos = 0 # backtrace
        end
    end

    # parse flags
    while c in "#0- + '"
        c, k = next_or_die(s,k)
    end
    flags = String(s[j:prevind(s,k)-1]) # exploiting that all flags are one-byte.
    # parse width
    while '0' <= c <= '9'
        width = 10*width + c-'0'
        c, k = next_or_die(s,k)
    end
    # parse precision
    if c == '.'
        c, k = next_or_die(s,k)
        if '0' <= c <= '9'
            precision = 0
            while '0' <= c <= '9'
                precision = 10*precision + c-'0'
                c, k = next_or_die(s,k)
            end
        end
    end
    # parse length modifer (ignored)
    if c == 'h' || c == 'l'
        prev = c
        c, k = next_or_die(s,k)
        if c == prev
            c, k = next_or_die(s,k)
        end
    elseif c in "Ljqtz"
        c, k = next_or_die(s,k)
    end
    # validate conversion
    if !(c in "diouxXDOUeEfFgGaAcCsSp")
        throw(ArgumentError("invalid printf format string: $(repr(s))"))
    end
    # warn about silly flag/conversion combinations
    warn_unused(flags, width, precision, c)
    flags, width, precision, c, k, pos
end

function warn_unused(flags::String, width::Int, prec::Int, c::Char)
    c = lowercase(c)
    unc = c in "diucp" && '#' in flags
    unf = c in "csp" && flags != "" &&
        !( '-' in flags  || (c == 'c' && '0' in flags ) || (c == 's' && '#' in flags))
    unp = c in "cp" && prec >= 0

    if unc || unf || unp
        wt = width > 0 ? string(width) : ""
        pt = prec >= 0 ? string('.', prec) : ""
        tc = unf ? "flags " : ""
        tp = unp ? "precision " : ""
        text = string("unused ", tc, tp, "in format spec %", flags, wt, pt, c)
        throw(ArgumentError(text)) # or @warn text
    end
end

# generate code for macro expansion
# as format allow plain string literals (no interpolation) and raw string literals
#
function printf_expr(macroname, line, io, fmt, args)
    if isa(fmt, Expr) && fmt.head == :macrocall && fmt.args[1] == Symbol("@raw_str")
        fmt = fmt.args[end]
    end
    if !isa(fmt, AbstractString)
        throw(ArgumentError("$macroname: format must be a plain static string (no interpolation or macro prefix)"))
    end

    sargs, blk, perm = gen(fmt)

    bind_vars!(macroname, true, io, sargs, blk, perm, args)
end

## check relationship between arguments in format and provided data
# 1. for each format specifier ther must by exactly one data value
# 2. each data item has to be output by the format
# 3. the same data item may be re-used for several format specifiers
#
function check_args(mana, argl, argmax, perm)
    if argl != argmax
        throw(ArgumentError("$mana: wrong number of arguments ($argl) should be ($argmax)"))
    end
    if !( 1:argl ⊆ perm )
        throw(ArgumentError("$mana: invalid permutation '$perm'"))
    end
end

## Add code to pair symbolic args with data items sarg[i] := arg[perm[i]] by assignment
# macroname: identification for error messages
# ism: is called from macro
# io: output IO
# sargs: contain variable names and data types as provided by gen-subroutines
# blk: generated code
# perm: list of argument-indices associated to the variable names
# args: arguments from macro call or (already evaluated) from format-call
#
# note: The statements are added in reverse order of later execution
#       art the beginning of the given block.
#
function bind_vars!(macroname, ism::Bool, io, sargs, blk, perm, args)

    has_splatting = ism && any(arg -> isa(arg, Expr) && arg.head == :..., args)
    argmax = max(0, perm...)
    #
    #  Immediately check for corresponding arguments if there is no splatting
    #
    blkst = blk.args # the executable statements of the block
    if !has_splatting
        argl = length(args)
        check_args(macroname, argl, argmax, perm)
    end

    for i = length(sargs):-1:1
        var = sargs[i].args[1]
        if has_splatting
            # add assignment statement with values stored in runtime variable G
            pushfirst!(blkst, :($var = G[$(perm[i])]))
        else
            rhs = args[perm[i]]
            rhs = ism ? esc(rhs) : rhs
            if rhs != var
                # add assignment statement with constant values
                pushfirst!(blkst, :($var = $rhs))
            end
        end
    end

    #
    #  Delay generation of argument list and check until evaluation time instead of macro
    #  expansion time if there is splatting.
    #
    if has_splatting
       x = Expr(:call, :tuple, args...) # all arguments expanded into one tuple
       pushfirst!(blkst,
          quote
             G = $(esc(x))  # evaluate at runtime
             check_args($macroname, length(G), $argmax, $perm)
          end
       )
    end

    io != :out && pushfirst!(blkst, :(out = $io))
    Expr(:let, Expr(:block), blk)
end

"""
    @printf([io::IO,] "%Fmt", args...)

Print `args` using C `printf` style format specification string, with some caveats:
`Inf` and `NaN` are printed consistently as `Inf` and `NaN` for flags `%a`, `%A`,
`%e`, `%E`, `%f`, `%F`, `%g`, and `%G`. Furthermore, if a floating point number is
equally close to the numeric values of two possible output strings, the output
string further away from zero is chosen.
The variable length-specifiers `*` are not supported.
The argument-numbering syntax of C and C++ boost is supported, where the number delimiter
`\$` may be replaced by `&`.

Optionally, an `IO` stream of buffer
may be passed as the first argument to redirect output.

# Examples
```jldoctest
julia> @printf("%f %F %f %F\\n", Inf, Inf, NaN, NaN)
Inf Inf NaN NaN\n

julia> @printf "%.0f %.1f %f\\n" 0.5 0.025 -0.0078125
1 0.0 -0.007813
```
"""
macro printf(args...)
    isempty(args) && throw(ArgumentError("@printf: called with no arguments"))
    if isa(args[1], AbstractString) || is_str_expr(args[1])
        printf_expr("@printf", __source__, STDOUT, args[1], args[2:end])
    else
        (length(args) >= 2 && (isa(args[2], AbstractString) || is_str_expr(args[2]))) ||
            throw(ArgumentError("@printf: first or second argument must be a format string"))
        printf_expr("@printf", __source__, esc(args[1]), args[2], args[3:end])
    end
end

"""
    @sprintf("%Fmt", args...)

Return `@printf` formatted output as string.

# Examples
```jldoctest
julia> s = @sprintf "this is a %2% s %1&15.1f" 34.567 "test";

julia> println(s)
this is a test            34.6
```
"""
macro sprintf(args...)
    isempty(args) && throw(ArgumentError("@sprintf: called with zero arguments"))
    isa(args[1], AbstractString) || is_str_expr(args[1]) ||
        throw(ArgumentError("@sprintf: first argument must be a format string"))
    letexpr = printf_expr("@sprintf", __source__, :(IOBuffer()), args[1], args[2:end])
    push!(letexpr.args[2].args, :(String(take!(out))))
    letexpr
end
"""
    format" literal format string"

Generate a format function f, which can be called like f(io, args...).
That has the same effect as `@printf(io, "format string", args..)` while the generated
code is not inlined, but in a separate local function.
The format function is generated only once for each different format string.
"""
macro format_str(arg)
    format_expr(__source__, arg)
end

##### formatting functions
## kind of type intersection, base on the symbols :Any, :Real, :Integer, :Ptr
# note: gen_d returns :Real (not :Integer)
#       gen_c returns :Integer (not :Char)
#
function typeinter(a::Symbol, b::Symbol)
    a === :Any ? b :
    b === :Any ? a :
    a === b ?    a :
    a === :Real ? b :
    b === :Real ? a :
                :PIC
end

## map to the actual argument types according to format specs
#
function typearg(a::Symbol)
    a === :Any ? a :
    a === :Real ? a :
    a === :Integer ? Expr(:curly, :Union, :Char, a) :
    a === :Ptr ? Expr(:curly, :Union, :Integer, a) :
                :Integer
end

## for all arguments determine the formal argument types
# to be accepted by the function
# Re-use the varaible names of the symbolic arguments
# If multiple format specs are delivered by the same argument
# the argument type must be restricted reasonably
#
function typeconsolidate(sargs, perm)
    n = max(0, perm...)
    m = length(perm)
    fargs = Vector{Any}(nothing, n)

    for i = 1:m
        j = perm[i]
        faj = fargs[j]
        if faj === nothing
            fargs[j] = sargs[i] # :(var::type)
        else
            ft = faj.args[2]
            st = sargs[i].args[2]
            if st != :Any
                if ft == :Any
                    faj.args[2] = st
                else
                    faj.args[2] = typeinter(ft, st)
                end
            end
        end
    end
    for j = 1:n
        faj = fargs[j]
        if faj === nothing
            fargs[j] = :(::Any) # unused argument without variable name
        else
            faj.args[2] = typearg(faj.args[2]) # final type
        end
    end
    fargs
end

## generate function definition for format function
#
function format_expr(line::LineNumberNode, fmt::AbstractString)
    sargs, blk, perm = gen(fmt)
    fargs = typeconsolidate(sargs, perm)
    args = map(p->p.args[1], fargs)
    blk = bind_vars!("format", false, :out, sargs, blk, perm, args)

    fun = :($(genfun())(out::IO) = nothing)
    append!(fun.args[1].args, fargs)
    fun.args[2].args[2] = blk
    fun.args[2].args[1] = line
    fun
end

# generate anonymous function
function format_fun(line::LineNumberNode, fmt)
    f(io::IO, args...) = Base.invokelatest(eval(format_expr(line, fmt)), io, args...)
end
"""
    format(fmt::AbstractString) -> Function(::IO, args...)

Generate a format function f, which can be called like f(io, args...).
That has the same effect as `printf(io, "format string", args..)`.
The format function is generated only once for each different format string.

The variable length-specifiers `*` are not supported.
The argument-numbering syntax of C and C++ boost is supported, where the number delimiter
`\$` may be replaced by `&`.
The arguments are numberd starting with `1`. The arguments may be explicitly referred to
from the format string. In each format specification after the starting `%` and optional
argument number may be given in decimal digits, followed by `&` or `\$`. As an additional
feature copied from C++, `%<digits>%` behaves like `%<digits>&s`.

The number and type of the arguments is restricted by the format. All arguments need to
be used by the format. Each argument may be used multiple times.
The argument types are determined by the conversion character in the format:
```
    %c                              Char and Integer
    %d                              Real - printed value is rounded to Integer
    %a, %A, %e, %f, %F, %g, %G      Real
    %s                              Any
    %p                              Ptr and Integer
```
"""
function format(fmt::AbstractString)
    line = LineNumberNode(@__LINE__, Symbol(joinpath(@__DIR__, @__FILE__)))
    global ALL_FORMATS
    get!(ALL_FORMATS, fmt) do
        f = format_fun(line, fmt)
    end
end
"""
    printf([io::IO,] fmt, args...)
    sprintf(fmt, args) -> String

Print `args` using C `printf` style format specification string, with some caveats:
`Inf` and `NaN` are printed consistently as `Inf` and `NaN` for flags `%d`, `%s`, `%a`, `%A`,
`%e`, `%E`, `%f`, `%F`, `%g`, and `%G`. Furthermore, if a floating point number is
equally close to the numeric values of two possible output strings, the output
string further away from zero is chosen.
For further details of the format see [`format`](@ref).
The format `fmt` is not restricted to literal strings without interpolation syntax as is
the case for the corresponding macro calls.

Optionally, an `IO` stream or buffer
may be passed as the first argument to redirect output.

# Examples
```jldoctest
julia> printf(STDERR, "%f %F %f %F\\n", Inf, Inf, NaN, NaN)
Inf Inf NaN NaN\n

julia> printf("%2&.0f %1\\\$.1f %3&f\\n", 0.025, 0.5, -0.0078125)
1 0.0 -0.007813
```
"""
printf(io::IO, fmt::Function, args...) = fmt(io, args...)
printf(fmt::Function, args...) = printf(STDOUT, fmt, args...)
function sprintf(fmt::Function, args...)
    io = IOBuffer()
    printf(io, fmt, args...)
    String(take!(io))
end

printf(io::IO, fmt::AbstractString, args...) = printf(io, format(fmt), args...)
printf(fmt::AbstractString, args...) = printf(format(fmt), args...)
sprintf(fmt::AbstractString, args...) = sprintf(format(fmt), args...)

# generate a syntactic correct unique function name in this Module
genfun() = Symbol("fmt", replace(string(gensym()), '#'=>'_'))

"""
    ALL_FORMATS::Dict{String,Function}

Global variable storing all generated format functions generated during the first
invocation of `format`, `printf`, `sprintf`, or `@format`. Not for `@printf` and
`@sprintf`.
"""
const ALL_FORMATS = Dict{String, Function}()

end # module


text:
# Interface for openBuildNet's ports

type OBNInput{T} <: OBNInputAbstract
  name::String
  node::OBNNode
  portid::Cint
  strict::Bool
  rcv_cb::OBNCallback    # Callback for message receive

  function OBNInput(name::String, node::OBNNode, portid::Cint, strict::Bool)
    obj = new(name, node, portid, strict, OBNCallback())
    node.input_ports[portid] = obj    # Add the input port object to the node, so we can refer to it later on
    finalizer(obj, x -> delete!(x.node.input_ports, x.portid))  # Remove the port from node if it's deleted
    obj
  end
end

type OBNOutput{T} <: OBNOutputAbstract
  name::String
  node::OBNNode
  portid::Cint
end

typealias OBNEI_PortType Cuint
typealias OBNEI_ContainerType Cuint
typealias OBNEI_ElementType Cuint
typealias OBNEI_FormatType Cuint

# Type to pass port's info
type OBNEI_PortInfo
  port_type::OBNEI_PortType
  container_type::OBNEI_ContainerType
  element_type::OBNEI_ElementType
  strict::Cuchar
end

const OBNEI_Port_Input = convert(OBNEI_PortType, 0)
const OBNEI_Port_Output = convert(OBNEI_PortType, 1)
const OBNEI_Port_Data = convert(OBNEI_PortType, 2)

# const OBNEI_Container_Scalar = convert(OBNEI_ContainerType, 0)
# const OBNEI_Container_Vector = convert(OBNEI_ContainerType, 1)
# const OBNEI_Container_Matrix = convert(OBNEI_ContainerType, 2)
# const OBNEI_Container_Binary = convert(OBNEI_ContainerType, 3)

elementtypeid(::Type{Bool}) = convert(OBNEI_ElementType, 0)
elementtypeid(::Type{Float64}) = convert(OBNEI_ElementType, 1)
elementtypeid(::Type{Int32}) = convert(OBNEI_ElementType, 2)
elementtypeid(::Type{Int64}) = convert(OBNEI_ElementType, 3)
elementtypeid(::Type{UInt32}) = convert(OBNEI_ElementType, 4)
elementtypeid(::Type{UInt64}) = convert(OBNEI_ElementType, 5)

function obntypeid{T}(::Type{T})
  (convert(OBNEI_ContainerType, 0), elementtypeid(T))
end

function obntypeid{T}(::Type{Vector{T}})
  (convert(OBNEI_ContainerType, 1), elementtypeid(T))
end

function obntypeid{T}(::Type{Matrix{T}})
  (convert(OBNEI_ContainerType, 2), elementtypeid(T))
end

const OBNEI_Format_ProtoBuf = convert(OBNEI_FormatType, 0)

# Create an input / output port on a node
# Currently only ProtoBuf (format = :PB) is supported
function create_input(node::OBNNode, name::String, T::DataType, strict::Bool = false, format::Symbol = :PB)
  @assert node.valid "Not allowed to create ports on non-valid node."

  if format == :PB || format == :ProtoBuf
    # Get the container and element types
    container, elementtype = obntypeid(T)
    result = ccall(_api_createInputPort, Cint, (Csize_t, Cstring, OBNEI_FormatType, OBNEI_ContainerType, OBNEI_ElementType, Cuchar),
                   node.node_id, name, OBNEI_Format_ProtoBuf, container, elementtype, strict)
    if result < 0
      error("Error creating input port [$result]: ", lastErrorMessage())
    else
      return OBNInput{T}(name, node, result, strict)
    end
  else
    error("Unsupported format: ", format)
  end
end

function create_output(node::OBNNode, name::String, T::DataType, format::Symbol = :PB)
  @assert node.valid "Not allowed to create ports on non-valid node."

  if format == :PB || format == :ProtoBuf
    # Get the container and element types
    container, elementtype = obntypeid(T)
    result = ccall(_api_createOutputPort, Cint, (Csize_t, Cstring, OBNEI_FormatType, OBNEI_ContainerType, OBNEI_ElementType),
                   node.node_id, name, OBNEI_Format_ProtoBuf, container, elementtype)
    if result < 0
      error("Error creating output port [$result]: ", lastErrorMessage())
    else
      return OBNOutput{T}(name, node, result)
    end
  else
    error("Unsupported format: ", format)
  end
end

# Read a Vector port, with possibly a default value (in case a strict port has no pending value to read from)
# function get{T}(port::OBNInput{Vector{T}}, defval = nothing)
# end

_scalar_input_api(::Type{Float64}) = _api_inputScalarDoubleGet
_scalar_input_api(::Type{Int32}) = _api_inputScalarInt32Get
_scalar_input_api(::Type{Int64}) = _api_inputScalarInt64Get
_scalar_input_api(::Type{UInt32}) = _api_inputScalarUInt32Get
_scalar_input_api(::Type{UInt64}) = _api_inputScalarUInt64Get
_scalar_input_api(::Type{Bool}) = _api_inputScalarBoolGet

_vector_input_api(::Type{Float64}) = (_api_inputVectorDoubleGet, _api_inputVectorDoubleRelease)
_vector_input_api(::Type{Int32}) = (_api_inputVectorInt32Get, _api_inputVectorInt32Release)
_vector_input_api(::Type{Int64}) = (_api_inputVectorInt64Get, _api_inputVectorInt64Release)
_vector_input_api(::Type{UInt32}) = (_api_inputVectorUInt32Get, _api_inputVectorUInt32Release)
_vector_input_api(::Type{UInt64}) = (_api_inputVectorUInt64Get, _api_inputVectorUInt64Release)
_vector_input_api(::Type{Bool}) = (_api_inputVectorBoolGet, _api_inputVectorBoolRelease)

_matrix_input_api(::Type{Float64}) = (_api_inputMatrixDoubleGet, _api_inputMatrixDoubleRelease)
_matrix_input_api(::Type{Int32}) = (_api_inputMatrixInt32Get, _api_inputMatrixInt32Release)
_matrix_input_api(::Type{Int64}) = (_api_inputMatrixInt64Get, _api_inputMatrixInt64Release)
_matrix_input_api(::Type{UInt32}) = (_api_inputMatrixUInt32Get, _api_inputMatrixUInt32Release)
_matrix_input_api(::Type{UInt64}) = (_api_inputMatrixUInt64Get, _api_inputMatrixUInt64Release)
_matrix_input_api(::Type{Bool}) = (_api_inputMatrixBoolGet, _api_inputMatrixBoolRelease)

import Base.get

# Read a scalar port, with possibly a default value (in case a strict port has no pending value to read from)
function get{T}(port::OBNInput{T}, defval = nothing)
  @assert port.node.valid

  val = Ref{T}()
  result = ccall(_scalar_input_api(T), Cint, (Csize_t, Csize_t, Ref{T}), port.node.node_id, port.portid, val)

  if result < 0
    error("Error reading from input [$result]: ", lastErrorMessage())
  end
  (result==0)?val[]:defval
end

# function get(port::OBNInput{Bool}, defval = nothing)
#   @assert port.node.valid
#
#   val = Ref{Cuchar}()
#   result = ccall(_api_inputScalarBoolGet, Cint, (Csize_t, Csize_t, Ref{Cuchar}), port.node.node_id, port.portid, val)
#
#   if result < 0
#     error("Error reading from input [$result]", lastErrorMessage())
#   end
#   (result==0)?convert(Bool,val[]):defval
# end

# Read a vector port, with possibly a default value (in case a strict port has no pending value to read from)
function get{T}(port::OBNInput{Vector{T}}, defval = nothing)
  @assert port.node.valid

  (getapi, releaseapi) = _vector_input_api(T)

  manobj = Ref{Ptr{Void}}()
  nelems = Ref{Csize_t}(0)
  result = ccall(getapi, Cint, (Csize_t, Csize_t, Ref{Ptr{Void}}, Ptr{Ptr{T}}, Ref{Csize_t}), port.node.node_id, port.portid, manobj, C_NULL, nelems)

  if result < 0
    error("Error reading from input [$result]: ", lastErrorMessage())
  end

  if result == 0
    val = Vector{T}(nelems[])
    if nelems[] > 0
      ccall(releaseapi, Void, (Ptr{Void}, Ref{T}), manobj[], val)
    else
      ccall(releaseapi, Void, (Ptr{Void}, Ptr{T}), manobj[], C_NULL)
    end
    val
  else
    defval
  end
end

# function get(port::OBNInput{Vector{Bool}}, defval = nothing)
#   @assert port.node.valid
#
#   manobj = Ref{Ptr{Void}}()
#   nelems = Ref{Csize_t}(0)
#   result = ccall(_api_inputVectorBoolGet, Cint, (Csize_t, Csize_t, Ref{Ptr{Void}}, Ptr{Ptr{Cuchar}}, Ref{Csize_t}), port.node.node_id, port.portid, manobj, C_NULL, nelems)
#
#   if result < 0 || nelems[] < 0
#     error("Error reading from input [$result]", lastErrorMessage())
#   end
#
#   if result == 0
#     val = Vector{Bool}(nelems[])
#     if nelems[] > 0
#       ccall(_api_inputVectorBoolRelease, Void, (Ptr{Void}, Ref{Cuchar}), manobj[], val)
#     else
#       ccall(_api_inputVectorBoolRelease, Void, (Ptr{Void}, Ptr{Cuchar}), manobj[], C_NULL)
#     end
#     val
#   else
#     defval
#   end
# end

# Read a matrix port, with possibly a default value (in case a strict port has no pending value to read from)
function get{T}(port::OBNInput{Matrix{T}}, defval = nothing)
  @assert port.node.valid

  (getapi, releaseapi) = _matrix_input_api(T)

  manobj = Ref{Ptr{Void}}()
  nrows = Ref{Csize_t}(0)
  ncols = Ref{Csize_t}(0)
  result = ccall(getapi, Cint, (Csize_t, Csize_t, Ref{Ptr{Void}}, Ptr{Ptr{T}}, Ref{Csize_t}, Ref{Csize_t}), port.node.node_id, port.portid, manobj, C_NULL, nrows, ncols)

  if result < 0
    error("Error reading from input [$result]: ", lastErrorMessage())
  end

  if result == 0
    val = Matrix{T}(nrows[], ncols[])
    if nrows[] > 0 && ncols[] > 0
      ccall(releaseapi, Void, (Ptr{Void}, Ref{T}), manobj[], val)
    else
      ccall(releaseapi, Void, (Ptr{Void}, Ptr{T}), manobj[], C_NULL)
    end
    val
  else
    defval
  end
end

# Ask an output port to send (Synchronously)
# Returns true if successful; check lastErrorMessage() otherwise
function sendsync(port::OBNOutputAbstract)
  @assert port.node.valid
  result = ccall(_api_outputSendSync, Cint, (Csize_t, Csize_t), port.node.node_id, port.portid)
  result == 0
end

function sendsync(port::OBNOutputAbstract, val)
  set(port, val)
  sendsync(port)
end

_scalar_output_api(::Type{Float64}) = _api_outputScalarDoubleSet
_scalar_output_api(::Type{Int32}) = _api_outputScalarInt32Set
_scalar_output_api(::Type{Int64}) = _api_outputScalarInt64Set
_scalar_output_api(::Type{UInt32}) = _api_outputScalarUInt32Set
_scalar_output_api(::Type{UInt64}) = _api_outputScalarUInt64Set
_scalar_output_api(::Type{Bool}) = _api_outputScalarBoolSet

_vector_output_api(::Type{Float64}) = _api_outputVectorDoubleSet
_vector_output_api(::Type{Int32}) = _api_outputVectorInt32Set
_vector_output_api(::Type{Int64}) = _api_outputVectorInt64Set
_vector_output_api(::Type{UInt32}) = _api_outputVectorUInt32Set
_vector_output_api(::Type{UInt64}) = _api_outputVectorUInt64Set
_vector_output_api(::Type{Bool}) = _api_outputVectorBoolSet

_matrix_output_api(::Type{Float64}) = _api_outputMatrixDoubleSet
_matrix_output_api(::Type{Int32}) = _api_outputMatrixInt32Set
_matrix_output_api(::Type{Int64}) = _api_outputMatrixInt64Set
_matrix_output_api(::Type{UInt32}) = _api_outputMatrixUInt32Set
_matrix_output_api(::Type{UInt64}) = _api_outputMatrixUInt64Set
_matrix_output_api(::Type{Bool}) = _api_outputMatrixBoolSet

# Write to a scalar port
function set{T}(port::OBNOutput{T}, val::T)
  @assert port.node.valid
  result = ccall(_scalar_output_api(T), Cint, (Csize_t, Csize_t, T), port.node.node_id, port.portid, val)

  if result < 0
    error("Error writing to output [$result]: ", lastErrorMessage())
  end
  result
end

# function set(port::OBNOutput{Bool}, val::Bool)
#   @assert port.node.valid
#   result = ccall(_api_outputScalarBoolSet, Cint, (Csize_t, Csize_t, Cuchar), port.node.node_id, port.portid, val)
#
#   if result < 0
#     error("Error writing to output [$result]", lastErrorMessage())
#   end
# end

# Write to a vector port
function set{T}(port::OBNOutput{Vector{T}}, val::Vector{T})
  @assert port.node.valid
  result = ccall(_vector_output_api(T), Cint, (Csize_t, Csize_t, Ref{T}, Csize_t), port.node.node_id, port.portid, val, length(val))

  if result < 0
    error("Error writing to output [$result]: ", lastErrorMessage())
  end
  result
end

# Write to a matrix port
function set{T}(port::OBNOutput{Matrix{T}}, val::Matrix{T})
  @assert port.node.valid
  nr, nc = size(val)
  result = ccall(_matrix_output_api(T), Cint, (Csize_t, Csize_t, Ref{T}, Csize_t, Csize_t), port.node.node_id, port.portid, val, nr, nc)

  if result < 0
    error("Error writing to output [$result]: ", lastErrorMessage())
  end
  result
end

# Query if a port has a pending value
function pending(port::OBNInputAbstract)
  @assert port.node.valid
  result = ccall(_api_inputPending, Cint, (Csize_t, Csize_t), port.node.node_id, port.portid)
  result > 0
end

# Connect a given port to an input port
# Returns true if successful; check lastErrorMessage() otherwise
function connectfrom(port::OBNInputAbstract, other::String)
  @assert port.node.valid
  result = ccall(_api_portConnect, Cint, (Csize_t, Csize_t, Cstring), port.node.node_id, port.portid, other)
  result == 0
end

# Get port's information
function portinfo(port::OBNPort)
  @assert port.node.valid

  info = Ref{OBNEI_PortInfo}()
  result = ccall(_api_portInfo, Cint, (Csize_t, Csize_t, Ref{OBNEI_PortInfo}), port.node.node_id, port.portid, info)

  if result != 0
    error("Error getting port information [$result]: ", lastErrorMessage())
  end
  info[]
end

# Callback when a message is received on an input port
# Returns true if successful
function on_receive(f::Function, port::OBNInputAbstract, optargs...)
  @assert port.node.valid
  result = ccall(_api_portEnableRcvEvent, Cint, (Csize_t, Csize_t), port.node.node_id, port.portid)
  if result == 0
    port.rcv_cb = OBNCallback( (f, optargs) )
    return true
  else
    return false
  end
end


text:
using Bio.Seq
using BioBenchmarks

chrom = first(open(datafile("chr1.fa"), FASTA)).seq
@benchmark translate(RNASequence(chrom))


text:
using Images, ImageView
img = load("sample-images/cats-3061372_640.jpg");
img_area_range = 320:640

img_area = img[:, img_area_range]
img_area = imfilter(img_area, Kernel.gaussian(5))
img[:, img_area_range] = img_area

imshow(img)

# create a border of a fixed style based on a content from an image
using Images, ImageView

border_size = 50
gaussian_kernel_value = 10

img = load("sample-images/cats-3061372_640.jpg");

# add borders
img = padarray(img, Pad(:reflect, border_size, 0))
img = parent(img) # reset the indices after using padarray

# apply blurring to top border
img_area_top = 1:border_size
img[img_area_top, :] = imfilter(img[img_area_top, :], Kernel.gaussian(gaussian_kernel_value))

# apply blurring to bottom border
img_area_bottom = size(img, 1)-border_size:size(img, 1)
img[img_area_bottom, :] = imfilter(img[img_area_bottom, :], Kernel.gaussian(gaussian_kernel_value))

imshow(img)


text:
import CPLEX
import Plots
import PolyhedralRelaxations
using JuMP

const PR = PolyhedralRelaxations

opt = CPLEX.Optimizer

# This line changes the Plots backend to GR, which speeds up plotting.
# An alternative is `Plots.plotly()`, which uses Python's plotly library as the backend.
Plots.gr()

function plot_mip(f::Function, base_partition::Vector{<:Real}; f_dash::Union{Function,Nothing}=nothing, out_file="plots/mip.pdf")
    if isa(f_dash, Nothing)
        milp_data, function_data = PR.construct_milp_relaxation(f, base_partition)
    else
        milp_data, function_data = PR.construct_milp_relaxation(f, f_dash, base_partition)
    end
    lb, ub = PR.get_variable_bounds(milp_data)
    x_min, x_max = lb[milp_data.x_index], ub[milp_data.x_index]
    @info "x bounds: $x_min, $x_max"

    # Add variables
    milp = Model(opt)
    set_optimizer_attribute(milp, "CPXPARAM_ScreenOutput", 0)
    num_vars = PR.get_num_variables(milp_data)
    binary_indices = PR.get_binary_indices(milp_data)
    @variable(milp, lb[i] <= x[i=1:num_vars] <= ub[i])
    for k in binary_indices
        set_binary(x[k])
    end
    for k in 1:num_vars
        set_name(x[k], PR.get_variable_names(milp_data)[k])
    end

    # Add constraints
    A, b = PR.get_eq_constraint_matrices(milp_data)
    @constraint(milp, A * x .== b)
    A, b = PR.get_leq_constraint_matrices(milp_data)
    @constraint(milp, A * x .<= b)

    # Add objective.
    set_objective_function(milp, x[milp_data.y_index])

    # Collect solutions from different linear objectives.
    xs = collect(x_min:0.01:x_max)
    max_ys = Float64[]
    min_ys = Float64[]

    for p in xs
        set_lower_bound(x[milp_data.x_index], p)
        set_upper_bound(x[milp_data.x_index], p)

        set_objective_sense(milp, MOI.MAX_SENSE)
        optimize!(milp)
        @assert termination_status(milp) == MOI.OPTIMAL
        yval = value.(x[milp_data.y_index])
        obj = objective_value(milp)
        @assert isapprox(yval, obj, atol = 1e-5) "obj $obj and y $yval different in max problem"
        push!(max_ys, yval)

        set_objective_sense(milp, MOI.MIN_SENSE)
        optimize!(milp)
        @assert termination_status(milp) == MOI.OPTIMAL
        yval = value.(x[milp_data.y_index])
        obj = objective_value(milp)
        @assert isapprox(yval, obj, atol = 1e-5) "obj $obj and y $yval different in min problem"
        push!(min_ys, yval)
    end

    fs = [f(x) for x in xs]
    p = Plots.plot(xs, hcat(max_ys,fs,min_ys), legend=false)
    Plots.savefig(p, out_file)
    @info "finished plotting mip"
end

function get_lp_cut(m::Float64, c::Float64, x_min::Float64, x_max::Float64, y_min::Float64,
        y_max::Float64)::Pair{Vector{Float64},Vector{Float64}}
    xs, ys = Float64[], Float64[]
    for a in [x_min, x_max]
        y = (m * a) + c
        if y >= y_max
            x = (y_max - c) / m
        elseif y >= y_min
            x = a
        else
            x = (y_min - c) / m
        end
        push!(xs, x)
        push!(ys, (m*x)+c)
    end
    return Pair(xs,ys)
end

function plot_lp(f::Function, base_partition::Vector{<:Real}; f_dash::Union{Nothing,Function}=nothing, out_file="plots/lp.pdf")
    if isa(f_dash, Nothing)
        lp_data, fn_data = PR.construct_lp_relaxation(f, base_partition)
    else
        lp_data, fn_data = PR.construct_lp_relaxation(f, f_dash, base_partition)
    end
    lb, ub = PR.get_variable_bounds(lp_data)

    x_min, x_max = lb[lp_data.x_index], ub[lp_data.x_index]
    @info "x bounds: $x_min, $x_max"

    y_min, y_max = lb[lp_data.y_index], ub[lp_data.y_index]
    @info "y bounds: $y_min, $y_max"

    xs = collect(x_min:0.01:x_max)
    fs = [f(x) for x in xs]
    p = Plots.plot(xs,fs,legend=false,color=:orange)

    lp = Model(opt)
    set_optimizer_attribute(lp, "CPXPARAM_ScreenOutput", 0)

    num_vars = PR.get_num_variables(lp_data)
    @variable(lp, lb[i] <= y[i=1:num_vars] <= ub[i])
    A, b = PR.get_eq_constraint_matrices(lp_data)
    @constraint(lp, A * y .== b)

    for k in 1:num_vars
        set_name(y[k], PR.get_variable_names(lp_data)[k])
    end

    res = π/100
    for α ∈ collect(0:res:π)
        xvar = y[lp_data.x_index]
        yvar = y[lp_data.y_index]

        if (isapprox(α, π/2, atol=1e-5) || isapprox(α, -π/2, atol=1e-5))
            continue
        end

        if isapprox(α, 0, atol=1e-5)
            @objective(lp, Max, yvar)
        # elseif (isapprox(α, π/2, atol=1e-5) || isapprox(α, -π/2, atol=1e-5))
        #     @objective(lp, Max, xvar)
        else
            m = tan(α)
            @objective(lp, Max, yvar - (xvar * tan(α)))
        end

        optimize!(lp)
        @assert termination_status(lp) == MOI.OPTIMAL
        c = objective_value(lp)

        if isapprox(α, 0, atol=1e-5)
            line_xs, line_ys = [x_min, x_max], [c, c]
        # elseif (isapprox(α, π/2, atol=1e-5) || isapprox(α, -π/2, atol=1e-5))
        #     line_xs, line_ys = [c, c], [y_min, y_max]
        else
            line_xs, line_ys = get_lp_cut(m,c,x_min,x_max,y_min,y_max)
        end

        if !isempty(line_xs)
            Plots.plot!(line_xs, line_ys, color=:blue)
        end

        set_objective_sense(lp, MOI.MIN_SENSE)
        optimize!(lp)
        @assert termination_status(lp) == MOI.OPTIMAL
        c = objective_value(lp)
        if isapprox(α, 0, atol=1e-5)
            line_xs, line_ys = [x_min, x_max], [c, c]
        # elseif (isapprox(α, π/2, atol=1e-5) || isapprox(α, -π/2, atol=1e-5))
        #     line_xs, line_ys = [c, c], [y_min, y_max]
        else
            line_xs,line_ys = get_lp_cut(m,c,x_min,x_max,y_min,y_max)
        end
        Plots.plot!(line_xs, line_ys, color=:green)
    end

    Plots.savefig(p, out_file)
    @info "finished plotting lp"
end

function generate_plots()
    # f, bp = x->x^3, collect(-1.0:0.25:1.0)
    # f, bp = sin, collect(-(2*π):π/8:(2*π))
    # f, bp = cos, collect(-(2*π):π/8:(2*π))
    # f, bp = log, collect(1.0:1.0:5.0)
    # f, bp = exp, collect(0.0:0.5:2.0)
    # f, bp = cot, collect(π/16:π/32:(π-π/16))
    @info "plotting sin(x)"
    f, bp = sin, collect(0:π/8:(2*π))
    plot_mip(f, bp, out_file="plots/sin-mip.pdf")
    plot_lp(f, bp, out_file="plots/sin-lp.pdf")
    @info "plotting x ⋅ abs(x)"
    f, bp = x -> x * abs(x), collect(-1.0:0.25:1.0)
    plot_mip(f, bp, out_file="plots/xabsx-mip.pdf", f_dash=x->2*abs(x))
    plot_lp(f, bp, out_file="plots/xabsx-lp.pdf", f_dash=x->2*abs(x))
    @info "plotting x^4 - x^3"
    f, bp = x -> x^4 - x^3, collect(-0.5:0.1:1.0)
    plot_mip(f, bp, out_file="plots/poly-mip.pdf")
    plot_lp(f, bp, out_file="plots/poly-lp.pdf")
end

max_stars_repo_name: sujeevraja/univariate-experiments

max_issues_repo_name: sujeevraja/univariate-experiments

max_forks_repo_name: sujeevraja/univariate-experiments


text:
struct Problem
    C::Vector{Int64} # Vecteur des pofits de chaque variables
    A::Union{Vector{Vector{Int64}}, Array{Int64, 2}} # Matrice creuse/Liste d'adgacence des contraintes

    #Donées annexes

	nbVar::Int64
	nbConst::Int64
	equalWeights::Bool

    listOccur::Array{Int64, 2}
    occurVar :: Vector{Int64} # Nombre de variable que chaque variable penalise
    occurConst::Vector{Int64} #Nombre de fois qu'une variable apparait  dans une contraintes
end

function Problem(fname::String = "Data/test.txt")::Problem

	    f=open(fname)

	    # lecture du nbre de contraintes (m) et de variables (n)
	    m, n = parse.(Int, split(readline(f)) )

	    # lecture des n coefficients de la fonction economique et cree le vecteur d'entiers c
	    C = parse.(Int, split(readline(f)) )

		# Savoir si tout les profits sont égaux ou non
		equalWeights::Bool = true
		for var = 1:n
			equalWeights = equalWeights && (C[var] == C[1])
		end

	    # lecture des m contraintes et reconstruction de la matrice binaire A
	    A = zeros(Int64, m, n)

		# Matrice décrivant la corrélation entre les variables
	    listOccur = zeros(Int64, n, n)

		# Liste du nombre de variable correlée à une variable
	    occurVar = zeros(Int64, n)

		# Liste du nombre de contraintes où apparaissent une variable
	    occurConst = zeros(Int64, n)

	    for i=1:m
	        # lecture du nombre d'elements non nuls sur la contrainte i (non utilise)
	        nbVarConst = parse.(Int, split(readline(f)))

	        # lecture des indices des elements non nuls sur la contrainte i
	        line = parse.(Int, split(readline(f)))

	        for var in line
	            A[i, var] = 1
	            occurConst[var] += 1
	            for var2 in line
	                listOccur[var, var2] = 1
	                occurVar[var] += 1
	            end
	        end
	    end
	    close(f)

	    return Problem(C, A, n, m, equalWeights, listOccur, occurVar, occurConst)
end

struct Solution
    x::Vector{Int64}
    z::Int64

	# Données annexes

	var1::Vector{Int64} # Indice des variables misent à 1
	var0::Vector{Int64} # Indice des variables misent à 0

	nbVarLibre::Int64
	estLibre::BitArray
end

function Solution(x::Vector{Int64}, z:: Int64, estLibre::BitArray)::Solution
	nbVar::Int64 = length(x)

	num1 = sum(x)
	var1 = Vector{Int64}(undef, num1)
	var0 = Vector{Int64}(undef, nbVar-num1)

	num0 = 0
	num1 = 0
	for ind = 1:length(x)
		if x[ind] == 0
			num0 += 1
			var0[num0] = ind
		else
			num1 += 1
			var1[num1] = ind
		end
	end

	nbVarLibre::Int64 = 0
	for ind = 1:nbVar
		if estLibre[ind] == 1
			nbVarLibre += 1
		end
	end
	return Solution(x, z, var1, var0, nbVarLibre, estLibre)
end

import Base.copy

function Base.copy(sol::Solution)::Solution
	return Solution(
		copy(x),
		z,
		copy(var1),
		copy(var0),
		nbVarLibre,
		copy(estLibre)
	)
end


text:
export BBCache, transformByDim!
import AQFED.Random: AbstractSeq, skipTo, nextn! #for vectorized way

mutable struct Entry{K,V}
    key::K
    value::V
    previous::Union{Nothing,Entry{K,V}}
    next::Union{Nothing,Entry{K,V}}
end

#LRU cache implementation (not threadsafe)
mutable struct BBCache{K,V} <: AbstractDict{K,V}
    map::Dict{K,Entry{K,V}}
    head::Union{Nothing,Entry{K,V}}
    tail::Union{Nothing,Entry{K,V}}
    maxsize::Int
    function BBCache{K,V}(size::Int) where {K,V}
        new{K,V}(Dict{K,Entry{K,V}}(), nothing, nothing, size)
    end
end

Base.haskey(cache::BBCache{K,V}, key::K) where {K,V} = haskey(cache.map, key)

function Base.getindex(cache::BBCache{K,V}, key::K) where {K,V}
    entry = cache.map[key]
    removeNode(cache, entry)
    addFirst(cache, entry)
    return entry.value
end

function Base.get(cache::BBCache{K,V}, key::K, default::Union{Nothing,V}) where {K,V}
    if haskey(cache.map, key)
        return Base.getindex(cache, key)
    else
        return default
    end
end

function Base.setindex!(cache::BBCache{K,V}, value::V, key::K) where {K,V}
    if haskey(cache.map, key)
        entry = cache.map[key]
        entry.value = value
        removeNode(cache, entry)
        addFirst(cache, entry)
    else
        newEntry = Entry(key, value, nothing, nothing)
        if length(cache.map) > cache.maxsize
            delete!(cache, cache.tail)
            addFirst(cache, newEntry)
        else
            addFirst(cache, newEntry)
        end
        cache.map[key] = newEntry
    end
end

Base.length(cache::BBCache) = Base.length(cache.map)

function Base.delete!(cache::BBCache, entry::Entry)
    delete!(cache.map, entry.key)
    removeNode(cache, entry)
end

Base.isempty(cache::BBCache) = isempty(cache.map)
function Base.sizehint!(cache::BBCache, n::Integer)
    sizehint!(cache.map, n)
    return cache
end

function addFirst(cache::BBCache, entry::Entry)
    entry.next = cache.head
    entry.previous = nothing
    if (cache.head != nothing)
        cache.head.previous = entry
    end
    cache.head = entry
    if (cache.tail == nothing)
        cache.tail = cache.head
    end
end

function removeNode(cache::BBCache, entry::Entry)
    if entry.previous != nothing
        entry.previous.next = entry.next
    else
        cache.head = entry.next
    end
    if entry.next != nothing
        entry.next.previous = entry.previous
    else
        cache.tail = entry.previous
    end
end

function Base.empty!(cache::BBCache)
    empty!(cache.map)
    return cache
end

Base.show(io::IO, cache::BBCache{K,V}) where {K,V} =
    print(io, "BBCache{$K, $V}(; maxsize = $(cache.maxsize))")

#===========================================================================
 Vectorized Brownian Bridge algorithm

 Copyright (C) 2017 Jherek Healy "Applied Quantitative Finance for Equity Derivatives"
 All right reserved.

 Permission to use, copy, modify, and distribute this software is freely
 granted, provided that this notice is preserved.
===========================================================================#
function transformRecursive(
    bb::BrownianBridgeConstruction,
    qrng::AbstractSeq,
    start::Int,
    length::Int,
    dimIndex::Int,
    cache::AbstractDict,
)
    #println("transformRecursive ",dimIndex)
    ent = get(cache, dimIndex, nothing)
    if ent != nothing
        return ent
    end #else...
    outputl = Vector{Float64}(undef, length)
    if dimIndex == Base.length(bb.stdDev)
        skipTo(qrng, 1, start)
        nextn!(qrng, 1, outputl)
        #println("transformRecursive dim=",dimIndex, "fetching sobol d=",1)
        @. outputl *= bb.stdDev[1]
        cache[dimIndex] = outputl
    else
        i = bb.pmap[dimIndex] #bijective
        j = bb.leftIndex[i]
        k = bb.rightIndex[i]
        if j != 1
            outputj = transformRecursive(bb, qrng, start, length, j - 1, cache)
            outputk = transformRecursive(bb, qrng, start, length, k, cache)
            skipTo(qrng, i, start)
            nextn!(qrng, i, outputl)
        #    println("transformRecursive dim=",dimIndex, " fetching sobol d=",i," and use points ",j," and ",k)
            wl = bb.leftWeight[i]
            wr = bb.rightWeight[i]
            sd = bb.stdDev[i]
            @. outputl = wl * outputj + wr * outputk + sd * outputl
            cache[dimIndex] = outputl
        else
            outputk = transformRecursive(bb, qrng, start, length, k, cache)
            skipTo(qrng, i, start)
            nextn!(qrng, i, outputl)
            # println("transformRecursive dim=",dimIndex, "fetching sobol d=",i," and use right point ",k)
            wr = bb.rightWeight[i]
            sd = bb.stdDev[i]
            @. outputl = wr * outputk + sd * outputl
            cache[dimIndex] = outputl
        end
    end
    return outputl
end

#Vectorized BB algorithm as described in Jherek Healy "Applied Quantitative Finance for Equity Derivatives"
function transformByDim!(
    bb::BrownianBridgeConstruction,
    qrng::AbstractSeq,
    start::Int,
    dimIndex::Int,
    output::AbstractVector{<:AbstractFloat},
    cache::AbstractDict,
)
    outputDim = transformRecursive(bb, qrng, start, length(output), dimIndex, cache)
    if dimIndex == 1
        output .= outputDim
        return
    end
    outputm = transformRecursive(bb, qrng, start, length(output), dimIndex - 1, cache)
    @. output = outputDim - outputm
end

#Vectorized BB algorithm for a d-dimensional Brownian path.
function transformByDim!(
    bb::BrownianBridgeConstruction,
    qrng::AbstractSeq,
    start::Int,
    dimIndex::Int,
    output::Array{Float64,2}, #size (n,d) retrieve column(di) = out[:,di] (column major order)
    cache::AbstractDict,
)
    outputDim = transformRecursive(bb, qrng, start, size(output), dimIndex, cache)
    if dimIndex == 1
        output .= outputDim
        return
    end
    outputm = transformRecursive(bb, qrng, start, size(output), dimIndex - 1, cache)
    @. output = outputDim - outputm
end

function transformRecursive(
    bb::BrownianBridgeConstruction,
    qrng::AbstractSeq,
    start::Int,
    osize::Tuple{Int,Int},
    dimIndex::Int,
    cache::AbstractDict,
)
    ent = get(cache, dimIndex, nothing)
    if ent != nothing
        return ent
    end #else...
    outputl = Array{Float64}(undef, osize)  #size (n,d) retrieve column(di) = out[:,di] (column major order)
    d = osize[2]
    if dimIndex == Base.length(bb.stdDev)
        @inbounds for di = 1:d
            skipTo(qrng, di, start)
            nextn!(qrng, di, @view outputl[:, di])
        end
        @. outputl *= bb.stdDev[1]
        cache[dimIndex] = outputl
    else
        i = bb.pmap[dimIndex] #bijective
        j = bb.leftIndex[i]
        k = bb.rightIndex[i]
        if j != 1
            outputj = transformRecursive(bb, qrng, start, osize, j - 1, cache)
            outputk = transformRecursive(bb, qrng, start, osize, k, cache)
            @inbounds for di = 1:d
                skipTo(qrng, d * (i - 1) + di, start) #i
                nextn!(qrng, d * (i - 1) + di, @view outputl[:, di])
            end
            wl = bb.leftWeight[i]
            wr = bb.rightWeight[i]
            sd = bb.stdDev[i]
            @. outputl = wl * outputj + wr * outputk + sd * outputl
            cache[dimIndex] = outputl
        else
            outputk = transformRecursive(bb, qrng, start, osize, k, cache)
            @inbounds for di = 1:d
                skipTo(qrng, d * (i - 1) + di, start) #i
                nextn!(qrng, d * (i - 1) + di, @view outputl[:, di])
            end
            wr = bb.rightWeight[i]
            sd = bb.stdDev[i]
            #println("sizes ",Base.size(outputl)," ",Base.size(outputk))
            @. outputl = wr * outputk + sd * outputl
            cache[dimIndex] = outputl
        end
    end
    return outputl
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module ConjugatePriors

using Compat
using PDMats
using Distributions

import Base.LinAlg: Cholesky
import Base: scale, mean

import Distributions:
    BernoulliStats,
    BinomData,
    BinomialStats,
    CategoricalData,
    CategoricalStats,
    DirichletCanon,
    ExponentialStats,
    IncompleteFormulation,
    MultinomialStats,
    MvNormalStats,
    MvNormalKnownCovStats,
    NormalStats,
    NormalKnownSigma,
    NormalKnownMu,
    NormalKnownSigmaStats,
    NormalKnownMuStats,

    add!,
    add_categorical_counts!,
    shape,
    rate,
    mode,
    rand,
    pdf,
    logpdf

include("fallbacks.jl")
include("beta_binom.jl")
include("dirichlet_multi.jl")
include("gamma_exp.jl")

include("normalgamma.jl")
include("normalinversegamma.jl")
include("normalwishart.jl")
include("normalinversewishart.jl")
include("normal.jl")
include("mvnormal.jl")

end # module


text:
function sieve(limit)

end


text:
using MonteCarloX

function test_histogram(;verbose=false)
  pass = true

  list_ints = [i for i = 1:100]
  hist = MonteCarloX.Histogram(list_ints , dx="none")

  add!(hist,1)
  add!(hist,2)

  sum_values = 0.0
  for (x, H) in hist
    sum_values += H
  end

  pass &= sum_values == 102
  pass &= hist[1] == 2
  pass &= hist[2] == 2
  pass &= hist[3] == 1
  return pass

end

function test_distribution(;verbose=false)
end


text:
abstract type EoSModel end
export EoSModel

"""
    eos(model::EoSModel, V, T, z=SA[1.0])

basic Clapeyron function, returns the total Helmholtz free energy.

# Inputs:
- `model::EoSModel` Thermodynamic model to evaluate
- `V` Total volume, in [m³]
- `T` Temperature, in [K]
- `z` mole amounts, in [mol], by default is `@SVector [1.0]`

# Outputs:

- Total Helmholtz free energy, in [J]

by default, it calls `R̄*T*∑(z)*(a_ideal(ideal_model,V,T,z) + a_res(model,V,T,z))` where `ideal_model == idealmodel(model)`, where `a_res` is the reduced residual Helmholtz energy and `a_ideal` is the reduced ideal Helmholtz energy.
You can mix and match ideal models if you provide:
- `idealmodel(model)`: extracts the ideal model from your Thermodynamic model
- `a_res(model,V,T,z)`: residual reduced Helmholtz free energy

"""
function eos(model::EoSModel, V, T, z=SA[1.0])
    return N_A*k_B*sum(z)*T * (a_ideal(idealmodel(model),V,T,z)+a_res(model,V,T,z))
end
"""
    idealmodel(model::EoSModel)

retrieves the ideal model from the input's model.

# Examples:

```julia-repl
julia> pr = PR(["water"],idealmodel=IAPWS95Ideal)
PR{IAPWS95Ideal} with 1 component:
 "water"
Contains parameters: a, b, acentricfactor, Tc, Mw
julia> Clapeyron.idealmodel(pr)
IAPWS95Ideal()
```

"""
idealmodel(model::EoSModel) = model.idealmodel

"""
    eos_res(model::EoSModel, V, T, z=SA[1.0])

basic Clapeyron function, returns the residual Helmholtz free energy.

# Inputs:
- `model::EoSModel` Thermodynamic model to evaluate
- `V` Total volume, in [m³]
- `T` Temperature, in [K]
- `z` mole amounts, in [mol], by default is `@SVector [1.0]`

# Outputs:

- Residual Helmholtz free energy, in [J]

by default, it calls `R̄*T*∑(z)*(a_res(model,V,T,z))` where `a_res` is the reduced residual Helmholtz energy.
"""
function eos_res(model::EoSModel, V, T, z=SA[1.0])
    return N_A*k_B*sum(z)*T*(a_res(model,V,T,z))
end

Base.broadcastable(model::EoSModel) = Ref(model)

"""
    @comps

This macro is an alias to

    1:length(model)

The caveat is that `model` has to exist in the local namespace.
`model` is expected to be an EoSModel type that contains the `icomponents` field.
`icomponents` is an iterator that goes through all component indices.
"""
macro comps()
    return quote
        1:length(model)
    end |> esc
end

has_sites(::Type{<:EoSModel}) = false
has_groups(::Type{<:EoSModel}) = false
has_sites(::T) where T<:EoSModel = has_sites(T)
has_groups(::T) where T<:EoSModel = has_groups(T)


text:
using Optim

function PipeFlow(D,P,L;Fluid="Water",
    ΔT=40,
     FlowType="Turbulent",
     AbsRoughness=1.5e-6)

    if Fluid == "Air"
        DynVisc = 20e-6; # Pa*seconds ***Air***
        Dens = 1.225; # ***Air***
        C = 1000;# J/(kg*k) AIR
    elseif Fluid == "Oil"
        DynVisc = 5e-3; # Pa*seconds ***Oil***
        Dens = 900; # ***Oil***
        C = 1800;# J/(kg*k) AIR
    else()
        C = 4187;# J/(kg*k) WATER
        Dens = 1000;# kg/m^3
        DynVisc = 7.97e-4; # Pa*seconds ***Water***
    end

    P = P * 6894.76 # converting from PSI to Pa
    Re = 1000 # Initial guess
    ϵ = AbsRoughness / D;
    if FlowType == "Laminar"

        f, Re, LinVel = Frict(Re, P, L, D, Dens, DynVisc)

    elseif FlowType == "Turbulent"
        f, Re, LinVel = Frict(ϵ, Re, P, L, D, Dens, DynVisc)
    end

    println("Friction factor:$(f)")

    println("Re Num:$(Re)")

    Area = π / 4 * D.^2
    MassFlow = Area * LinVel * Dens
    PowerDiss = MassFlow * C * ΔT
    println("Mass Flow:$MassFlow)")
    println("Max power:$(PowerDiss)")

    return MassFlow, LinVel, PowerDiss
end

function Frict(Re, P, L, D, Dens, DynVisc)

    f = 64 / Re
    LinVel = sqrt(2 * P / (L / D * f * Dens))

    for i = 1:10
        Re = Dens * D * LinVel / DynVisc
        f = 64 / Re
        LinVel = sqrt(2 * P / (L / D * f * Dens))

    end

    return f, Re, LinVel

end

function Frict(ϵ, Re, P, L, D, Dens, DynVisc)
        # A = -2*log((ϵ/3.7)-5.02/Re*log(ϵ/3.7+13/Re)); %"A Review of Explicit Friction Factor Equations"  Zigrang and Sylvester eq. 13
        # FFac = A.^(-2);

        # FFacEq(F)  = 1.14-2*log10(ϵ+9.35/(Re.*sqrt(F)))-1/sqrt(F)
    FFacEq1(F)  = -2 * log10((ϵ / 3.7) - 5.02 / Re * log10(ϵ / 3.7 + 13 / Re)) - 1 / √(F)
    OptimCost1(F) = abs.(FFacEq1(F))
    f = optimize(OptimCost1, 0., 1.).minimizer
    LinVel = sqrt(2 * P / (L / D * f * Dens))

    for i = 1:10
        Re = Dens * D * LinVel / DynVisc
        FFacEq(F)  = -2 * log10((ϵ / 3.7) - 5.02 / Re * log10(ϵ / 3.7 + 13 / Re)) - 1 / √(F)
        OptimCost(F) = abs.(FFacEq(F))
        f = optimize(OptimCost, 0., 1.).minimizer
        LinVel = sqrt(2 * P / (L / D * f * Dens))

    end

    return f, Re, LinVel

end


text:
using BinaryBuilder

name = "Kaleido"
version = v"0.1.0"

url_prefix = "https://github.com/plotly/Kaleido/releases/download/v$(version)/kaleido"
sources = [
    ArchiveSource("$(url_prefix)_linux_x64.zip", "6214b5e3082f315ead32133d13c3aacd385d1dab0da8eca9c6a6febf1669e9c8"; unpack_target = "x86_64-linux-gnu"),
    ArchiveSource("$(url_prefix)_linux_arm64.zip", "a66d0ad6da9edb0ea00508c2eaf79b386f6f22cb67b40812537a8aa05ac2e746"; unpack_target = "aarch64-linux-gnu"),
    ArchiveSource("$(url_prefix)_mac.zip", "4a583ee9363a9feed3ed6b7308ffeb1f11e8da57a978583bb7e9b4591dc55e38"; unpack_target = "x86_64-apple-darwin14"),
    ArchiveSource("$(url_prefix)_win_x64.zip", "58b57a973e660a4e7395e595a277492aa9008e0a25c3446656a52f7652b5f8a6"; unpack_target = "x86_64-w64-mingw32"),
    ArchiveSource("$(url_prefix)_win_x86.zip", "aa87ac1957e84f4f66f70de1c2dd9a2355ef7737cdf16b3e4125350836fe0484"; unpack_target = "i686-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
cp -r ${target}/* ${prefix}/.
chmod +x ${bindir}/*
if [[ "${target}" == *-mingw* ]]; then
    chmod +x ${prefix}/kaleido.cmd
fi
if [[ "${target}" != *-apple-darwin* ]]; then
    chmod +x ${bindir}/swiftshader/*.${dlext}
fi
LIC_DIR="${prefix}/share/licenses/${SRC_NAME}"
mkdir -p "${LIC_DIR}"
mv "${prefix}/LICENSE.txt" "${LIC_DIR}/LICENSE.txt"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"),
    Platform("aarch64", "linux"),
    Platform("x86_64", "macos"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("kaleido", :kaleido),
]

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# DataLoaders - Basic input/output operations
#
# Lucas Ondel 2020

"""
    save(path, obj)

Write `obj` to file `path` in the [BSON format](http://bsonspec.org/).
The intermediate directories are created if they do not exists.
If `path` does not end with the extension ".bson", the extension is
appended to the output path. The function returns the type of the
object saved. See [`load`](@ref) to load this file again.
"""
function save(path, obj)
    if ! endswith(path, ".bson")
        path *= ".bson"
    end
    mkpath(dirname(path))
    T = typeof(obj)
    bson(path, data = obj, type = T)
    T
end

"""
    load(path)

Load a julia object saved in `path` with the function [`save`](@ref).
If `path` does not end with thex extension ".bson", the extension is
appended to input path.
"""
function load(path)
    if ! endswith(path, ".bson")
        path *= ".bson"
    end
    dict = BSON.load(path)
    convert(dict[:type], dict[:data])
end

max_stars_repo_name: lucasondel/DataLoaders

max_issues_repo_name: lucasondel/DataLoaders

max_forks_repo_name: lucasondel/DataLoaders


text:
struct DuplicatedAliasError <: Exception
    name::String
end

function Base.showerror(io::IO, err::DuplicatedAliasError)
    print(io, "DuplicatedAliasError: duplicated alias name: ")
    return printstyled(io, err.name; color=:cyan)
end

struct InvalidKeyError <: Exception
    got::Union{Symbol, String}
    keys::Union{Vector{Symbol}, Vector{String}}
end

function Base.showerror(io::IO, err::InvalidKeyError)
    print(io, "InvalidKeyError: invalid key ")
    printstyled(io, err.got; color=:light_blue)
    print(io, ", possible keys are: ")

    if length(err.keys) > 8
        for idx in 1:8
            printstyled(io, err.keys[idx]; color=:light_blue)
            if idx != 8
                print(io, ", ")
            end
        end
        print(io, "... please check documentation for other valid keys")
    else
        for idx in eachindex(err.keys)
            printstyled(io, err.keys[idx]; color=:light_blue)
            if idx != lastindex(err.keys)
                print(io, ", ")
            end
        end
    end
    return nothing
end

"""
    DuplicatedFieldError(name, type)

A field with `name` of given option `type` is duplicated in the subfields option type.
Thus one cannot use the field keyword convention when seeing this error.
"""
struct DuplicatedFieldError <: Exception
    name::Symbol
    type
end

function Base.showerror(io::IO, err::DuplicatedFieldError)
    print(io, "DuplicatedFieldError: duplicated field ")
    printstyled(io, err.name; color=:light_blue)
    print(io, " in type ")
    printstyled(io, err.type; color=:green)
    return print(io, " and its sub-fields")
end


text:
module RecipeUtils

include("process_recipes.jl")

end # module


text:
using DiffEqJump: AbstractAggregatorAlgorithm, NullAggregator

PDMPProblem(prob,jumps::ConstantRateJump;kwargs...) = PDMPProblem(prob,JumpSet(jumps);kwargs...)
PDMPProblem(prob,jumps::VariableRateJump;kwargs...) = PDMPProblem(prob,JumpSet(jumps);kwargs...)
PDMPProblem(prob,jumps::RegularJump;kwargs...) = PDMPProblem(prob,JumpSet(jumps);kwargs...)
PDMPProblem(prob,jumps::MassActionJump;kwargs...) = PDMPProblem(prob,JumpSet(jumps);kwargs...)
PDMPProblem(prob,jumps::DiffEqJump.AbstractJump...;kwargs...) = PDMPProblem(prob,JumpSet(jumps...);kwargs...)

PDMPProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::ConstantRateJump;kwargs...) = PDMPProblem(prob,aggregator,JumpSet(jumps);kwargs...)
PDMPProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::VariableRateJump;kwargs...) = PDMPProblem(prob,aggregator,JumpSet(jumps);kwargs...)
PDMPProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::RegularJump;kwargs...) = PDMPProblem(prob,aggregator,JumpSet(jumps);kwargs...)
PDMPProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::MassActionJump;kwargs...) = PDMPProblem(prob,aggregator,JumpSet(jumps);kwargs...)
PDMPProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::DiffEqJump.AbstractJump...;kwargs...) = PDMPProblem(prob,aggregator,JumpSet(jumps...);kwargs...)
PDMPProblem(prob,jumps::JumpSet;kwargs...) = PDMPProblem(prob,NullAggregator(),jumps;kwargs...)

struct DiffeqJumpWrapper
	diffeqpb
	jumps
	u
end

# encode the vector field
function (wrap::DiffeqJumpWrapper)(ẋ, xc, xd, p, t::Float64)
	wrap.diffeqpb.f(ẋ, xc, p, t)
	nothing
end

# encode the rate function
function (wrap::DiffeqJumpWrapper)(rate, xc, xd, p, t::Float64, issum::Bool)
	for ii in eachindex(rate)
		rate[ii] = wrap.jumps.variable_jumps[ii].rate(xc, p, t)
	end
	return sum(rate)
end

# encode the jump function
function (wrap::DiffeqJumpWrapper)(xc, xd, p, t::Float64, ind_reaction::Int64)
	# this is a hack to be able to call affect! from DiffEqJump which requires an integrator as an argument. But if the type of affect! is not enforced, it should work
	@inbounds for ii=1:length(wrap.u)
		wrap.u[ii] = xc[ii]
	end
	wrap.jumps.variable_jumps[ind_reaction].affect!(wrap)
	@inbounds for ii=1:length(wrap.u)
		xc[ii] = wrap.u[ii]
	end
	nothing
end

function PDMPProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::JumpSet;
	save_positions = typeof(prob) <: DiffEqBase.AbstractDiscreteProblem ? (false,true) : (true, true), kwargs...)

	@assert isinplace(prob) "The current interface requires the ODE to be written inplace"
	@assert jumps.regular_jump == nothing
	@assert jumps.massaction_jump == nothing

	pb_wrapper = DiffeqJumpWrapper(prob, jumps, copy(prob.u0))

	# get PDMP characteristics
	F = (xdot,xc,xd,p,t::Float64) -> pb_wrapper(xdot,xc,xd,p,t)
	R = (rate,xc,xd,p,t::Float64,issum::Bool) -> pb_wrapper(rate,xc,xd,p,t,issum)
	Delta = (xc,xd,p,t::Float64,ind_reaction::Int) -> pb_wrapper(xc,xd,p,t,ind_reaction)

	xc0 = copy(prob.u0)

	xd0 = [0]
	tspan = prob.tspan
	p = prob.p

	# determine the number of reactions
	nb_reactions = length(jumps.variable_jumps)

	return PDMPProblem(F,R,Delta,nb_reactions,xc0,xd0,p,tspan)
end


text:
import Pkg

if haskey(Pkg.installed(),"MosekTools")
    using MosekTools
    const dpkf_ok = true
else
    dpkf_ok = false
    println("WARNING: MosekTools.jl not installed, the functions
    for differentially private Kalman filtering cannot be used!")
end

#=
if haskey(Pkg.installed(),"SCS")
    using SCS  # does not work with SCS currently for some reason
    const dpkf_ok = true
else
    dpkf_ok = false
    println("WARNING: SCS.jl not installed, the functions for differentially
    private Kalman filtering cannot be used!")
end
=#

using JuMP
using ControlSystems
using LinearAlgebra

"""    (D, P_val, X_val, Ω_val, M) =
            staticInputBlock_DPKF_ss(Ls, As, Cs, Winvs, Vs, Vinvs, k_priv, ρ,
                                    m=size(As,1), p=size(Cs,1), r=size(Ls,1))

Compute an input matrix D for the two-block differentially private
steady-state Kalman filter mechanism. D takes linear combinations of
individual signals before privacy-preserving Gaussian noise injection.

Note that we do not guarantee observability of the pair (A,D*C). If
observability is lost and needed, call the dfactor function with the
returned matrix M and a truncation threshold lower than 1e-6 in order
to compute a new matrix D with more rows.

Inputs: one matrix per individual

- Ls: size(r,nusers*m), or size (r,m,nusers) if want estimator of sum_i L_i x_i
- As: size (m,m,nusers)
- Cs: size (p,m,nusers)
- Winvs: size (m,m,nusers) -- inverses of process noise cov. matrices
- Vs: size (p,p,nusers) -- measurement noise cov. matrices
- Vinvs: inverses of Vs (computed externally for modularity/efficiency)
- ρ: vector of size nusers, for the adjacency definition: ||yᵢ-yᵢ'||₂ <= ρᵢ
for user i.
- k_priv: proportionality constant for Gaussian privacy-preserving noise
injection. Compute it with k_priv = gaussianMechConstant(ϵ,δ) for your
choice of ϵ, δ.
"""

function staticInputBlock_DPKF_ss(Ls, As, Cs, Winvs, Vs, Vinvs, ρ, k_priv,
                                  nusers=size(As,3), m=size(As,1),
                                  p=size(Cs,1), r=size(Ls,1))

    modl = Model(Mosek.Optimizer)
    #modl = Model(with_optimizer(Mosek.Optimizer))
    #modl = Model(with_optimizer(SCS.Optimizer))

    Nm = nusers * m
    Np = nusers * p

    A = zeros(Nm, Nm)
    C = zeros(Np, Nm)
    V = zeros(Np, Np)
    Winv = zeros(Nm, Nm)
    Vinv = zeros(Np, Np)

    for i=1:nusers
        A[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = As[:,:,i]
        C[(i-1)*p+1:i*p, (i-1)*m+1:i*m] = Cs[:,:,i]
        V[(i-1)*p+1:i*p, (i-1)*p+1:i*p] = Vs[:,:,i]
        Winv[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = Winvs[:,:,i]
        Vinv[(i-1)*p+1:i*p, (i-1)*p+1:i*p] = Vinvs[:,:,i]
    end

    if size(Ls,3) > 1
        L = zeros(r, Nm)
        for i=1:nusers
            L[:, (i-1)*m+1:i*m] = Ls[:,:,i]
        end
    else
        L = copy(Ls)
    end

    if r == 1
        @variable(modl, X >= 0)
        @objective(modl, Min, X)
    else
        @variable(modl, X[1:r, 1:r], PSD)
        @objective(modl, Min, tr(X))
    end
    if Nm == 1
        @variable(modl, Ω >= 0)
    else
        @variable(modl, Ω[1:Nm, 1:Nm], PSD)
    end
    if Np == 1
        @variable(modl, P >= 0)
    else
        @variable(modl, P[1:Np, 1:Np], PSD)
    end

    @SDconstraint(modl, hvcat((2,2), X, L, L', Ω) >= zeros(r+Nm,r+Nm))

    @SDconstraint(modl, hvcat((2,2),
        C'*P*C-Ω+Winv, Winv*A, A'*Winv, Ω+A'*Winv*A) >= zeros(2*Nm,2*Nm))

    for i=1:nusers
        e = zeros(p,Np); e[:,(i-1)*p+1:i*p]=Matrix{Float64}(I,p,p)
        @SDconstraint(modl,
          hvcat((2,2), Matrix{Float64}(I,p,p)/(k_priv^2*ρ[i]^2)+Vinvs[:,:,i], e, e', V-V*P*V)
          >= zeros(Np+p,Np+p))
    end

    optimize!(modl)
    status = termination_status(modl)

    println("============================================")
    println("Solution status: ", status)
    println("Objective value: ", objective_value(modl))
    println("============================================")

    P_val = value.(P)

    X_val = value.(X)
    Ω_val = value.(Ω)

    #= # Debug
    Σ_val = inv(Ω_val)
    println("Norm of the difference between Σ^{-1} and (A Σ A' + W)^{-1}+C' P C")
    E1 = C'*P_val*C+inv(A*Σ_val*A'+inv(Winv))-inv(Σ_val)
    println(norm(E1))
    =#

    # Matrix to factorize
    M = k_priv^2 * (inv(V-V*P_val*V)-Vinv)
    # Compute D matrix
    D = dfactor(M; svaltol=1e-6)

    # We return M also, so that the factorization can be redone if needed
    return (D, P_val, X_val, Ω_val, M)
end

"""    (D, P_val, X_val, Ω_val, M) =
            staticInputBlock_DPLQG_ss(Ls, As, Bs, Cs, Vs, Vinvs, Winvs, k_priv, ρ,
                                    m=size(As,1), p=size(Cs,1), r=size(Ls,1))

Compute an input matrix D for the two-block differentially private
steady-state LQG mechanism. D takes linear combinations of
individual signals before privacy-preserving Gaussian noise injection.

Note that we do not guarantee observability of the pair (A,D*C). If
observability is lost and needed, call the dfactor function with the
returned matrix M and a truncation threshold lower than 1e-6 in order
to compute a new matrix D with more rows.

Inputs: one matrix per individual

- Q: size(nusers*m,nusers*m) (m = dim of individual state-space)
- R: size(du,du) (du = dim of control input signal)
- As: size (m,m,nusers)
- Bs: size (m,du,nsusers) - All B matrices pre-multiply the same input u
- Cs: size (p,m,nusers)
- Vs: size (p,p,nusers)
- Vinvs: inverses of Vs (computed externally for modularity/efficiency)
- Winvs: size (m,m,nusers)
- ρ: vector of size nusers, for the adjacency definition: ||yᵢ-yᵢ'||₂ <= ρᵢ
for user i.
- k_priv: proportionality constant for Gaussian privacy-preserving noise
injection. Compute it with k_priv = gaussianMechConstant(ϵ,δ) for your
choice of ϵ, δ.
"""
function staticInputBlock_DPLQG_ss(Q, R, As, Bs, Cs, Ws, Winvs, Vs, Vinvs,
                                ρ, k_priv, nusers=size(As,3), m=size(As,1),
                                p=size(Cs,1), du=size(Bs,2))

    Nm = nusers * m

    A = zeros(Nm, Nm)
    B = zeros(Nm, du)  # the same control input is shared by all agents
    W = zeros(Nm, Nm)

    for i=1:nusers
        A[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = As[:,:,i]
        B[(i-1)*m+1:i*m, :] = Bs[:,:,i]  # stack the Bs in one column
        W[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = Ws[:,:,i]
    end

    P = dare(A, B, Q, R)  #  computes ss LQR cost
    N = A'*P*A+Q-P
    L = dfactor(N, svaltol=1e-5)  # factorize L' L = P
    r = size(L,1)

    (D, perf_val, X_val, Ω_val, M) =
        staticInputBlock_DPKF_ss(L, As, Cs, Winvs, Vs, Vinvs, ρ, k_priv, nusers, m, p, r)
    return (D, tr(P*W)+tr(X_val), P, X_val, M)
end

"""    dfactor(M; svaltol=1e-4)

Compute D such that D'D = M for M positive semidefinite, via SVD. Tries to
output a wide D matrix by neglecting the singular values of M smaller than
svaltol*(max svals).
"""
function dfactor(M; svaltol=1e-4)
    (U,S,V) = svd(M)
    svalmax = maximum(S)
    threshold = svaltol*svalmax
    cols = []
    for i=1:length(S)
        if S[i] >= threshold
            push!(cols, i)
        end
    end
    D = zeros(size(M,1), length(cols))
    for j=1:length(cols)
        D[:,j] = U[:,cols[j]] * sqrt(S[cols[j]])
    end
    return sign(D[1,1])*D'
end

"""    evaluateKFperf(D, Ls, As, Cs, Vs, Ws, ρ, k_priv)

Compute the steady-state MSE of the two-block differentially private
Kalman filter, with a static matrix D to combine input signals.
The performance is computed by solving an algebraic Riccati equation.
It corresponds to the MSE at the end of a measurement update state in
the Kalman filter.

Inputs:
- D: size (q,nusers*p)  (q can be arbitrary, p = individual output signal dim.)
one matrix per individual
- Ls: size (r,m,nusers)  (if want estimator of sum_i L_i x_i)
- As: size (m,m,nusers)
- Cs: size (p,m,nusers)
- Ws: size (m,m,nusers) -- process noise cov. matrix
- Vs: size (p,p,nusers) -- measurement noise cov. matrix
- ρ: vector of size nusers, for the adjacency definition: ||yᵢ-yᵢ'||₂ <= ρᵢ
for user i.
- k_priv: proportionality constant for Gaussian privacy-preserving noise
injection. Compute it with k_priv = gaussianMechConstant(ϵ,δ) for your
choice of ϵ, δ.

Outputs: (performance, Σ, Σupdate)
where Σ and Σupdate are the error covariance matrices after the
prediction and measurement update step respectively
"""
function evaluateKFperf(D, Ls, As, Cs, Ws, Vs, ρ, k_priv)
    nusers=size(As,3)
    m=size(As,1)
    p=size(Cs,1)
    r=size(Ls,1)

    Nm = nusers * m
    Np = nusers * p

    L = zeros(r, Nm)
    A = zeros(Nm, Nm)
    C = zeros(Np, Nm)
    W = zeros(Nm, Nm)
    V = zeros(Np, Np)

    for i=1:nusers
        L[:, (i-1)*m+1:i*m] = Ls[:,:,i]
        A[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = As[:,:,i]
        C[(i-1)*p+1:i*p, (i-1)*m+1:i*m] = Cs[:,:,i]
        W[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = Ws[:,:,i]
        V[(i-1)*p+1:i*p, (i-1)*p+1:i*p] = Vs[:,:,i]
        #Winv[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = Winvs[:,:,i]
        #Vinv[(i-1)*p+1:i*p, (i-1)*p+1:i*p] = Vinvs[:,:,i]
    end

    # compute sensitivity
    Δ₂ = 0
    for i=1:nusers
        Δ₂ = max(Δ₂, ρ[i] * maximum(svdvals(D[:,(i-1)*p+1:i*p])))
    end

    V₁ = D*V*D' + k_priv^2  * Δ₂^2 * Matrix{Float64}(I,size(D,1),size(D,1))
    C₁ = D*C
    V₁ = 0.5*(V₁+V₁')  # might have lost perfect symmetry due to numerical issues
    Σ = dare(A', C₁', W, V₁)  #  computes ss cov. after time update step
    Σupdate = Σ - Σ*C₁'*inv(C₁*Σ*C₁'+V₁)*C₁*Σ  # cov. afer meas. update step

    #return (trace(L*Σupdate*L'), Δ₂, Σupdate)
    return tr(L*Σupdate*L'), Σ, Σupdate
end

"""    evaluateLQGperf(D, Q, R, As, Bs, Cs, Ws, Vs, ρ, k_priv)

Compute the steady-state MSE of the two-block differentially private
LQG controller, with a static matrix D to combine input signals.
The performance is computed by solving two algebraic Riccati equations.

Inputs:
- D: size (q,nusers*p)  (q can be arbitrary, p = individual output signal dim.)
one matrix per individual
- Q: size(nusers*m,nusers*m) (m = dim of individual state-space)
- R: size(du,du) (du = dim of control input signal)
- As: size (m,m,nusers)
- Bs: size (m,du,nsusers) - All B matrices pre-multiply the same input u
- Cs: size (p,m,nusers)
- Ws: size (m,m,nusers) -- Process noise
- Vs: size (p,p,nusers) -- Measurement noise
- ρ: vector of size nusers, for the adjacency definition: ||yᵢ-yᵢ'||₂ <= ρᵢ
for user i.
- k_priv: proportionality constant for Gaussian privacy-preserving noise
injection. Compute it with k_priv = gaussianMechConstant(ϵ,δ) for your
choice of ϵ, δ.

Outputs: (performance, Σ, Σupdate, P)
where Σ and Σupdate are the error covariance matrices after the
prediction and measurement update step respectively, P is the solution
of the control algebraic Riccati equation
"""
function evaluateLQGperf(D, Q, R, As, Bs, Cs, Ws, Vs, ρ, k_priv)
    nusers=size(As,3)
    m=size(As,1)
    p=size(Cs,1)
    du=size(R,1)

    Nm = nusers * m
    Np = nusers * p

    A = zeros(Nm, Nm)
    B = zeros(Nm, du)  # the same control input is shared by all agents
    C = zeros(Np, Nm)
    W = zeros(Nm, Nm)
    V = zeros(Np, Np)

    for i=1:nusers
        A[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = As[:,:,i]
        B[(i-1)*m+1:i*m, :] = Bs[:,:,i]  # stack the Bs in one column
        C[(i-1)*p+1:i*p, (i-1)*m+1:i*m] = Cs[:,:,i]
        W[(i-1)*m+1:i*m, (i-1)*m+1:i*m] = Ws[:,:,i]
        V[(i-1)*p+1:i*p, (i-1)*p+1:i*p] = Vs[:,:,i]
    end

    P = dare(A, B, Q, R)  #  computes ss LQR cost matrix

    # compute sensitivity
    Δ₂ = 0
    for i=1:nusers
        Δ₂ = max(Δ₂, ρ[i] * maximum(svdvals(D[:,(i-1)*p+1:i*p])))
    end

    V₁ = D*V*D' + k_priv^2  * Δ₂^2 * Matrix{Float64}(I,size(D,1),size(D,1))
    C₁ = D*C
    V₁ = 0.5*(V₁+V₁')  # might have lost perfect symmetry due to numerical issues
    Σ = dare(A', C₁', W, V₁)  #  computes ss cov. after time update step
    Σupdate = Σ - Σ*C₁'*inv(C₁*Σ*C₁'+V₁)*C₁*Σ  # cov. afer meas. update step

    N = Q+A'*P*A-P
    return tr(P*W + N*Σupdate), P, Σ, Σupdate
end


text:
# Test if elements in the sparse array (S) are present
# in the dense array (D)
function SinD(S, D, verbose = false)
    Si = S[1]
    Sx = S[2]
    Sy = S[3]
    Sz = S[4]

    Dxs = D[:,:,:,:,1]
    Dys = D[:,:,:,:,2]
    Dzs = D[:,:,:,:,3]
    for i in eachindex(Si)
        idx = Si[i] .+ 1
        x = Sx[i]
        y = Sy[i]
        z = Sz[i]
        Dx = Dxs[idx...]
        Dy = Dys[idx...]
        Dz = Dzs[idx...]

        for (k,Dk) in [(x,Dx), (y,Dy), (z,Dz)]
            if abs(k - Dk) > 1e-10
                println("$idx diff found: $k against $Dk")
                return false
            end
        end
        verbose ? println("$idx OK") : nothing
    end
    return true
end

atoms = Molecules.parse_string("""
C   -2.131551124300    2.286168823700    0.000000000000
H   -1.061551124300    2.286168823700    0.000000000000
H   -2.488213906200    1.408104616400    0.496683911300
H   -2.488218762100    2.295059432700   -1.008766153900
H   -2.488220057000    3.155340844300    0.512081313000""")

bs = BasisSet("cc-pvdz", atoms)

@testset "∂S/∂X" begin
    for iA = 1:5
        dS = ∇overlap(bs, iA)
        for k = 1:3
            @test dS[:,:,k] ≈ ∇FD_overlap(bs, iA, k)
        end
    end
end

@testset "∂T/∂X" begin
    for iA = 1:5
        dT = ∇kinetic(bs, iA)
        for k = 1:3
            @test dT[:,:,k] ≈ ∇FD_kinetic(bs, iA, k)
        end
    end
end

@testset "∂V/∂X" begin
    for iA = 1:5
        dV = ∇nuclear(bs, iA)
        for k = 1:3
            @test dV[:,:,k] ≈ ∇FD_nuclear(bs, iA, k)
        end
    end
end

@testset "∂[ij|kl]/∂X" begin
    @testset "Dense" begin
        for iA = 1:5
            dERI = ∇ERI_2e4c(bs, iA)
            for k = 1:3
                @test dERI[:,:,:,:,k] ≈ ∇FD_ERI_2e4c(bs, iA, k)
            end
        end
    end

    @testset "Sparse" begin
        for iA = 1:5
            sparse = ∇sparseERI_2e4c(bs, iA)
            dense = ∇ERI_2e4c(bs, iA)
            @test SinD(sparse, dense)
        end
    end
end