text:
module ContentTrees

export extract_tree, verify_tree, repack_tree, patch_tree

import Logging
import Pkg.TOML
import Random: randstring
import SHA
import Tar
import BSDiff

## PathNode type ##

mutable struct PathNode
    type::Symbol
    hash::NTuple{2,String}
    link::String
    copy::Union{PathNode,Nothing}
    parent::PathNode
    children::Dict{String,PathNode}
    extra::Dict{String,Any}
    function PathNode(type::Symbol)
        node = new(type)
        if type == :directory
            node.children = Dict{String,PathNode}()
        end
        return node
    end
end

function PathNode(type::Symbol, parent::PathNode)
    node = PathNode(type)
    node.parent = parent
    return node
end

function Base.show(io::IO, node::PathNode)
    print(io, "PathNode(")
    show(io, node.type)
    if isdefined(node, :hash)
        print(io, ", hash = ")
        if node.hash[1] == node.hash[2]
            show(io, node.hash[1])
        else
            show(io, node.hash)
        end
    end
    if isdefined(node, :link)
        print(io, ", link = ")
        show(io, node.link)
    end
    if isdefined(node, :children)
        print(io, ", $(length(node.children)) children")
    end
    print(io, ")")
end

## main API ##

function extract_tree(
    tarball::AbstractString,
    root::AbstractString,
    hash::Union{AbstractString, Nothing} = nothing;
    can_symlink::Union{Bool, Nothing} = nothing,
    HashType::DataType = SHA.SHA1_CTX,
)
    open(tarball) do tar
        extract_tree(tar, root, hash; can_symlink, HashType)
    end
end

function extract_tree(
    tar::IO,
    root::AbstractString,
    hash::Union{AbstractString, Nothing} = nothing;
    can_symlink::Union{Bool, Nothing} = nothing,
    HashType::DataType = SHA.SHA1_CTX,
)
    # remove destination first if it exists
    ispath(root) && @warn "path already exists, replacing" path=root

    # create tree info structure
    tree = PathNode(:directory)
    temp, can_symlink = temp_path(root, can_symlink)

    # extract tarball, recording contents
    Tar.extract(tar, temp) do hdr
        executable = hdr.type == :file && (hdr.mode & 0o100) != 0
        node = path_node!(tree, hdr.path, executable ? :executable : hdr.type)
        hdr.type == :symlink && (node.link = hdr.link)
        hdr.type != :symlink || can_symlink
    end
    resolve_symlinks!(tree)
    compute_hashes!(tree, temp; HashType)

    # simulate simlinks with copies
    !can_symlink && copy_symlinks(temp, tree)

    # verify the tree has the expected hash
    if hash !== nothing && hash ∉ tree.hash
        msg  = "Tree hash mismatch!\n"
        msg *= "  Expected: $hash\n"
        if tree.hash[1] == tree.hash[2]
            msg *= "  Computed: $(tree.hash[1])"
        else
            msg *= "  Computed: $(tree.hash[1]) (excluding empty trees)\n"
            msg *= "  Computed: $(tree.hash[2]) (including empty trees)"
        end
        chmod(temp, 0o700, recursive=true)
        rm(temp, recursive=true)
        error(msg)
    end

    # are we skipping empty trees?
    skip_empty = hash !== nothing && hash == tree.hash[2]
    @assert hash === nothing || hash == tree.hash[skip_empty+1]

    # if tree_info path exists, remove it
    tree_info_file = joinpath(temp, ".tree_info.toml")
    if haskey(tree.children, ".tree_info.toml")
        set_extra!(tree, "diffable", false)
        @warn "overwriting extracted `.tree_info.toml`" path=tree_info_file
        chmod(tree_info_file, 0o700, recursive=true)
        rm(tree_info_file, recursive=true)
    end

    # construct & write tree_info to file
    tree_info = to_toml(tree; skip_empty)
    open(tree_info_file, write=true) do io
        TOML.print(io, sorted=true, tree_info)
    end

    # move temp dir to right place
    mv(temp, root, force=true)
    return tree.hash[skip_empty+1]
end

function verify_tree(
    root::AbstractString,
    hash::Union{AbstractString, Nothing} = nothing;
    HashType::DataType = SHA.SHA1_CTX,
    tree::PathNode = tree_info(root; HashType),
)
    errors = Dict{String,Vector{String}}()
    verify_hashes!(tree, root; HashType) do node, path, msg
        haskey(errors, msg) || (errors[msg] = String[])
        push!(errors[msg], relpath(path, root))
    end
    isempty(errors) && return
    err = "Content tree $(repr(root)) has been modified:"
    for (msg, paths) in sort!(collect(errors))
        for path in sort!(paths)
            err *= "\n - $(repr(path)): $msg"
        end
    end
    error(err)
end

function repack_tree(
    tarball::AbstractString,
    root::AbstractString;
    HashType::DataType = SHA.SHA1_CTX,
    tree::PathNode = tree_info(root; HashType),
)
    open(tarball, write=true) do tar
        repack_tree(tar, root; HashType, tree)
    end
    return tarball
end

function repack_tree(
    tar::IO,
    root::AbstractString;
    HashType::DataType = SHA.SHA1_CTX,
    tree::PathNode = tree_info(root; HashType),
)
    Tar.write_tarball(tar, tree) do node, tar_path
        hdr, sys_path = node_header(root, node, tar_path)
        hdr, node.type == :directory ? node.children : sys_path
    end
    return tar
end

function patch_tree(
    old_root::AbstractString,
    new_root::AbstractString,
    patch_file::AbstractString;
    can_symlink::Union{Bool, Nothing} = nothing,
    HashType::DataType = SHA.SHA1_CTX,
    old_tree::PathNode = tree_info(old_root; HashType),
    new_hash::Union{AbstractString, Nothing} = nothing,
)
    old_data = sprint() do old_io
        repack_tree(old_io, old_root; HashType, tree=old_tree)
    end |> codeunits
    open(patch_file) do patch_io
        format = BSDiff.detect_format(patch_io, true)
        patch_type = BSDiff.patch_type(format)
        pipe = Base.BufferStream()
        @sync begin
            @async BSDiff.bspatch_core(
                patch_type,
                old_data,
                new_root, # TODO: no file is correct here
                pipe,
                patch_io,
            )
            @async extract_tree(
                pipe,
                new_root,
                new_hash;
                can_symlink,
                HashType,
            )
        end
    end
end

## loading and validating `.tree_info.toml` as PathNode tree ##

function tree_info(
    root::AbstractString;
    HashType::DataType = SHA.SHA1_CTX,
)
    # look for the `.tree_info.toml` file
    file = joinpath(root, ".tree_info.toml")
    isdir(root) || error("no directory found at $root")
    ispath(file) || error("no tree info file found at $file")

    # load data & convert to PathNode tree
    data = TOML.parsefile(file)
    tree = from_toml(data; HashType)
    resolve_symlinks!(tree)

    return tree
end

## manipulating a tree of path nodes ##

function set_extra!(node::PathNode, key::String, value::Any)
    if !isdefined(node, :extra)
        node.extra = Dict{String,Any}()
    end
    node.extra[key] = value
end

function path_node!(tree::PathNode, path::AbstractString, type::Symbol)
    node = tree
    parts = split_tar_path(path)
    for (i, name) in enumerate(parts)
        if node.type ≠ :directory
            here = join(parts[1:i-1], '/')
            error("non-directory $(repr(here)) has contents: $(repr(path))")
        end
        if i < length(parts)
            node = get!(node.children, name) do
                PathNode(:directory, node)
            end
        else
            # overwrite whatever is already there
            node = node.children[name] = PathNode(type, node)
        end
    end
    return node
end

function resolve_symlinks!(node::PathNode)
    if node.type == :symlink
        node.copy = symlink_target(node.parent, node.link)
    elseif node.type == :directory
        foreach(resolve_symlinks!, values(node.children))
    end
    if !isdefined(node, :parent)
        follow_symlinks!(node::PathNode)
    end
end

function follow_symlinks!(node::PathNode)
    if node.type == :symlink && node.copy !== nothing
        node.copy = follow_symlinks(node.copy, [node])
        if node.copy !== nothing && is_anscestor(node.copy, node)
            node.copy = nothing
        end
    elseif node.type == :directory
        foreach(follow_symlinks!, values(node.children))
    end
end

function compute_hashes!(
    node::PathNode,
    path::AbstractString;
    HashType::DataType,
)
    if node.type == :directory
        for (name, child) in node.children
            compute_hashes!(child, joinpath(path, name); HashType)
        end
    end
    node.hash = (
        git_hash(node, path; HashType, skip_empty = false),
        git_hash(node, path; HashType, skip_empty = true),
    )
end

function verify_hashes!(
    handler::Function,
    node::PathNode,
    path::AbstractString;
    HashType::DataType,
)
    # helper to invoke the error handler callback
    err(msg::AbstractString) = handler(node, path, msg)

    # check external consistency
    stat = lstat(path)
    if node.type == :directory
        if isdir(stat)
            for (name, child) in node.children
                verify_hashes!(handler, child, joinpath(path, name); HashType)
            end
        else
            err("missing directory")
        end
    elseif node.type == :symlink
        if islink(stat)
            link = readlink(path)
            if node.link != link
                err("incorrect symlink ($(repr(link)))")
            end
        elseif node.copy !== nothing
            if ispath(stat)
                verify_hashes!(handler, node.copy, path; HashType)
            else
                err("missing symlink")
            end
        elseif ispath(stat)
            err("should be a symlink")
        end
    else # file/executable
        if isfile(stat)
            hash = git_hash(node, path; HashType)
            if hash ∉ node.hash
                err("file modified ($hash)")
            end
        elseif ispath(stat)
            err("should be a file")
        else
            err("missing file")
        end
    end
end

function git_hash(
    node::PathNode,
    path::Union{AbstractString, Nothing} = nothing;
    HashType::DataType,
    skip_empty::Bool = false,
)
    if node.type == :directory
        # collect and sort children in git-tree order
        children = Pair{String,PathNode}[
            name => child for (name, child) in node.children
        ]
        let by((name, child)) = child.type == :directory ? "$name/" : name
            sort!(children; by)
        end

        # compute and return the tree hash
        return git_object_hash("tree"; HashType) do io
            for (name, child) in children
                skip_empty && is_empty_directory(child) && continue
                mode = child.type == :directory  ?  "40000" :
                       child.type == :executable ? "100755" :
                       child.type == :file       ? "100644" :
                       child.type == :symlink    ? "120000" : @assert false
                print(io, mode, ' ', name, '\0')
                write(io, hex2bytes(child.hash[skip_empty+1]))
            end
        end
    elseif node.type == :symlink
        return git_object_hash("blob"; HashType) do io
            write(io, node.link)
        end
    else # file/executable
        path === nothing &&
            error("git_hash called on file without a path")
        return git_file_hash(path; HashType)
    end
end

function git_object_hash(
    emit::Function,
    kind::AbstractString;
    HashType::DataType,
)
    ctx = HashType()
    body = codeunits(sprint(emit))
    SHA.update!(ctx, codeunits("$kind $(length(body))\0"))
    SHA.update!(ctx, body)
    return bytes2hex(SHA.digest!(ctx))
end

function git_file_hash(
    path::AbstractString;
    HashType::DataType,
    buf::Vector{UInt8} = Vector{UInt8}(undef, 4096),
)
    ctx = HashType()
    size = filesize(path)
    SHA.update!(ctx, codeunits("blob $size\0"))
    open(path) do io
        while size > 0
            n = min(size, length(buf))
            r = readbytes!(io, buf, n)
            r < n && eof(io) &&
                error("file too small (size changed during hashing)")
            v = view(buf, 1:min(r, size))
            SHA.update!(ctx, v)
            size -= length(v)
        end
        eof(io) ||
            error("file too large (size changed during hashing)")
    end
    @assert size == 0
    return bytes2hex(SHA.digest!(ctx))
end

function copy_symlinks(root::AbstractString, node::PathNode)
    if node.type == :symlink
        path = node_path(root, node)
        if node.copy === nothing
            @warn("skipping copy of broken, circular or external symlink",
                path, node.link)
        elseif !ispath(path)
            copy_symlinks(root, node.copy)
            cp(node_path(root, node.copy), path)
        end
    elseif node.type == :directory
        for child in values(node.children)
            copy_symlinks(root, child)
        end
    end
end

symlink_target(node::PathNode, path::AbstractString) =
    startswith(path, '/') ? nothing : symlink_target(node, split(path, r"/+"))

function symlink_target(node::PathNode, parts::Vector{<:AbstractString}, i::Int=1)
    if i > length(parts)
        node
    elseif parts[i] in ("", ".")
        if node.type == :directory
            symlink_target(node, parts, i+1)
        end
    elseif parts[i] == ".."
        if node.type == :directory && isdefined(node, :parent)
            symlink_target(node.parent, parts, i+1)
        end
    elseif node.type == :directory
        child = get(node.children, parts[i], nothing)
        if child !== nothing
            symlink_target(child, parts, i+1)
        end
    end
end

function follow_symlinks(node::PathNode, seen::Vector{PathNode})
    node.type != :symlink && return node
    (node ∈ seen || node.copy === nothing) && return nothing
    return follow_symlinks(node.copy, push!(seen, node))
end

node_path(node::PathNode) = something(node_path(nothing, node), ".")

function node_path(root::Union{AbstractString, Nothing}, node::PathNode)
    isdefined(node, :parent) || return root
    for (name, sibling) in node.parent.children
        sibling === node || continue
        parent = node_path(root, node.parent)
        return root !== nothing ? joinpath(parent, name) :
            parent === nothing ? name : "$parent/$name"
    end
    error("internal error: node doesn't appear in parent's children")
end

function node_header(
    root::AbstractString,
    node::PathNode,
    tar_path::AbstractString = node_path(node),
)
    sys_path = node_path(root, node)
    hdr = Tar.path_header(sys_path, tar_path)
    matching =
        node.type == :file       ? hdr.type == :file && hdr.mode == 0o644 :
        node.type == :executable ? hdr.type == :file && hdr.mode == 0o755 :
        node.type == hdr.type
    matching && return hdr, sys_path
    error("node type ($(node.type)) ≠ path type ($(hdr.type)) " *
          "for path $(repr(tar_path)) in $root")
end

is_empty_directory(node::PathNode) =
    node.type == :directory && all(is_empty_directory, values(node.children))

is_anscestor(a::PathNode, b::PathNode) =
    a === b || isdefined(b, :parent) && is_anscestor(a, b.parent)

const METADATA_KEYS = split("type link hash")

name_to_key(name::AbstractString)::String =
    name in METADATA_KEYS ? "./$name" : name

key_to_name(key::AbstractString)::String =
    startswith(key, "./") ? chop(key, head=2, tail=0) : key

function to_toml(node::PathNode; skip_empty::Bool=false)
    dict = Dict{String,Any}()
    if node.type == :directory
        if !isdefined(node, :parent) || isempty(node.children)
            dict["type"] = node.type
        end
        if !isdefined(node, :parent)
            dict["hash"] = node.hash[skip_empty+1]
        end
        for (name, child) in node.children
            dict[name_to_key(name)] = to_toml(child)
        end
        if isdefined(node, :extra) && !isempty(node.extra)
            dict["."] = copy(node.extra)
        end
    else
        dict["type"] = node.type
        if node.type == :symlink
            dict["link"] = node.link
        else
            dict["hash"] = node.hash[skip_empty+1]
        end
        isdefined(node, :extra) && merge!(dict, node.extra)
    end
    return dict
end

function from_toml(
    data::Dict{<:AbstractString};
    HashType::DataType = SHA.SHA1_CTX,
)
    type = get(data, "type", "directory")
    type in ("directory", "symlink", "file", "executable") ||
        error("invalid node type: $(repr(type))")
    type = Symbol(type)
    node = PathNode(type)
    link = get(data, "link", nothing)
    if link !== nothing
        link isa AbstractString ||
            error("invalid value for `link` key: $(repr(link))")
        node.link = link
    end
    if type == :symlink
        isdefined(node, :link) ||
            error("symlink entry without `link`: $(repr(data))")
    end
    if type == :directory
        for (key, value) in data
            key in METADATA_KEYS && continue
            if key == "."
                value isa AbstractDict ||
                    error("entry for `.` is not a dict: $(repr(value))")
                node.extra = copy(value)
            else
                child = from_toml(value)
                child.parent = node
                node.children[key_to_name(key)] = child
            end
        end
    else
        for (key, value) in  data
            key in METADATA_KEYS && continue
            set_extra!(node, key, value)
        end
    end
    if type == :file || type == :executable
        "hash" in keys(data) ||
            error("file entry without hash value: $(repr(data))")
    end
    hash = get(data, "hash", nothing)
    if hash !== nothing
        hash isa AbstractString ||
            error("invalid value for `hash` key: $(repr(hash))")
        hash = normalize_hash(hash)
    end
    if type == :file || type == :executable
        node.hash = (hash, hash)
    else
        node.hash = (
            git_hash(node; HashType, skip_empty = false),
            git_hash(node; HashType, skip_empty = true),
        )
        hash === nothing || hash in node.hash ||
            error("internally inconsistent hash: $hash ∉ $(node.hash)")
    end
    return node
end

## helper functions ##

function split_tar_path(path::AbstractString)
    path == "." && return String[]
    isempty(path) &&
        error("invalid empty tar path")
    parts = split(path, r"/+")
    isempty(parts[end]) && pop!(parts)
    for part in parts
        part in (".", "..") &&
            error("invalid tar path contains $(repr(part)): $path")
    end
    return parts
end

function temp_path(path::AbstractString, can_symlink::Union{Bool, Nothing})
    temp = "$path.$(randstring(8)).tmp"
    mkdir(temp)
    Base.Filesystem.temp_cleanup_later(temp)
    if can_symlink === nothing
        link_path = joinpath(temp, "link")
        loglevel = Logging.min_enabled_level(Logging.current_logger())
        can_symlink = try
            Logging.disable_logging(Logging.Warn)
            symlink("target", link_path)
            true
        catch err
            err isa Base.IOError || rethrow()
            false
        finally
            Logging.disable_logging(loglevel-1)
            rm(link_path; force=true)
        end
    end
    return temp, can_symlink
end

function normalize_hash(hash::AbstractString, bits::Integer=160)
    bits % 16 == 0 ||
        throw(ArgumentError("Invalid number of bits for a hash: $bits"))
    len = bits >> 2
    len_ok = length(hash) == len
    chars_ok = occursin(r"^[0-9a-f]*$"i, hash)
    if !len_ok || !chars_ok
        msg = "Hash value must be $len hexadecimal characters ($bits bits); "
        msg *= "Given hash value "
        if !chars_ok
            if isascii(hash)
                msg *= "contains non-hexadecimal characters"
            else
                msg *= "is non-ASCII"
            end
        end
        if !chars_ok && !len_ok
            msg *= " and "
        end
        if !len_ok
            msg *= "has the wrong length ($(length(hash)))"
        end
        msg *= ": $(repr(hash))"
        throw(ArgumentError(msg))
    end
    return lowercase(hash)
end

end # module


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Affine decision function #
# ========================== #
struct AffineDecisionFunction{T} <: MOI.AbstractScalarFunction
    variable_part::MOI.ScalarAffineFunction{T}
    decision_part::MOI.ScalarAffineFunction{T}
end

function AffineDecisionFunction{T}(f::MOI.SingleVariable) where T
    AffineDecisionFunction(MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(one(T), f.variable)], zero(T)),
                           convert(MOI.ScalarAffineFunction{T}, zero(T)))
end

function AffineDecisionFunction{T}(f::SingleDecision) where T
    AffineDecisionFunction(convert(MOI.ScalarAffineFunction{T}, zero(T)),
                           MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(one(T), f.decision)], zero(T)))
end

# Affine decision vector function #
# ========================== #
struct VectorAffineDecisionFunction{T} <: MOI.AbstractVectorFunction
    variable_part::MOI.VectorAffineFunction{T}
    decision_part::MOI.VectorAffineFunction{T}
end
MOI.output_dimension(f::VectorAffineDecisionFunction) = MOI.output_dimension(f.variable_part)

# Base overrides #
# ========================== #
function Base.copy(f::F) where F <: Union{AffineDecisionFunction,
                                          VectorAffineDecisionFunction}
    return F(copy(f.variable_part),
             copy(f.decision_part))
end

function Base.convert(::Type{AffineDecisionFunction{T}}, α::T) where T
    return AffineDecisionFunction(convert(MOI.ScalarAffineFunction{T}, α),
                                  convert(MOI.ScalarAffineFunction{T}, zero(T)))
end

function Base.convert(::Type{AffineDecisionFunction{T}},
                      f::MOI.SingleVariable) where T
    return AffineDecisionFunction(convert(MOI.ScalarAffineFunction{T}, f),
                                  convert(MOI.ScalarAffineFunction{T}, zero(T)))
end

function Base.convert(::Type{AffineDecisionFunction{T}},
                      f::SingleDecision) where T
    return AffineDecisionFunction{T}(f)
end

function Base.convert(::Type{AffineDecisionFunction{T}},
                      f::MOI.ScalarAffineFunction{T}) where T
    return AffineDecisionFunction(f, convert(MOI.ScalarAffineFunction{T}, zero(T)))
end

function Base.convert(::Type{AffineDecisionFunction{T}},
                      f::AffineDecisionFunction{T}) where T
    return f
end

function Base.convert(::Type{AffineDecisionFunction{T}},
                      f::AffineDecisionFunction) where T
    return AffineDecisionFunction(convert(MOI.ScalarAffineFunction{T}, f.variable_part),
                                  convert(MOI.ScalarAffineFunction{T}, f.decision_part))
end

function Base.convert(::Type{MOI.SingleVariable}, f::AffineDecisionFunction)
    if !iszero(f.variable_part.constant) || !isone(length(f.variable_part.terms)) || !isone(f.variable_part.terms[1].coefficient)
        throw(InexactError(:convert, MOI.SingleVariable, f))
    end
    return MOI.SingleVariable(f.variable_part.terms[1].variable_index)
end

function Base.convert(::Type{SingleDecision}, f::AffineDecisionFunction)
    if !iszero(f.decision_part.constant) || !isone(length(f.decision_part.terms)) || !isone(f.decision_part.terms[1].coefficient)
        throw(InexactError(:convert, SingleDecision, f))
    end
    return SingleDecision(f.decision_part.terms[1].variable_index)
end

function Base.convert(::Type{VectorAffineDecisionFunction{T}}, α::T) where T
    return VectorAffineDecisionFunction(MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], [α]),
                                        MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], [zero(T)]))
end

function Base.convert(::Type{VectorAffineDecisionFunction{T}},
                      f::MOI.SingleVariable) where T
    return VectorAffineDecisionFunction(
        MOI.VectorAffineFunction{T}([MOI.VectorAffineTerm{T}(1, MOI.ScalarAffineTerm(one(T), f.variable))], [zero(T)]),
        MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], [zero(T)]))
end

function Base.convert(::Type{VectorAffineDecisionFunction{T}},
                      f::SingleDecision) where T
    return VectorAffineDecisionFunction(
        MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], [zero(T)]),
        MOI.VectorAffineFunction{T}([MOI.VectorAffineTerm{T}(1, MOI.ScalarAffineTerm(one(T), f.decision))], [zero(T)]))
end

function Base.convert(::Type{VectorAffineDecisionFunction{T}},
                      f::MOI.VectorOfVariables) where T
    n = length(f.variables)
    terms = map(1:n) do i
        MOI.VectorAffineTerm{T}(i, MOI.ScalarAffineTerm(one(T), f.variables[i]))
    end
    variable_part = MOI.VectorAffineFunction{T}(
        terms,
        zeros(T, n))
    return VectorAffineDecisionFunction(variable_part,
                                        MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], zeros(T, n)))
end

function Base.convert(::Type{VectorAffineDecisionFunction{T}},
                      f::VectorOfDecisions) where T
    n = length(f.decisions)
    terms = map(1:n) do i
        MOI.VectorAffineTerm{T}(i, MOI.ScalarAffineTerm(one(T), f.decisions[i]))
    end
    decision_part = MOI.VectorAffineFunction{T}(
        terms,
        zeros(T, n))
    return VectorAffineDecisionFunction(MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], zeros(T, n)),
                                        decision_part)
end

function Base.convert(::Type{VectorAffineDecisionFunction{T}},
                      f::MOI.ScalarAffineFunction{T}) where T
    variable_part = MOI.VectorAffineFunction{T}([MOI.VectorAffineTerm{T}(1, t) for t in f.terms],
                                                [f.constant])
    return VectorAffineDecisionFunction(variable_part,
                                        MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], [zero(T)]))
end

function Base.convert(::Type{VectorAffineDecisionFunction{T}},
                      f::AffineDecisionFunction{T}) where T
    variable_part = MOI.VectorAffineFunction{T}([MOI.VectorAffineTerm{T}(1, t) for t in f.variable_part.terms],
                                                [f.variable_part.constant])
    decision_part = MOI.VectorAffineFunction{T}([MOI.VectorAffineTerm{T}(1, t) for t in f.decision_part.terms],
                                                [f.decision_part.constant])
    return VectorAffineDecisionFunction(variable_part,
                                        decision_part)
end

function Base.convert(::Type{VectorAffineDecisionFunction{T}},
                      f::MOI.VectorAffineFunction{T}) where T
    n = MOI.output_dimension(f)
    return VectorAffineDecisionFunction(f,
                                        MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], zeros(T, n)))
end

function Base.zero(F::Type{AffineDecisionFunction{T}}) where T
    return convert(F, zero(T))
end

function Base.iszero(f::AffineDecisionFunction)
    return iszero(MOI.constant(f)) && MOIU._is_constant(MOIU.canonical(f))
end

function Base.isone(f::AffineDecisionFunction)
    return isone(MOI.constant(f)) && MOIU._is_constant(MOIU.canonical(f))
end

function Base.one(F::Type{AffineDecisionFunction{T}}) where T
    return convert(F, one(T))
end

function MOI._dicts(f::Union{AffineDecisionFunction, VectorAffineDecisionFunction})
    return (MOI.sum_dict(MOI.term_pair.(f.variable_part.terms)),
            MOI.sum_dict(MOI.term_pair.(f.decision_part.terms)))
end

function Base.isapprox(f::F, g::G; kwargs...) where {
    F<:Union{
        AffineDecisionFunction,
        VectorAffineDecisionFunction,
    },
    G<:Union{
        AffineDecisionFunction,
        VectorAffineDecisionFunction,
    },
}
    return isapprox(MOI.constant(f), MOI.constant(g); kwargs...) && all(
        MOI.dict_compare.(
            MOI._dicts(f),
            MOI._dicts(g),
            (α, β) -> isapprox(α, β; kwargs...),
        ),
    )
end

# JuMP overrides #
# ========================== #
function AffineDecisionFunction(aff::DAE)
    JuMP._assert_isfinite(aff)
    decision_terms = Vector{MOI.ScalarAffineTerm{Float64}}()
    for (coef, dvar) in linear_terms(aff.decisions)
        # Any fixed decision value is set in the decision bridge
        push!(decision_terms, MOI.ScalarAffineTerm(coef, index(dvar)))
    end
    return AffineDecisionFunction(JuMP.moi_function(aff.variables),
                                  MOI.ScalarAffineFunction(decision_terms, 0.0))
end
JuMP.moi_function(aff::DAE) = AffineDecisionFunction(aff)
function JuMP.moi_function_type(::Type{DecisionAffExpr{T}}) where T
    return AffineDecisionFunction{T}
end

is_decision_type(::Type{<:AffineDecisionFunction}) = true

AffineDecisionFunction(aff::_DAE) = AffineDecisionFunction(convert(DAE, aff))

function VectorAffineDecisionFunction(affs::Vector{DAE})
    # Decision part
    dlength = sum(aff -> length(linear_terms(aff.decisions)), affs)
    decision_terms = Vector{MOI.VectorAffineTerm{Float64}}(undef, dlength)
    offset = 0
    for (i, aff) in enumerate(affs)
        offset = JuMP._fill_vaf!(decision_terms, offset, i, aff.decisions)
    end
    VectorAffineDecisionFunction(JuMP.moi_function([aff.variables for aff in affs]),
                                 MOI.VectorAffineFunction(decision_terms, zeros(length(affs))))
end
JuMP.moi_function(affs::Vector{<:DecisionAffExpr}) = VectorAffineDecisionFunction(affs)
function JuMP.moi_function_type(::Type{<:Vector{<:DecisionAffExpr{T}}}) where {T}
    return VectorAffineDecisionFunction{T}
end

is_decision_type(::Type{<:VectorAffineDecisionFunction}) = true

JuMP.moi_function(aff::_DAE) = AffineDecisionFunction(aff)
function JuMP.moi_function_type(::Type{_DecisionAffExpr{T}}) where T
    return AffineDecisionFunction{T}
end

function _DecisionAffExpr(m::Model, f::MOI.ScalarAffineFunction)
    aff = _DAE()
    for t in f.terms
        JuMP.add_to_expression!(aff, t.coefficient, DecisionRef(m, t.variable_index))
    end
    # There should be not any constants in the decision terms
    return aff
end

function DAE(model::Model, f::AffineDecisionFunction)
    return DAE(AffExpr(model, f.variable_part),
               _DecisionAffExpr(model, f.decision_part))
end

function JuMP.jump_function_type(::Model,
                                 ::Type{AffineDecisionFunction{T}}) where T
    return DecisionAffExpr{T}
end
function JuMP.jump_function(model::Model, f::AffineDecisionFunction{T}) where T
    return DecisionAffExpr{T}(model, f)
end

function JuMP.jump_function_type(::Model,
                                 ::Type{VectorAffineDecisionFunction{T}}) where T
    return Vector{DecisionAffExpr{T}}
end
function JuMP.jump_function(model::Model, f::VectorAffineDecisionFunction{T}) where T
    return DecisionAffExpr{T}[
        DecisionAffExpr{T}(model, f) for f in MOIU.eachscalar(f)]
end

# MOI Function interface #
# ========================== #
MOI.constant(f::AffineDecisionFunction) = MOI.constant(f.variable_part)
MOI.constant(f::AffineDecisionFunction, T::DataType) = MOI.constant(f)
MOI.constant(f::VectorAffineDecisionFunction) = MOI.constant(f.variable_part)
MOI.constant(f::VectorAffineDecisionFunction, T::DataType) = MOI.constant(f)

function MOIU.eval_variables(varval::Function, f::AffineDecisionFunction)
    var_value = MOIU.eval_variables(varval, f.variable_part)
    dvar_value = MOIU.eval_variables(varval, f.decision_part)
    return var_value + dvar_value
end
function MOIU.eval_variables(varval::Function, f::VectorAffineDecisionFunction)
    var_value = MOIU.eval_variables(varval, f.variable_part)
    dvar_value = MOIU.eval_variables(varval, f.decision_part)
    return var_value + dvar_value
end

function MOIU.map_indices(index_map::Function, f::Union{AffineDecisionFunction, VectorAffineDecisionFunction})
    return typeof(f)(MOIU.map_indices(index_map, f.variable_part),
                     MOIU.map_indices(index_map, f.decision_part))
end

Base.eltype(it::MOIU.ScalarFunctionIterator{VectorAffineDecisionFunction{T}}) where T = AffineDecisionFunction{T}

function MOIU.ScalarFunctionIterator(f::VectorAffineDecisionFunction)
    return MOIU.ScalarFunctionIterator(
        f,
        (
            MOIU.output_index_iterator(f.variable_part.terms, MOI.output_dimension(f)),
            MOIU.output_index_iterator(f.decision_part.terms, MOI.output_dimension(f))
        ),
    )
end

function Base.getindex(it::MOIU.ScalarFunctionIterator{VectorAffineDecisionFunction{T}}, output_index::Integer) where T
    variable_part = MOI.ScalarAffineFunction{T}(
        MOI.ScalarAffineTerm{T}[
            it.f.variable_part.terms[i].scalar_term for
            i in MOIU.ChainedIteratorAtIndex(it.cache[1], output_index)
        ],
        it.f.variable_part.constants[output_index],
    )
    decision_part = MOI.ScalarAffineFunction{T}(
        MOI.ScalarAffineTerm{T}[
            it.f.decision_part.terms[i].scalar_term for
            i in MOIU.ChainedIteratorAtIndex(it.cache[2], output_index)
        ],
        it.f.decision_part.constants[output_index],
    )
    AffineDecisionFunction(variable_part, decision_part)
end
function Base.getindex(it::MOIU.ScalarFunctionIterator{VectorAffineDecisionFunction{T}}, output_indices::AbstractVector) where T
    variable_terms = MOI.VectorAffineTerm{T}[]
    decision_terms = MOI.VectorAffineTerm{T}[]
    for (i, output_index) in enumerate(output_indices)
        for j in MOIU.ChainedIteratorAtIndex(it.cache[1], output_index)
            push!(variable_terms, MOI.VectorAffineTerm(i, it.f.variable_part.terms[j].scalar_term))
        end
        for j in MOIU.ChainedIteratorAtIndex(it.cache[2], output_index)
            push!(decision_terms, MOI.VectorAffineTerm(i, it.f.decision_part.terms[j].scalar_term))
        end
    end
    return VectorAffineDecisionFunction(MOIU.VAF(variable_terms, it.f.variable_part.constants[output_indices]),
                                        MOIU.VAF(decision_terms, it.f.decision_part.constants[output_indices]))
end

function MOIU.zero_with_output_dimension(::Type{VectorAffineDecisionFunction{T}}, n::Integer) where T
    return MOI.VectorAffineDecisionFunction{T}(MOIU.zero_with_output_dimension(MOI.VAF{T}, n),
                                               MOIU.zero_with_output_dimension(MOI.VAF{T}, n))
end

function MOIU.substitute_variables(variable_map::Function, f::AffineDecisionFunction{T}) where T
    g = AffineDecisionFunction(convert(MOI.ScalarAffineFunction{T}, MOI.constant(f.variable_part)),
                               convert(MOI.ScalarAffineFunction{T}, zero(T)))
    # Substitute variables
    for term in f.variable_part.terms
        func::AffineDecisionFunction{T} = variable_map(term.variable_index)
        new_term = MOIU.operate(*, T, term.coefficient, func)::AffineDecisionFunction{T}
        MOIU.operate!(+, T, g, new_term)::AffineDecisionFunction{T}
    end
    # Substitute decisions
    for term in f.decision_part.terms
        func::AffineDecisionFunction{T} = variable_map(term.variable_index)
        new_term = MOIU.operate(*, T, term.coefficient, func)::AffineDecisionFunction{T}
        MOIU.operate!(+, T, g, new_term)::AffineDecisionFunction{T}
    end
    return g
end

function MOIU.substitute_variables(variable_map::Function, f::VectorAffineDecisionFunction{T}) where T
    n = MOI.output_dimension(f)
    g = VectorAffineDecisionFunction(MOI.VectorAffineFunction(MOI.VectorAffineTerm{T}[], copy(MOI.constant(f.variable_part))),
                                     MOI.VectorAffineFunction(MOI.VectorAffineTerm{T}[], zeros(T, n)))
    # Substitute variables
    for term in f.variable_part.terms
        func::AffineDecisionFunction{T} = variable_map(term.scalar_term.variable_index)
        new_term = MOIU.operate(*, T, term.scalar_term.coefficient, func)::AffineDecisionFunction{T}
        MOIU.operate_output_index!(+, T, term.output_index, g, new_term)::typeof(g)
    end
    # Substitute decisions
    for term in f.decision_part.terms
        func::AffineDecisionFunction{T} = variable_map(term.scalar_term.variable_index)
        new_term = MOIU.operate(*, T, term.scalar_term.coefficient, func)::AffineDecisionFunction{T}
        MOIU.operate_output_index!(+, T, term.output_index, g, new_term)::VectorAffineDecisionFunction{T}
    end
    return g
end

MOIU.constant_vector(f::AffineDecisionFunction) = [f.variable_part.constant]
MOIU.constant_vector(f::VectorAffineDecisionFunction) = f.variable_part.constants

MOIU.scalar_type(::Type{VectorAffineDecisionFunction{T}}) where T = AffineDecisionFunction{T}

MOIU.is_canonical(f::Union{AffineDecisionFunction, VectorAffineDecisionFunction}) =
    MOIU.is_canonical(f.variable_part) &&
    MOIU.is_canonical(f.decision_part)

MOIU.canonical(f::Union{AffineDecisionFunction, VectorAffineDecisionFunction}) = MOIU.canonicalize!(copy(f))

function MOIU.canonicalize!(f::Union{AffineDecisionFunction, VectorAffineDecisionFunction})
    MOIU.canonicalize!(f.variable_part)
    MOIU.canonicalize!(f.decision_part)
    return f
end

MOIU._is_constant(f::AffineDecisionFunction) =
    MOIU._is_constant(f.variable_part) &&
    MOIU._is_constant(f.decision_part)

MOIU.all_coefficients(p::Function, f::AffineDecisionFunction) =
    MOIU.all_coefficients(p, f.variable_part) &&
    MOIU.all_coefficients(p, f.decision_part)

MOIU.isapprox_zero(f::AffineDecisionFunction, tol) =
    MOIU.isapprox_zero(f.variable_part, tol) &&
    MOIU.isapprox_zero(f.decision_part, tol)

function MOIU.filter_variables(keep::Function, f::Union{AffineDecisionFunction, VectorAffineDecisionFunction})
    # Only filter variable part and decision part
    return typeof(f)(MOIU.filter_variables(keep, f.variable_part),
                     MOIU.filter_variables(keep, f.decision_part))
end

function MOIU.operate_coefficients(f, func::Union{AffineDecisionFunction, VectorAffineDecisionFunction})
    return typeof(func)(MOIU.operate_coefficients(f, func.variable_part),
                        MOIU.operate_coefficients(f, func.decision_part))
end

function MOIU.vectorize(funcs::AbstractVector{AffineDecisionFunction{T}}) where T
    return VectorAffineDecisionFunction{T}(
        MOIU.vectorize([f.variable_part for f in funcs]),
        MOIU.vectorize([f.decision_part for f in funcs]))
end

function MOIU.scalarize(f::VectorAffineDecisionFunction{T}, ignore_constants::Bool = false) where T
    variable_part = MOIU.scalarize(f.variable_part, ignore_constants)
    decision_part = MOIU.scalarize(f.decision_part, ignore_constants)
    return map(zip(variable_part, decision_part)) do (var_part, dvar_part)
        AffineDecisionFunction(var_part, dvar_part)
    end
end

function MOIU.convert_approx(::Type{SingleDecision},
                             func::AffineDecisionFunction{T};
                             tol = MOIU.tol_default(T)) where {T}
    f = MOIU.canonical(func)
    i = findfirst(t -> isapprox(t.coefficient, one(T), atol = tol), f.decision_part.terms)
    if abs(MOI.constant(f)) > tol ||
        i === nothing ||
        any(j -> j != i && abs(f.decision_part.terms[j].coefficient) > tol, eachindex(f.decision_part.terms)) ||
        any(j -> abs(f.variable_part.terms[j].coefficient) > tol, eachindex(f.variable_part.terms))
        throw(InexactError(:convert_approx, SingleDecision, func))
    end
    return SingleDecision(f.decision_part.terms[i].variable_index)
end

function MOIU.convert_approx(::Type{VectorOfDecisions},
                             func::VectorAffineDecisionFunction{T};
                             tol = MOIU.tol_default(T)) where {T}
    return VectorOfDecisions([
        MOIU.convert_approx(SingleDecision, f, tol = tol).decision for
        f in MOIU.scalarize(func)
    ])
end

function modify_coefficient!(terms::Vector{MOI.ScalarAffineTerm{T}},
                             index::MOI.VariableIndex,
                             new_coefficient::Number) where T
    i = something(findfirst(t -> t.variable_index == index,
                            terms), 0)
    if iszero(i)
        # The variable was not already included in the terms
        if !iszero(new_coefficient)
            # Add it
            push!(terms,
                  MOI.ScalarAffineTerm(T(new_coefficient), index))
        end
    else
        # The variable is included in the terms
        if iszero(new_coefficient)
            # Remove it
            deleteat!(terms, i)
        else
            # Update coefficient
            terms[i] = MOI.ScalarAffineTerm(T(new_coefficient), index)
        end
    end
    return nothing
end

function modify_coefficients!(terms::Vector{MOI.VectorAffineTerm{T}},
                              index::MOI.VariableIndex,
                              new_coefficients::AbstractVector) where T
    rowmap = Dict(c[1]=>i for (i,c) in enumerate(new_coefficients))
    del = Int[]
    for i in findall(t -> t.variable_index == index, terms)
        row = terms[i].output_index
        j = Base.get(rowmap, row, 0)
        if !iszero(j)
            if iszero(new_coefficients[j][2])
                push!(del, i)
            else
                terms[i] = MOI.VectorAffineTerm(row, MOI.ScalarAffineTerm(new_coefficients[j][2], index))
            end
            rowmap[row] = 0
        end
    end
    deleteat!(terms, del)
    for (row, j) in rowmap
        new_coefficient = new_coefficients[j][2]
        if !iszero(j) && !iszero(new_coefficient)
            push!(terms, MOI.VectorAffineTerm(row, MOI.ScalarAffineTerm(new_coefficient, index)))
        end
    end
    return nothing
end

function add_term!(terms::Vector{MOI.ScalarAffineTerm{T}},
                   term::MOI.ScalarAffineTerm{T}) where T
    index = term.variable_index
    coefficient = term.coefficient
    i = something(findfirst(t -> t.variable_index == index,
                            terms), 0)
    if iszero(i)
        if !iszero(coefficient)
            # The variable was not already included in the terms
            push!(terms, MOI.ScalarAffineTerm(coefficient, index))
        end
    else
        # Variable already included, increment the coefficient
        new_coeff = terms[i].coefficient + coefficient
        if iszero(new_coeff)
            deleteat!(terms, i)
        else
            terms[i] = MOI.ScalarAffineTerm(new_coeff, index)
        end
    end
    return nothing
end

function remove_term!(terms::Vector{MOI.ScalarAffineTerm{T}},
                      term::MOI.ScalarAffineTerm{T}) where T
    index = term.variable_index
    coefficient = term.coefficient
    i = something(findfirst(t -> t.variable_index == index,
                            terms), 0)
    if iszero(i) && !iszero(coefficient)
        # The variable was not already included in the terms
        push!(terms, MOI.ScalarAffineTerm(-coefficient, index))
    else
        # Variable already included, increment the coefficient
        new_coeff = terms[i].coefficient - coefficient
        if iszero(new_coeff)
            deleteat!(terms, i)
        else
            terms[i] = MOI.ScalarAffineTerm(new_coeff, index)
        end
    end
    return nothing
end

ext: jl


text:
using Test, DiscreteExteriorCalculus, DiscretePDEs
using UniqueVectors: UniqueVector
using AdmittanceModels: lossless_modes_dense, apply_transform

geo_code = GeoCode(gdspy_open(joinpath(@__DIR__, "cpw_resonator.gds"), "main"),
    layer_translation=Dict(23 => 50),
    layer_extrusion=Dict(21 => 10, 22 => 10, 23 => -100, 24 => 50),
    layer_volume_difference=[(24, 22, true), (23, 24, true), (23, 22, true)],
    layer_volume_groups=Dict("Superconductor" => [22], "Substrate" => [23], "Vacuum" => [24]),
    layer_mesh_size=Dict(21 => 10))
file_name = joinpath(@__DIR__, "cpw_resonator.geo")
geo_write!(file_name, geo_code, characteristic_length_factor=1.5, random_factor=1e-8,
    random_factor_3D=1e-8)

initialize!()
gmsh_open!(file_name)
N, K = 3, 4
mesh!(K)

node_tags, points, tcomp = get_triangulated_complex(N, K, 1e-6)
group_dict = get_physical_groups(node_tags, points)
comp = tcomp.complex
orient!(comp)
m = Metric(N)
mesh = Mesh(tcomp, circumcenter(m))

μ⁻_form = get_material(comp, 1/DiscretePDEs.μ₀, 3)
Λ⁻_form = get_material(comp, 0, 2)
σ_form = get_material(comp, 0, 2)
ϵr = 11.9 # silicon
ϵ_form = DiscretePDEs.ϵ₀ * (get_material(comp, 1, 2) +
    get_material(comp, group_dict["Substrate"], ϵr-1, 2))

pso, null_basis = electrodynamics_pso(m, mesh, Vector{Cell{N}}[],
    append!(subcomplex(comp, group_dict["Superconductor"]), boundary(comp)),
    μ⁻_form, Λ⁻_form, σ_form, ϵ_form)
constrained_pso = apply_transform(pso, null_basis)

λs, vs = lossless_modes_dense(constrained_pso, min_freq=1e9)
freqs = imag.(λs)/(2π)

hfss_answers = [47.0946, 140.140] * 1e9
@testset "cpw resonator" begin
    @test isapprox(freqs[1:2], hfss_answers, rtol=1.5e-1)
end

if false # set to true to display plots with Gmsh
    comp_points = UniqueVector([c.points[1] for c in comp.cells[1]])
    ordering = [findfirst(isequal(p), comp_points) for p in points]
    for i in 1:2
        v = null_basis * vs[:,i]
        v /= maximum(abs.(v))
        vec_A = sharp(m, comp, v)
        add_field!("Vector potential mode $i", node_tags, vec_A[ordering])
    end
    gui!()
end
rm(file_name)


text:
module Test_Signal4Assemblies

using  Modia3D
using  Modia3D.StaticArrays
import Modia3D.ModiaMath

vgreen  = Modia3D.Material(color="Green", transparency=0.5)
vyellow = Modia3D.Material(color="Yellow", transparency=0.5)
vred    = Modia3D.Material(color="Red", transparency=0.5)
cmat    = "Steel"

Dx = 0.2
Lx = 1.0
Ly = 0.3
Lz = 0.3

# a few signal macros
@signal SinWave(;A=5.0, w=2.0) begin
   y = ModiaMath.RealScalar("y", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::SinWave, sim::ModiaMath.SimulationState)
    signal.y.value = signal.A*sin(signal.w*sim.time)
end

@signal CosWave(;A=5.0, w=2.0) begin
   y = ModiaMath.RealScalar("y", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::CosWave, sim::ModiaMath.SimulationState)
    signal.y.value = signal.A*cos(sim.time)
end

@signal LinearSig(;A=5.0, B=2.0) begin
   y = ModiaMath.RealScalar("y", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::LinearSig, sim::ModiaMath.SimulationState)
    signal.y.value = signal.A*sim.time + signal.B
end

# 3 sub assemblies
@assembly greenAssembly() begin
  green     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vgreen),[ [-Lx/2,0,0],[ Lx/2,0,0] ])
  green2    = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vgreen),[ [-Lx/2,0,0],[ Lx/2,0,0] ])
  rev12     = Modia3D.Revolute(green.frames[2], green2.frames[1])

  sinWave     = SinWave()
  cosWave     = CosWave()
  sigGreen    = Modia3D.SignalToFlangeAngle(sinWave.y)
  sigGreen2   = Modia3D.SignalToFlangeAngle(cosWave.y)
  Modia3D.connect(sigGreen2, rev12)
end

@assembly yellowAssembly() begin
  gelb    = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vyellow),[ [-Lx/2,0,0],[ Lx/2,0,0] ])

  linearSig   = LinearSig()
  sigGelb     = Modia3D.SignalToFlangeAngle(linearSig.y)
end

# main assembly
@assembly Pendulum() begin
world    = Modia3D.Object3D()
world_f1 = Modia3D.Object3D(world , r=[0.5*Lx, 0.0, Lz/2])
gruen    = greenAssembly()
gelb     = yellowAssembly()

rev1     = Modia3D.Revolute(world, gruen.green.frames[1])
rev2     = Modia3D.Revolute(world_f1, gelb.gelb.frames[1])

Modia3D.connect(gruen.sigGreen, rev1)
Modia3D.connect(gelb.sigGelb,   rev2)
end

pendulum = Pendulum(sceneOptions=Modia3D.SceneOptions(visualizeFrames=true,defaultFrameLength=0.3, enableContactDetection=false))

model    = Modia3D.SimulationModel(pendulum, analysis=ModiaMath.KinematicAnalysis)
# ModiaMath.print_ModelVariables(model)
result   = ModiaMath.simulate!(model, stopTime=5.0)
ModiaMath.plot(result, ["gruen.sinWave.y", "rev1.phi", "gruen.cosWave.y", "gruen.rev12.phi", "gelb.linearSig.y", "rev2.phi"] )

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


text:
"""

# StructuralCausalModels

StructuralCausalModels.jl provides functionality to analyse directed acyclic
graph (DAG) based causal models as described in `StatisticalRethinking`, `Causal
Inference in Statistics` and `Cause and Correlation in Biology`.

My initial goal for this package is to have a way to apply SCM ideas to the
examples in the context of StatisticalRethinking, i.e. a working version of
`basis_set()`, `d_separation()`, `m_separation()` and `adjustment_sets()`.

SCM can be used as an alias to StructuralCausalModels.
"""
module StructuralCausalModels

using Reexport

@reexport using DataStructures, LinearAlgebra, Statistics
@reexport using Distributions, NamedArrays, DataFrames, CSV
@reexport using Combinatorics

src_path = @__DIR__

using DocStringExtensions: SIGNATURES, FIELDS, TYPEDEF

const scm_src_path = @__DIR__

"""

SCM

Alias for StructuralCausalModels
"""
const SCM = StructuralCausalModels

"""

# scm_path

Relative path using the StructuralCausalModels.jl src/ directory.

### Example to get access to the data subdirectory
```julia
scm_path("..", "data")
```

Part of the API, exported.
"""
scm_path(parts...) = normpath(joinpath(scm_src_path, parts...))

SymbolList = Union{Symbol, Vector{Symbol}}
SymbolListOrNothing = Union{SymbolList, Nothing}
OrderedDictOrNothing = Union{OrderedDict, Nothing}
NamedArrayOrNothing = Union{NamedArray, Nothing}
DataFrameOrNothing = Union{DataFrame, Nothing}
ModelDefinition = Union{OrderedDict, AbstractString, NamedArray}

include("types/DAG.jl")
include("types/Path.jl")
include("types/BasisSet.jl")
include("types/ConditionalIndependence.jl")

include("methods/dag_methods.jl")
include("methods/basis_set.jl")
include("methods/shipley_test.jl")
include("methods/d_separation.jl")
include("methods/pcor.jl")
include("methods/pcor_test.jl")
include("methods/induced_covariance_graph.jl")
include("methods/all_paths.jl")
include("methods/open_paths.jl")
include("methods/backdoor_paths.jl")
include("methods/adjustment_sets.jl")
include("methods/graph_utils.jl")
include("methods/ancestral_graph.jl")
include("methods/ribbon_graph.jl")
include("methods/maximize.jl")
include("methods/m_separation.jl")

include("utils/show_dag_path.jl")
include("utils/ggm_conversions.jl")
include("utils/dagitty_conversions.jl")
include("utils/graphviz_conversions.jl")

export
  SCM,
  scm_path

end # module


text:
# using TraitWrappers

# abstract type AbstractJLBoostPredictTraitWrapper <: AbstractTraitWrapper end
#
# struct JLBoostHasPredictTraitWrapper <: AbstractJLBoostPredictTraitWrapper
#     object::T
# end
#
# struct IterableTraitWrapper{T} <: AbstractTraitWrapper
#     object::T
# end
#
# struct ColumnBangAccessible{T} <: AbstractTraitWrapper
#     object::T
#     ColumnBangAccessible(t::T) = new{T}(t)
# end
#
# # checks whether T satisfies the IterableTraitWrapper
# satisfy(::Type{T}) where T = begin
#     hasmethod(iterate, Tuple{T})
# end
#
# satisfy(::Type{W}) where W <: AbstractVector = true
#
# abstract type ABC end
#
# ABC(x) = begin
#     satisfy(ABC, x)
# end
#
# satisfy(::Type{<:ABC}, obj) = begin
#     println("ok")
#     true
# end
#
# struct DEF <: ABC
#     object
# end
#
# DEF(1)
#
# satisfy(DEF, 1)


text:
# Function for placing bounds on a variable:

function log_bounds_upper(x::T,x1::Float64,x2::Float64) where {T <: Real}
@assert(x2 > x1)
# Places bounds on variable, gradually decreasing probability
# from x1 to x2:
if x <= x1
  # Return prior probability and derivative:
  return zero(T),zero(T)
elseif x < x2
  # Cubic interpolation:
  y = (x2-x)/(x2-x1)
  return log(3y^2-2y^3),6/(x2-x1)*(y^2-y)
end
return -Inf,zero(T)
end

function log_bounds_lower(x::T,x1::Float64,x2::Float64) where {T <: Real}
@assert(x2 > x1)
# Places bounds on variable, gradually decreasing probability
# from x2 to x1:
if x <= x1
  # Return prior probability of zero below x1:
  return -Inf,zero(T)
elseif x < x2
  # Cubic interpolation between x1 and x2:
  y = (x-x1)/(x2-x1)
  return log(3y^2-2y^3),6/(x2-x1)*(y-y^2)
end
# Above x2, return probability of unity:
return zero(T),zero(T)
end


text:
## CPU

# manually vectorized version of trilinear kernel
function _vscattering(p::NTuple{3,A}, xi::NTuple{3,B}, F::Array{Float64,3}) where {A,B}
    # unpack tuples
    dx, dy, dz = @. 1 / (x[2] - x[1], y[2] - y[1], z[2] - z[1])
    px, py, pz = p
    x, y, z = xi

    # indices of lowermost-left corner of the cell
    # containing the particle
    idx_x, idx_y, idx_z = parent_cell(p, dxi)

    # distance from particle to lowermost-left corner of the cell
    dx_particle = px - x[idx_x]
    dy_particle = py - y[idx_y]
    dz_particle = pz - z[idx_z]

    a1 = VectorizationBase.Vec(
        (
            F[idx_x + 1, idx_y, idx_z], # 3
            F[idx_x + 1, idx_y + 1, idx_z], # 4
            F[idx_x, idx_y + 1, idx_z], # 2
            F[idx_x, idx_y, idx_z], # 1
        )...,
    )
    a2 = VectorizationBase.Vec(
        (
            F[idx_x + 1, idx_y, idx_z + 1], # 3
            F[idx_x + 1, idx_y + 1, idx_z + 1], # 4
            F[idx_x, idx_y + 1, idx_z + 1], # 2
            F[idx_x, idx_y, idx_z + 1], # 1
        )...,
    )

    b1 = VectorizationBase.Vec(
        (
            (1 - dx_particle * dx) * (dy_particle * dy) * (1 - dz_particle * dz),     # 3
            (dx_particle * dx) * (dy_particle * dy) * (1 - dz_particle * dz),       # 4
            (dx_particle * dx) * (1 - dy_particle * dy) * (1 - dz_particle * dz),     # 2
            (1 - dx_particle * dx) * (1 - dy_particle * dy) * (1 - dz_particle * dz),   # 1
        )...,
    )
    b2 = VectorizationBase.Vec(
        (
            (1 - dx_particle * dx) * (dy_particle * dy) * (dz_particle * dz),   # 3
            (dx_particle * dx) * (dy_particle * dy) * (dz_particle * dz),     # 4
            (dx_particle * dx) * (1 - dy_particle * dy) * (dz_particle * dz),   # 2
            (1 - dx_particle * dx) * (1 - dy_particle * dy) * (dz_particle * dz), # 1
        )...,
    )

    # Interpolate field F onto particle
    Fp =
        VectorizationBase.vsum(VectorizationBase.vmul(a1, b1)) +
        VectorizationBase.vsum(VectorizationBase.vmul(a2, b2))

    return Fp
end

function _vscattering(p::NTuple{3,A}, xi::NTuple{3,B}, F::Array{Float32,3}) where {A,B}
    # unpack tuples
    dx, dy, dz = @. 1 / (x[2] - x[1], y[2] - y[1], z[2] - z[1])
    px, py, pz = p
    x, y, z = xi

    # indices of lowermost-left corner of the cell
    # containing the particle
    idx_x, idx_y, idx_z = parent_cell(p, dxi)

    # distance from particle to lowermost-left corner of the cell
    dx_particle = px - x[idx_x]
    dy_particle = py - y[idx_y]
    dz_particle = pz - z[idx_z]

    a = VectorizationBase.Vec(
        (
            F[idx_x + 1, idx_y, idx_z], # 3
            F[idx_x + 1, idx_y + 1, idx_z], # 4
            F[idx_x, idx_y + 1, idx_z], # 2
            F[idx_x, idx_y, idx_z], # 1
            F[idx_x + 1, idx_y, idx_z + 1], # 3
            F[idx_x + 1, idx_y + 1, idx_z + 1], # 4
            F[idx_x, idx_y + 1, idx_z + 1], # 2
            F[idx_x, idx_y, idx_z + 1], # 1
        )...,
    )

    b = VectorizationBase.Vec(
        (
            (1 - dx_particle * dx) * (dy_particle * dy) * (1 - dz_particle * dz),     # 3
            (dx_particle * dx) * (dy_particle * dy) * (1 - dz_particle * dz),       # 4
            (dx_particle * dx) * (1 - dy_particle * dy) * (1 - dz_particle * dz),     # 2
            (1 - dx_particle * dx) * (1 - dy_particle * dy) * (1 - dz_particle * dz),   # 1
            (1 - dx_particle * dx) * (dy_particle * dy) * (dz_particle * dz),       # 3
            (dx_particle * dx) * (dy_particle * dy) * (dz_particle * dz),         # 4
            (dx_particle * dx) * (1 - dy_particle * dy) * (dz_particle * dz),       # 2
            (1 - dx_particle * dx) * (1 - dy_particle * dy) * (dz_particle * dz),     # 1
        )...,
    )

    # Interpolate field F onto particle
    Fp = VectorizationBase.vsum(VectorizationBase.vmul(a, b))

    return Fp
end

function _scattering(
    p::NTuple{3,A}, dxi::NTuple{3,A}, xi::NTuple{3,B}, F::Array{C,3}
) where {A,B,C}
    # unpack tuples
    dx, dy, dz = dxi
    px, py, pz = p
    x, y, z = xi

    # check that the particle is inside the grid
    isinside(px, py, pz, x, y, z)

    # indices of lowermost-left corner of the cell
    # containing the particle
    idx_x, idx_y, idx_z = parent_cell(p, dxi)

    # distance from particle to lowermost-left corner of the cell
    tx = (px - x[idx_x]) / dx
    ty = (py - y[idx_y]) / dy
    tz = (pz - z[idx_z]) / dz

    # Interpolate field F onto particle
    Fp = trilinear(
        tx,
        ty,
        tz,
        F[idx_x, idx_y, idx_z],   # v000
        F[idx_x + 1, idx_y, idx_z],   # v100
        F[idx_x, idx_y, idx_z + 1], # v001
        F[idx_x + 1, idx_y, idx_z + 1], # v101
        F[idx_x, idx_y + 1, idx_z],   # v010
        F[idx_x + 1, idx_y + 1, idx_z],   # v110
        F[idx_x, idx_y + 1, idx_z + 1], # v011
        F[idx_x + 1, idx_y + 1, idx_z + 1], # v111
    )
    return Fp
end

function scattering(xi, F::Array{T,3}, particle_coords) where {T}
    # unpack tuples
    dxi = grid_size(xi)
    np = length(particle_coords[1])
    Fp = zeros(T, np)

    Threads.@threads for i in 1:np
        @inbounds Fp[i] = _scattering(
            (particle_coords[1][i], particle_coords[2][i], particle_coords[3][i]),
            dxi,
            xi,
            F,
        )
    end

    return Fp
end

function scattering!(Fp, xi, F::Array{T,3}, particle_coords) where {T}
    # unpack tuples
    x, y, z = xi
    dxi = (x[2] - x[1], y[2] - y[1], z[2] - z[1])
    np = length(particle_coords[1])

    Threads.@threads for i in 1:np
        @inbounds Fp[i] = _scattering(
            (particle_coords[1][i], particle_coords[2][i], particle_coords[3][i]),
            dxi,
            xi,
            F,
        )
    end
end

## CUDA

function _scattering!(
    Fp::CuDeviceVector{T,1},
    p::NTuple{3,CuDeviceVector{T,1}},
    dxi::NTuple{3,T},
    xi::NTuple{3,A},
    F::CuDeviceArray{T,3},
) where {A,T}
    ix = (blockIdx().x - 1) * blockDim().x + threadIdx().x

    # unpack tuples
    dx, dy, dz = dxi
    px, py, pz = p
    x, y, z = xi

    @inbounds if ix ≤ length(px)

        # check that the particle is inside the grid
        isinside(px[ix], py[ix], pz[ix], x, y, z)

        # indices of lowermost-left corner of the cell
        # containing the particle
        idx_x, idx_y, idx_z = parent_cell((px[ix], py[ix], pz[ix]), dxi)

        # # distance from particle to lowermost-left corner of the cell
        tx = (px[ix] - x[idx_x]) / dx
        ty = (py[ix] - y[idx_y]) / dy
        tz = (pz[ix] - z[idx_z]) / dz

        # Interpolate field F onto particle
        Fp[ix] = trilinear(
            tx,
            ty,
            tz,
            F[idx_x, idx_y, idx_z],   # v000
            F[idx_x + 1, idx_y, idx_z],   # v100
            F[idx_x, idx_y, idx_z + 1], # v001
            F[idx_x + 1, idx_y, idx_z + 1], # v101
            F[idx_x, idx_y + 1, idx_z],   # v010
            F[idx_x + 1, idx_y + 1, idx_z],   # v110
            F[idx_x, idx_y + 1, idx_z + 1], # v011
            F[idx_x + 1, idx_y + 1, idx_z + 1], # v111
        )
    end

    return nothing
end

function scattering(
    xi, Fd::CuArray{T,3}, particle_coords::NTuple{3,CuArray}; nt=512
) where {T}
    x, y, z = xi
    dxi = (x[2] - x[1], y[2] - y[1], z[2] - z[1])
    N = length(particle_coords[1])
    Fpd = CuArray{T,1}(undef, N)

    numblocks = ceil(Int, N / nt)
    CUDA.@sync begin
        @cuda threads = nt blocks = numblocks _scattering!(
            Fpd, particle_coords, dxi, xi, Fd
        )
    end

    return Fpd
end

function scattering!(
    Fpd, xi, Fd::CuArray{T,3}, particle_coords::NTuple{3,CuArray}; nt=512
) where {T}
    x, y, z = xi
    dxi = (x[2] - x[1], y[2] - y[1], z[2] - z[1])
    N = length(particle_coords[1])
    numblocks = ceil(Int, N / nt)
    CUDA.@sync begin
        @cuda threads = nt blocks = numblocks _scattering!(
            Fpd, particle_coords, dxi, xi, Fd
        )
    end
end


text:
using ImagineAnalyses, IntervalSets, Unitful, Test
import Unitful:μm, s, Hz

#find_circular
samps = repeat([1:9...; 10:-1:2...], 3)
circ_list = [2; 4]
pad_nsamps = 0
i_s = find_circular(samps, circ_list, pad_nsamps)
@test i_s[1] == findall(x->x==circ_list[1], samps)
@test i_s[2] == findall(x->x==circ_list[2], samps)
@test all(find_circular(samps, circ_list, 1) .== i_s)
@test_throws Exception find_circular(samps, circ_list, 2)
circ_list = [4; 2]
@test_throws Exception find_circular(samps, circ_list, pad_nsamps) #because it doesn't begin decreasing
samps = repeat([10:-1:2...; 1:9...], 3) #begins decreasing
i_s = find_circular(samps, circ_list, pad_nsamps)
@test i_s[1] == findall(x->x==circ_list[1], samps)
@test i_s[2] == findall(x->x==circ_list[2], samps)

samps = repeat([1:9...; 10:-1:2...], 3)
circ_list = [2.5; 4.0] #fractional
i_s = find_circular(samps, circ_list, pad_nsamps)
@test i_s[1][1:2:end] == findall(x->x==3, samps)[1:2:end] #on increasing sweep
@test i_s[1][2:2:end] == findall(x->x==2, samps)[2:2:end] #on decreasing sweep
@test i_s[2] == findall(x->x==circ_list[2], samps) #unaffected

#again with units
samps *= Unitful.μm
circ_list *= Unitful.μm
i_s = find_circular(samps, circ_list, pad_nsamps)
@test i_s[1][1:2:end] == findall(x->x==3μm, samps)[1:2:end] #on increasing sweep
@test i_s[1][2:2:end] == findall(x->x==2μm, samps)[2:2:end] #on decreasing sweep
@test i_s[2] == findall(x->x==circ_list[2], samps) #unaffected

samps = fill(1.0μm, 30)
@test all(mean_cycle(samps, 10; start_idx = 1).==samps[1:10])
samps[1] = 2.0μm
@test all(mean_cycle(samps, 10; start_idx = 2).==samps[2:11])
samps[1:2] .= 2.0μm
@test all(mean_cycle(samps, 8; start_idx = 3).==samps[3:10])

pulse_mod = Float64.(ImagineInterface.gen_pulses(20, [10..12]))
pulse_mon = Float64.(ImagineInterface.gen_pulses(20, [12..14]))
@test mon_delay(pulse_mod, pulse_mon) == 2
@test mon_delay(pulse_mon, pulse_mod) == -2
@test mon_delay(convert(Vector{Float64}, pulse_mod).*Unitful.μm, convert(Vector{Float64}, pulse_mon).*Unitful.μm) == 2

pos_fwd, pos_back = ImagineInterface.gen_bidi_pos(0.0μm, 20.0μm, 1.0s, 20.0Hz)
mod_cyc = ustrip.(vcat(pos_fwd, pos_back))
mod_cyc = round.(Int, mod_cyc) * unit(pos_fwd[1]) #inexact nature of floats can cause make tests nonintuitive
mon_cyc = mod_cyc
targets = [5μm; 10μm; 15μm]
pad_nsamps = 0
pt = pulse_timings(mod_cyc, mon_cyc, targets, pad_nsamps)
for i=1:length(targets)
    @test all(mon_cyc[pt[i]].==targets[i])
end

flash_cyc, cam_cyc, offset = flash_cam_cycs(mod_cyc, mon_cyc, targets, pad_nsamps, 0.05s, 0.1s, 20.0Hz; is_bidi = true)
@test length(findall(flash_cyc)) == sum([length(x) for x in pt])
@test length(findall(cam_cyc)) == sum([length(x) for x in pt]) * 2
for i=1:length(targets)
    for ii in pt[i]
        @test flash_cyc[ii] == true
        @test all(cam_cyc[ii-1:ii].== true) #Exposure lasts two samples, begins the sample previous to the flash
    end
end

#now with a shift
mon_cyc = circshift(mod_cyc, 1)
pt = pulse_timings(mod_cyc, mon_cyc, targets, pad_nsamps)
for i=1:length(targets)
    @test all(mon_cyc[pt[i]].==targets[i])
end
flash_cyc, cam_cyc, offset = flash_cam_cycs(mod_cyc, mon_cyc, targets, pad_nsamps, 0.05s, 0.1s, 20.0Hz; is_bidi = true)
@test offset == 1
flash_cyc = circshift(flash_cyc, offset)
cam_cyc = circshift(cam_cyc, offset)
@test length(findall(flash_cyc)) == sum([length(x) for x in pt])
@test length(findall(cam_cyc)) == sum([length(x) for x in pt]) * 2
for i=1:length(targets)
    for ii in pt[i]
        @test flash_cyc[ii] == true
        @test all(cam_cyc[ii-1:ii].== true) #Exposure lasts two samples, begins the sample previous to the flash
    end
end

#function cam_flash_cycs(flash_ctr_is, flash_time::HasTimeUnits, exp_time::HasTimeUnits, sample_rate::HasInverseTimeUnits)

v = [1.0:1.0:10.0...;9.0:-1.0:2.0...].*Unitful.μm
zs = [2.0; 4.0; 6.0] .* Unitful.μm
sr = 1.0Hz

@test max_exp(v, sr, zs) == 1.0s
eff = framerate_efficiency(v, sr, zs)
@test eff == 1.0
@test framerate_efficiency(v, sr) == 1.0 #version that doesn't consider slice locations (worst case)

zs = [2.0; 5.0] .* Unitful.μm #last slice is limiting
@test max_exp(v, sr, zs) == 1.0s
eff = framerate_efficiency(v, sr, zs)
@test eff == 1.0

v = [1.0;3.0;5.0;6.0:10.0...;9.0:-1.0:2.0...].*Unitful.μm #nonlinear in the beginning of the forward sweep
zs = [3.0; 6.0] .* Unitful.μm #max exp should be 2.0 in the linear case, 1.0 due to nonlinearity: efficiency of 0.5
@test max_exp(v, sr, zs) == 1.0s
eff = framerate_efficiency(v, sr, zs)
@test eff == 0.5


text:
module ADMPS

export Ising, TFIsing, Heisenberg
export Z,magnetisation, energy
export hamiltonian, model_tensor, mag_tensor, energy_tensor
export num_grad, optimisemps

include("cuda_patch.jl")
include("hamiltonianmodels.jl")
include("exampletensors.jl")
include("exampleobs.jl")
include("environment.jl")
include("autodiff.jl")
include("variationalmps.jl")

end


text:
@testset "demodata" begin
    x, A = demodata()

    @plottest begin
        plot(x, A[1, :], label="1")
        for n in 4:3:13
            plot!(x, A[n, :], label="$n")
        end
        xlabel!("Peak position")
        ylabel!("I")
    end "demodata.png"
end


text:
@testset "data.jl" begin
    @testset "dataset(name)" begin
        # nonexistent file throws
        @test_throws ErrorException Extremes.dataset("nonexistant")

        # portpirie loading
        df = Extremes.dataset("portpirie")
        @test size(df, 1) == 65
        @test size(df, 2) == 2
    end

end


text:
abstract type AbstractIBPoissonSolution{TU,TF} end # subtypes require ψ and f field. Change this eventually to getψ and getf methods. Similarly for RHS types.

struct PoissonSolution{TU}
    ψ::TU
end

struct IBPoissonSolution{TU,TF} <: AbstractIBPoissonSolution{TU,TF}
    ψ::TU
    f::TF
end

struct ConstrainedIBPoissonSolution{TU,TF,T} <: AbstractIBPoissonSolution{TU,TF}
    ψ::TU
    f::TF
    ψ₀::Vector{T}
    δΓ_vec::Vector{T}
end

function _scaletophysicalspace!(sol::PoissonSolution, Δx::Real)
    sol.ψ .*= Δx # The computed ψ is approximately equal to the continuous streamfunction divided by ∆x.
end

function _scaletophysicalspace!(sol::IBPoissonSolution, Δx::Real)
    sol.ψ .*= Δx # The computed ψ is approximately equal to the continuous streamfunction divided by ∆x.
    sol.f .*= Δx # Because Δψ + Rf = -w, f also has to be scaled back. The result is f = γ*Δs
end

function _scaletophysicalspace!(sol::ConstrainedIBPoissonSolution, Δx::Real)
    sol.ψ .*= Δx # The computed ψ is approximately equal to the continuous streamfunction divided by ∆x.
    sol.f .*= Δx # Because Δψ + Rf = -w, f also has to be scaled back. The result is f = γ*Δs
    sol.ψ₀ .*= Δx
    sol.δΓ_vec .*= Δx
end


text:
export IRDFT

"""
`IRDFT([domainType=Float64::Type,] dim_in::Tuple, d::Int, [,dims=1])`

`IRDFT(x::AbstractArray, d::Int, [,dims=1])`

Creates a `LinearOperator` which, when multiplied with a complex array `x`, returns the IDFT over the dimension `dims`, exploiting Hermitian symmetry. Like in the function `BASE.irfft`, `d` must satisfy `div(d,2)+1 == size(x,dims)`.

```julia
julia> A = IRDFT(Complex{Float64},(10,),19)
ℱ⁻¹  ℂ^10 -> ℝ^19

julia> A = IRDFT((5,10,8),19,2)
ℱ⁻¹  ℂ^(5, 10, 8) -> ℝ^(5, 19, 8)

```

"""
struct IRDFT{T <:Number,
	     N,
	     D,
	     T1<:AbstractFFTs.Plan,
	     T2<:AbstractFFTs.Plan,
	     T3<:NTuple{N,Any}
	     } <: LinearOperator
	dim_in::NTuple{N,Int}
	dim_out::NTuple{N,Int}
	A::T1
	At::T2
	idx::T3
end

# Constructors
#standard constructor

function IRDFT(x::AbstractArray{Complex{T},N}, d::Int, dims::Int=1) where {T<:Number,N}
	A = plan_irfft(x,d,dims)
	dim_in = size(x)
	dim_out = ()
	idx = ()
	for i = 1:N
		dim_out = i == dims ? (dim_out...,               d) : (dim_out...,dim_in[i]   )
        idx     = i == dims ? (idx...    , 2:ceil(Int,d/2)) : (idx...    ,Colon()     )
	end
	At = plan_rfft(zeros(dim_out),dims)
	IRDFT{T,N,dims,typeof(A),typeof(At),typeof(idx)}(dim_in,dim_out,A,At,idx)
end

IRDFT(T::Type,dim_in::NTuple{N,Int},d::Int,dims::Int=1) where {N} = IRDFT(zeros(T,dim_in),d,dims)
IRDFT(dim_in::NTuple{N,Int},d::Int,dims::Int=1) where {N} = IRDFT(zeros(Complex{Float64},dim_in),d,dims)

# Mappings

function mul!(y::C1,
              L::IRDFT{T,N,D,T1,T2,T3},
              b::C2) where {N,T,D,T1,T2,T3,C1<:AbstractArray{T,N},
				C2<:AbstractArray{Complex{T},N}}
	mul!(y,L.A,b)
end

function mul!(y::C2,
              L::AdjointOperator{IRDFT{T,N,D,T1,T2,T3}},
              b::C1) where {N,T,D,T1,T2,T3,C1<:AbstractArray{T,N},
				 C2<:AbstractArray{Complex{T},N}}

    A = L.A
	mul!(y,A.At,b)
	y ./= size(b,D)
	@views y[A.idx...] .*= 2
	return y

end

# Properties

size(L::IRDFT) = (L.dim_out,L.dim_in)

fun_name(A::IRDFT) = "ℱ⁻¹"

domainType(L::IRDFT{T,N}) where {T,N} = Complex{T}
codomainType(L::IRDFT{T,N}) where {T,N} = T

is_AAc_diagonal(    L::IRDFT) = false #TODO but might be true?
is_invertible(      L::IRDFT) = true
is_full_row_rank(   L::IRDFT) = true


text:
module TestAqua

import Aqua
import UnionArrays
using Test

# Default `Aqua.test_all(UnionArrays)` does not work due to ambiguities
# in upstream packages.
function test_aqua()
    Aqua.test_all(UnionArrays; ambiguities = false)
end

function test_method_ambiguity()
    Aqua.test_ambiguities(UnionArrays)
end

end  # module


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using FinEtoolsDeforLinear

# Pkg.test()
# include("examples/statics/3-d/bend_hex_spectrum_examples.jl"); bend_hex_spectrum_examples.allrun()
# include("examples/statics/3-d/comp_hex_spectrum_examples.jl"); comp_hex_spectrum_examples.allrun()


text:
# # Flock model
# ![](flock.gif)
#
# The flock model illustrates how flocking behavior can emerge when each bird follows three simple rules:
#
# * maintain a minimum distance from other birds to avoid collision
# * fly towards the average position of neighbors
# * fly in the average direction of neighbors

# ## Defining the core structures

# We begin by calling the required packages and defining an agent type representing a bird.

using Agents, LinearAlgebra
using Random # hide

mutable struct Bird <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    speed::Float64
    cohere_factor::Float64
    separation::Float64
    separate_factor::Float64
    match_factor::Float64
    visual_distance::Float64
end

# The fields `id` and `pos` are required for every agent.
# The field `vel` is required for using [`move_agent!`](@ref) in `ContinuousSpace`.
# `speed` defines how far the bird travels in the direction defined by `vel` per `step`.
# `seperation` defines the minimum distance a bird must maintain from its neighbors.
# `visual_distance` refers to the distance a bird can see and defines a radius of neighboring birds.
# The contribution of each rule defined above recieves an importance weight: `cohere_factor`
# is the importance of maintaining the average position of neighbors,
# `match_factor` is the importance of matching the average trajectory of neighboring birds,
# and `separate_factor` is the importance of maining the minimum
# distance from neighboring birds.

# The function `initialize_model` generates birds and returns a model object using default values.
function initialize_model(;
    n_birds = 100,
    speed = 1.0,
    cohere_factor = 0.25,
    separation = 4.0,
    separate_factor = 0.25,
    match_factor = 0.01,
    visual_distance = 5.0,
    dims = (100, 100),
)
    space2d = ContinuousSpace(2; periodic = true, extend = dims)
    model = ABM(Bird, space2d, scheduler = random_activation)
    for _ in 1:n_birds
        vel = Tuple(rand(2) * 2 .- 1)
        add_agent!(
            model,
            vel,
            speed,
            cohere_factor,
            separation,
            separate_factor,
            match_factor,
            visual_distance,
        )
    end
    index!(model)
    return model
end
nothing # hide

# ## Defining the agent_step!
# `agent_step!` is the primary function called for each step and computes velocity
# according to the three rules defined above.
function agent_step!(bird, model)
    ## Obtain the ids of neighbors within the bird's visual distance
    ids = space_neighbors(bird, model, bird.visual_distance)
    ## Compute velocity based on rules defined above
    bird.vel =
        (
            bird.vel .+ cohere(bird, model, ids) .+ separate(bird, model, ids) .+
            match(bird, model, ids)
        ) ./ 2
    bird.vel = bird.vel ./ norm(bird.vel)
    ## Move bird according to new velocity and speed
    move_agent!(bird, model, bird.speed)
end

distance(a1, a2) = sqrt(sum((a1.pos .- a2.pos) .^ 2))

get_heading(a1, a2) = a1.pos .- a2.pos
nothing # hide

# `cohere` computes the average position of neighboring birds, weighted by importance
function cohere(bird, model, ids)
    N = max(length(ids), 1)
    birds = model.agents
    coherence = (0.0, 0.0)
    for id in ids
        coherence = coherence .+ get_heading(birds[id], bird)
    end
    return coherence ./ N .* bird.cohere_factor
end
nothing # hide

# `separate` repels the bird away from neighboring birds
function separate(bird, model, ids)
    seperation_vec = (0.0, 0.0)
    N = max(length(ids), 1)
    birds = model.agents
    for id in ids
        neighbor = birds[id]
        if distance(bird, neighbor) < bird.separation
            seperation_vec = seperation_vec .- get_heading(neighbor, bird)
        end
    end
    return seperation_vec ./ N .* bird.separate_factor
end
nothing # hide

# `match` computes the average trajectory of neighboring birds, weighted by importance
function match(bird, model, ids)
    match_vector = (0.0, 0.0)
    N = max(length(ids), 1)
    birds = model.agents
    for id in ids
        match_vector = match_vector .+ birds[id].vel
    end
    return match_vector ./ N .* bird.match_factor
end
nothing # hide

# ## Running the model
Random.seed!(23182) # hide
n_steps = 500
model = initialize_model()
step!(model, agent_step!, n_steps)

# ## Plotting the birds
# The great thing about [`plotabm`](@ref) is its flexibility. We can incorporate the
# direction of the birds when plotting them, by making the "marker" function `am`
# create a `Shape`: a triangle with same orientation as the bird's velocity.
# It is as simple as defining the following function:
function bird_triangle(b::Bird)
    φ = atan(b.vel[2], b.vel[1])
    xs = [(i ∈ (0, 3) ? 2 : 1) * cos(i * 2π / 3 + φ) for i in 0:3]
    ys = [(i ∈ (0, 3) ? 2 : 1) * sin(i * 2π / 3 + φ) for i in 0:3]
    Shape(xs, ys)
end
nothing # hide

# And here is the animation
using AgentsPlots
gr() # hide
Random.seed!(23182) # hide
cd(@__DIR__) #src
model = initialize_model()
e = model.space.extend
anim = @animate for i in 0:100
    i > 0 && step!(model, agent_step!, 1)
    p1 = plotabm(
        model;
        am = bird_triangle,
        as = 10,
        showaxis = false,
        grid = false,
        xlims = (0, e[1]),
        ylims = (0, e[2]),
    )
    title!(p1, "step $(i)")
end
gif(anim, "flock.gif", fps = 30)


text:
Pkg.add("DecisionTree")

using RDatasets: dataset
using DecisionTree

iris = dataset("datasets", "iris")
features = convert(Array, iris[:, 1:4]);
labels = convert(Array, iris[:, 5]);

# train full-tree classifier
model = build_tree(labels, features)
# prune tree: merge leaves having >= 90% combined purity (default: 100%)
model = prune_tree(model, 0.9)
# pretty print of the tree, to a depth of 5 nodes (optional)
print_tree(model, 5)
# apply learned model
apply_tree(model, [5.9,3.0,5.1,1.9])
# get the probability of each label
apply_tree_proba(model, [5.9,3.0,5.1,1.9], ["setosa", "versicolor", "virginica"])
# run n-fold cross validation for pruned tree,
# using 90% purity threshold pruning, and 3 CV folds
accuracy = nfoldCV_tree(labels, features, 0.9, 3)

# train random forest classifier
# using 2 random features, 10 trees, 0.5 portion of samples per tree (optional), and a maximum tree depth of 6 (optional)
model = build_forest(labels, features, 2, 10, 0.5, 6)
# apply learned model
apply_forest(model, [5.9,3.0,5.1,1.9])
# get the probability of each label
apply_forest_proba(model, [5.9,3.0,5.1,1.9], ["setosa", "versicolor", "virginica"])
# run n-fold cross validation for forests
# using 2 random features, 10 trees, 3 folds and 0.5 of samples per tree (optional)
accuracy = nfoldCV_forest(labels, features, 2, 10, 3, 0.5)

# train adaptive-boosted stumps, using 7 iterations
model, coeffs = build_adaboost_stumps(labels, features, 7);
# apply learned model
apply_adaboost_stumps(model, coeffs, [5.9,3.0,5.1,1.9])
# get the probability of each label
apply_adaboost_stumps_proba(model, coeffs, [5.9,3.0,5.1,1.9], ["setosa", "versicolor", "virginica"])
# run n-fold cross validation for boosted stumps, using 7 iterations and 3 folds
accuracy = nfoldCV_stumps(labels, features, 7, 3)

n, m = 10^3, 5 ;
features = randn(n, m);
weights = rand(-2:2, m);
labels = features * weights;

# train regression tree, using an averaging of 5 samples per leaf (optional)
model = build_tree(labels, features, 5)
# apply learned model
apply_tree(model, [-0.9,3.0,5.1,1.9,0.0])
# run n-fold cross validation, using 3 folds, averaging of 5 samples per leaf (optional)
# returns array of coefficients of determination (R^2)
r2 = nfoldCV_tree(labels, features, 3, 5)

# train regression forest, using 2 random features, 10 trees,
# averaging of 5 samples per leaf (optional), 0.7 of samples per tree (optional)
model = build_forest(labels,features, 2, 10, 5, 0.7)
# apply learned model
apply_forest(model, [-0.9,3.0,5.1,1.9,0.0])
# run n-fold cross validation on regression forest
# using 2 random features, 10 trees, 3 folds, averaging of 5 samples/leaf (optional),
# and 0.7 porition of samples per tree (optional)
# returns array of coefficients of determination (R^2)
r2 = nfoldCV_forest(labels, features, 2, 10, 3, 5, 0.7)

max_stars_repo_path: Chapter 8/ch08-ensemblelearning.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 8/ch08-ensemblelearning.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 8/ch08-ensemblelearning.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
@testset "processors" begin
    @testset "StackFrames" begin
        cb = CircularArrayBuffer{Float32}(2, 3, 4)
        s = StackFrames(2, 3, 2)
        push!(cb, s)
        @test size(cb) == (2, 3, 1)

        s(ones(Float32, 2, 3))
        @test s[:, :, 1] == zeros(2, 3)
        @test s[:, :, 2] == ones(2, 3)

        push!(cb, s)
        @test size(cb) == (2, 3, 2)

        s = StackFrames(2, 3)  # one dimension lower
        s(ones(2))
        s(2 * ones(2))
        s(3 * ones(2))

        push!(cb, s)
        @test cb[:, :, end] == [1 2 3; 1 2 3]
    end
end


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

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

function SCIPcreateConsPseudobooleanWithConss(scip, cons, name, lincons, linconstype, andconss, andcoefs, nandconss, indvar, weight, issoftcons, intvar, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsPseudobooleanWithConss, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_CONS}, SCIP_LINEARCONSTYPE, Ptr{Ptr{SCIP_CONS}}, Ptr{Cdouble}, Cint, Ptr{SCIP_VAR}, Cdouble, UInt32, Ptr{SCIP_VAR}, Cdouble, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, lincons, linconstype, andconss, andcoefs, nandconss, indvar, weight, issoftcons, intvar, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsPseudoboolean(scip, cons, name, linvars, nlinvars, linvals, terms, nterms, ntermvars, termvals, indvar, weight, issoftcons, intvar, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsPseudoboolean, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{Cdouble}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{SCIP_VAR}, Cdouble, UInt32, Ptr{SCIP_VAR}, Cdouble, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, linvars, nlinvars, linvals, terms, nterms, ntermvars, termvals, indvar, weight, issoftcons, intvar, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicPseudoboolean(scip, cons, name, linvars, nlinvars, linvals, terms, nterms, ntermvars, termvals, indvar, weight, issoftcons, intvar, lhs, rhs)
    ccall((:SCIPcreateConsBasicPseudoboolean, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{Cdouble}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{SCIP_VAR}, Cdouble, UInt32, Ptr{SCIP_VAR}, Cdouble, Cdouble), scip, cons, name, linvars, nlinvars, linvals, terms, nterms, ntermvars, termvals, indvar, weight, issoftcons, intvar, lhs, rhs)
end

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

function SCIPaddTermPseudoboolean(scip, cons, vars, nvars, val)
    ccall((:SCIPaddTermPseudoboolean, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_VAR}}, Cint, Cdouble), scip, cons, vars, nvars, val)
end

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

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

function SCIPgetLinearConsTypePseudoboolean(scip, cons)
    ccall((:SCIPgetLinearConsTypePseudoboolean, libscip), SCIP_LINEARCONSTYPE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

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

function SCIPgetLinDatasWithoutAndPseudoboolean(scip, cons, linvars, lincoefs, nlinvars)
    ccall((:SCIPgetLinDatasWithoutAndPseudoboolean, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{Cint}), scip, cons, linvars, lincoefs, nlinvars)
end

function SCIPgetAndDatasPseudoboolean(scip, cons, andconss, andcoefs, nandconss)
    ccall((:SCIPgetAndDatasPseudoboolean, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_CONS}}, Ptr{Cdouble}, Ptr{Cint}), scip, cons, andconss, andcoefs, nandconss)
end

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

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

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

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

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