text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GTK3"
version = v"3.24.29"

# Collection of sources required to build GTK
sources = [
    ArchiveSource("http://ftp.gnome.org/pub/gnome/sources/gtk+/$(version.major).$(version.minor)/gtk+-$(version).tar.xz",
                  "f57ec4ade8f15cab0c23a80dcaee85b876e70a8823d9105f067ce335a8268caa"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gtk+-*/

# We need to run some commands with a native Glib
apk add glib-dev gtk+3.0

if [[ "${target}" == *-linux-* ]]; then
    # We need to run `wayland-scanner` on the build system only when Wayland is
    # available, i.e. only on Linux
    apk add wayland-dev
fi

# This is awful, I know
ln -sf /usr/bin/glib-compile-resources ${prefix}/bin/glib-compile-resources
ln -sf /usr/bin/glib-compile-schemas ${prefix}/bin/glib-compile-schemas
ln -sf /usr/bin/gdk-pixbuf-pixdata ${prefix}/bin/gdk-pixbuf-pixdata
# Remove wayland-scanner when present, so that we can call the native one
rm -f ${prefix}/bin/wayland-scanner

# Patch upstreamed: https://gitlab.gnome.org/GNOME/gtk/-/merge_requests/3526
atomic_patch -p1 $WORKSPACE/srcdir/patches/gdkwindow-quartz_c.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/meson_build.patch
# Patch upstreamed: https://gitlab.gnome.org/GNOME/gtk/-/merge_requests/3527
atomic_patch -p1 $WORKSPACE/srcdir/patches/c_standard_c99.patch

FLAGS=()
if [[ "${target}" == *-apple-* ]]; then
    FLAGS+=(-Dx11_backend=false -Dwayland_backend=false)
elif [[ "${target}" == *-freebsd* ]]; then
    FLAGS+=(-Dwayland_backend=false)
fi

if [[ "${target}" == *-mingw* ]]; then
    # Despite what's declared in their meson.build, they require Meson 0.54:
    # https://gitlab.gnome.org/GNOME/gtk/-/issues/3927
    apk add meson --repository=http://dl-cdn.alpinelinux.org/alpine/edge/testing
fi

mkdir build-gtk && cd build-gtk
meson .. \
    -Dintrospection=false \
    -Ddemos=false \
    -Dexamples=false \
    -Dtests=false \
    "${FLAGS[@]}" \
    --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install

# Remove temporary links
rm ${prefix}/bin/gdk-pixbuf-pixdata ${prefix}/bin/glib-compile-{resources,schemas}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libgailutil-3", :libgailutil3),
    LibraryProduct("libgdk-3", :libgdk3),
    LibraryProduct("libgtk-3", :libgtk3),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Glib_jll", v"2.59.0"; compat="2.59.0"),
    Dependency("Cairo_jll"),
    Dependency("Pango_jll", v"1.42.4"; compat="1.42.4"),
    Dependency("FriBidi_jll"),
    Dependency("FreeType2_jll"),
    Dependency("gdk_pixbuf_jll"),
    Dependency("Libepoxy_jll"),
    # Gtk 3.24.29 requires ATK 2.35.1
    Dependency("ATK_jll", v"2.36"; compat="2.35.1"),
    Dependency("HarfBuzz_jll"),
    Dependency("xkbcommon_jll"),
    Dependency("iso_codes_jll"),
    Dependency("Wayland_jll"),
    Dependency("Xorg_libXrandr_jll"),
    Dependency("Xorg_libX11_jll"),
    Dependency("Xorg_libXrender_jll"),
    Dependency("Xorg_libXi_jll"),
    Dependency("Xorg_libXext_jll"),
    Dependency("Xorg_libXcursor_jll"),
    Dependency("Xorg_libXdamage_jll"),
    Dependency("Xorg_libXfixes_jll"),
    Dependency("Xorg_libXcomposite_jll"),
    Dependency("Xorg_libXinerama_jll"),
    Dependency("Fontconfig_jll"),
    Dependency("at_spi2_atk_jll"),
]

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


text:
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/JuliaOpt/ECOS.jl
#############################################################################
# ECOSSolverInterface.jl
# MathProgBase.jl interface for the ECOS.jl solver wrapper
#############################################################################

require(joinpath(Pkg.dir("MathProgBase"),"src","MathProgSolverInterface.jl"))
importall MathProgSolverInterface

#############################################################################
# Define the MPB Solver and Model objects
export ECOSSolver
immutable ECOSSolver <: AbstractMathProgSolver
end

type ECOSMathProgModel <: AbstractMathProgModel
    nvar::Int                           # Number of variables
    nineq::Int                          # Number of inequalities Gx <=_K h
    neq::Int                            # Number of equalities Ax = b
    npos::Int                           # Number of ???
    ncones::Int                         # Number of SO cones
    conedims::Vector{Int}               # ?
    G::SparseMatrixCSC{Float64,Int}     # The G matrix (inequalties)
    A::SparseMatrixCSC{Float64,Int}     # The A matrix (equalities)
    c::Vector{Float64}                  # The objective coeffs (always min)
    orig_sense::Symbol                  # Original objective sense
    h::Vector{Float64}                  # RHS for inequality
    b::Vector{Float64}                  # RHS for equality
    # Post-solve
    solve_stat::Symbol
    obj_val::Float64
    primal_sol::Vector{Float64}
    fwd_map::Vector{Int}                # To reorder solution if we solved
end                                     # using the conic interface
ECOSMathProgModel() = ECOSMathProgModel(0,0,0,0,0,
                                        Int[],
                                        spzeros(0,0),
                                        spzeros(0,0),
                                        Float64[], :Min,
                                        Float64[], Float64[],
                                        :NotSolved, 0.0, Float64[], Int[])

#############################################################################
# Begin implementation of the MPB low-level interface
# Implements
# - model
# - loadproblem!
# - optimize!
# - status
# http://mathprogbasejl.readthedocs.org/en/latest/lowlevel.html

model(s::ECOSSolver) = ECOSMathProgModel()

# Loads the provided problem data to set up the linear programming problem:
# min c'x
# st  lb <= Ax <= ub
#      l <=  x <= u
# where sense = :Min or :Max
function loadproblem!(m::ECOSMathProgModel, A, collb, colub, obj, rowlb, rowub, sense)
    (nvar = length(collb)) == length(colub) || error("Unequal lengths for column bounds")
    (nrow = length(rowlb)) == length(rowub) || error("Unequal lengths for row bounds")

    # Turn variable bounds into constraints
    # Inefficient, because keeps allocating memory!
    # Would need to batch, get tricky...
    for j = 1:nvar
        if collb[j] != -Inf
            # Variable has lower bound
            newrow = zeros(1, nvar)
            newrow[j] = -1.0
            A = vcat(A, newrow)
            rowlb = vcat(rowlb, -Inf)
            rowub = vcat(rowub, -collb[j])
            nrow += 1
        end
        if colub[j] != +Inf
            # Variable has upper bound
            newrow = zeros(1, nvar)
            newrow[j] = 1.0
            A = vcat(A, newrow)
            rowlb = vcat(rowlb, -Inf)
            rowub = vcat(rowub, colub[j])
            nrow += 1
        end
    end

    eqidx   = Int[]      # Equality row indicies
    ineqidx = Int[]      # Inequality row indicies
    eqbnd   = Float64[]  # Bounds for equality rows
    ineqbnd = Float64[]  # Bounds for inequality row
    for it in 1:nrow
        # Equality constraint
        if rowlb[it] == rowub[it]
            push!(eqidx, it)
            push!(eqbnd, rowlb[it])
        # Range constraint - not supported
        elseif rowlb[it] != -Inf && rowub[it] != Inf
            error("Ranged constraints unsupported!")
        # Less-than constraint
        elseif rowlb[it] == -Inf
            push!(ineqidx, it)
            push!(ineqbnd, rowub[it])
        # Greater-than constraint - flip sign so only have <= constraints
        else
            push!(ineqidx, it)
            push!(ineqbnd, -rowlb[it])
            A[it,:] *= -1 # flip signs so we have Ax<=b
        end
    end

    m.nvar      = nvar                  # Number of variables
    m.nineq     = length(ineqidx)       # Number of inequalities Gx <=_K h
    m.neq       = length(eqidx)         # Number of equalities Ax = b
    m.npos      = length(ineqidx)       # Number of ???
    m.ncones    = 0                     # Number of SO cones
    m.conedims  = Int[]                 # ???
    m.G         = sparse(A[ineqidx,:])  # The G matrix (inequalties)
    m.A         = sparse(A[eqidx,:])    # The A matrix (equalities)
    m.c         = (sense == :Max) ? obj * -1 : obj[:]
                                        # The objective coeffs (always min)
    m.orig_sense = sense                # Original objective sense
    m.h         = ineqbnd               # RHS for inequality
    m.b         = eqbnd                 # RHS for equality
    m.fwd_map   = [1:nvar]              # Identity mapping
end

function optimize!(m::ECOSMathProgModel)
    ecos_prob_ptr = setup(
        m.nvar, m.nineq, m.neq,
        m.npos, m.ncones, m.conedims,
        m.G, m.A,
        m.c[:],  # Seems to modify this
        m.h, m.b)

    flag = solve(ecos_prob_ptr)
    if flag == ECOS_OPTIMAL
        m.solve_stat = :Optimal
    elseif flag == ECOS_PINF
        m.solve_stat = :Infeasible
    elseif flag == ECOS_DINF  # Dual infeasible = primal unbounded, probably
        m.solve_stat = :Unbounded
    elseif flag == ECOS_MAXIT
        m.solve_stat = :UserLimit
    else
        m.solve_stat = :Error
    end
    # Extract solution
    ecos_prob = pointer_to_array(ecos_prob_ptr, 1)[1]
    m.primal_sol = pointer_to_array(ecos_prob.x, m.nvar)[:]
    m.obj_val = dot(m.c, m.primal_sol) * (m.orig_sense == :Max ? -1 : +1)
    cleanup(ecos_prob_ptr, 0)
end

status(m::ECOSMathProgModel) = m.solve_stat
getobjval(m::ECOSMathProgModel) = m.obj_val
getsolution(m::ECOSMathProgModel) = m.primal_sol[m.fwd_map]

#############################################################################
# Begin implementation of the MPB conic interface
# Implements
# - loadconicproblem!
# http://mathprogbasejl.readthedocs.org/en/latest/conic.html

function loadconicproblem!(m::ECOSMathProgModel, c, A, b, cones)
    # TODO (if it matters): make this more efficient for sparse A

    # We don't support SOCRotated, SDP, or Exp*
    bad_cones = [:SOCRotated, :SDP, :ExpPrimal, :ExpDual]
    for cone_vars in cones
        cone_vars[1] in bad_cones && error("Cone type $(cone_vars[1]) not supported")
    end

    # MathProgBase form             ECOS form
    # min c'x                       min c'x
    # st  A x = b                   st  A x = b
    #       x in K                      h - Gx in K

    # Expand out the cones info
    # The cones can come in any order, so we need to build a mapping
    # from the variable indices in the input to the internal ordering
    # we will use.

    # In the first past we'll just count up the number of variables
    # of each type.
    num_vars = 0
    for (cone_type, idxs) in cones
        num_vars += length(idxs)
    end
    fwd_map = Array(Int,    num_vars)  # Will be used for SOCs
    rev_map = Array(Int,    num_vars)  # Need to restore sol. vec.
    idxcone = Array(Symbol, num_vars)  # We'll uses this for non-SOCs

    # Now build the mapping
    pos = 1
    for (cone, idxs) in cones
        for i in idxs
            fwd_map[i]   = pos   # fwd_map = orig idx -> internal idx
            rev_map[pos] = i     # rev_map = internal idx -> orig idx
            idxcone[pos] = cone
            pos += 1
        end
    end

    # Rearrange data into the internal ordering
    ecos_c = c[rev_map]
    ecos_A = A[:,rev_map]
    ecos_b = b[:]

    # For all variables in the :Zero cone, fix at 0 with an
    # equality constraint. TODO: Don't even include them

    for j = 1:num_vars
        idxcone[j] != :Zero && continue

        new_row    = zeros(1,num_vars)
        new_row[j] = 1.0
        ecos_A     = vcat(ecos_A, new_row)
        ecos_b     = vcat(ecos_b, 0.0)
    end

    # Build G matrix
    # There will be one row for every :NonNeg and :NonPos cone
    # and an additional row for every variable in a :SOC cone
    # Or in other words, everything that isn't a :Free or :Zero
    # gets a row in G and h
    num_G_row = 0
    for j = 1:num_vars
        idxcone[j] == :Free && continue
        idxcone[j] == :Zero && continue
        num_G_row += 1
    end
    ecos_G = zeros(num_G_row,num_vars)
    ecos_h = zeros(num_G_row)

    # First, handle the :NonNeg, :NonPos cases
    num_pos_orth = 0
    G_row = 1
    for j = 1:num_vars
        if idxcone[j] == :NonNeg
            ecos_G[G_row,j] = -1.0
            G_row += 1
            num_pos_orth += 1
        elseif idxcone[j] == :NonPos
            ecos_G[G_row,j] = +1.0
            G_row += 1
            num_pos_orth += 1
        end
    end
    @assert G_row == num_pos_orth + 1
    # Now handle the SOCs
    # The MPB unput form is basically just says a vector of
    # variables (y,x) lives in the SOC  || x || <= y
    # ECOS wants somethings in the form h - Gx in Q so we
    # will prove 0 - Ix \in Q
    num_SOC_cones = 0
    SOC_conedims  = Int[]
    for (cone, idxs) in cones
        cone != :SOC && continue
        # Found a new SOC
        num_SOC_cones += 1
        push!(SOC_conedims, length(idxs))
        # Add the entries (carrying on from pos. orthant)
        for j in idxs
            ecos_G[G_row,fwd_map[j]] = -1.0
            G_row += 1
        end
    end
    @assert G_row == num_G_row + 1

    # Store in the ECOS structure
    m.nvar          = num_vars
    m.nineq         = num_G_row
    m.neq           = length(ecos_b)
    m.npos          = num_pos_orth
    m.ncones        = num_SOC_cones
    m.conedims      = SOC_conedims
    m.G             = ecos_G
    m.A             = ecos_A
    m.c             = ecos_c
    m.orig_sense    = :Min
    m.h             = ecos_h
    m.b             = ecos_b
    m.fwd_map       = fwd_map
end

ext: jl


text:
module PerfMon

using ..LIKWID:
    LibLikwid,
    perfmon_initialized,
    topo_initialized,
    numa_initialized,
    timer_initialized,
    init_topology,
    init_numa,
    GroupInfoCompact
using OrderedCollections

function init(cpus::AbstractVector{Int32}=[Int32(i - 1) for i in 1:Threads.nthreads()])
    perfmon_initialized[] && finalize()

    if !topo_initialized[]
        init_topology() || error("Couldn't init topology.")
    end
    if !numa_initialized[]
        init_numa() || error("Couldn't init numa.")
    end

    nthreads = length(cpus)
    ret = LibLikwid.perfmon_init(nthreads, cpus)
    if ret == 0
        perfmon_initialized[] = true
        timer_initialized[] = true
        return true
    end
    return false
end

init(cpus::AbstractVector{<:Integer}) = init(convert(Vector{Int32}, cpus))
init(cpu::Integer) = init(Int32[cpu])

function finalize()
    LibLikwid.perfmon_finalize()
    perfmon_initialized[] = false
    return nothing
end

_check_groupid(gid) = 0 ≤ gid < get_number_of_groups()
_check_eventidx(gid, eidx) = 0 ≤ eidx < get_number_of_events(gid)
_check_metricidx(gid, eidx) = 0 ≤ eidx < get_number_of_metrics(gid)
_check_threadidx(tidx) = 0 ≤ tidx < get_number_of_threads()

"""
Return the number of threads initialized in the perfmon module.
"""
get_number_of_threads() = LibLikwid.perfmon_getNumberOfThreads()

"""
Return the number of groups currently registered in the perfmon module.
"""
get_number_of_groups() = LibLikwid.perfmon_getNumberOfGroups()

"""
Return the amount of events in the group.
"""
get_number_of_events(groupid::Integer) = LibLikwid.perfmon_getNumberOfEvents(groupid)

"""
Return the amount of metrics in the group.
Always zero for custom event sets.
"""
get_number_of_metrics(groupid::Integer) = LibLikwid.perfmon_getNumberOfMetrics(groupid)

"""
Return a list of all available perfmon groups.

# Examples
```jldoctest
julia> LIKWID.PerfMon.get_groups()
30-element Vector{LIKWID.GroupInfoCompact}:
 MEM_SP => Overview of arithmetic and main memory performance
 CYCLE_ACTIVITY => Cycle Activities
 ENERGY => Power and Energy consumption
 UOPS_RETIRE => UOPs retirement
 PMM => Intel Optane DC bandwidth in MBytes/s
 TLB_INSTR => L1 Instruction TLB miss rate/ratio
 DATA => Load to store ratio
 UOPS_ISSUE => UOPs issueing
 L2CACHE => L2 cache miss rate/ratio
 L2 => L2 cache bandwidth in MBytes/s
 ⋮
 TMA => Top down cycle allocation
 CLOCK => Power and Energy consumption
 FLOPS_DP => Double Precision MFLOP/s
 CYCLE_STALLS => Cycle Activities (Stalls)
 L3 => L3 cache bandwidth in MBytes/s
 UPI => UPI traffic
 L3NEW => L3 cache bandwidth in MBytes/s
 L3PF => L3 cache bandwidth in MBytes/s
 L2L3 => L3 cache bandwidth in MBytes/s
```
"""
function get_groups()
    if !topo_initialized[]
        init_topology() || error("Couldn't init topology.")
    end
    # refs to char**
    groups_ref = Ref{Ptr{Ptr{Cchar}}}()
    shorts_ref = Ref{Ptr{Ptr{Cchar}}}()
    longs_ref = Ref{Ptr{Ptr{Cchar}}}()

    ret = LibLikwid.perfmon_getGroups(groups_ref, shorts_ref, longs_ref)
    ret <= 0 && return nothing

    groups_vec = unsafe_wrap(Array, groups_ref[], ret)
    shorts_vec = unsafe_wrap(Array, shorts_ref[], ret)
    longs_vec = unsafe_wrap(Array, longs_ref[], ret)
    res = Vector{GroupInfoCompact}(undef, ret)
    for i in 1:ret
        res[i] = GroupInfoCompact(
            unsafe_string(groups_vec[i]),
            unsafe_string(shorts_vec[i]),
            unsafe_string(longs_vec[i]),
        )
    end
    LibLikwid.perfmon_returnGroups(ret, groups_ref[], shorts_ref[], longs_ref[])
    return res
end

"""
    add_event_set(estr) -> groupid
Add a performance group or a custom event set to the perfmon module.
Returns a `groupid` which is required to later specify the event set.
"""
function add_event_set(estr::AbstractString)
    perfmon_initialized[] || return nothing
    groupid = LibLikwid.perfmon_addEventSet(estr)
    return Int(groupid)
end

"""
Return the name of the group identified by `groupid`.
If it is a custom event set, the name is set to `Custom`.
"""
function get_name_of_group(groupid::Integer)
    perfmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    name = unsafe_string(LibLikwid.perfmon_getGroupName(groupid))
    return name
end

"""
Return the short information about a performance group.
"""
function get_shortinfo_of_group(groupid::Integer)
    perfmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    sinfo = unsafe_string(LibLikwid.perfmon_getGroupInfoShort(groupid))
    return sinfo
end

"""
Return the (long) description of a performance group.
"""
function get_longinfo_of_group(groupid::Integer)
    perfmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    linfo = unsafe_string(LibLikwid.perfmon_getGroupInfoLong(groupid))
    return linfo
end

"""
Return the name of the event identified by `groupid` and `eventidx`.
"""
function get_name_of_event(groupid::Integer, eventidx::Integer)
    if !perfmon_initialized[] ||
       !_check_groupid(groupid) ||
       !_check_eventidx(groupid, eventidx)
        return nothing
    end
    name = unsafe_string(LibLikwid.perfmon_getEventName(groupid, eventidx))
    return name
end

"""
Return the name of the counter register identified by `groupid` and `eventidx`.
"""
function get_name_of_counter(groupid::Integer, eventidx::Integer)
    if !perfmon_initialized[] ||
       !_check_groupid(groupid) ||
       !_check_eventidx(groupid, eventidx)
        return nothing
    end
    name = unsafe_string(LibLikwid.perfmon_getCounterName(groupid, eventidx))
    return name
end

"""
Return the name of a derived metric identified by `groupid` and `metricidx`.
"""
function get_name_of_metric(groupid::Integer, metricidx::Integer)
    if !perfmon_initialized[] ||
       !_check_groupid(groupid) ||
       !_check_metricidx(groupid, metricidx)
        return nothing
    end
    name = unsafe_string(LibLikwid.perfmon_getMetricName(groupid, metricidx))
    return name
end

"""
Program the counter registers to measure all events in group `groupid`. Returns `true` on success.
"""
function setup_counters(groupid::Integer)
    perfmon_initialized[] || return nothing
    _check_groupid(groupid) || return false
    ret = LibLikwid.perfmon_setupCounters(groupid)
    return ret == 0
end

"""
Start the counter registers. Returns `true` on success.
"""
function start_counters()
    perfmon_initialized[] || return nothing
    ret = LibLikwid.perfmon_startCounters()
    return ret == 0
end

"""
Stop the counter registers. Returns `true` on success.
"""
function stop_counters()
    perfmon_initialized[] || return nothing
    ret = LibLikwid.perfmon_stopCounters()
    return ret == 0
end

"""
Read the counter registers.
To be executed after `start_counters` and before `stop_counters`.
Returns `true` on success.
"""
function read_counters()
    perfmon_initialized[] || return nothing
    ret = LibLikwid.perfmon_readCounters()
    return ret == 0
end

"""
Return the `groupid` of the currently activate group.
"""
function get_id_of_active_group()
    perfmon_initialized[] || return nothing
    return LibLikwid.perfmon_getIdOfActiveGroup()
end

"""
Switch currently active group to `groupid`. Returns `true` on success.
"""
function switch_group(groupid::Integer)
    perfmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    groupid == get_id_of_active_group() && return true
    ret = LibLikwid.perfmon_switchActiveGroup(groupid)
    return ret == 0
end

"""
Return the raw counter register result of the last measurement cycle identified by group `groupid` and the indices for event `eventidx` and thread `threadidx`.
"""
function get_last_result(groupid::Integer, eventidx::Integer, threadidx::Integer)
    perfmon_initialized[] || return nothing
    _check_eventidx(groupid, eventidx) || return nothing
    _check_threadidx(threadidx) || return nothing
    res = LibLikwid.perfmon_getLastResult(groupid, eventidx, threadidx)
    return res
end

"""
Return the raw counter register result of all measurements identified by group `groupid` and the indices for event `eventidx` and thread `threadidx`.
"""
function get_result(groupid::Integer, eventidx::Integer, threadidx::Integer)
    perfmon_initialized[] || return nothing
    _check_eventidx(groupid, eventidx) || return nothing
    _check_threadidx(threadidx) || return nothing
    res = LibLikwid.perfmon_getResult(groupid, eventidx, threadidx)
    return res
end

"""
Return the derived metric result of all measurements identified by group `groupid` and the indices for metric `metricidx` and thread `threadidx`.
"""
function get_metric(groupid::Integer, metricidx::Integer, threadidx::Integer)
    perfmon_initialized[] || return nothing
    _check_metricidx(groupid, metricidx) || return nothing
    _check_threadidx(threadidx) || return nothing
    res = LibLikwid.perfmon_getMetric(groupid, metricidx, threadidx)
    return res
end

"""
Return the derived metric result of the last measurement cycle identified by group `groupid` and the indices for metric `metricidx` and thread `threadidx`.
"""
function get_last_metric(groupid::Integer, metricidx::Integer, threadidx::Integer)
    perfmon_initialized[] || return nothing
    _check_metricidx(groupid, metricidx) || return nothing
    _check_threadidx(threadidx) || return nothing
    res = LibLikwid.perfmon_getLastMetric(groupid, metricidx, threadidx)
    return res
end

"""
Return the measurement time for group identified by `groupid`.
"""
function get_time_of_group(groupid::Integer)
    perfmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    time = LibLikwid.perfmon_getTimeOfGroup(groupid)
    return time
end

"""
List all the metrics of a given group.
"""
function list_metrics(groupid::Integer)
    perfmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    nmetrics = get_number_of_metrics(groupid)
    return get_name_of_metric.(Ref(groupid), 0:nmetrics-1)
end

"""
Get the name and results of all metrics of a given group for a given cpu thread.
"""
function get_metric_results(groupid::Integer, threadid::Integer)
    perfmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    _check_threadidx(threadid) || return nothing
    nmetrics = get_number_of_metrics(groupid)
    d = OrderedDict{String, Float64}()
    for m in 1:nmetrics
        metricid = m-1
        metric = get_name_of_metric(groupid, metricid)
        d[metric] = get_last_metric(groupid, metricid, threadid)
    end
    return d
end

"""
Get the name and results of all events of a given group for a given cpu thread.
"""
function get_event_results(groupid::Integer, threadid::Integer)
    perfmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    _check_threadidx(threadid) || return nothing
    nevents = get_number_of_events(groupid)
    d = OrderedDict{String, Float64}()
    for i in 1:nevents
        eventid = i-1
        event = get_name_of_event(groupid, eventid)
        d[event] = get_last_result(groupid, eventid, threadid)
    end
    return d
end

end # module


text:
function are_type_and_fields_in_output(obj::T) where {T <: Component}
    match = true
    normal = repr(obj)
    io = IOBuffer()
    show(io, "text/plain", obj)
    custom = String(take!(io))
    fields = fieldnames(T)

    # Type must always be present. name should be also, if the type defines it.
    #for text in (normal, custom)
    for text in (custom,)
        if !occursin(string(T), text)
            @error "type name is not in output" string(T) text
            match = false
        end
        if :name in fields
            if !occursin(obj.name, text)
                @error "name is not in output" name text
                match = false
            end
        end
    end

    for (name, type) in zip(fields, fieldtypes(T))
        val = getfield(obj, name)
        if val === nothing || type <: IS.InfrastructureSystemsInternal
            continue
        end

        # Account for the fact that type may be abstract.
        actual_type = typeof(val)
        if actual_type <: IS.InfrastructureSystemsType ||
           actual_type <: IS.TimeSeriesContainer ||
           actual_type <: Vector{<:Service}
            expected = string(actual_type)
        else
            expected = string(val)
        end

        if !occursin(expected, custom)
            @error "field's value is not in custom output" name custom
            match = false
        end
    end

    return match
end

@testset "Test printing of system and components" begin
    sys = PSB.build_system(PSITestSystems, "test_RTS_GMLC_sys")
    @test are_type_and_fields_in_output(iterate(get_components(Bus, sys))[1])
    @test are_type_and_fields_in_output(iterate(get_components(Generator, sys))[1])
    @test are_type_and_fields_in_output(iterate(get_components(ThermalGen, sys))[1])
    @test are_type_and_fields_in_output(iterate(get_components(Branch, sys))[1])
    @test are_type_and_fields_in_output(iterate(get_components(ElectricLoad, sys))[1])

    # Just make sure nothing blows up.
    for component in iterate_components(sys)
        print(devnull, component)
        print(devnull, MIME"text/plain")
        @test !isempty(summary(component))
    end
    for time_series in get_time_series_multiple(sys)
        show(devnull, time_series)
        show(devnull, MIME"text/plain")
        @test !isempty(summary(time_series))
    end

    @test !isempty(summary(sys))

    @test isnothing(
        show(
            IOBuffer(),
            "text/plain",
            PowerSystemTableData(RTS_GMLC_DIR, 100.0, DESCRIPTORS),
        ),
    )

    @test isnothing(show(IOBuffer(), "text/plain", sys))
end

@testset "Test printing of non-PowerSystems struct" begin
    struct MyComponent <: Component
        name::String
        internal::IS.InfrastructureSystemsInternal
    end

    PSY.get_internal(x::MyComponent) = x.internal
    PSY.get_name(x::MyComponent) = string(x.name)

    component = MyComponent("component1", IS.InfrastructureSystemsInternal())
    @test isnothing(show(IOBuffer(), component))
    @test isnothing(show(IOBuffer(), "text/plain", component))
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 463264a9-0f9b-43ee-930f-74730e158372
# we write a function that caculates the product of two matrices A and B
function mult(A,B)

	# we find the sizes of A and B
	aRows, aCols = size(A)
	bRows, bCols = size(B)

	# we check if the matrix sizes are consistent with multiplication
	if aCols != bRows
		error("matrix sizes don't agree")
	end

	# we initialize C to a zero matrix of the appropriate size
	C = zeros(aRows, bCols)

	# we loop over all pairs of rows and columns
	for i = 1:aRows
		for j = 1:bCols
			# and for each pair we store their dot pruduct in C
			for k = 1:aCols
				C[i,j] += A[i,k] * B[k,j];
			end
		end
	end

	# we return C as the answer
	return C
end

# ╔═╡ c089bff8-8b2b-420e-95d4-f40026afb897
md"Simple test"

# ╔═╡ 5d274440-f3b4-4b98-953c-340ca2fb8b8e
A = [1 2; 3 4];

# ╔═╡ a8dc36db-ffc4-4068-af75-890dd3a98819
B = [2 1; 1 3];

# ╔═╡ ab386751-af17-41f9-b032-3735fd3492f4
mult(A, B)

# ╔═╡ 1073d7c1-cd7f-41eb-9531-b7a63c2d6286
md"Check with MATLAB built-in multiplication"

# ╔═╡ eb028ed8-2ff5-4c97-8b98-156f45402c7c
A*B

# ╔═╡ 507e8447-6446-44b3-a230-940254bedb6b
md"Two $100\times 100$ matrices of random numbers"

# ╔═╡ 413585fe-dbe5-4351-88b0-57cb415c31e5
time100 = @elapsed mult(rand(100,100), rand(100,100))

# ╔═╡ f014193b-04d6-430f-83c1-2bee5864ea99
md"Two $200\times 200$ matrices of random numbers"

# ╔═╡ 0faabc1c-a618-4351-9876-9fce4ceab769
time200 = @elapsed mult(rand(200,200), rand(200,200))

# ╔═╡ 31917d5e-28cf-4caf-80a4-50e630094901
md"Two $1000\times 1000$ matrices of random numbers"

# ╔═╡ 7c242e1a-c41c-4e1b-b24b-7343417fbfdd
time1000 = @elapsed mult(rand(1000,1000), rand(1000,1000))

# ╔═╡ 695be13b-3dd7-4943-a339-17872c1d128e
md"""
Now compute the ratio of the times taken for the two different sizes
and compare with the theroetical complexity $O(n^3)$, where $n$ is the size of the matrices. Note the prediction does not match exactly because it is valid only as $n$ tends to infinity, but it is a reasonable guide.
"""

# ╔═╡ 99bdb4c4-4af0-410b-a3f8-3e3e8e13d80f
ratio = time1000/time200

# ╔═╡ 25ae1cca-d886-46ca-9e57-7abd83c7ed86
predictedRatio = (1000/200)^3

# ╔═╡ 23559250-2979-40e1-87b1-f48b7eeee335
md"""
The last part shows that MATLAB is rather efficient at dealing with matrices...
"""

# ╔═╡ 426e472d-c146-41de-bf1f-bbe635a4f8cd
time = @elapsed rand(1000,1000)*rand(1000,1000)

# ╔═╡ Cell order:
# ╠═463264a9-0f9b-43ee-930f-74730e158372
# ╟─c089bff8-8b2b-420e-95d4-f40026afb897
# ╠═5d274440-f3b4-4b98-953c-340ca2fb8b8e
# ╠═a8dc36db-ffc4-4068-af75-890dd3a98819
# ╠═ab386751-af17-41f9-b032-3735fd3492f4
# ╟─1073d7c1-cd7f-41eb-9531-b7a63c2d6286
# ╠═eb028ed8-2ff5-4c97-8b98-156f45402c7c
# ╟─507e8447-6446-44b3-a230-940254bedb6b
# ╠═413585fe-dbe5-4351-88b0-57cb415c31e5
# ╟─f014193b-04d6-430f-83c1-2bee5864ea99
# ╠═0faabc1c-a618-4351-9876-9fce4ceab769
# ╟─31917d5e-28cf-4caf-80a4-50e630094901
# ╠═7c242e1a-c41c-4e1b-b24b-7343417fbfdd
# ╟─695be13b-3dd7-4943-a339-17872c1d128e
# ╠═99bdb4c4-4af0-410b-a3f8-3e3e8e13d80f
# ╠═25ae1cca-d886-46ca-9e57-7abd83c7ed86
# ╟─23559250-2979-40e1-87b1-f48b7eeee335
# ╠═426e472d-c146-41de-bf1f-bbe635a4f8cd

max_stars_repo_path: 1a/Matrix Multiplication/mult_test.jl

max_issues_repo_path: 1a/Matrix Multiplication/mult_test.jl

max_forks_repo_path: 1a/Matrix Multiplication/mult_test.jl


text:
# Filter prototypes, transformations, and transforms

using ..Windows

abstract type FilterType end

#
# Butterworth prototype
#

function Butterworth(::Type{T}, n::Integer) where {T<:Real}
    n > 0 || error("n must be positive")

    poles = zeros(Complex{T}, n)
    for i = 1:div(n, 2)
        w = convert(T, 2i-1)/2n
        pole = complex(-sinpi(w), cospi(w))
        poles[2i-1] = pole
        poles[2i] = conj(pole)
    end
    if isodd(n)
        poles[end] = -1
    end
    ZeroPoleGain{:s}(Complex{T}[], poles, one(T))
end

"""
    Butterworth(n)

``n`` pole Butterworth filter.
"""
Butterworth(n::Integer) = Butterworth(Float64, n)

#
# Chebyshev type I and II prototypes
#

function chebyshev_poles(::Type{T}, n::Integer, ε::Real) where {T<:Real}
    p = zeros(Complex{T}, n)
    μ = asinh(convert(T, 1)/ε)/n
    b = -sinh(μ)
    c = cosh(μ)
    for i = 1:div(n, 2)
        w = convert(T, 2i-1)/2n
        pole = complex(b*sinpi(w), c*cospi(w))
        p[2i-1] = pole
        p[2i] = conj(pole)
    end
    if isodd(n)
        w = convert(T, 2*div(n, 2)+1)/2n
        pole = b*sinpi(w)
        p[end] = pole
    end
    p
end

function Chebyshev1(::Type{T}, n::Integer, ripple::Real) where {T<:Real}
    n > 0 || error("n must be positive")
    ripple >= 0 || error("ripple must be non-negative")

    ε = sqrt(10^(convert(T, ripple)/10)-1)
    p = chebyshev_poles(T, n, ε)
    k = one(T)
    for i = 1:div(n, 2)
        k *= abs2(p[2i])
    end
    if iseven(n)
        k /= sqrt(1+abs2(ε))
    else
        k *= real(-p[end])
    end
    ZeroPoleGain{:s}(Complex{T}[], p, k)
end

"""
    Chebyshev1(n, ripple)

`n` pole Chebyshev type I filter with `ripple` dB ripple in
the passband.
"""
Chebyshev1(n::Integer, ripple::Real) = Chebyshev1(Float64, n, ripple)

function Chebyshev2(::Type{T}, n::Integer, ripple::Real) where {T<:Real}
    n > 0 || error("n must be positive")
    ripple >= 0 || error("ripple must be non-negative")

    ε = 1/sqrt(10^(convert(T, ripple)/10)-1)
    p = chebyshev_poles(T, n, ε)
    for i = 1:length(p)
        p[i] = inv(p[i])
    end

    z = zeros(Complex{T}, n-isodd(n))
    k = one(T)
    for i = 1:div(n, 2)
        w = convert(T, 2i-1)/2n
        ze = complex(zero(T), -inv(cospi(w)))
        z[2i-1] = ze
        z[2i] = conj(ze)
        k *= abs2(p[2i])/abs2(ze)
    end
    isodd(n) && (k *= -real(p[end]))

    ZeroPoleGain{:s}(z, p, k)
end

"""
    Chebyshev2(n, ripple)

`n` pole Chebyshev type II filter with `ripple` dB ripple in
the stopband.
"""
Chebyshev2(n::Integer, ripple::Real) = Chebyshev2(Float64, n, ripple)

#
# Elliptic prototype
#
# See Orfanidis, S. J. (2007). Lecture notes on elliptic filter design.
# Retrieved from http://www.ece.rutgers.edu/~orfanidi/ece521/notes.pdf

# Compute Landen sequence for evaluation of elliptic functions
function landen(k::Real)
    niter = 7
    kn = Vector{typeof(k)}(undef, niter)
    # Eq. (50)
    for i = 1:niter
        kn[i] = k = abs2(k/(1+sqrt(1-abs2(k))))
    end
    kn
end

# cde computes cd(u*K(k), k)
# sne computes sn(u*K(k), k)
# Both accept the Landen sequence as generated by landen above
for (fn, init) in ((:cde, :(cospi(u/2))), (:sne, :(sinpi(u/2))))
    @eval begin
        function $fn(u::Number, landen::Vector{T}) where T<:Real
            winv = inv($init)
            # Eq. (55)
            for i = length(landen):-1:1
                oldwinv = winv
                winv = 1/(1+landen[i])*(winv+landen[i]/winv)
            end
            w = inv(winv)
        end
    end
end

# sne inverse
function asne(w::Number, k::Real)
    oldw = NaN
    while w != oldw
        oldw = w
        kold = k
        # Eq. (50)
        k = abs2(k/(1+sqrt(1-abs2(k))))
        # Eq. (56)
        w = 2*w/((1+k)*(1+sqrt(1-abs2(kold)*w^2)))
    end
    2*asin(w)/π
end

function Elliptic(::Type{T}, n::Integer, rp::Real, rs::Real) where {T<:Real}
    n > 0 || error("n must be positive")
    rp > 0 || error("rp must be positive")
    rp < rs || error("rp must be less than rs")

    # Eq. (2)
    εp = sqrt(10^(convert(T, rp)/10)-1)
    εs = sqrt(10^(convert(T, rs)/10)-1)

    # Eq. (3)
    k1 = εp/εs
    k1 >= 1 && error("filter order is too high for parameters")

    # Eq. (20)
    k1′² = 1 - abs2(k1)
    k1′ = sqrt(k1′²)
    k1′_landen = landen(k1′)

    # Eq. (47)
    k′ = one(T)
    for i = 1:div(n, 2)
        k′ *= sne(convert(T, 2i-1)/n, k1′_landen)
    end
    k′ = k1′²^(convert(T, n)/2)*k′^4

    k = sqrt(1 - abs2(k′))
    k_landen = landen(k)

    # Eq. (65)
    v0 = -im/convert(T, n)*asne(im/εp, k1)

    z = Vector{Complex{T}}(undef, 2*div(n, 2))
    p = Vector{Complex{T}}(undef, n)
    gain = one(T)
    for i = 1:div(n, 2)
        # Eq. (43)
        w = convert(T, 2i-1)/n

        # Eq. (62)
        ze = complex(zero(T), -inv(k*cde(w, k_landen)))
        z[2i-1] = ze
        z[2i] = conj(ze)

        # Eq. (64)
        pole = im*cde(w - im*v0, k_landen)
        p[2i] = pole
        p[2i-1] = conj(pole)

        gain *= abs2(pole)/abs2(ze)
    end

    if isodd(n)
        pole = im*sne(im*v0, k_landen)
        p[end] = pole
        gain *= abs(pole)
    else
        gain *= 10^(-convert(T, rp)/20)
    end

    ZeroPoleGain{:s}(z, p, gain)
end

"""
    Elliptic(n, rp, rs)

`n` pole elliptic (Cauer) filter with `rp` dB ripple in the
passband and `rs` dB attentuation in the stopband.
"""
Elliptic(n::Integer, rp::Real, rs::Real) = Elliptic(Float64, n, rp, rs)

#
# Prototype transformation types
#

# returns frequency in half-cycles per sample ∈ (0, 1)
function normalize_freq(w::Real, fs::Real)
    w <= 0 && error("frequencies must be positive")
    f = 2*w/fs
    f >= 1 && error("frequencies must be less than the Nyquist frequency $(fs/2)")
    f
end

struct Lowpass{T} <: FilterType
    w::T
end

"""
    Lowpass(Wn[; fs])

Low pass filter with cutoff frequency `Wn`. If `fs` is not
specified, `Wn` is interpreted as a normalized frequency in
half-cycles/sample.
"""
Lowpass(w::Real; fs::Real=2) = Lowpass{typeof(w/1)}(normalize_freq(w, fs))

struct Highpass{T} <: FilterType
    w::T
end

"""
    Highpass(Wn[; fs])

High pass filter with cutoff frequency `Wn`. If `fs` is not
specified, `Wn` is interpreted as a normalized frequency in
half-cycles/sample.
"""
Highpass(w::Real; fs::Real=2) = Highpass{typeof(w/1)}(normalize_freq(w, fs))

struct Bandpass{T} <: FilterType
    w1::T
    w2::T
end

"""
    Bandpass(Wn1, Wn2[; fs])

Band pass filter with normalized pass band (`Wn1`, `Wn2`). If
`fs` is not specified, `Wn1` and `Wn2` are interpreted as
normalized frequencies in half-cycles/sample.
"""
function Bandpass(w1::Real, w2::Real; fs::Real=2)
    w1 < w2 || error("w1 must be less than w2")
    Bandpass{Base.promote_typeof(w1/1, w2/1)}(normalize_freq(w1, fs), normalize_freq(w2, fs))
end

struct Bandstop{T} <: FilterType
    w1::T
    w2::T
end

"""
    Bandstop(Wn1, Wn2[; fs])

Band stop filter with normalized stop band (`Wn1`, `Wn2`). If
`fs` is not specified, `Wn1` and `Wn2` are interpreted as
normalized frequencies in half-cycles/sample.
"""
function Bandstop(w1::Real, w2::Real; fs::Real=2)
    w1 < w2 || error("w1 must be less than w2")
    Bandstop{Base.promote_typeof(w1/1, w2/1)}(normalize_freq(w1, fs), normalize_freq(w2, fs))
end

#
# Prototype transformation implementations
#
# The formulas implemented here come from the documentation for the
# corresponding functionality in Octave, available at
# https://staff.ti.bfh.ch/sha1/Octave/index/f/sftrans.html
# The Octave implementation was not consulted in creating this code.

# Create a lowpass filter from a lowpass filter prototype
transform_prototype(ftype::Lowpass, proto::ZeroPoleGain{:s, Z, P, K}) where {Z, P, K} =
    ZeroPoleGain{:s, Z, P, K}(ftype.w * proto.z, ftype.w * proto.p,
              proto.k * ftype.w^(length(proto.p)-length(proto.z)))

# Create a highpass filter from a lowpass filter prototype
function transform_prototype(ftype::Highpass, proto::ZeroPoleGain{:s})
    z = proto.z
    p = proto.p
    k = proto.k
    nz = length(z)
    np = length(p)
    TR = Base.promote_eltype(z, p)
    newz = zeros(TR, max(nz, np))
    newp = zeros(TR, max(nz, np))

    num = one(eltype(z))
    for i = 1:nz
        num *= -z[i]
        newz[i] = ftype.w / z[i]
    end

    den = one(eltype(p))
    for i = 1:np
        den *= -p[i]
        newp[i] = ftype.w / p[i]
    end

    abs(real(num) - 1) < np*eps(real(num)) && (num = 1)
    abs(real(den) - 1) < np*eps(real(den)) && (den = 1)
    ZeroPoleGain{:s}(newz, newp, oftype(k, k * real(num)/real(den)))
end

# Create a bandpass filter from a lowpass filter prototype
function transform_prototype(ftype::Bandpass, proto::ZeroPoleGain{:s})
    z = proto.z
    p = proto.p
    k = proto.k
    nz = length(z)
    np = length(p)
    ncommon = min(nz, np)
    TR = Base.promote_eltype(z, p)
    newz = zeros(TR, 2*nz+np-ncommon)
    newp = zeros(TR, 2*np+nz-ncommon)
    for (oldc, newc) in ((p, newp), (z, newz))
        for i = 1:length(oldc)
            b = oldc[i] * ((ftype.w2 - ftype.w1)/2)
            pm = sqrt(b^2 - ftype.w2 * ftype.w1)
            newc[2i-1] = b + pm
            newc[2i] = b - pm
        end
    end
    ZeroPoleGain{:s}(newz, newp, oftype(k, k * (ftype.w2 - ftype.w1) ^ (np - nz)))
end

# Create a bandstop filter from a lowpass filter prototype
function transform_prototype(ftype::Bandstop, proto::ZeroPoleGain{:s})
    z = proto.z
    p = proto.p
    k = proto.k
    nz = length(z)
    np = length(p)
    npairs = nz+np-min(nz, np)
    TR = Base.promote_eltype(z, p)
    newz = Vector{TR}(undef, 2*npairs)
    newp = Vector{TR}(undef, 2*npairs)

    num = one(eltype(z))
    for i = 1:nz
        num *= -z[i]
        b = (ftype.w2 - ftype.w1)/2/z[i]
        pm = sqrt(b^2 - ftype.w2 * ftype.w1)
        newz[2i-1] = b - pm
        newz[2i] = b + pm
    end

    den = one(eltype(p))
    for i = 1:np
        den *= -p[i]
        b = (ftype.w2 - ftype.w1)/2/p[i]
        pm = sqrt(b^2 - ftype.w2 * ftype.w1)
        newp[2i-1] = b - pm
        newp[2i] = b + pm
    end

    # Any emaining poles/zeros are real and not cancelled
    npm = sqrt(-complex(ftype.w2 * ftype.w1))
    for (n, newc) in ((np, newp), (nz, newz))
        for i = n+1:npairs
            newc[2i-1] = -npm
            newc[2i] = npm
        end
    end

    abs(real(num) - 1) < np*eps(real(num)) && (num = 1)
    abs(real(den) - 1) < np*eps(real(den)) && (den = 1)
    ZeroPoleGain{:s}(newz, newp, oftype(k, k * real(num)/real(den)))
end

transform_prototype(ftype, proto::FilterCoefficients{:s}) =
    transform_prototype(ftype, convert(ZeroPoleGain, proto))

"""
    analogfilter(responsetype, designmethod)

Construct an analog filter. See below for possible response and
filter types.
"""
analogfilter(ftype::FilterType, proto::FilterCoefficients) =
    transform_prototype(ftype, proto)

# Bilinear transform
bilinear(f::FilterCoefficients{:s}, fs::Real) = bilinear(convert(ZeroPoleGain, f), fs)
function bilinear(f::ZeroPoleGain{:s,Z,P,K}, fs::Real) where {Z,P,K}
    ztype = typeof(0 + zero(Z)/fs)
    z = fill(convert(ztype, -1), max(length(f.p), length(f.z)))

    ptype = typeof(0 + zero(P)/fs)
    p = Vector{typeof(zero(P)/fs)}(undef, length(f.p))

    num = one(one(fs) - one(Z))
    for i = 1:length(f.z)
        z[i] = (2 + f.z[i] / fs)/(2 - f.z[i] / fs)
        num *= (2 * fs - f.z[i])
    end

    den = one(one(fs) - one(P))
    for i = 1:length(f.p)
        p[i] = (2 + f.p[i] / fs)/(2 - f.p[i]/fs)
        den *= (2 * fs - f.p[i])
    end

    ZeroPoleGain{:z}(z, p, f.k * real(num)/real(den))
end

# Pre-warp filter frequencies for digital filtering
prewarp(ftype::Union{Lowpass, Highpass}) = (typeof(ftype))(prewarp(ftype.w))
prewarp(ftype::Union{Bandpass, Bandstop}) = (typeof(ftype))(prewarp(ftype.w1), prewarp(ftype.w2))
# freq in half-samples per cycle
prewarp(f::Real) = 4*tan(pi*f/2)

# Digital filter design
"""
    digitalfilter(responsetype, designmethod)

Construct a digital filter. See below for possible response and
filter types.
"""
digitalfilter(ftype::FilterType, proto::FilterCoefficients) =
    bilinear(transform_prototype(prewarp(ftype), proto), 2)

#
# Special filter types
#

"""
    iirnotch(Wn, bandwidth[; fs])

Second-order digital IIR notch filter [^Orfandis] at frequency `Wn` with
bandwidth `bandwidth`. If `fs` is not specified, `Wn` is
interpreted as a normalized frequency in half-cycles/sample.

[^Orfandis]: Orfanidis, S. J. (1996). Introduction to signal processing. Englewood Cliffs, N.J: Prentice Hall, p. 370.
"""
function iirnotch(w::Real, bandwidth::Real; fs=2)
    w = normalize_freq(w, fs)
    bandwidth = normalize_freq(bandwidth, fs)

    # Eq. 8.2.23
    b = 1/(1+tan(pi*bandwidth/2))
    # Eq. 8.2.22
    cosw0 = cospi(w)
    Biquad(b, -2b*cosw0, b, -2b*cosw0, 2b-1)
end

#
# FIR filter design
#

# Get length and alpha for Kaiser window filter with specified
# transition band width and stopband attenuation in dB
function kaiserord(transitionwidth::Real, attenuation::Real=60)
    n = ceil(Int, (attenuation - 7.95)/(π*2.285*transitionwidth))+1

    if attenuation > 50
        β = 0.1102*(attenuation - 8.7)
    elseif attenuation >= 21
        β = 0.5842*(attenuation - 21)^0.4 + 0.07886*(attenuation - 21)
    else
        β = 0.0
    end

    return n, β/π
end

struct FIRWindow{T}
    window::Vector{T}
    scale::Bool
end

"""
    FIRWindow(window; scale=true)

FIR filter design using window `window`, a vector whose length
matches the number of taps in the resulting filter.

If `scale` is `true` (default), the designed FIR filter is
scaled so that the following holds:

- For [`Lowpass`](@ref) and [`Bandstop`](@ref) filters, the frequency
  response is unity at 0 (DC).
- For [`Highpass`](@ref) filters, the frequency response is unity
  at the Nyquist frequency.
- For [`Bandpass`](@ref) filters, the frequency response is unity
  in the center of the passband.
"""
FIRWindow(window::Vector; scale::Bool=true) = FIRWindow(window, scale)

# FIRWindow(n::Integer, window::Function, args...) = FIRWindow(window(n, args...))
"""
    FIRWindow(; transitionwidth, attenuation=60, scale=true)

Kaiser window FIR filter design. The required number of taps is
calculated based on `transitionwidth` (in half-cycles/sample)
and stopband `attenuation` (in dB). `attenuation` defaults to
60 dB.
"""
FIRWindow(; transitionwidth::Real=throw(ArgumentError("must specify transitionwidth")),
          attenuation::Real=60, scale::Bool=true) =
    FIRWindow(kaiser(kaiserord(transitionwidth, attenuation)...), scale)

# Compute coefficients for FIR prototype with specified order
function firprototype(n::Integer, ftype::Lowpass)
    w = ftype.w

    [w*sinc(w*(k-(n-1)/2)) for k = 0:(n-1)]
end

function firprototype(n::Integer, ftype::Bandpass)
    w1 = ftype.w1
    w2 = ftype.w2

    [w2*sinc(w2*(k-(n-1)/2)) - w1*sinc(w1*(k-(n-1)/2)) for k = 0:(n-1)]
end

function firprototype(n::Integer, ftype::Highpass)
    w = ftype.w
    isodd(n) || throw(ArgumentError("FIRWindow highpass filters must have an odd number of coefficients"))

    out = [-w*sinc(w*(k-(n-1)/2)) for k = 0:(n-1)]
    out[div(n, 2)+1] += 1
    out
end

function firprototype(n::Integer, ftype::Bandstop)
    w1 = ftype.w1
    w2 = ftype.w2
    isodd(n) || throw(ArgumentError("FIRWindow bandstop filters must have an odd number of coefficients"))

    out = [w1*sinc(w1*(k-(n-1)/2)) - w2*sinc(w2*(k-(n-1)/2)) for k = 0:(n-1)]
    out[div(n, 2)+1] += 1
    out
end

scalefactor(coefs::Vector, ::Union{Lowpass, Bandstop}) = sum(coefs)
function scalefactor(coefs::Vector, ::Highpass)
    c = zero(coefs[1])
    for k = 1:length(coefs)
        c += ifelse(isodd(k), coefs[k], -coefs[k])
    end
    c
end
function scalefactor(coefs::Vector, ftype::Bandpass)
    n = length(coefs)
    freq = middle(ftype.w1, ftype.w2)
    c = zero(coefs[1])
    for k = 0:n-1
        c += coefs[k+1]*cospi(freq*(k-(n-1)/2))
    end
    c
end

function digitalfilter(ftype::FilterType, proto::FIRWindow)
    coefs = firprototype(length(proto.window), ftype)
    @assert length(proto.window) == length(coefs)
    out = coefs .* proto.window
    proto.scale ? rmul!(out, 1/scalefactor(out, ftype)) : out
end

# Compute FIR coefficients necessary for arbitrary rate resampling
function resample_filter(rate::AbstractFloat, Nϕ = 32, rel_bw = 1.0, attenuation = 60)
    f_nyq       = rate >= 1.0 ? 1.0/Nϕ : rate/Nϕ
    cutoff      = f_nyq * rel_bw
    trans_width = cutoff * 0.2

    # Determine resampling filter order
    hLen, α = kaiserord(trans_width, attenuation)

    # Round the number of taps up to a multiple of Nϕ.
    # Otherwise the missing taps will be filled with 0.
    hLen = Nϕ * ceil(Int, hLen/Nϕ)

    # Ensure that the filter is an odd length
    if (iseven(hLen))
        hLen += 1
    end

    # Design filter
    h = digitalfilter(Lowpass(cutoff), FIRWindow(kaiser(hLen, α)))
    rmul!(h, Nϕ)
end

# Compute FIR coefficients necessary for rational rate resampling
function resample_filter(rate::Union{Integer,Rational}, rel_bw = 1.0, attenuation = 60)
    Nϕ          = numerator(rate)
    decimation  = denominator(rate)
    f_nyq       = min(1/Nϕ, 1/decimation)
    cutoff      = f_nyq * rel_bw
    trans_width = cutoff * 0.2

    # Determine resampling filter order
    hLen, α = kaiserord(trans_width, attenuation)

    # Round the number of taps up to a multiple of Nϕ (same as interpolation factor).
    # Otherwise the missing taps will be filled with 0.
    hLen = Nϕ * ceil(Int, hLen/Nϕ)

    # Ensure that the filter is an odd length
    if (iseven(hLen))
        hLen += 1
    end

    # Design filter
    h = digitalfilter(Lowpass(cutoff), FIRWindow(kaiser(hLen, α)))
    rmul!(h, Nϕ)
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module NonNegLeastSquares

using Compat
using Distributed
using LinearAlgebra
import SparseArrays

export nonneg_lsq

## Algorithms
include("nnls.jl")
using .NNLS: nnls
include("fnnls.jl")
include("pivot.jl")
include("pivot_comb.jl")
include("pivot_cache.jl")
include("admm.jl")

## Common interface to algorithms
include("interface.jl")

## Helper functions
include("cssls.jl") # combinatorial subspace least squares (CSSLS)

end


text:
###
### Counting
###
### Counting operations on biological sequence types.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

###
### Naive counting
###

function count_naive(pred, seq::BioSequence)
    n = 0
    @inbounds for i in eachindex(seq)
        n += pred(seq[i])::Bool
    end
    return n
end

function count_naive(pred, seqa::BioSequence, seqb::BioSequence)
    n = 0
    @inbounds for i in 1:min(length(seqa), length(seqb))
        n += pred(seqa[i], seqb[i])::Bool
    end
    return n
end

"""
Count how many positions in a sequence satisfy a condition (i.e. f(seq[i]) -> true).

The first argument should be a function which accepts an element of the sequence
as its first parameter, additional arguments may be passed with `args...`.
"""
Base.count(pred, seq::BioSequence) = count_naive(pred, seq)

Base.count(pred, seqa::BioSequence, seqb::BioSequence) = count_naive(pred, seqa, seqb)

BioSymbols.isambiguous(x::T, y::T) where {T<:NucleicAcid} = isambiguous(x) | isambiguous(y)
BioSymbols.isgap(x::T, y::T) where {T<:NucleicAcid} = isgap(x) | isgap(y)
BioSymbols.iscertain(x::T, y::T) where {T<:NucleicAcid} = iscertain(x) & iscertain(y)

Base.count(::typeof(isambiguous), seqa::S, seqb::S) where {S<:BioSequence{<:NucleicAcidAlphabet{2}}} = 0
Base.count(::typeof(isgap), seqa::S, seqb::S) where {S<:BioSequence{<:NucleicAcidAlphabet{2}}} = 0
Base.count(::typeof(iscertain), seqa::S, seqb::S) where {S<:BioSequence{<:NucleicAcidAlphabet{2}}} = min(length(seqa), length(seqb))

###
### Aliases for various uses of `count`.
###

"""
    gc_content(seq::BioSequence)

Calculate GC content of `seq`.
"""
gc_content(seq::NucleotideSeq) = isempty(seq) ? 0.0 : count(isGC, seq) / length(seq)

n_ambiguous(seq) = count(isambiguous, seq)
n_ambiguous(seqa::BioSequence, seqb::BioSequence) = count(isambiguous, seqa, seqb)

n_certain(seq) = count(iscertain, seq)
n_certain(seqa::BioSequence, seqb::BioSequence) = count(iscertain, seqa, seqb)

n_gaps(seq::BioSequence) = count(isgap, seq)
n_gaps(seqa::BioSequence, seqb::BioSequence) = count(isgap, seqa, seqb)

mismatches(seqa::BioSequence, seqb::BioSequence) = count(!=, seqa, seqb)
matches(seqa::BioSequence, seqb::BioSequence) = count(==, seqa, seqb)


text:
#!/usr/bin/env julia

# Given some model parameters, synthesize the images. Need to run DJInitialize.jl first.

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskJockey.constants
using DiskJockey.model

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]

pars = convert_dict(config["parameters"], config["model"])

npix = config["npix"] # number of pixels

grid = Grid(config["grid"])
(sizeau_desired, sizeau_command) = size_au(config["size_arcsec"], pars.dpc, grid) # [AU]

# lambda should have already been written by DJInitialize.jl
tic()
run(`radmc3d image incl $(pars.incl) posang $(pars.PA) npix $npix loadlambda sizeau $sizeau_command`)
println("Synthesis runtime")
toc()

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
"""
$(TYPEDEF)

Epigraph of real symmetric or complex Hermitian matrix spectral norm (i.e.
maximum absolute value of eigenvalues) cone of dimension `dim` in svec format.

    $(FUNCTIONNAME){T, R}(dim::Int, use_dual::Bool = false)
"""
mutable struct EpiNormSpectralTri{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    d::Int
    is_complex::Bool
    rt2::T

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    hess_aux_updated::Bool
    inv_hess_aux_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    s::Vector{T}
    V::Matrix{R}
    mu::Vector{T}
    zeta::Vector{T}
    Vrzi::Matrix{R}
    Vmrzi::Matrix{R}
    cu::T
    Zu::T
    sumzi::Vector{T}
    umzzzi::Matrix{T}
    zumziz::Matrix{T}

    w1::Matrix{R}
    w2::Matrix{R}
    w3::Matrix{R}
    w4::Matrix{R}
    s1::Vector{T}
    s2::Vector{T}

    function EpiNormSpectralTri{T, R}(
        dim::Int;
        use_dual::Bool = false,
        ) where {T <: Real, R <: RealOrComplex{T}}
        @assert dim >= 2
        cone = new{T, R}()
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        cone.rt2 = sqrt(T(2))
        cone.is_complex = (R <: Complex)
        cone.d = svec_side(R, dim - 1)
        return cone
    end
end

reset_data(cone::EpiNormSpectralTri) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_aux_updated =
    cone.inv_hess_aux_updated = cone.hess_fact_updated = false)

function setup_extra_data!(
    cone::EpiNormSpectralTri{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    d = cone.d
    cone.V = zeros(R, d, d)
    cone.zeta = zeros(T, d)
    cone.mu = zeros(T, d)
    cone.Vrzi = zeros(R, d, d)
    cone.Vmrzi = zeros(R, d, d)
    cone.sumzi = zeros(T, d)
    cone.zumziz = zeros(T, d, d)
    cone.w1 = zeros(R, d, d)
    cone.w2 = zeros(R, d, d)
    cone.w3 = zeros(R, d, d)
    cone.w4 = zeros(R, d, d)
    cone.s1 = zeros(T, d)
    cone.s2 = zeros(T, d)
    return cone
end

get_nu(cone::EpiNormSpectralTri) = 1 + cone.d

function set_initial_point!(
    arr::AbstractVector{T},
    cone::EpiNormSpectralTri{T},
    ) where {T <: Real}
    arr .= 0
    arr[1] = sqrt(T(get_nu(cone)))
    return arr
end

function update_feas(cone::EpiNormSpectralTri{T}) where {T <: Real}
    @assert !cone.feas_updated
    u = cone.point[1]
    cone.feas_updated = true
    cone.is_feas = false
    (u > eps(T)) || return false
    @views W = svec_to_smat!(cone.V, cone.point[2:end], cone.rt2)
    copytri!(W, 'U', true)

    # fast bounds (note opinf = op1):
    # spec <= frob <= spec * rtd
    # op1 / rtd <= spec <= op1
    frob = norm(W, 2)
    op1 = opnorm(W, 1)
    rtd = sqrt(T(cone.d))

    # lower bounds
    lb = max(frob, op1) / rtd
    (u - lb > eps(T)) || return false

    # upper bounds
    ub = min(frob, op1)
    if u - ub < eps(T)
        # use fast Cholesky-based feasibility check, rescale W*W' by inv(u)
        Z = mul!(cone.w1, W, W', -inv(u), false)
        @inbounds for i in 1:cone.d
            Z[i, i] += u
        end
        isposdef(cholesky!(Hermitian(Z, :U), check = false)) || return false
    end

    # compute eigendecomposition and final feasibility check
    cone.s = update_eigen!(cone.V)
    cone.is_feas = (u - maximum(abs, cone.s) > eps(T))

    return cone.is_feas
end

function is_dual_feas(cone::EpiNormSpectralTri{T}) where {T <: Real}
    u = cone.dual_point[1]
    (u > eps(T)) || return false
    @views W = svec_to_smat!(cone.w1, cone.dual_point[2:end], cone.rt2)
    copytri!(W, 'U', true)

    # fast bounds: frob <= nuc <= frob * rtd
    frob = norm(W, 2)
    (u - sqrt(T(cone.d)) * frob > eps(T)) && return true
    (u - frob > eps(T)) || return false

    # final feasibility check
    return (u - sum(abs, eigvals!(Hermitian(W, :U))) > eps(T))
end

function update_grad(cone::EpiNormSpectralTri{T}) where T
    @assert cone.is_feas
    u = cone.point[1]
    V = cone.V
    s = cone.s
    mu = cone.mu
    zeta = cone.zeta
    s1 = cone.s1
    w1 = cone.w1
    w2 = cone.w2
    g = cone.grad

    @. mu = s / u
    @. zeta = T(0.5) * (u - mu * s)
    cone.cu = (cone.d - 1) / u

    g[1] = cone.cu - sum(inv, zeta)

    @. s1 = mu / zeta
    mul!(w1, V, Diagonal(s1))
    mul!(w2, w1, V')
    @views smat_to_svec!(g[2:end], w2, cone.rt2)

    cone.grad_updated = true
    return cone.grad
end

function update_hess_aux(cone::EpiNormSpectralTri{T}) where T
    @assert !cone.hess_aux_updated
    @assert cone.grad_updated
    s1 = cone.s1
    s2 = cone.s2

    @. s1 = sqrt(cone.zeta)
    @. s2 = cone.mu / s1
    @. cone.Vmrzi = s2 * cone.V'
    @. s2 = inv(s1)
    @. cone.Vrzi = s2 * cone.V'

    cone.hess_aux_updated = true
end

function update_hess(cone::EpiNormSpectralTri{T}) where T
    cone.hess_aux_updated || update_hess_aux(cone)
    isdefined(cone, :hess) || alloc_hess!(cone)
    d = cone.d
    isdefined(cone, :umzzzi) || (cone.umzzzi = zeros(T, d, d))
    u = cone.point[1]
    mu = cone.mu
    zeta = cone.zeta
    umzzzi = cone.umzzzi
    s1 = cone.s1
    w1 = cone.w1
    w2 = cone.w2
    H = cone.hess.data
    ui = inv(u)

    # u, w
    @. s1 = -inv(zeta)
    mul!(w2, Diagonal(s1), cone.Vmrzi)
    mul!(w1, cone.Vrzi', w2)
    @views smat_to_svec!(H[1, 2:end], w1, cone.rt2)

    # w, w
    # TODO or write faster symmetric spectral kron
    @inbounds for j in 1:cone.d
        z_j = zeta[j]
        zi_j = inv(z_j)
        mzi_j = mu[j] / z_j
        for i in 1:(j - 1)
            umzzzi[i, j] = umzzzi[j, i] = T(0.5) *
                (zi_j + mzi_j * mu[i]) / zeta[i]
        end
        umzzzi[j, j] = (zi_j - ui) / z_j
    end
    @views Hww = H[2:end, 2:end]
    eig_dot_kron!(Hww, umzzzi, cone.V, w1, w2, cone.w3, cone.rt2)

    # u, u
    H[1, 1] = tr(umzzzi) - cone.cu / u

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::EpiNormSpectralTri{T},
    ) where T
    cone.hess_aux_updated || update_hess_aux(cone)
    d = cone.d
    u = cone.point[1]
    zeta = cone.zeta
    Vrzi = cone.Vrzi
    Vmrzi = cone.Vmrzi
    sim = r = w1 = cone.w1
    w2 = cone.w2
    S1 = cone.w3
    @views S1diag = S1[diagind(S1)]

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views svec_to_smat!(r, arr[2:end, j], cone.rt2)

        pui = p / u
        mul!(w2, Vrzi, Hermitian(r, :U))
        mul!(sim, w2, Vmrzi')
        @. S1 = T(0.5) * (sim + sim')
        @. S1diag -= p / zeta

        prod[1, j] = -sum((pui + real(S1[i, i])) / zeta[i] for i in 1:d) -
            cone.cu * pui

        mul!(w2, Hermitian(S1, :U), Vmrzi, true, inv(u))
        mul!(w1, Vrzi', w2)
        @views smat_to_svec!(prod[2:end, j], w1, cone.rt2)
    end

    return prod
end

function update_inv_hess_aux(cone::EpiNormSpectralTri{T}) where T
    @assert !cone.inv_hess_aux_updated
    @assert cone.grad_updated
    u = cone.point[1]
    s = cone.s
    zeta = cone.zeta
    zumziz = cone.zumziz

    cone.Zu = -cone.cu + sum(inv, u - z_i for z_i in zeta)

    @. cone.sumzi = s / (u - zeta)

    @inbounds for j in 1:cone.d
        mu_j = cone.mu[j]
        z_j = zeta[j]
        for i in 1:(j - 1)
            zumziz[i, j] = zumziz[j, i] = 2 * zeta[i] / (u + mu_j * s[i]) * z_j
        end
        zumziz[j, j] = z_j / (u - z_j) * z_j
    end

    cone.inv_hess_aux_updated = true
end

function update_inv_hess(cone::EpiNormSpectralTri)
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    u = cone.point[1]
    V = cone.V
    s1 = cone.s1
    w1 = cone.w1
    w2 = cone.w2
    Hi = cone.inv_hess.data

    # u, u
    hiuu = Hi[1, 1] = u / cone.Zu

    # w, w
    @views Hiww = Hi[2:end, 2:end]
    eig_dot_kron!(Hiww, cone.zumziz, V, w1, w2, cone.w3, cone.rt2)

    # u, w and w, w
    rthiuu = sqrt(hiuu)
    @. s1 = rthiuu * cone.sumzi
    mul!(w2, V, Diagonal(s1))
    mul!(w1, w2, V')
    @views Hiuwvec = Hi[1, 2:end]
    smat_to_svec!(Hiuwvec, w1, cone.rt2)
    mul!(Hiww, Hiuwvec, Hiuwvec', true, u)
    Hiuwvec .*= rthiuu

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormSpectralTri,
    )
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    d = cone.d
    u = cone.point[1]
    V = cone.V
    sumzi = cone.sumzi
    sim = r = w1 = cone.w1
    w2 = cone.w2
    @views simdiag = sim[diagind(sim)]

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views svec_to_smat!(r, arr[2:end, j], cone.rt2)

        mul!(w2, Hermitian(r, :U), V)
        mul!(sim, V', w2)

        c1 = (p + sum(sumzi[i] * real(sim[i, i]) for i in 1:d)) / cone.Zu
        prod[1, j] = u * c1

        sim .*= cone.zumziz
        @. simdiag += c1 * sumzi
        mul!(w2, V, Hermitian(sim, :U))
        mul!(w1, w2, V', u, false)
        @views smat_to_svec!(prod[2:end, j], w1, cone.rt2)
    end

    return prod
end

function dder3(cone::EpiNormSpectralTri{T}, dir::AbstractVector{T}) where T
    cone.hess_aux_updated || update_hess_aux(cone)
    u = cone.point[1]
    zeta = cone.zeta
    Vrzi = cone.Vrzi
    Vmrzi = cone.Vmrzi
    r = w1 = cone.w1
    simU = w2 = cone.w2
    sim = S2 = cone.w3
    S1 = cone.w4
    @views S1diag = S1[diagind(S1)]
    @views S2diag = S2[diagind(S2)]
    dder3 = cone.dder3

    p = dir[1]
    @views svec_to_smat!(r, dir[2:end], cone.rt2)

    pui = p / u
    mul!(simU, Vrzi, Hermitian(r, :U))
    mul!(sim, simU, Vmrzi')
    @. S1 = T(-0.5) * (sim + sim')
    @. S1diag += p / zeta

    mul!(S2, simU, simU', T(-0.5) / u, false)
    @. S2diag += T(0.5) * p / zeta * pui
    mul!(S2, Hermitian(S1, :U), S1, -1, true)

    @inbounds dder3[1] = -sum((real(S1[i, i]) * pui + real(S2[i, i])) / zeta[i]
        for i in 1:cone.d) - cone.cu * abs2(pui)

    mul!(w1, Hermitian(S2, :U), Vmrzi)
    mul!(w1, Hermitian(S1, :U), simU, inv(u), true)
    mul!(w2, Vrzi', w1)
    @views smat_to_svec!(dder3[2:end], w2, cone.rt2)

    return dder3
end


text:
using BinaryBuilder

name = "Openspecfun"
version = v"0.5.3"

# Collection of sources required to build openspecfun
sources = [
    "https://github.com/JuliaMath/openspecfun/archive/v0.5.3.tar.gz" =>
    "1505c7a45f9f39ffe18be36f7a985cb427873948281dbcd376a11c2cd15e41e7",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd openspecfun-*/

# It needs to be told it's on Windows
if [[ ${target} == *mingw* ]]; then
    OS=WINNT
    libdir=$prefix/bin
elif [[ ${target} == *darwin* ]]; then
    OS=Darwin
fi

# Build it
make OS=${OS} CC="$CC" CXX="$CXX" FC="$FC" -j${nproc}

# Install it
make install OS=${OS} prefix=$prefix
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line. Since openspecfun uses
# Fortran for AMOS, we need the combinatorial explosion of platforms
# and GCC versions.
platforms = expand_gcc_versions(supported_platforms())

# The products that we will ensure are always built
products = prefix -> [
    LibraryProduct(prefix, "libopenspecfun", :libopenspecfun)
]

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

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


text:
using GLVisualize, FileIO, Colors
using GLBooks, Reactive

GLBooks.init()

surf(i, N) = Float32[sin(x*i) * cos(y*i) for x = linspace(0, 2pi, N), y = linspace(0, 2pi, N)]

@block "Overview" begin
    loadasset("doge.png")
end

@block "Surface" begin
    s1 = slider(linspace(0.1f0, 2pi, 100), "surface function")
    drawpage(map(x-> surf(x, 128), s1), :surface)
end

@block "Surface with Image" begin
    video = loadasset("kittens-look.gif")
    play_s = playbutton("play")
    idx = foldp(1, fpswhen(play_s, 30.0)) do v0, t
        mod1(v0 + 1, size(video, 3))
    end
    imstream = map(i-> video[:, :, i], idx)
    s1 = slider(linspace(0.5f0, 2pi, 100), "surface function")
    drawpage(map(x-> surf(x, 128), s1), :surface, color = imstream)
end


text:
qp_solver = optimizer_with_attributes(
    Ipopt.Optimizer,
    "print_level" => 0,
    "warm_start_init_point" => "yes",
)
optimizer_solver = optimizer_with_attributes(
    SQP.Optimizer,
    "external_optimizer" => qp_solver,
    "algorithm" => "SQP-TR",
    "OutputFlag" => 0,
)

model = Model(optimizer_solver)

@variable(model, X);
@variable(model, Y);
@objective(model, Min, X^2 + X);
@NLconstraint(model, X^2 - X == 2);
@NLconstraint(model, X * Y == 1);
@NLconstraint(model, X * Y >= 0);
@constraint(model, X >= -2);

JuMP.optimize!(model);

xsol = JuMP.value.(X)
ysol = JuMP.value.(Y)
status = termination_status(model)


text:
_KGS_COLUMNS = "ABCDEFGHJKLMNOPQRST"
_SGF_COLUMNS = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"

# Converts from a board coordinate to a flattened coordinate
to_flat(coord, pos::GomokuPosition) = to_flat(coord, pos.env)
to_flat(coord, env::GomokuEnv) = coord == nothing ? env.N * env.N + 1 :
                                  env.N * (coord[2] - 1) + coord[1]

# Converts from a flattened coordinate to a board coordinate
from_flat(f, pos::GomokuPosition) = from_flat(f, pos.env)
from_flat(f, env::GomokuEnv) = f == env.N * env.N + 1 ? nothing :
                        tuple([i + 1 for i in reverse(divrem(f - 1, env.N))]...)
#=
function from_sgf(sgfc)
  # Interprets coords. aa is top left corner; sa is top right corner
  if sgfc == nothing || sgfc == ""
    return nothing
  end
  return searchindex(_SGF_COLUMNS, sgfc[2]), searchindex(_SGF_COLUMNS, sgfc[1])
end

# Converts from a board coordinate to an SGF coordinate
to_sgf(coord) =  coord ==  nothing ? "" : _SGF_COLUMNS[coord[2]] * _SGF_COLUMNS[coord[1]]

function from_kgs(kgsc, env)
  # Interprets coords. A1 is bottom left; A9 is top left.
  if kgsc == "pass"
    return nothing
  end
  kgsc = uppercase(kgsc)
  col = searchindex(_KGS_COLUMNS, kgsc[1])
  row_from_bottom = parse(Int, kgsc[2:end])
  return env.N - row_from_bottom + 1, col
end

# Converts from a board coordinate to a KGS coordinate.
to_kgs(coord, env) = coord == nothing ? "pass" : "$(_KGS_COLUMNS[coord[2]])$(env.N - coord[1] + 1)"=#


text:
using Random, Omega, Test

function testrng()
  x = normal(0, 1)
  rng = Random.MersenneTwister(345679)
  s1 = rand(rng, x, x > 0, alg = RejectionSample)
  rng = Random.MersenneTwister(345679)
  s2 = rand(rng, x, x > 0, alg = RejectionSample)
  @test s1 == s2
end

testrng()


text:
using CompScienceMeshes, BEAST

Γ = meshcuboid(1.0, 1.0, 1.0, 0.15)
X = raviartthomas(Γ)

κ,  η  = 1.0, 1.0
κ′, η′ = 1.4κ, η/1.4

γ, γ′ = im*κ, im*κ′
T,  K  = MWSingleLayer3D(γ),  MWDoubleLayer3D(γ)
T′, K′ = MWSingleLayer3D(γ′), MWDoubleLayer3D(γ′)

o, x, y, z = euclidianbasis(3)
d = normalize(x+y+z); p = normalize(d × z)
E = planewavemw3d(direction=d, polarization=p, wavenumber=κ)
H = -1/(im*κ*η)*curl(E)

e, h = (n × E) × n, (n × H) × n

@hilbertspace j m
@hilbertspace k l

α, α′ = 1/η, 1/η′
pmchwt = @discretise(
    (η*T+η′*T′)[k,j] +      (K+K′)[k,m] -
         (K+K′)[l,j] + (α*T+α′*T′)[l,m] == h[k] + e[l],
    j∈X, m∈X, k∈X, l∈X)

u = solve(pmchwt)

## Post-processing
u1 = u[1:numfunctions(X)]
fcr, geo = facecurrents(u1,X)

include(Pkg.dir("CompScienceMeshes","examples","plotlyjs_patches.jl"))
p = patch(geo, real.(norm.(fcr)))
PlotlyJS.plot(p)


text:
using CartesianGrids
using Statistics

@testset "Histories" begin

  a = 1.0
  b = 0.0
  c = 5.0

  h = History(collect(1:5))
  @test mean(h) == 3

  h = History(Nodes(Dual,(5,5)))
  d = Nodes(Dual,(5,5))
  fill!(d,a)
  push!(h,deepcopy(d))
  fill!(d,b)
  push!(h,deepcopy(d))
  fill!(d,c)
  push!(h,deepcopy(d))

  @test length(h.r) == length(h)

  dh = diff(h)

  @test typeof(dh) <: History

  @test dh[2][1,1] == c-b

  # test periodic history
  hp = History(h.vec,htype=PeriodicHistory)
  @test hp[4] == hp[1]

  # test the differencing of periodic history
  dhp = diff(hp)
  @test typeof(dhp) <: History
  @test dhp[1][1,1] == b-a
  @test dhp[5] == dhp[2]

  # test staggered indexing
  hp2 = History(hp[1:2:5],htype=PeriodicHistory)
  @test hp2[2][1,1] == c

  # test circular shifting
  hpshift = circshift(hp,1)
  @test hpshift[1][1,1] == c

  # testing arithmetic
  h3 = 3*h
  @test h3[3][1,1] == 3*c

  h4 = h+h3
  @test h4[1][1,1] == 4*a

  h5 = -h
  @test h5[3][1,1] == -c

  h6 = h/2
  @test h6[1][1,1] == a/2

  # test setting ghosts on regular history
  h_pre = deepcopy(h)
  h_post = deepcopy(h)
  set_first_ghost!(h,h_pre)
  @test h[1][1,1] == c

  set_last_ghost!(h,h_post)
  @test h[end][1,1] == a

  @test length(h) == length(h.r)+2
  @test h.r.start == 2
  @test h.r.stop == 4

end

@testset "Storage" begin

  u = ones(5,5)
  t = 0.0

  tsample = 0.05
  S = StorePlan(tsample,"state" => u, "time" => t)

  data_history = initialize_storage(S)

  store_data!(data_history,0.05,S,"state"=>deepcopy(u),"time"=>t)
  store_data!(data_history,0.1,S,"state"=>2*deepcopy(u),"time"=>2*t)
  store_data!(data_history,0.12,S,"state"=>3*deepcopy(u),"time"=>3*t)

  # test that we cannot pass the wrong type of data for one of the entries
  @test_throws AssertionError store_data!(data_history,0.15,S,"state"=>t,"time"=>t)

  # test that we cannot pass the wrong name of data for one of the entries
  @test_throws AssertionError store_data!(data_history,0.15,S,"another"=>u,"time"=>t)

  @test length(data_history) == 2

  @test data_history["state"][2] == 2*u
  @test data_history["time"][2] == 2*t

end

@testset "Interpolations" begin

  u0 = ones(5,5)
  xg, yg = 0:1:4, 0:1:4
  t = 0.0
  u = deepcopy(u0)

  t_sample = 0.05

  S = StorePlan(t_sample,"state" => u0, "time" => t)

  data_history = initialize_storage(S)

  store_data!(data_history,t,S,"state"=>deepcopy(u),"time"=>t)

  for i in 1:100
    u .+= u0 # simple model, just for making it dynamic
    t += 0.01
    store_data!(data_history,t,S,"state"=>u,"time"=>t)
  end
  tr = data_history["time"][1]:t_sample:data_history["time"][end]

  ui = interpolated_history(data_history["state"],xg,yg,tr)

  @test ui(0.5,0.4,0.54) ≈ 55.0

end


text:
# # Multiple ionospheric species
#
# This example demonstrates an ionosphere with multiple constituents---not only electrons,
# but ions as well.
#
# ## Background
#
# As discussed in the background of [Interpreting h′ and β](@ref interpreting-hpbeta),
# the propagation of VLF waves depends on the number density and collision frequency of
# charged species. Together, these two quantities describe the conductivity profile of the
# ionosphere. Not only electrons, but massive ions, influence the conductivity of the
# D-region. The influence of multiple species are combined by summing the susceptibility
# tensors for each species.
#
# ## Implementation
#
# There's nothing special about electrons in LongwaveModePropagator.
# In all of the other examples, the only ionospheric constituent is electrons, but they
# are represented as a [`Species`](@ref) type.
# LongwaveModePropagator let's you define multiple ionospheric `Species` and pass them as a
# `Tuple` or `Vector` to most functions that take the argument `species`.
#
# Let's compare electrons-only and multiple-species ionospheres.

using Plots
using LongwaveModePropagator

# The [`Species`](@ref) type defines `charge`, `mass`, `numberdensity` as a function
# of altitude, and neutral `collisionfrequency` as a function of altitude.
#
# We'll import the internally defined electron charge `QE` and mass `ME` for convenience.

using LongwaveModePropagator: QE, ME

# `QE` is literally the charge of an electron and is therefore negative.

QE

# LMP also exports the function [`waitprofile`](@ref) for the Wait and Spies (1964) electron
# density, [`electroncollisionfrequency`](@ref) for the electron-neutral collision frequency,
# and [`ioncollisionfrequency`](@ref) for the ion-neutral collision frequency.
#
# The Long Wavelength Propagation Capability (LWPC) supports built-in electrons, a positive
# ion, and a negative ion.
# Each of the ions have a mass of 58,000 electrons (approximately the mass of O₂).
#
# Here we'll use `waitprofile` for the electron number density profile, but make up a
# not-very-realistic exponential profile for the positive ion.

Ne(z) = waitprofile(z, 75, 0.32)
Np(z) = 2e6*exp(1e-4*z)
nothing  #hide

# We define density of the negative ions to conserve charge neutrality.

Nn(z) = Np(z) - Ne(z)
nothing  #hide

# For plotting, we'll replace densities below 1 m⁻³ with `NaN`.

mask(x) = x < 1 ? NaN : x

z = 0:1e3:110e3
plot(mask.(Ne.(z)), z/1000, label="Ne", linewidth=1.5)
plot!(mask.(Np.(z)), z/1000, label="Np", linewidth=1.5)
plot!(mask.(Nn.(z)), z/1000, label="Nn", linewidth=1.5)
plot!(xlabel="Density (m⁻³)", ylabel="Altitude (km)",
    xscale=:log10, legend=:topleft)
#md savefig("multiplespecies_density.png"); nothing # hide
#md # ![](multiplespecies_density.png)

# Here are the electron and ion collision frequencies:

plot(electroncollisionfrequency.(z), z/1000, label="νe", linewidth=1.5)
plot!(ioncollisionfrequency.(z), z/1000, label="νi", linewidth=1.5)
plot!(xlabel="Collision frequency (s⁻¹)", ylabel="Altitude (km)", xscale=:log10)
#md savefig("multiplespecies_collisionfreq.png"); nothing # hide
#md # ![](multiplespecies_collisionfreq.png)

# The [`Species`](@ref):

electrons = Species(QE, ME, Ne, electroncollisionfrequency)
posions = Species(abs(QE), 58000*ME, Np, ioncollisionfrequency)
negions = Species(QE, 58000*ME, Nn, ioncollisionfrequency)
nothing  #hide

# We'll compare an electrons-only and electrons-ions ionosphere.

tx = Transmitter(24e3)
rx = GroundSampler(0:5e3:2000e3, Fields.Ez)

bfield = BField(50e-6, π/2, 0)
ground = GROUND[10]

ewvg = HomogeneousWaveguide(bfield, electrons, ground)
eiwvg = HomogeneousWaveguide(bfield, (electrons, posions, negions), ground)

Ee, ae, pe = propagate(ewvg, tx, rx)
Eei, aei, pei  = propagate(eiwvg, tx, rx)

p1 = plot(rx.distance/1000, ae, label="electrons", ylabel="Amplitude (dB μV/m)")
plot!(p1, rx.distance/1000, aei, label="electrons & ions")
p2 = plot(rx.distance/1000, aei-ae,
    ylims=(-0.5, 0.5), xlabel="Range (km)", ylabel="Δ", legend=false)
plot(p1, p2, layout=grid(2,1,heights=[0.7, 0.3]))
#md savefig("multiplespecies_amplitude.png"); nothing # hide
#md # ![](multiplespecies_amplitude.png)

# The influence here is minor; the difference is hardly above the noise floor of many VLF
# receivers.
# Also, running with 3 species increases the runtime over 1 species by ~50%.


text:
function foobar(x, nrun, statistic)
    model = LinearModel
    formula = Surv ~ Dose
    yy = :Surv
    n = nrow(x)
    r0 = randn(n)
    r1 = copy(r0)
    data1 = deepcopy(x)
    t1 = zeros(nrun)
    f1 = fit(model, formula, data1)
    for i in 1:nrun
        sample!(r0, r1)
        data1[:,yy] = r0 + r1
        f1 = fit(model, formula, data1)
        t1[i] = statistic(f1)
    end
    return t1
end

@time foobar(surv, 1000, slope);
## n: 1000
## full: 0.165652 seconds (339.05 k allocations: 31.047 MB, 1.61% gc time)
## w/o sample!: 0.174539 seconds (339.05 k allocations: 31.047 MB, 1.61% gc time)
## w/o data1: 0.197007 seconds (333.05 k allocations: 30.711 MB, 1.58% gc time)
## w/o fit: 0.006245 seconds (8.38 k allocations: 417.047 KB) !!!!!

## glm general
srand(1)
X = rand(10, 2)
Y = logistic(X * [3; -3])
@time gm3 = glm(X, Y, Binomial(); dofit = false)

fit!(gm3)

function foobar2(X, Y, nrun, statistic)
    n = size(X, 1)
    r0 = randn(n)
    r1 = copy(r0)
    t1 = zeros(nrun)
    f0 = glm(X, Y, Normal(); dofit = false)
    fit!(f0)
    for i in 1:nrun
        sample!(r0, r1)
        f0.rr.y = Y + r1
        f0.fit = false
        fit!(f0)
        t1[i] = statistic(f0)
    end
    return t1
end

function foobar3(X, Y, nrun, statistic)
    n = size(X, 1)
    r0 = randn(n)
    r1 = copy(r0)
    t1 = zeros(nrun)
    f0 = glm(X, Y, Normal())
    for i in 1:nrun
        sample!(r0, r1)
        f0 = glm(X, Y + r1, Normal())
        t1[i] = statistic(f0)
    end
    return t1
end

@time foobar2(X, Y, 10000, slope);

@time foobar3(X, Y, 10000, slope);


text:
using Revise
using FGES
using Random, Statistics

Random.seed!(314)
#Generate a dataset
numFeatures = 100
numObservations = 1000
data = zeros(numObservations, numFeatures)

#number of previous features to use in caluclating the next feature
n=5

for i in 1:numFeatures
    if i ≤ n
        data[:,i] = randn(numObservations)
    else
        data[:,i] = sum(rand()*data[:,i-j] for j∈1:n) + randn(numObservations)
        #Stop the data from exploding
        data[:,i] ./= mean(data[:,i])
    end
end

g = fges(data)


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Solids (Modia3D/Solids/_module.jl)
#

#------------------------------------------------------------------------------------
# supportPoint(geo,r_abs,R_abs,e) - Return support point of geometry geo at position r_abs/R_abs in direction of unit vector e
# function is called from mpr algorithm

"""
    r = supportPoint(geo, r_abs, R_abs, e)

Return the absolute position vector `r` from world frame to *support point* of
solid `geo`, resolved in world frame in [m]. The support point is the point of solid `geo`
that is the most extreme in direction of unit vector `e`.

# Arguments
- `geo::Modia3D.AbstractSolidGeometry`: Solid geometry object.
- `r_abs::AbstractVector`: Absolute position vector of `geo` reference frame.
- `R_abs::AbstractMatrix`: Rotation matrix to rotate world frame in `geo` reference frame.
- `e::AbstractVector`: Unit vector pointing into the desired direction.
"""
supportPoint(geo::Modia3D.AbstractSolidGeometry, r_abs::AbstractVector, R_abs::AbstractMatrix, e::AbstractVector) =
            r_abs + R_abs'*centroid(geo) + R_abs'*supportPoint_abs(geo, R_abs*e) + e*geo.rsmall

supportPoint(geo::SolidSphere, r_abs::AbstractVector, R_abs::AbstractMatrix, e::AbstractVector) = r_abs + (geo.Dx/2)*e

supportPoint(geo::SolidEllipsoid, r_abs::AbstractVector, R_abs::AbstractMatrix, e::AbstractVector) =
            r_abs + R_abs'*supportPoint_abs(geo, R_abs*e)

supportPoint(geo::SolidFileMesh, r_abs::AbstractVector, R_abs::AbstractMatrix, e::AbstractVector) =
            r_abs + R_abs'*supportPoint_abs(geo, R_abs*e) # + e*geo.rsmall

#------------------------------------------------------------------------------------
# Otter: Nomenclature seems not good. Better:
#          supportPoint_abs -> supportPoint_ref
#          e_abs            -> e_ref
#
# supportPoint_abs(geo,e_abs) - Return support point of geometry geo at the reference frame of geo in direction of unit vector e_abs (resolved in geo reference frame)
# see e.g. [Gino v.d. Bergen, Collision Detection in Interactive 3D Environments, 2003]
# [G. Snethen, Xenocollide: Complex collision made simple. In Game Programming Gems 7, 2008]

# [Gino v.d. Bergen, p. 135]
function supportPoint_abs(geo::SolidEllipsoid, e_abs::AbstractVector)
   return [(geo.Dx/2)^2*e_abs[1],(geo.Dy/2)^2*e_abs[2],(geo.Dz/2)^2*e_abs[3]]/norm([(geo.Dx/2)*e_abs[1],(geo.Dy/2)*e_abs[2],(geo.Dz/2)*e_abs[3]])
end

# [Gino v.d. Bergen, p. 135]

supportPoint_abs(geo::SolidBox, e_abs::AbstractVector) = [Basics.sign_eps(e_abs[1])*(geo.Lx/2-geo.rsmall), Basics.sign_eps(e_abs[2])*(geo.Ly/2-geo.rsmall), Basics.sign_eps(e_abs[3])*(geo.Lz/2-geo.rsmall)]

# [Gino v.d. Bergen, p. 136, XenoCollide, p. 168, 169]
supportPoint_abs(geo::SolidCylinder, e_abs::AbstractVector) = norm([e_abs[1],e_abs[2]]) <= Modia3D.neps ?
                Basics.sign_eps(e_abs[3])*[0.0,0.0,geo.Lz/2] :
                Basics.sign_eps(e_abs[3])*[0.0,0.0,geo.Lz/2] + [(geo.Dx/2)*e_abs[1],(geo.Dy/2)*e_abs[2],0.0]/norm([e_abs[1],e_abs[2]])

# A. Neumayr: Cylinder + Ellipsoid as bottom and top
function supportPoint_abs(geo::SolidCapsule, e_abs::AbstractVector)
  if norm([(geo.Dx/2)*e_abs[1],(geo.Dy/2)*e_abs[2],(geo.Dy/2)*e_abs[3]]) <= Modia3D.neps
    return Basics.sign_eps(e_abs[3])*[0.0,0.0,geo.Lz/2]
  else
    return Basics.sign_eps(e_abs[3])*[0.0,0.0,geo.Lz/2] + [(geo.Dx/2)^2*e_abs[1],(geo.Dy/2)^2*e_abs[2],(geo.Dy/2)^2*e_abs[3]]/norm([(geo.Dx/2)*e_abs[1],(geo.Dy/2)*e_abs[2],(geo.Dy/2)*e_abs[3]])
  end
end

# A. Neumayr: Cuboid + arc of a circle in x dimension
function supportPoint_abs(geo::SolidBeam, e_abs::AbstractVector)
  return [Basics.sign_eps(e_abs[1])*geo.Lx/2 , Basics.sign_eps(e_abs[2])*geo.Dy/2, Basics.sign_eps(e_abs[3])*geo.Lz/2] + [ (geo.Dy/2)*Basics.sign_eps(e_abs[1]), 0.0, 0.0] #/norm([e_abs[1]])
end

# for cone: [Gino v.d. Bergen, p. 136]]
# for frustum of a cone: A. Neumayr
function supportPoint_abs(geo::SolidCone, e_abs::AbstractVector)
  radius_x = geo.Dx/2
  eta = geo.Lz/2
  sin_phi_x = radius_x/sqrt(radius_x^2 + (2*eta)^2)
  radius_y = geo.Dy/2
  sin_phi_y = radius_y/sqrt(radius_y^2 + (2*eta)^2)

  value = e_abs[3]/norm([e_abs[1],e_abs[2],e_abs[3]])

  if geo.extras[1] == 0.0 # total cone
    if      value >= sin_phi_x  && value >= sin_phi_y  # apex is support point
      return [0.0,0.0,eta]
    end
  else # frustum of a cone
    if      value >= sin_phi_x && value >= sin_phi_y && norm([e_abs[1],e_abs[2]]) <= Modia3D.neps   # apex is support point
      return [0.0,0.0,eta]
    elseif  value >= sin_phi_x && value >= sin_phi_y && norm([e_abs[1],e_abs[2]]) > Modia3D.neps # top surface is a smaller ellipse
      return ([0.0,0.0,eta]  + [(geo.Dx/2*geo.extras[1])*e_abs[1],(geo.Dy/2*geo.extras[1])*e_abs[2],0.0]/norm([e_abs[1],e_abs[2]]))
    end
  end
  if value < sin_phi_x && value < sin_phi_y && norm([e_abs[1],e_abs[2]]) > Modia3D.neps        # bottom surface is a bigger ellipse
    return ([0.0,0.0,-eta] + [(geo.Dx/2)*e_abs[1],(geo.Dy/2)*e_abs[2],0.0]/norm([e_abs[1],e_abs[2]]))
  else
    return  [0.0,0.0,-eta]
  end
end

# see Cylinder
supportPoint_abs(geo::SolidPipe, e_abs::AbstractVector) = norm([e_abs[1],e_abs[2]]) <= Modia3D.neps ?
                Basics.sign_eps(e_abs[3])*[0.0,0.0,geo.Lz/2] :
                Basics.sign_eps(e_abs[3])*[0.0,0.0,geo.Lz/2] + [(geo.Dx/2)*e_abs[1],(geo.Dy/2)*e_abs[2],0.0]/norm([e_abs[1],e_abs[2]])

# [Gino v.d. Bergen, p. 131]
function supportPoint_abs(geo::SolidFileMesh, e_abs::AbstractVector)
  # error("SolidFileMesh: function missing das haben wir noch nicht")
  max_value = Float64
  position = Float64
  for i in eachindex(geo.objPoints)
    if i == 1
      max_value = dot(geo.objPoints[i],e_abs)
      position = i
    else
      act_value = dot(geo.objPoints[i],e_abs)
      if act_value > max_value
        max_value = act_value
        position = i
      end
    end
  end
  return geo.objPoints[position]
end

#------------------------------------------------------------------------------------
# boundingBox!(..) calculates the AABB (= Axis Aligned Bounding Box) for the geometries
# it also uses support points for finding the points which are furthest away in each x,y,z - direction
# it calles supportPoint_i
"""
    boundingBox!(geo, AABB, r_abs, R_abs; tight=true, scaleFactor=0.01)

Returns the *Axis Aligned Bounding Box* of solid `geo` in argument `AABB`.

# Arguments
- `geo::Modia3D.AbstractSolidGeometry`: Solid geometry object.
- `AABB::Modia3D.BoundingBox`: On return, update AABB with the actual *Axis Aligned Bounding Box* of solid `geo`.
- `r_abs::AbstractVector`: Absolute position vector of `geo` reference frame.
- `R_abs::AbstractMatrix`: Rotation matrix to rotate world frame in `geo` reference frame.
- `tight::Bool`: If true, return the tightest AABB. If false return an AABB that is
                scaleFactor bigger than the best fitting AABB
                (for example, scaleFactor=0.1 means that the returned AABB is 10 percent bigger
                than the best fitting AABB).
- `scaleFactor::Float64`: If `tight=false`, the returned AABB is `scaleFactor` bigger than the
                          best fitting AABB.
"""
function boundingBox!(geo::Modia3D.AbstractSolidGeometry, AABB::Basics.BoundingBox, r_abs::AbstractVector, R_abs::AbstractMatrix; tight::Bool=true, scaleFactor::Float64=0.01)
  xmin = supportPoint_i(geo, r_abs[1], R_abs[:,1], -1)
  xmax = supportPoint_i(geo, r_abs[1], R_abs[:,1], +1)
  ymin = supportPoint_i(geo, r_abs[2], R_abs[:,2], -1)
  ymax = supportPoint_i(geo, r_abs[2], R_abs[:,2], +1)
  zmin = supportPoint_i(geo, r_abs[3], R_abs[:,3], -1)
  zmax = supportPoint_i(geo, r_abs[3], R_abs[:,3], +1)

  if tight    # best fitting AABB which is possible
   AABB.x_min = xmin
   AABB.x_max = xmax
   AABB.y_min = ymin
   AABB.y_max = ymax
   AABB.z_min = zmin
   AABB.z_max = zmax
 else
   # scales each axis with scaleFactor 0.1 means 10 percent bigger than best fitting AABB and so on
   @assert(scaleFactor > 0.0)
   scale = max(abs(xmax - xmin)*scaleFactor/2,abs(ymax - ymin)*scaleFactor/2,abs(zmax - zmin)*scaleFactor/2)
   AABB.x_min = xmin - scale
   AABB.x_max = xmax + scale
   AABB.y_min = ymin - scale
   AABB.y_max = ymax + scale
   AABB.z_min = zmin - scale
   AABB.z_max = zmax + scale
 end
 return AABB
 #=
  # other idea for AABB
  # scales each axis with scaleFactor 0.1 means 10 percent bigger than best fitting AABB and so on
  scaleX = abs(xmax - xmin)*scaleFactor/2
  scaleY = abs(ymax - ymin)*scaleFactor/2
  scaleZ = abs(zmax - zmin)*scaleFactor/2
  AABB.x_min = xmin - scaleX
  AABB.x_max = xmax + scaleX
  # so on

  # FileMesh
   longestEdge = max(geo.scaleFactor[1], geo.scaleFactor[2], geo.scaleFactor[3]) * geo.longestEdge * sqrt(2)
   x_min = -shape.longestEdge + r_rel[1]
   x_max =  shape.longestEdge + r_rel[1]
   y_min = -shape.longestEdge + r_rel[2]
   y_max =  shape.longestEdge + r_rel[2]
   z_min = -shape.longestEdge + r_rel[3]
   z_max = shape.longestEdge + r_rel[3]
 =#
end

function boundingBox!(geo::SolidSphere, AABB::Basics.BoundingBox, r_abs::AbstractVector, R_abs::AbstractMatrix; tight::Bool=true, scaleFactor::Float64=0.01)
  r = geo.Dx/2
  if tight    # best fitting AABB which is possible
    AABB.x_min = r_abs[1] - r
    AABB.x_max = r_abs[1] + r
    AABB.y_min = r_abs[2] - r
    AABB.y_max = r_abs[2] + r
    AABB.z_min = r_abs[3] - r
    AABB.z_max = r_abs[3] + r
  else
    # scales each axis with scaleFactor 0.1 means 10 percent bigger than best fitting AABB and so on
    @assert(scaleFactor > 0.0)
    scale = r*scaleFactor/2
    AABB.x_min = r_abs[1] - r - scale
    AABB.x_max = r_abs[1] + r + scale
    AABB.y_min = r_abs[2] - r - scale
    AABB.y_max = r_abs[2] + r + scale
    AABB.z_min = r_abs[3] - r - scale
    AABB.z_max = r_abs[3] + r + scale
  end
  return AABB
end

#------------------------------------------------------------------------------------
# supportPoint_i are utility functions for boundingBox!(..)
# it is based on supportPoint_abs (its nearly the same than supportPoint(..))
supportPoint_i(geo::Modia3D.AbstractSolidGeometry, r_absi::Float64, R_absi::AbstractVector, isign::Int) =
               r_absi + R_absi'*centroid(geo) + R_absi'*supportPoint_abs(geo, isign*R_absi) + isign*geo.rsmall

supportPoint_i(geo::SolidEllipsoid, r_absi::Float64, R_absi::AbstractVector, isign::Int) =
               r_absi + dot(R_absi, supportPoint_abs(geo, isign*R_absi))

supportPoint_i(geo::SolidFileMesh, r_absi::Float64, R_absi::AbstractVector, isign::Int) =
               r_absi + dot(R_absi, supportPoint_abs(geo, isign*R_absi))


text:
using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

@testset "2d curved block with frictionless finite sliding contact using forwarddiff" begin
    # FIXME: needs verification of some other fem software
    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_curved.med"
    mesh = aster_read_mesh(meshfile)

    upper = Problem(Elasticity, "upper", 2)
    upper.properties.formulation = :plane_stress
    upper.properties.finite_strain = true
    upper.properties.geometric_stiffness = true
    upper.elements = create_elements(mesh, "UPPER")
    update!(upper, "youngs modulus", 96.0)
    update!(upper, "poissons ratio", 1/3)

    lower = Problem(Elasticity, "lower", 2)
    lower.properties.formulation = :plane_stress
    lower.properties.finite_strain = true
    lower.properties.geometric_stiffness = true
    lower.elements = create_elements(mesh, "LOWER")
    update!(lower, "youngs modulus", 96.0)
    update!(lower, "poissons ratio", 1/3)

    bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement")
    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
    update!(bc_upper, "displacement 1", 0.0)
    update!(bc_upper, "displacement 2", -0.15)

    bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement")
    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
    update!(bc_lower, "displacement 1", 0.0)
    update!(bc_lower, "displacement 2", 0.0)

    contact = Problem(Contact, "contact between upper and lower block", 2, "displacement")
    contact.properties.rotate_normals = true
    contact.properties.finite_sliding = true
    contact.properties.friction = false
    contact.properties.use_forwarddiff = true
    contact_slave_elements = create_elements(mesh, "LOWER_TOP")
    contact_master_elements = create_elements(mesh, "UPPER_BOTTOM")
    update!(contact_slave_elements, "master elements", contact_master_elements)
    contact.elements = [contact_master_elements; contact_slave_elements]

    solver = NonlinearSolver(upper, lower, bc_upper, bc_lower, contact)
    solver()
    normu = norm(contact.assembly.u)
    info("displacement vector norm = $normu")

    # while accurate solution is unknown this is very close to linear solution
    # sqrt( ((Stress 11 - Stress 22)^2 + (Stress 22 - Stress 33)^2 + (Stress 33-Stress 11)^2 + 6*(Stress 12^2 + Stress 23^2 + Stress 13^2))/2 )
#   @test isapprox(normu, 0.49745873784105105)
    @test isapprox(normu, 0.49745872893844145)
end

#= TODO: Fix test, this is not converging
@testset "project from master to slave" begin
    el = Element(Seg2, [1, 2])
    x1 = DVTI(Vector{Float64}[
             [ 0.07406987526791842, 0.6628967239474994],
             [-0.24633092752656838, 0.4732606367688589]])
    n1 = DVTI(Vector{Float64}[
             [0.40398625635635355, 0.9147650543583191],
             [-0.5093430176405901, 0.860563588807229]])
    x2 = [0.5049198709043257, 0.27765317280577695]
    xi = project_from_master_to_slave(el, x1, n1, x2; debug=true)
    info("xi = $xi")
end
=#


text:
"""
    fixhorizon(m::Union{MDP,POMDP}, horizon::Int)

Wrap infinite horizon (PO)MDP `m` and `horizon` to the new structure creating Finite Horizon (PO)MDP.
"""
fixhorizon(m::MDP, horizon::Int) = FixedHorizonMDPWrapper(m, horizon)
fixhorizon(m::POMDP, horizon::Int) = FixedHorizonPOMDPWrapper(m, horizon)

##################
# Wrapper types
#################
struct FixedHorizonMDPWrapper{S,A,M<:MDP} <: MDP{Tuple{S, Int}, A}
    m::M
    horizon::Int
end

FixedHorizonMDPWrapper(m::MDP, d::Int) = FixedHorizonMDPWrapper{statetype(m), actiontype(m), typeof(m)}(m, d)

# TODO: Make it optional to include the stage in the observation
struct FixedHorizonPOMDPWrapper{S,A,O,M<:POMDP} <: POMDP{Tuple{S, Int}, A, Tuple{O,Int}}
    m::M
    horizon::Int
end

FixedHorizonPOMDPWrapper(m::POMDP, d::Int) = FixedHorizonPOMDPWrapper{statetype(m), actiontype(m), obstype(m), typeof(m)}(m, d)

const FHWrapper = Union{FixedHorizonMDPWrapper,FixedHorizonPOMDPWrapper}

###############################
# FiniteHorizonPOMDPs interface
###############################
HorizonLength(::Type{<:FHWrapper}) = FiniteHorizon()
horizon(w::FHWrapper) = w.horizon

###################################
# changed elements of POMDPs interface
###################################
"""
   POMDPs.states(w::HFWrapper)

Create a product of Infinite Horizon MDP's states with all stages (`1:horizon(w) + 1`).
"""
POMDPs.states(w::FHWrapper) = Iterators.product(states(w.m), 1:horizon(w)+1)

"""
    POMDPS.stateindex(w::FHWrapper, ss::Tuple{<:Any, Int})

Compute the index of the given state in Finite Horizon state space (meaning in state space of all stages).
"""
function POMDPs.stateindex(w::FHWrapper, ss::Tuple{<:Any, Int})
    s, k = ss
    return (k-1)*length(stage_states(w, 1)) + stateindex(w.m, s)
end

"""
    POMDPs.isterminal(w::FHWrapper, ss::Tuple{<:Any,Int})

Mark the state as terminal if its stage number if greater than horizon, else let the Infinite Horizon MDP's `isterminal` method decide.
"""
POMDPs.isterminal(w::FHWrapper, ss::Tuple{<:Any,Int}) = stage(w, ss) > horizon(w) || isterminal(w.m, first(ss))

"""
    POMDPs.gen(w::FHWrapper, ss::Tuple{<:Any,Int}, a, rng::AbstractRNG)

Implement the entire MDP/POMDP generative model by returning a NamedTuple.
"""
function POMDPs.gen(w::FHWrapper, ss::Tuple{<:Any,Int}, a, rng::AbstractRNG)
    out = gen(w.m, first(ss), a, rng)
    if haskey(out, :sp)
        return merge(out, (sp=(out.sp, stage(w, ss)+1),))
    else
        return out
    end
end

"""
    POMDPs.transition(w::FHWrapper, ss::Tuple{<:Any,Int}, a)

Wrap the transition result of Infinite Horizon MDP with stage number.
"""
POMDPs.transition(w::FHWrapper, ss::Tuple{<:Any,Int}, a) = InStageDistribution(transition(w.m, first(ss), a), stage(w, ss)+1)
# TODO: convert_s

"""
    POMDPs.actions(w::FHWrapper, ss::Tuple{<:Any,Int})

Return the actions of Infinite Horizon (PO)MDP.
This method assumes similar actions for all stages.
"""
POMDPs.actions(w::FHWrapper, ss::Tuple{<:Any,Int}) = actions(w.m, first(ss))

"""
    POMDPs.initialstate(w::FHWrapper)

Wrap Infinite Horizon MDP's `initialstate(mdp)` with stage 1.
"""
POMDPs.initialstate(w::FHWrapper) = InStageDistribution(initialstate(w.m), 1)

"""
    POMDPs.observations(w::FixedHorizonPOMDPWrapper)

Create a product of Infinite Horizon MDP's observations with all not-terminal stages (`1:horizon(w)`).
"""
POMDPs.observations(w::FixedHorizonPOMDPWrapper) = Iterators.product(observations(w.m), 1:horizon(w))

"""
    POMDPs.obsindex(w::FixedHorizonPOMDPWrapper, o::Tuple{<:Any, Int})::Int

Compute the index of the given observation in the Finite Horizon observation space (meaning in observation space of all stages).
"""
function POMDPs.obsindex(w::FixedHorizonPOMDPWrapper, o::Tuple{<:Any, Int})::Int
    s, k = o
    return (k-1)*length(stage_observations(w, 1)) + obsindex(w.m, s)
end

"""
    POMDPs.observation(w::FixedHorizonPOMDPWrapper[, ss::Tuple{<:Any,Int}], a, ssp::Tuple{<:Any, Int})

Create a product of Infinite Horizon MDP's observations given destination state and action (and original state) with original state's stage.
"""
POMDPs.observation(w::FixedHorizonPOMDPWrapper, ss::Tuple{<:Any,Int}, a, ssp::Tuple{<:Any, Int}) = InStageDistribution(observation(w.m, first(ss), a, first(ssp)), stage(w, ss))
POMDPs.observation(w::FixedHorizonPOMDPWrapper, a, ssp::Tuple{<:Any, Int}) = InStageDistribution(observation(w.m, a, first(ssp)), stage(w, ssp)-1)

"""
    POMDPs.initialstate(w::FHWrapper)

Return Infinite Horizon MDP's initial observations.
"""
POMDPs.initialobs(w::FixedHorizonPOMDPWrapper, ss::Tuple{<:Any,Int}) = initialobs(w.m, first(ss))
# TODO: convert_o

###############################
# FiniteHorizonPOMDPs interface
###############################
stage(w::FHWrapper, ss::Tuple{<:Any,Int}) = last(ss)
stage_states(w::FHWrapper, stage::Int) = Iterators.product(states(w.m), stage)
stage_stateindex(w::FHWrapper, ss::Tuple{<:Any,Int}) = stateindex(w.m, first(ss))
stage_observations(w::FixedHorizonPOMDPWrapper, stage::Int) = Iterators.product(observations(w.m), stage)
stage_obsindex(w::FixedHorizonPOMDPWrapper, o::Tuple{<:Any,Int}) = obsindex(w.m, first(o))
ordered_stage_states(w::FHWrapper, stage::Int) = POMDPModelTools.ordered_vector(statetype(typeof(w)), s->stage_stateindex(w,s), stage_states(w, stage), "stage_state")
ordered_stage_observations(w::FHWrapper, stage::Int) = POMDPModelTools.ordered_vector(obstype(typeof(w)), o->stage_obsindex(w,o), stage_observations(w, stage), "stage_observation")

###############################
# Forwarded parts of POMDPs interface
###############################
POMDPs.reward(w::FHWrapper, ss::Tuple{<:Any,Int}, a, ssp::Tuple{<:Any,Int}) = reward(w.m, first(ss), a, first(ssp))
POMDPs.reward(w::FixedHorizonPOMDPWrapper, ss, a, ssp, so) = reward(w.m, first(ss), a, first(ssp), first(so))
POMDPs.reward(w::FHWrapper, ss::Tuple{<:Any,Int}, a) = reward(w.m, first(ss), a)
POMDPs.actions(w::FHWrapper) = actions(w.m)
POMDPs.actionindex(w::FHWrapper, a) = actionindex(w.m, a)
POMDPs.discount(w::FHWrapper) = discount(w.m)
# TODO: convert_a

#################################
# distribution with a fixed stage
#################################
"""
    InStageDistribution{D}

Wrap given distribution with a given stage
"""
struct InStageDistribution{D}
    d::D
    stage::Int
end

"""
    distribution(d::InStageDistribution{D})::D

Return distrubution wrapped in InStageDistribution without stage
"""
function distribution(d::InStageDistribution)
    return d.d
end

"""
    stage(d::InStageDistribution)

Return stage of InStageDistribution
"""
stage(d::InStageDistribution) = d.stage

Base.rand(rng::AbstractRNG, s::Random.SamplerTrivial{<:InStageDistribution}) = (rand(rng, s[].d), s[].stage)

function POMDPs.pdf(d::InStageDistribution, ss::Tuple{<:Any, Int})
    s, k = ss
    if k == stage(d)
        return pdf(distribution(d), s)
    else
        return 0.0
    end
end

POMDPs.mean(d::InStageDistribution) = (mean(distribution(d)), stage(d))
POMDPs.mode(d::InStageDistribution) = (mode(distribution(d)), stage(d))
POMDPs.support(d::InStageDistribution) = Iterators.product(support(distribution(d)), stage(d))
POMDPs.rand(r::AbstractRNG, d::FiniteHorizonPOMDPs.InStageDistribution) = (rand(r, distribution(d)), stage(d))

#################################
# POMDPModelTools ordered actions
#################################
POMDPModelTools.ordered_actions(w::FHWrapper) = ordered_actions(w.m)