text:
# Layer that integrates with k8s
# - deploys pipeline stages on k8s
# - monitors the status of the pipeline stages

using Kuber

const projectdir = dirname(@__DIR__)

const ApiImpl = Kuber.ApiImpl
const CoreV1 = ApiImpl.Typedefs.CoreV1

const Pod = CoreV1.Pod
const Service = CoreV1.Service
const PodList = CoreV1.PodList
const WatchEvent = CoreV1.WatchEvent

"""
    k8s_run

Runs either a Pod or a Service on the k8s cluster
"""
k8s_run(ctx::KuberContext, typ::Symbol, entity) = put!(ctx, entity)
k8s_run(ctx::KuberContext, entity::Pod) = k8s_run(ctx, :Pod, entity)
k8s_run(ctx::KuberContext, entity::Service) = k8s_run(ctx, :Service, entity)

"""
    k8s_delete

Deletes a Pod or Service from the k8s cluster
"""
function k8s_delete(ctx::KuberContext, typ::Symbol, name::String)
    @debug("deleting", typ, name)
    try
        delete!(ctx, typ, name)
    catch
        # ignore not found exceptions
    end
end

"""
    k8s_run_pipeline_command(ctx, name, method)

Runs one of the commands that constitute the pipeline for
preparing the search server for deployment. Each command
deploys a pod on the cluster that invokes a certain method
provided by the search server.

Commands can be:
- fetch_sources
- extract_sources
- index_sources
"""
function k8s_run_pipeline_command(ctx, name, method)
    pod = kuber_obj(ctx, """{
        "apiVersion": "v1",
        "kind": "Pod",
        "metadata":{
            "name": "$name",
            "namespace": "default",
            "labels": {
                "name": "$name"
            }
        },
        "spec": {
            "containers": [{
                "name": "demo",
                "image": "demopipeline:v1",
                "args": ["$method"],
                "volumeMounts": [{
                    "name": "project",
                    "mountPath": "/project"
                }]
            }],
            "volumes": [{
                "name": "project",
                "hostPath": {
                    "path": "$projectdir"
                }
            }],
            "restartPolicy": "Never"
        }
    }""")
    k8s_run(ctx, pod)
end

k8s_fetch_sources(ctx) = k8s_run_pipeline_command(ctx, "fetch", "fetch_sources")
k8s_extract_sources(ctx) = k8s_run_pipeline_command(ctx, "extract", "extract_sources")
k8s_index_sources(ctx) = k8s_run_pipeline_command(ctx, "index", "index_sources")

"""
    k8s_run_search_server(ctx)

Deploy both
- the server as a pod
- a service that exposes it over a port
"""
function k8s_run_search_server(ctx)
    service = kuber_obj(ctx, """{
        "apiVersion": "v1",
        "kind": "Service",
        "metadata": {
            "name": "search",
            "namespace": "default",
            "labels": {"name": "search"}
        },
        "spec": {
            "type": "NodePort",
            "ports": [{
                "port": 5555,
                "nodePort": 30005
            }],
            "selector": {"name": "search"}
        }
    }""")
    k8s_run(ctx, service)

    pod = kuber_obj(ctx, """{
        "apiVersion": "v1",
        "kind": "Pod",
        "metadata":{
            "name": "search",
            "namespace": "default",
            "labels": {
                "name": "search"
            }
        },
        "spec": {
            "containers": [{
                "name": "demo",
                "image": "demopipeline:v1",
                "args": ["run_search_server"],
                "volumeMounts": [{
                    "name": "project",
                    "mountPath": "/project"
                }],
                "ports": [{"containerPort": 5555}]
            }],
            "volumes": [{
                "name": "project",
                "hostPath": {
                    "path": "$projectdir"
                }
            }],
            "restartPolicy": "Never"
        }
    }""")
    k8s_run(ctx, pod)
end

"""
    k8s_delete_search_server

Deletes both the search server pod and the service
"""
function k8s_delete_search_server(ctx::KuberContext=KuberContext())
    k8s_delete(ctx, :Pod, "search")
    k8s_delete(ctx, :Service, "search")
end

"""
    advance

Advances the pipeline between stages
"""
function advance(ctx::KuberContext, stage::Symbol)
    try
        if stage == :init
            stage = :fetching
            k8s_delete_search_server(ctx)
            k8s_fetch_sources(ctx)
        elseif stage == :fetching
            stage = :extracting
            k8s_delete(ctx, :Pod, "fetch")
            k8s_extract_sources(ctx)
        elseif stage == :extracting
            stage = :indexing
            k8s_delete(ctx, :Pod, "extract")
            k8s_index_sources(ctx)
        elseif stage == :indexing
            stage = :running
            k8s_delete(ctx, :Pod, "index")
            k8s_run_search_server(ctx)
        end
    catch ex
        @error("error handling event", exception=(ex,catch_backtrace()))
        # TODO: better error handling here
        # for the demo we just log
    end

    return stage
end

function can_advance(pod, stage)
    # We are interested in knowing when a Pod state is modified to "Succeeded".
    # In reality we should also monitor failure states for error handling.
    (pod.status.phase == "Succeeded") &&
    ((pod.metadata.name == "fetch"   && stage == :fetching) ||
     (pod.metadata.name == "extract" && stage == :extracting) ||
     (pod.metadata.name == "index"   && stage == :indexing))
end

"""
    k8s_update_search_server()

Run this periodically to refresh the search server.

Actions:
- fetch inputs
- extract inputs
- stop the search server if it is running
- index
- start the search server back up
"""
function k8s_update_search_server()
    ctx = KuberContext()
    Kuber.set_api_versions!(ctx)
    stage = :init
    println("starting")

    watch(ctx, list, :Pod) do stream
        for event in stream
            if isa(event, PodList)
                # watch is set up
                # the first event of a watch is where it returns the existing list of pods
                # we trigger the first step of the pipeline here
                stage = advance(ctx, stage)
            elseif isa(event, WatchEvent)
                # subsequent events are WatchEvents on individual pods
                # we are interestes only on pod state modification events
                pod = kuber_obj(ctx, event.object)
                if event.type == "MODIFIED"
                    print("$(pod.metadata.name) $(pod.status.phase)                    \r")
                    if can_advance(pod, stage)
                        stage = advance(ctx, stage)
                        if stage == :running
                            # we are done when our application is updated and reaches the running stage
                            close(stream)
                        end
                    end
                end
            end
        end
        println("")
        print("update complete!")
    end
end

max_stars_repo_name: tanmaykm/KubernetesIntegratedJulia

max_issues_repo_name: tanmaykm/KubernetesIntegratedJulia

max_forks_repo_name: tanmaykm/KubernetesIntegratedJulia


text:
cd(@__DIR__)

using Clang.Generators
using ExprTools, MacroTools, JSON3, JuliaFormatter

const CIMGUI_INCLUDE_DIR = joinpath(@__DIR__,"cimgui-pack","cimgui")
const CIMPLOT_INCLUDE_DIR = joinpath(@__DIR__,"cimgui-pack", "cimplot")
const CIMPLOT_H = normpath(CIMPLOT_INCLUDE_DIR, "cimplot.h")

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
pushfirst!(args, "-DCIMGUI_DEFINE_ENUMS_AND_STRUCTS")
pushfirst!(args, "-isystem$CIMGUI_INCLUDE_DIR")
push!(args, "-I$CIMPLOT_INCLUDE_DIR")

# add definitions
@add_def ImVec2
@add_def ImVec4
@add_def ImGuiMouseButton
@add_def ImGuiKeyModFlags
@add_def ImS8
@add_def ImU8
@add_def ImS16
@add_def ImU16
@add_def ImS32
@add_def ImU32
@add_def ImS64
@add_def ImU64
@add_def ImTextureID
@add_def ImGuiCond
@add_def ImGuiDragDropFlags
@add_def ImDrawList
@add_def ImGuiContext

@add_def ImGuiStyleVar
@add_def ImGuiStyleMod
@add_def ImGuiCol
@add_def ImGuiColorMod
@add_def ImGuiID
@add_def ImGuiStoragePair
@add_def ImGuiTextBuffer
@add_def ImGuiStorage

@add_def ImVector_float
@add_def ImVector_ImU32
@add_def ImVector_ImGuiStyleMod
@add_def ImRect
@add_def ImPoolIdx
@add_def ImVector_ImGuiColorMod

const METADATA_DIR = joinpath(@__DIR__, "cimgui-pack","cimplot","generator", "output")

json_enums = read(joinpath(METADATA_DIR, "structs_and_enums.json"), String);
json_defs = read(joinpath(METADATA_DIR, "definitions.json"), String);
json_typedefs = read(joinpath(METADATA_DIR, "typedefs_dict.json"), String);
enums = JSON3.read(json_enums);
types = JSON3.read(json_typedefs);
FUNCTION_METADATA = JSON3.read(json_defs);
const ENUMS = Symbol.(chop.(string.(propertynames(enums.enums))))
const DESPECIALIZE = ["LinkNextPlotLimits"]

# Check if originally from implot_internal.h; this filters most internal functions
function internal_check(x)
    if hasproperty(x, :location)
        startswith(x.location, "implot_internal") && return true
    end
    return false
end

# Find and extract metadata for specific cimplot function
for objfield in propertynames(FUNCTION_METADATA)
    indexes = nothing
    objvec = getproperty(FUNCTION_METADATA, objfield)
    indexes = findall(internal_check, objvec)
    if !isnothing(indexes)
        append!(options["general"]["printer_blacklist"], getproperty.(objvec[indexes],:ov_cimguiname))
    end
end
#
#for (obj, location) in enums.locations
#    if startswith(location, "implot_internal")
#        push!(options["general"]["printer_blacklist"], String(obj))
#    end
#end

imdatatypes = [:Cfloat, :Cdouble, :ImS8, :ImU8, :ImS16, :ImU16, :ImS32, :ImU32, :ImS64, :ImU64]
jldatatypes = [:Float32, :Float64, :Int8, :UInt8, :Int16, :UInt16, :Int32, :UInt32, :Int64, :UInt64]
imtojl_lookup = Dict(zip(imdatatypes, jldatatypes))
const PRIMITIVE_TYPES = [:ImPlotPoint, :ImPlotRange, :ImVec2, :ImVec4, :ImPlotLimits]
ctx = create_context(CIMPLOT_H, args, options)
build!(ctx, BUILDSTAGE_NO_PRINTING)

function split_ccall(body)
    local funsymbol, rettype, argtypes, argnames
    for ex in body.args
        @capture(ex, ccall((funsymbol_, libcimgui), rettype_, (argtypes__,), argnames__)) && break
    end
    return (funsymbol, rettype, argtypes, argnames)
end

function parse_default(T::DataType, str, ptr_type = :notparsed)
    str == "((void*)0)" && return :C_NULL
    T <: Integer && return (startswith(str, "sizeof") ? :(sizeof($ptr_type)) : Meta.parse(str))
    T <: AbstractFloat && return Meta.parse(str)
    T <: Cstring && return Meta.parse(str)
    T <: Bool && return Meta.parse(str)
    T <: Symbol && return Symbol(str)
    return @warn "Not parsing default value of: $str"
end

function revise_arg(def, metadata, i, sym, jltype, ptr_type = :notparsed)
    if jltype ∈ (:Cint, :Clong, :Cshort, :Cushort, :Culong, :Cuchar, :Cchar,
                 :Int8, :UInt8, :Int16, :UInt16, :Int32, :UInt32, :Int64, :UInt64)
        if hasproperty(metadata.defaults, sym)
            val = parse_default(eval(jltype), getproperty(metadata.defaults,sym), ptr_type)
            def[:args][i] = :($( Expr(:kw, :($sym::Integer), val)) )
        else
            def[:args][i] = :($sym::Integer)
        end
        return
    elseif jltype ∈ (:Cfloat, :Cdouble, :Float64, :Float32)
        if hasproperty(metadata.defaults, sym)
            val = parse_default(eval(jltype), getproperty(metadata.defaults,sym))
            def[:args][i] = :($( Expr(:kw, :($sym::Real), val)) )
        else
            def[:args][i] = :($sym::Real)
        end
        return
    elseif jltype ∈ (:Cstring,:Bool)
        # Don't annotate string arguments--we want to be able to pass C_NULL
        if hasproperty(metadata.defaults, sym)
            val = parse_default(eval(jltype), getproperty(metadata.defaults, sym))
            def[:args][i] = :($(Expr(:kw, sym, val)))
        end

        return
    elseif startswith(string(jltype), "Im")
        if hasproperty(metadata.defaults, sym)
            raw_val = getproperty(metadata.defaults,sym)
            if  raw_val == "((void*)0)"
                val = :C_NULL
                def[:args][i] = :($(Expr(:kw, sym, val)) )
                return
            elseif startswith(raw_val, "Im")
                #if jltype in PRIMITIVE_TYPES && endswith(raw_val, r"\(.+\)")
                #    rx = match(r"\(.+\)",raw_val)
                #    tupex = Meta.parse(rx.match)
                #    def[:args][i] = :($( Expr(:kw, :($sym::$jltype), :($jltype($(tupex.args...))) )))
                #else
                #    def[:args][i] =  :($( Expr(:kw, sym, :($(Symbol(raw_val))))) )
                #end
                if jltype in PRIMITIVE_TYPES
                    def[:args][i] = :($(Expr(:kw, :($sym::$jltype), Meta.parse(raw_val))))
                else
                    def[:args][i] = :($(Expr(:kw, sym, Meta.parse(raw_val))))
                end
                return
            elseif jltype in ENUMS
                def[:args][i] = :($(Expr(:kw, :($sym::Union{$(Symbol(string(jltype)*"_")),Integer}), Meta.parse(raw_val))))
                return
            else
                def[:args][i] = :($(Expr(:kw, sym, Meta.parse(raw_val))))
                return
            end
        elseif jltype in ENUMS
            def[:args][i] = :($sym::Union{$(Symbol(string(jltype)*"_")),Integer})
            return
        else
                #println("Raw value: $raw_val not processed for $sym :: $jltype in $(def[:name])")
            def[:args][i] = :($sym::$jltype)
            return
        end
    elseif @capture(jltype, Ptr{ptrtype_})
        if hasproperty(metadata.defaults, sym)
            raw_val = getproperty(metadata.defaults,sym)
            if raw_val !== "((void*)0)" && endswith(raw_val, r"\(.+\)")
                rx = match(r"\(.+\)",raw_val)
                println("matched ptr default: $(rx.match)")
                tupex = Meta.parse(rx.match)
                def[:args][i] = :($(Expr(:kw,
                :($sym::Union{$ptrtype,AbstractArray{$ptrtype}}), :($ptrtype($(tupex.args...))))))
            end
        elseif ptrtype == :Cstring
            def[:args][i] = :($sym::Union{Ptr{Nothing},String,AbstractArray{String}})
        elseif ptrtype == :Cvoid
            return
        else
           def[:args][i] = :($sym::Union{$ptrtype,AbstractArray{$ptrtype}})
        end
        return
    end
        @info "Not processing argument: $sym::$jltype from $(def[:name])"
        return
end

function make_plotmethod(def, metadata)
    def[:name] = Symbol(metadata.funcname)
    (funsymbol, rettype, argtypes, argnames) = split_ccall(def[:body])
    datatype = :notparsed
    for (i, argtype) in enumerate(argtypes)
        sym = argnames[i]
        jltype = argtype ∈ imdatatypes ? imtojl_lookup[argtype] : argtype
        if @capture(jltype, Ptr{ptrtype_}) && ptrtype ∈ imdatatypes
            datatype = ptrtype
            def[:args][i] = :($sym::Union{Ptr{$ptrtype},Ref{$ptrtype},AbstractArray{$ptrtype}})
        else
            revise_arg(def, metadata, i, sym, jltype, datatype)
        end
    end
    def[:body] = Expr(:block,
                      :(ccall(($funsymbol, libcimgui), $rettype, ($(argtypes...),), $(argnames...))))
end

function make_finalizer!(def, metadata)
    def[:name] = :(Base.finalizer)
    (funsymbol, rettype, argtypes, argnames) = split_ccall(def[:body])
    argtype, argname = only(argtypes), only(argnames)
    @capture(argtype, Ptr{ptrtype_})
    def[:args] = [:($argname::Union{$argtype,$ptrtype})]
    new_ccall = :(ccall(($funsymbol, libcimgui), $rettype, ($argtype,), $argname))
    new_body = Expr(:block, :(ptr = pointer_from_objref($argname)), :(GC.@preserve $argname $new_ccall))
    def[:body] = MacroTools.prewalk(rmlines, new_body)
end

function make_constructor!(def, metadata)
    def[:name] = Symbol(metadata.stname)
    (funsymbol, rettype, argtypes, argnames) = split_ccall(def[:body])
    new_ccall = :(ccall(($funsymbol, libcimgui), $rettype, ($(argtypes...),), $(argnames...)))
    def[:body] = Expr(:block, new_ccall)
end

function make_objmethod!(def, metadata)
    def[:name] = Symbol(metadata.funcname)
    (funsymbol, rettype, argtypes, argnames) = split_ccall(def[:body])
    sym, argtype = first(argnames), first(argtypes)

    @capture(argtype, Ptr{ptr_type_})
    def[:args][1] = :($sym::Union{$ptr_type,$argtype})
    for (i, argtype) in enumerate(argtypes)
        i == 1 && continue
        sym = argnames[i]
        jltype = argtype ∈ imdatatypes ? imtojl_lookup[argtype] : argtype
        # Skip pointer types
        if @capture(jltype, Ptr{ptrtype_})
            ptrtype ∉ vcat(PRIMITIVE_TYPES, imdatatypes) && continue
            if ptrtype in (imdatatypes..., :Cstring)
                if ptrtype == :Cstring
                    def[:args][i] = :($sym::Union{Ptr{Nothing},String,AbstractArray{String}})
                else
                    def[:args][i] = :($sym::Union{Ptr{$ptrtype},Ref{$ptrtype},AbstractArray{$ptrtype}})
                end
                continue
            end
        end
        revise_arg(def, metadata, i, sym, jltype) # offset bc we pop off first arg above
     end
end

function make_nonudt(def, metadata)
    def[:name] = Symbol(metadata.funcname)
    (funsymbol, rettype, argtypes, argnames) = split_ccall(def[:body])
    sym = popfirst!(def[:args])
    @capture(first(argtypes), Ptr{ptr_type_})
    argtypes[1] = :(Ref{$ptr_type})
    def[:body] = Expr(:block,
                      :($sym = Ref{$ptr_type}()),
                      :(ccall(($funsymbol, libcimgui), $rettype, ($(argtypes...),), $(argnames...))),
                      ptr_type in PRIMITIVE_TYPES ? :($sym[]) : :($sym))

   for (i, argtype) in enumerate(argtypes)
       i == 1 && continue
       sym = argnames[i]
       jltype = argtype ∈ imdatatypes ? imtojl_lookup[argtype] : argtype
       # Skip pointer types
       if @capture(jltype, Ptr{ptrtype_})
            ptrtype ∉ vcat(PRIMITIVE_TYPES, imdatatypes) && continue
            if ptrtype in (imdatatypes..., :Cstring)
                if ptrtype == :Cstring
                    def[:args][i] = :($sym::Union{Ptr{Nothing},String,AbstractArray{String}})
                else
                    def[:args][i] = :($sym::Union{Ptr{$ptrtype},Ref{$ptrtype},AbstractArray{$ptrtype}})
                end
                continue
            end
       end
       revise_arg(def, metadata, i-1, sym, jltype) # offset bc we pop off first arg above
    end
end

function make_generic(def, metadata)
    def[:name] = Symbol(metadata.funcname)
    (funsymbol, rettype, argtypes, argnames) = split_ccall(def[:body])
    def[:body] = Expr(:block,
                      :(ccall(($funsymbol, libcimgui), $rettype, ($(argtypes...),), $(argnames...))))
    for (i, argtype) in enumerate(argtypes)
        sym = argnames[i]
        jltype = argtype ∈ imdatatypes ? imtojl_lookup[argtype] : argtype

        # Remove all type annotations from marked functions
        if metadata.funcname in DESPECIALIZE
            def[:args][i] = :($sym)
            continue
        end

        # Skip pointer types
        if @capture(jltype, Ptr{ptrtype_})
            ptrtype ∉ vcat(PRIMITIVE_TYPES, imdatatypes) && continue
            if ptrtype in (imdatatypes..., :Cstring)
                if ptrtype == :Cstring
                    def[:args][i] = :($sym::Union{Ptr{Nothing},String,AbstractArray{String}})
                else
                    def[:args][i] = :($sym::Union{Ptr{$ptrtype},Ref{$ptrtype},AbstractArray{$ptrtype}})
                end
                continue
            end
        end
        revise_arg(def, metadata, i, sym, jltype)
   end
end

function revise_function(ex::Expr, all_metadata, options)
    def = ExprTools.splitdef(ex)
    # Skip Expr function names (e.g. :(Base.getproperty))
    # FIXME: sometimes getting functors??!!
    def[:name] isa Symbol || return ex
    fun_name = string(def[:name])
    # Skip functions not in the JSON metadata
    any(startswith.(fun_name,string.(propertynames(all_metadata)))) || return ex
    local metadata
    # Find and extract metadata for specific cimplot function
    for objfield in propertynames(all_metadata)
        objvec = getproperty(all_metadata, objfield)
        idx = findfirst(x -> isequal(x.ov_cimguiname, fun_name), objvec)
        if !isnothing(idx)
            metadata = objvec[idx]
            break
        end
    end
    if !@isdefined(metadata)
        @warn "Could not find function: $fun_name in JSON metadata"
        return ex #throw("Could not find cimgui function in JSON metadata")
    end
    # Check if it's for a type
    if metadata.stname !== ""
        # Skip constructors/destructors for primitive types--we can handle these with Julia
        if metadata.stname ∉ String.(PRIMITIVE_TYPES)
            if hasproperty(metadata, :destructor)
                make_finalizer!(def, metadata)
                return ExprTools.combinedef(def)
            elseif hasproperty(metadata, :constructor)
                # write contructor...
                make_constructor!(def, metadata)
                return ExprTools.combinedef(def)
            end
        end
        # Fall through to object method
        if !(hasproperty(metadata, :destructor) || hasproperty(metadata, :constructor))
            make_objmethod!(def,metadata)
            return ExprTools.combinedef(def)
        else
            return ex
        end
    elseif startswith(metadata.funcname, "Plot") && !startswith(metadata.funcname, "PlotToPixels")
        # Since Plot functions are templated, dispatch on pointer (data input) arguments
        make_plotmethod(def, metadata)
        return ExprTools.combinedef(def)
    elseif hasproperty(metadata, :nonUDT)
        # Pop the pOut argument and insert a Ref creation and unload...
        make_nonudt(def, metadata)
        return ExprTools.combinedef(def)
    else
        make_generic(def, metadata)
        out = ExprTools.combinedef(def)
        return out
    end
    @warn "function $(def[:name]) not parsed"
        return ex
end

function rewrite!(dag::ExprDAG, metadata, options)
    for node in get_nodes(dag)
        expressions = get_exprs(node)
        for (i, expr) in enumerate(expressions)
            if Meta.isexpr(expr, :function)
                    expressions[i] = revise_function(expr, metadata, options)
            end
        end
    end
end

ctx = create_context(CIMPLOT_H, args, options)
build!(ctx, BUILDSTAGE_NO_PRINTING)
rewrite!(ctx.dag, FUNCTION_METADATA, options)
build!(ctx, BUILDSTAGE_PRINTING_ONLY)
#format(normpath(@__DIR__,"..","src"), YASStyle())

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    localanisotropies(Gradients, grid, propname, w)

Extract `LocalAnisotropy` from a reference scenario. The `propname` variable
from the cartesian `grid` object is scanned and local gradients are extracted
from it. These gradients are smoothed within a squared/cubic window sized
2`w` x 2`w` (x 2`w`) in order to return an ellipse/ellipsoid along the
directions with more/less variability (`w` must be an `Int` value and corresponds
to the number of grid cells). A preliminar magnitude is also extracted from it
and can be later rescaled to proper limits using [`rescale_magnitude`](@ref)
function.

## Example

```julia
lpars = localanisotropies(Gradients, grid, :CO2, 5)
```
"""
function localanisotropies(::Type{Gradients}, obj, prop, window)
    # get dimensions
    dims  = domain(obj).dims
    N     = length(dims)
	propv = Tables.getcolumn(Tables.columns(values(obj)), prop)

    # extract gradients (maybe need an extra method to deal with big datasets)
    factor = maximum(propv) - minimum(propv)
    factor = factor < 1000 ? 100 : 1
    img = round.(Int,(factor .* reshape(propv, Size(dims)) )) # temp solution
    g = imgradients(img, KernelFactors.sobel, "replicate")

    quat = Array{Quaternion}(undef,size(img)) # make some better way to store it
    m = Array{Vector}(undef,size(img)) # make some better way to store it

    Threads.@threads for i in CartesianIndices(img)
        tensor = zeros(Float64,N,N)
        for ng in gridneighbors(img, i, window)
            for j in CartesianIndices((1:N,1:N))
                tensor[j] += g[j[1]][ng]*g[j[2]][ng]
            end
        end

        T = eigen(Symmetric(SMatrix{N,N}(tensor)))
        V = T.vectors[:, sortperm(T.values)]'
        if N==3
            eigv = V
			det(V) < 0 && (eigv = Diagonal(SVector{3}([-1,1,1])) * eigv)
        else
            eigv = zeros(Float64,3,3)
            eigv[1:2,1:2] = V
			eigv[1,1] ≉ eigv[2,2] && (eigv[1,:] .*= -1)
            eigv[3,3] = 1.0
            eigv = SMatrix{3,3}(eigv)
        end
        q = dcm_to_quat(eigv)
        λ = sort(T.values, rev=true)
        quat[i] = q
        m[i] = λ/λ[1]
    end

    LocalAnisotropy(vec(quat), reduce(hcat,vec(m)))
end


text:
import PyPlot.spy
function spy(A::Assoc)
   spy(Adj(logical(A)) );
   X = Col(A);
   Y = Row(A);
   xticks(0:(length(X)-1), X,rotation="vertical")
   yticks(0:(length(Y)-1), Y)

end


text:
#
println("Ap) Test of restricted-active-space (RAS) expansion.")
#

name        = "Beryllium 1s^2 2s^2 ^1S_0 ground state"
refConfigs  = [Configuration("[He] 2s^2")]
rasSettings = RasSettings([1], 24, 1.0e-6, false, true, [1,2,3] )
from        = [Shell("2s")]
#
frozen      = [Shell("1s")]
to          = [Shell("2s"), Shell("2p")]
step1       = RasStep(RasStep(), seFrom=from, seTo=deepcopy(to), deFrom=from, deTo=deepcopy(to), frozen=deepcopy(frozen))
#
append!(frozen, [Shell("2s"), Shell("2p")])
append!(to,     [Shell("3s"), Shell("3p"), Shell("3d")])
step2       = RasStep(step1; seTo=deepcopy(to), deTo=deepcopy(to), frozen=deepcopy(frozen))
#
append!(frozen, [Shell("3s"), Shell("3p"), Shell("3d")])
append!(to,     [Shell("4s"), Shell("4p"), Shell("4d"), Shell("4f")])
step3       = RasStep(step2, seTo=deepcopy(to), deTo=deepcopy(to), frozen=deepcopy(frozen))
#
wa          = Representation(name, Nuclear.Model(4.), Radial.Grid(true), refConfigs,
                             RasExpansion(LevelSymmetry(0, Basics.plus), 4, [step1, step2, step3], rasSettings) )
println("wa = $wa")

wb = generate(wa, output=true)


text:
"""
A type representing a conjunctive normal form.
"""
struct CNF
    N::Int
    M::Int
    clauses::Vector{Vector{Int}}
end

function CNF(clauses::Vector{Vector{Int}})
    M = length(clauses)
    N = maximum(maximum(abs.(c)) for c in clauses)
    return CNF(N, M, clauses)
end

"""
    randomcnf(; N=100, k=3, α=0.1, seed=-1, planted = Vector{Vector{Int}}())

Generates a random instance of the k-SAT problem, with `N` variables
and `αN` clauses.

Any configuration in `planted` is guaranteed to be a solution of the problem.
"""
function randomcnf(; N::Int = 100, k::Int = 3, α::Float64 = 0.1, seed::Int=-1,
                    planted = Vector{Vector{Int}}())
    seed > 0 && Random.seed!(seed)
    M = round(Int, N*α)
    clauses = Vector{Vector{Int}}()
    for p in planted
        @assert length(p) == N   "Wrong size for planted configurations ($N != $(lenght(p)) )"
    end
    for a=1:M
        while true
            c = rand(1:N, k)
            length(union(c)) != k && continue
            c = c .* rand([-1,1], k)

            # reject if not satisfies the planted solutions
            sat = Bool[any(i -> i>0, sol[abs.(c)] .* c) for sol in planted]
            !all(sat) && continue

            push!(clauses, c)
            break
        end
    end
    return CNF(N, M, clauses)
end

"""
    readcnf(fname::String)

Reads a CNF from file `fname`.
"""
function readcnf(fname::String)
    f = open(fname, "r")
    head = split(readline(f))
    N, M = parse(Int, head[3]), parse(Int, head[4])
    clauses = Vector{Vector{Int}}()
    for i=1:M
        line = readline(f)
        c = [parse(Int64, e) for e in split(line)[1:end-1]]
        push!(clauses, c)
    end
    return CNF(N, M, clauses)
end

"""
    writecnf(fname::String, cnf::CNF)

Writes `cnf` to file `fname`.
"""
function writecnf(fname::String, cnf::CNF)
    f = open(fname, "w")
    println(f, "p cnf $(cnf.N) $(cnf.M)")
    for c in cnf.clauses
        for i in c
            print(f, i, " ")
        end
        print(f, "0\n")
    end
    close(f)
end

"""
    energy(cnf, σ)

Counts the number of violated clauses.
"""
function energy(cnf::CNF, σ::Vector{Int})
    E = 0
    for c in cnf.clauses
        issatisfied = false
        for i in c
            if sign(i) == σ[abs(i)]
                issatisfied = true
                break
            end
        end
        E += issatisfied ? 0 : 1
    end
    E
end

"""
    adjlist(cnf::CNF)

Returns a vector containing for each variable a vector of the adjacent clause indexes.
"""
function adjlist(cnf::CNF)
    adj = Vector{Vector{Int}}(cnf.N)
    for (μ, c) in enumerate(cnf.clauses)
        for i in c
            push!(adj[abs(i)], μ)
        end
    end
end


text:
## Unsafe scalar indexing

# Extract relevant fields of a DataArray to a tuple
# The extracted tuple can be passed to `unsafe_ismissing`,
# `unsafe_getindex_notmissing`, `unsafe_setmissing!`, `unsafe_setnotmissing!`, and
# `unsafe_dasetindex!`. This has a meaningful performance impact within
# very tight loops.
daextract(da::DataArray) = (da.data, da.na.chunks)
daextract(pda::PooledDataArray) = (pda.refs, pda.pool)
daextract(a) = nothing

# Check for missing
unsafe_ismissing(da::DataArray, extr, idx::Real) = Base.unsafe_bitgetindex(extr[2], idx)
unsafe_ismissing(pda::PooledDataArray, extr, idx::Real) = extr[1][idx] == 0
unsafe_ismissing(a, extr, idx::Real) = false
unsafe_getindex_notmissing(da::DataArray, extr, idx::Real) = getindex(extr[1], idx)
unsafe_getindex_notmissing(pda::PooledDataArray, extr, idx::Real) = getindex(extr[2], extr[1][idx])
unsafe_getindex_notmissing(a, extr, idx::Real) = Base.unsafe_getindex(a, idx)

# Set missing or data portion of DataArray

unsafe_bitsettrue!(chunks::Vector{UInt64}, idx::Real) =
    chunks[Base._div64(Int(idx)-1)+1] |= (UInt64(1) << Base._mod64(Int(idx)-1))
unsafe_bitsetfalse!(chunks::Vector{UInt64}, idx::Real) =
    chunks[Base._div64(Int(idx)-1)+1] &= ~(UInt64(1) << Base._mod64(Int(idx)-1))

unsafe_setmissing!(da::DataArray, extr, idx::Real) = unsafe_bitsettrue!(extr[2], idx)
unsafe_setmissing!(da::PooledDataArray, extr, idx::Real) = setindex!(extr[1], 0, idx)
unsafe_setnotmissing!(da::DataArray, extr, idx::Real) = unsafe_bitsetfalse!(extr[2], idx)
unsafe_setnotmissing!(da::PooledDataArray, extr, idx::Real) = nothing

# Fast setting of missing values in DataArrays
# These take the data and chunks (extracted as da.data and
# da.na.chunks), a value, and a linear index. They assume
# a certain initialization pattern:
#
# - For DataArrays, da.na should be falses
# - For PooledDataArrays, pda.refs should be zeros
unsafe_dasetindex!(data::Array, na_chunks::Vector{UInt64}, val::Missing, idx::Real) =
    unsafe_bitsettrue!(na_chunks, idx)
unsafe_dasetindex!(data::Array, na_chunks::Vector{UInt64}, val, idx::Real) =
    setindex!(data, val, idx)
unsafe_dasetindex!(da::DataArray, extr, val::Missing, idx::Real) =
    unsafe_setmissing!(da, extr, idx)
unsafe_dasetindex!(da::PooledDataArray, extr, val::Missing, idx::Real) = nothing
unsafe_dasetindex!(da::DataArray, extr, val, idx::Real) = setindex!(extr[1], val, idx)
unsafe_dasetindex!(pda::PooledDataArray, extr, val, idx::Real) =
    setindex!(extr[1], getpoolidx(pda, val), idx)
unsafe_dasetindex!(a::AbstractArray, extr, val, idx::Real) = setindex!(a, val, idx)

## PooledDataArray helper functions

# Append newpool to pool. Return indices of newpool in pool.
function combine_pools!(pool, newpool)
    seen = Dict{eltype(pool),Int}()
    sizehint!(seen, length(pool)+length(newpool))

    # Create mapping from pool elements to indices
    i = 0
    for elem in pool
        seen[elem] = (i += 1)
    end

    # Find pool elements in existing array, or add them
    poolidx = Vector{Int}(uninitialized, length(newpool))
    for j = 1:length(newpool)
        poolidx[j] = Base.@get!(seen, newpool[j], (push!(pool, newpool[j]); i += 1))
    end
    poolidx
end

## General indexing functions

# Indexing with missing throws an error
function Base.to_index(A::DataArray)
    any(A.na) && throw(MissingException("cannot convert missing value to index"))
    Base.to_index(A.data)
end

if isdefined(Base, :checkindex) && isdefined(Base, :AbstractUnitRange)
    Base.checkindex(::Type{Bool}, ::AbstractUnitRange, ::Missing) =
        throw(MissingException("missing values are not allowed in indices"))
elseif isdefined(Base, :checkindex)
    Base.checkindex(::Type{Bool}, ::UnitRange, ::Missing) =
        throw(MissingException("missing values are not allowed in indices"))
else
    Base.checkbounds(::Type{Bool}, sz::Int, I::AbstractDataVector{Bool}) = length(I) == sz
    function Base.checkbounds(::Type{Bool}, sz::Int, I::AbstractDataArray{T}) where {T<:Real}
        any(ismissing, I) && throw(MissingException("missing values are not allowed in indices"))
        extr = daextract(I)
        b = true
        for i = 1:length(I)
            @inbounds v = unsafe_getindex_notmissing(I, extr, i)
            b &= Base.checkbounds(Bool, sz, v)
        end
        b
    end
end

import Base: index_shape, index_lengths, setindex_shape_check

if isdefined(Base, :OneTo)
    _index_shape(x...) = Base.to_shape(index_shape(x...))
else
    _index_shape = index_shape
end

# Fallbacks to avoid ambiguity
Base.setindex!(t::AbstractDataArray, x, i::Real) =
    throw(MethodError(setindex!, typeof(t), typeof(x), typeof(i)))
Base.getindex(t::AbstractDataArray, i::Real) =
    throw(MethodError(getindex, typeof(t), typeof(i)))

## getindex: DataArray

Base.IndexStyle(::Type{<:AbstractDataArray}) = Base.IndexLinear()

# Scalar case
function Base.getindex(da::DataArray, I::Real)
    if getindex(da.na, I)
        return missing
    else
        return getindex(da.data, I)
    end
end

@generated function Base._unsafe_getindex!(dest::DataArray, src::DataArray, I::Union{Real, AbstractArray}...)
    N = length(I)
    quote
        $(Expr(:meta, :inline))
        flipbits!(dest.na) # similar initializes with missings
        @nexprs $N d->(J_d = I[d])
        srcextr = daextract(src)
        destextr = daextract(dest)
        srcsz = size(src)
        D = eachindex(dest)
        Ds = start(D)
        @nloops $N j d->J_d begin
            offset_0 = @ncall $N sub2ind srcsz j
            d, Ds = next(D, Ds)
            if unsafe_ismissing(src, srcextr, offset_0)
                unsafe_dasetindex!(dest, destextr, missing, d)
            else
                unsafe_dasetindex!(dest, destextr, unsafe_getindex_notmissing(src, srcextr, offset_0), d)
            end
        end
        dest
    end
end

## getindex: PooledDataArray

# Scalar case
function Base.getindex(pda::PooledDataArray, I::Real)
    if getindex(pda.refs, I) == 0
        return missing
    else
        return pda.pool[getindex(pda.refs, I)]
    end
end

@inline function Base.getindex(pda::PooledDataArray, I::Integer...)
    if getindex(pda.refs, I...) == 0
        return missing
    else
        return pda.pool[getindex(pda.refs, I...)]
    end
end

# Vector case
@inline function Base.getindex(A::PooledDataArray, I::Union{AbstractVector,Colon}...)
    PooledDataArray(RefArray(getindex(A.refs, I...)), copy(A.pool))
end

## setindex!: DataArray

function Base.setindex!(da::DataArray, val::Missing, i::Real)
    da.na[i] = true
    return da
end

function Base.setindex!(da::DataArray, val, ind::Real)
    da.data[ind] = val
    da.na[ind] = false
    return da
end

## setindex!: PooledDataArray

function Base.setindex!(pda::PooledDataArray, val::Missing, ind::Real)
    pda.refs[ind] = 0
    return pda
end

function Base.setindex!(x::PooledDataArray, val, ind::Real)
    x.refs[ind] = getpoolidx(x, val)
    return x
end

## setindex!: both DataArray and PooledDataArray

@generated function Base.setindex!(A::AbstractDataArray, x, J::Union{Real,Colon,AbstractArray}...)
    N = length(J)
    quote
        if !isa(x, AbstractArray) && isa(A, PooledDataArray)
            # Only perform one pool lookup when assigning a scalar value in
            # a PooledDataArray
            setindex!(A.refs, getpoolidx(A, x), J...)
            return A
        end

        Aextr = daextract(A)
        @nexprs $N d->(I_d = Base.to_indices(A, J)[d])
        @ncall $N checkbounds A I
        stride_1 = 1
        @nexprs $N d->(stride_{d+1} = stride_d*size(A,d))
        @nexprs $N d->(offset_d = 1)  # really only need offset_$N = 1
        if !isa(x, AbstractArray)
            @nloops $N i d->I_d d->(@inbounds offset_{d-1} = offset_d + (i_d - 1)*stride_d) begin
                if isa(x, Missing)
                    @inbounds unsafe_setmissing!(A, Aextr, offset_0)
                else
                    @inbounds unsafe_setnotmissing!(A, Aextr, offset_0)
                    @inbounds unsafe_dasetindex!(A, Aextr, x, offset_0)
                end
            end
        else
            X = x
            idxlens = @ncall $N index_lengths I
            @ncall $N setindex_shape_check X (d->idxlens[d])
            k = 1
            if isa(A, PooledDataArray) && isa(X, PooledDataArray)
                # When putting one PDA into another, first unify the pools
                # and then translate the references
                poolmap = combine_pools!(A.pool, X.pool)
                Arefs = A.refs
                Xrefs = X.refs
                @nloops $N i d->I_d d->(@inbounds offset_{d-1} = offset_d + (i_d - 1)*stride_d) begin
                    @inbounds Arefs[offset_0] = Xrefs[k] == 0 ? 0 : poolmap[Xrefs[k]]
                    k += 1
                end
            else
                Xextr = daextract(X)
                @nloops $N i d->I_d d->(@inbounds offset_{d-1} = offset_d + (i_d - 1)*stride_d) begin
                    @inbounds if isa(X, AbstractDataArray) && unsafe_ismissing(X, Xextr, k)
                        unsafe_setmissing!(A, Aextr, offset_0)
                    else
                        unsafe_setnotmissing!(A, Aextr, offset_0)
                        unsafe_dasetindex!(A, Aextr, unsafe_getindex_notmissing(X, Xextr, k), offset_0)
                    end
                    k += 1
                end
            end
        end
        A
    end
end


text:
using LinearAlgebra
using Test

using CompScienceMeshes
using BEAST

import CompScienceMeshes.SComplex2D

fn = joinpath(@__DIR__, "assets", "thick_cladding.msh")
G01 = SComplex2D(CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma01"))
G02 = SComplex2D(CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma02"))
G12 = SComplex2D(CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma12"))

G = SComplex2D(CompScienceMeshes.read_gmsh_mesh(fn))

b = boundary(G02)
G = weld(G01,G02,seam=b)
@test numcells(G01) + numcells(G02) == numcells(G)
@test numcells(skeleton(G01,1)) + numcells(skeleton(G02,1)) == numcells(skeleton(G,1)) + numcells(b)

G = weld(G,G12,seam=b)
@test numcells(G01) + numcells(G02) + numcells(G12) == numcells(G)
@test numcells(skeleton(G01,1)) + numcells(skeleton(G02,1)) + numcells(skeleton(G12,1)) == numcells(skeleton(G,1)) + 2*numcells(b)

G1 = weld(-G01, -G12, seam=b)
@test CompScienceMeshes.isoriented(G1)
G2 = weld(-G02, G12, seam=b)
@test CompScienceMeshes.isoriented(G2)

isclosed(m) = (numcells(boundary(m)) == 0)

@test isclosed(G1)
@test isclosed(G2)

E = skeleton(G,1)
E1 = skeleton(G1,1)
A1 = CompScienceMeshes.embedding(E1,E)

nd = BEAST.nedelec(G,E)
@test numfunctions(nd) == numcells(E)
nd1 = BEAST.nedelec(G1,E1)
@test numfunctions(nd1) == numcells(E1)
Z_direct = assemble(BEAST.Identity(), nd1, nd)
Z_expand = assemble(BEAST.Identity(), nd1, nd1) * A1
@test norm(Z_direct - Z_expand, Inf) ≤ √(eps(1.0))


text:
export BIFM, BIFMMeta, functional_dependencies

@doc raw"""
The BIFM node is a node that can be used as a substitution for a state space model. It includes all factor of the time slice to perform efficient inference.
This node needs to be used in conjuction with the BIFMHelper node for efficient inference.

```julia
out ~ BIFM(in, zprev, znext)
```

Interfaces:
1. out - latent output (observation) of the BIFM node
2. in - latent input of the BIFM node
3. zprev - previous latent state of the BIFM node
4. znext - next latent state of the BIFM node

*Note: When performing inference, first subscribe to the marginals (in the order: z, out, in) and then to the free energy score function.*

## Example

```julia
# set priors
z_prior ~ MvNormalMeanPrecision(zeros(latent_dim), diagm(ones(latent_dim)))
z_tmp   ~ BIFMHelper(z_prior) where { q = MeanField()}

# update last/previous hidden state
z_prev = z_tmp

# loop through observations
for i in 1:nr_samples

    # specify input as random variable
    u[i]   ~ MvNormalMeanPrecision(μu, Wu)

    # specify observation
    xt[i]  ~ BIFM(u[i], z_prev, z[i]) where { meta = BIFMMeta(A, B, C) }
    x[i]   ~ MvNormalMeanPrecision(xt[i], Wx)

    # update last/previous hidden state
    z_prev = z[i]

end
```
"""
struct BIFM end

@node BIFM Deterministic [ out, in, zprev, znext ]

@doc raw"""
The BIFMMeta structure contains all characterizing information of the BIFM node. In addition, it stores intermediate variables for efficient computations.
"""
mutable struct BIFMMeta{T}
    "Transition matrix A: converting the previous to the next latent state."
    A       :: Matrix{T}
    "Transition matrix B: converting the input to the latent state."
    B       :: Matrix{T}
    "Transition matrix C: converting the latent state to the output."
    C       :: Matrix{T}
    H       :: Union{Matrix{T}, Nothing}
    BHBt    :: Union{Matrix{T}, Nothing}
    ξz      :: Union{Vector{T}, Nothing}
    Λz      :: Union{Matrix{T}, Nothing}
    ξztilde :: Union{Vector{T}, Nothing}
    Λztilde :: Union{Matrix{T}, Nothing}
    μu      :: Union{Vector{T}, Nothing}
    Σu      :: Union{Matrix{T}, Nothing}
end

@doc raw"""
Initialization of the BIFMMeta object can be performed by calling
```julia
meta = BIFMMeta(A, B, C)
```
where `A`, `B` and `C` are the transition matrices in the model.
"""
function BIFMMeta(A::Array{T1, 2}, B::Array{T2, 2}, C::Array{T3, 2}) where { T1, T2, T3 }
    T = promote_type(T1, T2, T3)
    # check whether the dimensionality of transition matrices makes sense
    @assert size(A,1) == size(B,1)
    @assert size(A,1) == size(C,2)

    # return default metadata for BIFM node
    return BIFMMeta{T}(A, B, C, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing)
end
@doc raw"""
Initialization of the BIFMMeta object can be performed by calling
```julia
meta = BIFMMeta(A, B, C, μu, Σu)
```
where `A`, `B` and `C` are the transition matrices in the model. `μu` and `Σu` are the mean vector and covariance matrix of the input. Importantly, in this setting we assume that these are known and do not change due to an external model.
"""
function BIFMMeta(A::Array{T1, 2}, B::Array{T2, 2}, C::Array{T3, 2}, μu::Array{T4, 1}, Σu::Array{T5, 2}) where { T1, T2, T3, T4, T5 }
    T = promote_type(T1, T2, T3, T4, T5)
    # check whether the dimensionality of transition matrices makes sense
    @assert size(A,1)  == size(B,1)
    @assert size(A,1)  == size(C,2)
    @assert size(μu,1) == size(Σu,1)
    @assert size(μu,1) == size(B,2)

    # return default Meta data for BIFM node
    return BIFMMeta{T}(A, B, C, nothing, nothing, nothing, nothing, nothing, nothing, μu, Σu)
end

getA(meta::BIFMMeta)                = meta.A
getB(meta::BIFMMeta)                = meta.B
getC(meta::BIFMMeta)                = meta.C
getH(meta::BIFMMeta)                = meta.H
getBHBt(meta::BIFMMeta)             = meta.BHBt
getξz(meta::BIFMMeta)               = meta.ξz
getΛz(meta::BIFMMeta)               = meta.Λz
getξztilde(meta::BIFMMeta)          = meta.ξztilde
getΛztilde(meta::BIFMMeta)          = meta.Λztilde
getμu(meta::BIFMMeta)               = meta.μu
getΣu(meta::BIFMMeta)               = meta.Σu

function setH!(meta::BIFMMeta, H)
    meta.H = H
end

function setBHBt!(meta::BIFMMeta, BHBt)
    meta.BHBt = BHBt
end

function setξztilde!(meta::BIFMMeta, ξztilde)
    meta.ξztilde = ξztilde
end

function setΛztilde!(meta::BIFMMeta, Λztilde)
    meta.Λztilde = Λztilde
end

function setξz!(meta::BIFMMeta, ξz)
    meta.ξz = ξz
end

function setΛz!(meta::BIFMMeta, Λz)
    meta.Λz = Λz
end

function setμu!(meta::BIFMMeta, μu)
    meta.μu = μu
end

function setΣu!(meta::BIFMMeta, Σu)
    meta.Σu = Σu
end

default_meta(::Type{ BIFM }) = error("BIFM node requires meta flag explicitly specified")


text:
# capturing MSGibbs' final label selection per sample

##

using ApproxManifoldProducts
# using KernelDensityEstimate

##

# build _BiDictMap between variable products

##

X1 = kde!([1;2;3.0],[1.0;]);
X2 = kde!([0.5;1.5;2.5],[1.0;]);
X3 = kde!([4;5;6.0],[1.0;]);

##

dpe = AMP._buildDensityProductElements([X1;X2;X3], inNames=[:X1, :X2, :X3], outName = :Y)

##

AMP._recalcProductKernel(dpe, 1)

##

AMP._recalcProductKernel(dpe, 2)

AMP._recalcProductKernel(dpe, 3)

##

AMP._listOutElementLabelSelections(dpe, 1)
AMP._listOutElementLabelSelections(dpe, 2)
AMP._listOutElementLabelSelections(dpe, 3)

##

using LightGraphs
using TensorCast

## new basic example

X_ = kde!(5 .+ randn(5), [1.0;]);
# Y_ = kde!(randn(5), [1.0;]);
Z_ = kde!(5 .+ randn(5), [1.0;]);

##

# Y_ = X_*Z_

dpeY = AMP._buildDensityProductElements([X_; Z_], outName = :Y, inNames=[:X, :Z], inFctNames=[:XYf1; :YZf1])

@cast pts[i,j] := dpeY.outElements[j][i]
Y_ = kde!(pts, dpeY.outBW[1])

##

X = kde!(pts .- 5)
dpeX = AMP._buildDensityProductElements([X;], outName = :X, inNames=[:Y;], inFctNames=[:XYf1;])

##

Z = kde!(pts .+ 5)
dpeZ = AMP._buildDensityProductElements([Z;], outName = :Z, inNames=[:Y;], inFctNames=[:YZf1;])

##

struct _AMPDiGraph{G,B}
  dg::G
  bd::B
end

function buildDiGraphKernelProduct!(dpel::AMP.DensityProductElements,
                                    abg::_AMPDiGraph = _AMPDiGraph(DiGraph(), AMP._BiDictMap(sizehint=100))  )
  #

  dg = abg.dg
  _bd = abg.bd

  vari = dpel.inElemNames
  N = length(dpel.outElements)

  for s in vari
    # one vertex to represent the variable
    if !haskey(_bd, s)
      add_vertex!(dg)
      _bd[s] = nv(dg)
    end
    for i in 1:N
      # one vertex for each kernel in the vertex belief
      kers = Symbol(s,"_$i")
      if !haskey(_bd, kers)
        add_vertex!(dg)
        _bd[kers] = nv(dg)
        add_edge!(dg, _bd[kers], _bd[s])
      end
    end
  end

  ## add who is getting the product

  s = dpel.outElemName
  if !haskey(_bd, s)
    add_vertex!(dg)
    _bd[s] = nv(dg)

    for i in 1:length(dpel.outElements)
      # one vertex for each kernel in the vertex belief
      add_vertex!(dg)
      kers = Symbol(s,"_$i")
      _bd[kers] = nv(dg)
      add_edge!(dg, _bd[kers], _bd[s])
    end
  end

  ## go through product and add edges from product of proposals

  # dId incoming density ID
  # sId selected label ID
  # eId product element ID
  # kCmb vector of label combinations from incoming densities
  for (eId, kCmb) in enumerate(dpel.lblCombinations), (dId,sId) in enumerate(kCmb)
    # @show dId, eId, sId, kCmb
    # @show dpel.inElemNames
    from = Symbol(dpel.inElemNames[dId], "_$sId")
    to = Symbol(dpel.outElemName, "_$eId")
    # add if the edge does not yet exist
    if !(_bd[to] in outneighbors(dg, _bd[from]) )
      add_edge!(dg, _bd[from], _bd[to])
    end
  end

  return abg
end

# from dependent to independent variable, `d/d ind  dep = dz/dx`
# used for TAF gradient calculations
function findPathsOfKernels(abg::_AMPDiGraph, dep::Symbol, ind::Symbol)
  #
  allpths = Vector{Vector{Symbol}}()

  # get number of independent variables
  allsyms = string.([abg.bd[n] for n in keys(abg.bd)])
  @show Ni = sum(occursin.("$(ind)_", allsyms) )
  @show Nj = sum(occursin.("$(dep)_", allsyms) )

  for i in 1:Ni, j in 1:Nj
    # get the path
    spth = dijkstra_shortest_paths(abg.dg, abg.bd[Symbol(dep,"_$i")])
    # trim to first hit on independent variable
    pth1 = enumerate_paths(spth, abg.bd[Symbol(ind, "_$j")]) .|> x->getindex(abg.bd, x)
    strp = string.(pth1)
    len = occursin.("$ind", strp)
    len_ = 0 < length(len) ? findfirst(len) : 0
    # trim start for first hit on dependent variable
    # pth1 = enumerate_paths(spth, abg.bd[Symbol(dep, "_$j")]) .|> x->getindex(abg.bd, x)
    stt = occursin.("$dep", strp)
    stt_ = 0 < length(stt) ? findlast(stt) : 0

    # assemble all paths
    stt_ < len_ ? push!(allpths, pth1[stt_:len_]) : nothing
  end

  # return all unique paths
  unique(allpths)
end

##

abg = buildDiGraphKernelProduct!( dpeY )

##

abg = buildDiGraphKernelProduct!( dpeX, abg )
abg = buildDiGraphKernelProduct!( dpeZ, abg )

##

# lets see what happened
using GraphPlot, Compose, Cairo

##

pth = "/tmp/gplot.pdf"
gplot(abg.dg, nodelabel=(i->abg.bd[i]).(keys(abg.bd)) ) |> PDF(pth)
@async run(`evince $pth`)

##

using Gadfly
Gadfly.set_default_plot_size(25cm,20cm)

##

spy(adjacency_matrix(abg.dg))

## list all paths worth find gradients for

unqpths = findPathsOfKernels(abg, :X, :Z)

unqpths = findPathsOfKernels(abg, :Z, :X)

## Get value and factor pairings

function _getKernelVariablesOnPath(unqpath::AbstractVector{Symbol})
  varnames = Vector{Symbol}()
  for ker in unqpath
    chplen = length(split(string(ker), '_')[end]) + 1
    nnm = Symbol(string(ker)[1:end-chplen])
    push!(varnames, nnm)
  end
  varnames
end

function _getKernelIdInPath(unqpath::AbstractVector{Symbol}, idx::Int)
  m = match(r"_[0-9]+", unqpath[idx] |> string)
  parse(Int, split(m.match, '_')[end])
end

# get values and measurements associated with each selection

function _rebuildVariableKernelPairs(uqph::AbstractVector{Symbol})
  allpairs = Vector{Tuple{Symbol, Int}}()
  pthvars = _getKernelVariablesOnPath(uqph)
  for i in 1:length(pthvars)
    sId = _getKernelIdInPath(uqph,i)
    push!(allpairs, (pthvars[i], sId) )
  end

  return allpairs
end

##

_rebuildVariableKernelPairs(unqpths[3])

# when proposal is done in IIF, must also store the associated factor connected variables that were use for n-ary case

##

xval = dpeX.outElements[sId]

sId = _getKernelIdInPath(uqph,3)
zval = dpeZ.outElements[sId]

## get measurements for each factor along path

using IncrementalInference

##

varList = [:X, :Y, :Z]

whatvar = unqpths[1][1]

onei = outneighbors(abg.dg, abg.bd[whatvar]) .|> x->getindex(abg.bd, x)
intersect(varList, onei)

## get the measurement for each factor

## each path is a minimization problem that must be solved to get gradients

##

findFactorsBetweenFrom()

##

## ===============================================================================================

##

dpe.inElements[3][1]

# idea is not to `sequence` a change on own dpe, but rather inform other clique variable of a change
# put on the other variable changes to those proposals
_sequenceUpdateProposalElement(dpe, [4.0;], 3, 1)
# _updateProposalElement(dpe, 4.0, 3, 1)

##

# on the other variable now take the updated proposal element (incoming) and
# find which outgoing Element needs to be updated
_takeOutElementFromSequence(dpe)

## test case

## the subgraph
# X--Y--Z

## find Y
# localProduct!(Y)
# this should list dpeY. with inProductLabel containing random X, Z labels
# update source factor label selections

## do single factor links last
# localProduct!(X)
# update source factor label selections

# localProduct!(Z)
# update source factor label selections

#

#

##

# gradient from kernel perturbation

## DEV CODE BELOW

glbs = KDE.makeEmptyGbGlb();
glbs.recordChoosen = true

p123 = *([X1;X2;X3], glbs=glbs, addEntropy=false );

## Analyze what happened

lc = glbs.labelsChoosen

Ndens = 3
XX = [X1;X2;X3]

for i in 1:3
  μ = zeros(Ndens)
  for densid in 1:Ndens
    level = length(lc[i][densid])
    μ[densid] = getPoints(XX[densid])[lc[i][densid][level] - 3]
  end

  Λμ = sum(μ)

  @show μ = 1/3*Λμ
end

##

#


text:
import NLPModels: obj, grad, grad!, objgrad, objgrad!, objcons, objcons!, cons, cons!,
                  jth_con, jth_congrad, jth_congrad!, jth_sparse_congrad,
                  jac_structure!, jac_structure, jac_coord!, jac_coord,
                  jac, jprod, jprod!, jtprod, jtprod!, jac_op, jac_op!,
                  jth_hprod, jth_hprod!, ghjvprod, ghjvprod!,
                  hess_structure!, hess_structure, hess_coord!, hess_coord,
                  hess, hess_op, hess_op!, hprod, hprod!

"""
Convert an MPCCModel to a parametric NLPModels as follows.
Definit le type RlxMPCC :
min 	f(x)
s.t. 	l <= x <= u
	 lcon <= cnl(x,t,tb) <= ucon

with

cnl(x) := c(x),G(x)+tb,H(x)+tb,Phi(G(x).*H(x),t)
"""
mutable struct RlxMPCC <: AbstractNLPModel

 meta     :: NLPModelMeta
 counters :: Counters
 mod      :: AbstractMPCCModel

 r        :: Float64
 s        :: Float64
 t        :: Float64
 tb       :: Float64

 function RlxMPCC(mod    :: AbstractMPCCModel,
                  r      :: Float64,
                  s      :: Float64,
                  t      :: Float64,
                  tb     :: Float64;
				  name   :: String = "Generic",
				  lin    :: AbstractVector{<: Integer}=Int[])
  x = mod.meta.x0

  nvar, ncc = mod.meta.nvar, mod.meta.ncc
  new_lcon = vcat(mod.meta.lcon,
                  mod.meta.lccG,
                  mod.meta.lccH,
                  -Inf*ones(ncc))
  new_ucon = vcat(mod.meta.ucon,
                  Inf*ones(2*ncc),
                  zeros(ncc))
  ncon = maximum([length(new_lcon); length(new_ucon)])
  y0 = vcat(mod.meta.y0, zeros(3*ncc))

  nnzh = nvar * (nvar + 1) / 2
  nnzj = nvar * ncon
  nln = setdiff(1:ncon, lin)

  meta = NLPModelMeta(nvar, x0 = x, lvar = mod.meta.lvar, uvar = mod.meta.uvar,
                                    ncon = mod.meta.ncon+3*ncc, y0=y0,
                                    lcon = new_lcon, ucon = new_ucon,
								    nnzj=nnzj, nnzh=nnzh, lin=lin, nln=nln,
								    minimize=true, islp=false, name=name)

  if minimum([r,s,t,tb])<0.0
   throw(error("Domain error: (r,s,t) must be non-negative"))
  end

  return new(meta, Counters(), mod, r, s, t, tb)
 end
end

"""
Return the number of classical nonlinear constraints
"""
function get_ncon(rlxmpcc :: RlxMPCC) return rlxmpcc.meta.ncon - 3*rlxmpcc.ncc end

"""
update_rlx!(:: RlxMPCC, :: Float64, :: Float64, :: Float64, :: Float64)
update the parameters r,s,t,tb of the RlxMPCC.
"""
function update_rlx!(rlx :: RlxMPCC, r :: Float64, s :: Float64, t :: Float64, tb :: Float64)
 rlx.r, rlx.s, rlx.t, rlx.tb = r, s, t, tb
 return rlx
end

function obj(nlp :: RlxMPCC, x ::  AbstractVector)
 increment!(nlp, :neval_obj)
 return obj(nlp.mod, x)
end

function grad!(nlp :: RlxMPCC, x :: Vector, gx :: AbstractVector)
 increment!(nlp, :neval_grad)
 return grad!(nlp.mod, x, gx)
end

function cons!(nlp :: RlxMPCC, x :: AbstractVector, c :: AbstractVector)
  increment!(nlp, :neval_cons)
  if nlp.mod.meta.ncc > 0
   Gx, Hx = consG(nlp.mod, x), consH(nlp.mod, x)
   c[1:nlp.meta.ncon] .= vcat(cons_nl(nlp.mod, x),
                              Gx .+ nlp.tb,
	  						  Hx .+ nlp.tb,
							  phi(Gx, Hx, nlp.r, nlp.s, nlp.t))
  else
   c[1:nlp.meta.ncon] .= cons_nl(nlp.mod, x)
  end

  return c
end

function jac(nlp :: RlxMPCC, x :: AbstractVector)
 increment!(nlp, :neval_jac)
 ncc, mod = nlp.mod.meta.ncc, nlp.mod

 if nlp.mod.meta.ncc > 0
  JGx, JHx = jacG(mod,x),   jacH(mod,x)
  Gx,  Hx  = consG(mod, x), consH(mod, x)
  dp = dphi(Gx, Hx, nlp.r, nlp.s, nlp.t)
  A = vcat(jac_nl(mod,x),
           JGx, JHx,
           diagm(0 => dp[1:ncc]) * JGx + diagm(0 => dp[ncc+1:2*ncc]) * JHx)
  else
   A = jac_nl(mod,x)
  end

 return A
end

function jac_coord!(nlp :: RlxMPCC, x :: AbstractVector, vals ::AbstractVector)
 Jx = jac(nlp, x) #findnz(jac(nlp.mod, x))
 m, n = size(Jx)
 vals[1 : n*m] .= Jx[:]
 return vals
end

function jac_structure!(nlp :: RlxMPCC,
	                   rows :: AbstractVector{<:Integer},
					   cols :: AbstractVector{<:Integer})
 m, n = nlp.meta.ncon, nlp.meta.nvar
 I = ((i,j) for i = 1:m, j = 1:n)
 rows[1 : n*m] .= getindex.(I, 1)[:]
 cols[1 : n*m] .= getindex.(I, 2)[:]
 return rows, cols
end

function jprod!(nlp :: RlxMPCC,
	            x   :: AbstractVector,
				v   :: AbstractVector,
				Jv  :: AbstractVector)
  increment!(nlp, :neval_jprod)
  Jv .= jac(nlp, x) * v
  return Jv
end

function jtprod!(nlp :: RlxMPCC,
	             x   :: AbstractVector,
				 v   :: AbstractVector,
				 Jtv :: AbstractVector)
  increment!(nlp, :neval_jtprod)
  Jtv .= jac(nlp,x)' * v
  return Jtv
end

function hess(nlp        :: RlxMPCC,
	          x          :: AbstractVector;
			  obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  Hx = hess(nlp.mod, x, obj_weight = obj_weight)
  return tril(Hx)
end

function hess(nlp        :: RlxMPCC,
	          x          :: AbstractVector,
			  y          :: AbstractVector;
			  obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  ncc, nlc, nvar = nlp.mod.meta.ncc, nlp.mod.meta.ncon, nlp.mod.meta.nvar

  if ncc > 0
   JGx, JHx = jacG(nlp.mod, x),   jacH(nlp.mod, x)
   Gx,  Hx  = consG(nlp.mod, x), consH(nlp.mod, x)

   lG, lH = y[nlc+1:nlc+ncc], y[nlc+ncc+1:nlc+2*ncc]
   lphi   = y[nlc+2*ncc+1:nlc+3*ncc]

   dgg, dgh, dhh = ddphi(Gx, Hx, nlp.r, nlp.s, nlp.t)
   dp            =  dphi(Gx, Hx, nlp.r, nlp.s, nlp.t)

   ny = vcat(y[1:nlc], lG - lphi .* dp[1:ncc], lH - lphi .* dp[ncc+1:2*ncc])
  else
   ny = y[1:nlc]
  end

  Hx = hess(nlp.mod, x, ny, obj_weight = obj_weight)

  if ncc > 0
   hess_P = JGx' * (diagm(0 => (dgg .* lphi)) * JGx)
          + JHx' * (diagm(0 => (dgh .* lphi)) * JHx)
	 	  + JGx' * (diagm(0 => (dgh .* lphi)) * JHx)
	  	  + JHx' * (diagm(0 => (dgh .* lphi)) * JGx)
   #for i=1:ncc
   # Hx += lphi[i] * (JGx[i,:] * JHx[i,:]' + JHx[i,:] * JGx[i,:]')
   #end
   Hx += hess_P
  end

  return tril(Hx)
end

function hess_structure!(nlp  :: RlxMPCC,
	                     rows :: AbstractVector{<: Integer},
						 cols :: AbstractVector{<: Integer})
  n = nlp.meta.nvar
  I = ((i,j) for i = 1:n, j = 1:n if i ≥ j)
  rows .= getindex.(I, 1)
  cols .= getindex.(I, 2)
  return rows, cols
end

function hess_coord!(nlp        :: RlxMPCC,
	                 x          :: AbstractVector,
					 vals       :: AbstractVector;
					 obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  Hx = hess(nlp, x, obj_weight=obj_weight)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hess_coord!(nlp        :: RlxMPCC,
	                 x          :: AbstractVector,
					 y          :: AbstractVector,
					 vals       :: AbstractVector;
					 obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hess)
  Hx = hess(nlp, x, y, obj_weight=obj_weight)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function hprod!(nlp        :: RlxMPCC,
	            x          :: AbstractVector,
				v          :: AbstractVector,
				Hv         :: AbstractVector;
				obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hprod)
  Hx = hess(nlp, x, obj_weight = obj_weight)
  Hv .= (Hx + Hx' - diagm(0 => diag(Hx))) * v
  return Hv
end

function hprod!(nlp        :: RlxMPCC,
	            x          :: AbstractVector,
				y          :: AbstractVector,
				v          :: AbstractVector,
				Hv         :: AbstractVector;
				obj_weight :: Real = one(eltype(x)))
  increment!(nlp, :neval_hprod)
  Hx = hess(nlp, x, y, obj_weight = obj_weight)
  Hv .= (Hx + Hx' - diagm(0 => diag(Hx))) * v
  return Hv
end

"""
Return the violation of the constraints
lb   <= x    <= ub,
lc   <= c(x) <= uc,
lccG <= G(x) + tb,
lccH <= H(x) + tb,
Phi(G(x),H(x),r s, t) <= 0.
"""
function viol(nlp :: RlxMPCC, x :: AbstractVector)

 mod, n, ncc = nlp.mod, nlp.mod.meta.nvar, nlp.mod.meta.ncc
 x = length(x) == n ? x : x[1:n]

 feas_x = vcat(max.(mod.meta.lvar-x, 0), max.(x-mod.meta.uvar, 0))

 if mod.meta.ncon > 0

  c = cons_nl(mod, x)
  feas_c = vcat(max.(mod.meta.lcon-c, 0), max.(c-mod.meta.ucon, 0))

 else

  feas_c = Float64[]

 end

 if ncc > 0

  Gx, Hx = consG(mod, x), consH(mod, x)

  feas_cp = vcat(max.(mod.meta.lccG - Gx .- nlp.tb, 0),
                 max.(mod.meta.lccH - Hx .- nlp.tb, 0))
  feas_cc = max.(phi(Gx, Hx, nlp.r, nlp.s, nlp.t), 0)
 else
  feas_cp = Float64[]
  feas_cc = Float64[]
 end

 return vcat(feas_x, feas_c, feas_cp, feas_cc)
end


text:
# # Incompressible Elasticity
#
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`incompressible_elasticity.ipynb`](@__NBVIEWER_ROOT_URL__/examples/incompressible_elasticity.ipynb)

# ## Introduction
#
# Mixed elements can be used to overcome locking when the material becomes
# incompressible. However, for an element to be stable, it needs to fulfill
# the LBB condition.
# In this example we will consider two different element formulations
# - linear displacement with linear pressure approximation (does *not* fulfill LBB)
# - quadratic displacement with linear pressure approximation (does fulfill LBB)
# The quadratic/linear element is also known as the Taylor-Hood element.
# We will consider Cook's Membrane with an applied traction on the right hand side.
#-
# ## Commented Program
#
#md # The full program, without comments, can be found in the next
#md # [section](@ref incompressible_elasticity-plain-program).
using Ferrite
using BlockArrays, SparseArrays, LinearAlgebra

# First we generate a simple grid, specifying the 4 corners of Cooks membrane.
function create_cook_grid(nx, ny)
    corners = [Vec{2}((0.0,   0.0)),
               Vec{2}((48.0, 44.0)),
               Vec{2}((48.0, 60.0)),
               Vec{2}((0.0,  44.0))]
    grid = generate_grid(Triangle, (nx, ny), corners);
    ## facesets for boundary conditions
    addfaceset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0);
    addfaceset!(grid, "traction", x -> norm(x[1]) ≈ 48.0);
    return grid
end;

# Next we define a function to set up our cell- and facevalues.
function create_values(interpolation_u, interpolation_p)
    ## quadrature rules
    qr      = QuadratureRule{2,RefTetrahedron}(3)
    face_qr = QuadratureRule{1,RefTetrahedron}(3)

    ## geometric interpolation
    interpolation_geom = Lagrange{2,RefTetrahedron,1}()

    ## cell and facevalues for u
    cellvalues_u = CellVectorValues(qr, interpolation_u, interpolation_geom)
    facevalues_u = FaceVectorValues(face_qr, interpolation_u, interpolation_geom)

    ## cellvalues for p
    cellvalues_p = CellScalarValues(qr, interpolation_p, interpolation_geom)

    return cellvalues_u, cellvalues_p, facevalues_u
end;

# We create a DofHandler, with two fields, `:u` and `:p`,
# with possibly different interpolations
function create_dofhandler(grid, ipu, ipp)
    dh = DofHandler(grid)
    push!(dh, :u, 2, ipu) # displacement
    push!(dh, :p, 1, ipp) # pressure
    close!(dh)
    return dh
end;

# We also need to add Dirichlet boundary conditions on the `"clamped"` faceset.
# We specify a homogeneous Dirichlet bc on the displacement field, `:u`.
function create_bc(dh)
    dbc = ConstraintHandler(dh)
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "clamped"), (x,t) -> zero(Vec{2}), [1,2]))
    close!(dbc)
    t = 0.0
    update!(dbc, t)
    return dbc
end;

# The material is linear elastic, which is here specified by the shear and bulk moduli
struct LinearElasticity{T}
    G::T
    K::T
end

# Now to the assembling of the stiffness matrix. This mixed formulation leads to a blocked
# element matrix. Since Ferrite does not force us to use any particular matrix type we will
# use a `PseudoBlockArray` from `BlockArrays.jl`.
function doassemble(cellvalues_u::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim},
                    facevalues_u::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,
                    dh::DofHandler, mp::LinearElasticity) where {dim}

    f = zeros(ndofs(dh))
    assembler = start_assemble(K, f)
    nu = getnbasefunctions(cellvalues_u)
    np = getnbasefunctions(cellvalues_p)

    fe = PseudoBlockArray(zeros(nu + np), [nu, np]) # local force vector
    ke = PseudoBlockArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix

    ## traction vector
    t = Vec{2}((0.0, 1/16))
    ## cache ɛdev outside the element routine to avoid some unnecessary allocations
    ɛdev = [zero(SymmetricTensor{2, dim}) for i in 1:getnbasefunctions(cellvalues_u)]

    for cell in CellIterator(dh)
        fill!(ke, 0)
        fill!(fe, 0)
        assemble_up!(ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)
        assemble!(assembler, celldofs(cell), fe, ke)
    end

    return K, f
end;

# The element routine integrates the local stiffness and force vector for all elements.
# Since the problem results in a symmetric matrix we choose to only assemble the lower part,
# and then symmetrize it after the loop over the quadrature points.
function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)

    n_basefuncs_u = getnbasefunctions(cellvalues_u)
    n_basefuncs_p = getnbasefunctions(cellvalues_p)
    u▄, p▄ = 1, 2
    reinit!(cellvalues_u, cell)
    reinit!(cellvalues_p, cell)

    ## We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
    @inbounds for q_point in 1:getnquadpoints(cellvalues_u)
        for i in 1:n_basefuncs_u
            ɛdev[i] = dev(symmetric(shape_gradient(cellvalues_u, q_point, i)))
        end
        dΩ = getdetJdV(cellvalues_u, q_point)
        for i in 1:n_basefuncs_u
            divδu = shape_divergence(cellvalues_u, q_point, i)
            δu = shape_value(cellvalues_u, q_point, i)
            for j in 1:i
                Ke[BlockIndex((u▄, u▄), (i, j))] += 2 * mp.G * ɛdev[i] ⊡ ɛdev[j] * dΩ
            end
        end

        for i in 1:n_basefuncs_p
            δp = shape_value(cellvalues_p, q_point, i)
            for j in 1:n_basefuncs_u
                divδu = shape_divergence(cellvalues_u, q_point, j)
                Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ
            end
            for j in 1:i
                p = shape_value(cellvalues_p, q_point, j)
                Ke[BlockIndex((p▄, p▄), (i, j))] += - 1/mp.K * δp * p * dΩ
            end

        end
    end

    symmetrize_lower!(Ke)

    ## We integrate the Neumann boundary using the facevalues.
    ## We loop over all the faces in the cell, then check if the face
    ## is in our `"traction"` faceset.
    @inbounds for face in 1:nfaces(cell)
        if onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "traction")
            reinit!(facevalues_u, cell, face)
            for q_point in 1:getnquadpoints(facevalues_u)
                dΓ = getdetJdV(facevalues_u, q_point)
                for i in 1:n_basefuncs_u
                    δu = shape_value(facevalues_u, q_point, i)
                    fe[i] += (δu ⋅ t) * dΓ
                end
            end
        end
    end
end

function symmetrize_lower!(K)
    for i in 1:size(K,1)
        for j in i+1:size(K,1)
            K[i,j] = K[j,i]
        end
    end
end;

# Now we have constructed all the necessary components, we just need a function
# to put it all together.

function solve(ν, interpolation_u, interpolation_p)
    ## material
    Emod = 1.
    Gmod = Emod / 2(1 + ν)
    Kmod = Emod * ν / ((1+ν) * (1-2ν))
    mp = LinearElasticity(Gmod, Kmod)

    ## grid, dofhandler, boundary condition
    n = 50
    grid = create_cook_grid(n, n)
    dh = create_dofhandler(grid, interpolation_u, interpolation_p)
    dbc = create_bc(dh)

    ## cellvalues
    cellvalues_u, cellvalues_p, facevalues_u = create_values(interpolation_u, interpolation_p)

    ## assembly and solve
    K = create_sparsity_pattern(dh);
    K, f = doassemble(cellvalues_u, cellvalues_p, facevalues_u, K, grid, dh, mp);
    apply!(K, f, dbc)
    u = Symmetric(K) \ f;

    ## export
    filename = "cook_" * (isa(interpolation_u, Lagrange{2,RefTetrahedron,1}) ? "linear" : "quadratic") *
                         "_linear"
    vtk_grid(filename, dh) do vtkfile
        vtk_point_data(vtkfile, dh, u)
    end
    return u
end

# All that is left is to solve the problem. We choose a value of Poissons
# ratio that is near incompressibility -- $ν = 0.5$ -- and thus expect the
# linear/linear approximation to return garbage, and the quadratic/linear
# approximation to be stable.
linear    = Lagrange{2,RefTetrahedron,1}()
quadratic = Lagrange{2,RefTetrahedron,2}()

u1 = solve(0.4999999, linear, linear)
u2 = solve(0.4999999, quadratic, linear);

## test the result                 #src
using Test                         #src
@test norm(u2) ≈ 919.2122668839389 #src

#md # ## [Plain Program](@id incompressible_elasticity-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [incompressible_elasticity.jl](incompressible_elasticity.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
#these tests are used for test DenseMultiVector under both serial and MPI comms

function denseMultiVectorTests(comm::Comm)
    #number of elements in vectors
    n = 8

    pid = myPid(comm)
    nProcs = numProc(comm)

    curMap = BlockMap(nProcs*n, n, comm)

    # test basic construction with setting data to zeros
    vect = DenseMultiVector{Float64}(curMap, 3, true)
    @test n == localLength(vect)
    @test nProcs*n == globalLength(vect)
    @test 3 == numVectors(vect)
    @test curMap == getMap(vect)
    @test zeros(Float64, (n, 3)) == vect.data

    # test basic construction without setting data to zeros
    vect = DenseMultiVector{Float64}(curMap, 3, false)
    @test n == localLength(vect)
    @test nProcs*n == globalLength(vect)
    @test 3 == numVectors(vect)
    @test curMap == getMap(vect)

    # test wrapper constructor
    arr = Array{Float64, 2}(undef, n, 3)
    vect = DenseMultiVector(curMap, arr)
    @test n == localLength(vect)
    @test nProcs*n == globalLength(vect)
    @test 3 == numVectors(vect)
    @test curMap == getMap(vect)
    @test arr === vect.data

    # test copy
    vect2 = copy(vect)
    @test n == localLength(vect)
    @test nProcs*n == globalLength(vect)
    @test 3 == numVectors(vect)
    @test curMap == getMap(vect)
    @test vect.data == vect2.data
    @test vect.data !== vect2.data #ensure same contents, but different address

    vect2 = DenseMultiVector{Float64}(curMap, 3, false)
    @test vect2 === copyto!(vect2, vect)
    @test localLength(vect) == localLength(vect2)
    @test globalLength(vect) == globalLength(vect2)
    @test numVectors(vect) == numVectors(vect2)
    @test getMap(vect) == getMap(vect2)
    @test vect.data == vect2.data
    @test vect.data !== vect2.data

    # test scale and scale!
    vect = DenseMultiVector(curMap, ones(Float64, n, 3))
    @test vect === scale!(vect, pid*5.0)
    @test pid*5*ones(Float64, (n, 3)) == vect.data

    increase = pid*nProcs

    vect = DenseMultiVector(curMap, ones(Float64, n, 3))
    @test vect isa DenseMultiVector
    @test vect === scale!(vect, increase .+ [2.0, 3.0, 4.0])
    @test hcat( (increase+2)*ones(Float64, n),
                (increase+3)*ones(Float64, n),
                (increase+4)*ones(Float64, n))  == vect.data

    for i = 1:3
        act = i .+ 1 .+ repeat(Float64[increase], inner=n)
        @test act == getVectorView(vect, i)
        @test act == getVectorCopy(vect, i)
    end

    #test dot
    vect = DenseMultiVector(curMap, ones(Float64, n, 3))
    @test fill(n*nProcs, (1, 3)) == dot(vect, vect)

    #test fill!
    fill!(vect, 8)
    @test 8*ones(Float64, (n, 3)) == vect.data

    #test broadcast
    vect = DenseMultiVector(curMap, ones(Float64, n, 3))
    fours = fill(4, (n, 3))
    vect2 = vect .+ fours
    @test fill(5.0, (n, 3)) == vect2.data
    @test fill(1.0, (n, 3)) == vect.data

    @. vect2 = vect * 3
    @test fill(3.0, (n, 3)) == vect2.data
    @test fill(1.0, (n, 3)) == vect.data

    #test comm reduce
    arr = (10^pid)*ones(Float64, n, 3)
    vect = DenseMultiVector(BlockMap(n, n, comm), arr)
    commReduce(vect)
    @test sum(10^i for i in 1:nProcs)*ones(Float64, n, 3) == vect.data

    #test norm
    arr = ones(Float64, n, 3)
    vect = DenseMultiVector(curMap, arr)
    @test [sqrt(n*nProcs), sqrt(n*nProcs), sqrt(n*nProcs)]' == norm(vect, 2)
    @test [(n*nProcs)^(1/3), (n*nProcs)^(1/3), (n*nProcs)^(1/3)]' == norm(vect, 3)

    arr = 2*ones(Float64, n, 3)
    vect = DenseMultiVector(curMap, arr)
    @test [sqrt(4*n*nProcs), sqrt(4*n*nProcs), sqrt(4*n*nProcs)]' == norm(vect, 2)
    @test [(8*n*nProcs)^(1/3), (8*n*nProcs)^(1/3), (8*n*nProcs)^(1/3)]' == norm(vect, 3)

    #test imports/exports
    source = DenseMultiVector(curMap,
        Array{Float64, 2}(reshape(collect(1:(3*n)), (n, 3))))
    target = DenseMultiVector{Float64}(curMap, 3, false)
    impor = Import(curMap, curMap)
    doImport(source, target, impor, REPLACE)
    @test reshape(Array{Float64, 1}(collect(1:(3*n))), (n, 3)) == target.data

    source = DenseMultiVector(curMap,
        Array{Float64, 2}(reshape(collect(1:(3*n)), (n, 3))))
    target = DenseMultiVector{Float64}(curMap, 3, false)
    expor = Export(curMap, curMap)
    doExport(source, target, expor, REPLACE)
    @test reshape(Array{Float64, 1}(collect(1:(3*n))), (n, 3)) == target.data

    source = DenseMultiVector(curMap,
        Array{Float64, 2}(reshape(collect(1:(3*n)), (n, 3))))
    target = DenseMultiVector{Float64}(curMap, 3, false)
    impor = Import(curMap, curMap)
    doExport(source, target, impor, REPLACE)
    @test reshape(Array{Float64, 1}(collect(1:(3*n))), (n, 3)) == target.data

    source = DenseMultiVector(curMap,
        Array{Float64, 2}(reshape(collect(1:(3*n)), (n, 3))))
    target = DenseMultiVector{Float64}(curMap, 3, false)
    expor = Export(curMap, curMap)
    doImport(source, target, expor, REPLACE)
    @test reshape(Array{Float64, 1}(collect(1:(3*n))), (n, 3)) == target.data

    #TODO create import expor tests to test non trivial case

end


text:
function is_leap_year(year::Int)
    if year % 400 == 0
        return true
    elseif year % 100 == 0
        return false
    elseif year % 4 == 0
        return true
    end
end

is_leap_year(2015)

max_stars_repo_name: Devlious/LearningJulia

max_issues_repo_name: Devlious/LearningJulia

max_forks_repo_name: Devlious/LearningJulia


text:
# Code for working with complexes,
# and complex-balanced systems

# (species,stoichiometric coefficients)
const SpeciesComplex = Tuple{IntVec,IntVec}

# composition matrix has dimensions n_species, n_complexes
function composition_matrix(n_species::Integer, cmplx::Vector{SpeciesComplex})
    Γ = spzeros(Int64,n_species, length(cmplx))
    for i in 1:length(cmplx)
        spec   = cmplx[i][1]
        stoich = cmplx[i][2]
        for j = 1 : length(spec)
            Γ[spec[j],i] += stoich[j]
        end
    end
    return Γ
end

"""
    complex_decomposition(rn::ReactionNetwork)

Return composition and incidence matrices of reaction net.
composition matrix has dimensions (`n_species`, `n_complexes`).
incidence matrix has dimensionality (`n_complex`, `n_reactions`),
the entries indicate if a complex is the reactant or product.
"""
function complex_decomposition(rn::ReactionNetwork)
    n_species = rn.n_species
    reactions = rn.reactions

    c = SpeciesComplex[]
    inc_mat = spzeros(Int64,2*length(reactions),length(reactions))  # we can have at most 2*rho complexes
    n_complex = 0

    idx = 0

    for i = 1:length(reactions)
        r = reactions[i]

        # Reactant complex
        rc = (r.reactants, r.stoichr)

        # Find complex in list of complexes
        idx = findfirst(x->x==rc,c)
        if idx==nothing  # If not already in, put it in
            append!(c, [rc])
            n_complex += 1
            idx = n_complex
        end

        inc_mat[idx, i] += -1

        # Product complex
        pc = (r.products, r.stoichp)

        # Find complex in list of complexes
        idx = findfirst(x->x==pc,c)
        if idx==nothing  # If not already in, put it in
            append!(c, [pc])
            n_complex += 1
            idx = n_complex
        end

        inc_mat[idx, i] +=  1
    end

    return composition_matrix(n_species, c), inc_mat[1:n_complex,:]
end

# Compute deficiency of reaction network
function deficiency(rn::ReactionNetwork)
    comp_mat, inc_mat = complex_decomposition(rn)

    return n_cycles(rn) - (size(inc_mat,2)-rank(full(inc_mat)))
end

# Take the species concentration vector,
# And return a vector of concentrations of complexes
# z: species concentration
complex_concentration(z, composition_matrix) = vec(exp((log(z)'*composition_matrix)))

# Create a reaction system for the complexes
function complex_basis_reactions(reactions::Unsigned, incidence_matrix)
    r = Reaction[]
    n_complexes = size(incidence_matrix,1)
    names = ["(C$x)" for x=1:n_complexes]
    for i = 1 : length(reactions)
        nzind,nzval = findnz(incidence_matrix[:,i])

        reactants = nzval[1] < 0 ? [nzind[1]] : [nzind[2]]
        products  = nzval[2] > 0 ? [nzind[2]] : [nzind[1]]

        append!(r, [Reaction(reactants,[1],products,[1], reactions[i].kf, reactions[i].kr, names)])
    end
    return r
end