text:
# ## Data

import CSV
import Random
import TuringModels

using DataFrames

Random.seed!(1)

delim = ';'
data_path = joinpath(TuringModels.project_root, "data", "Kline.csv")
df = CSV.read(data_path, DataFrame; delim)

df.log_pop = log.(df.population)
df.contact_high = [contact == "high" ? 1 : 0 for contact in df.contact]
df

# ## Model

using Turing

@model m10_10stan(total_tools, log_pop, contact_high) = begin
    α ~ Normal(0, 100)
    βp ~ Normal(0, 1)
    βc ~ Normal(0, 1)
    βpc ~ Normal(0, 1)

    for i ∈ 1:length(total_tools)
        λ = exp(α + βp*log_pop[i] + βc*contact_high[i] +
            βpc*contact_high[i]*log_pop[i])
        total_tools[i] ~ Poisson(λ)
    end
end;

# ## Output

chns = sample(
    m10_10stan(df.total_tools, df.log_pop, df.contact_high),
    NUTS(),
    1000
)

# \defaultoutput{}

# ## Original output

m_10_10t_result = "
     mean   sd  5.5% 94.5% n_eff Rhat
 a    0.94 0.37  0.36  1.53  3379    1
 bp   0.26 0.04  0.21  0.32  3347    1
 bc  -0.08 0.84 -1.41  1.23  2950    1
 bpc  0.04 0.09 -0.10  0.19  2907    1
";


text:
#=
Run this file to regenerate `geos_c.jl` and `geos_common.jl`.

It expects a GEOS install in the deps folder, run `build LibGEOS` in Pkg mode
if these are not in place.

The wrapped GEOS version and provided GEOS version should be kept in sync.
So when updating the GEOSBuilder provided version, also rerun this wrapper.
This way we ensure that the provided library has the same functions available
as the wrapped one. Furthermore this makes sure constants in `geos_common.jl`
like `GEOS_VERSION`, which are just strings, are correct.
=#

using Clang
using MacroTools

includedir = normpath(joinpath(@__DIR__, "..", "deps", "usr", "include"))
headerfiles = [joinpath(includedir, "geos_c.h")]

"""
Custom rewriter for Clang.jl's C wrapper

Gets called with all expressions in a header file, or all expressions in a common file.
The expressions get sent to `rewriter(::Expr)`` for
further treatment.
"""
function rewriter(xs::Vector)
    rewritten = Any[]
    for x in xs
        # Clang.jl inserts strings like "# Skipping MacroDefinition: X"
        # keep these to get a sense of what we are missing
        if x isa String
            push!(rewritten, x)
            continue
        end
        @assert x isa Expr

        x2 = rewriter(x)
        push!(rewritten, x2)
    end
    rewritten
end

"Functions that return a Cstring are wrapped in unsafe_string to return a String"
function rewriter(x::Expr)
    if @capture(x,
        function f_(fargs__)
            ccall(fname_, rettype_, argtypes_, argvalues__)
        end
    )
        # it is a function wrapper around a ccall

        # bind the ccall such that we can easily wrap it
        cc = :(ccall($fname, $rettype, $argtypes, $(argvalues...)))

        cc2 = if rettype == :Cstring
            :(unsafe_string($cc))
        else
            cc
        end

        # stitch the modified function expression back together
        x2 = :(function $f($(fargs...))
            $cc2
        end) |> prettify
        x2
    else
        # do not modify expressions that are no ccall function wrappers
        x
    end
end

wc = init(; headers = headerfiles,
            output_file = joinpath(@__DIR__, "..", "src", "geos_c.jl"),
            common_file = joinpath(@__DIR__, "..", "src", "geos_common.jl"),
            clang_includes = [includedir, CLANG_INCLUDE],
            clang_args = ["-I", includedir],
            header_wrapped = (root, current) -> root == current,
            header_library = x -> "libgeos",
            clang_diagnostics = true,
            rewriter = rewriter,
            )

run(wc)

# delete Clang.jl helper files
rm(joinpath(@__DIR__, "..", "src", "LibTemplate.jl"))
rm(joinpath(@__DIR__, "..", "src", "ctypes.jl"))

ext: jl


text:
function sortslicesperm(A::AbstractArray; dims::Union{Integer, Tuple{Vararg{Integer}}}, kws...)
    _sortslicesperm(A, Val{dims}(); kws...)
end

function _sortslicesperm(A::AbstractArray, d::Val{dims}; kws...) where dims
    itspace = Base.compute_itspace(A, d)
    vecs = map(its->view(A, its...), itspace)
    p = sortperm(vecs; kws...)
    if ndims(A) == 2 && isa(dims, Integer) && isa(A, Array)
        # At the moment, the performance of the generic version is subpar
        # (about 5x slower). Hardcode a fast-path until we're able to
        # optimize this.
        # return dims == 1 ? A[p, :] : A[:, p]
        return p
    else
        B = similar(A)
        for (x, its) in zip(p, itspace)
            B[its...] = vecs[x]
        end
        B
    end
end

function calculate_hv(points, bounds, k, n_sample)

    # This function is modified from the code in
    # http://www.tik.ee.ethz.ch/sop/download/supplementary/hype/

    N, M = size(points)
    if M > 2
        # Use the estimated method for three or more objectives
        alpha = zeros(1,N);
        for i = 1:k
            J = 1:i-1
            alpha[i] = prod((k .- J)./(N .- J ))./i
        end

        f_min = ideal(points)
        S = f_min' .+ (bounds - f_min)' .* rand(n_sample, M)

        PdS  = zeros(Bool, N,n_sample)
        dS   = zeros(Int, n_sample)
        for i = 1:N
            x = sum(points[i,:]' .- S .<= 0, dims = 2) .== M
            mask = view(x, :,1)
            PdS[i, mask] .= true
            dS[mask]    = dS[mask] .+ 1
        end

        F = zeros(N)
        for i = 1:N
            mask = view(dS, view(PdS, i,:))
            F[i] = sum(alpha[mask])
        end

        # ΔS
        return F .* prod(bounds - f_min) / n_sample
    end

    return zeros(0)

    #= ignored since  computation is carried out in `compute_contribution!(...) when M == 2`
    #
    #
    # Use the accurate method for two objectives
    pvec  = 1:size(points,1);
    alpha = zeros(1,k);
    for i = 1 : k
        j = 1 : i-1;
        alpha[i] = prod((k-j)./(N-j))./i;
    end

    return hypesub(N,points,M,bounds,pvec,alpha,k);
    =#
end

#=
function hypesub(l,A,M,bounds,pvec,alpha,k)
    # The recursive function for the accurate method

    h     = zeros(1,l);
    #[S,mask] = sortrows(A,M);
    _A = view(A, :, M)
    maks = sortslicesperm( _A, dims=1 )
    S = A[mask,:]

    pvec  = pvec[mask];
    for i = 1 : size(S,1)
        if i < size(S,1)
            extrusion = S[i+1,M] - S[i,M];
        else
            extrusion = bounds[M] - S[i,M];
        end
        if M == 1
            if i > k
                break;
            end
            if alpha >= 0
                h[pvec[1:i]] = h[pvec[1:i]] + extrusion*alpha[i];
            end
        elseif extrusion > 0
            h += extrusion*hypesub(l,S[1:i,:],M-1,bounds,pvec[1:i],alpha,k);
        end
    end

    return h
end
=#


text:
module ITensorNetworks

using ITensors

include("ising_classical_2d.jl")
include("tensor_networks.jl")
include("boundary_mps_projectors.jl")

end


text:
# Julia wrapper for header: /usr/include/libavformat/avio.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    avio_find_protocol_name,
    avio_check,
    avpriv_io_move,
    avpriv_io_delete,
    avio_open_dir,
    avio_read_dir,
    avio_close_dir,
    avio_free_directory_entry,
    avio_alloc_context,
    avio_context_free,
    avio_w8,
    avio_write,
    avio_wl64,
    avio_wb64,
    avio_wl32,
    avio_wb32,
    avio_wl24,
    avio_wb24,
    avio_wl16,
    avio_wb16,
    avio_put_str,
    avio_put_str16le,
    avio_put_str16be,
    avio_write_marker,
    avio_seek,
    avio_skip,
    avio_tell,
    avio_size,
    avio_feof,
    avio_flush,
    avio_read,
    avio_read_partial,
    avio_r8,
    avio_rl16,
    avio_rl24,
    avio_rl32,
    avio_rl64,
    avio_rb16,
    avio_rb24,
    avio_rb32,
    avio_rb64,
    avio_get_str,
    avio_get_str16le,
    avio_get_str16be,
    avio_open,
    avio_open2,
    avio_close,
    avio_closep,
    avio_open_dyn_buf,
    avio_get_dyn_buf,
    avio_close_dyn_buf,
    avio_enum_protocols,
    avio_pause,
    avio_seek_time,
    avio_read_to_bprint,
    avio_accept,
    avio_handshake

function avio_find_protocol_name(url)
    ccall((:avio_find_protocol_name, libavformat), Cstring, (Cstring,), url)
end

function avio_check(url, flags::Integer)
    ccall((:avio_check, libavformat), Cint, (Cstring, Cint), url, flags)
end

function avpriv_io_move(url_src, url_dst)
    ccall((:avpriv_io_move, libavformat), Cint, (Cstring, Cstring), url_src, url_dst)
end

function avpriv_io_delete(url)
    ccall((:avpriv_io_delete, libavformat), Cint, (Cstring,), url)
end

function avio_open_dir(s, url, options)
    ccall((:avio_open_dir, libavformat), Cint, (Ptr{Ptr{AVIODirContext}}, Cstring, Ref{Ptr{AVDictionary}}), s, url, options)
end

function avio_read_dir(s, next)
    ccall((:avio_read_dir, libavformat), Cint, (Ptr{AVIODirContext}, Ptr{Ptr{AVIODirEntry}}), s, next)
end

function avio_close_dir(s)
    ccall((:avio_close_dir, libavformat), Cint, (Ptr{Ptr{AVIODirContext}},), s)
end

function avio_free_directory_entry(entry)
    ccall((:avio_free_directory_entry, libavformat), Cvoid, (Ptr{Ptr{AVIODirEntry}},), entry)
end

function avio_alloc_context(buffer, buffer_size::Integer, write_flag::Integer, opaque, read_packet, write_packet, seek)
    ccall((:avio_alloc_context, libavformat), Ptr{AVIOContext}, (Ptr{Cuchar}, Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}), buffer, buffer_size, write_flag, opaque, read_packet, write_packet, seek)
end

function avio_context_free(s)
    ccall((:avio_context_free, libavformat), Cvoid, (Ptr{Ptr{AVIOContext}},), s)
end

function avio_w8(s, b::Integer)
    ccall((:avio_w8, libavformat), Cvoid, (Ptr{AVIOContext}, Cint), s, b)
end

function avio_write(s, buf, size::Integer)
    ccall((:avio_write, libavformat), Cvoid, (Ptr{AVIOContext}, Ptr{Cuchar}, Cint), s, buf, size)
end

function avio_wl64(s, val::UInt64)
    ccall((:avio_wl64, libavformat), Cvoid, (Ptr{AVIOContext}, UInt64), s, val)
end

function avio_wb64(s, val::UInt64)
    ccall((:avio_wb64, libavformat), Cvoid, (Ptr{AVIOContext}, UInt64), s, val)
end

function avio_wl32(s, val::Integer)
    ccall((:avio_wl32, libavformat), Cvoid, (Ptr{AVIOContext}, UInt32), s, val)
end

function avio_wb32(s, val::Integer)
    ccall((:avio_wb32, libavformat), Cvoid, (Ptr{AVIOContext}, UInt32), s, val)
end

function avio_wl24(s, val::Integer)
    ccall((:avio_wl24, libavformat), Cvoid, (Ptr{AVIOContext}, UInt32), s, val)
end

function avio_wb24(s, val::Integer)
    ccall((:avio_wb24, libavformat), Cvoid, (Ptr{AVIOContext}, UInt32), s, val)
end

function avio_wl16(s, val::Integer)
    ccall((:avio_wl16, libavformat), Cvoid, (Ptr{AVIOContext}, UInt32), s, val)
end

function avio_wb16(s, val::Integer)
    ccall((:avio_wb16, libavformat), Cvoid, (Ptr{AVIOContext}, UInt32), s, val)
end

function avio_put_str(s, str)
    ccall((:avio_put_str, libavformat), Cint, (Ptr{AVIOContext}, Cstring), s, str)
end

function avio_put_str16le(s, str)
    ccall((:avio_put_str16le, libavformat), Cint, (Ptr{AVIOContext}, Cstring), s, str)
end

function avio_put_str16be(s, str)
    ccall((:avio_put_str16be, libavformat), Cint, (Ptr{AVIOContext}, Cstring), s, str)
end

function avio_write_marker(s, time::Int64, _type::AVIODataMarkerType)
    ccall((:avio_write_marker, libavformat), Cvoid, (Ptr{AVIOContext}, Int64, AVIODataMarkerType), s, time, _type)
end

function avio_seek(s, offset::Int64, whence::Integer)
    ccall((:avio_seek, libavformat), Int64, (Ptr{AVIOContext}, Int64, Cint), s, offset, whence)
end

function avio_skip(s, offset::Int64)
    ccall((:avio_skip, libavformat), Int64, (Ptr{AVIOContext}, Int64), s, offset)
end

function avio_tell(s)
    ccall((:avio_tell, libavformat), Int64, (Ptr{AVIOContext},), s)
end

function avio_size(s)
    ccall((:avio_size, libavformat), Int64, (Ptr{AVIOContext},), s)
end

function avio_feof(s)
    ccall((:avio_feof, libavformat), Cint, (Ptr{AVIOContext},), s)
end

function avio_flush(s)
    ccall((:avio_flush, libavformat), Cvoid, (Ptr{AVIOContext},), s)
end

function avio_read(s, buf, size::Integer)
    ccall((:avio_read, libavformat), Cint, (Ptr{AVIOContext}, Ptr{Cuchar}, Cint), s, buf, size)
end

function avio_read_partial(s, buf, size::Integer)
    ccall((:avio_read_partial, libavformat), Cint, (Ptr{AVIOContext}, Ptr{Cuchar}, Cint), s, buf, size)
end

function avio_r8(s)
    ccall((:avio_r8, libavformat), Cint, (Ptr{AVIOContext},), s)
end

function avio_rl16(s)
    ccall((:avio_rl16, libavformat), UInt32, (Ptr{AVIOContext},), s)
end

function avio_rl24(s)
    ccall((:avio_rl24, libavformat), UInt32, (Ptr{AVIOContext},), s)
end

function avio_rl32(s)
    ccall((:avio_rl32, libavformat), UInt32, (Ptr{AVIOContext},), s)
end

function avio_rl64(s)
    ccall((:avio_rl64, libavformat), UInt64, (Ptr{AVIOContext},), s)
end

function avio_rb16(s)
    ccall((:avio_rb16, libavformat), UInt32, (Ptr{AVIOContext},), s)
end

function avio_rb24(s)
    ccall((:avio_rb24, libavformat), UInt32, (Ptr{AVIOContext},), s)
end

function avio_rb32(s)
    ccall((:avio_rb32, libavformat), UInt32, (Ptr{AVIOContext},), s)
end

function avio_rb64(s)
    ccall((:avio_rb64, libavformat), UInt64, (Ptr{AVIOContext},), s)
end

function avio_get_str(pb, maxlen::Integer, buf, buflen::Integer)
    ccall((:avio_get_str, libavformat), Cint, (Ptr{AVIOContext}, Cint, Cstring, Cint), pb, maxlen, buf, buflen)
end

function avio_get_str16le(pb, maxlen::Integer, buf, buflen::Integer)
    ccall((:avio_get_str16le, libavformat), Cint, (Ptr{AVIOContext}, Cint, Cstring, Cint), pb, maxlen, buf, buflen)
end

function avio_get_str16be(pb, maxlen::Integer, buf, buflen::Integer)
    ccall((:avio_get_str16be, libavformat), Cint, (Ptr{AVIOContext}, Cint, Cstring, Cint), pb, maxlen, buf, buflen)
end

function avio_open(s, url, flags::Integer)
    ccall((:avio_open, libavformat), Cint, (Ptr{Ptr{AVIOContext}}, Cstring, Cint), s, url, flags)
end

function avio_open2(s, url, flags::Integer, int_cb, options)
    ccall((:avio_open2, libavformat), Cint, (Ptr{Ptr{AVIOContext}}, Cstring, Cint, Ptr{AVIOInterruptCB}, Ref{Ptr{AVDictionary}}), s, url, flags, int_cb, options)
end

function avio_close(s)
    ccall((:avio_close, libavformat), Cint, (Ptr{AVIOContext},), s)
end

function avio_closep(s)
    ccall((:avio_closep, libavformat), Cint, (Ptr{Ptr{AVIOContext}},), s)
end

function avio_open_dyn_buf(s)
    ccall((:avio_open_dyn_buf, libavformat), Cint, (Ptr{Ptr{AVIOContext}},), s)
end

function avio_get_dyn_buf(s, pbuffer)
    ccall((:avio_get_dyn_buf, libavformat), Cint, (Ptr{AVIOContext}, Ptr{Ptr{UInt8}}), s, pbuffer)
end

function avio_close_dyn_buf(s, pbuffer)
    ccall((:avio_close_dyn_buf, libavformat), Cint, (Ptr{AVIOContext}, Ptr{Ptr{UInt8}}), s, pbuffer)
end

function avio_enum_protocols(opaque, output::Integer)
    ccall((:avio_enum_protocols, libavformat), Cstring, (Ptr{Ptr{Cvoid}}, Cint), opaque, output)
end

function avio_pause(h, pause::Integer)
    ccall((:avio_pause, libavformat), Cint, (Ptr{AVIOContext}, Cint), h, pause)
end

function avio_seek_time(h, stream_index::Integer, timestamp::Int64, flags::Integer)
    ccall((:avio_seek_time, libavformat), Int64, (Ptr{AVIOContext}, Cint, Int64, Cint), h, stream_index, timestamp, flags)
end

function avio_read_to_bprint(h, pb, max_size::Csize_t)
    ccall((:avio_read_to_bprint, libavformat), Cint, (Ptr{AVIOContext}, Ptr{AVBPrint}, Csize_t), h, pb, max_size)
end

function avio_accept(s, c)
    ccall((:avio_accept, libavformat), Cint, (Ptr{AVIOContext}, Ptr{Ptr{AVIOContext}}), s, c)
end

function avio_handshake(c)
    ccall((:avio_handshake, libavformat), Cint, (Ptr{AVIOContext},), c)
end


text:
using DWARF
using DataStructures
using Cxx

"""
Show the five instructions surrounding ipoffset (two before, two after), with an
indicator which instruction is at ipoffset.
"""
function disasm_around_ip(io, insts, ipoffset; ipbase = 0, circular = true, triple = getDefaultTargetTriple())
    Offset = 0
    ctx = DisAsmContext(triple)
    const InstInfo = Tuple{Int,Bool,AbstractString}
    buf = circular ? CircularBuffer{InstInfo}(5) : Vector{InstInfo}()
    targetn = typemax(Int64)
    ipinstoffset = 0
    while Offset < sizeof(insts) && targetn > 0
        (Inst, InstSize) = getInstruction(insts, Offset; ctx = ctx)
        lastoffset = Offset
        Offset += InstSize
        iobuf = IOContext(IOBuffer(),:disasmctx=>ctx)
        print(iobuf, Inst)
        push!(buf, (lastoffset, mayAffectControlFlow(Inst,ctx),
            takebuf_string(iobuf.io)))
        targetn -= 1
        if circular && Offset > ipoffset && targetn > 2 # Two more
            ipinstoffset = lastoffset
            targetn = 2
        end
        free(Inst)
    end
    for i = 1:length(buf)
        off, branch, inst = buf[i]
        attarget = off == ipinstoffset
        p = string(attarget ? "=> " : "   ",
            "0x",hex(UInt64(ipbase+off),2sizeof(UInt64)),"<+",off,">:",
            inst)
        if attarget
            print_with_color(:yellow, io, p)
        elseif branch
            print_with_color(:red, io, p)
        else
            print(io, p)
        end
        println(io)
    end
end

# From DIDebug
cxx"""
#include "llvm/MC/SubtargetFeature.h"
#include "llvm/MC/MCObjectFileInfo.h"
#include "llvm/MC/MCContext.h"
#include "llvm/MC/MCDisassembler.h"
#include "llvm/MC/MCInst.h"
#include "llvm/MC/MCInstPrinter.h"
#include "llvm/MC/MCAsmInfo.h"
#include "llvm/Support/TargetSelect.h"
#include "llvm/Support/TargetRegistry.h"
#include "llvm/MC/MCInstrInfo.h"
using namespace llvm;
"""

immutable DisAsmContext
  MAI::pcpp"llvm::MCAsmInfo"
  MRI::pcpp"llvm::MCRegisterInfo"
  MII::pcpp"llvm::MCInstrInfo"
  MOFI::pcpp"llvm::MCObjectFileInfo"
  MCtx::pcpp"llvm::MCContext"
  MSTI::pcpp"llvm::MCSubtargetInfo"
  DisAsm::pcpp"llvm::MCDisassembler"
  MIP::pcpp"llvm::MCInstPrinter"
end

getDefaultTargetTriple() = icxx"sys::getDefaultTargetTriple();"
function DisAsmContext(TripleName=getDefaultTargetTriple())
  icxx"""
    llvm::InitializeNativeTargetAsmParser();
    llvm::InitializeNativeTargetDisassembler();
  """

  TheTarget = icxx"""
    std::string err;
    TargetRegistry::lookupTarget($TripleName, err);
  """

  MAI  = icxx" $TheTarget->createMCAsmInfo(*$TheTarget->createMCRegInfo($TripleName),$TripleName); "
  MII  = icxx" $TheTarget->createMCInstrInfo(); "
  MRI  = icxx" $TheTarget->createMCRegInfo($TripleName); "
  MOFI = icxx" new MCObjectFileInfo; "
  MCtx = icxx" new MCContext($MAI, $MRI, $MOFI); "
  MSTI = icxx"""
    Triple TheTriple(Triple::normalize($TripleName));
    $MOFI->InitMCObjectFileInfo(TheTriple, Reloc::PIC_, CodeModel::Default, *$MCtx);
    SubtargetFeatures Features;
    Features.getDefaultSubtargetFeatures(TheTriple);
    std::string MCPU = sys::getHostCPUName();
    $TheTarget->createMCSubtargetInfo($TripleName, MCPU, Features.getString());
  """

  DisAsm = icxx" $TheTarget->createMCDisassembler(*$MSTI, *$MCtx); "

  MIP = icxx"""
      int AsmPrinterVariant = $MAI->getAssemblerDialect();
      $TheTarget->createMCInstPrinter(
          Triple($TripleName), AsmPrinterVariant, *$MAI, *$MII, *$MRI);
  """

  DisAsmContext(MAI, MRI, MII, MOFI, MCtx, MSTI, DisAsm, MIP)
end

function getInstruction(data::Vector{UInt8}, offset; ctx = DisAsmContext())
  Inst = icxx" new MCInst; "
  InstSize = icxx"""
    uint8_t *Base = (uint8_t*)$(convert(Ptr{Void},pointer(data)));
    uint64_t Total = $(sizeof(data));
    uint64_t Offset = $offset;
    uint64_t LoadAddress = 0;
    uint64_t InstSize;
    MCDisassembler::DecodeStatus S;
    S = $(ctx.DisAsm)->getInstruction(*$Inst, InstSize,
        ArrayRef<uint8_t>(Base+Offset, Total-Offset),
          LoadAddress + Offset, /*REMOVE*/ nulls(), nulls());
    switch (S) {
      case MCDisassembler::SoftFail:
      case MCDisassembler::Fail:
        delete $Inst;
        return (uint64_t)0;
      case MCDisassembler::Success:
        return InstSize;
    }
  """
  if InstSize == 0
    error("Invalid Instruction")
  end
  (Inst, InstSize)
end

function mayAffectControlFlow(inst, ctx)
  icxx"""
    $(ctx.MII)->get($inst->getOpcode()).mayAffectControlFlow(*$inst,*$(ctx.MRI));
  """
end

function mayBranch(inst, ctx)
  icxx"""
    $(ctx.MII)->get($inst->getOpcode()).isBranch();
  """
end

#=
function disassembleInstruction(f::Function, data::Vector{UInt8}, offset; ctx = DisAsmContext())
  icxx"""
    uint8_t *Base = (uint8_t*)$(convert(Ptr{Void},pointer(data)));
    uint64_t Total = $(sizeof(data));
    uint64_t Offset = $offset;
    uint64_t LoadAddress = 0;
    uint64_t InstSize;
    MCInst Inst;
    MCDisassembler::DecodeStatus S;
    S = $(ctx.DisAsm)->getInstruction(Inst, InstSize,
        ArrayRef<uint8_t>(Base+Offset, Total-Offset),
          LoadAddress + Offset, /*REMOVE*/ nulls(), nulls());
    switch (S) {
      case MCDisassembler::SoftFail:
      case MCDisassembler::Fail:
        return 0;
      case MCDisassembler::Success:
        :(f(icxx"&Inst;"); nothing);
        return InstSize;
    }
  """
end
=#

function Base.show(io::IO, Inst::pcpp"llvm::MCInst")
  ctx = isa(io,IOContext) && haskey(io, :disasmctx) ?
    get(io, :disasmctx, nothing) : DisAsmContext()
  print(io,String(icxx"""
  std::string O;
  raw_string_ostream OS(O);
  ($(ctx.MIP))->printInst(($Inst), OS, "", *$(ctx.MSTI));
  OS.flush();
  O;
  """))
end

free(x::pcpp"llvm::MCInst") = icxx" delete $x; "

function loclist2rangelist{T}(list::DWARF.LocationList{T})
  rangelist = Array(UnitRange{T},0)
  for entry in list.entries
    push!(rangelist, entry.first:entry.last)
  end
  rangelist
end

function rangelists(dbgs)
  seek(dbgs.oh, sectionoffset(dbgs.debug_loc))
  lists = Array(Vector{UnitRange{UInt64}},0)
  while position(dbgs.oh) < sectionoffset(dbgs.debug_loc)+sectionsize(dbgs.debug_loc)
      push!(lists,loclist2rangelist(read(dbgs.oh, DWARF.LocationList{UInt64})))
  end
  lists
end

function rangelist(dbgs, offset)
    seek(dbgs.oh, sectionoffset(dbgs.debug_loc)+offset)
    loclist2rangelist(read(dbgs.oh, DWARF.LocationList{UInt64}))
end

default_colors = [:blue, :red, :green, :yellow, :purple]
function disassemble2(data::Vector{UInt8}, instrange = nothing; io = STDOUT,
    ctx = DisAsmContext(), rloffset = 0, rangelists = [], colors = default_colors, offset = 0)
  Offset = instrange !== nothing ? first(instrange) - offset : 0
  while Offset < sizeof(data) && (instrange === nothing || Offset <= last(instrange)-offset)
    (Inst, InstSize) = getInstruction(data, Offset)
    print(io,"0x",hex(Offset+offset,2*sizeof(Offset)))
    print(io,":")
    str = repr(Inst; ctx = ctx)
    # This is bad, but I need things to line up
    lastt = findlast(str, '\t')
    rest = str[(lastt+1):end]
    print(io, str[1:lastt])
    if lastt == 1
      print(io,rest,'\t')
      rest = ""
    end
    printfield(io, rest, 25; align = :left)
    # Print applicable range lists
    for (i,rangelist) in enumerate(rangelists)
      print(io," ")
      found = false
      color = colors[mod1(i,length(colors))]
      for range in rangelist
        # If this is the start of a range
        if Offset+offset == first(range)+rloffset
          print_with_color(color,io,"x")
        # Or the last
        elseif last(range)+rloffset == Offset+offset+InstSize
          print_with_color(color,io,"x")
        elseif first(range)+rloffset <= Offset + offset < last(range)+rloffset
          print_with_color(color,io,"|")
        else
          continue
        end
        found = true
        break
      end
      found || print(io," ")
    end
    println(io)
    free(Inst)
    Inst = nothing
    Offset += InstSize
  end
end
function disassemble2(base, size; kwargs...)
    data = pointer_to_array(convert(Ptr{UInt8},base), (size,), false)
    disassemble2(data; kwargs...)
end

function extract_next_inst(insts::Vector{UInt8}, ctx = DisAsmContext())
  size = icxx"""
    uint64_t InstSize;
    uint8_t *Base = (uint8_t*)$(pointer(insts));
    uint64_t Total = $(sizeof(insts));
    uint64_t LoadAddress = 0;
    for (uint64_t Offset = 0; Offset < Total; Offset += InstSize)
    {
      MCInst Inst;
      MCDisassembler::DecodeStatus S;

      S = $(ctx.DisAsm)->getInstruction(Inst, InstSize,
          ArrayRef<uint8_t>(Base+Offset, Total-Offset),
            LoadAddress + Offset, /*REMOVE*/ nulls(), nulls());

      switch (S) {
      case MCDisassembler::SoftFail:
      case MCDisassembler::Fail:
        return (uint64_t)0;
      case MCDisassembler::Success:
          $(ctx.MIP)->printInst(&Inst, outs(), "", *$(ctx.MSTI));
        return InstSize;
      }
    }
    return (uint64_t)0;
  """
  insts[1:size]
end

function disassemble(insts::Vector{UInt8}, ctx = DisAsmContext())
  icxx"""
    uint64_t InstSize;
    uint8_t *Base = (uint8_t*)$(pointer(insts));
    uint64_t Total = $(sizeof(insts));
    uint64_t LoadAddress = 0;
    for (uint64_t Offset = 0; Offset < Total; Offset += InstSize)
    {
      MCInst Inst;
      MCDisassembler::DecodeStatus S;

      S = $(ctx.DisAsm)->getInstruction(Inst, InstSize,
          ArrayRef<uint8_t>(Base+Offset, Total-Offset),
            LoadAddress + Offset, /*REMOVE*/ nulls(), nulls());

      switch (S) {
      case MCDisassembler::SoftFail:
      case MCDisassembler::Fail:
        return false;
      case MCDisassembler::Success:
          $(ctx.MIP)->printInst(&Inst, outs(), "", *$(ctx.MSTI));
          outs() << "\n";
        break;
      }
    }
    return true;
  """
end


text:
using PiMC
using BenchmarkTools
using Random
using Printf

Random.seed!(100);

n_samples = 10^8;

Random.seed!(100);
pi_est = pi_estimation(n_samples);
@printf("π estimate with %.1e samples = %g on %d threads\n", n_samples,pi_est, Threads.nthreads());

@btime pi_estimation(n_samples);


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

function write_opwrap_lds_stor_init(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	## Extract data frames from input dictionary
	dfGen = inputs["dfGen"]
	G = inputs["G"]

	# Initial level of storage in each modeled period
	NPeriods = size(inputs["Period_Map"])[1]
	dfStorageInit = DataFrame(Resource = inputs["RESOURCES"], Zone = dfGen[!,:Zone])
	socw = zeros(G,NPeriods)
	for i in 1:G
		if i in inputs["STOR_LONG_DURATION"]
			socw[i,:] = value.(EP[:vSOCw])[i,:]
		end
	end
	dfStorageInit = hcat(dfStorageInit, DataFrame(socw, :auto))
	auxNew_Names=[Symbol("Resource");Symbol("Zone");[Symbol("n$t") for t in 1:NPeriods]]
	rename!(dfStorageInit,auxNew_Names)
	CSV.write(joinpath(path, "StorageInit.csv"), dftranspose(dfStorageInit, false), writeheader=false)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using Rebugger
using Rebugger: StopException
using Test, UUIDs, InteractiveUtils, REPL, Pkg, HeaderREPLs
using REPL.LineEdit
using Revise, Colors

if !isdefined(Main, :RebuggerTesting)
    includet("testmodule.jl")   # so the source code here gets loaded
end

const empty_kwvarargs = Rebugger.kwstasher()
uuidextractor(str) = UUID(match(r"getstored\(\"([a-z0-9\-]+)\"\)", str).captures[1])

struct ErrorsOnShow end
Base.show(io::IO, ::ErrorsOnShow) = throw(ArgumentError("no show"))

@testset "Rebugger" begin
    id = uuid1()
    @test uuidextractor("vars = getstored(\"$id\") and more stuff") == id
    @testset "Debug core" begin
        @testset "Deepcopy" begin
            args = (3.2, rand(3,3), Rebugger, [Rebugger], "hello", sum, (2,3))
            argc = Rebugger.safe_deepcopy(args...)
            @test argc == args
        end
        @testset "Signatures" begin
            @test Rebugger.signature_names!(:(f(x::Int, @nospecialize(y::String)))) == (:f, (:x, :y), (), ())
            @test Rebugger.signature_names!(:(f(x::Int, $(Expr(:meta, :nospecialize, :(y::String)))))) ==
                (:f, (:x, :y), (), ())
            ex = :(f(::Type{T}, ::IndexStyle, x::Int, ::IndexStyle) where T)
            @test Rebugger.signature_names!(ex) == (:f, (:T, :__IndexStyle_1, :x, :__IndexStyle_2), (), ())
            @test ex == :(f(::Type{T}, __IndexStyle_1::IndexStyle, x::Int, __IndexStyle_2::IndexStyle) where T)
            ex = :(f(Tuseless::Type{T}, ::IndexStyle, x::Int) where T)
            @test Rebugger.signature_names!(ex) == (:f, (:Tuseless, :__IndexStyle_1, :x), (), (:T,))
            @test ex == :(f(Tuseless::Type{T}, __IndexStyle_1::IndexStyle, x::Int) where T)
            # issue #34
            ex = :(_mapreduce_dim(f, op, ::NamedTuple{()}, A::AbstractArray, ::Colon))
            @test Rebugger.signature_names!(ex) == (:_mapreduce_dim, (:f, :op, :__NamedTuple_1, :A, :__Colon_1), (), ())
            @test ex == :(_mapreduce_dim(f, op, __NamedTuple_1::NamedTuple{()}, A::AbstractArray, __Colon_1::Colon))
        end
        @testset "Caller buffer capture and insertion" begin
            function run_insertion(str, atstr)
                RebuggerTesting.cbdata1[] = RebuggerTesting.cbdata2[] = Rebugger.stashed[] = nothing
                io = IOBuffer()
                idx = findfirst(atstr, str)
                print(io, str)
                seek(io, first(idx)-1)
                callexpr = Rebugger.prepare_caller_capture!(io)
                capstring = String(take!(io))
                capexpr   = Meta.parse(capstring)
                try
                    Core.eval(RebuggerTesting, capexpr)
                catch err
                    isa(err, StopException) || rethrow(err)
                end
            end

            str = """
            for i = 1:5
                cbdata1[] = i
                foo(12, 13; akw="modified")
                cbdata2[] = i
            end
            """
            @test run_insertion(str, "foo")
            @test RebuggerTesting.cbdata1[] == 1
            @test RebuggerTesting.cbdata2[] == nothing
            @test Rebugger.stashed[] == (RebuggerTesting.foo, (12, 13), Rebugger.kwstasher(;akw="modified"))
            str = """
            for i = 1:5
                error("not caught")
                foo(12, 13; akw="modified")
            end
            """
            @test_throws ErrorException("not caught") run_insertion(str, "foo")
            @test_throws Rebugger.StepException("Rebugger can only step into expressions, got 77") run_insertion("x = 77", "77")

            # Module-scoped calls
            io = IOBuffer()
            cmdstr = "Scope.func(x, y, z)"
            print(io, cmdstr)
            seek(io, 0)
            callexpr = Rebugger.prepare_caller_capture!(io)
            @test callexpr == :(Scope.func(x, y, z))
            take!(io)

            # getindex and setindex! expressions
            cmdstr = "x = a[2,3]"
            print(io, cmdstr)
            seek(io, first(findfirst("a", cmdstr))-1)
            callexpr = Rebugger.prepare_caller_capture!(io)
            @test callexpr == :(getindex(a, 2, 3))
            take!(io)

            cmdstr = "a[2,3] = x"
            print(io, cmdstr)
            seek(io, 0)
            callexpr = Rebugger.prepare_caller_capture!(io)
            @test callexpr == :(setindex!(a, x, 2, 3))
            take!(io)

            # Expressions that go beyond "user intention".
            # More generally we should support marking, but in the case of && and || it's
            # handled by lowering, so there is nothing to step into anyway.
            for cmdstr in ("f1(x) && f2(z)", "f1(x) || f2(z)")
                print(io, cmdstr)
                seek(io, 0)
                callexpr = Rebugger.prepare_caller_capture!(io)
                @test callexpr == :(f1(x))
                take!(io)
            end

            # issue #5
            cmdstr = "abs(abs(x))"
            print(io, cmdstr)
            seek(io, 4)
            callexpr = Rebugger.prepare_caller_capture!(io)
            @test callexpr == :(abs(x))
            take!(io)

            # splat expressions
            cmdstr = "foo(bar(x)..., 1)"
            print(io, cmdstr)
            idx = findfirst("bar", cmdstr)
            seek(io, first(idx)-1)
            callexpr = Rebugger.prepare_caller_capture!(io)
            @test callexpr == :(bar(x))
        end

        @testset "Callee variable capture" begin
            def = quote
                function complexargs(x::A, y=1, str="1.0"; kw1=Float64, kw2=7, kwargs...) where A<:AbstractArray{T} where T
                    return (x .+ y, parse(kw1, str), kw2)
                end
            end
            f = Core.eval(RebuggerTesting, def)
            @test f([8,9]) == ([9,10], 1.0, 7)
            m = collect(methods(f))[end]
            uuid = Rebugger.method_capture_from_callee(m, def)
            @test  Rebugger.method_capture_from_callee(m, def) == uuid  # calling twice returns the previously-defined objects
            fc = Rebugger.storefunc[uuid]
            @test_throws StopException fc([8,9], 2, "13"; kw1=Int, kw2=0)
            @test Rebugger.stored[uuid].varnames == (:x, :y, :str, :kw1, :kw2, :kwargs, :A, :T)
            @test Rebugger.stored[uuid].varvals  == ([8,9], 2, "13", Int, 0, empty_kwvarargs, Vector{Int}, Int)
            @test_throws StopException fc([8,9]; otherkw=77)
            @test Rebugger.stored[uuid].varnames == (:x, :y, :str, :kw1, :kw2, :kwargs, :A, :T)
            @test Rebugger.stored[uuid].varvals  == ([8,9], 1, "1.0", Float64, 7, pairs((otherkw=77,)), Vector{Int}, Int)

            uuid2 = Rebugger.method_capture_from_callee(m, def; overwrite=true)
            @test uuid2 != uuid
            # note overwriting methods are not stored in storefunc, but our old `f` will call the new method
            @test f([8,9], 2, "13"; kw1=Int, kw2=0) == ([10,11], 13, 0)
            Core.eval(RebuggerTesting, def)
            @test Rebugger.stored[uuid2].varnames == (:x, :y, :str, :kw1, :kw2, :kwargs, :A, :T)
            @test Rebugger.stored[uuid2].varvals  == ([8,9], 2, "13", Int, 0, empty_kwvarargs, Vector{Int}, Int)

            def = quote
                @inline modifies!(x) = (x[1] += 1; x)
            end
            f = Core.eval(RebuggerTesting, def)
            @test f([8,9]) == [9,9]
            m = collect(methods(f))[end]
            uuid = Rebugger.method_capture_from_callee(m, def)
            fc = Rebugger.storefunc[uuid]
            @test_throws StopException fc([8,9])
            @test Rebugger.stored[uuid].varnames == (:x,)
            @test Rebugger.stored[uuid].varvals  == ([8,9],)

            # Extensions of functions from other modules
            m = @which RebuggerTesting.foo()
            uuid = Rebugger.method_capture_from_callee(m)
            fc = Rebugger.storefunc[uuid]
            @test_throws StopException fc()
            @test Rebugger.stored[uuid].varnames == Rebugger.stored[uuid].varvals == ()
        end

        @testset "Step in" begin
            function run_stepin(str, atstr)
                io = IOBuffer()
                idx = findfirst(atstr, str)
                @test !isempty(idx)
                print(io, str)
                seek(io, first(idx)-1)
                Rebugger.stepin(io)
            end

            str = "RebuggerTesting.snoop0()"
            uuidref, cmd = run_stepin(str, str)
            uuid1 = uuidextractor(cmd)
            @test uuid1 == uuidref
            @test cmd == """
            @eval Main.RebuggerTesting let () = Main.Rebugger.getstored("$uuid1")
            begin
                snoop1("Spy")
            end
            end"""
            _, cmd = run_stepin(cmd, "snoop1")
            uuid2 = uuidextractor(cmd)
            @test cmd == """
            @eval Main.RebuggerTesting let (word,) = Main.Rebugger.getstored("$uuid2")
            begin
                snoop2(word, "on")
            end
            end"""
            _, cmd = run_stepin(cmd, "snoop2")
            uuid3 = uuidextractor(cmd)
            @test cmd == """
            @eval Main.RebuggerTesting let (word1, word2) = Main.Rebugger.getstored("$uuid3")
            begin
                snoop3(word1, word2, "arguments")
            end
            end"""
            @test Rebugger.getstored(string(uuid1)) == ()
            @test Rebugger.getstored(string(uuid2)) == ("Spy",)
            @test Rebugger.getstored(string(uuid3)) == ("Spy", "on")

            str = "RebuggerTesting.kwvarargs(1)"
            _, cmd = run_stepin(str, str)
            uuid = uuidextractor(cmd)
            @test cmd == """
            @eval Main.RebuggerTesting let (x, kw1, kwargs) = Main.Rebugger.getstored("$uuid")
            begin
                kwvarargs2(x; kw1=kw1, kwargs...)
            end
            end"""
            @test Rebugger.getstored(string(uuid)) == (1, 1, empty_kwvarargs)
            cmd = run_stepin(cmd, "kwvarargs2")

            str = "RebuggerTesting.kwvarargs(1; passthrough=false)"
            _, cmd = run_stepin(str, str)
            uuid = uuidextractor(cmd)
            @test cmd == """
            @eval Main.RebuggerTesting let (x, kw1, kwargs) = Main.Rebugger.getstored("$uuid")
            begin
                kwvarargs2(x; kw1=kw1, kwargs...)
            end
            end"""
            @test Rebugger.getstored(string(uuid)) == (1, 1, pairs((passthrough=false,)))
            _, cmd = run_stepin(cmd, "kwvarargs2")

            # Step in to call-overloading methods
            str = "RebuggerTesting.hv_test(\"hi\")"
            _, cmd = run_stepin(str, str)
            uuid = uuidextractor(cmd)
            @test cmd == """
            @eval Main.RebuggerTesting let (hv, str) = Main.Rebugger.getstored("$uuid")
            begin
                hv.x
            end
            end"""
            @test Rebugger.getstored(string(uuid)) == (RebuggerTesting.hv_test, "hi")

            # Step in to methods that do tuple-destructuring of arguments
            str = "RebuggerTesting.destruct(1, (2,3), 4)"
            @test eval(Meta.parse(str)) == 2
            _, cmd = run_stepin(str, str)
            uuid = uuidextractor(cmd)
            @test cmd == """
            @eval Main.RebuggerTesting let (x, (a, b), y) = Main.Rebugger.getstored("$uuid")
            begin
                a
            end
            end"""
            @test Rebugger.getstored(string(uuid)) == (1, (2,3), 4)

            # Step in to a broadcast call
            str = "sum.([[1,2], (3,5)])"
            uuid, cmd = run_stepin(str, str)
            s = Rebugger.stored[uuid]
            @test s.method.name == :broadcast
            @test cmd == """
            @eval Base.Broadcast let (f, As, Tf) = Main.Rebugger.getstored("$uuid")
            begin
                materialize(broadcasted(f, As...))
            end
            end"""
            @test Rebugger.getstored(string(uuid)) == (sum, (Any[[1,2], (3,5)],), typeof(sum))
            Core.eval(Main, Meta.parse(cmd)) == [3,8]

            str = "max.([1,5], [2,-3])"
            uuid, cmd = run_stepin(str, str)
            s = Rebugger.stored[uuid]
            @test s.method.name == :broadcast
            @test cmd == """
            @eval Base.Broadcast let (f, As, Tf) = Main.Rebugger.getstored("$uuid")
            begin
                materialize(broadcasted(f, As...))
            end
            end"""
            @test Rebugger.getstored(string(uuid)) == (max, ([1,5], [2,-3]), typeof(max))
            Core.eval(Main, Meta.parse(cmd)) == [2,5]

            # Step in to a do block
            str = "RebuggerTesting.calldo()"
            uuidref, cmd = run_stepin(str, str)
            uuid1 = uuidextractor(cmd)
            @test uuid1 == uuidref
            @test cmd == """
            @eval Main.RebuggerTesting let () = Main.Rebugger.getstored("$uuid1")
            begin
                apply(2, 3, 4) do x, y, z
                    snoop3(x, y, z)
                end
            end
            end"""
            uuidref, cmd = run_stepin(cmd, "apply")
            uuid1 = uuidextractor(cmd)
            @test uuid1 == uuidref
            @test cmd == """
            @eval Main.RebuggerTesting let (f, args) = Main.Rebugger.getstored("$uuid1")
            begin
                kwvarargs(f)
                f(args...)
            end
            end"""
        end

        @testset "Capture stacktrace" begin
            uuids = nothing
            mktemp() do path, iostacktrace
                redirect_stderr(iostacktrace) do
                    uuids = Rebugger.capture_stacktrace(RebuggerTesting, :(snoop0()))
                end
                flush(iostacktrace)
                str = read(path, String)
                @test occursin("snoop3", str)
            end
            @test Rebugger.stored[uuids[1]].varvals == ()
            @test Rebugger.stored[uuids[2]].varvals == ("Spy",)
            @test Rebugger.stored[uuids[3]].varvals == ("Spy", "on")
            @test Rebugger.stored[uuids[4]].varvals == ("Spy", "on", "arguments", "simply", empty_kwvarargs, String)
            @test_throws ErrorException("oops") RebuggerTesting.snoop0()

            st = try RebuggerTesting.kwfunctop(3) catch; stacktrace(catch_backtrace()) end
            usrtrace, defs = Rebugger.pregenerated_stacktrace(st; topname=Symbol("macro expansion"))
            @test length(unique(usrtrace)) == length(usrtrace)
            m = @which RebuggerTesting.kwfuncmiddle(1,1)
            @test m ∈ usrtrace

            # A case that tests inlining and several other aspects of argument capture
            ex = :([1, 2, 3] .* [1, 2])
            # Capture the actual stack trace, trimming it to avoid
            # anything involving the `eval` itself
            trace = try
                Core.eval(Main, ex)
            catch
                stacktrace(catch_backtrace())
            end
            i = 1
            while i <= length(trace)
                t = trace[i]
                if t.func == Symbol("top-level scope")
                    deleteat!(trace, i:length(trace))
                end
                i += 1
            end
            # Get the capture from Rebugger
            uuids = mktemp() do path, iostacktrace
                redirect_stderr(iostacktrace) do
                    Rebugger.capture_stacktrace(Main, ex)
                end
            end
            @test length(uuids) == length(trace)
            for (uuid, t) in zip(reverse(uuids), trace)
                @test Rebugger.stored[uuid].method.name == t.func
            end

            # Try capturing a method from Core. On binaries this would throw
            # if we didn't catch it.
            # Because the first entry is "top-level scope", and that terminates
            # processing in Rebugger.pregenerated_stacktrace, we have to intervene a bit.
            mod, ex = Main, :(Core.throw(ArgumentError("oops")))
            trace = try Core.eval(mod, command) catch err stacktrace(catch_backtrace()) end
            usrtrace, defs = Rebugger.pregenerated_stacktrace(trace[2:3])
            @test usrtrace isa Vector
        end
    end

    @testset "User interface" begin
        @testset "Printing header" begin
            h = Rebugger.RebugHeader()
            h.uuid = uuid = uuid1()
            meth = @which RebuggerTesting.foo(1,2)
            h.current_method = meth
            Rebugger.stored[uuid] = Rebugger.Stored(meth, (:x, :y), (1, ErrorsOnShow()))
            h.warnmsg = "This is a warning"
            h.errmsg  = "You will not have a second chance"
            io = IOBuffer()
            Rebugger.print_header(io, h)
            str = String(take!(io))
            @test startswith(str, """
            This is a warning
            You will not have a second chance
            foo(x, y) in Main.RebuggerTesting at """) # skip the "upper" part of the file location
            @test endswith(str, "testmodule.jl:7\n  x = 1\n  y errors in its show method")
        end

        @testset "Demos" begin
            function prepare_step_command(cmd, atstr)
                LineEdit.edit_clear(mistate)
                idx = findfirst(atstr, cmd)
                @test !isempty(idx)
                LineEdit.replace_line(mistate, cmd)
                buf = LineEdit.buffer(mistate)
                seek(buf, first(idx)-1)
                return mistate
            end

            function do_capture_stacktrace(cmd)
                l = length(hist.history)
                LineEdit.replace_line(mistate, cmd)
                Rebugger.capture_stacktrace(mistate)
                LineEdit.transition(mistate, julia_prompt)
                return l+1:length(hist.history)
            end

            if isdefined(Base, :active_repl)
                repl = Base.active_repl
                mistate = repl.mistate
                julia_prompt = find_prompt(mistate, "julia")
                LineEdit.transition(mistate, julia_prompt)
                hist = julia_prompt.hist
                header = Rebugger.rebug_prompt_ref[].repl.header
                histdel = 0

                @testset "show demo" begin  # this is a demo that appears in the documentation
                    cmd1 = "show([1,2,4])"
                    s = prepare_step_command(cmd1, cmd1)
                    Rebugger.stepin(s)
                    histdel += 1
                    uuid = header.uuid
                    @test Rebugger.getstored(string(uuid)) == ([1,2,4],)
                    cmd2 = LineEdit.content(s)
                    s = prepare_step_command(cmd2, "show(stdout::IO, x)")
                    Rebugger.stepin(s)
                    histdel += 1
                    uuid = header.uuid
                    @test Rebugger.getstored(string(uuid))[2] == [1,2,4]
                    cmd3 = LineEdit.content(s)
                    s = prepare_step_command(cmd3, "_show_empty")
                    Rebugger.stepin(s)
                    histdel += 1
                    @test header.warnmsg == "Execution did not reach point"
                end

                @testset "Colors demo" begin  # another demo that appears in the documentation
                    desc = "hsl(80%, 20%, 15%)"
                    cmd = "colorant\"hsl(80%, 20%, 15%)\""
                    local idx
                    mktemp() do path, io
                        redirect_stderr(io) do
                            logs, _ = Test.collect_test_logs() do
                                idx = do_capture_stacktrace(cmd)
                            end
                        end
                        flush(io)
                        seek(io, 0)
                        @test countlines(io) >= 4
                    end
                    histdel += length(idx)
                    @test length(idx) >= 5
                    @test hist.history[idx[1]] == cmd
                    @test occursin("error", hist.history[idx[end]])
                end

                @testset "Pkg demo" begin
                    updated = Pkg.UPDATED_REGISTRY_THIS_SESSION[]
                    Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true
                    uuids = Rebugger.capture_stacktrace(Pkg, :(add("NoPkg")))
                    @test length(uuids) >= 2
                    Pkg.UPDATED_REGISTRY_THIS_SESSION[] = updated
                end

                @testset "Empty stacktraces" begin
                    cmd = "ccall(:jl_throw, Nothing, (Any,), ArgumentError(\"oops\"))"
                    mktemp() do path, io
                        redirect_stderr(io) do
                            LineEdit.replace_line(mistate, cmd)
                            @test Rebugger.capture_stacktrace(mistate) === nothing
                            LineEdit.transition(mistate, julia_prompt)
                        end
                        flush(io)
                        str = read(path, String)
                        @test occursin("failed to capture", str)
                    end
                end

                LineEdit.edit_clear(mistate)
                l = length(hist.history)
                deleteat!(hist.history, l-histdel+1:l)
                deleteat!(hist.modes, l-histdel+1:l)
                hist.cur_idx = length(hist.history)+1
            end
        end
    end
end


text:
# Invertible HINT coupling layer from Kruse et al. (2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export CouplingLayerHINT, CouplingLayerHINT3D

"""
    H = CouplingLayerHINT(n_in, n_hidden; logdet=false, permute="none", k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, ndims=2) (2D)

    H = CouplingLayerHINT(n_in, n_hidden; logdet=false, permute="none", k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, ndims=3) (3D)

    H = CouplingLayerHINT3D(n_in, n_hidden; logdet=false, permute="none", k1=3, k2=3, p1=1, p2=1, s1=1, s2=1) (3D)

 Create a recursive HINT-style invertible layer based on coupling blocks.

 *Input*:

 - `n_in`, `n_hidden`: number of input and hidden channels

 - `logdet`: bool to indicate whether to return the log determinant. Default is `false`.

 - `permute`: string to specify permutation. Options are `"none"`, `"lower"`, `"both"` or `"full"`.

 - `k1`, `k2`: kernel size of convolutions in residual block. `k1` is the kernel of the first and third
    operator, `k2` is the kernel size of the second operator.

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`)

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`)

 - `ndims` : number of dimensions

 *Output*:

 - `H`: Recursive invertible HINT coupling layer.

 *Usage:*

 - Forward mode: `Y = H.forward(X)`

 - Inverse mode: `X = H.inverse(Y)`

 - Backward mode: `ΔX, X = H.backward(ΔY, Y)`

 *Trainable parameters:*

 - None in `H` itself

 - Trainable parameters in coupling layers `H.CL`

 See also: [`CouplingLayerBasic`](@ref), [`ResidualBlock`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
mutable struct CouplingLayerHINT <: NeuralNetLayer
    CL::AbstractArray{CouplingLayerBasic, 1}
    C::Union{Conv1x1, Nothing}
    logdet::Bool
    permute::String
    is_reversed::Bool
end

@Flux.functor CouplingLayerHINT

# Get layer depth for recursion
function get_depth(n_in)
    count = 0
    nc = n_in
    while nc > 4
        nc /= 2
        count += 1
    end
    return count +1
end

# Constructor for given coupling layer and 1 x 1 convolution
CouplingLayerHINT(CL::AbstractArray{CouplingLayerBasic, 1}, C::Union{Conv1x1, Nothing};
    logdet=false, permute="none") = CouplingLayerHINT(CL, C, logdet, permute, false)

# 2D Constructor from input dimensions
function CouplingLayerHINT(n_in::Int64, n_hidden::Int64; logdet=false, permute="none",
                           k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, ndims=2)

    # Create basic coupling layers
    n = get_depth(n_in)
    CL = Array{CouplingLayerBasic}(undef, n)
    for j=1:n
        CL[j] = CouplingLayerBasic(Int(n_in/2^j), n_hidden; k1=k1, k2=k2, p1=p1, p2=p2,
                                   s1=s1, s2=s2, logdet=logdet, ndims=ndims)
    end

    # Permutation using 1x1 convolution
    if permute == "full" || permute == "both"
        C = Conv1x1(n_in)
    elseif permute == "lower"
        C = Conv1x1(Int(n_in/2))
    else
        C = nothing
    end

    return CouplingLayerHINT(CL, C, logdet, permute, false)
end

CouplingLayerHINT3D(args...;kw...) = CouplingLayerHINT(args...; kw..., ndims=3)

# Input is tensor X
function forward(X, H::CouplingLayerHINT; scale=1, permute=nothing, logdet=nothing)
    isnothing(logdet) ? logdet = (H.logdet && ~H.is_reversed) : logdet = logdet
    isnothing(permute) ? permute = H.permute : permute = permute

    # Permutation
    if permute == "full" || permute == "both"
        X = H.C.forward(X)
    end
    Xa, Xb = tensor_split(X)
    permute == "lower" && (Xb = H.C.forward(Xb))

    # Determine whether to continue recursion
    recursive = false
    if typeof(X) <: AbstractArray{Float32, 4} && size(X, 3) > 4
        recursive = true
    elseif typeof(X) <: AbstractArray{Float32, 5} && size(X, 4) > 4
        recursive = true
    end

    # HINT coupling
    if recursive
        # Call function recursively
        Ya, logdet1 = forward(Xa, H; scale=scale+1, permute="none")
        Y_temp, logdet2 = forward(Xb, H; scale=scale+1, permute="none")
        if logdet
            Yb, logdet3 = H.CL[scale].forward(Xa, Y_temp)[[2,3]]
        else
            Yb = H.CL[scale].forward(Xa, Y_temp)[2]
            logdet3 = 0f0
        end
        logdet_full = logdet1 + logdet2 + logdet3
    else
        # Finest layer
        Ya = copy(Xa)
        if logdet
            Yb, logdet_full = H.CL[scale].forward(Xa, Xb)[[2,3]]
        else
            Yb = H.CL[scale].forward(Xa, Xb)[2]
            logdet_full = 0f0
        end
    end

    Y = tensor_cat(Ya, Yb)
    permute == "both" && (Y = H.C.inverse(Y))
    if scale == 1
        logdet ? (return Y, logdet_full) : (return Y)
    else
        return Y, logdet_full
    end
end

# Input is tensor Y
function inverse(Y::AbstractArray{Float32, N} , H::CouplingLayerHINT; scale=1, permute=nothing, logdet=nothing) where N
    isnothing(logdet) ? logdet = (H.logdet && H.is_reversed) : logdet = logdet
    isnothing(permute) ? permute = H.permute : permute = permute

    # Permutation
    permute == "both" && (Y = H.C.forward(Y))
    Ya, Yb = tensor_split(Y)

    # Check for recursion
    recursive = (size(Y, N-1) > 4)

    # Coupling layer
    if recursive
        Xa, logdet1 = inverse(Ya, H; scale=scale+1, permute="none")
        if logdet
            Y_temp, logdet2 = H.CL[scale].inverse(Xa, Yb; logdet=true)[[2,3]]
        else
            Y_temp = H.CL[scale].inverse(Xa, Yb)[2]
            logdet2 = 0f0
        end
        Xb, logdet3 = inverse(Y_temp, H; scale=scale+1, permute="none")
        logdet_full = logdet1 + logdet2 + logdet3
    else
        Xa = copy(Ya)
        if logdet
            Xb, logdet_full = H.CL[scale].inverse(Ya, Yb)[[2,3]]
        else
            Xb = H.CL[scale].inverse(Ya, Yb)[2]
            logdet_full = 0f0
        end
    end

    # Initial permutation
    permute == "lower" && (Xb = H.C.inverse(Xb))
    X = tensor_cat(Xa, Xb)
    if permute == "full" || permute == "both"
        X = H.C.inverse(X)
    end

    if scale == 1
        logdet ? (return X, logdet_full) : (return X)
    else
        return X, logdet_full
    end
end

# Input are two tensors ΔY, Y
function backward(ΔY, Y::AbstractArray{Float32, N}, H::CouplingLayerHINT; scale=1, permute=nothing, set_grad::Bool=true) where N
    isnothing(permute) ? permute = H.permute : permute = permute

    # Initializing output parameter array
    if !set_grad
        nscales_tot = length(H.CL)
        nparams = 5*(nscales_tot-scale+1)
        (permute != "none") && (nparams += 3)
        Δθ = Array{Parameter, 1}(undef, nparams)
        H.logdet && (∇logdet = Array{Parameter, 1}(undef, nparams))
    end

    if permute == "both"
        if set_grad
            ΔY, Y = H.C.forward((ΔY, Y))
        else
            ΔY, Δθ_C, Y = H.C.forward((ΔY, Y); set_grad=set_grad)
            Δθ[end-2:end] .= Δθ_C
            H.logdet && (∇logdet[end-2:end] .= [Parameter(zeros(Float32, size(H.C.v1))), Parameter(zeros(Float32, size(H.C.v2))), Parameter(zeros(Float32, size(H.C.v3)))])
        end
    end
    Ya, Yb = tensor_split(Y)
    ΔYa, ΔYb = tensor_split(ΔY)

    # Determine whether to continue recursion
    recursive = (size(Y, N-1) > 4)

    # HINT coupling
    if recursive
        if set_grad
            ΔXa, Xa = backward(ΔYa, Ya, H; scale=scale+1, permute="none")
            ΔXa_temp, ΔXb_temp, X_temp = H.CL[scale].backward(ΔXa.*0f0, ΔYb, Xa, Yb)[[1,2,4]]
            ΔXb, Xb = backward(ΔXb_temp, X_temp, H; scale=scale+1, permute="none")
        else
            if H.logdet
                ΔXa, Δθa, Xa, ∇logdet_a = backward(ΔYa, Ya, H; scale=scale+1, permute="none", set_grad=set_grad)
                ΔXa_temp, ΔXb_temp, Δθ_scale, _, X_temp, ∇logdet_scale = H.CL[scale].backward(ΔXa.*0f0, ΔYb, Xa, Yb; set_grad=set_grad)
                ΔXb, Δθb, Xb, ∇logdet_b = backward(ΔXb_temp, X_temp, H; scale=scale+1, permute="none", set_grad=set_grad)
                ∇logdet[1:5] .= ∇logdet_scale
                ∇logdet[6:5+length(∇logdet_a)] .= ∇logdet_a+∇logdet_b
            else
                ΔXa, Δθa, Xa = backward(ΔYa, Ya, H; scale=scale+1, permute="none", set_grad=set_grad)
                ΔXa_temp, ΔXb_temp, Δθ_scale, _, X_temp = H.CL[scale].backward(ΔXa.*0f0, ΔYb, Xa, Yb; set_grad=set_grad)
                ΔXb, Δθb, Xb = backward(ΔXb_temp, X_temp, H; scale=scale+1, permute="none", set_grad=set_grad)
            end
            Δθ[1:5] .= Δθ_scale
            Δθ[6:5+length(Δθa)] .= Δθa+Δθb
        end
        ΔXa += ΔXa_temp
    else
        Xa = copy(Ya)
        ΔXa = copy(ΔYa)
        if set_grad
            ΔXa_, ΔXb, Xb = H.CL[scale].backward(ΔYa.*0f0, ΔYb, Ya, Yb)[[1,2,4]]
        else
            if H.logdet
                ΔXa_, ΔXb, Δθ_scale, _, Xb, ∇logdet_scale = H.CL[scale].backward(ΔYa.*0f0, ΔYb, Ya, Yb; set_grad=set_grad)
                ∇logdet[1:5] .= ∇logdet_scale
            else
                ΔXa_, ΔXb, Δθ_scale, _, Xb = H.CL[scale].backward(ΔYa.*0f0, ΔYb, Ya, Yb; set_grad=set_grad)
            end
            Δθ[1:5] .= Δθ_scale
        end
        ΔXa += ΔXa_
    end
    if permute == "lower"
        if set_grad
            ΔXb, Xb = H.C.inverse((ΔXb, Xb))
        else
            ΔXb, Δθ_C, Xb = H.C.inverse((ΔXb, Xb); set_grad=set_grad)
            H.logdet && (∇logdet[end-2:end] .= [Parameter(zeros(Float32, size(H.C.v1))), Parameter(zeros(Float32, size(H.C.v2))), Parameter(zeros(Float32, size(H.C.v3)))])
        end
    end
    ΔX = tensor_cat(ΔXa, ΔXb)
    X = tensor_cat(Xa, Xb)
    if permute == "full" || permute == "both"
        if set_grad
            ΔX, X = H.C.inverse((ΔX, X))
        else
            ΔX, Δθ_C, X = H.C.inverse((ΔX, X); set_grad=set_grad)
            if permute == "full"
                Δθ[end-2:end] .= Δθ_C
                H.logdet && (∇logdet[end-2:end] .= [Parameter(zeros(Float32, size(H.C.v1))), Parameter(zeros(Float32, size(H.C.v2))), Parameter(zeros(Float32, size(H.C.v3)))])
            else
                Δθ[end-2:end] += Δθ_C
                H.logdet && (∇logdet[end-2:end] += [Parameter(zeros(Float32, size(H.C.v1))), Parameter(zeros(Float32, size(H.C.v2))), Parameter(zeros(Float32, size(H.C.v3)))])
            end
        end
    end

    if set_grad
        return ΔX, X
    else
        H.logdet ? (return ΔX, Δθ, X, ∇logdet) : (return ΔX, Δθ, X)
    end
end

# Input are two tensors ΔX, X
function backward_inv(ΔX, X::AbstractArray{Float32, N}, H::CouplingLayerHINT; scale=1, permute=nothing) where N
    isnothing(permute) ? permute = H.permute : permute = permute

    # Permutation
    if permute == "full" || permute == "both"
        (ΔX, X) = H.C.forward((ΔX, X))
    end
    ΔXa, ΔXb = tensor_split(ΔX)
    Xa, Xb = tensor_split(X)
    permute == "lower" && ((ΔXb, Xb) = H.C.forward((ΔXb, Xb)))

    # Check whether to continue recursion
    recursive = (size(X, N-1) > 4)

    # Coupling layer backprop
    if recursive
        ΔY_temp, Y_temp = backward_inv(ΔXb, Xb, H; scale=scale+1, permute="none")
        ΔYa_temp, ΔYb, Yb = backward_inv(0f0.*ΔXa, ΔY_temp, Xa, Y_temp, H.CL[scale])[[1,2,4]]
        ΔYa, Ya = backward_inv(ΔXa+ΔYa_temp, Xa, H; scale=scale+1, permute="none")
    else
        ΔYa = copy(ΔXa)
        Ya = copy(Xa)
        ΔYa_temp, ΔYb, Yb = backward_inv(0f0.*ΔYa, ΔXb, Xa, Xb, H.CL[scale])[[1,2,4]]
        ΔYa += ΔYa_temp
    end
    ΔY = tensor_cat(ΔYa, ΔYb)
    Y = tensor_cat(Ya, Yb)
    permute == "both" && ((ΔY, Y) = H.C.inverse((ΔY, Y)))
    return ΔY, Y
end

## Jacobian-related functions

function jacobian(ΔX, Δθ::Array{Parameter, 1}, X::AbstractArray{Float32, N}, H::CouplingLayerHINT;
                  scale=1, permute=nothing, logdet=nothing) where N
    isnothing(logdet) ? logdet = (H.logdet && ~H.is_reversed) : logdet = logdet
    isnothing(permute) ? permute = H.permute : permute = permute

    # Selecting parameters
    if permute == "none"
        Δθscale = Δθ[1:5]
        Δθ = Δθ[6:end]
    else
        Δθscale = Δθ[1:5]
        Δθ_C = Δθ[end-2:end]
        Δθ = Δθ[6:end-3]
    end

    # Permutation
    if permute == "full" || permute == "both"
        ΔX, X = H.C.jacobian(ΔX, Δθ_C, X)
    end
    Xa, Xb = tensor_split(X)
    ΔXa, ΔXb = tensor_split(ΔX)
    if permute == "lower"
        ΔXb, Xb = H.C.jacobian(ΔXb, Δθ_C, Xb)
    end

    # Initialize Gauss-Newton approx of logdet term
    if logdet
        nscales_tot = length(H.CL)
        nparams = 5*(nscales_tot-scale+1)
        (permute != "none") && (nparams += 3)
        GNΔθ_full = Array{Parameter, 1}(undef, nparams)
        (permute != "none") && (GNΔθ_full[end-2:end] = [Parameter(zeros(Float32, size(H.C.v1))), Parameter(zeros(Float32, size(H.C.v2))), Parameter(zeros(Float32, size(H.C.v3)))])
    end

    # Determine whether to continue recursion
    recursive = (size(X, N-1) > 4)

    # HINT coupling
    # idx_Δθ_scale = (scale-1)*5+1:scale*5
    if recursive
        # Call function recursively
        if logdet
            ΔYa, Ya, logdet1, GNΔθ1 = jacobian(ΔXa, Δθ, Xa, H; scale=scale+1, permute="none")
            ΔY_temp, Y_temp, logdet2, GNΔθ2 = jacobian(ΔXb, Δθ, Xb, H; scale=scale+1, permute="none")
            _, ΔYb, _, Yb, logdet3, GNΔθ3 = H.CL[scale].jacobian(ΔXa, ΔY_temp, Δθscale, Xa, Y_temp)
            logdet_full = logdet1 + logdet2 + logdet3
            GNΔθ_full[1:5] .= GNΔθ3
            (permute != "none") ? (GNΔθ_full[6:end-3] .= GNΔθ1 + GNΔθ2) : (GNΔθ_full[6:end] .= GNΔθ1 + GNΔθ2)
        else
            ΔYa, Ya = jacobian(ΔXa, Δθ, Xa, H; scale=scale+1, permute="none")
            ΔY_temp, Y_temp = jacobian(ΔXb, Δθ, Xb, H; scale=scale+1, permute="none")
            _, ΔYb, _, Yb = H.CL[scale].jacobian(ΔXa, ΔY_temp, Δθscale, Xa, Y_temp)
        end
    else
        # Coarsest scale
        Ya = copy(Xa)
        ΔYa = copy(ΔXa)
        if logdet
            _, ΔYb, _, Yb, logdet_full, GNΔθ3 = H.CL[scale].jacobian(ΔXa, ΔXb, Δθscale, Xa, Xb)
            GNΔθ_full[1:5] .= GNΔθ3
        else
            _, ΔYb, _, Yb = H.CL[scale].jacobian(ΔXa, ΔXb, Δθscale, Xa, Xb)
        end
    end

    Y = tensor_cat(Ya, Yb)
    ΔY = tensor_cat(ΔYa, ΔYb)
    if permute == "both"
        ΔY, Y = H.C.jacobianInverse(ΔY, Δθ_C, Y)
    end
    logdet ? (return ΔY, Y, logdet_full, GNΔθ_full) : (return ΔY, Y)

end

adjointJacobian(ΔY, Y, H::CouplingLayerHINT; scale=1, permute=nothing) = backward(ΔY, Y, H; scale=scale, permute=permute, set_grad=false)

## Other utils

# Clear gradients
function clear_grad!(H::CouplingLayerHINT)
    for j=1:length(H.CL)
        clear_grad!(H.CL[j])
    end
    ~isnothing(H.C) && clear_grad!(H.C)
end

# Get parameters
function get_params(H::CouplingLayerHINT)
    nlayers = length(H.CL)
    p = get_params(H.CL[1])
    if nlayers > 1
        for j=2:nlayers
            p = cat(p, get_params(H.CL[j]); dims=1)
        end
    end
    ~isnothing(H.C) && (p = cat(p, get_params(H.C); dims=1))
    return p
end

# Set is_reversed flag in full network tree
function tag_as_reversed!(H::CouplingLayerHINT, tag::Bool)
    H.is_reversed = tag
    nlayers = length(H.CL)
    for j=1:nlayers
        H.CL[j].is_reversed = tag
    end
    return H
end


text:
module EmbeddedFunctions

#--- Dependencies

import Base: resize!

#--- Exports

# Exported types
export StateVar, ZeroFunction, MonitorFunction, ClosedEmbeddedFunctions

# Exported functions
export getu, getu!, setu!, getmu, getmu!, setmu!, getvars, getvars!,
    setvars!, muidx, active, active!, dim_u, dim_mu, dim_phi, dim_psi,
    mu_idx, mu_name, eval_efuncs!, eval_mfuncs!

#--- State variables

struct StateVar{T <: Number}
    name::String  # name of the state variable
    u::Vector{T}  # initial state (may be changed during continuation)
    t::Vector{T}  # initial tangent (may be changed during continuation)
end
StateVar(name::String, u::Vector{T}) where {T <: Number} = StateVar(name, u, Vector{T}())

#--- Zero functions (otherwise known as zero problems)

struct ZeroFunction{T <: Number, F}
    f::F  # underlying function
    u::Vector{StateVar{T}}  # underlying state variables
    res::Vector{T}  # output vector
end
ZeroFunction(f, u::Vector{StateVar{T}}, m::Integer) where T <: Number = ZeroFunction(f, u, Vector{T}(undef, m))

#--- Monitor functions

struct MonitorFunction{T <: Number, F}
    f::F  #underlying function
    u::Vector{StateVar{T}}  # underlying state variables
    μ_name::String  # name of the continuation parameter
    active::Bool  # whether the continuation parameter is active initially
end

#--- Closed functions (for computational speed)

const SimpleView{T} = SubArray{T, 1, Array{T, 1}, Tuple{UnitRange{Int}}, true} where T

struct ClosedEmbeddedFunctions{T <: Number,
                               F <: Tuple,
                               FU <: Tuple,
                               G <: Tuple,
                               GU <: Tuple}
    u::Vector{StateVar{T}}
    uᵢ::Vector{Tuple{Int, Int}}
    uᵥ::Vector{SimpleView{T}}  # view cache; should be same size as u
    Φ::F
    Φᵤ::FU
    Φᵢ::Vector{Tuple{Int, Int}}
    μ::Vector{T}
    μᵢ::Vector{Int}
    μs::Dict{String, Int}
    𝕁::Vector{Bool}  # 𝕁 denotes active continuation parameters (i.e., continuation parameters that vary)
    Ψ::G
    Ψᵤ::GU
    Ψᵢ::Vector{Int}
end

"""
    ClosedEmbeddedFunctions(Φ::Vector{ZeroFunction}, Ψ::Vector{MonitorFunction})

A ClosedEmbeddedFunctions structure is an optimised structure that allows fast
execution of the embedded zero and monitor functions within the problem,
avoiding dynamic dispatch.

This function generates a ClosedEmbeddedFunctions structure from a vector of
zero functions and a vector of monitor functions. Any state variables referenced
will be automatically included in the closed problem.
"""
function ClosedEmbeddedFunctions(Φ::Vector{<: ZeroFunction{T}}, Ψ::Vector{<: MonitorFunction{T}}) where T <: Number
    # Check for uniqueness
    if !allunique(Φ)
        throw(ArgumentError("Some zero functions are included multiple times"))
    end
    if !allunique(Ψ)
        throw(ArgumentError("Some monitor functions are included multiple times"))
    end
    # Get the necessary state variables
    u = Vector{StateVar{T}}()
    Φᵤ = []  # a vector of tuples of different lengths
    for ϕ in Φ
        ϕᵤ = Int[]
        for uu in ϕ.u
            if !(uu in u)
                push!(u, uu)
            end
            push!(ϕᵤ, findfirst(isequal(uu), u))
        end
        push!(Φᵤ, (ϕᵤ...,))
    end
    Φᵢ = Vector{Tuple{Int, Int}}(undef, length(Φ))
    Ψᵤ = []  # a vector of tuples of different lengths so needs to be Any
    𝕁 = Bool[]
    μs = Dict{String, Int}()
    for i = eachindex(Ψ)
        ψ = Ψ[i]
        ψᵤ = Int[]
        for uu in ψ.u
            if !(uu in u)
                push!(u, uu)
            end
            push!(ψᵤ, findfirst(isequal(uu), u))
        end
        push!(𝕁, ψ.active)
        push!(Ψᵤ, (ψᵤ...,))
        push!(μs, ψ.μ_name => i)
    end
    Ψᵢ = Vector{Int}(undef, length(Ψ))
    uᵢ = Vector{Tuple{Int, Int}}(undef, length(u))
    uᵥ = Vector{SimpleView{T}}(undef, length(u))
    μ = Vector{T}(undef, length(Ψ))
    μᵢ = Vector{Int}(undef, length(μ))
    closed = ClosedEmbeddedFunctions(u, uᵢ, uᵥ, (Φ...,), (Φᵤ...,), Φᵢ, μ, μᵢ, μs, 𝕁, (Ψ...,), (Ψᵤ...,), Ψᵢ)
    resize!(closed)
end

"""
    resize!(closed::ClosedEmbeddedFunctions)

Update the indices within `closed` to reflect changes in the sizes of the
underlying state variables and the dimensions of the zero functions.
"""
function resize!(closed::ClosedEmbeddedFunctions)
    # State variables
    idx = 1
    for i = eachindex(closed.u)
        n = length(closed.u[i].u)
        closed.uᵢ[i] = (idx, idx + n - 1)
        idx += n
    end
    # Continuation parameters
    for i = eachindex(closed.μ)
        if closed.𝕁[i]
            closed.μᵢ[i] = idx
            idx += 1
        else
            closed.μᵢ[i] = 0
        end
    end
    # Zero and monitor functions
    idx = 1
    for i = eachindex(closed.Φ)
        m = length(closed.Φ[i].res)
        closed.Φᵢ[i] = (idx, idx + m - 1)
        idx += m
    end
    for i = eachindex(closed.Ψ)
        closed.Ψᵢ[i] = idx
        idx += 1
    end
    closed
end

@generated function eval_efuncs!(
        res::AbstractVector{T},
        closed::ClosedEmbeddedFunctions{T, F, FU, G, GU},
        prob,
        u::AbstractVector{T}
        ) where {T <: Number, F, FU, G, GU}
    body = quote
        # Views on the state variables get reused, so precompute them
        for i = eachindex(closed.uᵢ)
            (i0, i1) = closed.uᵢ[i]
            closed.uᵥ[i] = view(u, i0:i1)
        end
        # Copy any active continuation parameter values into the μ variable
        for i = eachindex(closed.μ)
            if closed.𝕁[i]
                closed.μ[i] = u[closed.μᵢ[i]]
            end
        end
    end
    for i in eachindex(FU.parameters)
        # Construct function calls of the form Φ[i](resᵥ[i], uᵥ[Φᵤ[i][1]], ..., uᵥ[Φᵤ[i][n]])
        if length(FU.parameters[i].parameters) == 0
            # No dependencies means pass everything
            push!(body.args, :(closed.Φ[$i].f(view(res, closed.Φᵢ[$i][1]:closed.Φᵢ[$i][2]), prob, u)))
        else
            push!(body.args, :(closed.Φ[$i].f(view(res, closed.Φᵢ[$i][1]:closed.Φᵢ[$i][2]), prob, $((:(closed.uᵥ[closed.Φᵤ[$i][$j]]) for j in eachindex(FU.parameters[i].parameters))...))))
        end
    end
    for i in eachindex(GU.parameters)
        # Construct function calls of the form res[Ψᵢ[i]] = Ψ[i](uᵥ[Ψᵤ[i][1]], ..., uᵥ[Ψᵤ[i][n]]) - μ[i]
        # Uses the return value of Ψ in contrast to Φ since it is assumed to be ℝ rather than ℝⁿ
        if length(GU.parameters[i].parameters) == 0
            # No dependencies means pass everything
            push!(body.args, :(res[closed.Ψᵢ[$i]] = closed.Ψ[$i].f(prob, u) - closed.μ[$i]))
        else
            push!(body.args, :(res[closed.Ψᵢ[$i]] = closed.Ψ[$i].f(prob, $((:(closed.uᵥ[closed.Ψᵤ[$i][$j]]) for j in eachindex(GU.parameters[i].parameters))...)) - closed.μ[$i]))
        end
    end
    push!(body.args, :res)
    body
end

@generated function eval_mfuncs!(
        res::AbstractVector{T},
        closed::ClosedEmbeddedFunctions{T, F, FU, G, GU},
        prob,
        u::AbstractVector{T}
        ) where {T <: Number, F, FU, G, GU}
    body = quote
        # Views on the state variables get reused, so precompute them
        for i = eachindex(closed.uᵢ)
            (i0, i1) = closed.uᵢ[i]
            closed.uᵥ[i] = view(u, i0:i1)
        end
    end
    for i in eachindex(GU.parameters)
        # Construct function calls of the form res[i] = Ψ[i](uᵥ[Ψᵤ[i][1]], ..., uᵥ[Ψᵤ[i][n]])
        # Uses the return value of Ψ in contrast to Φ since it is assumed to be ℝ rather than ℝⁿ
        if length(GU.parameters[i].parameters) == 0
            # No dependencies means pass everything
            push!(body.args, :(res[$i] = closed.Ψ[$i].f(prob, u)))
        else
            push!(body.args, :(res[$i] = closed.Ψ[$i].f(prob, $((:(closed.uᵥ[closed.Ψᵤ[$i][$j]]) for j in eachindex(GU.parameters[i].parameters))...))))
        end
    end
    push!(body.args, :res)
    body
end

"""
    mu_idx(closed::ClosedEmbeddedFunctions, μ::String)

Return the index in the continuation parameter vector of the specified
continuation parameter `μ`.
"""
mu_idx(closed::ClosedEmbeddedFunctions, μ::String) = closed.μs[μ]
mu_idx(closed::ClosedEmbeddedFunctions, μ::Integer) = μ

"""
    mu_name(closed::ClosedEmbeddedFunctions, μ::Integer)

Return the name of the specified continuation parameter `μ`.
"""
mu_name(closed::ClosedEmbeddedFunctions, μ::Integer) = closed.Ψ[μ].μ_name

"""
    active(closed::ClosedEmbeddedFunctions)

Return the number of active continuation parameters.
"""
active(closed::ClosedEmbeddedFunctions) = sum(closed.𝕁)

"""
    active(closed::ClosedEmbeddedFunctions, μ)

Return whether the specified continuation parameter `μ` is active or not.
"""
active(closed::ClosedEmbeddedFunctions, μ) = closed.𝕁[muidx(closed, μ)]

"""
    active!(closed::ClosedEmbeddedFunctions, μ, active::Bool)

Set the specified continuation parameter `μ` to be active or not.
"""
active!(closed::ClosedEmbeddedFunctions, μ, active::Bool) = (closed.𝕁[muidx(closed, μ)] = active)

"""
    dim_u(closed::ClosedEmbeddedFunctions)

Return the total number of state variables.
"""
dim_u(closed::ClosedEmbeddedFunctions) = closed.uᵢ[end][end]

"""
    dim_mu(closed::ClosedEmbeddedFunctions)

Return the total number of continuation parameters.
"""
dim_mu(closed::ClosedEmbeddedFunctions) = length(closed.μ)

"""
    dim_phi(closed::ClosedEmbeddedFunctions)

Return the number of output dimensions of the set of zero functions.
"""
dim_phi(closed::ClosedEmbeddedFunctions) = closed.Φᵢ[end][end]

"""
    dim_psi(closed::ClosedEmbeddedFunctions)

Return the number of output dimensions of the set of monitor functions.
"""
dim_psi(closed::ClosedEmbeddedFunctions) = length(closed.Ψ)

function getu!(u::AbstractVector{T}, closed::ClosedEmbeddedFunctions{T}) where T <: Number
    for i = eachindex(closed.u)
        u[closed.uᵢ[i][1]:closed.uᵢ[i][2]] .= closed.u[i].u
    end
    u
end
getu(closed::ClosedEmbeddedFunctions{T}) where {T <: Number} = getu!(zeros(T, dim_u(closed)), closed)

function setu!(closed::ClosedEmbeddedFunctions{T}, u::AbstractVector{T}) where T <: Number
    for i = eachindex(closed.u)
        closed.u[i].u .= u[closed.uᵢ[i][1]:closed.uᵢ[i][2]]
    end
end

function getmu!(μ::AbstractVector{T}, closed::ClosedEmbeddedFunctions{T}; mu=:all) where T <: Number
    if mu == :all
        μ .= closed.μ
    elseif mu == :active
        i = 1
        for j = eachindex(closed.μ)
            if closed.𝕁[j]
                μ[i] = closed.μ[j]
                i += 1
            end
        end
    elseif mu == :inactive
        i = 1
        for j = eachindex(closed.μ)
            if !closed.𝕁[j]
                μ[i] = closed.μ[j]
                i += 1
            end
        end
    else
        throw(ArgumentError("Invalid option for mu; valid options are :all, :active, and :inactive"))
    end
    μ
end

function getmu(closed::ClosedEmbeddedFunctions{T}; mu=:all) where T <: Number
    if mu == :all
        return getmu!(zeros(T, length(closed.μ)), closed)
    elseif mu == :active
        return getmu!(zeros(T, sum(closed.𝕁)), closed)
    elseif mu == :inactive
        return getmu!(zeros(T, sum(.!closed.𝕁)), closed)
    else
        throw(ArgumentError("Invalid option for mu; valid options are :all, :active, and :inactive"))
    end
end

function setmu!(closed::ClosedEmbeddedFunctions{T}, μ::AbstractVector{T}; mu=:all) where T <: Number
    if mu == :all
        closed.μ .= μ
    elseif mu == :active
        i = 1
        for j = eachindex(closed.μ)
            if closed.𝕁[j]
                closed.μ[j] = μ[i]
                i += 1
            end
        end
        if i != length(μ) + 1
            throw(DimensionMismatch("μ is the wrong size"))
        end
    elseif mu == :inactive
        i = 1
        for j = eachindex(closed.μ)
            if !closed.𝕁[j]
                closed.μ[j] = μ[i]
                i += 1
            end
        end
        if i != length(μ) + 1
            throw(DimensionMismatch("μ is the wrong size"))
        end
    else
        throw(ArgumentError("Invalid option for mu; valid options are :all, :active, and :inactive"))
    end
    μ
end

function getvars!(v::AbstractVector{T}, closed::ClosedEmbeddedFunctions{T}) where T <: Number
    for i = eachindex(closed.u)
        v[closed.uᵢ[i][1]:closed.uᵢ[i][2]] .= closed.u[i].u
    end
    for i = eachindex(closed.μ)
        if closed.𝕁[i]
            v[closed.μᵢ[i]] = closed.μ[i]
        end
    end
    v
end

function getvars(closed::ClosedEmbeddedFunctions{T}) where T <: Number
    getvars!(zeros(T, closed.uᵢ[end][end] + sum(closed.𝕁)), closed)
end

function setvars!(closed::ClosedEmbeddedFunctions{T}, v::AbstractVector{T}) where T <: Number
    for i = eachindex(closed.u)
        closed.u[i].u .= v[closed.uᵢ[i][1]:closed.uᵢ[i][2]]
    end
    for i = eachindex(closed.μ)
        if closed.𝕁[i]
            closed.μ[i] = v[closed.μᵢ[i]]
        end
    end
    v
end

end  # module


text:
module NtToolBox

# package code goes here

using LinearAlgebra

export load_image, imageplot, rescale, clamP, snr, perform_wavelet_transf, plot_wavelet, compute_wavelet_filter, perform_wavortho_transf, plot_levelset, gaussian_blur, bilinear_interpolate, Grad, Div, perform_haar_transf, load_sound, perform_stft, plot_spectrogram, compute_max, perform_blurring, plot_vf, meshgrid, read_mesh, compute_boundary, plot_mesh, compute_normal, perform_linprog, plot_hufftree, perform_conjugate_gradient, perform_dijkstra_fm, cconvol, subsampling, upsampling, trim_dim, rand_discr, plot_multiclasses

## m must be of type Array{Float32, 3} or Array{Float64, 3}, it makes products between the vectors of the 3rd dimension and v

function Mdot(m, v::Array{Float64, 1})
  w = ones(size(m, 1), size(m, 2))
  for i in 1:size(m, 1)
    for j in 1:size(m, 2)
      w[i, j] = (m[i, j, :]'*v)[1]
    end
  end
  return w
end

##

## Rescale linearly the dynamic of a vector to fit within a range [a,b]

function rescale(f, a = 0, b = 1)
  v = maximum(f) - minimum(f)
  g = copy(f - minimum(f))
  if v > 0
    g = g / v
  end
  return a + g*(b - a)
end

##

## Function clamp

function clamP(x, a = [], b = [])

  """
   clamp - clamp a value

     y = clamp(x,a,b);

   Default is [a,b]=[0,1].

     Copyright (c) 2004 Gabriel Peyre
  """

  if a == []
    a = 0.0
  end
  if b == []
    b = 1.0
  end
  return min(max(x, a), b)
end

#function plot_levelset(Z, level = 0, f = [])
#    """
#        f is supposed to be of the same shape as Z
#    """
#     if length(f) == 0
#         f = copy(Z)
#     end
#
#     (n, p) = size(Z)
#     (X, Y) = meshgrid(collect(0 : n - 1), collect(0 : p - 1))
#     contour(X, Y, Z, [level], linewidths = 2, colors = "red")
#     imageplot(f)
# end

function gaussian_blur(f, sigma)

      """ gaussian_blur - gaussian blurs an image
      %
      %   M = perform_blurring(M, sigma, options);
      %
      %   M is the original data
      %   sigma is the std of the Gaussian blur (in pixels)
      %
      %   Copyright (c) 2007 Gabriel Peyre
      """
      if sigma <= 0
          return
      end
      n = maximum(size(f))
      t = [collect(0:n/2); collect(-n/2:-2)]
      (Y, X) = meshgrid(t, t)
      h = exp( -(X.^2 + Y.^2)./(2.0*float(sigma)^2) )
      h = h./sum(h)
      return real(plan_ifft((plan_fft(f)*f).*(plan_fft(h)*h))*((plan_fft(f)*f).*(plan_fft(h)*h)))
end

function plot_levelset(Z, level = 0, f = [])
    """
        f is supposed to be of the same shape as Z
    """
    if length(f) == 0
        f = copy(Z)
    end

    (n, p) = size(Z)
    (X, Y) = meshgrid(collect(0 : n - 1), collect(0 : p - 1))
    contour(X, Y, Z, [level], linewidths = 2, colors="red")
    imageplot(f)
end

function compute_max(X,d)
    """ compute_max - compute maximum along dimension d

       Y,I = compute_max(X,d);

       Copyright (c) 2008 Gabriel Peyre
    """
    if ndims(X)>=2
        Y = maximum(X,d)
        I = mapslices(argmax,X,d)[:]
    elseif ndims(X)==1
        Y = maximum(X)
        I = indmaw(X)
    end
    return Y,I
end

function compute_max(X)
    if ndims(X)>=2
        Y = maximum(X,1)
        I = mapslices(argmax,X,1)[:]
    elseif ndims(X)==1
        Y = maximum(X)
        I = indmaw(X)
    end
    return Y,I
end

# function circshift(x, p)
#     """
#         Circular shift of an array.
#     """
#     y = copy(x)
#     y = cat(1, y[p[1]+1:end, :], y[1:p[1]+1, :])
#     if (size(x)[2] > 0) & (length(p) > 1)
#         y = cat(2, y[:, p[1]+1:end], y[:, 1:p[1]+1])
#     end
#     return y
# end

include("ndgrid.jl")
include("signal.jl")
include("perform_wavelet_transf.jl")
include("compute_wavelet_filter.jl")
include("Grad.jl")
include("Div.jl")
include("perform_blurring.jl")
include("read_bin.jl")
#include("isosurface.jl")
include("perform_thresholding.jl")
include("load_sound.jl")
include("perform_stft.jl")
include("plot_spectrogram.jl")
include("plot_vf.jl")
include("read_mesh.jl")
include("compute_boundary.jl")
include("plot_mesh.jl")
include("compute_normal.jl")
include("perform_linprog.jl")
include("plot_hufftree.jl")
include("perform_conjugate_gradient.jl")
include("graph.jl")
include("perform_haar_transf.jl")
#include("perform_fast_marching.jl")
include("rand_discr.jl")
include("plot_multiclasses.jl")

end # module


text:
function PSquarePercentile(arg0::PSquarePercentile)
    return PSquarePercentile((PSquarePercentile,), arg0)
end

function PSquarePercentile(arg0::jdouble)
    return PSquarePercentile((jdouble,), arg0)
end

function clear(obj::PSquarePercentile)
    return jcall(obj, "clear", void, ())
end

function copy(obj::PSquarePercentile)
    return jcall(obj, "copy", PSquarePercentile, ())
end

function equals(obj::PSquarePercentile, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function evaluate(obj::StorelessUnivariateStatistic, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function get_n(obj::PSquarePercentile)
    return jcall(obj, "getN", jlong, ())
end

function get_quantile(obj::PSquarePercentile)
    return jcall(obj, "getQuantile", jdouble, ())
end

function get_result(obj::PSquarePercentile)
    return jcall(obj, "getResult", jdouble, ())
end

function hash_code(obj::PSquarePercentile)
    return jcall(obj, "hashCode", jint, ())
end

function increment(obj::PSquarePercentile, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end

function new_markers(::Type{PSquarePercentile}, arg0::List, arg1::jdouble)
    return jcall(PSquarePercentile, "newMarkers", PSquarePercentile_PSquareMarkers, (List, jdouble), arg0, arg1)
end

function quantile(obj::PSquarePercentile)
    return jcall(obj, "quantile", jdouble, ())
end

function to_string(obj::PSquarePercentile)
    return jcall(obj, "toString", JString, ())
end


text:
module MixedModeDebugger

using Compat
using Cassette
using JuliaInterpreter: JuliaInterpreter, BreakpointSignature, BreakpointRef, on_breakpoints_updated
using JuliaInterpreter: Frame, enter_call, get_return

using Debugger

#=============== Exports ==================================================#
export @run_mixedmode

# reexport some useful things from JuliaInterpreter
using JuliaInterpreter: @bp, @breakpoint, breakpoint
export @breakpoint, breakpoint
#==========================================================================#

Cassette.@context MixModeDebugCtx
include("reflection_utils.jl")
include("run.jl")
include("make_overdub.jl")
include("remove_overdub.jl")
include("hooking.jl")

#==========================================================================#
function __init__()
    # Get any breakpoints that were already defined.
    for bp in JuliaInterpreter._breakpoints
        make_interpretted_mode_overdub(bp)
    end

    # Add hook to get any that will be defined later
    on_breakpoints_updated(breakpoint_hook)
end

end # module


text:
using Test
using SmartShockFinder

# Various Tests for the NNPP. Tests are limited to only 1 mesh & polydegree

@testset "Test_1d_Networks" begin
    n_dims = 1

    # add mesh and polynomdegrees to then generate Traindata
    meshes = SmartShockFinder.get_mesh_1d(1)
    polydeg = collect(2:3)
    train_setup = Traindata_Settings(polydeg, meshes)

    # datasets path should be cleared before testing
    path_datsets = joinpath("datasets", "1d")
    rm(path_datsets, recursive=true, force=true)

    @testset "Test1dNNPP" begin

        @test_nowarn mkpath(joinpath("datasets", "1d"))
        @test_nowarn mkpath(joinpath("networks", "1d"))
        network_type = NNPP(n_dims)

        # run test that generates a Train dataset for the NNPP
        @testset "Test_generate_Traindata_1dNNPP" begin
            @test_nowarn generate_traindataset1d(network_type, train_setup)
            @test isfile(joinpath(path_datsets, "traindata1dNNPP.h5"))
        end

        # run test similar to traindata. Now for Validation dataset
        @testset "Test_generate_Validdata_1dNNPP" begin
            @test_nowarn generate_validdataset1d(network_type, train_setup)
            @test isfile(joinpath(path_datsets, "validdata1dNNPP.h5"))
        end

        # now tested the training process with only 2 epochs
        @testset "Test_TrainNetwork_1dNNPP" begin

            # Parameters to train the Net
            η = 0.001                # learning rate
            β = 0.0001               # regularization paramater
            number_epochs = 2
            Sb = 500                 # batch size
            L = 10                   # early stopping parameter
            hidden_units = [10, 10, 10, 10, 10]

            model_setup = Model_Settings(η, β, number_epochs, Sb, L, hidden_units)

            @test_nowarn train_network(network_type, model_setup)
            @test_nowarn length(readdir(joinpath("networks", "1d"))) > 0
            @test_nowarn rm(joinpath("networks", "1d"), recursive=true, force=true)
            @test_nowarn rm(joinpath("datasets", "1d"), recursive=true, force=true)
        end
    end

    @testset "Test1dNNPPh" begin

        @test_nowarn mkpath(joinpath("datasets", "1d"))
        @test_nowarn mkpath(joinpath("networks", "1d"))
        network_type = NNPPh(n_dims)

        # run test that generates a Train dataset for the NNPRPh
        @testset "Test_generate_Traindata_1dNNPPh" begin
            @test_nowarn generate_traindataset1d(network_type, train_setup)
            @test isfile(joinpath("datasets", "1d", "traindata1dNNPPh.h5"))
        end

        # run test similar to traindata. Now for Validation dataset
        @testset "Test_generate_Validdata_1dNNPPh" begin
            @test_nowarn generate_validdataset1d(network_type, train_setup)
            @test isfile(joinpath("datasets", "1d", "validdata1dNNPPh.h5"))
        end

        # now tested the training process with only 2 epochs
        @testset "Test_TrainNetwork_1dNNPPh" begin

            # Parameters to train the Net
            η = 0.001                # learning rate
            β = 0.0001               # regularization paramater
            number_epochs = 2
            Sb = 500                 # batch size
            L = 10                   # early stopping parameter
            hidden_units = [20, 20, 20, 20, 20]

            model_setup = Model_Settings(η, β, number_epochs, Sb, L, hidden_units)

            @test_nowarn train_network(network_type, model_setup)
            @test_nowarn length(readdir(joinpath("networks", "1d"))) > 0
            @test_nowarn rm(joinpath("networks", "1d"), recursive=true, force=true)
            @test_nowarn rm(joinpath("datasets", "1d"), recursive=true, force=true)
        end
    end

    @testset "Test1dNNRH" begin

        @test_nowarn mkpath(joinpath("datasets", "1d"))
        @test_nowarn mkpath(joinpath("networks", "1d"))
        network_type = NNRH(n_dims)

        # run test that generates a Train dataset for the NNRH
        @testset "Test_generate_Traindata_1dNNRH" begin
            @test_nowarn generate_traindataset1d(network_type, train_setup)
            @test isfile(joinpath("datasets", "1d", "traindata1dNNRH.h5"))
        end

        # run test similar to traindata. Now for Validation dataset
        @testset "Test_generate_Validdata_1dNNRH" begin
            @test_nowarn generate_validdataset1d(network_type, train_setup)
            @test isfile(joinpath("datasets", "1d", "validdata1dNNRH.h5"))
        end

        # now tested the training process with only 2 epochs
        @testset "Test_TrainNetwork_1dNNRH" begin

            # Parameters to train the Net
            η = 0.002             # learning rate
            β = 0.009             # regularization paramater
            number_epochs = 2
            Sb = 500                 # batch size
            L = 10                   # early stopping parameter
            hidden_units = [10, 10, 10, 10, 10]

            model_setup = Model_Settings(η, β, number_epochs, Sb, L, hidden_units)

            @test_nowarn train_network(network_type, model_setup)
            @test_nowarn length(readdir(joinpath("networks", "1d"))) > 0
            @test_nowarn rm(joinpath("networks", "1d"), recursive=true, force=true)
            @test_nowarn rm(joinpath("datasets", "1d"), recursive=true, force=true)
        end
    end
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DiscriminativeCircuits
using DataFrames: DataFrame

# This tests are supposed to test queries on the circuits
@testset "Logistic Circuit Class Conditional" begin
    # Uses a Logistic Circuit with 4 variables, and tests 3 of the configurations to
    # match with python version.

    EPS = 1e-5; # TODO; Cannot do more with Float32
    logistic_circuit = zoo_dc("little_4var.circuit");
    @test logistic_circuit isa LogisticCircuit;

    # Step 1. Check Probabilities for 3 samples
    data = DataFrame(Bool.([0 0 0 0; 0 1 1 0; 0 0 1 1]), :auto);

    # Testing values before taking sigmoid, so will compare with output of `class_weights_per_instance`
    true_prob = [3.43147972 4.66740416;
                4.27595352 2.83503504;
                3.67415087 4.93793472]
    true_prob

    CLASSES = 2
    calc_prob = class_weights_per_instance(logistic_circuit, data)

    for i = 1:3
        for j = 1:2
            @test true_prob[i,j] ≈ calc_prob[i,j] atol= EPS;
        end
    end
end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
export rotation2D, OEDTypeBParam, getOEDTypeBParam
export getTransmissionMatrix

"""
type OEDTypeBParam <: AbstractOEDParam

Fields:

	r          - number of rays
	l          - number of projections (design parameters)
	Ftrue      - true models
	Mimg       - image domain
	Mgrid      - rotating grid
	c          - center of rotation
	noiseLevel -
	Q          - temp storage for rotation2D
	xc         - cell-centers of Mgrid
	T          - transmission matrix
	M          - forward operator
	thetaOld   - previous theta (to check if M must be rebuilt)
"""
type OEDTypeBParam <: AbstractOEDParam
	r::Int
	l::Int
	Ftrue::Array{Float64,2}
	Mimg::RegularMesh
	Mgrid::RegularMesh
	c::Vector{Float64}
	noiseLevel::Float64
	xc::Array{Float64}
	T::AbstractArray{Float64}
	Q::Array{Float64}
	M::Array{Float64}
	thetaOld::Vector{Float64}
end
function Base.display(param::OEDTypeBParam)
	println("--- OED.OEDTypeBParam ---")
	println("number of rays   = $(param.r)")
	println("number of proj.  = $(param.l)")
	println("number of images = $(size(param.Ftrue,2))")
	println("image size       = $(param.Mimg.n)")
	println("image domain     = $(param.Mimg.domain)")
	println("grid size        = $(param.Mgrid.n)")
	println("gird domain      = $(param.Mgrid.domain)")
	println("noise level      = $(param.noiseLevel)")
	println("rotation center  = $(param.c)")

	print("temp variables   = ")
	(!isempty(param.xc)) && print("xc $(size(param.xc)), ")
	(!isempty(param.T)) && print("T $(size(param.T)), ")
	(!isempty(param.Q)) && print("Q $(size(param.Q)).")
	(!isempty(param.M)) && print("M $(size(param.M)).")
	print("\n")
end

export getOEDTypeBParam
"""
function getOEDTypeBParam

constructs OEDTypeBParam

Required Input:

	r          - number of rays
	l          - number of projections
	Ftrue      - true models
	Mimg       - image domain
	Mgrid      - rotating grid

Optional Inputs:

	c          - center of rotation
	noiseLevel -

"""
function getOEDTypeBParam(r::Int,l::Int,Ftrue,Mimg::RegularMesh,Mgrid::RegularMesh;
	c=(Mgrid.domain[2:2:end]+Mgrid.domain[1:2:end])/2, noiseLevel=1e-4)
	return OEDTypeBParam(r,l,Ftrue,Mimg,Mgrid,c,noiseLevel,[],[],[],[],[])
end

getProblemSize(param::OEDTypeBParam) = (param.r*param.l,param.Mimg.nc,size(param.Ftrue,2))
getNumberOfDesignParameters(pFor::OEDTypeBParam) = pFor.l

function getGrid(pFor::OEDTypeBParam)
	if isempty(pFor.xc)
		pFor.xc = getCellCenteredGrid(pFor.Mgrid)
	end
	return pFor.xc
end

function rotation2D(w::Float64,pFor::OEDTypeBParam;doDerivative::Bool=false)
	if isempty(pFor.Q) || size(pFor.Q,1) != length(pFor.xc)
		xc     = getGrid(pFor)
		pFor.Q = getQ(xc)
	end
	return rotation2D(w,pFor.xc,Q=pFor.Q,doDerivative=doDerivative,c=pFor.c)
end

function getTransmissionMatrix(pFor)
	if isempty(pFor.T)
		Mgrid = pFor.Mgrid
		r     = pFor.r;

		start = floor(Int,(Mgrid.n[1]-r)/2);
		v = ones(1,Mgrid.n[1]);
		T = zeros(r,prod(Mgrid.n));
		T[:, start*Mgrid.n[1]+1:start*Mgrid.n[1]+r*Mgrid.n[1]] = kron(eye(r),v);
		pFor.T = T;
	end
	return pFor.T
end

function getM(theta::Vector{Float64}, pFor::OEDTypeBParam,doClear::Bool=true)
	# This function gets the projection matrix by rotation and summation
	if doClear || isempty(pFor.thetaOld) || norm(theta-pFor.thetaOld,Inf)/norm(theta,Inf) > 1e-14
		Mgrid = pFor.Mgrid
		Mimg  = pFor.Mimg
		r     = pFor.r;
		l     = pFor.l

		# Get projection angles in radians
		psi = theta*(pi/180);
		if length(psi)!=l; error("getM theta must be length $l."); end;
		xc = getCellCenteredGrid(Mgrid);

		M     = zeros(r*l,Mimg.nc);

		T = getTransmissionMatrix(pFor)
		for i = 1:l
		  yc  = rotation2D(psi[i],pFor)
		  R   = getLinearInterMatrix(Mimg,yc)

		  M[(i-1)*r+1:i*r,:] = T*R
		end
		pFor.thetaOld = theta; pFor.M= M
	end
	return pFor.M
end
function getdpMx(theta::Vector{Float64},x::Vector{Float64},pFor::OEDTypeBParam)
#  This function gets the Jacobian matrix times a vector f

	Mgrid = pFor.Mgrid
	Mimg  = pFor.Mimg

	# Get projection angles in radians and create grid
	psi = theta*(pi/180);
	xc  = getCellCenteredGrid(Mgrid);

	# Build transmission matrix T
	T  = getTransmissionMatrix(pFor)

	# Build Jacobian matrix
	Jac = Array{Array}(pFor.l)
	for i = 1:pFor.l
	  yc,yp = rotation2D(psi[i],pFor,doDerivative=true)
	  Rc,dR = linearInter(reshape(x,tuple(Mimg.n...)),Mimg.domain,yc,doDerivative=true)
	  Jac[i] = (pi/180)*T*(dR*yp)
	end
	dpM = SparseMatrixCSC(sparse(Jac[1]))
	for i=2:pFor.l
		dpM = blkdiag(dpM,SparseMatrixCSC(sparse(Jac[i])))
	end
	return dpM
end

function getdpMTx(theta::Vector{Float64}, x::Vector{Float64}, pFor::OEDTypeBParam)
# This function gets the Jacobian matrix (with respect to the transpose)
#   times a vector x

	Mgrid = pFor.Mgrid
	Mimg  = pFor.Mimg

	# Get projection angles in radians
	psi = theta*(pi/180)
	xc  = getGrid(pFor)

	T   = getTransmissionMatrix(pFor)
	TTx = T'*reshape(x,pFor.r,pFor.l);
	Jac = zeros(Mimg.nc,0);

	for i = 1:pFor.l
	  yc,yp  = rotation2D(psi[i],pFor,doDerivative=true);
	  AT,dAT = getLinearInterMatrixTranspose(Mimg.domain,Mimg.n,yc,Img=TTx[:,i],doDerivative=true)
	  Jac = [Jac (pi/180)*dAT*yp];
	end
	return Jac
end

function getD(theta::Vector{Float64},param::OEDTypeBParam)
	D  = getM(theta,param) *param.Ftrue
	D += param.noiseLevel*randn(size(D))
	return D
end

function getdpD(theta::Vector{Float64},k::Int,param::OEDTypeBParam)
	return getdpMx(theta,param.Ftrue[:,k],param)
end