text:
mutable struct Hole
    _begin::vk.VkDeviceSize
    _end::vk.VkDeviceSize
end

struct MemoryBlock
    memory::MemoryChunk
    size::vk.VkDeviceSize
    holes::Vector{Hole}
end

mutable struct Suballocator
    mDevice
    mBufferImageGranularity::vk.VkDeviceSize
    mMemoryProperties::vk.VkPhysicalDeviceMemoryProperties
    mLargeHeapBlockSize::vk.VkDeviceSize
    mSmallHeapBlockSize::vk.VkDeviceSize
    mSmallHeapThreshold::vk.VkDeviceSize
    mTypeBlocks::Vector{Vector{MemoryBlock}}

    function Suballocator(device, biGranularity::vk.VkDeviceSize)
        this = new()
        this.mDevice = device
        this.mBufferImageGranularity = biGranularity
        this.mMemoryProperties = VkExt.getMemoryProperties(getPhysicalDevice(this.mDevice))

        this.mLargeHeapBlockSize = vk.VkDeviceSize(256 * 1024 * 1024)
        this.mSmallHeapBlockSize = vk.VkDeviceSize(64 * 1024 * 1024)
        this.mSmallHeapThreshold = vk.VkDeviceSize(512 * 1024 * 1024)

        this.mTypeBlocks = Vector{Vector{MemoryBlock}}(undef, vk.VK_MAX_MEMORY_TYPES)
        for i = 1 : length(this.mTypeBlocks)
            this.mTypeBlocks[i] = Vector{MemoryBlock}()
        end
        return this
    end
end

@class Suballocator

@enum MemoryType begin
    VRAM = 1
    RAM = 2
end

function allocate(this::Suballocator, req::vk.VkMemoryRequirements, type::MemoryType)::MemoryChunk
    return allocate(this, req, flagsForType(type))
end

function allocate(this::Suballocator, req::vk.VkMemoryRequirements, flags::vk.VkMemoryPropertyFlags)::MemoryChunk
    # the bufferImageGranularity is the size of memory pages that must not be
    # shared between linear resources (buffers) and non-linear resources
    # (images). In order to limit complexity, we (for now) assume that every
    # currently bound resource is of the "wrong" type
    # => increase the alignment to the page size
    newAlignment = max(req.alignment, this.mBufferImageGranularity)

    typeIdx = typeIndexFor(this, req, flags)
    type = this.mMemoryProperties.memoryTypes[typeIdx]
    heapIndex = type.heapIndex # start from 1
    heap = this.mMemoryProperties.memoryHeaps[heapIndex + 1]
    # allocate a large block if heap.size exceed the threshold
    blocksize = heap.size > this.mSmallHeapThreshold ? this.mLargeHeapBlockSize : this.mSmallHeapBlockSize

    if (req.size > blocksize)
        return allocateDedicated(this, req, flags)
    end

    toAlign = offset::vk.VkDeviceSize->begin
        remainder = offset % newAlignment
        if remainder == 0
            return offset
        end
        return offset + newAlignment - remainder;
    end
    deallocator = chunk::MemoryChunk->begin
        deallocate(this, chunk)
    end

    mappable = (type.propertyFlags & vk.VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT) != 0
    blocks = this.mTypeBlocks[typeIdx]
    logicalDevice = getLogicalDevice(this.mDevice)
    for block in blocks
        count = length(block.holes)
        for i = 1 : count
            hole = block.holes[i]
            diff = Int64(hole._end) - Int64(toAlign(hole._begin)) - Int64(req.size)
            if (diff == 0)
                # Allocation fits hole exactly, throw out the hole
                deleteat!(block.holes, i)
                return MemoryChunk(handle(block.memory), logicalDevice,
                                    hole._begin, toAlign(hole._begin), req.size,
                                    typeIdx, mappable, deallocator)
            elseif (diff > 0)
                # Put the allocation at the beginning of the hole
                newBegin = hole._begin
                # move the begin of the hole backward because the part is used for the new memory chunk
                hole._begin = toAlign(newBegin) + req.size;
                return MemoryChunk(handle(block.memory), logicalDevice,
                                    newBegin, toAlign(newBegin), req.size,
                                    typeIdx, mappable, deallocator)
            end
        end
    end

    # Still no hole found, add a new block
    createInfo = Ref(vk.VkMemoryAllocateInfo(
        vk.VK_STRUCTURE_TYPE_MEMORY_ALLOCATE_INFO, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        blocksize, #allocationSize::VkDeviceSize
        typeIdx - 1 #memoryTypeIndex::UInt32
    ))

    memory = VkExt.allocateMemory(logicalDevice, createInfo)
    blockchunk = MemoryChunk(memory, logicalDevice,
                            UInt64(0), UInt64(0), blocksize, typeIdx, mappable)
    push!(blocks, MemoryBlock(blockchunk, blocksize, [Hole(req.size, blocksize)]))

    return MemoryChunk(handle(blockchunk), logicalDevice,
                            UInt64(0), UInt64(0), req.size, typeIdx, mappable, deallocator)
end

function allocateDedicated(this::Suballocator, req::vk.VkMemoryRequirements, type::MemoryType)::MemoryChunk
    return allocateDedicated(this, req, flagsForType(type))
end

function allocateDedicated(this::Suballocator, req::vk.VkMemoryRequirements, flags::vk.VkMemoryPropertyFlags)::MemoryChunk
    typeIdx = typeIndexFor(this, req, flags)
    type = this.mMemoryProperties.memoryTypes[typeIdx]
    createInfo = Ref(vk.VkMemoryAllocateInfo(
        vk.VK_STRUCTURE_TYPE_MEMORY_ALLOCATE_INFO, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        req.size, #allocationSize::VkDeviceSize
        typeIdx - 1 #memoryTypeIndex::UInt32
    ))
    logicalDevice = getLogicalDevice(this.mDevice)
    memory = VkExt.allocateMemory(logicalDevice, createInfo)
    mappable = (type.propertyFlags & vk.VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT) != 0
    return MemoryChunk(memory, logicalDevice, UInt64(0), UInt64(0), req.size, typeIdx, mappable, nothing)
end

function flagsForType(type::MemoryType)::vk.VkMemoryPropertyFlags
    flags::vk.VkMemoryPropertyFlags = 0
    if type == VRAM
        flags = vk.VK_MEMORY_PROPERTY_DEVICE_LOCAL_BIT |
                vk.VK_MEMORY_PROPERTY_LAZILY_ALLOCATED_BIT
    elseif type == RAM
        flags = vk.VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT |
                vk.VK_MEMORY_PROPERTY_HOST_COHERENT_BIT |
                vk.VK_MEMORY_PROPERTY_HOST_CACHED_BIT
    end
    return flags
end

function typeIndexFor(this::Suballocator, req::vk.VkMemoryRequirements, flags::vk.VkMemoryPropertyFlags)::UInt32
    memtype = typemax(UInt32)
    best_popcount::Int = 0
    for i::UInt32 = 0 : this.mMemoryProperties.memoryTypeCount - 1
        if (req.memoryTypeBits & (1 << i)) != 0
            current_flags = this.mMemoryProperties.memoryTypes[i + 1].propertyFlags
            ones = count_ones(flags & current_flags)

            if (ones > best_popcount)
                best_popcount = ones
                memtype = i + 1
            end
        end
    end

    if (memtype == typemax(UInt32))
        error("Could not find appropriate memory class.")
    end
    return memtype
end

function deallocate(this::Suballocator, chunk::MemoryChunk)
    blocks = this.mTypeBlocks[chunk.mType]
    it = findfirst(block->begin return handle(block.memory) == handle(chunk) end, blocks)
    @assert (it != nothing) "Couldn't find the block this chunk belongs to."

    block = blocks[it]
    holes = block.holes

    directlyBefore = findfirst(hole->begin hole._end == getAllocationOffset(chunk) end, holes)
    directlyAfter = findfirst(hole->begin hole._begin == getOffset(chunk) + getSize(chunk) end, holes)

    realsize = getSize(chunk) + (getOffset(chunk) - getAllocationOffset(chunk))
    if (directlyBefore != nothing && directlyAfter != nothing)
        # the chunk connects two holes to one
        # => resize the first, delete the second
        holes[directlyBefore]._end = holes[directlyAfter]._end
        deleteat!(holes, directlyAfter)
    elseif (directlyBefore != nothing)
        holes[directlyBefore]._end += realsize
    elseif (directlyAfter != nothing)
        holes[directlyAfter]._begin = getAllocationOffset(chunk)
    else
        push!(holes, Hole(getAllocationOffset(chunk), getOffset(chunk) + getSize(chunk)))
    end

    if (length(holes) == 1 && holes[1]._begin == 0 && holes[1]._end == block.size)
        # give completely empty blocks back to the driver
        deleteat!(blocks, it)
    end
end

# TODO
# heldBlockCount


text:
module TestFloatingTableView

using FloatingTableView
using Test

const W = FloatingTableView.TableView.WebIO

@testset "FloatingTableView.jl" begin
	t = (a = [1], b = [2])

	@test showtable(t) isa W.Scope
	@test browse(t) === nothing

	t = [(a = 1, b = 2), (a = 3, b =4)]
	@test showtable(t) isa W.Scope
	@test browse(t) === nothing

	t = (a = ["a string"], b = ["another string"])
	@test showtable(t) isa W.Scope
	@test browse(t) === nothing

	t = [1 3; 3 4]
	@test showtable(t) isa W.Scope
    @test browse(t) === nothing

    t = [(a = 1, b = 2), (a = 3, b = 4)]
    @test browse(t; newwindow=true) === nothing
end

end

ext: jl


text:
using DataFrames, HTTP, CSV, Plots, StatsPlots, StatsBase, SplitApplyCombine, Missings

data_from_web = HTTP.request("GET", "https://raw.githubusercontent.com/saw191/"*"Stephanie_Walker_2504_2021_PROJECT3/main/data/Melbourne_housing_FULL.csv")
string_from_web = String(data_from_web.body);

df = CSV.read(IOBuffer(string_from_web), DataFrame)


text:
module CURANDkernel
using ..Codegen

export curandState_t, curandStateSobol32_t
export curand_init, curand, curand_uniform, curand_normal, curand_log_normal, curand_poisson

const __current_compiler__ = Codegen.__current_compiler__

cxx"""
#include <curand_kernel.h>
"""

"""
    curandState_t

"""
typealias curandState_t cxxt"curandState_t"{48}

export new, delete
new(::Type{curandState_t}) = icxx"return (curandState_t*) malloc(sizeof(curandState_t));"
delete(x) = icxx"free($x);"

"""
   curand_init(see, sequence, offset, state)

Initialisation for `curandState_t`

# Arguments
- `seed`
- `sequence`
- `offset`
- `state` Cxx pointer to a state
"""
function curand_init(seed, sequence, offset, state)
  icxx"""
    curand_init($seed, $sequence, $offset, $state);
  """
end

"""
   curand(state)

"""
curand(state) = icxx"curand($state);"

curand_uniform(::Type{Float32}, state) = icxx"curand_uniform($state);"
curand_uniform(::Type{Float64}, state) = icxx"curand_uniform_double($state);"

curand_normal(::Type{Float32}, state) = icxx"curand_normal($state);"
curand_normal(::Type{Float64}, state) = icxx"curand_normal_double($state);"

curand_log_normal(::Type{Float32}, state, mean, stddev) =
  icxx"curand_log_normal($state, $mean, $stddev);"
curand_log_normal(::Type{Float64}, state, mean, stddev) =
  icxx"curand_log_normal_double($state, $mean, $stddev);"

curand_poisson(state, lambda) = icxx"curand_poisson($state, $lambda);"

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Data buffer for pipes.
mutable struct Buffer
    data::Vector{UInt8}
    ptr::Ptr{UInt8}
    size::Int

    function Buffer(bufsize)
        data = Vector{UInt8}(undef, bufsize)
        return new(data, pointer(data), 0)
    end
end

Base.empty!(buffer::Buffer) = buffer.size = 0
Base.getindex(buffer::Buffer, i::Integer) = unsafe_load(buffer.ptr, i)
Base.setindex!(buffer::Buffer, v::UInt8, i::Integer) = unsafe_store!(buffer.ptr, v, i)
Base.firstindex(buffer::Buffer) = 1
Base.lastindex(buffer::Buffer) = buffer.size
Base.pointer(buffer::Buffer) = buffer.ptr
capacity(buffer::Buffer) = Int(pointer(buffer.data, lastindex(buffer.data) + 1) - buffer.ptr)

function consumed!(buffer::Buffer, n::Integer)
    @assert n ≤ buffer.size
    buffer.ptr += n
    buffer.size -= n
end

function read_to_buffer(io::IO, buffer::Buffer)
    offset = buffer.ptr - pointer(buffer.data)
    copyto!(buffer.data, 1, buffer.data, offset + 1, buffer.size)
    buffer.ptr = pointer(buffer.data)
    if !eof(io)
        n = min(bytesavailable(io), capacity(buffer) - buffer.size)
        unsafe_read(io, buffer.ptr + buffer.size, n)
        buffer.size += n
    end
    return
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
"""
```julia
calculate_tgrad(sys::AbstractSystem)
```

Calculate the time gradient of a system.

Returns a vector of [`Expression`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_tgrad end

"""
```julia
calculate_gradient(sys::AbstractSystem)
```

Calculate the gradient of a scalar system.

Returns a vector of [`Expression`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_gradient end

"""
```julia
calculate_jacobian(sys::AbstractSystem)
```

Calculate the jacobian matrix of a system.

Returns a matrix of [`Expression`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_jacobian end

"""
```julia
calculate_factorized_W(sys::AbstractSystem)
```

Calculate the factorized W-matrix of a system.

Returns a matrix of [`Expression`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_factorized_W end

"""
```julia
calculate_hessian(sys::AbstractSystem)
```

Calculate the hessian matrix of a scalar system.

Returns a matrix of [`Expression`](@ref) instances. The result from the first
call will be cached in the system object.
"""
function calculate_hessian end

"""
```julia
generate_tgrad(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
```

Generates a function for the time gradient of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_tgrad end

"""
```julia
generate_gradient(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
```

Generates a function for the gradient of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_gradient end

"""
```julia
generate_jacobian(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; sparse = false, kwargs...)
```

Generates a function for the jacobian matrix matrix of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_jacobian end

"""
```julia
generate_factorized_W(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; sparse = false, kwargs...)
```

Generates a function for the factorized W-matrix matrix of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_factorized_W end

"""
```julia
generate_hessian(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; sparse = false, kwargs...)
```

Generates a function for the hessian matrix matrix of a system. Extra arguments control
the arguments to the internal [`build_function`](@ref) call.
"""
function generate_hessian end

"""
```julia
generate_function(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
```

Generate a function to evaluate the system's equations.
"""
function generate_function end

function Base.getproperty(sys::AbstractSystem, name::Symbol)
    if name ∈ fieldnames(typeof(sys))
        return getfield(sys,name)
    elseif !isempty(sys.systems)
        i = findfirst(x->x.name==name,sys.systems)
        if i !== nothing
            return rename(sys.systems[i],renamespace(sys.name,name))
        end
    end
    i = findfirst(x->x.name==name,sys.states)
    if i !== nothing
        x = rename(sys.states[i],renamespace(sys.name,name))
        if :iv ∈ fieldnames(typeof(sys))
            return x(getfield(sys,:iv)())
        else
            return x()
        end
    end
    if :ps ∈ fieldnames(typeof(sys))
        i = findfirst(x->x.name==name,sys.ps)
        if i !== nothing
            return rename(sys.ps[i],renamespace(sys.name,name))()
        end
    end
    throw(error("Variable $name does not exist"))
end

renamespace(namespace,name) = Symbol(namespace,:₊,name)

function namespace_variables(sys::AbstractSystem)
    [rename(x,renamespace(sys.name,x.name)) for x in states(sys)]
end

function namespace_parameters(sys::AbstractSystem)
    [rename(x,renamespace(sys.name,x.name)) for x in parameters(sys)]
end

namespace_equations(sys::AbstractSystem) = namespace_equation.(equations(sys),sys.name,sys.iv.name)

function namespace_equation(eq::Equation,name,ivname)
    _lhs = namespace_operation(eq.lhs,name,ivname)
    _rhs = namespace_operation(eq.rhs,name,ivname)
    _lhs ~ _rhs
end

function namespace_operation(O::Operation,name,ivname)
    if O.op isa Variable && O.op.name != ivname
        Operation(rename(O.op,renamespace(name,O.op.name)),namespace_operation.(O.args,name,ivname))
    else
        Operation(O.op,namespace_operation.(O.args,name,ivname))
    end
end
namespace_operation(O::Constant,name,ivname) = O

independent_variable(sys::AbstractSystem) = sys.iv
states(sys::AbstractSystem) = isempty(sys.systems) ? sys.states : [sys.states;reduce(vcat,namespace_variables.(sys.systems))]
parameters(sys::AbstractSystem) = isempty(sys.systems) ? sys.ps : [sys.ps;reduce(vcat,namespace_parameters.(sys.systems))]

function equations(sys::AbstractSystem)
    isempty(sys.systems) ? sys.eqs : [sys.eqs;reduce(vcat,namespace_equations.(sys.systems))]
end

function states(sys::AbstractSystem,name::Symbol)
    x = sys.states[findfirst(x->x.name==name,sys.states)]
    rename(x,renamespace(sys.name,x.name))(sys.iv())
end

function parameters(sys::AbstractSystem,name::Symbol)
    x = sys.ps[findfirst(x->x.name==name,sys.ps)]
    rename(x,renamespace(sys.name,x.name))()
end

function states(sys::AbstractSystem,args...)
    name = last(args)
    extra_names = reduce(Symbol,[Symbol(:₊,x.name) for x in args[1:end-1]])
    newname = renamespace(extra_names,name)
    rename(x,renamespace(sys.name,newname))(sys.iv())
end

function parameters(sys::AbstractSystem,args...)
    name = last(args)
    extra_names = reduce(Symbol,[Symbol(:₊,x.name) for x in args[1:end-1]])
    newname = renamespace(extra_names,name)
    rename(x,renamespace(sys.name,newname))()
end

function islinear(sys::AbstractSystem)
    rhs = [eq.rhs for eq ∈ equations(sys)]

    iv = sys.iv()
    dvs = [dv(iv) for dv ∈ states(sys)]

    all(islinear(r, dvs) for r in rhs)
end

struct AbstractSysToExpr
    sys::AbstractSystem
    states::Vector{Variable}
end
AbstractSysToExpr(sys) = AbstractSysToExpr(sys,states(sys))
function (f::AbstractSysToExpr)(O::Operation)
    any(isequal(O), f.states) && return O.op.name  # variables
    if isa(O.op, Variable)
        isempty(O.args) && return O.op.name  # 0-ary parameters
        return build_expr(:call, Any[O.op.name; f.(O.args)])
    end
    return build_expr(:call, Any[O.op; f.(O.args)])
end
(f::AbstractSysToExpr)(x) = convert(Expr, x)


text:
using MPSKit,TensorKit,Test,OptimKit,MPSKitModels,TestExtras

println("------------------------------------")
println("|     States                       |")
println("------------------------------------")
@timedtestset "FiniteMPS ($D,$d,$elt)" for (D,d,elt) in [
        (ComplexSpace(10),ComplexSpace(2),ComplexF64),
        (Rep[SU₂](1=>1,0=>3),Rep[SU₂](0=>1)*Rep[SU₂](0=>1),ComplexF32)
        ]

    ts = FiniteMPS(rand,elt,rand(3:20),d,D);

    ovl = dot(ts,ts);

    @test ovl ≈ norm(ts.AC[1])^2

    for i in 1:length(ts)
        @test ts.AC[i] ≈ ts.AL[i]*ts.CR[i]
        @test ts.AC[i] ≈ MPSKit._permute_front(ts.CR[i-1]*MPSKit._permute_tail(ts.AR[i]))
    end

    @test elt == eltype(eltype(ts))

    ts = ts*3
    @test ovl*9 ≈ norm(ts)^2
    ts = 3*ts
    @test ovl*9*9 ≈ norm(ts)^2
end

@timedtestset "InfiniteMPS ($D,$d,$elt)" for (D,d,elt) in [
        (ComplexSpace(10),ComplexSpace(2),ComplexF64),
        (Rep[U₁](1=>3),Rep[U₁](0=>1),ComplexF64)
        ]
    tol = Float64(eps(real(elt))*100);

    ts = InfiniteMPS([TensorMap(rand,elt,D*d,D),TensorMap(rand,elt,D*d,D)],tol = tol);

    for i in 1:length(ts)
        @tensor difference[-1,-2,-3] := ts.AL[i][-1,-2,1]*ts.CR[i][1,-3]-ts.CR[i-1][-1,1]*ts.AR[i][1,-2,-3];
        @test norm(difference,Inf) < tol*10;

        @test transfer_left(l_LL(ts,i),ts.AL[i],ts.AL[i]) ≈ l_LL(ts,i+1)
        @test transfer_left(l_LR(ts,i),ts.AL[i],ts.AR[i]) ≈ l_LR(ts,i+1)
        @test transfer_left(l_RL(ts,i),ts.AR[i],ts.AL[i]) ≈ l_RL(ts,i+1)
        @test transfer_left(l_RR(ts,i),ts.AR[i],ts.AR[i]) ≈ l_RR(ts,i+1)

        @test transfer_right(r_LL(ts,i),ts.AL[i],ts.AL[i]) ≈ r_LL(ts,i+1)
        @test transfer_right(r_LR(ts,i),ts.AL[i],ts.AR[i]) ≈ r_LR(ts,i+1)
        @test transfer_right(r_RL(ts,i),ts.AR[i],ts.AL[i]) ≈ r_RL(ts,i+1)
        @test transfer_right(r_RR(ts,i),ts.AR[i],ts.AR[i]) ≈ r_RR(ts,i+1)
    end
end

@timedtestset "MPSMultiline ($D,$d,$elt)" for (D,d,elt) in [
        (ComplexSpace(10),ComplexSpace(2),ComplexF64),
        (Rep[U₁](1=>3),Rep[U₁](0=>1),ComplexF32)
        ]

    tol = Float64(eps(real(elt))*100);
    ts = MPSMultiline([TensorMap(rand,elt,D*d,D) TensorMap(rand,elt,D*d,D);TensorMap(rand,elt,D*d,D) TensorMap(rand,elt,D*d,D)],tol = tol);

    for i = 1:size(ts,1), j = 1:size(ts,2)
        @tensor difference[-1,-2,-3] := ts.AL[i,j][-1,-2,1]*ts.CR[i,j][1,-3]-ts.CR[i,j-1][-1,1]*ts.AR[i,j][1,-2,-3];
        @test norm(difference,Inf) < tol*10;

        @test transfer_left(l_LL(ts,i,j),ts.AL[i,j],ts.AL[i,j]) ≈ l_LL(ts,i,j+1)
        @test transfer_left(l_LR(ts,i,j),ts.AL[i,j],ts.AR[i,j]) ≈ l_LR(ts,i,j+1)
        @test transfer_left(l_RL(ts,i,j),ts.AR[i,j],ts.AL[i,j]) ≈ l_RL(ts,i,j+1)
        @test transfer_left(l_RR(ts,i,j),ts.AR[i,j],ts.AR[i,j]) ≈ l_RR(ts,i,j+1)

        @test transfer_right(r_LL(ts,i,j),ts.AL[i,j],ts.AL[i,j]) ≈ r_LL(ts,i,j+1)
        @test transfer_right(r_LR(ts,i,j),ts.AL[i,j],ts.AR[i,j]) ≈ r_LR(ts,i,j+1)
        @test transfer_right(r_RL(ts,i,j),ts.AR[i,j],ts.AL[i,j]) ≈ r_RL(ts,i,j+1)
        @test transfer_right(r_RR(ts,i,j),ts.AR[i,j],ts.AR[i,j]) ≈ r_RR(ts,i,j+1)
    end
end

@timedtestset "MPSComoving" begin
    ham = nonsym_ising_ham(lambda=4.0);
    (gs,_,_) = find_groundstate(InfiniteMPS([ℂ^2],[ℂ^10]),ham,Vumps(verbose=false));

    #constructor 1 - give it a plain array of tensors
    window_1 = MPSComoving(gs,copy.([gs.AC[1];[gs.AR[i] for i in 2:10]]),gs);

    #constructor 2 - used to take a "slice" from an infinite mps
    window_2 = MPSComoving(gs,10);

    # we should logically have that window_1 approximates window_2
    ovl = dot(window_1,window_2)
    @test ovl ≈ 1 atol=1e-8

    #constructor 3 - random initial tensors
    window = MPSComoving(rand,ComplexF64,10,ℂ^2,ℂ^10,gs,gs)
    normalize!(window);

    for i in 1:length(window)
        @test window.AC[i] ≈ window.AL[i]*window.CR[i]
        @test window.AC[i] ≈ MPSKit._permute_front(window.CR[i-1]*MPSKit._permute_tail(window.AR[i]))
    end

    @test norm(window) ≈ 1
    window = window*3
    @test 9 ≈ norm(window)^2
    window = 3*window
    @test 9*9 ≈ norm(window)^2
    normalize!(window)

    e1 = expectation_value(window,ham);

    v1 = variance(window,ham)
    (window,envs,_) = find_groundstate(window,ham,Dmrg(verbose=false));
    v2 = variance(window,ham)

    e2 = expectation_value(window,ham);

    @test v2<v1
    @test real(e2[2]) ≤ real(e1[2])

    (window,envs) = timestep(window,ham,0.1,Tdvp2(),envs)
    (window,envs) = timestep(window,ham,0.1,Tdvp(),envs)

    e3 = expectation_value(window,ham);

    @test e2[1] ≈ e3[1]
    @test e2[2] ≈ e3[2]
end

@timedtestset "Quasiparticle state" begin
    @timedtestset "Finite" for (th,D,d) in [
        (nonsym_ising_ham(),ComplexSpace(10),ComplexSpace(2)),
        (su2_xxx_ham(spin=1),Rep[SU₂](1=>1,0=>3),Rep[SU₂](1=>1))
        ]

        ts = FiniteMPS(rand,ComplexF64,rand(4:20),d,D);
        normalize!(ts);

        #rand_quasiparticle is a private non-exported function
        qst1 = MPSKit.LeftGaugedQP(rand,ts);
        qst2 = MPSKit.LeftGaugedQP(rand,ts);

        @test norm(axpy!(1,qst1,copy(qst2))) ≤ norm(qst1) + norm(qst2)
        @test norm(qst1)*3 ≈ norm(qst1*3)

        normalize!(qst1);

        qst1_f = convert(FiniteMPS,qst1);
        qst2_f = convert(FiniteMPS,qst2);

        ovl_f = dot(qst1_f,qst2_f)
        ovl_q = dot(qst1,qst2)
        @test ovl_f ≈ ovl_q atol=1e-5
        @test norm(qst1_f) ≈ norm(qst1) atol=1e-5

        ev_f = sum(expectation_value(qst1_f,th)-expectation_value(ts,th))
        ev_q = dot(qst1,effective_excitation_hamiltonian(th,qst1));
        @test ev_f ≈ ev_q atol=1e-5
    end

    @timedtestset "Infinite" for (th,D,d) in [
        (nonsym_ising_ham(),ComplexSpace(10),ComplexSpace(2)),
        (su2_xxx_ham(spin=1),Rep[SU₂](1=>1,0=>3),Rep[SU₂](1=>1))
        ]

        period = rand(1:4);
        ts = InfiniteMPS(fill(d,period),fill(D,period));

        #rand_quasiparticle is a private non-exported function
        qst1 = MPSKit.LeftGaugedQP(rand,ts);
        qst2 = MPSKit.LeftGaugedQP(rand,ts);

        @test norm(axpy!(1,qst1,copy(qst2))) ≤ norm(qst1) + norm(qst2)
        @test norm(qst1)*3 ≈ norm(qst1*3)

        @test dot(qst1,convert(MPSKit.LeftGaugedQP,convert(MPSKit.RightGaugedQP,qst1))) ≈ dot(qst1,qst1) atol=1e-10
    end

end

println("------------------------------------")
println("|     Operators                    |")
println("------------------------------------")

@timedtestset "mpoham $((pspace,Dspace))" for (pspace,Dspace) in [(ℂ^4,ℂ^10),
        (Rep[U₁](0=>2),Rep[U₁]((0=>20))),
        (Rep[SU₂](1=>1),Rep[SU₂](1//2=>10,3//2=>5,5//2=>1))]

    #generate a 1-2-3 body interaction
    n = TensorMap(rand,ComplexF64,pspace,pspace); n+= n';
    nn = TensorMap(rand,ComplexF64,pspace*pspace,pspace*pspace); nn+=nn';
    nnn = TensorMap(rand,ComplexF64,pspace*pspace*pspace,pspace*pspace*pspace); nnn+=nnn';

    #can you pass in a proper mpo?
    identity = complex(isomorphism(oneunit(pspace)*pspace,oneunit(pspace)*pspace));
    mpoified = MPSKit.decompose_localmpo(MPSKit.add_util_leg(nnn));
    d3 = Array{Union{Missing,typeof(identity)},3}(missing,1,4,4);
    d3[1,1,1] = identity;
    d3[1,end,end] = identity;
    d3[1,1,2] = mpoified[1];
    d3[1,2,3] = mpoified[2];
    d3[1,3,4] = mpoified[3];
    h1 = MPOHamiltonian(d3);

    #¢an you pass in the actual hamiltonian?
    h2 = MPOHamiltonian(nn);

    #can you generate a hamiltonian using only onsite interactions?
    d1 = Array{Any,3}(missing,2,3,3);
    d1[1,1,1] = 1; d1[1,end,end] = 1;
    d1[1,1,2] = n; d1[1,2,end] = n;
    d1[2,1,1] = 1; d1[2,end,end] = 1;
    d1[2,1,2] = n; d1[2,2,end] = n;
    h3 = MPOHamiltonian(d1);

    #make a teststate to measure expectation values for
    ts1 = InfiniteMPS([pspace],[Dspace]);
    ts2 = InfiniteMPS([pspace,pspace],[Dspace,Dspace]);

    e1 = expectation_value(ts1,h1);
    e2 = expectation_value(ts1,h2);

    h1 = 2*h1-[1];
    @test e1[1]*2-1 ≈ expectation_value(ts1,h1)[1] atol=1e-10

    h1 = h1 + h2;

    @test e1[1]*2+e2[1]-1 ≈ expectation_value(ts1,h1)[1] atol=1e-10

    h1 = repeat(h1,2);

    e1 = expectation_value(ts2,h1);
    e3 = expectation_value(ts2,h3);

    @test e1+e3 ≈ expectation_value(ts2,h1+h3) atol=1e-10

    h4 = h1+h3;
    h4 = h4*h4;
    @test real(sum(expectation_value(ts2,h4)))>=0;
end

@timedtestset "InfiniteMPO"  for ham in (nonsym_ising_ham(),su2_xxx_ham(spin=1))
    physical_space = ham.pspaces[1];
    ou = oneunit(physical_space);

    ts = InfiniteMPS([physical_space],[fuse(ou*physical_space)]);

    W = make_time_mpo(ham,1im*0.5,WII());

    @test abs(dot(W*(W*ts),(W*W)*ts))≈1.0 atol=1e-10
end

println("------------------------------------")
println("|     Algorithms                   |")
println("------------------------------------")

@timedtestset "find_groundstate $(ind)" for (ind,(state,alg,ham)) in enumerate([
        (InfiniteMPS([ℂ^2],[ℂ^10]),Vumps(tol_galerkin=1e-8,verbose=false),nonsym_ising_ham(lambda=2.0)),
        (InfiniteMPS([ℂ^2],[ℂ^10]), GradientGrassmann(tol=1e-8, verbosity=0), nonsym_ising_ham(lambda=2.0)),
        (InfiniteMPS([ℂ^2],[ℂ^10]), GradientGrassmann(method=LBFGS(6; gradtol=1e-8, verbosity=0)), nonsym_ising_ham(lambda=2.0)),
        (InfiniteMPS([ℂ^2],[ℂ^10]),Idmrg1(tol_galerkin=1e-8,maxiter=400,verbose=false),nonsym_ising_ham(lambda=2.0)),
        (InfiniteMPS([ℂ^2,ℂ^2],[ℂ^10,ℂ^10]),Idmrg2(trscheme = truncdim(12),tol_galerkin=1e-8,maxiter=400,verbose=false),repeat(nonsym_ising_ham(lambda=2.0),2)),
        (InfiniteMPS([ℂ^2], [ℂ^10]), Vumps(tol_galerkin=1e-5,verbose=false)&GradientGrassmann(tol=1e-8, verbosity=0), nonsym_ising_ham(lambda=2.0)),
        (FiniteMPS(rand,ComplexF64,10,ℂ^2,ℂ^10),Dmrg2(verbose=false,trscheme=truncdim(10)),nonsym_ising_ham(lambda=2.0)),
        (FiniteMPS(rand,ComplexF64,10,ℂ^2,ℂ^10),Dmrg(verbose=false),nonsym_ising_ham(lambda=2.0)),
        (FiniteMPS(rand,ComplexF64,10,ℂ^2,ℂ^10),GradientGrassmann(verbosity=0),nonsym_ising_ham(lambda=2.0))
        ])

    v1 = variance(state,ham);
    (ts,envs,delta) =  find_groundstate(state,ham,alg)
    v2 = variance(ts,ham);

    @test v1 > v2
    @test sum(delta) < 1e-6
    @test v2 < 1e-2 #is the ground state variance relatively low?
end

@timedtestset "timestep $(ind)" for (ind,(state,alg,opp)) in enumerate([
    (FiniteMPS(fill(TensorMap(rand,ComplexF64,ℂ^1*ℂ^2,ℂ^1),5)),Tdvp(),nonsym_xxz_ham(spin=1//2)),
    (FiniteMPS(fill(TensorMap(rand,ComplexF64,ℂ^1*ℂ^2,ℂ^1),7)),Tdvp2(),nonsym_xxz_ham(spin=1//2)),
    (InfiniteMPS([ℂ^3,ℂ^3],[ℂ^50,ℂ^50]),Tdvp(),repeat(nonsym_xxz_ham(spin=1),2))
    ])

    edens = expectation_value(state,opp);

    (state,_) = timestep(state,opp,rand()/10,alg)

    @test sum(expectation_value(state,opp)) ≈ sum(edens) atol = 1e-2
end

@timedtestset "leading_boundary $(ind)" for (ind,alg) in enumerate([
        Vumps(tol_galerkin=1e-5,verbose=false);
        GradientGrassmann(verbosity=0)])

    mpo = nonsym_ising_mpo();
    state = InfiniteMPS([ℂ^2],[ℂ^10]);
    (state,envs) = leading_boundary(state,mpo,alg);
    (state,envs) = changebonds(state,mpo,OptimalExpand(trscheme=truncdim(3)),envs)
    (state,envs) = leading_boundary(state,mpo,alg);

    @test dim(space(state.AL[1,1],1)) == 13
    @test expectation_value(state,envs)[1,1] ≈ 2.5337 atol=1e-3
end

@timedtestset "quasiparticle_excitation" begin
    @timedtestset "infinite" begin
        th = nonsym_xxz_ham()
        ts = InfiniteMPS([ℂ^3],[ℂ^48]);
        (ts,envs,_) = find_groundstate(ts,th,Vumps(maxiter=400,verbose=false));
        (energies,Bs) = excitations(th,QuasiparticleAnsatz(),Float64(pi),ts,envs);
        @test energies[1] ≈ 0.41047925 atol=1e-4
        @test variance(Bs[1],th) < 1e-8
    end

    @timedtestset "finite" begin
        th = nonsym_ising_ham()
        ts = InfiniteMPS([ℂ^2],[ℂ^12]);
        (ts,envs,_) = find_groundstate(ts,th,Vumps(maxiter=400,verbose=false));
        (energies,Bs) = excitations(th,QuasiparticleAnsatz(),0.0,ts,envs);
        inf_en = energies[1];

        fin_en = map([30,20,10]) do len
            ts = FiniteMPS(rand,ComplexF64,len,ℂ^2,ℂ^12)
            (ts,envs,_) = find_groundstate(ts,th,Dmrg(verbose=false));

            #find energy with quasiparticle ansatz
            (energies_QP,Bs) = excitations(th,QuasiparticleAnsatz(),ts,envs);
            @test variance(Bs[1],th)<1e-8

            #find energy with normal dmrg
            (energies_dm,_) = excitations(th,FiniteExcited(gsalg=Dmrg(verbose=false)),ts);
            @test energies_dm[1] ≈ energies_QP[1]+sum(expectation_value(ts,th,envs)) atol=1e-4

            energies_QP[1]
        end

        @test issorted(abs.(fin_en.-inf_en))
    end
end

@timedtestset "changebonds $((pspace,Dspace))" for (pspace,Dspace) in [(ℂ^4,ℂ^10),
        (Rep[SU₂](1=>1),Rep[SU₂](0=>10,1=>5,2=>1))]

    @timedtestset "mpo" begin
        #random nn interaction
        nn = TensorMap(rand,ComplexF64,pspace*pspace,pspace*pspace);
        nn += nn';

        mpo1 = periodic_boundary_conditions(make_time_mpo(MPOHamiltonian(nn),0.1,WII()),10);
        mpo2 = changebonds(mpo1,SvdCut(trscheme = truncdim(5)));

        @test dim(space(mpo2[5],1)) < dim(space(mpo1[5],1))
    end

    @timedtestset "infinite mps" begin
        #random nn interaction
        nn = TensorMap(rand,ComplexF64,pspace*pspace,pspace*pspace);
        nn += nn';

        state = InfiniteMPS([pspace,pspace],[Dspace,Dspace]);

        state_re = changebonds(state,RandExpand(trscheme = truncdim(dim(Dspace)*dim(Dspace))));
        @test dot(state,state_re) ≈ 1 atol=1e-8

        (state_oe,_) = changebonds(state,MPOHamiltonian(nn),OptimalExpand(trscheme = truncdim(dim(Dspace)*dim(Dspace))));
        @test dot(state,state_oe) ≈ 1 atol=1e-8

        (state_vs,_) = changebonds(state,MPOHamiltonian(nn),VumpsSvdCut(trscheme=notrunc()));
        @test dim(virtualspace(state,1)) < dim(virtualspace(state_vs,1))

        state_vs_tr = changebonds(state_vs,SvdCut(trscheme = truncdim(dim(Dspace))));
        @test dim(virtualspace(state_vs_tr,1)) < dim(virtualspace(state_vs,1))
    end

    @timedtestset "finite mps" begin
        #random nn interaction
        nn = TensorMap(rand,ComplexF64,pspace*pspace,pspace*pspace);
        nn += nn';

        state = FiniteMPS(10,pspace,Dspace);

        state_re = changebonds(state,RandExpand(trscheme = truncdim(dim(Dspace)*dim(Dspace))));
        @test dot(state,state_re) ≈ 1 atol=1e-8

        (state_oe,_) = changebonds(state,MPOHamiltonian(nn),OptimalExpand(trscheme = truncdim(dim(Dspace)*dim(Dspace))));
        @test dot(state,state_oe) ≈ 1 atol=1e-8

        state_tr = changebonds(state_oe,SvdCut(trscheme = truncdim(dim(Dspace))));
        @test dim(virtualspace(state_tr,5)) < dim(virtualspace(state_oe,5))
    end

    @timedtestset "MPSMultiline" begin
        o = TensorMap(rand,ComplexF64,pspace*pspace,pspace*pspace);
        mpo = MPOMultiline(o);

        t = TensorMap(rand,ComplexF64,Dspace*pspace,Dspace);
        state = MPSMultiline(fill(t,1,1));

        state_re = changebonds(state,RandExpand(trscheme = truncdim(dim(Dspace)*dim(Dspace))));
        @test dot(state,state_re) ≈ 1 atol=1e-8

        (state_oe,_) = changebonds(state,mpo,OptimalExpand(trscheme = truncdim(dim(Dspace)*dim(Dspace))));
        @test dot(state,state_oe) ≈ 1 atol=1e-8

        state_tr = changebonds(state_oe,SvdCut(trscheme = truncdim(dim(Dspace))));
        @test dim(virtualspace(state_tr,1,1)) < dim(virtualspace(state_oe,1,1))

    end
end

@timedtestset "dynamicaldmrg" begin
    ham = nonsym_ising_ham(lambda=4.0);
    (gs,_,_) = find_groundstate(InfiniteMPS([ℂ^2],[ℂ^10]),ham,Vumps(verbose=false));
    window = MPSComoving(gs,copy.([gs.AC[1];[gs.AR[i] for i in 2:10]]),gs);

    szd = TensorMap([1 0;0 -1],ℂ^2,ℂ^2);
    @test expectation_value(gs,szd)[1] ≈ expectation_value(window,szd)[1] atol=1e-10

    polepos = expectation_value(gs,ham,10)
    @test polepos ≈ expectation_value(window,ham)[2]

    vals = (-0.5:0.2:0.5).+polepos
    eta = 0.3im;

    predicted = [1/(v+eta-polepos) for v in vals];

    data = similar(predicted);
    for (i,v) in enumerate(vals)
        (data[i],_) = dynamicaldmrg(window,v+eta,ham,verbose=false)
    end

    @test data ≈ predicted atol=1e-8
end

@timedtestset "fidelity susceptibility" begin
    lambs = [1.05,2.0,4.0]

    for l in lambs

        (X,Y,Z) = nonsym_spintensors(1//2).*2;

        @tensor ZZ_data[-1 -2;-3 -4] := Z[-1 -3]*Z[-2 -4]
        @tensor X_data[-1;-2] := X[-1 -2]

        ZZham = MPOHamiltonian(ZZ_data);
        Xham = MPOHamiltonian(X_data);

        th = ZZham + l*Xham;

        ts = InfiniteMPS([ℂ^2],[ℂ^20]);
        (ts,envs,_) = find_groundstate(ts,th,Vumps(maxiter=1000,verbose=false));

        #test if the infinite fid sus approximates the analytical one
        num_sus = fidelity_susceptibility(ts,th,[Xham],envs,maxiter=10);
        ana_sus = abs.(1/(16*l^2*(l^2-1)));
        @test ana_sus ≈ num_sus[1,1] atol=1e-2

        #test if the finite fid sus approximates the analytical one with increasing system size
        fin_en = map([30,20,10]) do len
            ts = FiniteMPS(rand,ComplexF64,len,ℂ^2,ℂ^20)
            (ts,envs,_) = find_groundstate(ts,th,Dmrg(verbose=false));
            num_sus = fidelity_susceptibility(ts,th,[Xham],envs,maxiter=10);
            num_sus[1,1]/len
        end
        @test issorted(abs.(fin_en.-ana_sus))
    end

end

#stub tests
@timedtestset "correlation length / entropy" begin

    st = InfiniteMPS([ℂ^2],[ℂ^10]);
    th = nonsym_ising_ham();
    (st,_) = find_groundstate(st,th,Vumps(verbose=false))
    len_crit = correlation_length(st)[1]
    entrop_crit = entropy(st);

    th = nonsym_ising_ham(lambda=4);
    (st,_) = find_groundstate(st,th,Vumps(verbose=false))
    len_gapped = correlation_length(st)[1]
    entrop_gapped = entropy(st);

    @test len_crit > len_gapped;
    @test entrop_crit > entrop_gapped;
end

@timedtestset "expectation value" begin
    st = InfiniteMPS([ℂ^2],[ℂ^10]);
    th = nonsym_ising_ham(lambda=4);
    (st,_) = find_groundstate(st,th,Vumps(verbose=false))

    sz_mpo =TensorMap([1.0 0;0 -1],ℂ^1*ℂ^2,ℂ^1*ℂ^2)
    sz =TensorMap([1.0 0;0 -1],ℂ^2,ℂ^2)
    @tensor szsz[-1 -2;-3 -4]:=sz[-1 -3]*sz[-2 -4]

    @test isapprox(expectation_value(st, [sz_mpo], 1) ,  expectation_value(st,sz,1) ,atol = 1e-2)
    @test isapprox(expectation_value(st, [sz_mpo,sz_mpo], 1) , expectation_value(st,szsz,1) ,atol = 1e-2)
    @test isapprox(expectation_value(st, [sz_mpo,sz_mpo], 2) , expectation_value(st,szsz,1) ,atol = 1e-2)
end

@timedtestset "approximate" begin
    @timedtestset "mpo * infinite ≈ infinite" begin
        st = InfiniteMPS([ℂ^2],[ℂ^10]);
        th = nonsym_ising_ham(lambda=4);

        dt = 1e-3;
        W1 = make_time_mpo(th,dt,WI());
        W2 = make_time_mpo(th,dt,WII());

        (st1,_) = approximate(st,(W1,st),Vomps(verbose=false));
        (st2,_) = approximate(st,(W2,st),Vomps(verbose=false));
        (st3,_) = timestep(st,th,dt,Tdvp());
        st4 = changebonds(W1*st,SvdCut(trscheme=truncdim(10)))

        @test abs(dot(st1,st3)) ≈ 1.0 atol = dt
        @test abs(dot(st2,st3)) ≈ 1.0 atol = dt
        @test abs(dot(st4,st3)) ≈ 1.0 atol = dt

        nW1 = changebonds(W1,SvdCut(trscheme=truncerr(dt))); #this should be a trivial mpo now
        @test dim(space(nW1.opp[1,1],1)) == 1
    end

    @timedtestset "finitemps1 ≈ finitemps2" for alg in [Dmrg(verbose=false),Dmrg2(verbose=false,trscheme=truncdim(10))]
        a = FiniteMPS(10,ℂ^2,ℂ^10);
        b = FiniteMPS(10,ℂ^2,ℂ^20);

        before = abs(dot(a,b));

        a = first(approximate(a,b,alg));

        after = abs(dot(a,b));

        @test before < after
    end
end

@timedtestset "periodic boundary conditions" begin
    len = 10;

    #impose periodic boundary conditions on the hamiltonian (cirkel size 10)
    th = nonsym_ising_ham();
    th = periodic_boundary_conditions(th,len);

    ts = FiniteMPS(len,ℂ^2,ℂ^10);

    (gs,envs) = find_groundstate(ts,th,Dmrg(verbose=false));

    #translation mpo:
    @tensor bulk[-1 -2;-3 -4] := isomorphism(ℂ^2,ℂ^2)[-2,-3]*isomorphism(ℂ^2,ℂ^2)[-1,-4];
    translation = periodic_boundary_conditions(InfiniteMPO(bulk),len);

    #the groundstate should be translation invariant:
    ut = Tensor(ones,ℂ^1);
    @tensor leftstart[-1 -2;-3] := l_LL(gs)[-1,-3]*conj(ut[-2]);
    v = transfer_left(leftstart,translation[:],[gs.AC[1];gs.AR[2:end]],[gs.AC[1];gs.AR[2:end]])
    expval = @tensor v[1,2,3]*r_RR(gs)[3,1]*ut[2]

    @test expval ≈ 1 atol=1e-5

    (energies,values) = exact_diagonalization(th,which=:SR);
    @test energies[1] ≈ sum(expectation_value(gs,th)) atol=1e-5
end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 1f6d20c6-2a06-11ec-30b2-f5462ac0ec76
begin
	using Pkg
	Pkg.activate("../../")
	using Revise
end

# ╔═╡ 46cd6d61-7721-49be-909f-7c411a5562ff
using ImgProcMic, PlutoUI

# ╔═╡ a0b4af53-ded5-4d30-9d2b-7c5d1c2bff80
begin
	using ColorSchemes 	# for different colors maps
	using Plots 		# for heatmap
	using TestImages 	# for testimage
	using ImageShow 	# for better rendering of inline images in Pluto
	using Colors 		# for Gray
end

# ╔═╡ a7007d37-6f30-462b-8be3-70999c980c60
using PlutoTest

# ╔═╡ d9d77790-30d1-4ca2-83ca-d974b10a9431
begin
	using PoissonRandom # for poisson noise (`pois_rand`)
	using Statistics # for mean, var
end

# ╔═╡ 7a0e442c-9f32-4210-9483-63e0dc9fbe49
using FourierTools

# ╔═╡ dd4f1881-4570-4972-bd40-a4123ad4b6b1
# for fitting
using LsqFit

# ╔═╡ 85a15a0e-37d4-4b81-986c-d8a2fe4185c7
md"# Homework 02
In this homework we work on images, colormaps and sensor calibration.
"

# ╔═╡ 375801bb-7b43-4983-a0ca-211d106ab11c
begin
	# load a image and shift it to negative values as well
	img = Float64.(testimage("moonsurface")) .- 0.5;
	histogram(img[:], xlabel="intensity", ylabel="occurence", title="Histogram Plot", legend=nothing)
end

# ╔═╡ e2768bdd-3346-4f13-88a8-ab8a43094761
md"## 1. Colormaps
Very often in microscopy, we look at grayscale images.
We already learnt that we can use `Gray.(img)` for that.
However, that is only a very simple viewer.
In general, there are much more sophisticated ones like `View5D.jl` or `Napari.jl`.

As you can see below, `Gray` fails pretty quickly and also does not make use of our human color vision.
"

# ╔═╡ cea22309-272b-4696-850a-f49ad3120345
Gray.(img)

# ╔═╡ c66329a6-84f4-4e0e-87a1-8a62f97ba9f7
md" ### 1.1 Heatmap
`img` contains also negative values which `Gray` does not support by default.
A nicer tool is heatmap.
"

# ╔═╡ 450d3dea-cb97-482b-aaab-bb2919e4ee9e
heatmap(img)

# ╔═╡ 41f64a18-bd0d-49a4-8347-36818b45a576
md"## 1.2 Task - Equal Intensity Color Map
We now want to create a RGB color map which always has the same intensity.

Hence, `green_value + red_value + blue_value = const` for all different inputs.

"

# ╔═╡ f6379b6e-7030-4967-afea-13494ccad294
"""
	to_equal_intensity_tuple(value::T) where T

Returns for a given scalar input a tuple of equal intensity.
There is no unique solution, we only require that the sum is always 1.
For example, below would be a valid solution.

## Examples
```julia-repl
julia> to_equal_intensity_tuple(1.0)
(0.04465819873852073, 0.33333333333333326, 0.6220084679281465)

julia> to_equal_intensity_tuple(0.0)
(0.044658198738520505, 0.3333333333333333, 0.6220084679281462)

julia> to_equal_intensity_tuple(0.3f0)
(0.26402956f0, 0.6503522f0, 0.08561832f0)
```
"""
function to_equal_intensity_tuple(value::T) where T
	# with the function signature we can access T (which is the type of `value`)
	# 0 has type Int
	# 1.0 has type Float64
	# 1.0f0 has type Float32
	if value < 0 || value > 1
		error("value must be within [0,1]")
	end

	# fix those three functions such that the property is achieved
	# maybe think about periodic functions ;)
	red_map(x) = x # TODO
	green_map(x) = x # TODO
	blue_map(x) = x # TODO

	return (3.0, 4f0, 2) # TODO, return the correct tuple!
end

# ╔═╡ 2b3c2ab1-9bec-43f5-a008-795ac354e080
md"#### Maybe a visualization helps
"

# ╔═╡ df017bc1-a29c-4bbc-80c3-bc2466cd95f5
begin
	unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))
	x = 0:0.01:1
	r,g,b = unzip(to_equal_intensity_tuple.(x))
	plot(x, r, color=:red, label="red")
	plot!(x, g, color=:green, label="green")
	plot!(x, b, color=:blue, label="blue")
end

# ╔═╡ f05ca6de-b46c-4400-a141-c662a1a16d7c
md"## 1.2 Test"

# ╔═╡ f1d7bf0a-6806-45fb-8213-313416017c18
# cyclic
PlutoTest.@test all(to_equal_intensity_tuple(0f0) .≈ to_equal_intensity_tuple(1f0))

# ╔═╡ 18bb1f0f-5b80-43d5-bc8a-526947a70f2c
PlutoTest.@test to_equal_intensity_tuple(0.1) != to_equal_intensity_tuple(0.2)

# ╔═╡ 76e2f439-1dff-4f2a-9c05-14e10a13fda6
# intensity is 1
PlutoTest.@test 1 ≈ sum(to_equal_intensity_tuple(rand(0f0:1f0)))

# ╔═╡ 5315fbfd-8ade-46d5-9851-f3ee98c349f7
# test for many values
PlutoTest.@test all(1 .≈ sum.(to_equal_intensity_tuple.(rand(0.0:1.0, 100))))

# ╔═╡ 9d4337c8-6cd1-4112-9477-4a40d13fe7e3
md"##### Check if the output is always 0<=output<=1

_Side Note: it took quite a while to engineer this line :D_
"

# ╔═╡ c0debce4-d81b-4d5d-8907-53d567cf664a
PlutoTest.@test foldl(&, foldl.(Ref((acc, x) -> acc || 0<=x<=1), to_equal_intensity_tuple.(rand(0.0:1.0, 100)), init=false))

# ╔═╡ 4456a207-7f56-4900-95dc-77fe7de31c5d
md"#### Type Stability

The test might be wrong.
Check what happens with expresions like `2π / 3`. Which type do they have?
How can we convert them to different types?
"

# ╔═╡ a6fbfa34-33c6-4f09-9dcd-773d92a4db0d
PlutoTest.@test NTuple{3, Float32} == typeof(to_equal_intensity_tuple(rand(Float32.(0.0:1.0))))

# ╔═╡ b70be939-7f1d-4e80-9d7f-cc2947d57280
PlutoTest.@test NTuple{3, Float64} == typeof(to_equal_intensity_tuple(rand(Float64.(0.0:1.0))))

# ╔═╡ e9cc6238-5ad9-4a50-93a7-f7065648574a
md"## 1.3 Register Color Map
We are able to create valid outputs, but we still need to pass that to the colormap mechanism.
We need to register the colormap so that we can load it with `heatmap`.

You don't need to know the details, but if you'd like understand, you find those in [ColorSchemes.jl](https://juliagraphics.github.io/ColorSchemes.jl/stable/basics/#Make-your-own-colorscheme).
"

# ╔═╡ 633d295e-cc1f-4864-a36b-11905bc47a50
begin
	# register that with the color map mechanism
	function create_equal_intensity_colormap()
		values = range(0, 1, length=128)
		return map(x -> RGB(x...), to_equal_intensity_tuple.(values))
	end

	begin
		my_equal_intensity_colormap = create_equal_intensity_colormap()
		loadcolorscheme(:my_equal_intensity_colormap, my_equal_intensity_colormap)
	end
end

# ╔═╡ 6b65dc02-9a84-43c2-9bb5-32104d23f580
heatmap(img, c=:my_equal_intensity_colormap)

# ╔═╡ 160051e5-56b9-413a-a78d-b371a85ff18c
md" ## 1.4 Task
Now we want to create a colormap, which maps values between 0 and 0.5
to blue values and all values above 0.5  to red.
Don't let you confuse by negativ and positive intensities.
`loadcolorschemes` needs the values within the interval [0, 1] and will automatically scale it to the image.

"

# ╔═╡ 1483db2c-6374-4c93-ade4-db8894071ce7
"""
	to_negativ_positive_tuple(value::T) where T

## Examples
```
julia> to_negative_positive_tuple(0.0)
(0.0, 0.0, 1.0)

julia> to_negative_positive_tuple(0.25)
(0.5, 0.5, 1.0)

julia> to_negative_positive_tuple(0.5)
(1.0, 1.0, 1.0)

julia> to_negative_positive_tuple(0.75)
(1.0, 0.5, 0.5)

julia> to_negative_positive_tuple(1.0)
(1.0, 0.0, 0.0)
```
"""
function to_negative_positive_tuple(value::T) where T
	# TODO build the function in analogy to to_equal_intensity_tuple

	# TODO
	# TODO

	# TODO
	# TODO

	# TODO
	# TODO

	# TODO
	# TODO

	return (0, 0, 0) # TODO
end

# ╔═╡ 0afb1220-1fdf-447d-8e79-b1ad7c274dd3
md"
Try to show again a heatmap image. The procedure is exactly the same above, except that you need to exchange the core part of the color generation."

# ╔═╡ f6e5d3c8-172e-4fb6-89ce-a1827a6ddd6c
begin
	# try to get :my_negative_positive_colormap working
	# see above what to do!

	# TODO
	# TODO
	# TODO
end

# ╔═╡ bd71f8c1-e1e8-4f72-9ad8-8714e9e215f0
heatmap(img, c=:my_negative_positive_colormap)

# ╔═╡ 45d098a5-86e1-45ed-9e27-d59dbc7cfc42
md"## 1.4 Test"

# ╔═╡ a389a422-3946-49ba-ac2a-44ba0ebbb7fe
PlutoTest.@test to_negative_positive_tuple(0.0) == (0.0, 0.0, 1)

# ╔═╡ cf832e74-1029-4554-94ff-876e7ce622ab
PlutoTest.@test to_negative_positive_tuple(1.0) == (1.0, 0.0, 0.0)

# ╔═╡ 968529e7-c30c-4ba9-96cc-c3468bbbf7e3
PlutoTest.@test to_negative_positive_tuple(0.5) == (1.0, 1.0, 1.0)

# ╔═╡ 65aab672-9fc2-48f9-b669-e6f5299b1e99
PlutoTest.@test to_negative_positive_tuple(0.75) == (1.0, 0.5, 0.5)

# ╔═╡ 00ee1b88-41fc-41d6-b6a5-91a68560384d
PlutoTest.@test to_negative_positive_tuple(0.25) == (0.5, 0.5, 1)

# ╔═╡ b7e01e25-4b4c-4021-ac52-b5fa4d54c1fa
PlutoTest.@test typeof(to_negative_positive_tuple(0.25f0)) == Tuple{Float32, Float32, Float32}

# ╔═╡ ac0c999b-afb3-4b0c-b981-1bde3275e182
PlutoTest.@test typeof(to_negative_positive_tuple(0.0)) == Tuple{Float64, Float64, Float64}

# ╔═╡ 279cd4f9-9d0e-4bef-80c9-907652371f12
md"## 2 Sensor Simulation

In this part, we want to simulate a sensor. Later, we try to calibrate our artifical sensor.

In short: Each pixel in the sensor measures a number of photons (`n_photons`) which is affected by Poisson noise.
Afterwards, the measured value is altered by additive Gauss noise with a certain
standard deviation `σ`.
This value is then converted to ADU units (analog to digital units) with a certain `gain` (linear conversion factor).
Additionally, each sensor has an `offset` which is finally added.
"

# ╔═╡ ffbc76ee-5c5b-460f-a7dc-ac50500dede9
md"## 2.1 Task
Fill in the missing parts.
"

# ╔═╡ 0b6561db-6286-49b6-bbfb-98f56cbafd1b
"""
      simulate_pixel(n_photons; read_σ=5, offset=0.5, gain=0.1)

In this function we simulate a single pixel.
It transforms a photon number to a digital unit (ADU).

The input is a integer number of photons which is then altered with poisson noise (use `pois_rand` for that).
This number is changed with additive Gauss noise (see HW01).
Finally, the result of this operation is multiplied by the linear `gain`.
At the end, we add `offset`.

* `n_photons`: number of photons
* `read_σ`: additive read noise of the sensor (applied before `gain`)
* `gain`: how much digital output per photon
* `offset`: how much digital offset

## Example Outputs - can vary in your case!
```julia
julia> simulate_pixel(10)
0.8570431305838304

julia> simulate_pixel(10)
2.060753446522642

julia> simulate_pixel(10)
1.020461460534897

julia> simulate_pixel(100)
8.96276245560446

julia> simulate_pixel.(ones((3,3)))
3×3 Matrix{Float64}:
 0.396078  0.96849   -0.0387851
 1.09982   0.210089   0.685605
 0.606973  1.42573   -0.675297
```
"""
function simulate_pixel(n_photons; read_σ=5, offset=0.5, gain=0.1)
	# try to compose the different operations
	# TODO
	# TODO
	return 0 # TODO
end

# ╔═╡ 3b1bbda7-dfa0-4181-b0f7-fc920f166a4b
single_output = simulate_pixel(100)

# ╔═╡ 5df0c449-c46e-4719-bff5-d42220991e8f
md"## 2.1 Test"

# ╔═╡ 02b3bafb-fc03-42ff-a281-47216a04f6ce
# check offset
PlutoTest.@test 0.0123 ≈ simulate_pixel(0, read_σ=0, offset=0.0123, gain=1)

# ╔═╡ 1d3e3c6f-7920-4b2c-a432-105fa1f131f7
# check mean poisson
PlutoTest.@test 123 ≈ mean(simulate_pixel.(ones((1000,)) * 123, read_σ=0, offset=0.0123, gain=1)) rtol=0.1

# ╔═╡ b01359d9-c1b2-455f-897c-8ab6bed4f20f
# check variance poisson
PlutoTest.@test 123 ≈ var(simulate_pixel.(ones((1000,)) * 123, read_σ=0, offset=0.0123, gain=1)) rtol=0.1

# ╔═╡ a598763f-63f9-48c9-9e5b-ba6256775221
# check readnoise
PlutoTest.@test 12.3 ≈ std(simulate_pixel.(zeros((1000,)) * 123, read_σ=12.3, offset=0, gain=1)) rtol=0.1

# ╔═╡ 93567e80-181f-4629-8a27-2f37e5e20b62
# check gain
PlutoTest.@test 42 * 12.3 ≈ std(simulate_pixel.(zeros((1000,)) * 123, read_σ=12.3, offset=0, gain=42)) rtol=0.1

# ╔═╡ 9ad71b79-872d-4c96-b4eb-3c2e841f295d
md"## 2.2 Variance Mean Projection
Now we want to do a Variance Mean Projection to fit the values of `offset` and `gain` since we often don't know them for real sensors.

For this test, you take a sample and in practice you defocus your microscope heavily.
Via that trick, you have an image which covers all values between very dark and very bright.
Hence, we expect that our image sensors both measures a large dynamic range.
"

# ╔═╡ 70fa1c62-b5df-49dd-ae66-b49cb9bf6c0a
# that's already done for you
begin
	ideal = resample(Float32.(testimage("resolution_test_512"))[50:450, 50:450], (150, 150)) # random amount of photons in the interval 0:500
	x_cords = range(-10, 10, length=size(ideal, 1))
	gauss = Float32.(exp.(-(x_cords.^2 .+ x_cords'.^2)))
	gauss ./= sum(gauss)
	img_blurry = conv_psf(ideal, gauss)
	img_blurry .-= minimum(img_blurry)
	img_blurry ./= maximum(img_blurry)
	img_blurry = round.(Ref(Int), img_blurry .* 1000) # round to integers
	Gray.(img_blurry ./ 500)
end;

# ╔═╡ 1e3f86d7-a952-41fc-9f67-c8ce87d14315
md"###### That's an ideal image, where we have both very bright and very dark regions"

# ╔═╡ 120846a5-7892-4f79-8e95-bde9bf06df25
# our blurry sample which emits 1000 photons in maximum
gray_show(img_blurry)

# ╔═╡ 3a6757db-050c-4923-a6bd-5c7c33581a93
histogram(img_blurry[:], xlabel="intensity", ylabel="occurence", title="Histogram Plot", legend=nothing)

# ╔═╡ ee7b9e26-c81d-4e09-9553-e9cded4db31d
md"## 2.2 Task - Apply to image
Now apply the function `simulate_pixel` to `img_blurry`.
Use the automatic broadcasting of Julia for that.
"

# ╔═╡ 6fddf2a4-45f0-4243-aa0c-df5ce2c354a9
simulated_img = similar(img) # TODO, fix that line such that it uses `simulate_pixel`. Use broadcasting mechanism

# ╔═╡ 48fb4f01-a3fa-4c7b-abb0-6e6f2b63dac7
md"Apply the same trick (broadcasting) for a 2D image to this 3D image"

# ╔═╡ f9b327ef-e184-4e1f-af3c-184f408f04bd
md"## 2.3 Task - Calculate mean and variance
Calculate the mean and variance of those images, but only along the third dimension!
Either use a library (Statistics) or write it yourself :)
"

# ╔═╡ ec0289d5-427c-4e19-b63a-361ed7165d42
md"#### Finally, we produce a scatter plot"

# ╔═╡ 1df9763e-81be-4f24-b439-e01559d78a9b
md" ## 2.4 Fitting
Having our noisy dataset, we want to extract the `gain` of the detector since the gain translate from a photon count to an eletrical signal (ADU).

Why does plotting the variance $\sigma^2$ over the mean $\mu$ allows to fit the gain?
We know, that in a Poisson distribution the expected mean $\mu$ is always equal to the variance $\sigma^2$.

If we now take 10 captures of the same image, we always expect the same mean $\mu$ for a certain pixel. Of course, those are 10 evaluation of a random measurement process, therefore we observe fluctuations. Hence, for each pixel we can calculate the mean and the variance.

Since the detector multiplies the measured photons with a certain gain, we do not measure the photons but the digital value.
Hence the measured mean is

$\hat \mu = \text{gain} \cdot \mu$

where $\mu$ would be the expected photon number.

However, the measured variance is

$$\hat \sigma^2 = \sum_{i} (\hat \mu - \text{gain} \cdot x_i)^2 = \text{gain}^2 \sigma^2$$
where $\sigma^2$ would be the variance of the photon number.

The slope of our variance over mean is then

$$\frac{\Delta \hat \sigma^2}{\Delta \hat \mu} = \text{gain}$$

Dividing $\hat \mu$ by $\text{gain}$ returns the measured photon number.

"

# ╔═╡ d6ea97c0-42a2-4119-94cb-c59d07c84992
md" ## 2.5 LsqFit.jl
We now want to fit a linear function with a slope and a offset on our data.
Try to find out [here](https://github.com/JuliaNLSolvers/LsqFit.jl) how that works
"

# ╔═╡ 529cf3ab-090d-4ed3-b236-eed577554bbb
@. model(x, p) = nothing # TODO

# ╔═╡ ca78fb0c-3cd6-461e-a4c7-5bbd61589001
p0 = nothing # TODO (maybe you need to convert it to Float32, maybe not)

# ╔═╡ 709c3594-d979-4a09-876f-7bbf4ad4aa18
fit = nothing # TODO

# ╔═╡ b2b821b3-e95b-48ec-aee0-acea4e6b01dc
md"##  Final Check
If you slide those parameters, you change the simulated values.
"

# ╔═╡ e68aa618-f5b0-4834-b492-268962b4c5be
md"
`offset`=$(@bind offset Slider(0.5:1:100, show_value=true))

`read_σ`=$(@bind read_σ Slider(2:1:20, show_value=true))

`gain`=$(@bind gain Slider(0.01:0.01:1, show_value=true))
"

# ╔═╡ 1140ab5b-6e3a-4fac-baf0-a150a07d7293
begin
	imgs = repeat(img_blurry, outer=(1, 1, 10));
	images = simulate_pixel.(imgs; offset, read_σ, gain);
end;

# ╔═╡ 1e97f466-796e-4837-b0bf-2a3c4edd89c0
md" In practice we would need to take a few images (like 10)

Hence we repeat the image 10 times digitally

The output is an array with `size(imgs) ==` $(size(imgs))
"

# ╔═╡ 3a7737ac-7884-41ba-863c-b352dceb73c9
begin
	# mean
	μ_arr = randn(size(imgs)[1:2]) # calculate the mean of `images` along dim 3
	μ = μ_arr[:]
end

# ╔═╡ 501a5a71-9cff-4338-b575-0f0be06d064f
begin
	# variance
	variance_arr = randn(size(imgs)[1:2])# calculate the variance of `images` along dim 3
	variance = variance_arr[:]
end

# ╔═╡ 41bdbfb9-06b9-4f4b-bc97-31ed2e6ae695
scatter(μ, variance, xlabel="mean", ylabel="variance")

# ╔═╡ 3ffb57c7-825b-4213-a847-7cc4e9466875
begin
	scatter(μ, variance, label="data")
	plot!(μ, model(μ, fit.param), linewidth=5, label="fit")
end

# ╔═╡ 0adf654e-f238-432c-ace4-50ad5d896870
md"The final solution is gain = $(round(fit.param[1], sigdigits=4)) ± $(round(stderror(fit)[1], sigdigits=2))"

# ╔═╡ Cell order:
# ╠═1f6d20c6-2a06-11ec-30b2-f5462ac0ec76
# ╠═46cd6d61-7721-49be-909f-7c411a5562ff
# ╠═85a15a0e-37d4-4b81-986c-d8a2fe4185c7
# ╠═a0b4af53-ded5-4d30-9d2b-7c5d1c2bff80
# ╠═375801bb-7b43-4983-a0ca-211d106ab11c
# ╟─e2768bdd-3346-4f13-88a8-ab8a43094761
# ╠═cea22309-272b-4696-850a-f49ad3120345
# ╠═c66329a6-84f4-4e0e-87a1-8a62f97ba9f7
# ╠═450d3dea-cb97-482b-aaab-bb2919e4ee9e
# ╠═41f64a18-bd0d-49a4-8347-36818b45a576
# ╠═f6379b6e-7030-4967-afea-13494ccad294
# ╟─2b3c2ab1-9bec-43f5-a008-795ac354e080
# ╠═df017bc1-a29c-4bbc-80c3-bc2466cd95f5
# ╠═f05ca6de-b46c-4400-a141-c662a1a16d7c
# ╠═a7007d37-6f30-462b-8be3-70999c980c60
# ╠═f1d7bf0a-6806-45fb-8213-313416017c18
# ╠═18bb1f0f-5b80-43d5-bc8a-526947a70f2c
# ╠═76e2f439-1dff-4f2a-9c05-14e10a13fda6
# ╠═5315fbfd-8ade-46d5-9851-f3ee98c349f7
# ╟─9d4337c8-6cd1-4112-9477-4a40d13fe7e3
# ╠═c0debce4-d81b-4d5d-8907-53d567cf664a
# ╟─4456a207-7f56-4900-95dc-77fe7de31c5d
# ╠═a6fbfa34-33c6-4f09-9dcd-773d92a4db0d
# ╠═b70be939-7f1d-4e80-9d7f-cc2947d57280
# ╠═e9cc6238-5ad9-4a50-93a7-f7065648574a
# ╠═633d295e-cc1f-4864-a36b-11905bc47a50
# ╠═6b65dc02-9a84-43c2-9bb5-32104d23f580
# ╠═160051e5-56b9-413a-a78d-b371a85ff18c
# ╠═1483db2c-6374-4c93-ade4-db8894071ce7
# ╠═0afb1220-1fdf-447d-8e79-b1ad7c274dd3
# ╠═f6e5d3c8-172e-4fb6-89ce-a1827a6ddd6c
# ╠═bd71f8c1-e1e8-4f72-9ad8-8714e9e215f0
# ╠═45d098a5-86e1-45ed-9e27-d59dbc7cfc42
# ╠═a389a422-3946-49ba-ac2a-44ba0ebbb7fe
# ╠═cf832e74-1029-4554-94ff-876e7ce622ab
# ╠═968529e7-c30c-4ba9-96cc-c3468bbbf7e3
# ╠═65aab672-9fc2-48f9-b669-e6f5299b1e99
# ╠═00ee1b88-41fc-41d6-b6a5-91a68560384d
# ╠═b7e01e25-4b4c-4021-ac52-b5fa4d54c1fa
# ╠═ac0c999b-afb3-4b0c-b981-1bde3275e182
# ╠═279cd4f9-9d0e-4bef-80c9-907652371f12
# ╠═d9d77790-30d1-4ca2-83ca-d974b10a9431
# ╠═ffbc76ee-5c5b-460f-a7dc-ac50500dede9
# ╠═0b6561db-6286-49b6-bbfb-98f56cbafd1b
# ╠═3b1bbda7-dfa0-4181-b0f7-fc920f166a4b
# ╠═5df0c449-c46e-4719-bff5-d42220991e8f
# ╠═02b3bafb-fc03-42ff-a281-47216a04f6ce
# ╠═1d3e3c6f-7920-4b2c-a432-105fa1f131f7
# ╠═b01359d9-c1b2-455f-897c-8ab6bed4f20f
# ╠═a598763f-63f9-48c9-9e5b-ba6256775221
# ╠═93567e80-181f-4629-8a27-2f37e5e20b62
# ╠═9ad71b79-872d-4c96-b4eb-3c2e841f295d
# ╠═7a0e442c-9f32-4210-9483-63e0dc9fbe49
# ╠═70fa1c62-b5df-49dd-ae66-b49cb9bf6c0a
# ╠═1e3f86d7-a952-41fc-9f67-c8ce87d14315
# ╠═120846a5-7892-4f79-8e95-bde9bf06df25
# ╠═3a6757db-050c-4923-a6bd-5c7c33581a93
# ╠═ee7b9e26-c81d-4e09-9553-e9cded4db31d
# ╠═6fddf2a4-45f0-4243-aa0c-df5ce2c354a9
# ╠═1e97f466-796e-4837-b0bf-2a3c4edd89c0
# ╠═48fb4f01-a3fa-4c7b-abb0-6e6f2b63dac7
# ╠═1140ab5b-6e3a-4fac-baf0-a150a07d7293
# ╠═f9b327ef-e184-4e1f-af3c-184f408f04bd
# ╠═3a7737ac-7884-41ba-863c-b352dceb73c9
# ╠═501a5a71-9cff-4338-b575-0f0be06d064f
# ╠═ec0289d5-427c-4e19-b63a-361ed7165d42
# ╠═41bdbfb9-06b9-4f4b-bc97-31ed2e6ae695
# ╠═1df9763e-81be-4f24-b439-e01559d78a9b
# ╠═dd4f1881-4570-4972-bd40-a4123ad4b6b1
# ╠═d6ea97c0-42a2-4119-94cb-c59d07c84992
# ╠═529cf3ab-090d-4ed3-b236-eed577554bbb
# ╠═ca78fb0c-3cd6-461e-a4c7-5bbd61589001
# ╠═709c3594-d979-4a09-876f-7bbf4ad4aa18
# ╠═3ffb57c7-825b-4213-a847-7cc4e9466875
# ╠═b2b821b3-e95b-48ec-aee0-acea4e6b01dc
# ╠═e68aa618-f5b0-4834-b492-268962b4c5be
# ╠═0adf654e-f238-432c-ace4-50ad5d896870


text:
module js

import GR

@static if VERSION < v"0.7.0-DEV.4762"
    macro cfunction(f, rt, tup)
        :(Base.cfunction($(esc(f)), $(esc(rt)), Tuple{$(esc(tup))...}))
    end
end

id_count = 0
js_running = false

mutable struct JSTermWidget
    identifier::String
    visible::Bool
end

function inject_js()
  _js_fallback = "https://gr-framework.org/downloads/gr-latest.js"
  _gr_js = if isfile(joinpath(ENV["GRDIR"], "lib", "gr.js"))
    _gr_js = try
      _gr_js = open(joinpath(ENV["GRDIR"], "lib", "gr.js")) do f
        _gr_js = read(f, String)
        _gr_js
      end
    catch e
      nothing
    end
    _gr_js
  end
  _jsterm = """/* jshint esversion: 6 */
  /*eslint no-console: ["error", { allow: ["warn"] }] */
  /*eslint no-unused-vars: ["error", { "varsIgnorePattern": "^_" }]*/
    if (typeof grJSTermRunning === 'undefined') {
      function JSTerm() {
        this.onready = [];
        this.gr = [];
        this.args = [];
        this.canvas = [];
        this.comm = undefined;

        function saveLoad(url, callback, maxtime) {
            let script = document.createElement('script');
            script.onload = function() {
                callback();
            };
            script.onerror = function() {
                console.error(url + ' can not be loaded.');
            };
            script.src = url;
            document.head.appendChild(script);
            setTimeout(function() {
                if (!grJstermReady) {
                    console.error(url + ' can not be loaded.');
                }
            }, maxtime);
        }

        function jsLoaded() {
            grJstermReady = true;
            for (let i = 0; i < this.onready.length; i++) {
                this.onready[i]();
            }
            this.onready = [];
        }

        function canvas_removed(id) {
          this.comm.send({"type":"removed","content":id})
        }

        function draw(msg) {
            if (!grJstermReady) {
                this.onready.push(function() {
                    return draw(msg);
                });
            } else if (!GR.is_ready) {
                GR.ready(function() {
                    return draw(msg);
                });
            } else {
                if (typeof this.canvas['jsterm-' + msg.content.data.canvasid] === 'undefined') {
                  this.canvas['jsterm-' + msg.content.data.canvasid] = document.getElementById('jsterm-' + msg.content.data.canvasid);
                  this.canvas['jsterm-' + msg.content.data.canvasid].addEventListener('DOMNodeRemoved', function() {
                      canvas_removed(msg.content.data.canvasid);
                  });
                }
                if (typeof this.gr['jsterm-' + msg.content.data.canvasid] === 'undefined') {
                    this.gr['jsterm-' + msg.content.data.canvasid] = new GR('jsterm-' + msg.content.data.canvasid);
                }
                if (!this.args.hasOwnProperty('jsterm-' + msg.content.data.canvasid) || typeof this.args['jsterm-' + msg.content.data.canvasid] === 'undefined') {
                    this.args['jsterm-' + msg.content.data.canvasid] = this.gr['jsterm-' + msg.content.data.canvasid].newmeta();
                }
                if (document.getElementById('jsterm-' + msg.content.data.canvasid) == null) {
                    canvas_removed(msg.content.data.canvasid);
                    this.gr['jsterm-' + msg.content.data.canvasid] = undefined;
                    setTimeout(function(){refr_plot('jsterm-' + msg.content.data.canvasid, msg, 0);}, 5);
                } else {
                    this.gr['jsterm-' + msg.content.data.canvasid].readmeta(this.args['jsterm-' + msg.content.data.canvasid], msg.content.data.json);
                    this.gr['jsterm-' + msg.content.data.canvasid].select_canvas();
                    this.gr['jsterm-' + msg.content.data.canvasid].clearws();
                    this.gr['jsterm-' + msg.content.data.canvasid].plotmeta(this.args['jsterm-' + msg.content.data.canvasid]);
                }
            }
        }

        function refr_plot(id, msg, count) {
            if (count >= 500) {
              return;
            }
            if (document.getElementById(id) == null) {
                setTimeout(function(){refr_plot(id, msg, count + 1);}, 5);
            } else {
                draw(msg);
            }
        }

        function register_comm(kernel) {
            kernel.comm_manager.register_target('jsterm_comm', function(comm) {
                comm.on_msg(function(msg) {
                    draw(msg);
                });
                comm.on_close(function() {
                });
                window.addEventListener('beforeunload', function (e) {
                  comm.close();
                });
                this.comm = comm;
            });
        }

        function onLoad() {
            Jupyter.notebook.events.on('execution_request.Kernel', function() {
                for (var key in this.gr) {
                    if (this.args.hasOwnProperty(key)) {
                        this.gr[key].deletemeta(this.args[key]);
                    }
                }
                this.gr = [];
                this.args = [];
            });
            // TODO restart event
            let kernel = Jupyter.notebook.kernel;
            if (typeof kernel === 'undefined' || kernel == null) {
                console.error('JSTerm: No kernel detected');
                return;
            }
            register_comm(kernel);
            Jupyter.notebook.events.on('kernel_ready.Kernel', function() {
                kernel = IPython.notebook.kernel;
                register_comm(kernel);
            });
        }
        onLoad();
    }
    JSTerm();
  }
  var grJSTermRunning = true;"""
  if _gr_js === nothing
      _gr_js = string("""
        saveLoad('""", _js_fallback, """', jsLoaded, 10000);
        var grJstermReady = false;
      """)
  else
    _gr_js = string(_gr_js, "var grJstermReady = true;")
  end
  display(HTML(string("""
    <script type="text/javascript">
      """, _gr_js, """
      """, _jsterm, """
    </script>
  """)))
end

function JSTermWidget(name::String, id::Int64)
  global id_count, js_running
  if GR.isijulia()
    id_count += 1
    if !js_running
      inject_js()
      js_running = true
    end
    JSTermWidget(string(name, id), false)
  else
    error("JSTermWidget is only available in IJulia environments")
  end
end

function jsterm_display(widget::JSTermWidget)
  if GR.isijulia()
    display(HTML(string("<canvas id=\"jsterm-", widget.identifier, "\" width=\"500\" height=\"500\"></canvas>")))
    widget.visible = true
  else
    error("jsterm_display is only available in IJulia environments")
  end
end

function jsterm_send(widget::JSTermWidget, data::String)
  global js_running
  if GR.isijulia()
    comm = Main.IJulia.Comm("jsterm_comm")
    comm.on_close = function comm_close_callback(msg)
      global js_running
      js_running = false
    end
    comm.on_msg = function comm_msg_callback(msg)
      if haskey(msg.content["data"], "type")
        if msg.content["data"]["type"] == "removed"
          jswidgets[msg.content["data"]["content"]].visible = false
          jsterm_display(widget)
        end
      end
    end
    Main.IJulia.send_comm(comm, Dict("json" => data, "canvasid" => widget.identifier))
  else
    error("jsterm_send is only available in IJulia environments")
  end
end

function recv(name::Cstring, id::Int64, msg::Cstring)
    # receives string from C and sends it to JS via Comm
    global js_running
    if !js_running
      inject_js()
    end
    _name = unsafe_string(name)
    if haskey(jswidgets, string(_name, id))
        widget = jswidgets[string(_name, id)]
    else
        widget = JSTermWidget(_name, id)
        jswidgets[string(_name, id)] = widget
    end
    if !widget.visible
      jsterm_display(widget)
    end
    jsterm_send(widget, unsafe_string(msg))
    return convert(Int32, 1)
end

function send(name::Cstring, id::Int64)
    # Dummy function, not in use
    print(string("send", name, id, "\n"))
    return convert(Cstring, "String")
end

jswidgets = nothing
send_c = nothing
recv_c = nothing

function initjs()
    global jswidgets, send_c, recv_c

    jswidgets = Dict{String, JSTermWidget}()
    send_c = @cfunction(send, Cstring, (Cstring, Int64))
    recv_c = @cfunction(recv, Int32, (Cstring, Int64, Cstring))
    send_c, recv_c
end

end # module


text:
using DReal
import DReal: opensmt_context, opensmt_set_verbosity, opensmt_assert, opensmt_check, opensmt_print_expr, opensmt_mk_context,
              opensmt_mk_bool_var, opensmt_mk_int_var, opensmt_mk_minus, opensmt_mk_leq, opensmt_mk_num_from_string,
              opensmt_push, opensmt_pop, opensmt_reset, opensmt_del_context
using Base.Test
using Compat

begin
# Blank slate
reset_ctx!()

print( "Creating context\n" )
ctx = opensmt_mk_context(@compat UInt32(1))

 # Creating boolean variables
print( "Creating some boolean variables\n" )
a = opensmt_mk_bool_var( ctx, "a" )
b = opensmt_mk_bool_var( ctx, "b" )
c = opensmt_mk_bool_var( ctx, "c" )

 # Creating integer variables
print( "Creating some integer variables\n" )
x = opensmt_mk_int_var( ctx, "x" , -100, 100)
y = opensmt_mk_int_var( ctx, "y" , -100, 100)
z = opensmt_mk_int_var( ctx, "z" , -100, 100)
 # Creating inequality
print( "Creating x - y <= 0\n" )

print( "  Creating x - y\n" )
minus = opensmt_mk_minus( ctx, x, y )
print( "  Expression created: " )
opensmt_print_expr( minus )
print( "\n" )

print( "  Creating 0\n" )
zero = opensmt_mk_num_from_string( ctx, "0" )
print( "  Expression created: " )
opensmt_print_expr( zero )
print( "\n" )

print( "  Creating x - y <= 0\n" )
leq = opensmt_mk_leq( ctx, minus, zero )
print( "  Expression created: " )
opensmt_print_expr( leq )
print( "\n" )

 # Creating inequality 2
print( "Creating y - z <= 0\n" )
minus2 = opensmt_mk_minus( ctx, y, z )
leq2 = opensmt_mk_leq( ctx, minus2, zero )
print( "  Expression created: " )
opensmt_print_expr( leq2 )
print( "\n" )

 # Creating inequality 3
print( "Creating z - x <= -1\n" )
minus3 = opensmt_mk_minus( ctx, z, x )
mone = opensmt_mk_num_from_string( ctx, "-1" )
leq3 = opensmt_mk_leq( ctx, minus3, mone )
print( "  Expression created: " )
opensmt_print_expr( leq3 )
print( "\n" )

 # Creating inequality 4
print( "Creating z - x <= 0\n" )
minus4 = opensmt_mk_minus( ctx, z, x )
leq4 = opensmt_mk_leq( ctx, minus4, zero )
print( "  Expression created: " )
opensmt_print_expr( leq4 )
print( "\n" )

 # Asserting first inequality
print( "Asserting ")
opensmt_print_expr( leq )
print( "\n" )
opensmt_assert( ctx, leq )
opensmt_push( ctx )

println("Expression is")
opensmt_print_expr( leq )

# Checking for consistency
println( "\nChecking for consistency: " )
res = opensmt_check( ctx )
println( "res = $(res == -1 ? "unsat" : "sat")")
@test res == 1

# Asserting second inequality
print( "Asserting ")
opensmt_print_expr( leq2 )
print( "\n" )
opensmt_assert( ctx, leq2 )
opensmt_push( ctx )

 # Checking for consistency
print( "\nChecking for consistency: " )
res = opensmt_check( ctx )
println( "res = $(res == -1 ? "unsat" : "sat")")
@test res == 1

# Asserting third inequality
print( "\nAsserting ")
opensmt_print_expr( leq3 )
print( "\n" )
opensmt_push( ctx )
opensmt_assert( ctx, leq3 )

# Checking for consistency - it is UNSAT
print( "\nChecking for consistency: " )
res = opensmt_check( ctx )
@test res == -1

 # Backtracking one step - so the state is SAT again
opensmt_pop( ctx )
print( "Backtracking one step\n\n" )

# Asserting fourth inequality
print( "Asserting ")
opensmt_print_expr( leq4 )
print( "\n" )
opensmt_assert( ctx, leq4 )

 # Checking for consistency
print( "\nChecking for consistency: " )
res = opensmt_check( ctx )
println( "res = $(res == -1 ? "unsat" : "sat")")
@test res == 1

# Resetting context
print( "Resetting context\n" )
opensmt_reset( ctx )

# # Deleting context
# print( "Deleting context\n" )
# opensmt_del_context( ctx )
end


text:
"""
    Object

Overarching type for objects displayable on `Axes`.
"""
abstract type Object end

"""
    Object

Overarching type for objects displayable on `Axes2D`.
"""
abstract type Object2D <: Object end

"""
    Object

Overarching type for objects displayable on `Axes3D`.
"""
abstract type Object3D <: Object end

"""
    Text <: Object

Place text somewhere relative to current axes.
"""
@with_kw mutable struct Text2D <: Object2D
    text::String            # what to write
    anchor::T2F             # where the text is located (x,y) position, relative to the axes
    position::String = "cc" # position with respect to anchor (by default centered on anchor)
    textstyle::TextStyle = TextStyle()
end

"""
    StraightLine2D <: Object2D

Place either a vertical or horizontal straightline at given anchor (from axis to axis).
"""
@with_kw mutable struct StraightLine2D <: Object2D
    anchor::Float64
    horiz::Bool
    linestyle::LineStyle = LineStyle()
end

"""
    Box2D <: Object2D

Place a 2D filled box.
"""
@with_kw mutable struct Box2D <: Object2D
    anchor::T2F  # where the box is (in graph units)
    size::T2F    # width, hei (in graph units)
    #
    position::String     = "bl"        # position of the anchor with respect to the box
    nobox::Bool          = true        # show an edge or not
    linestyle::LineStyle = LineStyle() # style of the box edge
    fillstyle::Option{FillStyle} = ∅   # if none --> transparent box
    # text::String
    # textStyle::TextStyle = TextStyle()
end

"""
    Colorbar <: Object2D

Add a colorbar.
"""
@with_kw mutable struct Colorbar <: Object2D
    zmin::Float64
    zmax::Float64
    cmap::Vector{Color}
    # -- optional things
    ticks::Ticks # constructed
    # --
    size::Option{T2F} = ∅ # (width, height)
    # --
    pixels::Int      = 100        # resolution for the color bar
    nobox::Bool      = true       #
    position::String = "right"    # left, right, bottom, top
    offset::T2F      = (0.3, 0.0) # h, v distance from closest corner of the graph
end
Colorbar(zmin::Float64, zmax::Float64, cmap::Vector{<:Color}) =
    Colorbar(zmin=zmin, zmax=zmax, cmap=cmap,
             ticks=Ticks(places=collect(range(zmin, stop=zmax, length=5))[2:end-1]))


text:
export TwoDPoisson

struct TwoDPoisson
    nx::Int
    ny::Int
    degree::Int
    dx::Float64
    dy::Float64

    eig_values_dtm1d_1::Vector{Float64}
    eig_values_dtm1d_2::Vector{Float64}
    eig_values_d1::Vector{ComplexF64}
    eig_values_d2::Vector{ComplexF64}
    eig_values_mass_0_1::Vector{Float64}
    eig_values_mass_0_2::Vector{Float64}

    array1d_x::Vector{ComplexF64}
    array1d_y::Vector{ComplexF64}
    scratch::Array{ComplexF64,2}
    scratchx::Array{ComplexF64,2}
    scratchy::Array{ComplexF64,2}

    function TwoDPoisson(mesh, degree)
        nx, ny = mesh.nx, mesh.ny
        dx, dy = mesh.dx, mesh.dy

        mass_line_0 = zeros(degree + 1)
        mass_line_1 = zeros(degree)
        eig_values_mass_1_1 = zeros(nx)
        eig_values_mass_1_2 = zeros(ny)

        array1d_x = zeros(ComplexF64, nx)
        array1d_y = zeros(ComplexF64, ny)
        scratch = zeros(ComplexF64, nx, ny)
        scratchx = zeros(ComplexF64, nx, ny)
        scratchy = zeros(ComplexF64, nx, ny)

        # Eigenvalues of mass matrices
        mass_line_0 = spline_fem_mass_line(degree)
        mass_line_1 = spline_fem_mass_line(degree - 1)

        eig_values_mass_0_1 = spline_fem_compute_mass_eig(nx, degree, mass_line_0 * dx)
        eig_values_mass_1_1 = spline_fem_compute_mass_eig(nx, degree - 1, mass_line_1 * dx)

        eig_values_mass_0_2 = spline_fem_compute_mass_eig(ny, degree, mass_line_0 * dy)
        eig_values_mass_1_2 = spline_fem_compute_mass_eig(ny, degree - 1, mass_line_1 * dy)

        eig_values_d1 = zeros(ComplexF64, nx)
        eig_values_dtm1d_1 = zeros(nx)

        for j in 2:nx
            angle = 2π * (j - 1) / nx
            eig_values_d1[j] = (1 - cos(angle)) / dx + 1im * sin(angle) / dx
            eig_values_dtm1d_1[j] = 2 / dx^2 * (1 - cos(angle)) * eig_values_mass_1_1[j]
        end

        eig_values_d2 = zeros(ComplexF64, ny)
        eig_values_dtm1d_2 = zeros(ny)

        for j in 2:ny
            angle = 2π * (j - 1) / ny
            eig_values_d2[j] = (1 - cos(angle)) / dy + 1im * sin(angle) / dy
            eig_values_dtm1d_2[j] = 2 / dy^2 * (1 - cos(angle)) * eig_values_mass_1_2[j]
        end

        return new(
            nx,
            ny,
            degree,
            dx,
            dy,
            eig_values_dtm1d_1,
            eig_values_dtm1d_2,
            eig_values_d1,
            eig_values_d2,
            eig_values_mass_0_1,
            eig_values_mass_0_2,
            array1d_x,
            array1d_y,
            scratch,
            scratchx,
            scratchy,
        )
    end
end

function spline_fem_mass_line(degree::Int)
    n = degree + 1
    spline_val = zeros(degree + 1, n)

    x, w = gausslegendre(n)
    x .+= 1
    x .*= 0.5
    w .*= 0.5

    for j in 1:n
        spline_val[:, j] .= uniform_bsplines_eval_basis(degree, x[j])
    end

    mass_line = zeros(degree + 1)
    for j in 1:(degree + 1), i in j:(degree + 1), k in 1:n
        mass_line[j] += spline_val[i, k] * spline_val[i - j + 1, k] * w[k]
    end

    return mass_line
end

"""
    spline_fem_compute_mass_eig( n_cells, degree, mass_line, eig_values_mass )

Compute eigenvalues of mass matrix with given line @ mass_line
"""
function spline_fem_compute_mass_eig(n_cells, degree, mass_line)
    eig_values_mass = zeros(n_cells)

    factor = 2π / n_cells

    for k in 1:n_cells
        eig_values_mass[k] = mass_line[1]
        for j in 1:degree
            eig_values_mass[k] += mass_line[j + 1] * 2 * cos(factor * (k - 1) * j)
        end
    end

    return eig_values_mass
end

export compute_rhs_from_function

"""
    compute_rhs_from_function(solver, func, coefs_dofs)

Compute the FEM right-hand-side for a given function f and periodic splines of given degree
Its components are ``\\int f N_i dx`` where ``N_i`` is the B-spline starting at ``x_i``
- solver :: Maxwell solver
- func : function
- coefs_dofs[:]  : Finite Element right-hand-side
"""
function compute_rhs_from_function(solver::TwoDPoisson, f)
    nx, ny = solver.nx, solver.ny
    dx, dy = solver.dx, solver.dy
    d = solver.degree

    rhs = zeros(nx * ny)

    # take enough Gauss points so that projection is exact for splines of degree deg
    # rescale on [0,1] for compatibility with B-splines
    bspl_d1 = zeros(d + 1, d + 1)
    xg, wg = gausslegendre(d + 1)
    xg .+= 1
    xg ./= 2
    wg ./= 2
    # Compute bsplines at gauss_points
    for k in 1:(d + 1)
        bspl_d1[k, :] .= uniform_bsplines_eval_basis(d, xg[k])
    end

    counter = 1
    # Compute coefs_dofs = int f(x)N_i(x)
    for i2 in 1:ny, i1 in 1:nx
        coef = 0.0
        # loop over support of B spline
        for j1 in 1:(d + 1), j2 in 1:(d + 1)
            # loop over Gauss points
            for k1 in 1:(d + 1), k2 in 1:(d + 1)
                x = dx * (xg[k1] + (i1 + j1 - 2))
                y = dy * (xg[k2] + (i2 + j2 - 2))
                coef +=
                    wg[k1] * wg[k2] * f(x, y) *
                    bspl_d1[k1, d + 2 - j1] *
                    bspl_d1[k2, d + 2 - j2]
            end
        end
        # rescale by cell size
        rhs[counter] = coef * dx * dy
        counter += 1
    end

    return rhs
end

function ifft2d!(outval, solver, inval)
    nx, ny = solver.nx, solver.ny

    for i in 1:nx
        for j in 1:ny
            solver.array1d_y[j] = inval[i, j]
        end
        ifft!(solver.array1d_y)
        for j in 1:ny
            inval[i, j] = solver.array1d_y[j]
        end
    end

    k = 0

    for j in 1:ny
        for i in 1:nx
            solver.array1d_x[i] = inval[i, j]
        end
        ifft!(solver.array1d_x)
        for i in 1:nx
            k = k + 1
            outval[k] = real(solver.array1d_x[i])
        end
    end
end

function fft2d!(solver, rho)
    nx, ny = solver.nx, solver.ny

    k = 0
    for j in 1:ny
        for i in 1:nx
            k = k + 1
            solver.array1d_x[i] = rho[k]
        end
        fft!(solver.array1d_x)
        for i in 1:nx
            solver.scratch[i, j] = solver.array1d_x[i]
        end
    end

    for i in 1:nx
        for j in 1:ny
            solver.array1d_y[j] = solver.scratch[i, j]
        end
        fft!(solver.array1d_y)
        for j in 1:ny
            solver.scratch[i, j] = solver.array1d_y[j]
        end
    end
end

function compute_e_from_rho!(efield, solver::TwoDPoisson, rho)
    nx, ny = solver.nx, solver.ny
    # Compute Fourier transform
    fft2d!(solver, rho)

    # Apply inverse matrix of eigenvalues on mode
    for j in 1:ny, i in 1:nx
        if (i == 1 && j == 1)
            solver.scratch[i, j] = 0
        else
            eig_val = (
                solver.eig_values_dtm1d_1[i] * solver.eig_values_mass_0_2[j] +
                solver.eig_values_mass_0_1[i] * solver.eig_values_dtm1d_2[j]
            )
            solver.scratch[i, j] = solver.scratch[i, j] / eig_val
        end

        solver.scratchx[i, j] = -solver.scratch[i, j] * solver.eig_values_d1[i]
        solver.scratchy[i, j] = -solver.scratch[i, j] * solver.eig_values_d2[j]
    end

    # Compute inverse Fourier transfrom

    ifft2d!(efield[1], solver, solver.scratchx)
    return ifft2d!(efield[2], solver, solver.scratchy)
end


text:
using ControlBenchmarks
using ControlSystems

# Test the default constructor
let benchmarkProb = controlbenchmark( HydraulicPositioningSystem() )
    @test ControlSystems.nstates( benchmarkProb.sys ) == 3
    @test ControlSystems.ninputs( benchmarkProb.sys ) == 1
end

# Test setting a non-default value for a parameter
let hydPosSys = HydraulicPositioningSystem()
    hydPosSys.B = 10000

    benchmarkProb = controlbenchmark( hydPosSys )
    @test ControlSystems.nstates( benchmarkProb.sys ) == 3
    @test ControlSystems.ninputs( benchmarkProb.sys ) == 1
end

# Test setting a non-default value for a parameter that is out of bounds
let hydPosSys = HydraulicPositioningSystem()
    hydPosSys.B = 8000

    # Test without ignoring the bounds
    @test_throws DomainError controlbenchmark( hydPosSys; ignoreBounds = false )

    # Test while ignoring the bounds
    benchmarkProb = controlbenchmark( hydPosSys; ignoreBounds = true )
    @test ControlSystems.nstates( benchmarkProb.sys ) == 3
    @test ControlSystems.ninputs( benchmarkProb.sys ) == 1
end


text:
"""
    geno_freq(x::Array{Union{Missing, Tuple},1})
Calculate genotype frequencies of all loci in a `PopObj`. Returns a `Dict` of
genotypes and their frequencies.
"""
function geno_freq(x::Array{Union{Missing,Tuple},1})
    d = Dict()
    # conditional testing if all genos are missing
    ismissing.(x) |> unique == [true] && return missing
    for row in x
        # sum up missing
        if row === missing
            continue
        else
        # sum up non-missing genotypes
            d[row] = get!(d, row, 0) + 1
        end
    end
    total = values(d) |> sum    # sum of all non-missing genotypes
    [d[i] = d[i] / total for i in keys(d)] # genotype count/total
    return d
end

"""
    geno_freq(x::SubArray{Union{Missing, Tuple},1})
Calculate genotype frequencies of all loci in `PopObj` split
by population using `group()`. Returns a `Dict` of genotypes and their frequencies.
"""
function geno_freq(x::SubArray{Union{Missing,Tuple},1})
    d = Dict()
    # conditional testing if all genos are missing
    ismissing.(x) |> unique == [true] && return missing
    for row in x
        # sum up missing
        if row === missing
            continue
        else
        # sum up non-missing genotypes
            d[row] = get!(d, row, 0) + 1
        end
    end
    total = values(d) |> sum    # sum of all non-missing genotypes
    [d[i] = d[i] / total for i in keys(d)] # genotype count/total
    return d
end

#=============================================================
REMOVE?

"""
  allele_freq_alpha(x::PopObj)
Returns an array of `Dicts` of allele counts per locus
"""
function allele_freq_alpha(x::PopObj)
    y = PopOpt(x)
    tmp = names(y.loci)[1]  # restrict to single locus for testing
    d = Dict()
    for row in y.loci[!, tmp]
        if row === missing
            continue
        else
            for allele in row
                d[allele] = get!(d, allele, 0) + 1
            end
        end
    end
    total = values(d) |> sum
    [d[i] = d[i] / total for i in keys(d)]
    return d
end
==============================================================#

"""
    allele_freq_mini(x::Array{Union{Missing, Tuple},1})
Calculate allele counts for a single locus of a `PopObj`. Returns a `Dict` of
allele's and their frequencies.
"""
function allele_freq_mini(x::Array{Union{Missing,Tuple},1})
    d = Dict()
    for row in x
        # sum up missing
        if row === missing
            continue
        else
        # sum up alleles
            for allele in row
                d[allele] = get!(d, allele, 0) + 1
            end
        end
    end
    total = values(d) |> sum    # sum of all non-missing alleles
    [d[i] = d[i] / total for i in keys(d)]  # allele count / sum
    return d
end

"""
    allele_freq_mini(x::SubArray{Union{Missing, Tuple},1})
Calculate allele counts for a single locus of a `PopObj` split by population
using `group()`. Returns a `Dict` of allele's and their frequencies.
"""
function allele_freq_mini(x::SubArray{Union{Missing,Tuple},1})
    d = Dict()
    for row in x
        # sum up missing
        if row === missing
            continue
        else
        # sum up alleles
            for allele in row
                d[allele] = get!(d, allele, 0) + 1
            end
        end
    end
    total = values(d) |> sum    # sum of all non-missing alleles
    [d[i] = d[i] / total for i in keys(d)]  # allele count / sum
    return d
end


text:
function cartToKep(x::AbstractArray, μ::AbstractFloat)

    # Position and velocity norms
    r = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
    v = sqrt(x[4]^2 + x[5]^2 + x[6]^2)

    # Precompute repeated computations
    rinv  = 1.0 / r
    μinv  = 1.0 / μ
    vsq   = v^2
    rDotV = x[1]*x[4] + x[2]*x[5] + x[3]*x[6]
    μRinv = μ*rinv
    twoπ  = 2.0*π

    # Special case tolerance
    scTol = 1.0e-14;

    # Special case flag
    scFlag = 0

    # Angular momentum
    hv = SVector(x[2]*x[6] - x[3]*x[5],
                 x[3]*x[4] - x[1]*x[6],
                 x[1]*x[5] - x[2]*x[4])
    h  = sqrt(hv[1]^2 + hv[2]^2 + hv[3]^2)

    # Normal vector
    nv = SVector(-hv[2], hv[1], 0.0)
    n = sqrt(nv[1]^2 + nv[2]^2)

    # Eccentricity
    term1 = (vsq - μRinv)*μinv
    term2 = rDotV*μinv
    ev = SVector(term1*x[1] - term2*x[4],
                 term1*x[2] - term2*x[5],
                 term1*x[3] - term2*x[6])
    e = sqrt(ev[1]^2 + ev[2]^2 + ev[3]^2)

    # Specific energy
    enrgy = 0.5*vsq - μRinv

    # Semiparameter
    if e != 1.0
        a = -μ / (2.0*enrgy)
        p = a*(1 - e^2)
    else
        a = Inf
        p = h^2*μinv
    end

    # Inclination
    i = acos(hv[3]/h)

    # Check if elliptical
    if e > scTol

        # RAAN
        Ω = (nv[2] > 0.0) ? acos(nv[1]/n) :
            twoπ - acos(nv[1]/n)

        # True anomoly
        eDotR = ev[1]*x[1] + ev[2]*x[2] + ev[3]*x[3]
        ν = (rDotV > 0.0) ? acos(eDotR / (r*e)) :
                twoπ - acos(eDotR / (r*e))

        # Argument of periapsis/True longitude of periapsis
        if i > scTol && abs(i - π) > scTol
            nDotE = nv[1]*ev[1] + nv[2]*ev[2] + nv[3]*ev[3]
            ω = (ev[3] > 0.0) ? acos(nDotE / (n*e)) :
                    twoπ - acos(nDotE / (n*e))

            return (SVector(a, e, i, Ω, ω, ν), scFlag)

        else # Elliptical equatorial
            scFlat = 1
            ωt = (ev[2] > 0.0) ? acos(e[1] / e) :
                    twoπ - acos(e[1] / e)

            return (SVector(a, e, i, Ω, ωt, ν), scFlag)
        end

    else # Circular
        if i > scTol && abs(i - π) > scTol # inclined

            # RAAN
            Ω = (nv[2] > 0.0) ? acos(nv[1]/n) :
                twoπ - acos(nv[1]/n)

            # Argument of latitude
            nDotR = nv[1]*x[1] + nv[2]*x[2] + nv[3]*x[3]
            u = (r[3] > 0.0) ? acos(nDotR / (n*r)) :
                    twoπ - acos(nDotR / (n*r))

            scFlag = 2
            return (SVector(a, e, i, Ω, NaN, u), scFlag)

        else # equitorial

            # True longitude
            λt = (x[2] > 0.0) ? acos(r[1] * rinv) :
                    twoπ - acos(r[1] * rinv)

            scFlag = 3
            return (SVector(a, e, i, NaN, NaN, λt), scFlag)
        end
    end
end


text:
"""
   abstract type CellDofBasis end

Abstract type that represents a cell array of `Dof`. The main motivation for
its definition is to provide a trait that informs whether the `Dof` entries are
defined for functions in the reference or physical space
"""
abstract type CellDofBasis end

RefStyle(::Type{<:CellDofBasis}) = @abstractmethod
get_array(::CellDofBasis) = @abstractmethod

function test_cell_dof_basis(cf::CellDofBasis,f::CellField)
  ar = get_array(cf)
  @test isa(ar,AbstractArray)
  a = evaluate(cf,f)
  _ = collect(a)
  rs = RefStyle(cf)
  @test isa(get_val_parameter(rs),Bool)
end

"""
    evaluate(dof_array::CellDofBasis,field_array::AbstractArray)

Evaluates the `CellDofBasis` for the `Field` objects
at the array `field` element by element.

The result is numerically equivalent to

    map(evaluate_dof, dof_array.array, field_array)

but it is described with a more memory-friendly lazy type.
"""
function evaluate(cell_dofs::CellDofBasis,cell_field::CellField)
 _evaluate_cell_dofs(cell_dofs,cell_field,RefStyle(cell_dofs))
end

function  _evaluate_cell_dofs(cell_dofs,cell_field,ref_trait::Val{true})
  ReferenceFEs.evaluate_dof_array(get_array(cell_dofs),get_array(to_ref_space(cell_field)))
end

function  _evaluate_cell_dofs(cell_dofs,cell_field,ref_trait::Val{false})
  ReferenceFEs.evaluate_dof_array(get_array(cell_dofs),get_array(to_physical_space(cell_field)))
end

"""
"""
struct GenericCellDofBasis{R} <: CellDofBasis
   ref_trait::Val{R}
   array::AbstractArray{<:Dof}
end

RefStyle(::Type{<:GenericCellDofBasis{R}}) where R = Val{R}()

get_array(a::GenericCellDofBasis) = a.array

function reindex(cf::CellDofBasis,a::AbstractVector)
  similar_object(cf,reindex(get_array(cf),a))
end

function similar_object(cf::CellDofBasis,a::AbstractArray)
  GenericCellDofBasis(RefStyle(cf),a)
end


text:
typealias WordId Int32          # [1:nvocab]
typealias DepRel UInt8          # [1:ndeps]
typealias PosTag UInt8          # [1:npostag]
typealias Position Int32        # [1:nword]
typealias Cost Int32            # [0:nword]
typealias Move Int              # [1:nmove]
typealias SFtype Int32          # representing sparse feature in x
typealias WVtype Float32        # element type for word vectors
typealias Pvec Vector{Position} # used for stack, head
typealias Dvec Vector{DepRel}   # used for deprel
Pzeros(n::Integer...)=zeros(Position, n...)
Droots(n::Integer...)=ones(DepRel, n...) # ROOT=1

type Vocab
    words::Vector{String}
    postags::Vector{String}
    deprels::Vector{String}
    wvecs::Matrix{WVtype}
    fdict::Dict{Any,SFtype}
    Vocab()=new(String[],String[],String[],Array(WVtype,0,0),Dict{Any,SFtype}())
end

type Sentence
    word::Vector{WordId}
    postag::Vector{PosTag}
    head::Vector{Position}
    deprel::Vector{DepRel}
    vocab::Vocab
    # wvec::Array{Float32}
    # TODO: consider using lc(word), lemma, xpostag, feats, deps, misc, lang-specific deprels ext
    # TODO: add word vectors, vocab, how do we add context vectors?
    Sentence()=new(WordId[],PosTag[],Position[],DepRel[])
end

typealias Corpus AbstractVector{Sentence}

wdim(s::Sentence)=size(s.vocab.wvecs,1)
wcnt(s::Sentence)=length(s.word)
wtype(s::Sentence)=eltype(s.vocab.wvecs)
wdim(c::Corpus)=wdim(c[1])
wcnt(c::Corpus)=(n=0;for s in c; n+=wcnt(s); end; n)
wtype(c::Corpus)=wtype(c[1])

function Base.show(io::IO, s::Sentence)
    print(io, "[")
    d = s.vocab.words
    print(io, (wcnt(s) <= 6 ?
               join(d[s.word], " ") :
               join([d[s.word[1:3]]..., "…", d[s.word[end-2:end]]...], " ")))
    print(io, "]")
    print(io, map(Int, s.head))
end


text:
export bin_vector, bin_matrix, bin_value
export unbin_value, unbin_matrix, unbin_vector
export combine_binned_matrix, combine_binned_vector
export combine_and_relabel_binned_matrix
export unary_of_matrix
export relabel

bin_value(v::Float64, bins::Int64, min=-1.0, max=1.0) = round(Int64, minimum([1+floor((v-min) / (max - min) * (bins-0.1)), bins]))

bin_vector(vec::Vector{Float64}, min::Float64, max::Float64, bins::Int64) = map(v->bin_value(v, bins, min, max), vec)

unbin_vector(vec::Vector{Float64}, min::Float64, max::Float64, bins::Int64; mode="centre") = map(v->unbin_value(v, bins, min, max, mode=mode), vec)

function unbin_value(v::Int64, bins::Int64, min=-1.0, max=1.0; mode="centre")
  known_mode = (mode == "centre" ||
                mode == "lower"  ||
                mode == "upper")
  @assert known_mode "Mode may be any of the following: [\"centre\", \"lower\", \"upper\"]"

  delta = (max - min) / Float64(bins)
  u = (v - 1) * delta + min

  if     mode == "centre"
    return u + 0.5 * delta
  elseif mode == "upper"
    return u + delta
  end
  return u
end

function bin_matrix(m::Matrix{Float64}, min::Float64, max::Float64, bins::Int64)
  r = zeros(size(m))
  for i=1:size(m)[1]
    r[i,:] = bin_vector(squeeze(m[i,:],2), min, max, bins)
  end
  convert(Matrix{Int64}, r)
end

function unbin_matrix(m::Matrix{Float64}, min::Float64, max::Float64, bins::Int64; mode="centre")
  r = zeros(size(m))
  for i=1:size(m)[1]
    r[i,:] = unbin_vector(squeeze(m[i,:],1), min, max, bins, mode=mode)
  end
  convert(Matrix{Int64}, r)
end

function combine_binned_vector(v::Vector{Int64}, bins::Int64)
  convert(Int64, sum([v[i] * bins^(i-1) for i = 1:length(v)]))
end

function combine_binned_matrix(v::Matrix{Int64})
  l = size(v)[1]
  r = zeros(l)
  m = maximum(v)
  for i = 1:l
    r[i] = combine_binned_vector(squeeze(v[i,:],1), m)
  end
  convert(Vector{Int64}, r)
end

relabel(v::Vector{Int64}) = indexin(v, sort(unique(v)))

combine_and_relabel_binned_matrix(data::Matrix{Int64}) = relabel(combine_binned_matrix(data))

unary_of_matrix(data::Matrix{Float64}, min::Float64, max::Float64, bins::Int64) = combine_and_relabel_binned_matrix(bin_matrix(data, min, max, bins))


text:
using Scallop.Metrics

@assert 12 == xmetrics(6)


text:
using NumberUnions

if VERSION >= v"0.7-"
    using Test
else
    using Base.Test
end

@test Float64 <: IntFloat64
@test Int64 <: IntFloat64
@test !(Float64 <: IntFloat32)

@test SysInt >: Int16
@test Integer32 == Union{Int32, UInt32}

@test bytes2Int( sizeof(Int16) ) === Int16
@test bytes2UInt( sizeof(UInt32) ) === UInt32
@test bytes2Float( sizeof(Float64) ) === Float64