text:
@testset "native" begin

include("definitions/native.jl")

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

if VERSION >= v"1.8-"
    using Cthulhu
    include(joinpath(dirname(pathof(Cthulhu)), "..", "test", "FakeTerminals.jl"))
    using .FakeTerminals

    test_interactive = true
else
    test_interactive = false
end

cread1(io) = readuntil(io, '↩'; keep=true)
cread(io) = cread1(io) * cread1(io)

@testset "reflection" begin
    job, _ = native_job(identity, (Int,))

    @test only(GPUCompiler.code_lowered(job)) isa Core.CodeInfo

    CI, rt = only(GPUCompiler.code_typed(job))
    @test rt === Int

    IR = sprint(io->GPUCompiler.code_warntype(io, job))
    @test contains(IR, "MethodInstance for identity")

    IR = sprint(io->GPUCompiler.code_llvm(io, job))
    @test contains(IR, "julia_identity")

    ASM = sprint(io->GPUCompiler.code_native(io, job))
    @test contains(ASM, "julia_identity")

    if test_interactive
        fake_terminal() do term, in, out
            T = @async begin
                GPUCompiler.code_typed(job, interactive=true, interruptexc=false, terminal=term)
            end
            lines = cread(out)
            @test contains(lines, "identity(x)")
            write(in, 'q')
            wait(T)
        end
    end
end

@testset "Compilation" begin
    @testset "Callable structs" begin
        struct MyCallable end
        (::MyCallable)(a, b) = a+b

        (CI, rt) = native_code_typed(MyCallable(), (Int, Int), kernel=false)[1]
        @test CI.slottypes[1] == Core.Compiler.Const(MyCallable())

        (CI, rt) = native_code_typed(typeof(MyCallable()), (Int, Int), kernel=false)[1]
        @test CI.slottypes[1] == Core.Compiler.Const(MyCallable())
    end

    @testset "Compilation database" begin
        @noinline inner(x) = x+1
        function outer(x)
            return inner(x)
        end

        job, _ = native_job(outer, (Int,))
        JuliaContext() do ctx
            ir, meta = GPUCompiler.compile(:llvm, job; ctx)

            meth = only(methods(outer, (Int,)))

            mis = filter(mi->mi.def == meth, keys(meta.compiled))
            @test length(mis) == 1

            other_mis = filter(mi->mi.def != meth, keys(meta.compiled))
            @test length(other_mis) == 1
            @test only(other_mis).def in methods(inner)
        end
    end

    @testset "Advanced database" begin
        foo(x) = sum(exp(fill(x, 10, 10)))

        job, _ = native_job(foo, (Float64,))
        JuliaContext() do ctx
            # shouldn't segfault
            ir, meta = GPUCompiler.compile(:llvm, job; ctx)

            meth = only(methods(foo, (Float64,)))

            mis = filter(mi->mi.def == meth, keys(meta.compiled))
            @test length(mis) == 1

            expfloat = filter(keys(meta.compiled)) do mi
                mi.def in methods(Base.exp) && mi.specTypes == Tuple{typeof(Base.exp), Float64}
            end
            @test length(expfloat) == 1
        end
    end
end

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

@testset "IR" begin

@testset "basic reflection" begin
    valid_kernel() = return
    invalid_kernel() = 1

    ir = sprint(io->native_code_llvm(io, valid_kernel, Tuple{}; optimize=false, dump_module=true))

    # module should contain our function + a generic call wrapper
    @test occursin(r"define\ .* void\ @.*julia_valid_kernel.*\(\)"x, ir)
    @test !occursin("define %jl_value_t* @jlcall_", ir)

    # there should be no debug metadata
    @test !occursin("!dbg", ir)

    @test native_code_llvm(devnull, invalid_kernel, Tuple{}) == nothing
    @test_throws KernelError native_code_llvm(devnull, invalid_kernel, Tuple{}; kernel=true) == nothing
end

@testset "unbound typevars" begin
    invalid_kernel() where {unbound} = return
    @test_throws KernelError native_code_llvm(devnull, invalid_kernel, Tuple{})
end

@testset "child functions" begin
    # we often test using `@noinline sink` child functions, so test whether these survive
    @noinline child(i) = sink(i)
    parent(i) = child(i)

    ir = sprint(io->native_code_llvm(io, parent, Tuple{Int}))
    @test occursin(r"call .+ @julia_child_", ir)
end

@testset "sysimg" begin
    # bug: use a system image function

    function foobar(a,i)
        Base.pointerset(a, 0, mod1(i,10), 8)
    end

    ir = sprint(io->native_code_llvm(io, foobar, Tuple{Ptr{Int},Int}))
    @test !occursin("jlsys_", ir)
end

@testset "tracked pointers" begin
    function kernel(a)
        a[1] = 1
        return
    end

    # this used to throw an LLVM assertion (#223)
    native_code_llvm(devnull, kernel, Tuple{Vector{Int}}; kernel=true)
end

@testset "CUDAnative.jl#278" begin
    # codegen idempotency
    # NOTE: this isn't fixed, but surfaces here due to bad inference of checked_sub
    # NOTE: with the fix to print_to_string this doesn't error anymore,
    #       but still have a test to make sure it doesn't regress
    native_code_llvm(devnull, Base.checked_sub, Tuple{Int,Int}; optimize=false)
    native_code_llvm(devnull, Base.checked_sub, Tuple{Int,Int}; optimize=false)

    # breaking recursion in print_to_string makes it possible to compile
    # even in the presence of the above bug
    native_code_llvm(devnull, Base.print_to_string, Tuple{Int,Int}; optimize=false)
end

@testset "LLVM D32593" begin
    @eval struct D32593_struct
        foo::Float32
        bar::Float32
    end

    D32593(ptr) = unsafe_load(ptr).foo

    native_code_llvm(devnull, D32593, Tuple{Ptr{D32593_struct}})
end

@testset "slow abi" begin
    x = 2
    f = () -> x+1
    ir = sprint(io->native_code_llvm(io, f, Tuple{}, entry_abi=:func, dump_module=true))
    @test occursin(r"define nonnull {}\* @jfptr", ir)
    @test occursin(r"define internal fastcc .+ @julia", ir)
    @test occursin(r"call fastcc .+ @julia", ir)
end

end

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

@testset "assembly" begin

@testset "basic reflection" begin
    valid_kernel() = return
    invalid_kernel() = 1

    @test native_code_native(devnull, valid_kernel, Tuple{}) == nothing
    @test native_code_native(devnull, invalid_kernel, Tuple{}) == nothing
    @test_throws KernelError native_code_native(devnull, invalid_kernel, Tuple{}; kernel=true)
end

@testset "idempotency" begin
    # bug: generate code twice for the same kernel (jl_to_ptx wasn't idempotent)

    kernel() = return
    native_code_native(devnull, kernel, Tuple{})
    native_code_native(devnull, kernel, Tuple{})
end

@testset "compile for host after gpu" begin
    # issue #11: re-using host functions after GPU compilation
    @noinline child(i) = sink(i+1)

    function fromhost()
        child(10)
    end

    function fromptx()
        child(10)
        return
    end

    native_code_native(devnull, fromptx, Tuple{})
    @test fromhost() == 11
end

end

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

@testset "errors" begin

@eval Main begin
struct CleverType{T}
    x::T
end
Base.unsafe_trunc(::Type{Int}, x::CleverType) = unsafe_trunc(Int, x.x)
end

@testset "non-isbits arguments" begin
    foobar(i) = (sink(unsafe_trunc(Int,i)); return)

    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{BigInt})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("BigInt", msg)
    end

    # test that we can handle abstract types
    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{Any})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("Any", msg)
    end

    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{Union{Int32, Int64}})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("Union{Int32, Int64}", msg)
    end

    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{Union{Int32, Int64}})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("Union{Int32, Int64}", msg)
    end

    # test that we get information about fields and reason why something is not isbits
    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{CleverType{BigInt}})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("CleverType", msg) &&
        occursin("BigInt", msg)
    end
end

@testset "invalid LLVM IR" begin
    foobar(i) = println(i)

    @test_throws_message(InvalidIRError,
                         native_code_execution(foobar, Tuple{Int})) do msg
        occursin("invalid LLVM IR", msg) &&
        (occursin(GPUCompiler.RUNTIME_FUNCTION, msg) ||
         occursin(GPUCompiler.UNKNOWN_FUNCTION, msg) ||
         occursin(GPUCompiler.DYNAMIC_CALL, msg)) &&
        occursin("[1] println", msg) &&
        occursin(r"\[2\] .*foobar", msg)
    end
end

@testset "invalid LLVM IR (ccall)" begin
    foobar(p) = (unsafe_store!(p, ccall(:time, Cint, ())); nothing)

    @test_throws_message(InvalidIRError,
                         native_code_execution(foobar, Tuple{Ptr{Int}})) do msg
        occursin("invalid LLVM IR", msg) &&
        occursin(GPUCompiler.POINTER_FUNCTION, msg) &&
        occursin(r"\[1\] .*foobar", msg)
    end
end

@testset "delayed bindings" begin
    kernel() = (undefined; return)

    @test_throws_message(InvalidIRError,
                         native_code_execution(kernel, Tuple{})) do msg
        occursin("invalid LLVM IR", msg) &&
        occursin(GPUCompiler.DELAYED_BINDING, msg) &&
        occursin("use of 'undefined'", msg) &&
        occursin(r"\[1\] .*kernel", msg)
    end
end

@testset "dynamic call (invoke)" begin
    @eval @noinline nospecialize_child(@nospecialize(i)) = i
    kernel(a, b) = (unsafe_store!(b, nospecialize_child(a)); return)

    @test_throws_message(InvalidIRError,
                         native_code_execution(kernel, Tuple{Int,Ptr{Int}})) do msg
        occursin("invalid LLVM IR", msg) &&
        occursin(GPUCompiler.DYNAMIC_CALL, msg) &&
        occursin("call to nospecialize_child", msg) &&
        occursin(r"\[1\] .+kernel", msg)
    end
end

@testset "dynamic call (apply)" begin
    func() = println(1)

    @test_throws_message(InvalidIRError,
                         native_code_execution(func, Tuple{})) do msg
        occursin("invalid LLVM IR", msg) &&
        occursin(GPUCompiler.DYNAMIC_CALL, msg) &&
        occursin("call to println", msg) &&
        occursin("[2] func", msg)
    end
end

end

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

@testset "LazyCodegen" begin
    import .LazyCodegen: call_delayed

    f(A) = (A[] += 42; nothing)

    global flag = [0]
    function caller()
        call_delayed(f, flag::Vector{Int})
    end
    @test caller() === nothing
    @test flag[] == 42

    ir = sprint(io->native_code_llvm(io, caller, Tuple{}, dump_module=true))
    @test occursin(r"add i64 %\d+, 42", ir)
    # NOTE: can't just look for `jl_f` here, since it may be inlined and optimized away.

    add(x, y) = x+y
    function call_add(x, y)
        call_delayed(add, x, y)
    end

    @test call_add(1, 3) == 4

    incr(r) = r[] += 1
    function call_incr(r)
        call_delayed(incr, r)
    end
    r = Ref{Int}(0)
    @test call_incr(r) == 1
    @test r[] == 1

    function call_real(c)
        call_delayed(real, c)
    end

    @test call_real(1.0+im) == 1.0

    # Test ABI removal
    ir = sprint(io->native_code_llvm(io, call_real, Tuple{ComplexF64}))
    @test !occursin("alloca", ir)

    ghostly_identity(x, y) = y
    @test call_delayed(ghostly_identity, nothing, 1) == 1

    # tests struct return
    if Sys.ARCH != :aarch64
        @test call_delayed(complex, 1.0, 2.0) == 1.0+2.0im
    else
        @test_broken call_delayed(complex, 1.0, 2.0) == 1.0+2.0im
    end

    throws(arr, i) = arr[i]
    @test call_delayed(throws, [1], 1) == 1
    @test_throws BoundsError call_delayed(throws, [1], 0)

    struct Closure
        x::Int64
    end
    (c::Closure)(b) = c.x+b

    @test call_delayed(Closure(3), 5) == 8

    struct Closure2
        x::Integer
    end
    (c::Closure2)(b) = c.x+b

    @test call_delayed(Closure2(3), 5) == 8
end

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

@testset "overrides" begin
    # NOTE: method overrides do not support redefinitions, so we use different kernels

    mod = @eval module $(gensym())
        kernel() = child()
        child() = 0
    end

    ir = sprint(io->native_code_llvm(io, mod.kernel, Tuple{}))
    @test occursin("ret i64 0", ir)

    mod = @eval module $(gensym())
        using ..GPUCompiler
        import ..method_table

        kernel() = child()
        child() = 0

        GPUCompiler.@override method_table child() = 1
    end

    ir = sprint(io->native_code_llvm(io, mod.kernel, Tuple{}))
    @test occursin("ret i64 1", ir)
end

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

end


text:
"""Accepts rating as a Float64 and then creates a TwoPartCost."""
function TwoPartCost(variable_cost::T, args...) where {T <: VarCostArgs}
    return TwoPartCost(VariableCost(variable_cost), args...)
end

"""Accepts rating as a Float64 and then creates a ThreePartCost."""
function ThreePartCost(variable_cost::T, args...) where {T <: VarCostArgs}
    return ThreePartCost(VariableCost(variable_cost), args...)
end

"""Accepts rating as a Float64 and then creates a TechRenewable."""
function RenewableFix(name::String, available::Bool, bus::Bus,
                        activepower::Float64, reactivepower::Float64,
                        prime_mover::PrimeMovers, rating::Float64)
    tech = TechRenewable(rating, prime_mover, nothing, 1.0)
    RenewableFix(name, available, bus, activepower, reactivepower, tech)
end

"""Accepts rating as a Float64 and then creates a TechRenewable."""
function RenewableDispatch(name::String, available::Bool, bus::Bus,
                           activepower::Float64, reactivepower::Float64,
                           prime_mover::PrimeMovers, rating::Float64, op_cost::TwoPartCost)
    tech = TechRenewable(rating, prime_mover, nothing, 1.0)
    return RenewableDispatch(name, available, bus, activepower, reactivepower, tech, op_cost)
end

function PowerLoadPF(name::String, available::Bool, bus::Bus,
                     model::Union{Nothing, LoadModel}, activepower::Float64,
                     maxactivepower::Float64, power_factor::Float64)
    maxreactivepower = maxactivepower * sin(acos(power_factor))
    reactivepower = activepower * sin(acos(power_factor))
    return PowerLoad(name,
                     available,
                     bus,
                     model,
                     activepower,
                     reactivepower,
                     maxactivepower,
                     maxreactivepower)
end

function PowerLoadPF(::Nothing)
    return PowerLoadPF("init", true, Bus(nothing), nothing, 0.0, 0.0, 1.0)
end

"""Accepts anglelimits as a Float64."""
function Line(name, available::Bool, activepower_flow::Float64,
    reactivepower_flow::Float64, arc::Arc, r, x, b, rate, anglelimits::Float64)
    return Line(name, available, activepower_flow, reactivepower_flow, arc::Arc, r, x, b, rate,
                (min=-anglelimits, max=anglelimits))
end

"""Allows construction with bus type specified as a string for legacy code."""
function Bus(number, name, bustype::String, angle, voltage, voltagelimits, basevoltage)
    return Bus(number, name, get_enum_value(BusType, bustype), angle, voltage,
               voltagelimits, basevoltage, InfrastructureSystemsInternal())
end

ext: jl


text:
using BSON
using Dates
using Distributed
using Documenter
using FilePathsBase: SystemPath
using InteractiveUtils
using Memento
using Pkg
using Random
using Serialization
using Suppressor
using Test

using JLSO
using JLSO: JLSOFile, LOGGER, upgrade_jlso

# To test different types from common external packages
using DataFrames
using Distributions
using TimeZones

include("test_data_setup.jl")
@testset "JLSO" begin
    include("backwards_compat.jl")
    include("JLSOFile.jl")
    include("file_io.jl")

    # The doctests fail on x86, so only run them on 64-bit hardware & Julia 1.5
    Sys.WORD_SIZE == 64 && v"1.5" <= VERSION < v"1.6" && doctest(JLSO)
end

max_stars_repo_name: invenia/JLSO

max_issues_repo_name: invenia/JLSO

max_forks_repo_name: invenia/JLSO


text:
@doc raw"""
    sumBassUB(p,q,m,c0,alpha,F,T,Pi)

This function implements Udell an Boyd's algorithm to maximize a sum of Bass functions
and returns the maximizer and maximum.
"""
function sumBassUB(p,q,m,c0,alpha,F,T,Pi::Real)
    n = length(p)
    # Generate model for solution
    f = Vector{Function}(undef,n)
    df = Vector{Function}(undef,n)
    z = Vector{Float64}(undef,n)
    for i in 1:n
        f[i] = t -> alpha[i]*Bass(t,p[i],q[i],m[i],c0[i])
        df[i] = t -> alpha[i]*Bass_derivative(t,p[i],q[i],m[i],c0[i])
        z[i] = Bass_inflection(p[i],q[i],m[i],c0[i])
    end

    C = ones(1,n)
    d = [Pi]

    A = zeros(0,n)
    b = zeros(0)

    problem = LinearSP(f, df, z, A, b, C, d)

    pq, bestnodes, lbs, ubs = solve_sp(F, T, problem,TOL=1e-6)

    node = dequeue!(pq)

    maximizer = Vector{Float64}(undef,n)
    maximizer .= node.x
    maximum = sum(alpha[i]*Bass(maximizer[i],p[i],q[i],m[i],c0[i]) for i in 1:n)
    return maximizer, maximum
end

    # Maximizer
    # node.x
    # Lower bound of maximum
    # node.lb
    # Upper bound of maximum
    # node.ub

max_stars_repo_name: andres-alban/res-alloc-bass-demands

max_issues_repo_name: andres-alban/res-alloc-bass-demands

max_forks_repo_name: andres-alban/res-alloc-bass-demands


text:
# Paired end reads
# ================

abstract type Pairing end
struct Paired <: Pairing end
struct UnPaired <: Pairing end

struct Substitution
    pos::UInt64
    base::DNA
end

struct Reads{P<:Pairing,V<:AbstractSequencingView}
    genome::Vector{LongSequence{DNAAlphabet{2}}}
    views::Vector{V}
    substitutions::Dict{Int,Vector{Substitution}}
end

@inline views(x::Reads) = x.views
@inline genome(x::Reads) = x.genome
@inline substitutions(x::Reads) = x.substitutions
@inline nreads(x::Reads) = length(x.views)
@inline pick_read(reads::Reads) = rand(1:nreads(reads))
@inline function nsubstitutions(x::Reads)
    subs = substitutions(x)
    return isempty(subs) ? 0 : sum(length(x) for x in values(subs))
end

# I/O and printing functions
# --------------------------

summary(x::Reads{Paired,T}) where {T} = string(nreads(x), " paired reads:")
summary(x::Reads{UnPaired,T}) where {T} = string(nreads(x), " unpaired reads:")

function printlen(x::Reads)
    mx, av, mn = summarize_lengths(views(x))
    return string(" Maximum read size: $mx\n Average read size: $av\n Minimum read size: $mn")
end

function Base.show(io::IO, reads::Reads)
    println(summary(reads))
    println(printlen(reads))
    println(string(" Number of errors: ", nsubstitutions(reads)))
end

# Constructors
# ------------

function Reads{Paired,V}(p::Molecules{V}, flen::Int, rlen::Int = flen) where {V<:AbstractSequencingView}
    vs = views(p)
    vse = take_paired_ends(vs, flen, rlen)
    rds = Reads{Paired,V}(genome(p), vse, Dict{Int,Vector{Substitution}}())
    return rds
end

function Reads{UnPaired,V}(p::Molecules{V}, len::Union{Int,Nothing}) where {V<:AbstractSequencingView}
    vs = views(p)
    vse = take_single_ends(vs, len)
    rds = Reads{UnPaired,V}(genome(p), vse, Dict{Int, Vector{Substitution}}())
    return rds
end

"""
    paired_reads(p::MoleculePool, flen::Int, rlen::Int = flen)

Create a set of paired-end reads from a pool of DNA molecules `p`.

`flen` sets the length of forward read, and `rlen` sets the length of the
reverse read. If you only provide `flen`, then the function sets `rlen = flen`.

!!! note
    If a molecule in the pool is not long enough to create a forward and/or
    reverse read, then that molecule will simply be skipped.
"""
function paired_reads(p::Molecules{T}, flen::Int, rlen::Int = flen) where {T<:AbstractSequencingView}
    return Reads{Paired,T}(p, flen, rlen)
end

"""
    unpaired_reads(p::Molecules{T}, len::Int) where {T<:AbstractSequencingView}

Create a set of single-end reads from a pool of DNA molecules `p`.

`len` sets the length of the reads.

The end (strand) from which the reading begins for each DNA molecule in the
pool is determined at random for each molecule, with 50:50 probability.

If you don't provide a value for `len`, then the function will read each
DNA molecule in it's entirety.

!!! note
    If a molecule in the pool is not long enough to create a forward and/or
    reverse read, then that molecule will simply be skipped.
"""
function unpaired_reads(p::Molecules{T}, len::Union{Int,Nothing}) where {T<:AbstractSequencingView}
    return Reads{UnPaired,T}(p, len)
end

const posbases = (
    (DNA_C, DNA_G, DNA_T),
    (DNA_A, DNA_G, DNA_T),
    (DNA_N, DNA_N, DNA_N),
    (DNA_A, DNA_C, DNA_T),
    (DNA_N, DNA_N, DNA_N),
    (DNA_N, DNA_N, DNA_N),
    (DNA_N, DNA_N, DNA_N),
    (DNA_A, DNA_C, DNA_G)
)

"""
A simple function object that can be used with [`edit_substitutions`](@ref) or
[`edit_substitutions!`](@ref).

Simply clears all sequencing errors - substitutions - for all reads.
"""
struct ClearSubstitutions <: Function end

function (f::ClearSubstitutions)(output::Vector{Substitution}, readseq::LongSequence{DNAAlphabet{2}})
    return empty!(output)
end

"""
A simple function object that can be used with [`edit_substitutions`](@ref) or
[`edit_substitutions!`](@ref).

Randomly applies sequencing errors - substitutions - to all reads according to a
fixed, uniform, per-base probability.
"""
struct FixedProbSubstitutions <: Function
    prob::Float64
end

function (f::FixedProbSubstitutions)(output::Vector{Substitution}, readseq::LongSequence{DNAAlphabet{2}})
    p = f.prob
    dice = rand(length(readseq))
    dosub = dice .< p
    @inbounds for pos in 1:length(readseq)
        if dosub[pos]
            basenuc = reinterpret(UInt8, readseq[pos])
            push!(output, Substitution(pos, rand(posbases[basenuc])))
        end
    end
end

"""
    edit_substitutions!(f::Function, reads::Reads)

Add or remove sequencing errors to reads in the form of single base substitutions,
according to so some error generating or deleting function `f`.

!!! note
    This is the in-place, mutating version of [`edit_substitutions`](@ref), and
    since Pseudoseq's experted sequencing API is functional, this version is not
    exported.

The function provided must be a function that accepts two arguments:
    1. A vector of Substitution structs, which the function will mutate.
    2. The nucleotide sequence of the read.

The function needs not return anything, and anything it does return will not be
used.

Some acceptable functions already provided with Pseudoseq include [`FixedProbSubstitutions`](@ref)
and [`ClearSubstitutions`](@ref). If you want to develop your own function, their
implementations provide a simple example of how it can be done. You can of course
also just create an anonymous function.
"""
function edit_substitutions!(f::Function, reads::Reads)
    vs = views(reads)
    subs = substitutions(reads)
    subsbuf = Vector{Substitution}()
    for (i, v) in enumerate(vs)
        readseq = extract_sequence(genome(reads), v)
        readsubs = get(subs, i, subsbuf)
        f(readsubs, readseq)
        if isempty(readsubs)           # Read no longer has any subs.
            delete!(subs, i)
        elseif readsubs === subsbuf    # Read did not have subs but now does.
            subs[i] = copy(subsbuf)
            empty!(subsbuf)
        else                           # Read did have subs, and still does.
            subs[i] = readsubs
        end
    end
end

"""
    edit_substitutions(f::Function, reads::Reads)

Add or remove sequencing errors to reads in the form of single base substitutions,
according to so some error generating or deleting function `f`.

The function provided must be a function that accepts two arguments:
    1. A vector of Substitution structs, which the function will mutate.
    2. The nucleotide sequence of the read.

The function needs not return anything, and anything it does return will not be
used.

Some acceptable functions already provided with Pseudoseq include [`FixedProbSubstitutions`](@ref)
and [`ClearSubstitutions`](@ref). If you want to develop your own function, their
implementations provide a simple example of how it can be done. You can of course
also just create an anonymous function.
"""
function edit_substitutions(f::Function, reads::Reads)
    newreads = typeof(reads)(genome(reads), views(reads), copy(substitutions(reads)))
    edit_substitutions!(f, newreads)
    return newreads
end

# Read file generation
# --------------------

@inline function _add_subs!(seq::BioSequence, subs::Vector{Substitution})
    for sub in subs
        seq[sub.pos] = sub.base
    end
end

function generate(wtr::FASTQ.Writer, reads::Reads{UnPaired,BasicSequencingView})
    vs = views(reads)
    g = genome(reads)
    subs = substitutions(reads)
    n = 0
    for i in eachindex(vs)
        v = vs[i]
        # Extract subsequence from the reference and add errors.
        seq = extract_sequence(g, v)
        if haskey(subs, i)
            _add_subs!(seq, subs[i])
        end
        # Make the quality string for the read, (currently no real meaning!).
        qual = Vector{Int}(undef, length(seq))
        FASTQ.encode_quality_string!(FASTQ.SANGER_QUAL_ENCODING, fill(30, length(seq)), qual, 1, length(seq))
        # Create the name for the read...
        fragname = string("Refseq_", seqid(v))
        fqread = FASTQ.Record(fragname, seq, qual)
        write(wtr, fqread)
        n += 1
    end
    @info string("- ✔ Wrote ", n, " single end reads to FASTQ file")
    return reads
end

function generate(R1W::FASTQ.Writer, R2W::FASTQ.Writer, reads::Reads{Paired,BasicSequencingView})
    vs = views(reads)
    g = genome(reads)
    subs = substitutions(reads)
    n = 0
    for i in eachindex(vs)
        v = vs[i]
        # Extract subsequence from the reference and add errors.
        seq = extract_sequence(g, v)
        if haskey(subs, i)
            _add_subs!(seq, subs[i])
        end
        # Make the quality string for the read, (currently no real meaning!).
        qual = Vector{Int}(undef, length(seq))
        FASTQ.encode_quality_string!(FASTQ.SANGER_QUAL_ENCODING, fill(30, length(seq)), qual, 1, length(seq))
        # Create the name for the read...
        if isodd(i)
            fragname = string("readpair_", i)
            fqread = FASTQ.Record(fragname, seq, qual)
            write(R1W, fqread)
        else
            fragname = string("readpair_", i - 1)
            fqread = FASTQ.Record(fragname, seq, qual)
            write(R2W, fqread)
        end
        n += 1
    end
    @info string("- ✔ Wrote ", n, " paired end reads to FASTQ file")
    return reads
end

function prepare_tags(reads::Reads{Paired,TaggedSequencingView})
    tags = keys(summarize_tags(views(reads)))
    tagrange = 0x0000000000000000:0x00000000ffffffff
    tagseqs = DNAMer{16}.(sample_values(tagrange, length(tags)))
    tagdict = Dict(zip(tags, tagseqs))
    return tagdict
end

function generate(R1W::FASTQ.Writer, R2W::FASTQ.Writer, reads::Reads{Paired,TaggedSequencingView})
    vs = views(reads)
    g = genome(reads)
    subs = substitutions(reads)

    bufseq = LongSequence{DNAAlphabet{2}}(rand(ACGT, 7))
    tagdict = prepare_tags(reads)
    n = 0

    for i in eachindex(vs)
        v = vs[i]
        # Extract subsequence from the reference and add errors.
        seq = extract_sequence(g, v)
        if haskey(subs, i)
            _add_subs!(seq, subs[i])
        end

        if isodd(i) # We are dealing with an R1 read
            tagseq = LongSequence{DNAAlphabet{2}}(tagdict[tag(v)])
            seq = tagseq * bufseq * seq
        end

        # Make the quality string for the read, (currently no real meaning!).
        qual = Vector{Int}(undef, length(seq))
        FASTQ.encode_quality_string!(FASTQ.SANGER_QUAL_ENCODING, fill(30, length(seq)), qual, 1, length(seq))

        # Create the name for the read...
        pairnum = Int(ceil(i / 2))
        fragname = string("readpair_", pairnum)
        fqread = FASTQ.Record(fragname, seq, qual)
        current_writer = ifelse(isodd(i), R1W, R2W)
        write(current_writer, fqread)

        n += 1
    end
    @info string("- ✔ Wrote ", n, " tagged paired end reads to FASTQ file")
    return reads
end

"""
    generate(R1name::String, R2name::String, reads::Reads{Paired,<:AbstractSequencingView})

This method only works for paired reads. Instead of interleaving R1 and R2 reads
in a single FASTQ file, R1 and R2 reads are partitioned into two seperate FASTQ
files.
"""
function generate(R1name::String, R2name::String, reads::Reads{Paired,<:AbstractSequencingView})
    R1W = open(FASTQ.Writer, R1name)
    R2W = open(FASTQ.Writer, R2name)
    try
        generate(R1W, R2W, reads)
    finally
        close(R1W)
        close(R2W)
    end
    return reads
end

generate(wtr::FASTQ.Writer, reads::Reads{Paired,<:AbstractSequencingView}) = generate(wtr, wtr, reads)

"""
    generate(filename::String, reads::Reads)

Write the `reads` out to a FASTQ formatted file with the given `filename`.

If this method is used with a paired-end read type, then the FASTQ file will
be interleaved; all R1 reads will be odd records, and all R2 reads will be even
records in the file.

!!! note
    Reads are named according to the sequence in the input genome they came
    from. e.g. `@Reference_1_R1` means the first sequence in the genome, and
    `@Reference_2_R1` means the second sequence in the genome.
"""
function generate(filename::String, reads::Reads)
    open(FASTQ.Writer, filename) do wtr
        generate(wtr, reads)
    end
    return reads
end

function substitution_hist(reads::Reads{Paired,BasicSequencingView})
    hist = zeros(summarize_lengths(views(reads))[1])
    for i in Base.OneTo(nreads(reads))
        for s in get(substitutions(reads), i, Substitution[])
            hist[s.pos] += 1
        end
    end
    return hist
end


text:
#==============================================================================================

    AmpModelFit.jl

    structs containing parameters
    of fitted amplification models

    Author: Tom Price
    Date:   June 2019

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

import DataStructures.OrderedDict
using JuMP

#==============================================================================================
    structs and constructors
==============================================================================================#

abstract type AmpModelFit end

struct MAK2Fit <: AmpModelFit
    max_d_idx   ::Int_T
    coef_strs   ::Vector{String}
    coefs       ::Vector{Float_T}
    status      ::Symbol
    obj_val     ::Float_T
    jmp_model   ::JuMP.Model
end

## constructor
MAK2Fit(;
    max_d_idx   ::Int_T             = 0,
    coef_strs   ::Vector{String}    = Vector{String}(),
    coefs       ::Vector{Float_T}   = zeros(0),
    status      ::Symbol            = :not_fitted,
    obj_val     ::Float_T           = 0.0,
    jmp_model   ::JuMP.Model        = JuMP.Model(),
) =
    MAK2Fit(max_d_idx, coef_strs, coefs,
        status, obj_val, jmp_model)

struct MAK3Fit <: AmpModelFit
    max_d_idx   ::Int_T
    fb_start    ::Float_T
    bl_k_start  ::Float_T
    coef_strs   ::Vector{String}
    coefs       ::Vector{Float_T}
    status      ::Symbol
    obj_val     ::Float_T
    jmp_model   ::JuMP.Model
end

## constructor
MAK3Fit(;
    max_d_idx   ::Int_T             = 0,
    fb_start    ::Float_T           = 0.0,
    bl_k_start  ::Float_T           = 0.0,
    coef_strs   ::Vector{String}    = Vector{String}(),
    coefs       ::Vector{Float_T}   = zeros(0),
    status      ::Symbol            = :not_fitted,
    obj_val     ::Float_T           = 0.0,
    jmp_model   ::JuMP.Model        = JuMP.Model()
) =
    MAK3Fit(max_d_idx, fb_start, bl_k_start,
        coef_strs, coefs, status, obj_val, jmp_model)

struct MAKERGAUL3Fit <: AmpModelFit
    max_of_idx  ::Int_T
    coef_strs   ::Vector{String}
    coefs       ::Vector{Float_T}
    status      ::Symbol
    obj_val     ::Float_T
    jmp_model   ::JuMP.Model
end

## constructor
MAKERGAUL3Fit(;
    max_of_idx  ::Int_T             = 0,
    coef_strs   ::Vector{String}    = Vector{String}(),
    coefs       ::Vector{Float_T}   = zeros(0),
    status      ::Symbol            = :not_fitted,
    obj_val     ::Float_T           = 0.0,
    jmp_model   ::JuMP.Model        = JuMP.Model()
) =
    MAKERGAUL3Fit(max_of_idx, coef_strs, coefs,
        status, obj_val, jmp_model)

struct MAKERGAUL4Fit <: AmpModelFit
    max_of_idx  ::Int_T
    fb_start    ::Float_T
    bl_k_start  ::Float_T
    coef_strs   ::Vector{String}
    coefs       ::Vector{Float_T}
    status      ::Symbol
    obj_val     ::Float_T
    jmp_model   ::JuMP.Model
end

## constructor
MAKERGAUL4Fit(;
    max_of_idx  ::Int_T             = 0,
    fb_start    ::Float_T           = 0.0,
    bl_k_start  ::Float_T           = 0.0,
    coef_strs   ::Vector{String}    = Vector{String}(),
    coefs       ::Vector{Float_T}   = zeros(0),
    status      ::Symbol            = :not_fitted,
    obj_val     ::Float_T           = 0.0,
    jmp_model   ::JuMP.Model        = JuMP.Model()
) =
    MAKERGAUL4Fit(max_of_idx, fb_start, bl_k_start,
        coef_strs, coefs, status, obj_val, jmp_model)

struct SFCFit <: AmpModelFit
    coef_strs   ::Vector{String}
    coefs       ::Vector{Float_T}
    status      ::Symbol
    obj_val     ::Float_T
    jmp_model   ::JuMP.Model
    init_coefs  ::OrderedDict{String,Float_T}
end

## constructor
SFCFit(;
    coef_strs   ::Vector{String}                = Vector{String}(),
    coefs       ::Vector{Float_T}               = zeros(0),
    status      ::Symbol                        = :not_fitted,
    obj_val     ::Float_T                 = 0.0,
    jmp_model   ::JuMP.Model                    = JuMP.Model(),
    init_coefs  ::OrderedDict{String,Float_T}   = OrderedDict{String,Float_T}()
) =
    SFCFit(coef_strs, coefs, status, obj_val, jmp_model, init_coefs)

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
# Literal value.

mutable struct LiteralClause <: AbstractSQLClause
    val

    LiteralClause(; val) =
        new(val)
end

LiteralClause(val) =
    LiteralClause(val = val)

"""
    LIT(; val)
    LIT(val)

A SQL literal.

In a context of a SQL clause, `missing`, numbers, strings and datetime values
are automatically converted to SQL literals.

# Examples

```jldoctest
julia> c = LIT(missing);

julia> print(render(c))
NULL
```

```jldoctest
julia> c = LIT("SQL is fun!");

julia> print(render(c))
'SQL is fun!'
```
"""
LIT(args...; kws...) =
    LiteralClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(LIT), pats::Vector{Any}) =
    dissect(scr, LiteralClause, pats)

Base.convert(::Type{AbstractSQLClause}, val::SQLLiteralType) =
    LiteralClause(val)

PrettyPrinting.quoteof(c::LiteralClause, ::QuoteContext) =
    Expr(:call, nameof(LIT), c.val)


text:
function (bquad::AbstractBQ)(
    rng::Random.AbstractRNG,
    model::AbstractBQModel,
    sampler::AbstractMCMC.AbstractSampler;
    x_init=[],
    nsamples=200,
    callback=nothing,
)
    x = sample(rng, model, sampler, nsamples; callback=callback, bquad=bquad)
    return quadrature(bquad, model, x), x
end

function (bquad::AbstractBQ)(
    model::AbstractBQModel,
    sampler::AbstractMCMC.AbstractSampler;
    x_init=[],
    nsamples=200,
    callback=nothing,
)
    return bquad(
        GLOBAL_RNG, model, sampler; x_init=x_init, nsamples=nsamples, callback=callback
    )
end


text:
# In addition to illustrating transformations, we'll start loading some
# shaders from a file.  We'll also demo using Reactive for animations.
using ModernGL, GeometryTypes, GLAbstraction, GLWindow, Images, FileIO, Reactive

# Load our textures. See "downloads.jl" to get the images.
kitten = load(GLAbstraction.dir("tutorials", "images", "kitten.png"))
puppy  = load(GLAbstraction.dir("tutorials", "images", "puppy.png"))

# Create the window. This sets all the hints and makes the context current.
window = create_glcontext("Transformations 1", resolution=(600,600))

vao = glGenVertexArrays()
glBindVertexArray(vao)

# The positions of the vertices in our rectangle
vertex_positions = Point{2,Float32}[(-0.5,  0.5),     # top-left
                                    ( 0.5,  0.5),     # top-right
                                    ( 0.5, -0.5),     # bottom-right
                                    (-0.5, -0.5)]     # bottom-left

# The texture coordinates of each vertex
vertex_texcoords = Vec2f0[(0, 0),
                          (1, 0),
                          (1, 1),
                          (0, 1)]

# Specify how vertices are arranged into faces
elements = Face{3,UInt32,-1}[(0,1,2),          # the first triangle
                             (2,3,0)]          # the second triangle

vertex_shader = vert"""
#version 150

in vec2 position;
in vec3 color;
in vec2 texcoord;

out vec3 Color;
out vec2 Texcoord;

uniform mat4 trans;

void main()
{
    Color = color;
    Texcoord = texcoord;
    gl_Position = trans*vec4(position, 0.0, 1.0);
}
"""
fragment_shader = load(joinpath(dirname(@__FILE__), "shaders", "puppykitten.frag"))

# Define the rotation matrix (could also use rotationmatrix_z)
# By wrapping it in a Signal, we can easily update it.
trans = Signal(rotate(0f0, Vec((0,0,1f0))))

# Link everything together, using the corresponding shader variable as
# the Dict key
bufferdict = Dict(:position=>Buffer(vertex_positions),
                  :texcoord=>Buffer(vertex_texcoords),
                  :texKitten=>Texture(kitten'),
                  :texPuppy=>Texture(puppy'),
                  :trans=>trans,
                  :indexes=>indexbuffer(elements)) # special for element buffers

ro = std_renderobject(bufferdict,
                      LazyShader(vertex_shader, fragment_shader))
# Do the rendering
glClearColor(0,0,0,1)
while !GLFW.WindowShouldClose(window)
    push!(trans, rotationmatrix_z(time()*deg2rad(180)))
    Reactive.run_till_now()
    glClear(GL_COLOR_BUFFER_BIT)
    GLAbstraction.render(ro)
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end

GLFW.DestroyWindow(window)  # needed if you're running this from the REPL


text:
##########################################################################
#   	Basic functionality for doing MPS calculations with OBC     	#
##########################################################################
using LinearAlgebra
using Arpack

# Define sites as 3 legged tensors filled with elements of some type
Site{T} = Array{T,3}
# An MPS is an array of Sites of a certain type
MPS{T} = Vector{Site{T}}
# Define operators as 4 legged tensors filled with elements of some type
Operator{T} = Array{T,4}
# A MPO is an array of Operators
MPO{T} = Vector{Operator{T}}

"""
    random_mps_obc(N::Int, D::Int, d, tensortype::Type{T}=ComplexF64)::MPS{T} where T

Function to generate a random MPS with OBC
N: Number of sites
D: Bond dimension
d: Vector of physical dimensions, if a single integer dimension is given, it is assumed that all physical dimensions are the same, if varying dimensions are wanted a vector of integers of length N has to be supplied
The convention for the 3 index tensors is that the first two are the virtual ones, the last one is the physical one
   3
   |
1--A--2
"""
function random_mps_obc(N::Int, D::Int, d, tensortype::Type{T} = ComplexF64)::MPS{T} where {T}
    # Ensure the the system size is at least two
    @assert(N > 1)

    # Initialize
    mps = Array{Site{T}}(undef, N)

    if isa(d, Number)
        # Single number, I assume all dimensions are the same
        dim = d * ones(Int64, N)
    else
        # Make sure that the input is really a vector as julia distinguishes between nx1 matrices and column vectors
        dim = vec(d)
        # Make sure input is valid
        @assert(length(dim) == N)
    end

    # Left boundary tensor (row vector)
    mps[1] = rand(tensortype, 1, D, dim[1])

    # Right boundary tensor (column vector)
    mps[N] = rand(tensortype, D, 1, dim[N])

    # Tensors in between
    for i = 2:N-1
        mps[i] = rand(tensortype, D, D, dim[i])
    end

    return mps
end

"""
    basis_state_obc(configuration::Vector{<:Int}, d::Int=2)::MPS

Prepare the a product state corresponding |configuration> on N sites where configuration is an Array containing N elements from 1 to d. Each site is then initialized in the dth canonical basis state
"""
function basis_state_obc(configuration::Vector{<:Int}, d::Int = 2)::MPS{Float64}
    # Some error checking
    if any(x -> (x < 1 || x > d), configuration)
        throw(ArgumentError("configuration must contain integer elements in the range from 1 to d, got d=$(repr(d)), configuration=$(repr(configuration))"))
    end
    # Generate the MPS
    N = length(configuration)
    psi = MPS{Float64}(undef, N)
    tensors = Vector{Array{Float64,3}}(undef, d)
    for i = 1:d
        tmp = zeros(Float64, 1, 1, d)
        tmp[1, 1, i] = 1.0
        tensors[i] = tmp
    end
    for i = 1:N
        psi[i] = tensors[configuration[i]]
    end
    return psi
end

"""
    calculate_overlap(mps1::MPS,mps2::MPS)::Number

Given two mps, compute the overlap <mps1|mps2>.
"""
function calculate_overlap(mps1::MPS, mps2::MPS)::Number
    # Make sure, that the two MPS have the same length
    N1 = length(mps1)
    N2 = length(mps2)
    @assert(N1 == N2)

    # Now compute ther overlap
    overlap = ones(Float64, 1, 1)
    for i = 1:N1
        overlap = contract_tensors(overlap, [2], mps2[i], [1])
        overlap = contract_tensors(conj(mps1[i]), [1; 3], overlap, [1; 3])
    end

    return overlap[1]
end

"""
   expectation_value(mps::MPS, mpo::MPO)::Number

Given a MPS and a MPO compute the expectation value <mps|mpo|mps>.
"""
function expectation_value(mps::MPS, mpo::MPO)::Number
    # Make sure, that the two MPS have the same length
    N1 = length(mps)
    N2 = length(mpo)
    @assert(N1 == N2)

    # Contract zipper like
    val = ones(Float64, 1, 1, 1)
    for i = 1:N1
        val = contract_tensors(val, [1], conj(mps[i]), [1])
        val = contract_tensors(val, [1; 4], mpo[i], [1; 3])
        val = contract_tensors(val, [1; 4], mps[i], [1; 3])
    end

    return val[1]
end

"""
    gaugeMPS(mps::MPS{T}, direction::Symbol=:right, normalize::Bool=false)::MPS{T}

Function to bring an MPS with OBC in canonical form
direction: tells, whether it will be left or right canonical gauge
If normalize is set to true, the resulting state will be normalized.
"""
function gaugeMPS(mps::MPS{T}, direction::Symbol = :right, normalize::Bool = false)::MPS{T} where {T}
    mps_gauged = deepcopy(mps)
    gaugeMPS!(mps_gauged, direction, normalize)
    return mps_gauged
end

"""
    gaugeMPS!(mps::MPS,direction::Dir=right,)

Bring an MPS with OBC in canonical form. Direction tells, whether it will be left or right canonical gauge If normalize is set to true, the resulting state will be
normalized. This function overwrites the input MPS with its gauged version
"""
function gaugeMPS!(mps::MPS, direction::Symbol = :right, normalize::Bool = false)
    # Check that we got a meaningful direction
    if (direction != :left && direction != :right)
        throw(ArgumentError("direction must be :left or :right, got $(repr(direction))"))
    end
    # Start gauging
    N = length(mps)
    if (direction == :left)
        # Case that I want left canonical gauge
        M = mps[1]
        for i = 1:N-1
            mps[i], res = gauge_site(M, direction)
            M = contract_tensors(res, [2], mps[i+1], [1])
            if (i == (N - 1))
                if (normalize)
                    mps[N], _ = gauge_site(M, direction)
                else
                    mps[N] = M
                end
            end
        end
    else
        # Case that I want right canonical gauge
        M = mps[N]
        for i = N:-1:2
            mps[i], res = gauge_site(M, direction)
            M = contract_tensors(mps[i-1], [2], res, [1], [1; 3; 2])
            if (i == 2)
                if (normalize)
                    mps[1], _ = gauge_site(M, direction)
                else
                    mps[1] = M
                end
            end
        end
    end
end

"""
    gauge_site(A::Site{T}, direction::Dir)::Tuple{Site{T},Matrix{T}} where T

Bring a single tensor A of an MPS with OBC in canonical form
direction tells, whether it will be left or right normalized
The new tensor is M, the residual matrix which has to be multiplied
in the next tensor is stored in res
"""
function gauge_site(A::Site{T}, direction::Symbol)::Tuple{Site{T},Matrix{T}} where {T}
    if (direction == :left)
        Dl, Dr, d = size(A)
        M = permutedims(A, [3 1 2])
        M = reshape(M, (d * Dl, Dr))
        U, S, V = svd(M)
        dsv = length(S)
        U = reshape(U, (d, Dl, dsv))
        M = permutedims(U, [2 3 1])
        res = Matrix(Diagonal(S)) * V'
        return M, res
    else
        Dl, Dr, d = size(A)
        M = permutedims(A, [1 3 2])
        M = reshape(M, (Dl, Dr * d))
        U, S, V = svd(M)
        dsv = length(S)
        V = V'
        V = reshape(V, (dsv, d, Dr))
        M = permutedims(V, [1 3 2])
        res = U * Matrix(Diagonal(S))
        return M, res
    end
end

"""
    contract_virtual_indices(mps::MPS)::Vector{<:Number}

Given an MPS contract the virtual indices such that one obtains a dense vector. The indices are ordered such that they are compatible with the standard Julia kronecker product.

Warning: the object constructed will have exponential memory requirements in terms of the number of sites, use with care!
"""
function contract_virtual_indices(mps::MPS)::Vector{<:Number}
    N = length(mps)

    # Since we deal with open boundary conditions, we drop the dummy indices one on the left (right) boundary for the first (last) tensor manually. We start from the right to have the physical indices in the order compatible with Julia's kronecker product
    res = mps[N][:, 1, :]
    for i = N-1:-1:2
        res = contract_tensors(res, [ndims(res) - 1], mps[i], [2])
    end
    res = contract_tensors(res, [ndims(res) - 1], mps[1][1, :, :], [1])

    # Now reshape the result accordingly
    res = reshape(res, prod(size(res)))

    return res
end

"""
    contract_virtual_indices(mps::MPO)::Matrix{<:Number}

Given an MPO contract the virtual indices such that one obtains a dense matrix. The indices are ordered such that they are compatible with the standard Julia kronecker product.

Warning: the object constructed will have exponential memory requirements in terms of the number of sites, use with care!
"""
function contract_virtual_indices(mpo::MPO)::Matrix{<:Number}
    N = length(mpo)

    # Since we deal with open boundary conditions, we drop the dummy indices one on the left (right) boundary for the first (last) tensor manually. We start from the right to have the physical indices in the order compatible with Julia's kronecker product
    res = mpo[N][:, 1, :, :]
    for i = N-1:-1:2
        res = contract_tensors(res, [ndims(res) - 2], mpo[i], [2])
    end
    res = contract_tensors(res, [ndims(res) - 2], mpo[1][1, :, :, :], [1])

    # Reshuffle the indices to be compatible with Julia's standard kronecker product
    res = permutedims(res, [collect(1:2:ndims(res)); collect(2:2:ndims(res))])

    # Now reshape the result accordingly
    dims = size(res)
    dr = prod(dims[1:N])
    dc = prod(dims[N+1:end])
    res = reshape(res, (dr, dc))

    return res
end

"""
    find_groundstate(H::Array{Operator},D::Int64,acc::Float64,max_sweeps::Int64=-1)

Given a Hamiltonian MPO H, find an MPS approximation for its ground state with bond
dimension D converged to a relative accuracy acc. If a positive number for max_sweeps
is given, the maximum amount of iterations is limited to max_sweeps
"""
function find_groundstate(H::Vector{Operator{T}}, D::Int64, d::Int64, acc::Float64, max_sweeps::Int64) where {T}
    N = length(H)

    # Get the physical dimensions from the Hamiltonian
    d = Vector{Int64}(undef, length(H))
    for i = 1:length(H)
        d[i] = size(H[i], 3)
    end

    # Random MPS to start the calculation with
    mps = random_mps_obc(N, D, d)

    # Put it in right canonical gauge
    gaugeMPS!(mps, :right)
    # Precalculate the partial contractions
    LR = setup_R(H, mps)

    # Now start the sweeping
    num_of_sweeps = 0
    E0 = 0
    Eold = 1E5
    while (true)
        num_of_sweeps = num_of_sweeps + 1
        println("Sweep number ", num_of_sweeps)

        # From left to right starting by 1 up to N-1
        E_local = 0
        res = 0
        for i = 1:N-1
            Left = LR[i]
            Right = LR[i+1]
            E_local, M = solve_eigenvalue_problem(Left, Right, H[i])
            mps[i], _ = gauge_site(M, :left)
            LR[i+1] = update_left(Left, mps[i], mps[i], H[i])
        end

        # From left to right starting by 1 up to N-1
        for i = N:-1:2
            Left = LR[i]
            Right = LR[i+1]
            E_local, M = solve_eigenvalue_problem(Left, Right, H[i])
            mps[i], res = gauge_site(M, :right)
            LR[i] = update_right(Right, mps[i], mps[i], H[i])
        end

        # Check convergence criterion
        if (abs((Eold - E_local) / Eold) < acc)
            E0 = E_local
            # Restore normalisation (permutation is needed to bring the indices back into right order after contracting)
            mps[1] = contract_tensors(mps[1], [2], res, [1], [1; 3; 2])
            println("Convergence to desired accuracy achieved")
            break
        elseif (max_sweeps > 0 && num_of_sweeps >= max_sweeps)
            E0 = E_local
            # Restore normalisation (permutation is needed to bring the indices back into right order after contracting)
            mps[1] = contract_tensors(mps[1], [2], res, [1], [1; 3; 2])
            warn("Reached maximum number of iterations before convergence to desired accuracy")
            break
        end
        Eold = E_local
    end

    return E0, mps, num_of_sweeps
end

"""
    getHeff(Left,Right,W::Operator)

Construct the effective Hamiltonian
"""
function getHeff(Left, Right, W::Operator)
    Htemp = contract_tensors(W, [2], Right, [2])
    Htemp = contract_tensors(Left, [2], Htemp, [1])
    Htemp = permutedims(Htemp, [1; 5; 3; 2; 6; 4])
    dim1, dim2, dim3, dim4, dim5, dim6 = size(Htemp)
    Heff = reshape(Htemp, (dim1 * dim2 * dim3, dim4 * dim5 * dim6))

    return Heff
end

"""
    solve_eigenvalue_problem(Left,Right,W)

Function to construct and solve the eigenvalue problem which arises on each site for the effective Hamiltonian
"""
function solve_eigenvalue_problem(Left, Right, W)
    Htemp = contract_tensors(W, [2], Right, [2])
    Htemp = contract_tensors(Left, [2], Htemp, [1])
    Htemp = permutedims(Htemp, [1; 5; 3; 2; 6; 4])
    dim1, dim2, dim3, dim4, dim5, dim6 = size(Htemp)
    Heff = reshape(Htemp, (dim1 * dim2 * dim3, dim4 * dim5 * dim6))

    E_local, M = eigs(Heff, nev = 1, which = :SR)
    M = reshape(M, (dim1, dim5, dim3))

    return real(E_local[1]), M
end

"""
    setup_R(H::MPO{T1},mps::MPS{T2}) where {T1,T2}

Function to calculate the partial contractions needed to form the effective Hamiltonian
"""
function setup_R(H::MPO{T1}, mps::MPS{T2}) where {T1,T2}
    N = length(H)
    Tres = Base.return_types(*, (T1, T2))[1]
    LR = Vector{Array{Tres,3}}(undef, N + 1)

    # We need only N-1 partial contractions, however, we set the edges to dummy values 1 that we can recursively compute every contraction resuing the previous ones
    LR[1] = ones(Tres, 1, 1, 1)
    LR[N+1] = ones(Tres, 1, 1, 1)

    # Now compute the partial contractions starting from the right (as I start sweeping on the left in the ground state search)
    for i = N:-1:2
        tmp = contract_tensors(LR[i+1], [3], mps[i], [2])
        tmp = contract_tensors(H[i], [2; 4], tmp, [2; 4])
        LR[i] = contract_tensors(conj(mps[i]), [2; 3], tmp, [3; 2])
    end

    return LR
end

"""
    update_left(LR, M_left::Site, M_right::Site, W::Operator)

Function to perform an update the partial contractions required for iterative ground state search starting from the left stored in LR
"""
function update_left(LR, M_left::Site, M_right::Site, W::Operator)
    res = contract_tensors(LR, [3], M_right, [1])
    res = contract_tensors(W, [1; 4], res, [2; 4])
    res = contract_tensors(conj(M_left), [1; 3], res, [3; 2])

    return res
end

"""
    update_right(LR, M_left::Site, M_right::Site, W::Operator)

Function to perform an update the partial contractions required for iterative ground state search starting from the right stored in LR
"""
function update_right(LR, M_left::Site, M_right::Site, W::Operator)
    res = contract_tensors(LR, [3], M_right, [2])
    res = contract_tensors(W, [2; 4], res, [2; 4])
    res = contract_tensors(conj(M_left), [2; 3], res, [3; 2])

    return res
end

"""
    apply_operator(operator::MPO{T1}, mps::MPS{T2})::MPS where {T1,T2}

Apply an operator given as MPO to an MPS. The resulting MPS will have a bond dimension that is the product of the bond dimensions of the MPS and the MPO.
"""
function apply_operator(operator::MPO{T1}, mps::MPS{T2})::MPS where {T1,T2}
    N1 = length(mps)
    N2 = length(operator)
    @assert(N1 == N2)

    # Generate a new MPO of the correct type
    Tres = Base.return_types(*, (T1, T2))[1]
    res = MPS{Tres}(undef, N1)

    # Apply the MPO to the MPS and generate new MPS
    for i = 1:N1
        temp = contract_tensors(operator[i], [4], mps[i], [3])
        temp = permutedims(temp, (1, 4, 2, 5, 3))
        dim1, dim2, dim3, dim4, dim5 = size(temp)
        res[i] = reshape(temp, (dim1 * dim2, dim3 * dim4, dim5))
    end
    return res
end

"""
    apply_operator!(operator::MPO{T1}, mps::MPS{T2})::MPS where {T1,T2}

Apply an operator given as MPO to an MPS. The resulting MPS will have a bond dimension that is the product of the bond dimensions of the MPS and the MPO. The input will be overwritten by the result. This requires that the type of the elements of the input MPS is able to accomodate the result of multiplying the MPO tensors into the MPS tensors (e.g. applying a complex MPO to a real MPS cannot be done inpalce as the result will be complex)
"""
function apply_operator!(operator::MPO, mps::MPS)
    N1 = length(mps)
    N2 = length(operator)
    @assert(N1 == N2)

    # Contract the tensors for each site
    for i = 1:N1
        temp = contract_tensors(operator[i], [4], mps[i], [3])
        temp = permutedims(temp, (1, 4, 2, 5, 3))
        dim1, dim2, dim3, dim4, dim5 = size(temp)
        mps[i] = reshape(temp, (dim1 * dim2, dim3 * dim4, dim5))
    end
end

"""
    apply_operator(op1::MPO{T1}, op2::MPO{T2})::MPO where {T1,T2}

Multiply two MPOs together to get an expression for op2 * op1 in MPO form. The resulting MPS will have a bond dimension that is the product of the bond dimensions of both MPOs.
"""
function apply_operator(op1::MPO{T1}, op2::MPO{T2})::MPO where {T1,T2}
    N1 = length(op1)
    N2 = length(op2)
    @assert(N1 == N2)

    # Generate a new MPS of the correct type
    Tres = Base.return_types(*, (T1, T2))[1]
    res = MPO{Tres}(undef, N1)

    # Contract the tensors for each site
    for i = 1:N1
        temp = contract_tensors(op2[i], [4], op1[i], [3])
        temp = permutedims(temp, (1, 4, 2, 5, 3, 6))
        dim1, dim2, dim3, dim4, dim5, dim6 = size(temp)
        res[i] = reshape(temp, (dim1 * dim2, dim3 * dim4, dim5, dim6))
    end
    return res
end

"""
    sum_states(mps1::MPS{T1}, mps2::MPS{T2})::MPS where {T1,T2}

Add two MPSs together to get an expression for mps1 + mps2 in MPS form. The resulting MPS will have a bond dimension that is the sum of the bond dimensions of both MPOs.
"""
function sum_states(mps1::MPS{T1}, mps2::MPS{T2})::MPS where {T1,T2}
    N1 = length(mps1)
    N2 = length(mps2)
    @assert(N1 == N2)

    # Generate a new MPO of the correct type
    Tres = Base.return_types(+, (T1, T2))[1]
    res = MPS{Tres}(undef, N1)

    # The first site needs special treatment
    tensor1 = mps1[1]
    tensor2 = mps2[1]
    _, Dr1, d1 = size(tensor1)
    _, Dr2, d2 = size(tensor2)
    new_tensor = zeros(Tres, 1, Dr1 + Dr2, d1)
    for r = 1:d1
        new_tensor[1, :, r] = [tensor1[1:1, :, r] tensor2[1:1, :, r]]
    end
    res[1] = new_tensor
    # The tensors in between
    for i = 2:N1-1
        tensor1 = mps1[i]
        tensor2 = mps2[i]
        Dl1, Dr1, d1 = size(tensor1)
        Dl2, Dr2, d2 = size(tensor2)
        new_tensor = zeros(Tres, Dl1 + Dl2, Dr1 + Dr2, d1)
        for r = 1:d1
            new_tensor[1:Dl1, 1:Dr1, r] = tensor1[:, :, r]
            new_tensor[Dl1+1:end, Dr1+1:end, r] = tensor2[:, :, r]
        end
        res[i] = new_tensor
    end
    # The last site needs special treatment
    tensor1 = mps1[N1]
    tensor2 = mps2[N1]
    Dl1, _, d1 = size(tensor1)
    Dl2, _, d2 = size(tensor2)
    new_tensor = zeros(Tres, Dl1 + Dl2, 1, d1)
    for r = 1:d1
        new_tensor[:, 1, r] = [tensor1[:, 1, r]; tensor2[:, 1, r]]
    end
    res[N1] = new_tensor

    return res
end

"""
    sum_operators(op1::MPO{T1}, op2::MPO{T2})::MPO where {T1,T2}

Add two MPOs together to get an expression for op2 + op1 in MPO form. The resulting MPO will have a bond dimension that is the sum of the bond dimensions of both MPOs.
"""
function sum_operators(op1::MPO{T1}, op2::MPO{T2})::MPO where {T1,T2}
    N1 = length(op1)
    N2 = length(op2)
    @assert(N1 == N2)

    # Generate a new MPO of the correct type
    Tres = Base.return_types(+, (T1, T2))[1]
    res = MPO{Tres}(undef, N1)

    # The first site needs special treatment
    tensor1 = op1[1]
    tensor2 = op2[1]
    _, Dr1, dr1, dc1 = size(tensor1)
    _, Dr2, dr2, dc2 = size(tensor2)
    new_tensor = zeros(Tres, 1, Dr1 + Dr2, dr1, dc1)
    for r = 1:dr1
        for c = 1:dc1
            new_tensor[1, :, r, c] = [tensor1[1:1, :, r, c] tensor2[1:1, :, r, c]]
        end
    end
    res[1] = new_tensor
    # The tensors in between
    for i = 2:N1-1
        tensor1 = op1[i]
        tensor2 = op2[i]
        Dl1, Dr1, dr1, dc1 = size(tensor1)
        Dl2, Dr2, dr2, dc2 = size(tensor2)
        new_tensor = zeros(Tres, Dl1 + Dl2, Dr1 + Dr2, dr1, dc1)
        for r = 1:dr1
            for c = 1:dc1
                new_tensor[1:Dl1, 1:Dr1, r, c] = tensor1[:, :, r, c]
                new_tensor[Dl1+1:end, Dr1+1:end, r, c] = tensor2[:, :, r, c]
            end
        end
        res[i] = new_tensor
    end
    # The last site needs special treatment
    tensor1 = op1[N1]
    tensor2 = op2[N1]
    Dl1, _, dr1, dc1 = size(tensor1)
    Dl2, _, dr2, dc2 = size(tensor2)
    new_tensor = zeros(Tres, Dl1 + Dl2, 1, dr1, dc1)
    for r = 1:dr1
        for c = 1:dc1
            new_tensor[:, 1, r, c] = [tensor1[:, 1, r, c]; tensor2[:, 1, r, c]]
        end
    end
    res[N1] = new_tensor

    return res
end

"""
    compute_entropy(mps::MPS, n:Int)::Float64

Compute the von Neumann along for the bipartion of the sites into two subsets A={1,...,n} and B={n+1,...,N} where N is the length of the MPS. If n<=0 is supplied, a bipartion of into A. For the result to make sense, mps has to be a normalized quantum state.
"""
function compute_entropy(mps::MPS{T}, n::Int = 0)::Float64 where {T}
    # Extact the length and check if the given position is reasonable
    N = length(mps)
    @assert(n <= N)

    # The entropy of the entire state is simply zero since it is a pure state, so nothing to compute
    if n == N
        return 0.0
    end

    # In case a value n<=0 is given, we assume the bipartion is taken in the center
    if n <= 0
        n = Int(round(N / 2))
    end

    # Get a copy of the input
    mps_loc = deepcopy(mps)

    # Put the sites left to n into left canonical gauge
    M = mps_loc[1]
    for i = 1:n
        mps_loc[i], res = gauge_site(M, :left)
        M = contract_tensors(res, [2], mps_loc[i+1], [1])
    end
    mps_loc[n+1] = M

    # Now start contracting from the right boundary to obtain the reduced density operator
    rdm = ones(T, 1, 1)
    for i = N:-1:n+1
        rdm = contract_tensors(rdm, [2], mps_loc[i], [2])
        rdm = contract_tensors(conj(mps_loc[i]), [2; 3], rdm, [1; 3])
    end

    # Now diagonalize the reduced density matrix and compute the entropy
    ev = real(eigvals(rdm))
    # Eigenvalues that are numerically zero sometimes become -1E-16, to prevent problems with the logarithm we filter them
    ev = filter(x -> x > 0.0, ev)
    # The von Neumann entropy for the reduced density operator
    entropy = -sum(ev .* log2.(ev))

    return entropy
end

"""
    sample_from_mps!(mps::MPS, gauge_input::Bool=true)::Vector{Int64}

Generate a sample from the probability distribution of basis states defined by the MPS following A. Ferris, G. Vidal, PRB 85 165146 (2021). If the flag gague_input is set to true, the input MPS will be put in right canoncial gauge and normalized, which is a requirement for the algorithm to work. In case one is sure that the MPS is already properly gauged and normalized the gauging step can be spared by setting the flag to false. In this case the input will stay untouched.
"""
function sample_from_mps!(mps::MPS, gauge_input::Bool = true)::Vector{Int64}
    # Extact the length and check if the given position is reasonable
    N = length(mps)
    res = zeros(Int64, N)

    # Put the state into right canonical gauge and make sure it is normalized
    if gauge_input
        gaugeMPS!(mps, :right, true)
    end

    # Now sample from the MPS
    A = 0
    p = 0.0
    M = mps[1][1, :, :]
    for i = 1:N
        d = size(M, 2)
        pacc = 0
        r = rand()
        for l = 1:d
            # Prepare the basis state
            basis_state = zeros(d)
            basis_state[l] = 1.0
            # Contract it into the physical index of the tensor
            A = contract_tensors(M, [2], basis_state, [1])
            # Determine the probability for the basis vector e_l
            p = contract_tensors(conj(A), [1], A, [1])
            pacc += real(p[1])
            if r < pacc
                res[i] = l
                break
            end
        end
        if i < N
            # Contract the result into the next tensor
            M = contract_tensors(1 / sqrt(real(p[1])) * A, [1], mps[i+1], [1])
        end
    end
    return res
end

"""
    sample_from_mps(mps::MPS)::Vector{Int64}

Generate a sample from the probability distribution of basis states defined by the MPS following A. Ferris, G. Vidal, PRB 85 165146 (2021).
"""
function sample_from_mps(mps::MPS)::Vector{Int64}
    return sample_from_mps!(deepcopy(mps))
end

"""
    svd_compress_mps(mps::MPS, Dmax::Int,, tol::Real = 0.0)::MPS

Compress a given MPS applying an a singular value decomposition at each bond. If Dmax > 0 is supplied (and simultaneously tol = 0.0), a maximum of Dmax singular values is kept, thus truncating the MPS to one with maximum bond dimension Dmax. If tol > 0 is given (and simultaneously Dmax = 0) then all singular values > tol are kept. If both are specified then at most Dmax singular values > tol are kept.
"""
function svd_compress_mps(mps::MPS, Dmax::Int, tol::Real = 0.0)::MPS
    # One of the two parameters has to be larger than zero
    @assert((Dmax > 0) || (tol > 0))
    # Extract the length and prepare a result
    N = length(mps)
    res = deepcopy(mps)
    # Gauge in both directions that redundant dimensions at the boundaries are removed
    gaugeMPS!(res, :left)
    gaugeMPS!(res, :right)
    Dnew = 0
    for i = 1:N-1
        # Check if the given bond dimension is larger than Dmax
        if size(res[i], 2) > Dmax
            Dl1, _, d1 = size(res[i])
            _, Dr2, d2 = size(res[i+1])
            tmp = contract_tensors(res[i], [2], res[i+1], [1])
            tmp = reshape(tmp, (Dl1 * d1, Dr2 * d2))
            U, S, V = svd(tmp)
            M = diagm(S) * V'
            # Now truncate
            if Dmax > 0 && tol == 0.0
                U = U[:, 1:Dmax]
                M = M[1:Dmax, :]
                Dnew = Dmax
            else
                ind = findall(x -> x > tol, S)
                if Dmax > 0 && length(ind) > Dmax
                    ind = ind[1:Dmax]
                end
                U = U[:, ind]
                M = M[ind, :]
                Dnew = length(ind)
            end
            # Reshape and set new tensors
            res[i] = permutedims(reshape(U, (Dl1, d1, Dnew)), (1, 3, 2))
            res[i+1] = reshape(M, (Dnew, Dr2, d2))
        end
    end
    return res
end

"""
    decompose_into_mpo(M::Matrix{T}, d::Vector{Int}) where T <:Number

Given a many-body operator in form of a dense matrix, decompose it into an MPO. It is assumed that the many-body operator follows the index convention of Julia's built-in kronecker product. The vector d specifies the local dimensions of the Hilbert space.

# Examples
Decomposing a simple two-site operator made up from sum of Pauli terms into MPO form.
```julia-repl
julia>  Id = [1.0 0; 0.0 1.0]
julia>  X = [0.0 1.0; 1.0 0.0]
julia>  Z = [1.0 0.0; 0.0 -1.0]
julia>  H = kron(X,X) + kron(Id,Z) + kron(Z,Id)
julia>  mpo = decompose_into_mpo(H, 2)
```
"""
function decompose_into_mpo(M::Matrix{T}, d::Vector{Int})::MPO{T} where {T<:Number}
    N = length(d)
    dim = prod(d)
    if size(M) != (prod(d), prod(d))
        throw(ArgumentError("matrix not compatible with specified local dimensions, got size(M)=$(repr(size(M))), d=$(repr(d))"))
    end
    # The MPO holding the result
    res = MPO{T}(undef, N)
    # Reshape the matrix and rearrange the indices, such that the row and column index for each site are adjacent and add dummy indices 1 at the boundaries
    A = reshape(M, (1, reverse(d)..., reverse(d)..., 1))
    ind = collect(Iterators.flatten(zip(collect(N:-1:1), collect(2N:-1:N+1))))
    A = permutedims(A, (1, ind .+ 1..., 2 * length(d) + 2))
    # Now split it into tensors using an SVD
    Dl = 1
    for i = 1:N-1
        # Take the first three indices together
        dims = size(A)
        A = reshape(A, (prod(dims[1:3]), prod(dims[4:end])))
        # SVD the matrix representation
        U, S, V = svd(A)
        A = diagm(S) * V'
        Dr = size(U, 2)
        # Extract the tensor and the remaining part
        U = reshape(U, (Dl, d[i], d[i], Dr))
        res[i] = permutedims(U, (1, 4, 2, 3))
        A = reshape(A, (Dr, dims[4:end]...))
        Dl = Dr
    end
    res[N] = permutedims(A, (1, 4, 2, 3))

    return res
end

"""
    decompose_into_mpo(M::Matrix{T}, d::Int) where T <:Number

Simplified interface to the more general method assuming that all local dimensions are equal to d.
"""
function decompose_into_mpo(M::Matrix{T}, d::Int)::MPO{T} where {T<:Number}
    dl, dr = size(M)
    if dl != dr
        throw(ArgumentError("local dimensions must all be the same, obtained a matrix with dimensions $(repr((dl,dr)))"))
    end
    N = Int(round(log(d, dl)))
    return decompose_into_mpo(M, d * ones(Int64, N))
end


text:
using MiniFB
using Test

@testset "MiniFB.jl" begin
    # Can't really test this package automatically
    #    since the windows needs to be closed manually
    #    Best we do is to test the package loads.
    #    Look at examples directory for usage hints.
    @test 1+1 == 2
end


text:
"""
    get_value(s::Symbol)

Get access to the value that got saved in `s` by a previous action.
If you want to access a position or angle check out [`get_position`](@ref)
and [`get_angle`](@ref).

# Returns
- `Any`: the value stored by a previous action.
"""
function get_value(s::Symbol)
    is_internal = first(string(s)) == '_'
    if is_internal
        internal_sym = Symbol(string(s)[2:end])
        if hasfield(ActionSetting, internal_sym)
            return getfield(get_current_setting(), internal_sym)
        end
    end

    defs = CURRENT_VIDEO[1].defs
    if haskey(defs, s)
        return defs[s]
    else
        error("The symbol $s is not defined.")
    end
end

"""
    val(x)

`val` is just a short-hand for [`get_value`](@ref)
"""
val(x) = get_value(x)

get_position(p::Point) = p
get_position(t::Transformation) = t.p

"""
    get_position(s::Symbol)

Get access to the position that got saved in `s` by a previous action.

# Returns
- `Point`: the point stored by a previous action.
"""
get_position(s::Symbol) = get_position(val(s))

"""
    pos(x)

`pos` is just a short-hand for [`get_position`](@ref)
"""
pos(x) = get_position(x)

# As it is just the number tuple -> return it
get_scale(x::Tuple{<:Number,<:Number}) = x

# If just the number -> return it as a tuple
get_scale(x::Number) = (x, x)

"""
    get_scale(s::Symbol)

Get access to the scaling that got saved in `s` by a previous action.

# Returns
- `Scaling`: the scale stored by a previous action.
"""
get_scale(s::Symbol) = get_scale(val(s))

"""
    scl(x)

`scl` is just a short-hand for [`get_scale`](@ref)
"""
scl(x) = get_scale(x)

get_angle(t::Transformation) = t.angle

"""
    get_angle(s::Symbol)

Get access to the angle that got saved in `s` by a previous action.

# Returns
- `Float64`: the angle stored by a previous action i.e via `return Transformation(p, angle)`
"""
get_angle(s::Symbol) = get_angle(val(s))

"""
    ang(x)

`ang` is just a short-hand for [`get_angle`](@ref)
"""
ang(x) = get_angle(x)


text:
module Parsers
using Base.Meta
import Base.push!

export parsefile, defs, funcs, recurse,
    MetaCollector, FuncCollector, AbstractCollector,
    walk, inexpr, findfunc, findassign, replacevar,
    postwalk, prewalk, replace, inexpr

include("macrotools.jl")
include("findfunc.jl")

"""    parsefile(path)

read in a julia source file and parse it.

Note: If the top level is not a simple expression or module definition the file is wrapped in a Module named modprefix.
"""
function parsefile(path, modprefix="Modeling")
    s = read(path, String)
    # open(path) do fp
    #     s = read(String, fp)
    # end
    try
        expr = Base.Meta.parse(s)
        return expr
    catch
        s = "module $modprefix\n$s \nend #module $modprefix"
        expr = Base.Meta.parse(s)
        return expr
    end
end

"""    AbstractCollector

subtypes of AbstractCollector support extracting and collecting information
from input sources.
"""
abstract type AbstractCollector end

"""    FuncCollector{T} <: AbstractCollector

collects function definitions and names
"""
struct FuncCollector{T} <: AbstractCollector
    defs::T
end

function push!(fc::AbstractCollector, expr::LineNumberNode)
    return nothing
end

function push!(fc::FuncCollector, expr::Expr)
    if expr.head == :function
        push!(fc.defs, expr.args[1] => expr.args[2])
    end
end

"""   MetaCollector{T,U,V,W} <: AbstractCollector

collects multiple pieces of information such as

- exprs: expressions
- fc: functions
- vc: variable assignments
- modc: module imports
"""
struct MetaCollector{T,U,V,W} <: AbstractCollector
    exprs::V
    fc::T
    vc::U
    modc::W
end

function push!(mc::MetaCollector, expr::Expr)
    push!(mc.exprs, expr)
    push!(mc.fc, expr)
    if expr.head == :(=)
        @debug "pushing into vc" expr=expr
        push!(mc.vc, expr.args[1]=>expr.args[2])
    elseif expr.head == :using
        push!(mc.modc, expr.args[1].args)
    else
        @info "unknown expr type for metacollector"
        @show expr
    end
end

"""    funcs(body)

collect the function definitions from a module expression.
"""
function funcs(body)
    fs = FuncCollector([])
    for subexpr in body
        push!(fs, subexpr)
    end
    return fs
end

"""    defs(body)

collect the function definitions and variable assignments from a module expression.
"""
function defs(body)
    # fs = funcs(body)
    mc = MetaCollector(Any[], FuncCollector([]), Any[], Any[])
    for expr in body
        push!(mc, expr)
    end
    return mc
end

function recurse(mc::AbstractCollector)
    subdefs = Any[]
    funcdefs = mc.fc.defs
    @show funcdefs
    for def in funcdefs
        funcname = def[1]
        funcquote = def[2]
        push!(subdefs, funcname=>defs(funcquote.args))
    end
    return subdefs
end

end


text:
module Follow
using Plots
using BoundaryIntegralMethod
using BoundaryIntegralMethod.DirectSweep
using JLD

function followResonance(k0, b0, changeFn, steps; convergence = 1e-6, n=200, outputDir=nothing, fileBase="",
    maxDist=nothing, plotFn = nothing, startStep=1, namingFn = nothing, callbackFn = nothing)
  plots = []
  kfs = []
  bs = []
  b = b0
  k = k0
  if outputDir == nothing
    outputDir = "output/follow-$k0"
  end
  if fileBase != ""
    fileBase = fileBase+"-"
  end
  mkpath(outputDir)
  if plotFn == nothing
    plotFn = (b, xf, kf, baseName) -> exportTotalField(b, xf, kf; outFile="$baseName.h5", xRange=(-2,2), yRange=(-2,2), dx=.01, dy=.01);
  end

  for i in startStep:steps
    try
      (ks, xs) = directSweep(b, n, k)
      ks = ks[1:6]
      # Try closest resonance
      @debug ks
      p = sortperm(abs.(imag.(ks).-imag(k)))
      @info "Choosing from ks:"
      [@info "$i: $(ks[i])" for i in 1:length(ks)]
      # candidateK = ks[p[1]]
      # candidateX = xs[p[1],:]
      candidateK = ks[1]
      candidateX = xs[1,:]
      @info "Sortperm would be $(p[1]): $(ks[p[1]]), abs would be 1: $(ks[1])"
      @info "Choosing 1: $(ks[1])"

      (kf, xf) = residualInverseIteration(25, convergence, b, n, candidateK, candidateX)
      @info "At step $i, converged to $kf"
      paramName = if (namingFn != nothing)
        "-$(namingFn(b))"
      else
        ""
      end

      baseName= "$outputDir/$fileBase$i$paramName-$kf"
      field = plotFn(b, xf, kf, baseName)
      p = heatmap(clamp.(abs.(field), 0, .6), cbar=false, aspect_ratio=1, ticks=nothing, title="$i$paramName-k=$kf")
      savefig("$baseName.png")
      if maxDist != nothing && maxDist < abs(k-kf)
        @error "Stopping after $i steps, max distance between wavenumber exceeded in a step ($k -> $kf)"
        break
      end
      push!(kfs, kf)
      push!(plots, p)
      push!(bs, b)
      writeDatafile("$outputDir/$(fileBase)data.jld", kfs, bs)
      if callbackFn != nothing
        callbackFn(b, kf, xf, "$fileBase$i$paramName", baseName)
      end
      k = kf
      b = changeFn(b, i)
    catch e
      @error "Stopping after $i steps, caught error" exception=e
      break
    end
  end
  (kfs, plots)
end

function writeDatafile(name, kfs, bs)
  jldopen(name, "w") do file
    write(file, "ks/full", kfs)
    write(file, "ks/Re", real(kfs))
    write(file, "ks/Im", imag(kfs))
    write(file, "boundaries", bs)
  end
end

end


text:
#!/usr/bin/env julia

function assertmsg(assertion, msg)
  if !assertion
    println("assert-ERROR: $msg")
    exit(1)
  else
    println("assert-OK: $msg")
  end
end

## expects string, returns decoded string
function process(inp)
  rest = inp
  outp = ""
  loop = 0
  println("  input=$inp")
  m = match(r"^(.*?)\((\d+)x(\d+)\)(.*)$", rest)
  while m != nothing
    loop += 1
    println("regex matched for input $inp, loop=$loop")
    pre, len, rept, remain = m.captures
    leni = parse(Int64, len)
    repti = parse(Int64, rept)
    println("  found pre=$pre, cmd(len=$len, repeat=$rept), remain=$remain")
    outp = string(outp, pre)
    println("  outp=$outp; prepended")
    torepeat = remain[1:leni]
    insrt = repeat(torepeat, repti)
    println("  torepeat=$torepeat; insrt=$insrt")
    outp = string(outp, insrt)
    println("  outp=$outp; inserted")
    rest = remain[leni+1:end]
    println("  new rest=$rest")
    m = match(r"^(.*?)\((\d+)x(\d+)\)(.*)$", rest)
  end
  outp = string(outp, rest)
  outp
end

#** MAIN

assertmsg(true, "success message")
#assertmsg(false, "fail message")
#assertmsg(true, "success message")

inp = "ABC"
expct = "ABC"
res = process(inp)
assertmsg(res == expct, "result $res of input $inp should be $expct")

tests = [
  ["ADVENT" "ADVENT"];
  ["A(1x5)BC" "ABBBBBC"];
  ["(3x3)XYZ" "XYZXYZXYZ"];
  ["A(2x2)BCD(2x2)EFG" "ABCBCDEFEFG"];
  ["(6x1)(1x3)A" "(1x3)A"];
  ["X(8x2)(3x3)ABCY" "X(3x3)ABC(3x3)ABCY"]
]

println(typeof(tests))
println(string(tests[1,1], " => ", tests[1,2]))
for i = 1:size(tests,1)
  println(string(tests[i,1], " => ", tests[i,2]))

  inp = tests[i,1]
  expct = tests[i,2]
  res = process(inp)
  assertmsg(res == expct, "result $res of input $inp should be $expct")
end

open("day09_data.txt") do f
  global data = readstring(f)
  #println("data-file=>$data<")
end
data = strip(data)
println("INPUT data-str=>$data<")
res = process(data)
println("RESULT=>$res<")
println(string("RESult length = ", length(res)))


text:
module TestConstructors
using Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType

@testset "Type parameter constraints" begin
    # cannot use categorical value as level type
    @test_throws ArgumentError CategoricalPool{CategoricalValue{Int,UInt8}, UInt8, CategoricalValue{CategoricalValue{Int,UInt8},UInt8}}(
            CategoricalValue{Int,UInt8}[], Dict{CategoricalValue{Int,UInt8}, UInt8}(), UInt8[], false)
    # cannot use non-categorical value as categorical value type
    @test_throws ArgumentError CategoricalPool{Int, UInt8, Int}(Int[], Dict{Int, UInt8}(), UInt8[], false)
    # level type of the pool and categorical value should match
    @test_throws ArgumentError CategoricalPool{Int, UInt8, CategoricalValue{String, UInt8}}(Int[], Dict{Int, UInt8}(), UInt8[], false)
    # reference type of the pool and categorical value should match
    @test_throws ArgumentError CategoricalPool{Int, UInt8, CategoricalValue{Int, UInt16}}(Int[], Dict{Int, UInt8}(), UInt8[], false)
    # correct types combination
    @test CategoricalPool{Int, UInt8, CategoricalValue{Int, UInt8}}(Int[], Dict{Int, UInt8}(), UInt8[], false) isa CategoricalPool
end

@testset "empty CategoricalPool{String}" begin
    pool = CategoricalPool{String}()

    @test isa(pool, CategoricalPool{String})

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 0

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 0
end

@testset "empty CategoricalPool{Int}" begin
    pool = CategoricalPool{Int, UInt8}()

    @test isa(pool, CategoricalPool{Int, UInt8, CategoricalValue{Int, UInt8}})

    @test isa(pool.index, Vector{Int})
    @test length(pool.index) == 0

    @test isa(pool.invindex, Dict{Int, UInt8})
    @test length(pool.invindex) == 0
end

@testset "CategoricalPool{String, DefaultRefType}(a b c)" begin
    pool = CategoricalPool(["a", "b", "c"])

    @test isa(pool, CategoricalPool{String, UInt32, CategoricalValue{String, UInt32}})

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)
end

@testset "CategoricalPool{String, UInt8}(a b c)" begin
    pool = CategoricalPool{String, UInt8}(["a", "b", "c"])

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict{String, UInt8})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === UInt8(1)
    @test pool.invindex["b"] === UInt8(2)
    @test pool.invindex["c"] === UInt8(3)
end

@testset "CategoricalPool(a b c) with specified reference codes" begin
    pool = CategoricalPool(
        Dict(
            "a" => DefaultRefType(1),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(3),
        )
    )

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)
end

@testset "CategoricalPool(a b c) with specified Int ref codes" begin
    # TODO: Make sure that invindex input is exhaustive
    # Raise an error if map misses any entries
    pool = CategoricalPool(
        Dict(
            "a" => 1,
            "b" => 2,
            "c" => 3,
        )
    )

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)
end

@testset "CategoricalPool(c b a)" begin
    pool = CategoricalPool(["c", "b", "a"])

    @test isa(pool, CategoricalPool)

    @test length(pool.index) == 3
    @test pool.index[1] == "c"
    @test pool.index[2] == "b"
    @test pool.index[3] == "a"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)
end

@testset "CategoricalPool(a b c) with ref codes not matching the natural order" begin
    pool = CategoricalPool(
        Dict(
            "a" => DefaultRefType(3),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(1),
        )
    )

    @test isa(pool, CategoricalPool)

    @test length(pool.index) == 3
    @test pool.index[1] == "c"
    @test pool.index[2] == "b"
    @test pool.index[3] == "a"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)
end

@testset "CategoricalPool(a b c) with specified levels order" begin
    pool = CategoricalPool(["c", "b", "a"], ["c", "b", "a"])

    @test isa(pool, CategoricalPool)

    @test length(pool.index) == 3
    @test pool.index[1] == "c"
    @test pool.index[2] == "b"
    @test pool.index[3] == "a"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)
end

@testset "CategoricalPool(a b c) with specified index and levels order" begin
    pool = CategoricalPool(
        Dict(
            "a" => DefaultRefType(3),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(1),
        ),
        ["c", "b", "a"]
    )

    @test isa(pool, CategoricalPool)

    @test length(pool.index) == 3
    @test pool.index[1] == "c"
    @test pool.index[2] == "b"
    @test pool.index[3] == "a"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)
end

@testset "CategoricalPool{Float64, UInt8}()" begin
    pool = CategoricalPool{Float64, UInt8}([1.0, 2.0, 3.0])

    @test isa(pool, CategoricalPool{Float64, UInt8, CategoricalValue{Float64, UInt8}})
    @test CategoricalValue(1, pool) isa CategoricalValue{Float64, UInt8}
end

end