text:
# new Sandshark example
# add more julia processes
# nprocs() < 7 ? addprocs(7-nprocs()) : nothing

using KernelDensityEstimate, IncrementalInference
using Distributions

const TU = TransformUtils

## Step: Building the factor graph
fg = initfg()
# Add a central beacon with a prior
addVariable!(fg, :l1, ContinuousScalar)
# Pinger location is (0.6; -16)
addFactor!(fg, [:l1], IIF.Prior( Normal(0.6, 0.1) ))

epochs = 50:2:100

index=0
for ep in epochs
    curvar = Symbol("x$index")
    addVariable!(fg, curvar, ContinuousScalar)

    if ep != epochs[1]
      # Odo factor x(i-1) -> xi
      addFactor!(fg, [Symbol("x$(index-1)"); curvar], LinearConditional(Normal(1, 0.1)))
    else
      # Prior to the first pose location (a "GPS" prior)
      addFactor!(fg, [curvar], IIF.Prior(Normal(5,5) ))
    end
    # Heading partial prior
    # addFactor!(fg, [curvar], RoME.PartialPriorYawPose2(Normal(interp_yaw(ep), deg2rad(3))))
    index+=1
end

# Just adding the first one...
addFactor!(fg, [:x0; :l1], LinearConditional(Normal(1.0, 0.1)))

addFactor!(fg, [:x5; :l1], LinearConditional(Normal(1.0, 0.1)))

addFactor!(fg, [:x10; :l1], LinearConditional(Normal(1.0, 0.1)))

addFactor!(fg, [:x13; :l1], LinearConditional(Normal(1.0, 0.1)))

addFactor!(fg, [:x15; :l1], LinearConditional(Normal(1.0, 0.1)))

addFactor!(fg, [:x17; :l1], LinearConditional(Normal(1.0, 0.1)))
addFactor!(fg, [:x18; :l1], LinearConditional(Normal(1.0, 0.1)))
addFactor!(fg, [:x19; :l1], LinearConditional(Normal(1.0, 0.1)))
addFactor!(fg, [:x20; :l1], LinearConditional(Normal(1.0, 0.1)))
addFactor!(fg, [:x21; :l1], LinearConditional(Normal(1.0, 0.1))) # breaks it!
addFactor!(fg, [:x22; :l1], LinearConditional(Normal(1.0, 0.1)))
addFactor!(fg, [:x23; :l1], LinearConditional(Normal(1.0, 0.1)))
addFactor!(fg, [:x24; :l1], LinearConditional(Normal(1.0, 0.1)))
addFactor!(fg, [:x25; :l1], LinearConditional(Normal(1.0, 0.1)))

plotPose(fg, [:x21])
plotKDE(fg, [:x20, :x21, :x22, :x23, :x24])
# plotKDE(ppbrDict[epochs[26]].range)
#
# plotKDE(ppbrDict[epochs[23]].bearing)
plotKDE(ppbrDict[epochs[23]].range)

writeGraphPdf(fg, engine="dot")

ensureAllInitialized!(fg)
solveTree!(fg)

drawPosesLandms(fg)

# saveDFG(fg, file="presolve_$t.jld")
IIF.solveTree!(fg) #, N=100
# saveDFG(fg, file="postsolve_$t.jld")

# pl = drawPoses(fg, spscale=2.75) # Just for odo plot
# Roll again for inspiration check
## PLOT BEAM PATTERNS
Gadfly.push_theme(:default)
pla = drawPosesLandmarksAndOdo(fg, ppbrDict, navkeys, X, Y, lblkeys, lblX, lblY)
Gadfly.draw(PDF("sandshark-beacon_$t.pdf", 12cm, 15cm), pla)
Gadfly.draw(PNG("sandshark-beacon_$t.png", 12cm, 15cm), pla)

## resolve funky 21

# shows all the proposals based on the clique in the tree -- similar to plotLocalProduct on factor graph
plotTreeProduct(fg, tree, :x25)
#
# using IncrementalInference
#
#
# import RoMEPlotting: spyCliqMat
#
# function spyCliqMat(cliq::Graphs.ExVertex; showmsg=true)
#   mat = deepcopy(IIF.getCliqMat(cliq, showmsg=showmsg))
#   # TODO -- add improved visualization here, iter vs skip
#   mat = map(Float64, mat)*2.0-1.0
#   numlcl = size(IIF.getCliqAssocMat(cliq),1)
#   mat[(numlcl+1):end,:] *= 0.9
#   mat[(numlcl+1):end,:] -= 0.1
#   numfrtl1 = floor(Int,length(cliq.attributes["data"].frontalIDs)+1)
#   mat[:,numfrtl1:end] *= 0.9
#   mat[:,numfrtl1:end] -= 0.1
#   @show cliq.attributes["data"].itervarIDs
#   @show cliq.attributes["data"].directvarIDs
#   @show cliq.attributes["data"].msgskipIDs
#   @show cliq.attributes["data"].directFrtlMsgIDs
#   @show cliq.attributes["data"].directPriorMsgIDs
#   sp = Gadfly.spy(mat)
#   push!(sp.guides, Gadfly.Guide.title("$(cliq.attributes["label"]) || $(cliq.attributes["data"].frontalIDs) :$(cliq.attributes["data"].separatorIDs)"))
#   push!(sp.guides, Gadfly.Guide.xlabel("fmcmcs $(cliq.attributes["data"].itervarIDs)"))
#   push!(sp.guides, Gadfly.Guide.ylabel("lcl=$(numlcl) || msg=$(size(IIF.getCliqMsgMat(cliq),1))" ))
#   return sp
# end
# function spyCliqMat(bt::BayesTree, lbl::Symbol; showmsg=true)
#   spyCliqMat(IIF.getClique(bt,lbl), showmsg=showmsg)
# end
#

tree = wipeBuildNewTree!(fg, drawpdf=true, show=true)

import IncrementalInference: getCliqMat

sym = :x21
getClique(tree, sym).attributes["label"]
spyCliqMat(tree, sym)

syms = union(getCliqSymbols(tree, :x21)...)
findfirst(syms, :x21)

getData(getClique(tree, :x21))

0

####  DEBUG PPBRDict ====================

LL = [0.6; -16; 0]
function expectedLocalBearing(curPose, LL)
    v = [LL[1] - curPose[1]; LL[2] - curPose[2]]
    world = atan2(v[2], v[1])
    loc = rad2deg(TU.wrapRad(world - curPose[3]))
    while loc < 0
        loc += 360
    end
    return loc
end
# println("$(interp_x[ep]), $(interp_y[ep]), world yaw = $(interp_yaw[ep] * 180 / pi)")
# XX = [interp_x(ep); interp_y[ep]; interp_yaw[ep]]
# expectedLocalBearing(XX, LL)
# println("Expected local yaw = $(expectedLocalBearing(XX, LL)). Sonar local yaw = $(rad2deg(getKDEMax(ppbr.bearing)[1]))")
# xyt = se2vee(SE2(XX[1:3]) \ SE2(LL))
# bear = rad2deg(TU.wrapRad(atan2(xyt[2],xyt[1]) -XX[3]))

bearExp = []
bearAcoustics = []
for ep in epochs
    XX = [interp_x(ep); interp_y[ep]; interp_yaw[ep]]
    ppbr = ppbrDict[ep]
    push!(bearExp, expectedLocalBearing(XX, LL))
    push!(bearAcoustics, rad2deg(getKDEMax(ppbr.bearing)[1]))
    println("local beacon heading: Expected = $(expectedLocalYaw(XX, LL)). Sonar = $(rad2deg(getKDEMax(ppbr.bearing)[1]))")
end

layers = []
push!(layers, Gadfly.layer(x=1:length(bearExp), y=bearExp, Geom.path())[1], Theme(default_color=color("red")))
push!(layers, Gadfly.layer(x=1:length(bearAcoustics), y=bearAcoustics, Geom.path(), Theme(default_color=color("blue"))))
Gadfly.plot(layers...)
# end

####################################

######################################################

# pl = drawPoses(fg, spscale=2.75)
# drawPosesLandms(fg, spscale=0.75) #Means so we don't run into MM == Union() || Dict{} in
# You rolled 20!

# Roll again for debug check
## Debugging strange headings
# posePts = get2DPoseMeans(fg)
# landPts = get2DLandmMeans(fg)
## Debugging landmark bearing range

# ppbrDict[epoch_slice[1]]
# getSample(ppbrDict[epoch_slice[1]],100)
#
#
# addVariable!(fg, :l1, Point2)
# addFactor!(fg, [:x0; :l1], ppbrDict[epoch_slice[1]])
#
# ls(fg, :l1)
# pts = IIF.approxConv(fg, :x0l1f1, :l1)
#
# fct = getData(getVert(fg, :x0l1f1, nt=:fnc))
#
# fct.fnc.zDim

# dev
# ep = epochs[1]
# # azidata[ep][:,1]
# pll = layerBeamPatternRose(ppbrDict[ep].bearing)
# Gadfly.plot(pll...)
#
# pll = layerBeamPatternRose(ppbrDict[ep].bearing, wRr=TU.R(pi/2))
# Gadfly.plot(pll...)

#
# diff(epochs[[end;1]])

## GADFLY EXAMPLE
# plotLocalProduct(fg, :x49, dims=[1;2])
# stuff = IncrementalInference.localProduct(fg, :x49)
# plotKDE(stuff[2], dims=[1;2], levels=10)
#
# # STUFF
# fsym = :x49l1f1
# const TU = TransformUtils
# XX, LL = (KDE.getKDEMax.(IIF.getBelief.(fg, IIF.lsf(fg, fsym)))...)
# @show xyt = se2vee(SE2(XX[1:3]) \ SE2([LL[1:2];0.0]))
# bear= rad2deg(TU.wrapRad(atan2(-xyt[2],xyt[1]) -XX[3]))
# b = IncrementalInference.getData(fg, fsym, nt=:fnc).fnc.usrfnc!
# b.bearing
# p = plotKDE(b.range)
# Gadfly.draw(PDF("multimodal_bearing_range.pdf", 12cm, 15cm), p)
#
# # SAVE THIS PLOT
# # epochs = timestamps[11:2:200]
# g = getBelief.(fg, [:l1, :x49])
# m1 = KDE.marginal(g[1], [1;2])
# m2 = KDE.marginal(g[2], [1;2])
# norm(diff(KDE.getKDEMax.([m1; m2])))
#
# kde!(rand(Normal(4.0, 0.1),100))
#
# ###############
#
# PL = []
#
# push!(PL, Gadfly.layer(x=collect(1:10), y=randn(10), Geom.path, Theme(default_color=colorant"red"))[1])
# push!(PL, Gadfly.layer(x=collect(1:10), y=randn(10), Geom.path, Theme(default_color=colorant"green"))[1])
# push!(PL, Coord.Cartesian(xmin=-1.0,xmax=11, ymin=-5.0, ymax=5.0))
# push!(PL, Guide.ylabel("test"))
# push!(PL, Guide.ylabel("test"))
# pl = drawPoses(fg)
# push!(pl.layers, )
#
# pl = drawPoses(fg)
# pl = Gadfly.plot(pl.layers...)
#
#
#
# pl = Gadfly.plot(PL...)
#
# @show fieldnames(pl)
# push!(pl.layers, Gadfly.layer(x=collect(1:10), y=randn(10), Geom.path, Theme(default_color=colorant"magenta"))[1] )
#
# pl
# # vstasck, hstack
#
#
# Gadfly.push_theme(:default)
# # SVG, PDF, PNG
# Gadfly.draw(PDF("/tmp/testfig.pdf", 12cm, 8cm), pl)
# run(`evince /tmp/testfig.pdf`)
#
#
# ### PATH + Poses
#
# pl = drawPosesLandms(fg, spscale=2.5)
# # Add odo
# navdf = DataFrame(
#   ts = navkeys,
#   x = X,
#   y = Y
# )
# # pl = Gadfly.layer(navdf, x=:x, y=:y, Geom.path())
# push!(pl.layers, Gadfly.layer(navdf, x=:x, y=:y, Geom.path())[1])
# Gadfly.plot(pl.layers)

fsym = :x49l1f1
const TU = TransformUtils
XX, LL = (KDE.getKDEMax.(IIF.getBelief.(fg, IIF.lsf(fg, fsym)))...)
@show xyt = se2vee(SE2(XX[1:3]) \ SE2([LL[1:2];0.0]))
bear= rad2deg(TU.wrapRad(atan2(-xyt[2],xyt[1]) -XX[3]))

#


text:
# Julia wrapper for header: spirv_cross_c.h
# Automatically generated using Clang.jl

function spvc_get_version(major, minor, patch)
    ccall((:spvc_get_version, libspirv_cross), Cvoid, (Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}), major, minor, patch)
end

function spvc_get_commit_revision_and_timestamp()
    ccall((:spvc_get_commit_revision_and_timestamp, libspirv_cross), Cstring, ())
end

function spvc_msl_vertex_attribute_init(attr)
    ccall((:spvc_msl_vertex_attribute_init, libspirv_cross), Cvoid, (Ptr{spvc_msl_vertex_attribute},), attr)
end

function spvc_msl_resource_binding_init(binding)
    ccall((:spvc_msl_resource_binding_init, libspirv_cross), Cvoid, (Ptr{spvc_msl_resource_binding},), binding)
end

function spvc_msl_get_aux_buffer_struct_version()
    ccall((:spvc_msl_get_aux_buffer_struct_version, libspirv_cross), UInt32, ())
end

function spvc_msl_constexpr_sampler_init(sampler)
    ccall((:spvc_msl_constexpr_sampler_init, libspirv_cross), Cvoid, (Ptr{spvc_msl_constexpr_sampler},), sampler)
end

function spvc_msl_sampler_ycbcr_conversion_init(conv)
    ccall((:spvc_msl_sampler_ycbcr_conversion_init, libspirv_cross), Cvoid, (Ptr{spvc_msl_sampler_ycbcr_conversion},), conv)
end

function spvc_hlsl_resource_binding_init(binding)
    ccall((:spvc_hlsl_resource_binding_init, libspirv_cross), Cvoid, (Ptr{spvc_hlsl_resource_binding},), binding)
end

function spvc_context_create(context)
    ccall((:spvc_context_create, libspirv_cross), spvc_result, (Ptr{spvc_context},), context)
end

function spvc_context_destroy(context)
    ccall((:spvc_context_destroy, libspirv_cross), Cvoid, (spvc_context,), context)
end

function spvc_context_release_allocations(context)
    ccall((:spvc_context_release_allocations, libspirv_cross), Cvoid, (spvc_context,), context)
end

function spvc_context_get_last_error_string(context)
    ccall((:spvc_context_get_last_error_string, libspirv_cross), Cstring, (spvc_context,), context)
end

function spvc_context_set_error_callback(context, cb, userdata)
    ccall((:spvc_context_set_error_callback, libspirv_cross), Cvoid, (spvc_context, spvc_error_callback, Ptr{Cvoid}), context, cb, userdata)
end

function spvc_context_parse_spirv(context, spirv, word_count, parsed_ir)
    ccall((:spvc_context_parse_spirv, libspirv_cross), spvc_result, (spvc_context, Ptr{SpvId}, Cint, Ptr{spvc_parsed_ir}), context, spirv, word_count, parsed_ir)
end

function spvc_context_create_compiler(context, backend, parsed_ir, mode, compiler)
    ccall((:spvc_context_create_compiler, libspirv_cross), spvc_result, (spvc_context, spvc_backend, spvc_parsed_ir, spvc_capture_mode, Ptr{spvc_compiler}), context, backend, parsed_ir, mode, compiler)
end

function spvc_compiler_get_current_id_bound(compiler)
    ccall((:spvc_compiler_get_current_id_bound, libspirv_cross), UInt32, (spvc_compiler,), compiler)
end

function spvc_compiler_create_compiler_options(compiler, options)
    ccall((:spvc_compiler_create_compiler_options, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_compiler_options}), compiler, options)
end

function spvc_compiler_options_set_bool(options, option, value)
    ccall((:spvc_compiler_options_set_bool, libspirv_cross), spvc_result, (spvc_compiler_options, spvc_compiler_option, spvc_bool), options, option, value)
end

function spvc_compiler_options_set_uint(options, option, value)
    ccall((:spvc_compiler_options_set_uint, libspirv_cross), spvc_result, (spvc_compiler_options, spvc_compiler_option, UInt32), options, option, value)
end

function spvc_compiler_install_compiler_options(compiler, options)
    ccall((:spvc_compiler_install_compiler_options, libspirv_cross), spvc_result, (spvc_compiler, spvc_compiler_options), compiler, options)
end

function spvc_compiler_compile(compiler, source)
    ccall((:spvc_compiler_compile, libspirv_cross), spvc_result, (spvc_compiler, Ptr{Cstring}), compiler, source)
end

function spvc_compiler_add_header_line(compiler, line)
    ccall((:spvc_compiler_add_header_line, libspirv_cross), spvc_result, (spvc_compiler, Cstring), compiler, line)
end

function spvc_compiler_require_extension(compiler, ext)
    ccall((:spvc_compiler_require_extension, libspirv_cross), spvc_result, (spvc_compiler, Cstring), compiler, ext)
end

function spvc_compiler_flatten_buffer_block(compiler, id)
    ccall((:spvc_compiler_flatten_buffer_block, libspirv_cross), spvc_result, (spvc_compiler, spvc_variable_id), compiler, id)
end

function spvc_compiler_variable_is_depth_or_compare(compiler, id)
    ccall((:spvc_compiler_variable_is_depth_or_compare, libspirv_cross), spvc_bool, (spvc_compiler, spvc_variable_id), compiler, id)
end

function spvc_compiler_hlsl_set_root_constants_layout(compiler, constant_info, count)
    ccall((:spvc_compiler_hlsl_set_root_constants_layout, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_hlsl_root_constants}, Cint), compiler, constant_info, count)
end

function spvc_compiler_hlsl_add_vertex_attribute_remap(compiler, remap, remaps)
    ccall((:spvc_compiler_hlsl_add_vertex_attribute_remap, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_hlsl_vertex_attribute_remap}, Cint), compiler, remap, remaps)
end

function spvc_compiler_hlsl_remap_num_workgroups_builtin(compiler)
    ccall((:spvc_compiler_hlsl_remap_num_workgroups_builtin, libspirv_cross), spvc_variable_id, (spvc_compiler,), compiler)
end

function spvc_compiler_hlsl_set_resource_binding_flags(compiler, flags)
    ccall((:spvc_compiler_hlsl_set_resource_binding_flags, libspirv_cross), spvc_result, (spvc_compiler, spvc_hlsl_binding_flags), compiler, flags)
end

function spvc_compiler_hlsl_add_resource_binding(compiler, binding)
    ccall((:spvc_compiler_hlsl_add_resource_binding, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_hlsl_resource_binding}), compiler, binding)
end

function spvc_compiler_hlsl_is_resource_used(compiler, model, set, binding)
    ccall((:spvc_compiler_hlsl_is_resource_used, libspirv_cross), spvc_bool, (spvc_compiler, SpvExecutionModel, UInt32, UInt32), compiler, model, set, binding)
end

function spvc_compiler_msl_is_rasterization_disabled(compiler)
    ccall((:spvc_compiler_msl_is_rasterization_disabled, libspirv_cross), spvc_bool, (spvc_compiler,), compiler)
end

function spvc_compiler_msl_needs_aux_buffer(compiler)
    ccall((:spvc_compiler_msl_needs_aux_buffer, libspirv_cross), spvc_bool, (spvc_compiler,), compiler)
end

function spvc_compiler_msl_needs_swizzle_buffer(compiler)
    ccall((:spvc_compiler_msl_needs_swizzle_buffer, libspirv_cross), spvc_bool, (spvc_compiler,), compiler)
end

function spvc_compiler_msl_needs_buffer_size_buffer(compiler)
    ccall((:spvc_compiler_msl_needs_buffer_size_buffer, libspirv_cross), spvc_bool, (spvc_compiler,), compiler)
end

function spvc_compiler_msl_needs_output_buffer(compiler)
    ccall((:spvc_compiler_msl_needs_output_buffer, libspirv_cross), spvc_bool, (spvc_compiler,), compiler)
end

function spvc_compiler_msl_needs_patch_output_buffer(compiler)
    ccall((:spvc_compiler_msl_needs_patch_output_buffer, libspirv_cross), spvc_bool, (spvc_compiler,), compiler)
end

function spvc_compiler_msl_needs_input_threadgroup_mem(compiler)
    ccall((:spvc_compiler_msl_needs_input_threadgroup_mem, libspirv_cross), spvc_bool, (spvc_compiler,), compiler)
end

function spvc_compiler_msl_add_vertex_attribute(compiler, attrs)
    ccall((:spvc_compiler_msl_add_vertex_attribute, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_msl_vertex_attribute}), compiler, attrs)
end

function spvc_compiler_msl_add_resource_binding(compiler, binding)
    ccall((:spvc_compiler_msl_add_resource_binding, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_msl_resource_binding}), compiler, binding)
end

function spvc_compiler_msl_add_discrete_descriptor_set(compiler, desc_set)
    ccall((:spvc_compiler_msl_add_discrete_descriptor_set, libspirv_cross), spvc_result, (spvc_compiler, UInt32), compiler, desc_set)
end

function spvc_compiler_msl_set_argument_buffer_device_address_space(compiler, desc_set, device_address)
    ccall((:spvc_compiler_msl_set_argument_buffer_device_address_space, libspirv_cross), spvc_result, (spvc_compiler, UInt32, spvc_bool), compiler, desc_set, device_address)
end

function spvc_compiler_msl_is_vertex_attribute_used(compiler, location)
    ccall((:spvc_compiler_msl_is_vertex_attribute_used, libspirv_cross), spvc_bool, (spvc_compiler, UInt32), compiler, location)
end

function spvc_compiler_msl_is_resource_used(compiler, model, set, binding)
    ccall((:spvc_compiler_msl_is_resource_used, libspirv_cross), spvc_bool, (spvc_compiler, SpvExecutionModel, UInt32, UInt32), compiler, model, set, binding)
end

function spvc_compiler_msl_remap_constexpr_sampler(compiler, id, sampler)
    ccall((:spvc_compiler_msl_remap_constexpr_sampler, libspirv_cross), spvc_result, (spvc_compiler, spvc_variable_id, Ptr{spvc_msl_constexpr_sampler}), compiler, id, sampler)
end

function spvc_compiler_msl_remap_constexpr_sampler_by_binding(compiler, desc_set, binding, sampler)
    ccall((:spvc_compiler_msl_remap_constexpr_sampler_by_binding, libspirv_cross), spvc_result, (spvc_compiler, UInt32, UInt32, Ptr{spvc_msl_constexpr_sampler}), compiler, desc_set, binding, sampler)
end

function spvc_compiler_msl_remap_constexpr_sampler_ycbcr(compiler, id, sampler, conv)
    ccall((:spvc_compiler_msl_remap_constexpr_sampler_ycbcr, libspirv_cross), spvc_result, (spvc_compiler, spvc_variable_id, Ptr{spvc_msl_constexpr_sampler}, Ptr{spvc_msl_sampler_ycbcr_conversion}), compiler, id, sampler, conv)
end

function spvc_compiler_msl_remap_constexpr_sampler_by_binding_ycbcr(compiler, desc_set, binding, sampler, conv)
    ccall((:spvc_compiler_msl_remap_constexpr_sampler_by_binding_ycbcr, libspirv_cross), spvc_result, (spvc_compiler, UInt32, UInt32, Ptr{spvc_msl_constexpr_sampler}, Ptr{spvc_msl_sampler_ycbcr_conversion}), compiler, desc_set, binding, sampler, conv)
end

function spvc_compiler_msl_set_fragment_output_components(compiler, location, components)
    ccall((:spvc_compiler_msl_set_fragment_output_components, libspirv_cross), spvc_result, (spvc_compiler, UInt32, UInt32), compiler, location, components)
end

function spvc_compiler_msl_get_automatic_resource_binding(compiler, id)
    ccall((:spvc_compiler_msl_get_automatic_resource_binding, libspirv_cross), UInt32, (spvc_compiler, spvc_variable_id), compiler, id)
end

function spvc_compiler_msl_get_automatic_resource_binding_secondary(compiler, id)
    ccall((:spvc_compiler_msl_get_automatic_resource_binding_secondary, libspirv_cross), UInt32, (spvc_compiler, spvc_variable_id), compiler, id)
end

function spvc_compiler_msl_add_dynamic_buffer(compiler, desc_set, binding, index)
    ccall((:spvc_compiler_msl_add_dynamic_buffer, libspirv_cross), spvc_result, (spvc_compiler, UInt32, UInt32, UInt32), compiler, desc_set, binding, index)
end

function spvc_compiler_msl_add_inline_uniform_block(compiler, desc_set, binding)
    ccall((:spvc_compiler_msl_add_inline_uniform_block, libspirv_cross), spvc_result, (spvc_compiler, UInt32, UInt32), compiler, desc_set, binding)
end

function spvc_compiler_get_active_interface_variables(compiler, set)
    ccall((:spvc_compiler_get_active_interface_variables, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_set}), compiler, set)
end

function spvc_compiler_set_enabled_interface_variables(compiler, set)
    ccall((:spvc_compiler_set_enabled_interface_variables, libspirv_cross), spvc_result, (spvc_compiler, spvc_set), compiler, set)
end

function spvc_compiler_create_shader_resources(compiler, resources)
    ccall((:spvc_compiler_create_shader_resources, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_resources}), compiler, resources)
end

function spvc_compiler_create_shader_resources_for_active_variables(compiler, resources, active)
    ccall((:spvc_compiler_create_shader_resources_for_active_variables, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_resources}, spvc_set), compiler, resources, active)
end

function spvc_resources_get_resource_list_for_type(resources, type, resource_list, resource_size)
    ccall((:spvc_resources_get_resource_list_for_type, libspirv_cross), spvc_result, (spvc_resources, spvc_resource_type, Ptr{Ptr{spvc_reflected_resource}}, Ptr{Cint}), resources, type, resource_list, resource_size)
end

function spvc_compiler_set_decoration(compiler, id, decoration, argument)
    ccall((:spvc_compiler_set_decoration, libspirv_cross), Cvoid, (spvc_compiler, SpvId, SpvDecoration, UInt32), compiler, id, decoration, argument)
end

function spvc_compiler_set_decoration_string(compiler, id, decoration, argument)
    ccall((:spvc_compiler_set_decoration_string, libspirv_cross), Cvoid, (spvc_compiler, SpvId, SpvDecoration, Cstring), compiler, id, decoration, argument)
end

function spvc_compiler_set_name(compiler, id, argument)
    ccall((:spvc_compiler_set_name, libspirv_cross), Cvoid, (spvc_compiler, SpvId, Cstring), compiler, id, argument)
end

function spvc_compiler_set_member_decoration(compiler, id, member_index, decoration, argument)
    ccall((:spvc_compiler_set_member_decoration, libspirv_cross), Cvoid, (spvc_compiler, spvc_type_id, UInt32, SpvDecoration, UInt32), compiler, id, member_index, decoration, argument)
end

function spvc_compiler_set_member_decoration_string(compiler, id, member_index, decoration, argument)
    ccall((:spvc_compiler_set_member_decoration_string, libspirv_cross), Cvoid, (spvc_compiler, spvc_type_id, UInt32, SpvDecoration, Cstring), compiler, id, member_index, decoration, argument)
end

function spvc_compiler_set_member_name(compiler, id, member_index, argument)
    ccall((:spvc_compiler_set_member_name, libspirv_cross), Cvoid, (spvc_compiler, spvc_type_id, UInt32, Cstring), compiler, id, member_index, argument)
end

function spvc_compiler_unset_decoration(compiler, id, decoration)
    ccall((:spvc_compiler_unset_decoration, libspirv_cross), Cvoid, (spvc_compiler, SpvId, SpvDecoration), compiler, id, decoration)
end

function spvc_compiler_unset_member_decoration(compiler, id, member_index, decoration)
    ccall((:spvc_compiler_unset_member_decoration, libspirv_cross), Cvoid, (spvc_compiler, spvc_type_id, UInt32, SpvDecoration), compiler, id, member_index, decoration)
end

function spvc_compiler_has_decoration(compiler, id, decoration)
    ccall((:spvc_compiler_has_decoration, libspirv_cross), spvc_bool, (spvc_compiler, SpvId, SpvDecoration), compiler, id, decoration)
end

function spvc_compiler_has_member_decoration(compiler, id, member_index, decoration)
    ccall((:spvc_compiler_has_member_decoration, libspirv_cross), spvc_bool, (spvc_compiler, spvc_type_id, UInt32, SpvDecoration), compiler, id, member_index, decoration)
end

function spvc_compiler_get_name(compiler, id)
    ccall((:spvc_compiler_get_name, libspirv_cross), Cstring, (spvc_compiler, SpvId), compiler, id)
end

function spvc_compiler_get_decoration(compiler, id, decoration)
    ccall((:spvc_compiler_get_decoration, libspirv_cross), UInt32, (spvc_compiler, SpvId, SpvDecoration), compiler, id, decoration)
end

function spvc_compiler_get_decoration_string(compiler, id, decoration)
    ccall((:spvc_compiler_get_decoration_string, libspirv_cross), Cstring, (spvc_compiler, SpvId, SpvDecoration), compiler, id, decoration)
end

function spvc_compiler_get_member_decoration(compiler, id, member_index, decoration)
    ccall((:spvc_compiler_get_member_decoration, libspirv_cross), UInt32, (spvc_compiler, spvc_type_id, UInt32, SpvDecoration), compiler, id, member_index, decoration)
end

function spvc_compiler_get_member_decoration_string(compiler, id, member_index, decoration)
    ccall((:spvc_compiler_get_member_decoration_string, libspirv_cross), Cstring, (spvc_compiler, spvc_type_id, UInt32, SpvDecoration), compiler, id, member_index, decoration)
end

function spvc_compiler_get_member_name(compiler, id, member_index)
    ccall((:spvc_compiler_get_member_name, libspirv_cross), Cstring, (spvc_compiler, spvc_type_id, UInt32), compiler, id, member_index)
end

function spvc_compiler_get_entry_points(compiler, entry_points, num_entry_points)
    ccall((:spvc_compiler_get_entry_points, libspirv_cross), spvc_result, (spvc_compiler, Ptr{Ptr{spvc_entry_point}}, Ptr{Cint}), compiler, entry_points, num_entry_points)
end

function spvc_compiler_set_entry_point(compiler, name, model)
    ccall((:spvc_compiler_set_entry_point, libspirv_cross), spvc_result, (spvc_compiler, Cstring, SpvExecutionModel), compiler, name, model)
end

function spvc_compiler_rename_entry_point(compiler, old_name, new_name, model)
    ccall((:spvc_compiler_rename_entry_point, libspirv_cross), spvc_result, (spvc_compiler, Cstring, Cstring, SpvExecutionModel), compiler, old_name, new_name, model)
end

function spvc_compiler_get_cleansed_entry_point_name(compiler, name, model)
    ccall((:spvc_compiler_get_cleansed_entry_point_name, libspirv_cross), Cstring, (spvc_compiler, Cstring, SpvExecutionModel), compiler, name, model)
end

function spvc_compiler_set_execution_mode(compiler, mode)
    ccall((:spvc_compiler_set_execution_mode, libspirv_cross), Cvoid, (spvc_compiler, SpvExecutionMode), compiler, mode)
end

function spvc_compiler_unset_execution_mode(compiler, mode)
    ccall((:spvc_compiler_unset_execution_mode, libspirv_cross), Cvoid, (spvc_compiler, SpvExecutionMode), compiler, mode)
end

function spvc_compiler_set_execution_mode_with_arguments(compiler, mode, arg0, arg1, arg2)
    ccall((:spvc_compiler_set_execution_mode_with_arguments, libspirv_cross), Cvoid, (spvc_compiler, SpvExecutionMode, UInt32, UInt32, UInt32), compiler, mode, arg0, arg1, arg2)
end

function spvc_compiler_get_execution_modes(compiler, modes, num_modes)
    ccall((:spvc_compiler_get_execution_modes, libspirv_cross), spvc_result, (spvc_compiler, Ptr{Ptr{SpvExecutionMode}}, Ptr{Cint}), compiler, modes, num_modes)
end

function spvc_compiler_get_execution_mode_argument(compiler, mode)
    ccall((:spvc_compiler_get_execution_mode_argument, libspirv_cross), UInt32, (spvc_compiler, SpvExecutionMode), compiler, mode)
end

function spvc_compiler_get_execution_mode_argument_by_index(compiler, mode, index)
    ccall((:spvc_compiler_get_execution_mode_argument_by_index, libspirv_cross), UInt32, (spvc_compiler, SpvExecutionMode, UInt32), compiler, mode, index)
end

function spvc_compiler_get_execution_model(compiler)
    ccall((:spvc_compiler_get_execution_model, libspirv_cross), SpvExecutionModel, (spvc_compiler,), compiler)
end

function spvc_compiler_get_type_handle(compiler, id)
    ccall((:spvc_compiler_get_type_handle, libspirv_cross), spvc_type, (spvc_compiler, spvc_type_id), compiler, id)
end

function spvc_type_get_base_type_id(type)
    ccall((:spvc_type_get_base_type_id, libspirv_cross), spvc_type_id, (spvc_type,), type)
end

function spvc_type_get_basetype(type)
    ccall((:spvc_type_get_basetype, libspirv_cross), spvc_basetype, (spvc_type,), type)
end

function spvc_type_get_bit_width(type)
    ccall((:spvc_type_get_bit_width, libspirv_cross), UInt32, (spvc_type,), type)
end

function spvc_type_get_vector_size(type)
    ccall((:spvc_type_get_vector_size, libspirv_cross), UInt32, (spvc_type,), type)
end

function spvc_type_get_columns(type)
    ccall((:spvc_type_get_columns, libspirv_cross), UInt32, (spvc_type,), type)
end

function spvc_type_get_num_array_dimensions(type)
    ccall((:spvc_type_get_num_array_dimensions, libspirv_cross), UInt32, (spvc_type,), type)
end

function spvc_type_array_dimension_is_literal(type, dimension)
    ccall((:spvc_type_array_dimension_is_literal, libspirv_cross), spvc_bool, (spvc_type, UInt32), type, dimension)
end

function spvc_type_get_array_dimension(type, dimension)
    ccall((:spvc_type_get_array_dimension, libspirv_cross), SpvId, (spvc_type, UInt32), type, dimension)
end

function spvc_type_get_num_member_types(type)
    ccall((:spvc_type_get_num_member_types, libspirv_cross), UInt32, (spvc_type,), type)
end

function spvc_type_get_member_type(type, index)
    ccall((:spvc_type_get_member_type, libspirv_cross), spvc_type_id, (spvc_type, UInt32), type, index)
end

function spvc_type_get_storage_class(type)
    ccall((:spvc_type_get_storage_class, libspirv_cross), SpvStorageClass, (spvc_type,), type)
end

function spvc_type_get_image_sampled_type(type)
    ccall((:spvc_type_get_image_sampled_type, libspirv_cross), spvc_type_id, (spvc_type,), type)
end

function spvc_type_get_image_dimension(type)
    ccall((:spvc_type_get_image_dimension, libspirv_cross), SpvDim, (spvc_type,), type)
end

function spvc_type_get_image_is_depth(type)
    ccall((:spvc_type_get_image_is_depth, libspirv_cross), spvc_bool, (spvc_type,), type)
end

function spvc_type_get_image_arrayed(type)
    ccall((:spvc_type_get_image_arrayed, libspirv_cross), spvc_bool, (spvc_type,), type)
end

function spvc_type_get_image_multisampled(type)
    ccall((:spvc_type_get_image_multisampled, libspirv_cross), spvc_bool, (spvc_type,), type)
end

function spvc_type_get_image_is_storage(type)
    ccall((:spvc_type_get_image_is_storage, libspirv_cross), spvc_bool, (spvc_type,), type)
end

function spvc_type_get_image_storage_format(type)
    ccall((:spvc_type_get_image_storage_format, libspirv_cross), SpvImageFormat, (spvc_type,), type)
end

function spvc_type_get_image_access_qualifier(type)
    ccall((:spvc_type_get_image_access_qualifier, libspirv_cross), SpvAccessQualifier, (spvc_type,), type)
end

function spvc_compiler_get_declared_struct_size(compiler, struct_type, size)
    ccall((:spvc_compiler_get_declared_struct_size, libspirv_cross), spvc_result, (spvc_compiler, spvc_type, Ptr{Cint}), compiler, struct_type, size)
end

function spvc_compiler_get_declared_struct_size_runtime_array(compiler, struct_type, array_size, size)
    ccall((:spvc_compiler_get_declared_struct_size_runtime_array, libspirv_cross), spvc_result, (spvc_compiler, spvc_type, Cint, Ptr{Cint}), compiler, struct_type, array_size, size)
end

function spvc_compiler_get_declared_struct_member_size(compiler, type, index, size)
    ccall((:spvc_compiler_get_declared_struct_member_size, libspirv_cross), spvc_result, (spvc_compiler, spvc_type, UInt32, Ptr{Cint}), compiler, type, index, size)
end

function spvc_compiler_type_struct_member_offset(compiler, type, index, offset)
    ccall((:spvc_compiler_type_struct_member_offset, libspirv_cross), spvc_result, (spvc_compiler, spvc_type, UInt32, Ptr{UInt32}), compiler, type, index, offset)
end

function spvc_compiler_type_struct_member_array_stride(compiler, type, index, stride)
    ccall((:spvc_compiler_type_struct_member_array_stride, libspirv_cross), spvc_result, (spvc_compiler, spvc_type, UInt32, Ptr{UInt32}), compiler, type, index, stride)
end

function spvc_compiler_type_struct_member_matrix_stride(compiler, type, index, stride)
    ccall((:spvc_compiler_type_struct_member_matrix_stride, libspirv_cross), spvc_result, (spvc_compiler, spvc_type, UInt32, Ptr{UInt32}), compiler, type, index, stride)
end

function spvc_compiler_build_dummy_sampler_for_combined_images(compiler, id)
    ccall((:spvc_compiler_build_dummy_sampler_for_combined_images, libspirv_cross), spvc_result, (spvc_compiler, Ptr{spvc_variable_id}), compiler, id)
end

function spvc_compiler_build_combined_image_samplers(compiler)
    ccall((:spvc_compiler_build_combined_image_samplers, libspirv_cross), spvc_result, (spvc_compiler,), compiler)
end

function spvc_compiler_get_combined_image_samplers(compiler, samplers, num_samplers)
    ccall((:spvc_compiler_get_combined_image_samplers, libspirv_cross), spvc_result, (spvc_compiler, Ptr{Ptr{spvc_combined_image_sampler}}, Ptr{Cint}), compiler, samplers, num_samplers)
end

function spvc_compiler_get_specialization_constants(compiler, constants, num_constants)
    ccall((:spvc_compiler_get_specialization_constants, libspirv_cross), spvc_result, (spvc_compiler, Ptr{Ptr{spvc_specialization_constant}}, Ptr{Cint}), compiler, constants, num_constants)
end

function spvc_compiler_get_constant_handle(compiler, id)
    ccall((:spvc_compiler_get_constant_handle, libspirv_cross), spvc_constant, (spvc_compiler, spvc_constant_id), compiler, id)
end

function spvc_compiler_get_work_group_size_specialization_constants(compiler, x, y, z)
    ccall((:spvc_compiler_get_work_group_size_specialization_constants, libspirv_cross), spvc_constant_id, (spvc_compiler, Ptr{spvc_specialization_constant}, Ptr{spvc_specialization_constant}, Ptr{spvc_specialization_constant}), compiler, x, y, z)
end

function spvc_compiler_get_active_buffer_ranges(compiler, id, ranges, num_ranges)
    ccall((:spvc_compiler_get_active_buffer_ranges, libspirv_cross), spvc_result, (spvc_compiler, spvc_variable_id, Ptr{Ptr{spvc_buffer_range}}, Ptr{Cint}), compiler, id, ranges, num_ranges)
end

function spvc_constant_get_scalar_fp16(constant, column, row)
    ccall((:spvc_constant_get_scalar_fp16, libspirv_cross), Cfloat, (spvc_constant, UInt32, UInt32), constant, column, row)
end

function spvc_constant_get_scalar_fp32(constant, column, row)
    ccall((:spvc_constant_get_scalar_fp32, libspirv_cross), Cfloat, (spvc_constant, UInt32, UInt32), constant, column, row)
end

function spvc_constant_get_scalar_fp64(constant, column, row)
    ccall((:spvc_constant_get_scalar_fp64, libspirv_cross), Cdouble, (spvc_constant, UInt32, UInt32), constant, column, row)
end

function spvc_constant_get_scalar_u32(constant, column, row)
    ccall((:spvc_constant_get_scalar_u32, libspirv_cross), UInt32, (spvc_constant, UInt32, UInt32), constant, column, row)
end

function spvc_constant_get_scalar_i32(constant, column, row)
    ccall((:spvc_constant_get_scalar_i32, libspirv_cross), Cint, (spvc_constant, UInt32, UInt32), constant, column, row)
end

function spvc_constant_get_scalar_u16(constant, column, row)
    ccall((:spvc_constant_get_scalar_u16, libspirv_cross), UInt32, (spvc_constant, UInt32, UInt32), constant, column, row)
end

function spvc_constant_get_scalar_i16(constant, column, row)
    ccall((:spvc_constant_get_scalar_i16, libspirv_cross), Cint, (spvc_constant, UInt32, UInt32), constant, column, row)
end

function spvc_constant_get_scalar_u8(constant, column, row)
    ccall((:spvc_constant_get_scalar_u8, libspirv_cross), UInt32, (spvc_constant, UInt32, UInt32), constant, column, row)
end

function spvc_constant_get_scalar_i8(constant, column, row)
    ccall((:spvc_constant_get_scalar_i8, libspirv_cross), Cint, (spvc_constant, UInt32, UInt32), constant, column, row)
end

function spvc_constant_get_subconstants(constant, constituents, count)
    ccall((:spvc_constant_get_subconstants, libspirv_cross), Cvoid, (spvc_constant, Ptr{Ptr{spvc_constant_id}}, Ptr{Cint}), constant, constituents, count)
end

function spvc_constant_get_type(constant)
    ccall((:spvc_constant_get_type, libspirv_cross), spvc_type_id, (spvc_constant,), constant)
end

function spvc_compiler_get_binary_offset_for_decoration(compiler, id, decoration, word_offset)
    ccall((:spvc_compiler_get_binary_offset_for_decoration, libspirv_cross), spvc_bool, (spvc_compiler, spvc_variable_id, SpvDecoration, Ptr{UInt32}), compiler, id, decoration, word_offset)
end

function spvc_compiler_buffer_is_hlsl_counter_buffer(compiler, id)
    ccall((:spvc_compiler_buffer_is_hlsl_counter_buffer, libspirv_cross), spvc_bool, (spvc_compiler, spvc_variable_id), compiler, id)
end

function spvc_compiler_buffer_get_hlsl_counter_buffer(compiler, id, counter_id)
    ccall((:spvc_compiler_buffer_get_hlsl_counter_buffer, libspirv_cross), spvc_bool, (spvc_compiler, spvc_variable_id, Ptr{spvc_variable_id}), compiler, id, counter_id)
end

function spvc_compiler_get_declared_capabilities(compiler, capabilities, num_capabilities)
    ccall((:spvc_compiler_get_declared_capabilities, libspirv_cross), spvc_result, (spvc_compiler, Ptr{Ptr{SpvCapability}}, Ptr{Cint}), compiler, capabilities, num_capabilities)
end

function spvc_compiler_get_declared_extensions(compiler, extensions, num_extensions)
    ccall((:spvc_compiler_get_declared_extensions, libspirv_cross), spvc_result, (spvc_compiler, Ptr{Ptr{Cstring}}, Ptr{Cint}), compiler, extensions, num_extensions)
end

function spvc_compiler_get_remapped_declared_block_name(compiler, id)
    ccall((:spvc_compiler_get_remapped_declared_block_name, libspirv_cross), Cstring, (spvc_compiler, spvc_variable_id), compiler, id)
end

function spvc_compiler_get_buffer_block_decorations(compiler, id, decorations, num_decorations)
    ccall((:spvc_compiler_get_buffer_block_decorations, libspirv_cross), spvc_result, (spvc_compiler, spvc_variable_id, Ptr{Ptr{SpvDecoration}}, Ptr{Cint}), compiler, id, decorations, num_decorations)
end

ext: jl


text:
#--------------------------------------------------------------
# TensorTrains package
#--------------------------------------------------------------
# Basic functionality of the tensor train format
module TensorTrains

using LinearAlgebra, Logging

const Tensor3 = Array{Float64,3}
const Tensor4 = Array{Float64,4}

import Base.reverse, Base.transpose, LinearAlgebra.norm, LinearAlgebra.dot

export Tensor, TensorMatrix, Tensor3, Tensor4,
      dimpermute,
      ttones, tteye, ttdelta, ttrandn, ttdiagm,
      sizes, ranks, decompress, revdecompress, unfold, evaluate,
      scale, scale!, hadamard, add!, add,
      dot, norm!, norm, rightOrth!, leftOrth!,
      svd!, svdtrunc!, softthresh!,
      matvec, matvec!, matmat,
      reverse, transpose, ttkron

# tensor train representation
const Tensor = Vector{Tensor3}

# representations of matrices acting on tensor trains
const TensorMatrix = Vector{Tensor4}

module Condition
# Auxiliary routines for evaluating representation condition numbers
# (or respective bounds for matrix representations)
include("Condition.jl")
end  # module Condition

module Solvers
# Simple solvers
include("Solvers.jl")
end  # module Solvers

module TTPy
# Interface to python TT toolbox
# To install the required python packages via julia:
#
# using Conda, PyCall
# Conda.add("numpy")
# Conda.add("scipy")
# Conda.add("cython")
# Conda.add("pip")
# run(`$(PyCall.pyprogramname) -m pip install --user ttpy`)
include("TTPy.jl")
end  # module TTPy

# construct tensor train with mode sizes n and ranks 0
function Tensor(n::Vector{Int64})
  m = Vector{Tensor3}(undef, length(n))
  m[1] = Tensor3(1, n[1], 0)
  for i = 2:(length(m)-1)
    m[i] = Tensor3(0, n[i], 0)
  end
  m[end] = Tensor3(0, n[end], 1)
  return m
end

# construct uninitialized tensor train with mode sizes n and ranks r
function Tensor(n::Vector{Int64}, r::Vector{Int64})
  m = Vector{Tensor3}(undef, length(n))
  m[1] = Tensor3(1, n[1], r[1])
  for i in 2:(length(m)-1)
    m[i] = Tensor3(r[i-1], n[i], r[i])
  end
  m[end] = Tensor3(r[end], n[end], 1)
  return m
end

# construct uninitialized tensor train with order d
function Tensor(d::Int64)
  m = Vector{Tensor3}(undef, d)
  return m
end

# construct tensor train of order d with mode sizes n
function Tensor(n::Int64, d::Int64)
  m = Tensor3[Tensor3(1,n,0), fill(Tensor3(0,n,0), d-2)..., Tensor3(0,n,1)]
  return m
end

# construct uninitialized tensor train with order d, mode sizes n, ranks r
function Tensor(n::Int64, d::Int64, r::Int64)
  m = Vector{Tensor3}(undef, d)
  m[1] = Tensor3(1, n, r)
  for i in 2:(length(m)-1)
    m[i] = Tensor3(r, n, r)
  end
  m[d] = Tensor3(r, n, 1)
  return m
end

function Tensor3(n1::Int64, n2::Int64, n3::Int64)
  return Tensor3(undef, n1,n2,n3)
end

function Tensor4(n1::Int64, n2::Int64, n3::Int64, n4::Int64)
  return Tensor4(undef, n1,n2,n3,n4)
end

function TensorMatrix(n::Vector{Tuple{Int64,Int64}})
  m = Vector{Tensor4}(undef, length(n))
  m[1] = Tensor4(1, n[1]..., 0)
  for i = 2:(length(m)-1)
    m[i] = Tensor4(0, n[i]..., 0)
  end
  m[end] = Tensor4(0, n[end]..., 1)
  return m
end

function TensorMatrix(n::Vector{Tuple{Int64,Int64}}, r::Vector{Int64})
  m = Vector{Tensor4}(undef, length(n))
  m[1] = Tensor4(1, n[1]..., r[1])
  for i in 2:(length(m)-1)
    m[i] = Tensor4(r[i-1], n[i]..., r[i])
  end
  m[end] = Tensor4(r[end], n[end]..., 1)
  return m
end

function TensorMatrix(d::Int64)
  m = Vector{Tensor4}(undef, d)
  return m
end

function TensorMatrix(n::Int64, d::Int64)
  m = Tensor4[Tensor4(1,n,n,0), fill(Tensor4(0,n,n,0), d-2)..., Tensor4(0,n,n,1)]
  return m
end

function TensorMatrix(n::Int64, d::Int64, r::Int64)
  m = Vector{Tensor4}(undef, d)
  m[1] = Tensor4(1, n, n, r)
  for i in 2:(length(m)-1)
    m[i] = Tensor4(r, n, n, r)
  end
  m[d] = Tensor4(r, n, n, 1)
  return m
end

# convert tensor matrix to regular tensor
function Tensor(t::TensorMatrix)
  r = Tensor(undef, length(t))
  for i in eachindex(t)
    r[i] = reshape(t[i], (size(t[i],1),
          size(t[i],2)*size(t[i],3),size(t[i],4)))
  end
  return r
end

# convert tensor to tensor matrix with specified mode sizes
function TensorMatrix(t::Tensor, n::Vector{Tuple{Int64,Int64}})
  r = TensorMatrix(undef, length(t))
  for i in eachindex(t)
    r[i] = reshape(t[i], (size(t[i],1), n[i]..., size(t[i],3)))
  end
  return r
end

# returns vector of mode sizes
function sizes(t::Tensor)
  return [size(t[i],2) for i in eachindex(t)]
end

# returns vector of mode size tuples
function sizes(t::TensorMatrix)
  return [(size(t[i],2),size(t[i],3)) for i in eachindex(t)]
end

# returns row or column mode size vectors
function sizes(t::TensorMatrix, i::Int64)
  if i == 1
    return [ size(t[i],2) for i in eachindex(t)]
  elseif i == 2
    return [ size(t[i],3) for i in eachindex(t)]
  end
end

# substitute for permutedims
function dimpermute(T::Array, perm::Vector)
  d = ndims(T)
  Tp = zeros(size(T)[perm])
  idx = zeros(Int64, d)
  for i in 1:length(T)
    rem = i - 1
    for k = 1:d
      b = size(T,k)
      r = fld(rem, b)
      idx[k] = rem - b*r
      rem = r
    end
    ip = idx[perm[end]]
    for k = (d-1):-1:1
      b = size(T,perm[k])
      ip = ip*b + idx[perm[k]]
    end
    ip += 1
    Tp[ip] = T[i]
  end
  return Tp
end

# tensor with 1 at index E and 0 otherwise
function ttdelta(n::Int64, E::Vector{Int64})
  t = Tensor(2, length(E), 1)
  for i in eachindex(t)
    fill!(t[i], 0.)
    t[i][1,E[i],1] = 1.
  end
  return t
end

# tensor with all entries 1
function ttones(L::Int64)
  t = Tensor(2, L, 1)
  for i in eachindex(t)
    t[i][1,:,1] = [1., 1.]
  end
  return t
end

# tensor with normally distributed random cores
function ttrandn(L::Int64)
  t = Tensor(2, L, 1)
  for i in eachindex(t)
    t[i][1,:,1] = randn(2)
  end
  return t
end

# identity matrix
function tteye(L::Int64)
  t = TensorMatrix(2, L, 1)
  for i in eachindex(t)
    t[i][1,:,:,1] = Matrix{Float64}(I, 2, 2)
  end
  return t
end

# convert tensor to diagonal matrix
function ttdiagm(t::Tensor)
  L = length(t)
  T = TensorMatrix([(n,n) for n in sizes(t)], ranks(t))
  for k = 1:L
    for i = 1:size(T[k],1)
      for j = 1:size(T[k],4)
        T[k][i,:,:,j] = diagm(0=>t[k][i,:,j])
      end
    end
  end
  return T
end

# turn tt representation into full array
function decompress(t::Tensor)
  s1 = size(t[1])
  r = reshape(t[1], (s1[1]*s1[2], s1[3]))
  for i = 2:length(t)
    r1 = size(r, 1)
    si = size(t[i])
    r = reshape(r*reshape(t[i], (si[1], si[2]*si[3])), (r1*si[2], si[3]))
  end
  return reshape(r, (sizes(t)...,))
end

# turn tt representation into full array with reversed order of indices
function revdecompress(t::Tensor)
  tx = reverse(t)
  return decompress(tx)
end

# assemble tensor train with first and/or last ranks > 1
function unfold(t::Tensor)
  s1 = size(t[1])
  r = reshape(t[1], (s1[1]*s1[2], s1[3]))
  for i = 2:length(t)
    r1 = size(r, 1)
    si = size(t[i])
    r = reshape(r*reshape(t[i], (si[1], si[2]*si[3])), (r1*si[2], si[3]))
  end
  return reshape(r, (s1[1],prod(sizes(t)),size(t[end],3)))
end

# assemble matrix tensor train
function unfold(T::TensorMatrix)
  s = [sizes(T)[i][j] for j = 1:2, i in eachindex(sizes(T))]
  t = Tensor(T)
  L = length(t)
  dim = zeros(Int64,2L);
  dim[1:L] = (2L-1):-2:1;
  dim[(L+1):(2L)] = (2L):-2:2;
  return reshape(dimpermute(
    reshape(decompress(t), (s[:]...,)),
    dim), (prod(s[1,:]), prod(s[2,:])))
end

# evaluate one entry of tensor represented by tensor train
function evaluate(t::Tensor, idx)
  r = t[1][:,idx[1],:]
  for i = 2:length(t)
    r *= t[i][:,idx[i],:]
  end
  return r
end

# multiply tensor by scalar (in-place)
function scale!(t::Tensor, α::Float64)
  t[1] .*= α
  return t
end

# multiply tensor by scalar (creating copy)
function scale(t::Tensor, α::Float64)
  s = deepcopy(t)
  s[1] .*= α
  return s
end

# multiply tensor matrix by scalar (in-place)
function scale!(t::TensorMatrix, α::Float64)
  t[1] .*= α
  return t
end

# multiply tensor matrix by scalar (creating copy)
function scale(t::TensorMatrix, α::Float64)
  s = deepcopy(t)
  s[1] .*= α
  return s
end

# entry-wise product of tensor trains
function hadamard(t1::Tensor, t2::Tensor)
  r = Tensor(length(t1))
  for n in eachindex(t1)
    s1 = size(t1[n])
    s2 = size(t2[n])
    r[n] = Tensor3(s1[1]*s2[1], s1[2], s1[3]*s2[3])
    for i = 1:s1[1], j = 1:s2[1], k = 1:s1[3], l = 1:s2[3]
      r[n][i+s1[1]*(j-1),:,k+s1[3]*(l-1)] = t1[n][i,:,k].*t2[n][j,:,l]
    end
  end
  return r
end

# add two tensor trains (result in t)
function add!(t::Tensor, s::Tensor)
  return add!(t, 1., s)
end

# add α*s to t (with result in t)
function add!(t::Tensor, α::Float64, s::Tensor)
  for (i,(τ,σ)) in enumerate(zip(t,s))
    if i == 1
      c = Tensor3(1,size(τ,2),size(τ,3)+size(σ,3))
      c[:,:,1:size(τ,3)] = τ;
      c[:,:,(size(τ,3)+1):size(c,3)] = α*σ
    elseif i == length(t)
      c = zeros(size(τ,1)+size(σ,1),size(τ,2),1)
      c[1:size(τ,1),:,:] = τ;
      c[(size(τ,1)+1):size(c,1),:,:] = σ
    else
      c = zeros(size(τ,1)+size(σ,1),size(τ,2),size(τ,3)+size(σ,3))
      c[1:size(τ,1),:,1:size(τ,3)] = τ;
      c[(size(τ,1)+1):size(c,1),:,(size(τ,3)+1):size(c,3)] = σ
    end
    t[i] = c
  end
  return t
end

# compute αt + βs in tt format (returning new tensor)
function add(α::Float64, t::Tensor, β::Float64, s::Tensor)
  r = Tensor(sizes(t))
  for (i,(τ,σ)) in enumerate(zip(t,s))
    if i == 1
      c = Tensor3(1,size(τ,2),size(τ,3)+size(σ,3))
      c[:,:,1:size(τ,3)] = α*τ;
      c[:,:,(size(τ,3)+1):size(c,3)] = β*σ
    elseif i == length(t)
      c = zeros(size(τ,1)+size(σ,1),size(τ,2),1)
      c[1:size(τ,1),:,:] = τ;
      c[(size(τ,1)+1):size(c,1),:,:] = σ
    else
      c = zeros(size(τ,1)+size(σ,1),size(τ,2),size(τ,3)+size(σ,3))
      c[1:size(τ,1),:,1:size(τ,3)] = τ;
      c[(size(τ,1)+1):size(c,1),:,(size(τ,3)+1):size(c,3)] = σ
    end
    r[i] = c
  end
  return r
end

# add βs to tensor matrix t (with result in t)
function add!(t::TensorMatrix, β::Float64, s::TensorMatrix)
  for (i,(τ,σ)) in enumerate(zip(t,s))
    if i == 1
      c = Tensor4(1,size(τ,2),size(τ,3),size(τ,4)+size(σ,4))
      c[:,:,:,1:size(τ,4)] = τ;
      c[:,:,:,(size(τ,4)+1):size(c,4)] = β*σ
    elseif i == length(t)
      c = zeros(size(τ,1)+size(σ,1),size(τ,2),size(τ,3),1)
      c[1:size(τ,1),:,:,:] = τ;
      c[(size(τ,1)+1):size(c,1),:,:,:] = σ
    else
      c = zeros(size(τ,1)+size(σ,1),size(τ,2),size(τ,3),size(τ,4)+size(σ,4))
      c[1:size(τ,1),:,:,1:size(τ,4)] = τ;
      c[(size(τ,1)+1):size(c,1),:,:,(size(τ,4)+1):size(c,4)] = σ
    end
    t[i] = c
  end
  return t
end

# add tensor matrices t and s (result in t)
function add!(t::TensorMatrix, s::TensorMatrix)
  return add!(t, 1., s)
end

function add(α::Float64, t::TensorMatrix, β::Float64, s::TensorMatrix)
  r = TensorMatrix(sizes(t))
  for (i,(τ,σ)) in enumerate(zip(t,s))
    if i == 1
      c = Tensor4(1,size(τ,2),size(τ,3),size(τ,4)+size(σ,4))
      c[:,:,:,1:size(τ,4)] = α*τ;
      c[:,:,:,(size(τ,4)+1):size(c,4)] = β*σ
    elseif i == length(t)
      c = zeros(size(τ,1)+size(σ,1),size(τ,2),size(τ,3),1)
      c[1:size(τ,1),:,:,:] = τ;
      c[(size(τ,1)+1):size(c,1),:,:,:] = σ
    else
      c = zeros(size(τ,1)+size(σ,1),size(τ,2),size(τ,3),size(τ,4)+size(σ,4))
      c[1:size(τ,1),:,:,1:size(τ,4)] = τ;
      c[(size(τ,1)+1):size(c,1),:,:,(size(τ,4)+1):size(c,4)] = σ
    end
    r[i] = c
  end
  return r
end

# apply tensor matrix to tensor (result in new tensor)
function matvec(M::TensorMatrix, t::Tensor)
  r = Tensor(sizes(M, 1))
  matvec!(M, t, r)
  return r
end

# apply tensor matrix M to tensor t, with result in r
# (which needs to have correct order)
function matvec!(M::TensorMatrix, t::Tensor, r::Tensor)
  for m in 1:length(t)
    s = size(t[m])
    S = size(M[m])
    T = t[m]
    r[m] = Tensor3(s[1]*S[1], S[2], s[3]*S[4])
    for i = 1:s[1], j = 1:S[1], k = 1:s[3], l = 1:S[4]
      r[m][i+s[1]*(j-1),:,k+s[3]*(l-1)] = M[m][j,:,:,l]*T[i,:,k]
    end
  end
  for m in (length(t)+1):length(M)
    r[m] = reshape(M[m], (size(M[m],1), size(M[m],2), size(M[m],4)))
  end
  if size(r[end],2) == 1
    x = r[end-1]; y = r[end]
    z = reshape(reshape(x, (size(x,1)*size(x,2), size(x,3)))*reshape(y, size(y,1)),
          (size(x,1), size(x,2), 1) )
    r[end-1] = z
    pop!(r)
  end
end

# compute representation of S * T
function matmat(S::TensorMatrix, T::TensorMatrix)
  R = TensorMatrix(undef, length(S))
  for m in eachindex(S)
    R[m] = matmat(S[m], T[m])
  end
  return R
end

function matmat(S::Tensor4, T::Tensor4)
  σ = size(S)
  τ = size(T)
  R = Tensor4(σ[1]*τ[1], σ[2], τ[3], σ[4]*τ[4])
  for i = 1:τ[1], j = 1:σ[1], k = 1:τ[4], l = 1:σ[4]
    R[i+τ[1]*(j-1),:,:,k+τ[4]*(l-1)] = S[j,:,:,l] * T[i,:,:,k]
  end
  return R
end

# euclidean inner product of tt tensors
function dot(t::Tensor, s::Tensor)
  sizet = size(t[1])
  sizes = size(s[1])
  G = reshape(t[1], (sizet[1]*sizet[2],sizet[3]))' *
      reshape(s[1], (sizes[1]*sizes[2],sizes[3]))
  for i = 2:length(t)
    Gi = (t[i][:,1,:])' * G * s[i][:,1,:]
    for n = 2:size(t[i], 2)
      Gi += (t[i][:,n,:])' * G * s[i][:,n,:]
    end
    G = Gi
  end
  return G[1,1]
end

# rank vector of tt tensor (omitting first and last, assumed =1)
function ranks(t::Tensor)
  rvec = zeros(Int64, length(t)-1)
  for i = 1:(length(t)-1)
    rvec[i] = size(t[i], 3)
  end
  return rvec
end

# rank vector of tt matrix (omitting first and last, assumed =1)
function ranks(t::TensorMatrix)
  rvec = zeros(Int64, length(t)-1)
  for i = 1:(length(t)-1)
    rvec[i] = size(t[i], 4)
  end
  return rvec
end

# euclidean norm of tt tensor (t not modified)
function norm(t::Tensor)
  s = deepcopy(t)
  return norm!(s)
end

# euclidean norm of tt tensor (t is orthogonalized)
function norm!(t::Tensor)
  leftOrth!(t)
  return norm(t[1])
end

# orthogonalize to the right in core i
function rightOrth!(t::Tensor, i::Int64)
  s = size(t[i])
  F = qr(reshape(t[i],(s[1]*s[2],s[3])))
  Q = Matrix(F.Q)
  t[i] = reshape(Q, (s[1],s[2],size(Q,2)))
  s = size(t[i+1])
  t[i+1] = reshape(F.R*reshape(t[i+1],(s[1],s[2]*s[3])), (size(F.R,1),s[2],s[3]))
end

# orthogonalize entire tensor to the right
function rightOrth!(t::Tensor)
  for i = 1:(length(t)-1)
    rightOrth!(t, i)
  end
end

# orthogonalize to the left in core i
function leftOrth!(t::Tensor, i::Int64)
  s = size(t[i])
  F = qr(reshape(t[i],(s[1],s[2]*s[3]))')
  Q = Matrix(F.Q)
  t[i] = reshape(Q', (size(Q,2),s[2],s[3]))
  s = size(t[i-1])
  t[i-1] = reshape(reshape(t[i-1],(s[1]*s[2],s[3]))*F.R', (s[1],s[2],size(F.R,1)))
end

# orthogonalize entire tensor to the left
function leftOrth!(t::Tensor)
  for i = length(t):-1:2
    leftOrth!(t, i)
  end
end

# bring tensor into tt-svd form, return list of singular value vectors
# of matricizations
function svd!(t::Tensor)
  if minimum(ranks(t)) > 0
    rightOrth!(t)
    svlist = Array{Vector{Float64}}(undef, length(t)-1)
    for i = length(t):-1:2
      s = size(t[i])
      U, σ, V = svd(reshape(t[i], (s[1],s[2]*s[3])))
      svlist[i-1] = σ
      s = (size(V,2), s[2], s[3])
      t[i] = reshape(V', s)
      rmul!(U, Diagonal(σ))
      s = size(t[i-1])
      t[i-1] = reshape(reshape(t[i-1], (s[1]*s[2],s[3]))*U, (s[1],s[2],size(U,2)))
    end
  else
    svlist = fill(Float64[], length(t)-1)
    t .= Tensor(sizes(t))
  end
  return svlist
end

# truncate t given in tt-svd form according to singular values, up to
# euclidean norm error ε
# usage example for t of type Tensor:
#   svdtrunc!(t, svd!(t), 1e-6)
function svdtrunc!(t::Tensor, svlist, ɛ)
  rvec = [length(s) for s in svlist]
  svvec = Vector{Tuple{Int64,Int64,Float64}}()
  for i in eachindex(svlist)
    for (j,σ) in enumerate(svlist[i])
      push!(svvec, (i, j, σ))
    end
  end
  sort!(svvec, lt = (x,y)->(x[3]<y[3] || (x[3]==y[3] && x[2]>y[2])))
  s = 0.
  stopi = length(svvec)
  for (i,τ) in enumerate(svvec)
    if s + τ[3]^2 <= ɛ^2
      s += τ[3]^2
    else
      stopi = i - 1
      break
    end
  end
  if stopi > 0
    for τ in svvec[1:stopi]
      rvec[τ[1]] = τ[2]-1
    end
  end
  t[1] = t[1][:,:,1:rvec[1]]
  for i = 2:(length(t)-1)
    t[i] = t[i][1:rvec[i-1],:,1:rvec[i]]
  end
  t[end] = t[end][1:rvec[end],:,:]
  if minimum(ranks(t)) == 0
    t .= Tensor(sizes(t))
  end
  return t
end

# soft thresholding of singular values with threshold α
function softthresh!(t::Tensor, α::Float64)
  if minimum(ranks(t)) > 0
    rightOrth!(t)
    for i = length(t):-1:2
      s = size(t[i])
      U, σ, V = svd(reshape(t[i], (s[1],s[2]*s[3])))
      σ = σ[σ .> α] .- α
      r = length(σ)
      s = (r, s[2], s[3])
      t[i] = reshape(V[:,1:r]', s)
      U = U[:,1:r]
      rmul!(U, Diagonal(σ))
      s = size(t[i-1])
      t[i-1] = reshape(reshape(t[i-1], (s[1]*s[2],s[3]))*U, (s[1],s[2],size(U,2)))
    end
    if minimum(ranks(t)) == 0
      t .= Tensor(sizes(t))
    end
  else
    t .= Tensor(sizes(t))
  end
  return t
end

# reverse indices in tensor
function reverse(t::Tensor)
	L = length(t)
	s = Tensor(L)
	for l ∈ 1:L
		s[L-l+1] = permutedims(t[l], (3,2,1))
	end
	return s
end

# reverse indices in tensor matrix
function reverse(T::TensorMatrix)
	L = length(T)
	S = TensorMatrix(L)
	for l ∈ 1:L
		S[L-l+1] = permutedims(T[l], (4,2,3,1))
	end
	return S
end

# transpose of tensor matrix
function transpose(T::TensorMatrix)
	L = length(T)
	S = TensorMatrix(L)
	for l ∈ 1:L
		S[l] = permutedims(T[l], (1,3,2,4))
	end
	return S
end

transpose(t::Tensor4) = permutedims(t, [1,3,2,4])

function ttkron(A::Tensor4, B::Tensor4)
  return reshape(permutedims(reshape(A[:]*B[:]',
      (size(A)...,size(B)...)), [1,5,2,6,3,7,4,8]),
      ((size(A,i)*size(B,i) for i = 1:4)...,) )
end

function ttkron(A::TensorMatrix, B::TensorMatrix)
  C = TensorMatrix(length(A))
  for n = 1:length(C)
    C[n] = ttkron(A[n], B[n])
  end
  return C
end

function ttkron(A::Tensor3, B::Tensor3)
  return reshape(permutedims(reshape(A[:]*B[:]',
      (size(A)...,size(B)...)), [1,4,2,5,3,6]),
      ((size(A,i)*size(B,i) for i = 1:3)...,) )
end

function ttkron(A::Tensor, B::Tensor)
  C = Tensor(length(A))
  for n = 1:length(C)
    C[n] = ttkron(A[n], B[n])
  end
  return C
end

end # module


text:
function cauchymt(f::AbstractVector{T}, k::Int) where T
#
# generate a k-column Cauchy matrix of polynomial f
#
    n = length(f)
    m = n + k - 1
    A = zeros(T, m, k)
    for j = 1:k
        A[j:j+n-1,j] = f
    end
    A
end


text:
"""
This BellComm submodule provides methods for polyhedral analysis and transformations.
"""
module ConvexPolytope

using XPORTA, LinearAlgebra

include("./ConvexPolytope/adjacency_decomposition.jl")

# """
# facet_constraints(vertices, filename)
#
#     Uses xporta software to compute the facet inequalities given a set of
#     polytope vertices.
#
# Side-effect:
#     filename.poi and filename.poi.ieq files are saved to porta_tmp. porta_tmp
#     is optionally deleted at the end of function's execution.
#
# Inputs:
#     vertices: Array, contains polytope vertices.
#     filename: String (optional), name of .poi and .ieq files
#     dir: String (optional), directory where porta_tmp is created
#     cleanup: Boolean, if true tmp files are removed
#     suppress: Boolean, if true xporta does not print to STDOUT
#
# Outputs:
#     constraints: Tuple, contains arrays of dual vectors to vertices, (e.g. facets).
#       The constraint tupe is ordered in the following way
#           1: less than constraints <= 0
#           2: equal to constraints == 0
#           3: greater than constraints >= 0
# """
function facet_constraints(vertices; filename="facet_constraints_tmp", dir="./", cleanup=true, suppress=true)
    ieq = traf(POI( vertices=convert.(Int, vcat(map(v -> v', vertices)...))), filename=filename, dir=dir, cleanup=cleanup, verbose=!suppress)

    equalities = map(eq -> cat(-1*eq[end], eq[1:end-1]', dims=2) , eachrow(ieq.equalities))
    inequalities = map(ineq -> cat(-1*ineq[end], ineq[1:end-1]', dims=2) , eachrow(ieq.inequalities))

    (inequalities, equalities, [])
end

# """
# linear_interpolation(x1, x2; dw=0.1)
#
#     Creates an array of points which lie on the line connecting x1 to x2. Points
#     are interpolated by a step interval of dw where 0:dw:1.
#
# Inputs:
#     x1: Col Vector, point in convex polytope space
#     x2: Col Vector, point in convex polytope space
#     dw: float, 0<dw<1, sets the number of steps used in the interpolation
#
# Output:
#     interpolation: Array, contains points in convex polytope which scan across
#                    the line starting from x1 going to x2.
# """
function linear_interpolation(x1, x2; dw=0.1)
    ws = 0:dw:1
    dim = length(x1)
    interpolation = map(w -> reshape((1 - w)*x1 + w*x2,(dim, 1)), ws)

    interpolation
end

# """
# facet_contains(facet, x)
#
#     Returns true if the point x lies approximately on the facet.
#
# Inputs:
#     facet: Row Vector, contains hyperplane constant and normal vector of hyperplane.
#     x: Col Vector, a point the convex polytope space.
#
# Outputs:
#     Boolean, true if point x lies on the facet.
# """
function facet_contains(facet, x)
    isapprox(1 + facet_violation(facet, x), 1)
end

# """
# contains(facets, x):
#
#     Checks if the convex polytope contains the point x by checking each facet
#     for violation.
#
# Inputs:
#     facets: Array, contains facet inequalities
#     x: Col Vector, the point to test
#
# Output:
#     Boolean, true if the polytope contains point x
# """
function contains(facets, x)
    length(violated_facets(facets, x)) == 0
end

# """
# violated_facets(facets, x)
#
#     Returns the list of facets which are violated by the provided point, x. If
#     no facets are violated, an empty array is returned.
#
# Inputs:
#     facet: Row Vector, contains hyperplane constant and normal vector of hyperplane.
#     x: Col Vector, a point the convex polytope space.
#
# Outputs:
#     violations: Array, contains the facet inequalities which were violated by x.
# """
function violated_facets(facets, x)
    violations = filter(
        facet -> facet_violation(facet, x) > 0 && !facet_contains(facet, x),
        facets
    )

    violations
end

# """
# facets_contain(facets, x)
#
#     Returns true if point x satisfies all facet inequalities.
#
# Inputs:
#     facet: Row Vector, contains hyperplane constant and normal vector of hyperplane.
#     x: Col Vector, a point the convex polytope space.
#
# Outputs:
#     x_is_contained: Boolean, true if x lies within or along the polytope boundary.
# """
function facets_contain(facets, x)
    x_is_contained = all(map(
        facet -> facet_violation(facet, x) <= 0,
        facets
    ))

    x_is_contained
end

# """
# facet_violation(facet, x)
#
#     Computes the violation quantity of point x for a given facet. A violation > 0
#     indicates a point which does not satisfy the facet inequality.
#
# Inputs:
#     facet: Row Vector, contains hyperplane constant and normal vector of hyperplane.
#     x: Col Vector, a point the convex polytope space.
#
# Outputs:
#     violation: Float, Inequality is satisfied for violation <= 0.
# """
function facet_violation(facet, x)
    violation = NaN
    if length(x) == length(facet)
        violation = (facet*x)[1]
    elseif (length(x) + 1) == length(facet)
        violation = facet[1] + (facet[2:end]'*x)[1]
    else
        throw(ArgumentError("dimension mismatch: facet not 'dual' to coordinate"))
    end

    violation
end

# """
# facet_distance(facet, x)
#
#     Computes the (shortest) distance between point x and a given facet. A
#     distance < 0 indicates that the point lies within the polytope.
#
# Inputs:
#     facet: Row Vector, contains hyperplane constant and normal vector of hyperplane.
#     x: Col Vector, a point the convex polytope space.
#
# Outputs:
#     dist: Float, distance a point is from a bounding hyperplane.
# """
function facet_distance(facet, x)
    dist = ConvexPolytope.facet_violation(facet, x)/norm(facet[2:end])
    dist
end

# """
# distance(x1, x2)
#     Computes the distance between two points in behavior space.
#
# Inputs:
#     x1/2: Col Vector, represents a point in behavior space.
#
# Output:
#     dist: Integer, the Euclidean distance between the two point.
# """
function vertex_distance(x1, x2)
    dist = norm(x1 - x2)
    dist
end

# """
# collect_vertices(facets, vertices):
#
#     For each facet, find all vertices that lie on that facet and construct a
#     dictionary linking the facet to the set of vertices.
#
# Input:
#     facets: Array, contains valid facet inequalities
#     vertices: Array, contains valid vertices
#
# Output:
#     Array, contains dictionaries in the following form
#         {
#             "facet": [facet inequality],
#             "vertices": [[v1],[v2],...,[vN]]
#         }
# """
function collect_vertices(facets, vertices)
    fv_groups = map(
        (f)->Dict(
            "facet" => f,
            "vertices" => filter(v->ConvexPolytope.facet_contains(f, v), vertices)
        ), facets)

    fv_groups
end

# """
# collect_facets(vertices, facets):
#
#     For each vertex, find all facets which that vertex lies on and construct
#     a dictionary linking each vertex to all of its facets.
#
# Input:
#     vertices: Array, contains valid vertices
#     facets: Array, contains valid facet inequalities
#
# Output:
#     Array, contains dictionaries of the following form
#         {
#             "vertex": [vertex]
#             "facets": [[f1],[f2],...,[fN]]
#         }
# """
function collect_facets(vertices, facets)
    vf_groups = map(
        (v)->Dict(
            "vertex" => v,
            "facets" => filter(f->ConvexPolytope.facet_contains(f, v), facets)
        ), vertices)

    vf_groups
end

end


text:
module PerronFrobenius
    import CausalityToolsBase
    export RectangularBinning

    include("api.jl")
    include("boundary_conditions.jl")
    include("transfer_operators.jl")
    include("invariant_measures.jl")
    include("gridestimators/SingleGrid.jl")
    include("invariant_measures_grid.jl")

    using Requires
    function __init__()
        @require Simplices="d5428e67-3037-59ba-9ab1-57a04f0a3b6a" begin
            include("triangulationestimators/TriangulationEstimators.jl")
        end

        @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
            #include("v1/plot_recipes/invariant_measures.jl")
        end
    end
    #import .TriangulationEstimators: SimplexPoint, SimplexExact
    export SimplexPoint, SimplexExact
end


text:
import Base: +, -, *, /
import Base.Broadcast: broadcasted

export add, mul, neg, sub, pow, div

"""
    Add <: Func
"""
@func mutable struct Add end

forward(f::Add, A, B) = A .+ B

function backward(f::Add, gy)
    gx1, gx2 = gy, gy
    x_shape1, x_shape2 = size.(f.args)
    if x_shape1 != x_shape2
        gx1 = sumto(gx1, x_shape1)
        gx2 = sumto(gx2, x_shape2)
    end
    return gx1, gx2
end

add(A, B) = Add()(A, B)

+(A::Var, B::Var) = add(A, B)
+(A::Var, B) = add(A, B)
+(A, B::Var) = add(A, B)

broadcasted(::typeof(+), A::Var, B::Var) = add(A, B)
broadcasted(::typeof(+), A::Var, B) = add(A, B)
broadcasted(::typeof(+), A, B::Var) = add(A, B)

"""
    Mul <: Func
"""
@func mutable struct Mul end

forward(f::Mul, A, B) = A .* B

function backward(f::Mul, gy)
    x1, x2 = f.args
    gx1 = gy .* x2
    gx2 = gy .* x1
    x_shape1, x_shape2 = size.(f.args)
    if x_shape1 != x_shape2
        gx1 = sumto(gx1, x_shape1)
        gx2 = sumto(gx2, x_shape2)
    end
    return gx1, gx2
end

mul(A, B) = Mul()(A, B)

broadcasted(::typeof(*), A::Var, B::Var) = mul(A, B)
broadcasted(::typeof(*), A::Var, B) = mul(A, B)
broadcasted(::typeof(*), A, B::Var) = mul(A, B)

"""
    Neg <: Func
"""
@func mutable struct Neg end

forward(f::Neg, x) = -x

backward(f::Neg, gy) = -gy

neg(x) = Neg()(x)

-(x::Var) = neg(x)

"""
    Sub <: Func
"""
@func mutable struct Sub end

forward(f::Sub, A, B) = A .- B

function backward(f::Sub, gy)
    gx1, gx2 = gy, -gy
    x_shape1, x_shape2 = size.(f.args)
    if x_shape1 != x_shape2
        gx1 = sumto(gx1, x_shape1)
        gx2 = sumto(gx2, x_shape2)
    end
    return gx1, gx2
end

sub(A, B) = Sub()(A, B)

-(A::Var, B::Var) = sub(A, B)
-(A::Var, B) = sub(A, B)
-(A, B::Var) = sub(A, B)

broadcasted(::typeof(-), A::Var, B::Var) = sub(A, B)
broadcasted(::typeof(-), A::Var, B) = sub(A, B)
broadcasted(::typeof(-), A, B::Var) = sub(A, B)

"""
    Pow <: Func
"""
@func mutable struct Pow{T}
    c::T
    Pow(c::T) where T = new{T}(c)
end

forward(f::Pow, x) = x.^f.c

backward(f::Pow, gy) = f.c .* f.args[1]^(f.c - 1) .* gy

pow(x, c) = Pow(c)(x)

Base.:^(x::Var, c) = pow(x, c)

"""
    Div <: Func
"""
@func mutable struct Div end

forward(f::Div, A, B) = A ./ B

function backward(f::Div, gy)
    x1, x2 = f.args
    gx1 = gy / x2
    gx2 = gy .* (-x1 / x2^2)
    x_shape1, x_shape2 = size.(f.args)
    if x_shape1 != x_shape2
        gx1 = sumto(gx1, x_shape1)
        gx2 = sumto(gx2, x_shape2)
    end
    return gx1, gx2
end

div(A, B) = Div()(A, B)

/(A::Var, B::Var) = div(A, B)
/(A::Var, B) = div(A, B)
/(A, B::Var) = div(A, B)

broadcasted(::typeof(/), A::Var, B::Var) = div(A, B)
broadcasted(::typeof(/), A::Var, B) = div(A, B)
broadcasted(::typeof(/), A, B::Var) = div(A, B)


text:
@testset "TreeView" begin
    for node in (Node{Node{LeafNode{Float64, 2, 1}, 2, 1}, 2, 1}(),
                 HashNode{Node{LeafNode{Float64, 2, 1}, 2, 1}, 2, 1}(),
                 DynamicNode{Node{LeafNode{Float64, 2, 1}, 2, 1}, 2}(2,2),
                 DynamicHashNode{Node{LeafNode{Float64, 2, 1}, 2, 1}, 2}(2,2))
        A = TreeView(node)

        @test size(A) == (8, 8)

        A[1:3, 1:3] .= reshape(1:9, 3, 3)
        @test A[1:3, 1:3] == reshape(1:9, 3, 3)

        A[7:8, 7:8] .= reshape(11:14, 2, 2)
        @test A[7:8, 7:8] == reshape(11:14, 2, 2)

        # check activation
        mask = map(eachindex(A)) do i
            isactive(A, Tuple(i)...)
        end
        @test all(mask[1:3, 1:3])
        @test all(mask[7:8, 7:8])
        @test !any(mask[1:3, 4:6])
        @test !any(mask[7:8, 1:3])

        # nleaves
        @test nleaves(A) == 3*3 + 2*2

        # deactivating
        for i in CartesianIndices((1:2, 1:2))
            A[i] = nothing
            @test isactive(A, Tuple(i)...) == false
        end

        # cleanup!
        TreeArrays.cleanup!(rootnode(A))
        @test isallocated(rootnode(A), 1) == true
        @test isallocated(TreeArrays.unsafe_getindex(rootnode(A), 1), 1) == false

        # check allocation after deactivating
        A .= nothing
        @test isallocated(rootnode(A), 1) == true
        @test isallocated(rootnode(A), 2) == false
        @test isallocated(rootnode(A), 3) == false
        @test isallocated(rootnode(A), 4) == true
        @test isallocated(TreeArrays.unsafe_getindex(rootnode(A), 1), 1) == false
        @test isallocated(TreeArrays.unsafe_getindex(rootnode(A), 1), 2) == true
        @test isallocated(TreeArrays.unsafe_getindex(rootnode(A), 1), 3) == true
        @test isallocated(TreeArrays.unsafe_getindex(rootnode(A), 1), 4) == true
        @test isallocated(TreeArrays.unsafe_getindex(rootnode(A), 4), 1) == false
        @test isallocated(TreeArrays.unsafe_getindex(rootnode(A), 4), 2) == false
        @test isallocated(TreeArrays.unsafe_getindex(rootnode(A), 4), 3) == false
        @test isallocated(TreeArrays.unsafe_getindex(rootnode(A), 4), 4) == true

        # cleanup!
        TreeArrays.cleanup!(rootnode(A))
        @test isallocated(rootnode(A), 1) == false
        @test isallocated(rootnode(A), 2) == false
        @test isallocated(rootnode(A), 3) == false
        @test isallocated(rootnode(A), 4) == false
    end
    @testset "threads computations" begin
        for node in (Node{Node{LeafNode{Float64, 2, 3}, 2, 3}, 2, 3}(),
                     HashNode{Node{LeafNode{Float64, 2, 3}, 2, 3}, 2, 3}(),
                     DynamicNode{Node{LeafNode{Float64, 2, 3}, 2, 3}, 2}(8,8))
            n = 2^3 * 2^3 * 2^3
            A = TreeView(node)
            @test size(A) == (n, n)

            A[1:100, 1:100] .= 1
            A[200:end, 400:end] .= 1
            A[400:end, 500:end] .= nothing
            @test nleaves(A) == 100*100 + length(200:n)*length(400:n) - length(400:n)*length(500:n)

            A .= nothing
            @test nleaves(A) == 0
        end
    end
end


text:
module ShaleDrillingLikelihood_StateSpaceTest

using ShaleDrillingLikelihood
using Test

using ShaleDrillingLikelihood: state_space_vector,
    actionspace,
    state,
    sprime,
    state_if_never_drilled,
    _nstates_per_D,
    end_lrn,
    end_ex0,
    end_ex1,
    _dmax,
    _τrem,
    inf_fm_lrn,
    exploratory_learning,
    ind_lrn,
    maxlease,
    state_space_vector,
    _regime,
    _τ1,
    _τ0,
    _D ,
    _d1,
    ind_exp, ind_inf, end_inf,ind_ex1,ind_ex0,
    state_idx,
    _horizon,
    _nSexp,
    s_of_D

@testset "Drilling state space" begin

    @testset "Maxlease and state_if_never_drilled" begin
        for unitprob ∈ (LeasedProblem, LeasedProblemContsDrill,)
            for wpp in (unitprob(3,4, 5,3,2), unitprob(3,4, 2,3,2), unitprob(3,4, 1,3,5), unitprob(3,16, 5,10,3), unitprob(3,4,1,-1,0))
                SS = state_space_vector(wpp)
                @test maximum(_τrem(s) for s in SS) == maxlease(wpp)
                for sidx ∈ 1:length(wpp)
                    for d in actionspace(wpp,sidx)
                        @test sprime(wpp,sidx,d) ∈ 1:length(wpp)
                    end
                end

                # test state if never drilled...
                for state0 = 1:end_ex0(wpp)
                    statet = state0
                    for t = 1:end_ex0(wpp)
                        statet = sprime(wpp,statet,0)
                        ifno_drill = state_if_never_drilled(wpp,state0,t)
                        @test statet == ifno_drill
                        # println("t,state0 = ($t,$state0): recursive sprime = $statet, state_if_never_drilled = $ifno_drill")
                    end
                end

            end
        end
    end

    @testset "Drilling transitions are ok" begin
        for unitprob ∈ (LeasedProblem, LeasedProblemContsDrill,PerpetualProblem,)
            wp = unitprob(4, 4, 5, 3, 2)
            nS = length(wp)

            # infill: d1 == 0 cascades down
            for i = end_lrn(wp)+_nstates_per_D(wp) : _nstates_per_D(wp) : nS-1
                @test sprime(wp,i,0) == i
            end

            # infill: d1 == 1
            if isa(wp,LeasedProblemContsDrill)
                for i = end_lrn(wp)+1 : 2 : nS-2
                    @test sprime(wp,i,0) == i+1
                    for d in 1:_dmax(wp,i)
                        @test sprime(wp,i,d) == i+2*d
                    end
                end
            end

            # initial lease term
            if isa(wp,PerpetualProblem)
                @test sprime(wp,1,0) == 1
            else
                for i = 1 :  end_ex1(wp)-1
                    @test sprime(wp,i,0) == i+1
                end

                # last lease term
                for i = end_ex1(wp)+1 : end_ex0(wp)
                    @test sprime(wp,i,0) == i+1
                end

                # any drilling during a primary term
                for i = 0+1 : end_ex0(wp)
                    for d = 1:4
                        @test sprime(wp,i,d) == end_ex0(wp)+1+d
                    end
                end
            end
        end
    end

    @testset "sprime is in 1:length(wp)?" begin
        for unitprob ∈ (LeasedProblem, LeasedProblemContsDrill, PerpetualProblem,)
            for dmx in 3:4
                wp = unitprob(dmx, 4, 5, 3, 2)
                SS = state_space_vector(wp)

                @test length(inf_fm_lrn(wp)) == length(exploratory_learning(wp))

                for (i,s) in enumerate(SS)
                    for d in  actionspace(wp,i)
                        sp = sprime(wp,i,d)
                        s = SS[i]  # this state
                        t = SS[sp] # next state
                        sp ∈ 1:length(wp) || println("dmax = $dmx d = $d, i = $i with si = $(SS[i]), sprime = $(sprime(wp,i,d))")
                        @test sp ∈ 1:length(wp)
                        @test t.D == s.D + d

                        @test SS[i] == state(wp,i)
                        @test SS[sp] == state(wp,sp)

                        if  i < length(wp) && i ∉ ind_lrn(wp) && isa(wp,LeasedProblemContsDrill)
                            @test t.d1 == Int(d > 0)
                        end

                        if i ∈ ind_ex1(wp)
                            @test t.τ1 == ( d == 0 ? s.τ1 -1 : -1)
                            @test t.τ0 == ( d == 0 ? s.τ0 - (i == end_ex1(wp)) : -1)
                        end

                        if i ∈ ind_ex0(wp)
                            @test t.τ1 ==  s.τ1 == -1
                            if isa(wp, PerpetualProblem)
                                @test t.τ0 == -1
                            else
                                @test t.τ0 == ( d == 0 ? s.τ0 - 1 : -1)
                            end
                        end

                    end
                end
            end
        end
    end

    @testset "make State Space" begin

        let dmx = 4,
            Dmx = 4,
            τ0mx = 5,
            τ1mx = 3,
            emx = 2

            for unitprob ∈ (LeasedProblem, LeasedProblemContsDrill, PerpetualProblem,)
                wp = unitprob(dmx, Dmx, τ0mx, τ1mx, emx)

                SS = state_space_vector(wp)

                SS[ 1:end_ex1(wp) ]                                # Exploratory 1
                SS[ 1:end_ex0(wp) ]                                # Exploratory 0
                SS[  end_ex0(wp) .+ (1:dmx+1)]                     # Exploratory Terminal + learning update
                # SS[  (τ0mx + τ1mx + dmx+3) .+ (1:2*Dmx-2)]   # Infill drilling
                # SS[  (τ0mx + τ1mx + dmx+3)  +   2*Dmx-2+1]   # Terminal

                nS = length(SS)
                [actionspace(wp,i) for i in 1:nS]

                a = SS, [
                    (i,
                    _regime(wp,i),
                    _τ1(wp,i),
                    _τ0(wp,i),
                    _D(wp,i),
                    _d1(wp,i),
                    collect(sprime(wp,i,d) for d in actionspace(wp,i))
                    ) for i in 1:nS]

                idxs = [ind_exp(wp)..., ind_inf(wp)..., end_inf(wp)..., ind_lrn(wp)...]

                @test idxs ⊆ 1:length(wp)
                @test 1:length(wp) ⊆ idxs

                # test that we get back the state we want to get.
                for i in [ind_exp(wp)..., ind_inf(wp)..., end_inf(wp)...,]
                    st = SS[i]
                    i_of_st = state_idx(wp, st.τ1, st.τ0, st.D, st.d1)
                    @test i_of_st == i
                    @test i ∈ s_of_D(wp, st.D)
                end

                for s in ind_inf(wp)
                    @test _horizon(wp,s) ∈ (:Infinite, :Finite)
                end

                for s in ind_exp(wp)
                    if isa(wp,PerpetualProblem)
                        @test _horizon(wp,s) == :Infinite
                    else
                        @test _horizon(wp,s) == :Finite
                    end
                end

                for s in ind_lrn(wp)
                    @test _horizon(wp,s) ∈ (:Terminal, :Learning)
                end
            end
        end
    end

    @testset "Checking on dims of dEVσ" begin
        for unitprob ∈ (LeasedProblem, LeasedProblemContsDrill,)
            for wpp in (unitprob(3,4, 5,3,2), unitprob(3,4, 2,3,2), unitprob(3,4, 1,3,5), unitprob(3,16, 5,10,3), unitprob(3,4,1,-1,0))
                nsexp = _nSexp(wpp)
                for sidx ∈ 1:length(wpp)
                    if sidx <= end_ex0(wpp)
                        for d in actionspace(wpp,sidx)
                            @test sprime(wpp,sidx,d) ∈ 1:nsexp
                        end
                    end
                end
            end
        end
    end
end

end


text:
#jl #! format: off
# # Serializing PowerSystem Data

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSystems.jl supports serializing/deserializing data with JSON. This notebook
# provides an example of how to write and read a `System` to/from disk.

# ### Dependencies
# Let's use a dataset from the [tabular data parsing example](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/2_PowerSystems_examples/parse_matpower.ipynb)

using PowerSystems
using TimeSeries

base_dir = PowerSystems.download(PowerSystems.TestData; branch = "master");
sys = System(joinpath(base_dir, "matpower", "case5_re.m"))
sys

# ### Write data to a temporary directory

folder = mktempdir()
path = joinpath(folder, "system.json")
println("Serializing to $path")
to_json(sys, path)

filesize(path) / (1024 * 1024) #MiB

# ### Read the JSON file and create a new `System`
sys2 = System(path)

max_stars_repo_name: NREL-SIIP/SIIP-Tutorial

max_issues_repo_name: NREL-SIIP/SIIP-Tutorial

max_forks_repo_name: NREL-SIIP/SIIP-Tutorial


text:
using LinearAlgebra
using Statistics

abstract type ClassicalSystem end
abstract type Propagator end
abstract type Phase end

function kinetic end
function potential end

struct HamiltonianPhase <: Phase
    q::Vector
    p::Vector
end

"""
    Integrator(propagators...)

Constructs an integrator from a series of propagator, so that within a time step these propagators act sequentially on the phase.
"""
struct Integrator
    propagators::Tuple{Vararg{Propagator}}
    Integrator(propagators...) = begin
        new(propagators)
    end
end
(integrator!::Integrator)(phase::Phase, system::ClassicalSystem) = begin
    for propagator! in integrator!.propagators
        propagator!(phase, system)
    end
end

"""
    md(initPhase, system, integrator!, timetable, properties)

Launches an MD simulation of `system` with the initial conditions `initPhase` and collect `properties`. During the simulation, `integrator!` is used to integrate the equations of motion. With `timetable = (equilibrationSteps, sampleSteps)`, the system is first equilibrated for `equilibrationSteps` and then sampled for `sampleSteps`.
"""
md(initPhase::Phase, system::ClassicalSystem, integrator!::Integrator, timetable::Tuple{<:Integer, <:Integer}, properties::Vector{<:Function}) = begin
    phase = deepcopy(initPhase)
    nEquil, nSample = timetable
    for i in 1:nEquil
        integrator!(phase, system)
    end
    data = Dict(property => [] for property in properties)
    for i in 1:nSample
        integrator!(phase, system)
        for property in properties
            push!(data[property], property(phase, system))
        end
    end
    Dict(property => mean(data[property]) for property in properties)
end

statistics(properties::Vector{<:Function}, data::Vector{<:Dict}) = begin
    Dict(
        property => (
            mean=mean([item[property] for item in data]),
            std=std([item[property] for item in data]) / √(length(data))
        )
        for property in properties
    )
end


text:
using ForwardDiff

import Distributions: AbstractMvNormal

abstract type AbstractStateSpaceModel end
const AbstractSSM = AbstractStateSpaceModel
abstract type AbstractGaussianSSM <: AbstractStateSpaceModel end

# LinearGaussianSSM     LinearKalmanFilter
#                       NonlinearKalmanFilter
#                       NonlinearFilter
#
# NonlinearGaussianSSM  NonlinearKalmanFilter
#                       NonlinearFilter
#
# NonlinearSSM          NonlinearFilter

###########################################################################
# Linear Gaussian
###########################################################################

struct LinearGaussianSSM <: AbstractGaussianSSM
    # Process transition matrix, control matrix, and noise covariance
    F :: Function
    B :: Function
    V :: Function

    # Observation matrix and noise covariance
    G :: Function
    W :: Function

    # Model dimensions
    nx :: Int
    ny :: Int
    nu :: Int

    function LinearGaussianSSM( F :: Function,
                                B :: Function,
                                V :: Function,
                                G :: Function,
                                W :: Function )

        @assert ispossemidef(V(1))
        @assert ispossemidef(W(1))

        nx, ny, nu = confirm_matrix_sizes(F(1), B(1), V(1), G(1), W(1))
        new(F, B, V, G, W, nx, ny, nu)

    end

    ## Time-dependent constructor
    function LinearGaussianSSM( F :: Function,
                                V :: Function,
                                G :: Function,
                                W :: Function )

        B = _ -> begin zeros(T, (size(V(1), 1), 1)) end where T

        LinearGaussianSSM(F, B, V, G, W)

    end

    # Time-independent constructor
    function LinearGaussianSSM( F :: Matrix{T},
                                V :: Matrix{T},
                                G :: Matrix{T},
                                W :: Matrix{T}) where T

        B = zeros(T, size(F, 1), 1)

        LinearGaussianSSM(_->F, _->B, _->V, _->G, _->W)

    end

    function LinearGaussianSSM( F :: Matrix{T},
                                V :: Matrix{T},
                                G :: Matrix{T},
                                W :: Matrix{T},
                                B :: Matrix{T}) where T

        LinearGaussianSSM(_->F, _->B, _->V, _->G, _->W)

    end

end

import Base.show

function show(io::IO, mod::LinearGaussianSSM)

    dx, dy = mod.nx, mod.ny
    println("LinearGaussianSSM, $dx-D process x $dy-D observations")
    println("Process evolution matrix F:")
    show(mod.F(1))
    println("\n\nGontrol input matrix B:")
    show(mod.B(1))
    println("\n\nProcess error covariance V:")
    show(mod.V(1))
    println("\n\nObservation matrix G:")
    show(mod.G(1))
    println("\n\nObseration error covariance W:")
    show(mod.W(1))

end

## Core methods
process_matrix(m::LinearGaussianSSM, state::Vector, t::Real=0.0) = m.F(t)
process_matrix(m::LinearGaussianSSM, state::AbstractMvNormal, t::Real=0.0) = m.F(t)

observation_matrix(m::LinearGaussianSSM, state::Vector, t::Real=0.0) = m.G(t)
observation_matrix(m::LinearGaussianSSM, state::AbstractMvNormal, t::Real=0.0) = m.G(t)

control_input(m::LinearGaussianSSM, u, t::Real=0.0) = m.B(t) * u

###########################################################################
# Nonlinear Gaussian
###########################################################################

struct NonlinearGaussianSSM <: AbstractGaussianSSM
    # Process transition function, jacobian, and noise covariance matrix
    f :: Function
    V :: Function

    # Control input function
    b :: Function

    # Observation function,
    g :: Function # actual observation function
    W :: Function

    nx :: Int
    ny :: Int
    nu :: Int

    function NonlinearGaussianSSM( f  :: Function,
                                   V  :: Function,
                                   b  :: Function,
                                   g  :: Function,
                                   W  :: Function,
                                   nx :: Int,
                                   ny :: Int,
                                   nu :: Int)

        @assert ispossemidef(V(1))
        @assert ispossemidef(W(1))
        @assert (nx, nx) == size(V(1))
        @assert (ny, ny) == size(W(1))

        new(f, V, b, g, W, nx, ny, nu)

    end

    function NonlinearGaussianSSM( f  :: Function,
                                   V  :: Matrix{T},
                                   g  :: Function,
                                   W  :: Matrix{T};
                                   nu :: Int=0) where T

        b  = _ -> 0
        nx = size(V, 1)
        ny = size(W, 1)

        new(f, _-> V, b, g, _-> W, nx, ny, nu)

    end

end

## Core methods
function process_matrix( m :: NonlinearGaussianSSM,
                         x :: Vector,
                         t :: Real=0.0)

    ForwardDiff.jacobian(m.f, x)

end

function process_matrix( m :: NonlinearGaussianSSM,
                         x :: AbstractMvNormal,
                         t :: Real=0.0)

    process_matrix(m, mean(x), t)

end

function observation_matrix( m :: NonlinearGaussianSSM,
                             x :: Vector,
                             t :: Real=0.0)

    ForwardDiff.jacobian(m.g, x)

end

function observation_matrix( m :: NonlinearGaussianSSM,
                             x :: AbstractMvNormal,
                             t :: Real=0.0)

    process_matrix(m, mean(x), t)

end

function control_input( m :: NonlinearGaussianSSM,
                        u,
                        t :: Real=0.0 )
    m.b(u)

end


text:
# No copy since I do not modify anything and a copy is done when building a polyhedron
type LPHRepresentation{N, T} <: HRepresentation{N, T}
  # Ax <= b
  A::AbstractMatrix{T}
  l::AbstractVector{T}
  u::AbstractVector{T}
  colleqs::IntSet
  colgeqs::IntSet
  coleqs::IntSet
  lb::AbstractVector{T}
  ub::AbstractVector{T}
  rowleqs::IntSet
  rowgeqs::IntSet
  roweqs::IntSet

  function LPHRepresentation(A::AbstractMatrix{T}, l::AbstractVector{T}, u::AbstractVector{T}, lb::AbstractVector{T}, ub::AbstractVector{T})
    if length(l) != length(u) || size(A, 1) != length(l)
      error("The length of l and u must be equal to the number of rows of A")
    end
    if length(lb) != length(ub) || size(A, 2) != length(lb)
      error("The length of lb and ub must be equal to the number of columns of A")
    end
    if size(A, 2) != N
      error("Type dimension does not match the number of rows of A")
    end
    colleqs = IntSet()
    colgeqs = IntSet()
    coleqs = IntSet()
    for i in 1:N
      if l[i] > typemin(T)
        push!(colgeqs, i)
      end
      if u[i] < typemax(T)
        push!(colleqs, i)
      end
      if l[i] > typemin(T) && u[i] < typemax(T) && myeq(l[i], u[i])
        push!(coleqs, i)
      end
    end
    rowleqs = IntSet()
    rowgeqs = IntSet()
    roweqs = IntSet()
    for i in 1:size(A, 1)
      if lb[i] > typemin(T)
        push!(rowgeqs, i)
      end
      if ub[i] < typemax(T)
        push!(rowleqs, i)
      end
      if lb[i] > typemin(T) && ub[i] < typemax(T) && myeq(lb[i], ub[i])
        push!(roweqs, i)
      end
    end
    new(A, l, u, colleqs, colgeqs, coleqs, lb, ub, rowleqs, rowgeqs, roweqs)
  end
end

LPHRepresentation{T <: Real}(A::AbstractMatrix{T}, l::AbstractVector{T}, u::AbstractVector{T}, lb::AbstractVector{T}, ub::AbstractVector{T}) = LPHRepresentation{size(A,2),T}(A, l, u, lb, ub)
function LPHRepresentation(A::AbstractMatrix, l::AbstractVector, u::AbstractVector, lb::AbstractVector, ub::AbstractVector)
  T = promote_type(eltype(A), eltype(l), eltype(u), eltype(lb), eltype(ub))
  LPHRepresentation{size(A,2),T}(AbstractMatrix{T}(A), AbstractVector{T}(l), AbstractVector{T}(u), AbstractVector{T}(lb), AbstractVector{T}(ub))
end

function LPHRepresentation{N, T}(it::HRepIterator{N, T})
  A = Matrix{T}(length(it), N)
  lb = Vector{T}(length(it))
  ub = Vector{T}(length(it))
  MathProgBase.warn_no_inf(T)
  l = fill(typemin(T), N)
  u = fill(typemax(T), N)
  for (i, h) in enumerate(it)
    A[i,:] = h.a
    ub[i] = h.β
    if islin(h)
      lb[i] = ub[i]
    else
      lb[i] = typemin(T)
    end
  end
  LPHRepresentation{N, T}(A, l, u, lb, ub)
end
function LPHRepresentation{N, T}(;eqs::Nullable{EqIterator{N, T}}=nothing, ineq::Nullable{IneqIterator{N, T}}=nothing)
  neq = isnull(eqs) ? 0 : length(eqs)
  nineq = isnull(ineqs) ? 0 : length(ineqs)
  nhrep = neq + nineq
  A = Matrix{T}(nhrep, N)
  lb = Vector{T}(nhrep)
  ub = Vector{T}(nhrep)
  MathProgBase.warn_no_inf(T)
  l = fill(typemin(T), N)
  u = fill(typemax(T), N)
  if !(eqs === nothing)
    for (i, h) in enumerate(get(eqs))
      A[i,:] = h.a
      lb[i] = h.β
      ub[i] = h.β
    end
  end
  if !(ineqs === nothing)
    for (i, h) in enumerate(get(ineqs))
      A[neq+i,:] = h.a
      lb[neq+i] = typemin(T)
      ub[neq+i] = h.β
    end
  end
  new(A, l, u, lb, ub)
end

Base.copy{N,T}(lp::LPHRepresentation{N,T}) = LPHRepresentation{N,T}(copy(A), copy(l), copy(u), copy(colleqs), copy(colgeqs), copy(coleqs), copy(lb), copy(ub), copy(rowleqs), copy(rowgeqs), copy(roweqs))

function checknext(lp::LPHRepresentation, colrow, i, lgeq, allowed)
  lgeq += 1
  ok = false
  while colrow <= 2 && !ok
    if colrow == 1
      lgeqs = (lp.colleqs, lp.colgeqs, lp.coleqs)
    else
      lgeqs = (lp.rowleqs, lp.rowgeqs, lp.roweqs)
    end
    while i <= (colrow == 1 ? size(lp.A, 2) : size(lp.A, 1)) && !ok
      while lgeq <= 3 && !ok
        if allowed(i) && i in lgeqs[lgeq]
          ok = true
        else
          lgeq += 1
        end
      end
      if !ok
        i += 1
        lgeq = 1
      end
    end
    if !ok
      colrow += 1
      i = 1
      lgeq = 1
    end
  end
  (colrow, i, lgeq)
end

neqs(lp::LPHRepresentation) = length(lp.coleqs) + length(lp.roweqs)
nineqs(lp::LPHRepresentation) = length(lp.colleqs) + length(lp.colgeqs) + length(lp.rowleqs) + length(lp.rowgeqs)

starthrep(lp::LPHRepresentation) = checknext(1, 0, 3, (i) -> true)
donehrep(lp::LPHRepresentation, state) = state[1] > 2
function nexthrep{N,T}(lp::LPHRepresentation{N,T}, state)
  colrow, i, lgeq = state[1], state[2], state[3]
  if state[1] == 1
    a = spzeros(T, N)
    a[i] = lgeq == 1 ? -one(T) : one(T)
    β = lgeq == 1 ? -lp.lb[i,:] : lp.ub[i,:]
  elseif state[1] == 1
    a = lgeq == 1 ? -lp.A[i,:] : lp.A[i,:]
    β = lgeq == 1 ? -lp.lb[i,:] : lp.ub[i,:]
  else
    error("The iterator is done")
  end
  (lgeq == 3 ? HyperPlane(a, β) : HalfSpace(a, β), checknext(colrow, i, lgeq, (i) -> true))
end

starteq(lp::LPHRepresentation) = checknext(1, 0, 3, (i) -> i == 3)
doneeq(lp::LPHRepresentation, state) = state[1] > 2
function nexteq{N,T}(lp::LPHRepresentation{N,T}, state)
  colrow, i, lgeq = state[1], state[2], state[3]
  if state[1] == 1
    a = spzeros(T, N)
    a[i] = lgeq == 1 ? -one(T) : one(T)
    β = lgeq == 1 ? -lp.lb[i,:] : lp.ub[i,:]
  elseif state[1] == 1
    a = lgeq == 1 ? -lp.A[i,:] : lp.A[i,:]
    β = lgeq == 1 ? -lp.lb[i,:] : lp.ub[i,:]
  else
    error("The iterator is done")
  end
  (HyperPlane(a, β), checknext(colrow, i, lgeq, (i) -> i == 3))
end

startineq(lp::LPHRepresentation) = checknext(1, 0, 3, (i) -> i <= 2)
doneineq(lp::LPHRepresentation, state) = state[1] > 2
function nextineq{N,T}(lp::LPHRepresentation{N,T}, state)
  colrow, i, lgeq = state[1], state[2], state[3]
  if state[1] == 1
    a = spzeros(T, N)
    a[i] = lgeq == 1 ? -one(T) : one(T)
    β = lgeq == 1 ? -lp.lb[i,:] : lp.ub[i,:]
  elseif state[1] == 1
    a = lgeq == 1 ? -lp.A[i,:] : lp.A[i,:]
    β = lgeq == 1 ? -lp.lb[i,:] : lp.ub[i,:]
  else
    error("The iterator is done")
  end
  (HalfSpace(a, β), checknext(colrow, i, lgeq, (i) -> i <= 2))
end


text:
@testset "threadsafe.jl" begin
    @testset "constructor" begin
        vi = VarInfo(gdemo_default)
        threadsafe_vi = @inferred DynamicPPL.ThreadSafeVarInfo(vi)

        @test threadsafe_vi.varinfo === vi
        @test threadsafe_vi.logps isa Vector{typeof(Ref(getlogp(vi)))}
        @test length(threadsafe_vi.logps) == Threads.nthreads()
        @test all(iszero(x[]) for x in threadsafe_vi.logps)
    end

    # TODO: Add more tests of the public API
    @testset "API" begin
        vi = VarInfo(gdemo_default)
        threadsafe_vi = DynamicPPL.ThreadSafeVarInfo(vi)

        lp = getlogp(vi)
        @test getlogp(threadsafe_vi) == lp

        acclogp!!(threadsafe_vi, 42)
        @test threadsafe_vi.logps[Threads.threadid()][] == 42
        @test getlogp(vi) == lp
        @test getlogp(threadsafe_vi) == lp + 42

        resetlogp!!(threadsafe_vi)
        @test iszero(getlogp(vi))
        @test iszero(getlogp(threadsafe_vi))
        @test all(iszero(x[]) for x in threadsafe_vi.logps)

        setlogp!!(threadsafe_vi, 42)
        @test getlogp(vi) == 42
        @test getlogp(threadsafe_vi) == 42
        @test all(iszero(x[]) for x in threadsafe_vi.logps)
    end

    @testset "model" begin
        println("Peforming threading tests with $(Threads.nthreads()) threads")

        x = rand(10_000)

        @model function wthreads(x)
            global vi_ = __varinfo__
            x[1] ~ Normal(0, 1)
            Threads.@threads for i in 2:length(x)
                x[i] ~ Normal(x[i - 1], 1)
            end
        end

        vi = VarInfo()
        wthreads(x)(vi)
        lp_w_threads = getlogp(vi)
        if Threads.nthreads() == 1
            @test vi_ isa VarInfo
        else
            @test vi_ isa DynamicPPL.ThreadSafeVarInfo
        end

        println("With `@threads`:")
        println("  default:")
        @time wthreads(x)(vi)

        # Ensure that we use `ThreadSafeVarInfo` to handle multithreaded observe statements.
        DynamicPPL.evaluate_threadsafe!!(
            wthreads(x),
            vi,
            SamplingContext(Random.GLOBAL_RNG, SampleFromPrior(), DefaultContext()),
        )
        @test getlogp(vi) ≈ lp_w_threads
        @test vi_ isa DynamicPPL.ThreadSafeVarInfo

        println("  evaluate_threadsafe!!:")
        @time DynamicPPL.evaluate_threadsafe!!(
            wthreads(x),
            vi,
            SamplingContext(Random.GLOBAL_RNG, SampleFromPrior(), DefaultContext()),
        )

        @model function wothreads(x)
            global vi_ = __varinfo__
            x[1] ~ Normal(0, 1)
            for i in 2:length(x)
                x[i] ~ Normal(x[i - 1], 1)
            end
        end

        vi = VarInfo()
        wothreads(x)(vi)
        lp_wo_threads = getlogp(vi)
        if Threads.nthreads() == 1
            @test vi_ isa VarInfo
        else
            @test vi_ isa DynamicPPL.ThreadSafeVarInfo
        end

        println("Without `@threads`:")
        println("  default:")
        @time wothreads(x)(vi)

        @test lp_w_threads ≈ lp_wo_threads

        # Ensure that we use `VarInfo`.
        DynamicPPL.evaluate_threadunsafe!!(
            wothreads(x),
            vi,
            SamplingContext(Random.GLOBAL_RNG, SampleFromPrior(), DefaultContext()),
        )
        @test getlogp(vi) ≈ lp_w_threads
        @test vi_ isa VarInfo

        println("  evaluate_threadunsafe!!:")
        @time DynamicPPL.evaluate_threadunsafe!!(
            wothreads(x),
            vi,
            SamplingContext(Random.GLOBAL_RNG, SampleFromPrior(), DefaultContext()),
        )
    end
end

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
# This is an idea for how to package thermodynamic states.
# Here's my rationale for splittlng it up like this:
# There are a lot of possible thermodynamic and transport properties. At any moment,
# some of those may have been calculated by a program and others may not have been.
# That's why all of the fields in these structs are Union{T_, Missing}.
# The structs themselves group properties by groups that indicate "If I know at
# least one or two of these properties, I can reasonably expect that I would be
# able to calculate the others if I wanted to."

function iscontainer(::AbstractPhaseState)
    return false
end

struct MinimalState{T_<:AbstractFloat} <: AbstractPhaseState # where T_ <: AbstractFloat
    P::Union{T_, Missing}
    V::Union{T_, Missing}
    ρ::Union{T_, Missing}
    T::Union{T_, Missing}
    Cp::Union{T_, Missing}
end

MinimalState{T_}() where T_ = MinimalState{T_}(missing, missing, missing, missing, missing)

struct TherodynamicPotentialsState{T_<:AbstractFloat} <: AbstractPhaseState
    S::Union{T_, Missing}
    U::Union{T_, Missing}
    H::Union{T_, Missing}
    A::Union{T_, Missing}
    G::Union{T_, Missing}
end

TherodynamicPotentialsState{T_}() where T_ = TherodynamicPotentialsState{T_}(missing, missing, missing, missing, missing)

struct FullPotentialsState{T_<:AbstractFloat} <: AbstractPhaseState
    minstate::MinimalState{T_}
    potentials::TherodynamicPotentialsState{T_}
end

iscontainer(::FullPotentialsState) = true

MinimalState() = MinimalState(missing, missing, missing, missing, missing)

struct ChemicalState{T_<:AbstractFloat} <: AbstractPhaseState
    μ::Union{T_, Missing}
    xᵢ::Union{T_, Missing}
    f::Union{T_, Missing}
    ϕ::Union{T_, Missing}
    aᵢ::Union{T_, Missing}
    γᵢ::Union{T_, Missing}
end

ChemicalState{T_}() where T_ = ChemicalState{T_}(missing, missing, missing, missing, missing)

struct FullChemicalState{T_<:AbstractFloat} <: AbstractPhaseState
    minstate::MinimalState{T_}
    potentials::TherodynamicPotentialsState{T_}
    chemstate::ChemicalState{T_}
end

iscontainer(::FullChemicalState) = true

struct TransportPropertiesState{T_<:AbstractFloat} <: AbstractPhaseState
    λ::Union{T_, Missing}
    μ::Union{T_, Missing}
    ν::Union{T_, Missing}
    α::Union{T_, Missing}
    Pr::Union{T_, Missing}
end

TransportPropertiesState{T_}() where T_ = TransportPropertiesState{T_}(missing, missing, missing, missing, missing)

struct FullState{T_<:AbstractFloat} <: AbstractPhaseState
    minstate::MinimalState{T_}
    potentials::TherodynamicPotentialsState{T_}
    chemstate::ChemicalState{T_}
    transtate::TransportPropertiesState{T_}
end

iscontainer(::FullState) = true