text:
using Base.Iterators

abstract type Part1 end
abstract type Part2 end

function update!(::Type{Part1}, grid::AbstractMatrix{Char})
    previous = copy(grid)
    updated = false
    for j in 1:size(grid, 2), i in 1:size(grid, 1)
        middle = previous[i, j]
        switch = true
        numoccupied = 0
        for (δj, δi) in product([-1, 0, 1], [-1, 0, 1])
            δi == δj == 0 && continue
            try
                neighbor = previous[i + δi, j + δj]
                if neighbor == '#'
                    middle == 'L' && (switch = false)
                    middle == '#' && (numoccupied += 1)
                end
            catch
            end
        end
        middle == 'L' && switch && (grid[i, j] = '#'; updated = true)
        middle == '#' && numoccupied ≥ 4 && (grid[i, j] = 'L'; updated = true)
    end
    updated
end

const Site = Tuple{Int, Int}

function visibility(init::AbstractMatrix{Char})
    N, M = size(init)

    visible = Dict{Site, Vector{Site}}()

    for j in 1:M, i in 1:N
        if init[i, j] != '.'
            for (δj, δi) in product([-1, 0, 1], [-1, 0, 1])
                δi == δj == 0 && continue
                u, v = i + δi, j + δj
                try
                    while init[u, v] == '.'
                        u += δi
                        v += δj
                    end
                    if haskey(visible, (i, j))
                        push!(visible[(i, j)], (u, v))
                    else
                        visible[(i, j)] = [(u, v)]
                    end
                catch
                end
            end
        end
    end
    visible
end

function update!(::Type{Part2}, grid::AbstractMatrix{Char}, visible)
    previous = copy(grid)
    updated = false
    for (i, j) in keys(visible)
        if previous[i, j] == 'L'
            switch = true
            for (u, v) in visible[(i, j)]
                if previous[u, v] == '#'
                    switch = false
                end
            end
            switch && (grid[i, j] = '#'; updated = true)
        elseif previous[i, j] == '#'
            numoccupied = 0
            for (u, v) in visible[(i, j)]
                if previous[u, v] == '#'
                    numoccupied += 1
                end
            end
            numoccupied ≥ 5 && (grid[i, j] = 'L'; updated = true)
        end
    end
    updated
end

function part1(grid)
    while update!(Part1, grid) end
    sum(grid .== '#')
end

function part2(grid)
    visible = visibility(grid)
    while update!(Part2, grid, visible) end
    sum(grid .== '#')
end

let init = hcat(collect.(readlines("input.txt"))...)
    @info "Solution" part1(init) part2(init)
end

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
export stnKD, PointNet

stnKD(K::Int) = Chain(
    Conv((1,), K => 64, relu),
    BatchNorm(64),
    Conv((1,), 64 => 128, relu),
    BatchNorm(128),
    Conv((1,), 128 => 1024, relu),
    BatchNorm(1024),
    x -> maximum(x, dims = 1),
    x -> reshape(x, :, size(x, 3)),
    Dense(1024, 512, relu),
    #     BatchNorm(512),
    Dense(512, 256, relu),
    BatchNorm(256),
    Dense(256, K * K),
    x -> reshape(x, K, K, size(x, 2)),
    # x -> x .+ I, #TODO: add identity matrix compatible with gpu
    x -> PermutedDimsArray(x, (2, 1, 3)),
)

"""
    PointNet(num_classes::Int=10, hidden_dims::Int=64)

Flux implementation of PointNet classification model.

### Fields:

- `num_classes` - Number of classes in dataset.
- `hidden_dims` - Hiddem dimension in PointNet model.

"""
struct PointNet
    stn::Any
    fstn::Any
    conv_block1::Any
    feat::Any
    cls::Any
end

function PointNet(num_classes::Int = 10, K::Int = 64)
    stn = stnKD(3)
    fstn = stnKD(K)
    conv_block1 = conv_bn_block(3, 64, (1,))
    feat = Chain(
        Conv((1,), 64 => 128, relu),
        BatchNorm(128),
        Conv((1,), 128 => 1024),
        BatchNorm(1024),
        x -> maximum(x, dims = 1),
        x -> reshape(x, 1024, :),
        Dense(1024, 512, relu),
        BatchNorm(512),
        Dense(512, 256, relu),
        Dropout(0.4),
        BatchNorm(256),
    )
    cls = Dense(256, num_classes, relu)
    PointNet(stn, fstn, conv_block1, feat, cls)
end

function (m::PointNet)(X)

    # X: [3, N, B]

    X = permutedims(X, (2, 1, 3))
    # X: [N, 3, B]

    X = Flux.batched_mul(X, m.stn(X))
    # X: [3, 3, B]

    X = m.conv_block1(X)
    # X: [3, 64, B]

    X = batched_mul(X, m.fstn(X))
    # X: [3, 64, B]

    X = m.feat(X)
    # X: [256, B]

    X = m.cls(X)
    # X: [num_classes, B]

    return softmax(X, dims = 1)
end

@functor PointNet

ext: jl


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

function vedaArgsCreate(args)
    ccall((:vedaArgsCreate, libveda), VEDAresult, (Ptr{VEDAargs},), args)
end

function vedaArgsDestroy(args)
    ccall((:vedaArgsDestroy, libveda), VEDAresult, (VEDAargs,), args)
end

function vedaArgsSetF32(args, idx, value)
    ccall((:vedaArgsSetF32, libveda), VEDAresult, (VEDAargs, Cint, Cfloat), args, idx, value)
end

function vedaArgsSetF64(args, idx, value)
    ccall((:vedaArgsSetF64, libveda), VEDAresult, (VEDAargs, Cint, Cdouble), args, idx, value)
end

function vedaArgsSetI16(args, idx, value)
    ccall((:vedaArgsSetI16, libveda), VEDAresult, (VEDAargs, Cint, Int16), args, idx, value)
end

function vedaArgsSetI32(args, idx, value)
    ccall((:vedaArgsSetI32, libveda), VEDAresult, (VEDAargs, Cint, Int32), args, idx, value)
end

function vedaArgsSetI64(args, idx, value)
    ccall((:vedaArgsSetI64, libveda), VEDAresult, (VEDAargs, Cint, Int64), args, idx, value)
end

function vedaArgsSetI8(args, idx, value)
    ccall((:vedaArgsSetI8, libveda), VEDAresult, (VEDAargs, Cint, Int8), args, idx, value)
end

function vedaArgsSetPtr(args, idx, value)
    ccall((:vedaArgsSetPtr, libveda), VEDAresult, (VEDAargs, Cint, VEDAdeviceptr), args, idx, value)
end

function vedaArgsSetVPtr(args, idx, value)
    ccall((:vedaArgsSetVPtr, libveda), VEDAresult, (VEDAargs, Cint, VEDAdeviceptr), args, idx, value)
end

function vedaArgsSetHMEM(args, idx, value)
    ccall((:vedaArgsSetHMEM, libveda), VEDAresult, (VEDAargs, Cint, Ptr{Cvoid}), args, idx, value)
end

function vedaArgsSetStack(args, idx, ptr, intent, size)
    ccall((:vedaArgsSetStack, libveda), VEDAresult, (VEDAargs, Cint, Ptr{Cvoid}, VEDAargs_intent, Csize_t), args, idx, ptr, intent, size)
end

function vedaArgsSetU16(args, idx, value)
    ccall((:vedaArgsSetU16, libveda), VEDAresult, (VEDAargs, Cint, UInt16), args, idx, value)
end

function vedaArgsSetU32(args, idx, value)
    ccall((:vedaArgsSetU32, libveda), VEDAresult, (VEDAargs, Cint, UInt32), args, idx, value)
end

function vedaArgsSetU64(args, idx, value)
    ccall((:vedaArgsSetU64, libveda), VEDAresult, (VEDAargs, Cint, UInt64), args, idx, value)
end

function vedaArgsSetU8(args, idx, value)
    ccall((:vedaArgsSetU8, libveda), VEDAresult, (VEDAargs, Cint, UInt8), args, idx, value)
end

function vedaCtxCreate(pctx, mode, dev)
    ccall((:vedaCtxCreate, libveda), VEDAresult, (Ptr{VEDAcontext}, Cint, VEDAdevice), pctx, mode, dev)
end

function vedaCtxDestroy(ctx)
    ccall((:vedaCtxDestroy, libveda), VEDAresult, (VEDAcontext,), ctx)
end

function vedaCtxGet(ctx, device)
    ccall((:vedaCtxGet, libveda), VEDAresult, (Ptr{VEDAcontext}, VEDAdevice), ctx, device)
end

function vedaCtxGetApiVersion(ctx, version)
    ccall((:vedaCtxGetApiVersion, libveda), VEDAresult, (VEDAcontext, Ptr{UInt32}), ctx, version)
end

function vedaCtxGetCurrent(pctx)
    ccall((:vedaCtxGetCurrent, libveda), VEDAresult, (Ptr{VEDAcontext},), pctx)
end

function vedaCtxGetDevice(device)
    ccall((:vedaCtxGetDevice, libveda), VEDAresult, (Ptr{VEDAdevice},), device)
end

function vedaCtxPopCurrent(pctx)
    ccall((:vedaCtxPopCurrent, libveda), VEDAresult, (Ptr{VEDAcontext},), pctx)
end

function vedaCtxPushCurrent(ctx)
    ccall((:vedaCtxPushCurrent, libveda), VEDAresult, (VEDAcontext,), ctx)
end

function vedaCtxSetCurrent(ctx)
    ccall((:vedaCtxSetCurrent, libveda), VEDAresult, (VEDAcontext,), ctx)
end

function vedaCtxStreamCnt(cnt)
    ccall((:vedaCtxStreamCnt, libveda), VEDAresult, (Ptr{Cint},), cnt)
end

function vedaCtxSynchronize()
    ccall((:vedaCtxSynchronize, libveda), VEDAresult, ())
end

function vedaDeviceDistance(distance, devA, devB)
    ccall((:vedaDeviceDistance, libveda), VEDAresult, (Ptr{Cfloat}, VEDAdevice, VEDAdevice), distance, devA, devB)
end

function vedaDeviceGet(device, ordinal)
    ccall((:vedaDeviceGet, libveda), VEDAresult, (Ptr{VEDAdevice}, Cint), device, ordinal)
end

function vedaDeviceGetAVEOId(id, dev)
    ccall((:vedaDeviceGetAVEOId, libveda), VEDAresult, (Ptr{Cint}, VEDAdevice), id, dev)
end

function vedaDeviceGetAttribute(pi, attrib, dev)
    ccall((:vedaDeviceGetAttribute, libveda), VEDAresult, (Ptr{Cint}, VEDAdevice_attribute, VEDAdevice), pi, attrib, dev)
end

function vedaDeviceGetCount(count)
    ccall((:vedaDeviceGetCount, libveda), VEDAresult, (Ptr{Cint},), count)
end

function vedaDeviceGetCurrent(current, dev)
    ccall((:vedaDeviceGetCurrent, libveda), VEDAresult, (Ptr{Cfloat}, VEDAdevice), current, dev)
end

function vedaDeviceGetCurrentEdge(current, dev)
    ccall((:vedaDeviceGetCurrentEdge, libveda), VEDAresult, (Ptr{Cfloat}, VEDAdevice), current, dev)
end

function vedaDeviceGetNUMAId(id, dev)
    ccall((:vedaDeviceGetNUMAId, libveda), VEDAresult, (Ptr{Cint}, VEDAdevice), id, dev)
end

function vedaDeviceGetName(name, len, dev)
    ccall((:vedaDeviceGetName, libveda), VEDAresult, (Cstring, Cint, VEDAdevice), name, len, dev)
end

function vedaDeviceGetPhysicalId(id, dev)
    ccall((:vedaDeviceGetPhysicalId, libveda), VEDAresult, (Ptr{Cint}, VEDAdevice), id, dev)
end

function vedaDeviceGetPower(power, dev)
    ccall((:vedaDeviceGetPower, libveda), VEDAresult, (Ptr{Cfloat}, VEDAdevice), power, dev)
end

function vedaDeviceGetTemp(tempC, coreIdx, dev)
    ccall((:vedaDeviceGetTemp, libveda), VEDAresult, (Ptr{Cfloat}, Cint, VEDAdevice), tempC, coreIdx, dev)
end

function vedaDeviceGetVoltage(voltage, dev)
    ccall((:vedaDeviceGetVoltage, libveda), VEDAresult, (Ptr{Cfloat}, VEDAdevice), voltage, dev)
end

function vedaDeviceGetVoltageEdge(voltage, dev)
    ccall((:vedaDeviceGetVoltageEdge, libveda), VEDAresult, (Ptr{Cfloat}, VEDAdevice), voltage, dev)
end

function vedaDevicePrimaryCtxGetState(dev, flags, active)
    ccall((:vedaDevicePrimaryCtxGetState, libveda), VEDAresult, (VEDAdevice, Ptr{UInt32}, Ptr{Cint}), dev, flags, active)
end

function vedaDevicePrimaryCtxRelease(dev)
    ccall((:vedaDevicePrimaryCtxRelease, libveda), VEDAresult, (VEDAdevice,), dev)
end

function vedaDevicePrimaryCtxReset(dev)
    ccall((:vedaDevicePrimaryCtxReset, libveda), VEDAresult, (VEDAdevice,), dev)
end

function vedaDevicePrimaryCtxRetain(pctx, dev)
    ccall((:vedaDevicePrimaryCtxRetain, libveda), VEDAresult, (Ptr{VEDAcontext}, VEDAdevice), pctx, dev)
end

function vedaDevicePrimaryCtxSetFlags(dev, flags)
    ccall((:vedaDevicePrimaryCtxSetFlags, libveda), VEDAresult, (VEDAdevice, UInt32), dev, flags)
end

function vedaDeviceTotalMem(bytes, dev)
    ccall((:vedaDeviceTotalMem, libveda), VEDAresult, (Ptr{Csize_t}, VEDAdevice), bytes, dev)
end

function vedaDriverGetVersion(str)
    ccall((:vedaDriverGetVersion, libveda), VEDAresult, (Ptr{Cstring},), str)
end

function vedaExit()
    ccall((:vedaExit, libveda), VEDAresult, ())
end

function vedaGetErrorName(error, pStr)
    ccall((:vedaGetErrorName, libveda), VEDAresult, (VEDAresult, Ptr{Cstring}), error, pStr)
end

function vedaGetErrorString(error, pStr)
    ccall((:vedaGetErrorString, libveda), VEDAresult, (VEDAresult, Ptr{Cstring}), error, pStr)
end

function vedaGetVersion(str)
    ccall((:vedaGetVersion, libveda), VEDAresult, (Ptr{Cstring},), str)
end

function vedaInit(Flags)
    ccall((:vedaInit, libveda), VEDAresult, (UInt32,), Flags)
end

function vedaLaunchHostFunc(stream, fn, userData)
    ccall((:vedaLaunchHostFunc, libveda), VEDAresult, (VEDAstream, VEDAhost_function, Ptr{Cvoid}), stream, fn, userData)
end

function vedaLaunchHostFuncEx(stream, fn, userData, result)
    ccall((:vedaLaunchHostFuncEx, libveda), VEDAresult, (VEDAstream, VEDAhost_function, Ptr{Cvoid}, Ptr{UInt64}), stream, fn, userData, result)
end

function vedaLaunchKernel(f, stream, arg1)
    ccall((:vedaLaunchKernel, libveda), VEDAresult, (VEDAfunction, VEDAstream, VEDAargs), f, stream, arg1)
end

function vedaLaunchKernelEx(f, stream, arg1, destroyArgs, result)
    ccall((:vedaLaunchKernelEx, libveda), VEDAresult, (VEDAfunction, VEDAstream, VEDAargs, Cint, Ptr{UInt64}), f, stream, arg1, destroyArgs, result)
end

function vedaMemAlloc(ptr, size)
    ccall((:vedaMemAlloc, libveda), VEDAresult, (Ptr{VEDAdeviceptr}, Csize_t), ptr, size)
end

function vedaMemAllocAsync(ptr, size, stream)
    ccall((:vedaMemAllocAsync, libveda), VEDAresult, (Ptr{VEDAdeviceptr}, Csize_t, VEDAstream), ptr, size, stream)
end

function vedaMemAllocHost(pp, bytesiz)
    ccall((:vedaMemAllocHost, libveda), VEDAresult, (Ptr{Ptr{Cvoid}}, Csize_t), pp, bytesiz)
end

function vedaMemAllocPitch(ptr, pPitch, WidthInBytes, Height, ElementSizeByte)
    ccall((:vedaMemAllocPitch, libveda), VEDAresult, (Ptr{VEDAdeviceptr}, Ptr{Csize_t}, Csize_t, Csize_t, UInt32), ptr, pPitch, WidthInBytes, Height, ElementSizeByte)
end

function vedaMemAllocPitchAsync(ptr, pPitch, WidthInBytes, Height, ElementSizeByte, stream)
    ccall((:vedaMemAllocPitchAsync, libveda), VEDAresult, (Ptr{VEDAdeviceptr}, Ptr{Csize_t}, Csize_t, Csize_t, UInt32, VEDAstream), ptr, pPitch, WidthInBytes, Height, ElementSizeByte, stream)
end

function vedaMemFree(ptr)
    ccall((:vedaMemFree, libveda), VEDAresult, (VEDAdeviceptr,), ptr)
end

function vedaMemFreeAsync(ptr, stream)
    ccall((:vedaMemFreeAsync, libveda), VEDAresult, (VEDAdeviceptr, VEDAstream), ptr, stream)
end

function vedaMemFreeHost(ptr)
    ccall((:vedaMemFreeHost, libveda), VEDAresult, (Ptr{Cvoid},), ptr)
end

function vedaMemGetAddressRange(base, size, ptr)
    ccall((:vedaMemGetAddressRange, libveda), VEDAresult, (Ptr{VEDAdeviceptr}, Ptr{Csize_t}, VEDAdeviceptr), base, size, ptr)
end

function vedaMemGetDevice(dev, ptr)
    ccall((:vedaMemGetDevice, libveda), VEDAresult, (Ptr{VEDAdevice}, VEDAdeviceptr), dev, ptr)
end

function vedaMemGetInfo(free, total)
    ccall((:vedaMemGetInfo, libveda), VEDAresult, (Ptr{Csize_t}, Ptr{Csize_t}), free, total)
end

function vedaMemHMEM(ptr, vptr)
    ccall((:vedaMemHMEM, libveda), VEDAresult, (Ptr{Ptr{Cvoid}}, VEDAdeviceptr), ptr, vptr)
end

function vedaMemHMEMSize(ptr, size, vptr)
    ccall((:vedaMemHMEMSize, libveda), VEDAresult, (Ptr{Ptr{Cvoid}}, Ptr{Csize_t}, VEDAdeviceptr), ptr, size, vptr)
end

function vedaMemPtr(ptr, vptr)
    ccall((:vedaMemPtr, libveda), VEDAresult, (Ptr{Ptr{Cvoid}}, VEDAdeviceptr), ptr, vptr)
end

function vedaMemPtrSize(ptr, size, vptr)
    ccall((:vedaMemPtrSize, libveda), VEDAresult, (Ptr{Ptr{Cvoid}}, Ptr{Csize_t}, VEDAdeviceptr), ptr, size, vptr)
end

function vedaMemReport()
    ccall((:vedaMemReport, libveda), VEDAresult, ())
end

function vedaMemSize(size, ptr)
    ccall((:vedaMemSize, libveda), VEDAresult, (Ptr{Csize_t}, VEDAdeviceptr), size, ptr)
end

function vedaMemSwap(A, B)
    ccall((:vedaMemSwap, libveda), VEDAresult, (VEDAdeviceptr, VEDAdeviceptr), A, B)
end

function vedaMemSwapAsync(A, B, hStream)
    ccall((:vedaMemSwapAsync, libveda), VEDAresult, (VEDAdeviceptr, VEDAdeviceptr, VEDAstream), A, B, hStream)
end

function vedaMemcpy(dst, src, ByteCount)
    ccall((:vedaMemcpy, libveda), VEDAresult, (VEDAdeviceptr, VEDAdeviceptr, Csize_t), dst, src, ByteCount)
end

function vedaMemcpyAsync(dst, src, ByteCount, hStream)
    ccall((:vedaMemcpyAsync, libveda), VEDAresult, (VEDAdeviceptr, VEDAdeviceptr, Csize_t, VEDAstream), dst, src, ByteCount, hStream)
end

function vedaMemcpyDtoD(dstDevice, srcDevice, ByteCount)
    ccall((:vedaMemcpyDtoD, libveda), VEDAresult, (VEDAdeviceptr, VEDAdeviceptr, Csize_t), dstDevice, srcDevice, ByteCount)
end

function vedaMemcpyDtoDAsync(dstDevice, srcDevice, ByteCount, hStream)
    ccall((:vedaMemcpyDtoDAsync, libveda), VEDAresult, (VEDAdeviceptr, VEDAdeviceptr, Csize_t, VEDAstream), dstDevice, srcDevice, ByteCount, hStream)
end

function vedaMemcpyDtoH(dstHost, srcDevice, ByteCount)
    ccall((:vedaMemcpyDtoH, libveda), VEDAresult, (Ptr{Cvoid}, VEDAdeviceptr, Csize_t), dstHost, srcDevice, ByteCount)
end

function vedaMemcpyDtoHAsync(dstHost, srcDevice, ByteCount, hStream)
    ccall((:vedaMemcpyDtoHAsync, libveda), VEDAresult, (Ptr{Cvoid}, VEDAdeviceptr, Csize_t, VEDAstream), dstHost, srcDevice, ByteCount, hStream)
end

function vedaMemcpyHtoD(dstDevice, srcHost, ByteCount)
    ccall((:vedaMemcpyHtoD, libveda), VEDAresult, (VEDAdeviceptr, Ptr{Cvoid}, Csize_t), dstDevice, srcHost, ByteCount)
end

function vedaMemcpyHtoDAsync(dstDevice, srcHost, ByteCount, hStream)
    ccall((:vedaMemcpyHtoDAsync, libveda), VEDAresult, (VEDAdeviceptr, Ptr{Cvoid}, Csize_t, VEDAstream), dstDevice, srcHost, ByteCount, hStream)
end

function vedaMemsetD128(dstDevice, x, y, N)
    ccall((:vedaMemsetD128, libveda), VEDAresult, (VEDAdeviceptr, UInt64, UInt64, Csize_t), dstDevice, x, y, N)
end

function vedaMemsetD128Async(dstDevice, x, u, N, hStream)
    ccall((:vedaMemsetD128Async, libveda), VEDAresult, (VEDAdeviceptr, UInt64, UInt64, Csize_t, VEDAstream), dstDevice, x, u, N, hStream)
end

function vedaMemsetD16(dstDevice, us, N)
    ccall((:vedaMemsetD16, libveda), VEDAresult, (VEDAdeviceptr, UInt16, Csize_t), dstDevice, us, N)
end

function vedaMemsetD16Async(dstDevice, us, N, hStream)
    ccall((:vedaMemsetD16Async, libveda), VEDAresult, (VEDAdeviceptr, UInt16, Csize_t, VEDAstream), dstDevice, us, N, hStream)
end

function vedaMemsetD2D128(dstDevice, dstPitch, x, y, Width, Height)
    ccall((:vedaMemsetD2D128, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt64, UInt64, Csize_t, Csize_t), dstDevice, dstPitch, x, y, Width, Height)
end

function vedaMemsetD2D128Async(dstDevice, dstPitch, x, y, Width, Height, hStream)
    ccall((:vedaMemsetD2D128Async, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt64, UInt64, Csize_t, Csize_t, VEDAstream), dstDevice, dstPitch, x, y, Width, Height, hStream)
end

function vedaMemsetD2D16(dstDevice, dstPitch, us, Width, Height)
    ccall((:vedaMemsetD2D16, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt16, Csize_t, Csize_t), dstDevice, dstPitch, us, Width, Height)
end

function vedaMemsetD2D16Async(dstDevice, dstPitch, us, Width, Height, hStream)
    ccall((:vedaMemsetD2D16Async, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt16, Csize_t, Csize_t, VEDAstream), dstDevice, dstPitch, us, Width, Height, hStream)
end

function vedaMemsetD2D32(dstDevice, dstPitch, ui, Width, Height)
    ccall((:vedaMemsetD2D32, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt32, Csize_t, Csize_t), dstDevice, dstPitch, ui, Width, Height)
end

function vedaMemsetD2D32Async(dstDevice, dstPitch, ui, Width, Height, hStream)
    ccall((:vedaMemsetD2D32Async, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt32, Csize_t, Csize_t, VEDAstream), dstDevice, dstPitch, ui, Width, Height, hStream)
end

function vedaMemsetD2D64(dstDevice, dstPitch, ul, Width, Height)
    ccall((:vedaMemsetD2D64, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt64, Csize_t, Csize_t), dstDevice, dstPitch, ul, Width, Height)
end

function vedaMemsetD2D64Async(dstDevice, dstPitch, ul, Width, Height, hStream)
    ccall((:vedaMemsetD2D64Async, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt64, Csize_t, Csize_t, VEDAstream), dstDevice, dstPitch, ul, Width, Height, hStream)
end

function vedaMemsetD2D8(dstDevice, dstPitch, uc, Width, Height)
    ccall((:vedaMemsetD2D8, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt8, Csize_t, Csize_t), dstDevice, dstPitch, uc, Width, Height)
end

function vedaMemsetD2D8Async(dstDevice, dstPitch, uc, Width, Height, hStream)
    ccall((:vedaMemsetD2D8Async, libveda), VEDAresult, (VEDAdeviceptr, Csize_t, UInt8, Csize_t, Csize_t, VEDAstream), dstDevice, dstPitch, uc, Width, Height, hStream)
end

function vedaMemsetD32(dstDevice, ui, N)
    ccall((:vedaMemsetD32, libveda), VEDAresult, (VEDAdeviceptr, UInt32, Csize_t), dstDevice, ui, N)
end

function vedaMemsetD32Async(dstDevice, ui, N, hStream)
    ccall((:vedaMemsetD32Async, libveda), VEDAresult, (VEDAdeviceptr, UInt32, Csize_t, VEDAstream), dstDevice, ui, N, hStream)
end

function vedaMemsetD64(dstDevice, ui, N)
    ccall((:vedaMemsetD64, libveda), VEDAresult, (VEDAdeviceptr, UInt64, Csize_t), dstDevice, ui, N)
end

function vedaMemsetD64Async(dstDevice, ui, N, hStream)
    ccall((:vedaMemsetD64Async, libveda), VEDAresult, (VEDAdeviceptr, UInt64, Csize_t, VEDAstream), dstDevice, ui, N, hStream)
end

function vedaMemsetD8(dstDevice, uc, N)
    ccall((:vedaMemsetD8, libveda), VEDAresult, (VEDAdeviceptr, UInt8, Csize_t), dstDevice, uc, N)
end

function vedaMemsetD8Async(dstDevice, uc, N, hStream)
    ccall((:vedaMemsetD8Async, libveda), VEDAresult, (VEDAdeviceptr, UInt8, Csize_t, VEDAstream), dstDevice, uc, N, hStream)
end

function vedaModuleGetFunction(hfunc, hmod, name)
    ccall((:vedaModuleGetFunction, libveda), VEDAresult, (Ptr{VEDAfunction}, VEDAmodule, Cstring), hfunc, hmod, name)
end

function vedaModuleLoad(_module, fname)
    ccall((:vedaModuleLoad, libveda), VEDAresult, (Ptr{VEDAmodule}, Cstring), _module, fname)
end

function vedaModuleUnload(hmod)
    ccall((:vedaModuleUnload, libveda), VEDAresult, (VEDAmodule,), hmod)
end

function vedaStreamAddCallback(stream, callback, userData, flags)
    ccall((:vedaStreamAddCallback, libveda), VEDAresult, (VEDAstream, VEDAstream_callback, Ptr{Cvoid}, UInt32), stream, callback, userData, flags)
end

function vedaStreamGetFlags(hStream, flags)
    ccall((:vedaStreamGetFlags, libveda), VEDAresult, (VEDAstream, Ptr{UInt32}), hStream, flags)
end

function vedaStreamQuery(hStream)
    ccall((:vedaStreamQuery, libveda), VEDAresult, (VEDAstream,), hStream)
end

function vedaStreamSynchronize(hStream)
    ccall((:vedaStreamSynchronize, libveda), VEDAresult, (VEDAstream,), hStream)
end
# Julia wrapper for header: veda_types.h
# Automatically generated using Clang.jl

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

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

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

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


text:
@testset "DAQ" begin
  include("DummyDAQTest.jl")
  include("SimpleSimulatedDAQTest.jl")
end


text:
#
# types.jl --
#
# Definitions of types.
#
#-------------------------------------------------------------------------------
#
# This file is part of "AndorCameras.jl" released under the MIT license.
#
# Copyright (C) 2017-2019, Éric Thiébaut.
#

"""

`AndorError` is an exception sub-type to report errors occuring in calls to the
Sofware Develement Kit (SDK) of Andor cameras.

"""
struct AndorError <: Exception
    func::Symbol
    code::AT.STATUS
    AndorError(func::Symbol, code::Integer) = new(func, code)
    AndorError(func::AbstractString, code::Integer) = new(Symbol(func), code)
end

"""

`AndorCameraModel` is used to wrap a constant identifier to quickly identify a
specific camera model.

"""
struct AndorCameraModel
    id::Int
end
Base.:(==)(a::AndorCameraModel, b::AndorCameraModel) = a.id == b.id
const _UNKNOWN_MODEL  = AndorCameraModel(0)
const _SIM_CAM_MODEL  = AndorCameraModel(1)
const _ZYLA_MODEL     = AndorCameraModel(2)
const _ZYLA_USB_MODEL = AndorCameraModel(3)

# Camera structure.
mutable struct Camera <: ScientificCamera
    state::Int
    bufs::Vector{Vector{UInt8}}        # buffers for the frame grabber
    lastimg::Array{T,2} where {T}      # last image
    bytesperline::Int
    clockfrequency::Int
    model::AndorCameraModel
    handle::AT.HANDLE
    mono12packed::Bool
    Camera() = new(0,
                   Vector{Vector{UInt8}}(undef, 0),
                   Array{UInt8,2}(undef, 0, 0),
                   0, 0, _UNKNOWN_MODEL, AT.HANDLE_UNINITIALISED, false)
end

const AndorCamera = Camera

abstract type AbstractFeature end

struct CommandFeature <: AbstractFeature
    name::Vector{AT.WCHAR}
    CommandFeature(sym::Symbol) = new(AT.widestring(sym))
end

struct BooleanFeature <: AbstractFeature
    name::Vector{AT.WCHAR}
    BooleanFeature(sym::Symbol) = new(AT.widestring(sym))
end

struct EnumeratedFeature <: AbstractFeature
    name::Vector{AT.WCHAR}
    EnumeratedFeature(sym::Symbol) = new(AT.widestring(sym))
end

struct IntegerFeature <: AbstractFeature
    name::Vector{AT.WCHAR}
    IntegerFeature(sym::Symbol) = new(AT.widestring(sym))
end

struct FloatingPointFeature <: AbstractFeature
    name::Vector{AT.WCHAR}
    FloatingPointFeature(sym::Symbol) = new(AT.widestring(sym))
end

struct StringFeature <: AbstractFeature
    name::Vector{AT.WCHAR}
    StringFeature(sym::Symbol) = new(AT.widestring(sym))
end

struct BooleanOrEnumeratedFeature <: AbstractFeature
    name::Vector{AT.WCHAR}
    BooleanOrEnumeratedFeature(sym::Symbol) = new(AT.widestring(sym))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ dcab9a18-5d58-4632-a325-d16af9c79344
using Pkg, DrWatson

# ╔═╡ c95f3f02-f179-4f81-9c53-6e6b87f25613
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers, Symbolics, SymbolicUtils
	using StaticArrays
end

# ╔═╡ 75361cb8-3c45-4475-94dc-bf8093114613
md" ##### How to build a Lagrangian polynomial in expression form."

# ╔═╡ c9036632-adb9-11eb-33d1-c5193729e72a
function lagrangianpolynomial1(x::Vector{Float64},
		y::Vector{Float64}, xi::Vector{Float64})

	np = length(x)
 	yi = Vector{Float64}(undef, length(xi))
 	m = sort(hcat(x, y), dims=1)
	for l in 1:length(xi)
		for i in 1:np
			term = 1.0
			for j in 1:np
				j !== i && (term *= (xi[l] - m[j, 1])/(m[i, 1] - m[j, 1]))
			end
			yi[l] += term * m[i, 2]
		end
	end
	columns=[]
	push!(columns, m[:,1])
	push!(columns, m[:,2])
	df = DataFrame(columns, [:x, :y])

	columns=[]
	push!(columns, xi)
	push!(columns, yi)
	dfxi = DataFrame(columns, [:xi, :yi])

	(df, dfxi)
end

# ╔═╡ 3d35a281-d906-495f-8534-a6379493665a
md" ##### Expressions from Strings:"

# ╔═╡ 298bbd55-d7b6-4229-8848-f0ea11f6448e
s = "var" * '_' * "t1"

# ╔═╡ 822ad40e-c9c1-4de9-bec4-6533f0d7ea8f
e = Meta.parse(s)

# ╔═╡ d2ade179-b33b-4c73-bebe-b362007a0221
ex = :(a+b*c+1)

# ╔═╡ 8dab1a23-9db4-46b9-bbf8-106d76de570a
typeof(ex)

# ╔═╡ 50b46461-3838-4d95-ae64-a63bf2d29b05
:(a + b*c + 1)       ==
	Meta.parse("a + b*c + 1") ==
       Expr(:call, :+, :a, Expr(:call, :*, :b, :c), 1)

# ╔═╡ 141ed7a4-fbac-40b1-9ada-7b74842b5e75
md" ##### `From the manual:` The use of $ for expression interpolation is intentionally reminiscent of string interpolation and command interpolation."

# ╔═╡ ce6ecce6-eafc-4070-b14f-45415452bb84
begin
	a = 1;
	ex1 = :($a + b)
end

# ╔═╡ e0b8881d-9dec-4ea4-9084-b7e045eeb047
ex2 = :(a in $:((1,2,3)) )

# ╔═╡ 3bb91119-11ef-4ec8-9eff-b7cf4d336417
md" ##### Splatting arguments:"

# ╔═╡ 56895695-745d-46ec-9326-3ce444d0597b
args = [:x, :y, :z];

# ╔═╡ f0aa8168-6bc9-4702-b40a-689869258741
:(f(1, $(args...)))

# ╔═╡ 1181e81e-8603-4d80-a752-1d9d1c51640d
begin
	x = :(1 + 2)
	ex3 = quote quote $x end end
end

# ╔═╡ 309baaeb-8a45-4502-9840-9d71a2b1349f
eval(ex3)

# ╔═╡ 06f51cf6-cb6e-4fce-b914-fdbf4a65ef79
eval(eval(ex3))

# ╔═╡ acbd6c5f-1b3a-4a77-a1b8-1ebf47f563ee
md" ##### Output to Terminal:"

# ╔═╡ cde57703-141a-4bd6-b411-d36ddf5eaf48
dump(Meta.parse(":(1+2)"))

# ╔═╡ 771c13c0-07a9-416c-b254-6ed9b6ee37c0
ex4 = :(x1 = 1)

# ╔═╡ e6c624d0-c468-482d-8c42-de8f44b7e7d4
try
	x1
catch exc
	exc
end

# ╔═╡ 8229d3f5-ff3d-482f-8029-a521f62bc19b
eval(ex4);

# ╔═╡ 8dabc67e-72a9-4359-9a7e-9bbaecd9f1a8
x1

# ╔═╡ 43851fcd-bd0b-47aa-8b50-83242670dfba
md" ##### Functions on expressions:"

# ╔═╡ ccbd0884-79ff-49cd-82e5-3247665a2edc
function math_expr(op, op1, op2)
	ex5 = Expr(:call, op, op1, op2)
	return ex5
end

# ╔═╡ 6133a422-d5fa-4ab6-8eb7-18de0f0c904f
ex6 = math_expr(:+, 1, Expr(:call, :*, 4, 5))

# ╔═╡ 727f59e1-89e2-4407-a32a-c2a9640ab031
eval(ex6)

# ╔═╡ 2447682c-407c-44e0-889a-45a0ea3d9843
function make_expr2(op, opr1, opr2)
	opr1f, opr2f = map(x -> isa(x, Number) ? 2*x : x, (opr1, opr2))
	retexpr = Expr(:call, op, opr1f, opr2f)
	return retexpr
end

# ╔═╡ baffafa2-14e4-4461-93c7-beeebd254d34
make_expr2(:+, 1, 2)

# ╔═╡ 76e490a7-26c2-4bab-9d04-83fd3bf9aa6e
ex7 = make_expr2(:+, 1, Expr(:call, :*, 5, 8))

# ╔═╡ 58d9f409-5da7-473d-a933-2221ae6a9da1
eval(ex7)

# ╔═╡ c8ee6c92-f51f-4cfc-b2eb-2a2eadb73be5
ex8 = make_expr2(:(.+), 1, [2, 3])

# ╔═╡ 657560a4-b7ae-4fff-8ba8-c17bfbfe9179
eval(ex8)

# ╔═╡ 2ddc5871-0735-4d52-ae8b-720c606fef3f
function make_expr3(op, opr1, opr2)
	opr1f, opr2f =
		map(x -> isa(x, Number) || isa(x, Vector{Int64}) ? 2 .* x : x, (opr1, opr2))
	retexpr = Expr(:call, op, opr1f, opr2f)
	return retexpr
end

# ╔═╡ c30be5cb-8aba-4ed3-9356-45d64c4dd236
ex9 = make_expr3(:(.+), 1, [2, 3])

# ╔═╡ 6d0fb950-9c30-4d7b-9954-a2877a90ff5e
typeof([2, 3])

# ╔═╡ fcb65821-e585-4af9-860f-6003268fdfe8
eval(ex9)

# ╔═╡ Cell order:
# ╠═dcab9a18-5d58-4632-a325-d16af9c79344
# ╠═c95f3f02-f179-4f81-9c53-6e6b87f25613
# ╟─75361cb8-3c45-4475-94dc-bf8093114613
# ╠═c9036632-adb9-11eb-33d1-c5193729e72a
# ╟─3d35a281-d906-495f-8534-a6379493665a
# ╠═298bbd55-d7b6-4229-8848-f0ea11f6448e
# ╠═822ad40e-c9c1-4de9-bec4-6533f0d7ea8f
# ╠═d2ade179-b33b-4c73-bebe-b362007a0221
# ╠═8dab1a23-9db4-46b9-bbf8-106d76de570a
# ╠═50b46461-3838-4d95-ae64-a63bf2d29b05
# ╟─141ed7a4-fbac-40b1-9ada-7b74842b5e75
# ╠═ce6ecce6-eafc-4070-b14f-45415452bb84
# ╠═e0b8881d-9dec-4ea4-9084-b7e045eeb047
# ╟─3bb91119-11ef-4ec8-9eff-b7cf4d336417
# ╠═56895695-745d-46ec-9326-3ce444d0597b
# ╠═f0aa8168-6bc9-4702-b40a-689869258741
# ╠═1181e81e-8603-4d80-a752-1d9d1c51640d
# ╠═309baaeb-8a45-4502-9840-9d71a2b1349f
# ╠═06f51cf6-cb6e-4fce-b914-fdbf4a65ef79
# ╟─acbd6c5f-1b3a-4a77-a1b8-1ebf47f563ee
# ╠═cde57703-141a-4bd6-b411-d36ddf5eaf48
# ╠═771c13c0-07a9-416c-b254-6ed9b6ee37c0
# ╠═e6c624d0-c468-482d-8c42-de8f44b7e7d4
# ╠═8229d3f5-ff3d-482f-8029-a521f62bc19b
# ╠═8dabc67e-72a9-4359-9a7e-9bbaecd9f1a8
# ╟─43851fcd-bd0b-47aa-8b50-83242670dfba
# ╠═ccbd0884-79ff-49cd-82e5-3247665a2edc
# ╠═6133a422-d5fa-4ab6-8eb7-18de0f0c904f
# ╠═727f59e1-89e2-4407-a32a-c2a9640ab031
# ╠═2447682c-407c-44e0-889a-45a0ea3d9843
# ╠═baffafa2-14e4-4461-93c7-beeebd254d34
# ╠═76e490a7-26c2-4bab-9d04-83fd3bf9aa6e
# ╠═58d9f409-5da7-473d-a933-2221ae6a9da1
# ╠═c8ee6c92-f51f-4cfc-b2eb-2a2eadb73be5
# ╠═657560a4-b7ae-4fff-8ba8-c17bfbfe9179
# ╠═2ddc5871-0735-4d52-ae8b-720c606fef3f
# ╠═c30be5cb-8aba-4ed3-9356-45d64c4dd236
# ╠═6d0fb950-9c30-4d7b-9954-a2877a90ff5e
# ╠═fcb65821-e585-4af9-860f-6003268fdfe8


text:
using EzXML

function aasearch(filename::String)
    pdf = readpdf("$filename.pdf")
    xml = readxml("$filename.xml")
    tree
    leaves = filter(isleaf, tree)
    for leaf in leaves
        leaf.data
    end
end

function simsearch(doc::Vector{PDContent}, query::String)

end

function featvec()
end

export dictmatch
function dictmatch!(pdf::Vector{PDContent}, dictfile::String)
    patterns = Regex[]
    for line in open(readlines, dictfile)
        isempty(line) && continue
        if startswith(line,"/") && endswith(line,"/")
            push!(patterns, Regex(line[2:end-1]))
        else

        end
    end

    text = join(map(c -> c.value, pdf))
    for p in patterns
        for m in eachmatch(p, text)
            println(m.offset)
        end
    end
    pdf
end


text:
#=
helpers for sparse factorizations and linear solves
=#

import SparseArrays.SparseMatrixCSC
import SuiteSparse
const SuiteSparseInt = SuiteSparse.CHOLMOD.SuiteSparse_long

#=
nonsymmetric
=#

abstract type SparseNonSymCache{T <: Real} end

# only works with Float64
mutable struct UMFPACKNonSymCache{Float64} <: SparseNonSymCache{Float64}
    analyzed::Bool
    umfpack::SuiteSparse.UMFPACK.UmfpackLU
    function UMFPACKNonSymCache{Float64}()
        cache = new{Float64}()
        cache.analyzed = false
        return cache
    end
end

# restrict int type to SuiteSparse_long
int_type(::UMFPACKNonSymCache) = SuiteSparseInt

function update_fact(
    cache::UMFPACKNonSymCache{Float64},
    A::SparseMatrixCSC{Float64, SuiteSparseInt},
    )
    if !cache.analyzed
        cache.umfpack = lu(A) # symbolic and numeric factorization
        cache.analyzed = true
    else
        lu!(cache.umfpack, A, check = false)
    end
    return
end

function inv_prod(
    cache::UMFPACKNonSymCache{Float64},
    x::Vector{Float64},
    A::SparseMatrixCSC{Float64, SuiteSparseInt},
    b::Vector{Float64},
    )
    ldiv!(x, cache.umfpack, b) # does not repeat symbolic or numeric factorization
    return x
end

# default to UMFPACK
SparseNonSymCache{Float64}() = UMFPACKNonSymCache{Float64}()

#=
symmetric
=#

abstract type SparseSymCache{T <: Real} end

# helper for symmetric solvers that need nonzero diagonal
diag_min(::SparseSymCache{T}) where T = zero(T)

# only works with Float64
mutable struct CHOLMODSymCache{Float64} <: SparseSymCache{Float64}
    analyzed::Bool
    cholmod::SuiteSparse.CHOLMOD.Factor
    function CHOLMODSymCache{Float64}()
        cache = new{Float64}()
        cache.analyzed = false
        return cache
    end
end
int_type(::CHOLMODSymCache) = SuiteSparseInt

# CHOLMOD needs nonzero diagonal
diag_min(::SparseSymCache{Float64}) = sqrt(eps())

function update_fact(
    cache::CHOLMODSymCache{Float64},
    A::SparseMatrixCSC{Float64, SuiteSparseInt},
    )
    A_symm = Symmetric(A, :L)

    if !cache.analyzed
        cache.cholmod = SuiteSparse.CHOLMOD.ldlt(A_symm, check = false)
        cache.analyzed = true
    else
        ldlt!(cache.cholmod, A_symm, check = false)
    end

    if !issuccess(cache.cholmod)
        @warn("numerical failure: sparse factorization failed")
        ldlt!(cache.cholmod, A_symm, shift = 1e-4, check = false)
        if !issuccess(cache.cholmod)
            @warn("numerical failure: could not fix sparse factorization failure")
        end
    end
    return
end

function inv_prod(
    cache::CHOLMODSymCache{Float64},
    x::Vector{Float64},
    A::SparseMatrixCSC{Float64, SuiteSparseInt},
    b::Vector{Float64},
    )
    x .= cache.cholmod \ b # TODO try to make this in-place
    return x
end

# default to CHOLMOD
SparseSymCache{Float64}() = CHOLMODSymCache{Float64}()

#=
helpers
=#

free_memory(::Union{
    UMFPACKNonSymCache{Float64},
    CHOLMODSymCache{Float64},
    }) = nothing


text:
#=
example5:
- Julia version:
- Author: Ashley Sommer
- Date: 2020-01-03

This example gets some Sensor Observations from a given Cosmoz station, with observations aggregated to 24H steps
=#

using Printf
using JSON
using HTTP
using Dates

const COSMOZ_API_URL = "https://esoil.io/cosmoz-data-pipeline/rest/"  # Keep the trailing slash on here.
const STATION_NUMBER = 21
const PROCESSING_LEVEL = 4  # Choose processing level 1, 2, 3, 4, or 0 (for Raw)

# Endpoint to get a station's observations is "pipeline/rest/stations/{id}/observations"
stations_endpoint = COSMOZ_API_URL * "stations"
station_endpoint = stations_endpoint * "/" * string(STATION_NUMBER)
station_obs_endpoint = station_endpoint * "/observations"

# Time Period Start Date
start_date = Dates.DateTime(2019, 1, 1) # UTC
start_date_str = Dates.format(start_date, "Y-m-dTH:M:S.000Z")  # ISO8601 Format
# Time Period End Date
end_date = Dates.DateTime(2019, 1, 31) # UTC
end_date_str = Dates.format(end_date, "Y-m-dTH:M:S.000Z")  # ISO8601 Format
# Add request query params
query_params = Dict(
    "processing_level"=>PROCESSING_LEVEL,
    "startdate"=>start_date_str,
    "enddate"=>end_date_str,
    "aggregate"=>"24h"  # Set aggregation to 24H
)

# Add a header to specifically ask for JSON output
request_headers = Dict("Accept"=>"application/json",)
# Construct a GET request, with that URL and those headers
r = HTTP.request("GET", station_obs_endpoint, request_headers, query=query_params)
payload = JSON.parse(String(r.body))
println("Showing raw payload")
println(payload)
println()
count = payload["meta"]["count"]
println("Found $count Observations for site $(string(STATION_NUMBER))")
for (i, c) in enumerate(payload["observations"])
    println("Observation $i")
    for (k,v) in c
        println("\t$(repr(k)): $(repr(v))")
    end
end


text:
module Term
# general utils
include("__text_utils.jl")
include("_ansi.jl")
include("_utils.jl")

# don't import other modules
include("measure.jl")
include("color.jl")
include("theme.jl")
include("highlight.jl")
include("console.jl")

# rely on other modules
include("markup.jl")
include("style.jl")
include("segment.jl")
include("macros.jl")

# renderables, rely heavily on other modules
include("box.jl")
include("renderables.jl")
include("layout.jl")
include("panel.jl")
include("inspect.jl")
include("errors.jl")
include("logging.jl")

include("logo.jl")

export RenderableText, Panel, TextBox
export Spacer, vLine, hLine
export theme, highlight
export inspect
export @red, @black, @green, @yellow, @blue, @magenta, @cyan, @white, @default
export @bold, @dim, @italic, @underline, @style
export tprint, install_stacktrace
export install_term_logger

# ----------------------------------- base ----------------------------------- #
import .measure
using .measure: Measure
using .consoles: Console, console, err_console, console_height, console_width

# ----------------------------------- style ---------------------------------- #
using .markup: extract_markup, MarkupTag, pairup_tags

using .color: NamedColor, BitColor, RGBColor, get_color

using .style: MarkupStyle, extract_style, apply_style

using .segment: Segment

"""
    Measure(seg::Segment)

gives the measure of a segment
"""
measure.Measure(seg::Segment) = seg.measure

"""
    Measure(segments::Vector{Segment})

gives the measure of a vector of segments
"""
function measure.Measure(segments::Vector{Segment})
    return Measure(
        max([seg.measure.w for seg in segments]...),
        sum([seg.measure.h for seg in segments]),
    )
end

# -------------------------------- renderables ------------------------------- #
using .box

using .renderables: AbstractRenderable, Renderable, RenderableText

using .layout: Padding, vstack, hstack, Spacer, vLine, hLine

using .panel: Panel, TextBox

tprint(x::AbstractString) = (println ∘ apply_style)(x)
tprint(x::AbstractRenderable) = println(x)
tpritn(x) = tprint(string(x))

# ---------------------------------- others ---------------------------------- #
using .introspection: inspect

using .errors: install_stacktrace

using .logging: install_term_logger, TermLogger

end


text:
@testset "intrinsics" begin

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

@testset "indexing" begin
    @on_device threadIdx().x
    @on_device blockDim().x
    @on_device blockIdx().x
    @on_device gridDim().x

    @on_device threadIdx().y
    @on_device blockDim().y
    @on_device blockIdx().y
    @on_device gridDim().y

    @on_device threadIdx().z
    @on_device blockDim().z
    @on_device blockIdx().z
    @on_device gridDim().z

    @testset "range metadata" begin
        foobar() = threadIdx().x
        ir = sprint(io->CUDAnative.code_llvm(io, foobar, Tuple{}))

        @test occursin(r"call .+ @llvm.nvvm.read.ptx.sreg.tid.x.+ !range", ir)
    end
end

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

@testset "math" begin
    buf = CuTestArray(Float32[0])

    function kernel(a, i)
        a[1] = CUDAnative.log10(i)
        return
    end

    @cuda kernel(buf, Float32(100))
    val = Array(buf)
    @test val[1] ≈ 2.0
end

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

@testset "I/O" begin

@testset "printing" begin
    _, out = @grab_output @on_device @cuprintf("")
    @test out == ""

    endline = Sys.iswindows() ? "\r\n" : "\n"

    _, out = @grab_output @on_device @cuprintf("Testing...\n")
    @test out == "Testing...$endline"

    # narrow integer
    _, out = @grab_output @on_device @cuprintf("Testing %d...\n", Int32(42))
    @test out == "Testing 42...$endline"

    # wide integer
    _, out = @grab_output @on_device @cuprintf("Testing %ld...\n", Int64(42))
    @test out == "Testing 42...$endline"

    _, out = @grab_output if Int == Int32
        @on_device @cuprintf("Testing %lld %d...\n",
                             blockIdx().x, threadIdx().x)
    elseif Sys.iswindows()
        @on_device @cuprintf("Testing %lld %lld...\n",
                             blockIdx().x, threadIdx().x)
    else
        @on_device @cuprintf("Testing %ld %ld...\n",
                             blockIdx().x, threadIdx().x)
    end
    @test out == "Testing 1 1...$endline"

    _, out = @grab_output @on_device begin
        @cuprintf("foo")
        @cuprintf("bar\n")
    end
    @test out == "foobar$endline"

    # c argument promotions
    function kernel(A)
        @cuprintf("%f %f\n", A[1], A[1])
        return
    end
    x = CuTestArray(ones(2, 2))
    _, out = @grab_output begin
        @cuda kernel(x)
        synchronize()
    end
    @test out == "1.000000 1.000000$endline"
end

end

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

# a composite type to test for more complex element types
@eval struct RGB{T}
    r::T
    g::T
    b::T
end

@testset "shared memory" begin

n = 1024

@testset "constructors" begin
    # static
    @on_device @cuStaticSharedMem(Float32, 1)
    @on_device @cuStaticSharedMem(Float32, (1,2))
    @on_device @cuStaticSharedMem(Tuple{Float32, Float32}, 1)
    @on_device @cuStaticSharedMem(Tuple{Float32, Float32}, (1,2))
    @on_device @cuStaticSharedMem(Tuple{RGB{Float32}, UInt32}, 1)
    @on_device @cuStaticSharedMem(Tuple{RGB{Float32}, UInt32}, (1,2))

    # dynamic
    @on_device @cuDynamicSharedMem(Float32, 1)
    @on_device @cuDynamicSharedMem(Float32, (1, 2))
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, 1)
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, (1,2))
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, 1)
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, (1,2))

    # dynamic with offset
    @on_device @cuDynamicSharedMem(Float32, 1, 8)
    @on_device @cuDynamicSharedMem(Float32, (1,2), 8)
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, 1, 8)
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, (1,2), 8)
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, 1, 8)
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, (1,2), 8)
end

@testset "dynamic shmem" begin

@testset "statically typed" begin
    function kernel(d, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuDynamicSharedMem(Float32, n)
        s[t] = d[t]
        sync_threads()
        d[t] = s[tr]

        return
    end

    a = rand(Float32, n)
    d_a = CuTestArray(a)

    @cuda threads=n shmem=n*sizeof(Float32) kernel(d_a, n)
    @test reverse(a) == Array(d_a)
end

@testset "parametrically typed" for T in [Int32, Int64, Float32, Float64]
    function kernel(d::CuDeviceArray{T}, n) where {T}
        t = threadIdx().x
        tr = n-t+1

        s = @cuDynamicSharedMem(T, n)
        s[t] = d[t]
        sync_threads()
        d[t] = s[tr]

        return
    end

    a = rand(T, n)
    d_a = CuTestArray(a)

    @cuda threads=n shmem=n*sizeof(T) kernel(d_a, n)
    @test reverse(a) == Array(d_a)
end

@testset "alignment" begin
    # bug: used to generate align=12, which is invalid (non pow2)
    function kernel(v0::T, n) where {T}
        shared = @cuDynamicSharedMem(T, n)
        @inbounds shared[Cuint(1)] = v0
        return
    end

    n = 32
    typ = typeof((0f0, 0f0, 0f0))
    @cuda shmem=n*sizeof(typ) kernel((0f0, 0f0, 0f0), n)
end

end

@testset "static shmem" begin

@testset "statically typed" begin
    function kernel(d, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuStaticSharedMem(Float32, 1024)
        s2 = @cuStaticSharedMem(Float32, 1024)  # catch aliasing

        s[t] = d[t]
        s2[t] = 2*d[t]
        sync_threads()
        d[t] = s[tr]

        return
    end

    a = rand(Float32, n)
    d_a = CuTestArray(a)

    @cuda threads=n kernel(d_a, n)
    @test reverse(a) == Array(d_a)
end

@testset "parametrically typed" for typ in [Int32, Int64, Float32, Float64]
    function kernel(d::CuDeviceArray{T}, n) where {T}
        t = threadIdx().x
        tr = n-t+1

        s = @cuStaticSharedMem(T, 1024)
        s2 = @cuStaticSharedMem(T, 1024)  # catch aliasing

        s[t] = d[t]
        s2[t] = d[t]
        sync_threads()
        d[t] = s[tr]

        return
    end

    a = rand(typ, n)
    d_a = CuTestArray(a)

    @cuda threads=n kernel(d_a, n)
    @test reverse(a) == Array(d_a)
end

@testset "alignment" begin
    # bug: used to generate align=12, which is invalid (non pow2)
    function kernel(v0::T) where {T}
        shared = CUDAnative.@cuStaticSharedMem(T, 32)
        @inbounds shared[Cuint(1)] = v0
        return
    end

    @cuda kernel((0f0, 0f0, 0f0))
end

end

@testset "dynamic shmem consisting of multiple arrays" begin

# common use case 1: dynamic shmem consists of multiple homogeneous arrays
#                    -> split using `view`
@testset "homogeneous" begin
    function kernel(a, b, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuDynamicSharedMem(eltype(a), 2*n)

        sa = view(s, 1:n)
        sa[t] = a[t]
        sync_threads()
        a[t] = sa[tr]

        sb = view(s, n+1:2*n)
        sb[t] = b[t]
        sync_threads()
        b[t] = sb[tr]

        return
    end

    a = rand(Float32, n)
    d_a = CuTestArray(a)

    b = rand(Float32, n)
    d_b = CuTestArray(b)

    @cuda threads=n shmem=2*n*sizeof(Float32) kernel(d_a, d_b, n)
    @test reverse(a) == Array(d_a)
    @test reverse(b) == Array(d_b)
end

# common use case 2: dynamic shmem consists of multiple heterogeneous arrays
#                    -> construct using pointer offset
@testset "heterogeneous" begin
    function kernel(a, b, n)
        t = threadIdx().x
        tr = n-t+1

        sa = @cuDynamicSharedMem(eltype(a), n)
        sa[t] = a[t]
        sync_threads()
        a[t] = sa[tr]

        sb = @cuDynamicSharedMem(eltype(b), n, n*sizeof(eltype(a)))
        sb[t] = b[t]
        sync_threads()
        b[t] = sb[tr]

        return
    end

    a = rand(Float32, n)
    d_a = CuTestArray(a)

    b = rand(Int64, n)
    d_b = CuTestArray(b)

    @cuda threads=n shmem=(n*sizeof(Float32) + n*sizeof(Int64)) kernel(d_a, d_b, n)
    @test reverse(a) == Array(d_a)
    @test reverse(b) == Array(d_b)
end

end

end

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

@testset "data movement and conversion" begin

if capability(dev) >= v"3.0"
@testset "shuffle" begin

@eval struct AddableTuple
    x::Int32
    y::Int64
    AddableTuple(val) = new(val, val*2)
end
Base.:(+)(a::AddableTuple, b::AddableTuple) = AddableTuple(a.x+b.x)

n = 14

@testset "down" for T in [Int32, Int64, Float32, Float64, AddableTuple]
    function kernel(d::CuDeviceArray{T}, n) where {T}
        t = threadIdx().x
        if t <= n
            d[t] += shfl_down(d[t], n÷2)
        end
        return
    end

    a = T[T(i) for i in 1:n]
    d_a = CuTestArray(a)

    threads = nearest_warpsize(dev, n)
    @cuda threads=threads kernel(d_a, n)

    a[1:n÷2] += a[n÷2+1:end]
    @test a == Array(d_a)
end

end
end

end

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

@testset "parallel synchronization and communication" begin

@on_device sync_threads()

@on_device sync_warp()
@on_device sync_warp(0xffffffff)

@testset "voting" begin

@testset "ballot" begin
    d_a = CuTestArray(UInt32[0])

    function kernel(a, i)
        vote = vote_ballot(threadIdx().x == i)
        if threadIdx().x == 1
            a[1] = vote
        end
        return
    end

    len = 4
    for i in 1:len
        @cuda threads=len kernel(d_a, i)
        @test Array(d_a) == [2^(i-1)]
    end
end

@testset "any" begin
    d_a = CuTestArray(UInt32[0])

    function kernel(a, i)
        vote = vote_any(threadIdx().x >= i)
        if threadIdx().x == 1
            a[1] = vote
        end
        return
    end

    @cuda threads=2 kernel(d_a, 1)
    @test Array(d_a) == [1]

    @cuda threads=2 kernel(d_a, 2)
    @test Array(d_a) == [1]

    @cuda threads=2 kernel(d_a, 3)
    @test Array(d_a) == [0]
end

@testset "all" begin
    d_a = CuTestArray(UInt32[0])

    function kernel(a, i)
        vote = vote_all(threadIdx().x >= i)
        if threadIdx().x == 1
            a[1] = vote
        end
        return
    end

    @cuda threads=2 kernel(d_a, 1)
    @test Array(d_a) == [1]

    @cuda threads=2 kernel(d_a, 2)
    @test Array(d_a) == [0]

    @cuda threads=2 kernel(d_a, 3)
    @test Array(d_a) == [0]
end

end

end

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

end


text:
# These functions return an operator for use in `get_clims(::Seres, op)`
process_clims(lims::Tuple{<:Number,<:Number}) =
    (zlims -> ifelse.(isfinite.(lims), lims, zlims)) ∘ ignorenan_extrema
process_clims(s::Union{Symbol,Nothing,Missing}) = ignorenan_extrema
# don't specialize on ::Function otherwise python functions won't work
process_clims(f) = f

get_clims(sp::Subplot)::Tuple{Float64,Float64} =
    haskey(sp.attr, :clims_calculated) ? sp[:clims_calculated] : update_clims(sp)
get_clims(series::Series)::Tuple{Float64,Float64} =
    haskey(series.plotattributes, :clims_calculated) ?
    series[:clims_calculated]::Tuple{Float64,Float64} : update_clims(series)
get_clims(sp::Subplot, series::Series)::Tuple{Float64,Float64} =
    series[:colorbar_entry] ? get_clims(sp) : get_clims(series)

function update_clims(sp::Subplot, op = process_clims(sp[:clims]))::Tuple{Float64,Float64}
    zmin, zmax = Inf, -Inf
    for series in series_list(sp)
        if series[:colorbar_entry]::Bool
            zmin, zmax = _update_clims(zmin, zmax, update_clims(series, op)...)
        else
            update_clims(series, op)
        end
    end
    return sp[:clims_calculated] = zmin <= zmax ? (zmin, zmax) : (NaN, NaN)
end

"""
    update_clims(::Series, op=Plots.ignorenan_extrema)
Finds the limits for the colorbar by taking the "z-values" for the series and passing them into `op`,
which must return the tuple `(zmin, zmax)`. The default op is the extrema of the finite
values of the input. The value is stored as a series property, which is retrieved by `get_clims`.
"""
function update_clims(series::Series, op = ignorenan_extrema)::Tuple{Float64,Float64}
    zmin, zmax = Inf, -Inf

    # keeping this unrolled has higher performance
    if series[:seriestype] ∈ _z_colored_series && series[:z] !== nothing
        zmin, zmax = update_clims(zmin, zmax, series[:z], op)
    end
    if series[:line_z] !== nothing
        zmin, zmax = update_clims(zmin, zmax, series[:line_z], op)
    end
    if series[:marker_z] !== nothing
        zmin, zmax = update_clims(zmin, zmax, series[:marker_z], op)
    end
    if series[:fill_z] !== nothing
        zmin, zmax = update_clims(zmin, zmax, series[:fill_z], op)
    end
    return series[:clims_calculated] = zmin <= zmax ? (zmin, zmax) : (NaN, NaN)
end

update_clims(zmin, zmax, vals::AbstractSurface, op)::Tuple{Float64,Float64} =
    update_clims(zmin, zmax, vals.surf, op)
update_clims(zmin, zmax, vals::Any, op)::Tuple{Float64,Float64} =
    _update_clims(zmin, zmax, op(vals)...)
update_clims(zmin, zmax, ::Nothing, ::Any)::Tuple{Float64,Float64} = zmin, zmax

_update_clims(zmin, zmax, emin, emax) = NaNMath.min(zmin, emin), NaNMath.max(zmax, emax)

@enum ColorbarStyle cbar_gradient cbar_fill cbar_lines

function colorbar_style(series::Series)
    colorbar_entry = series[:colorbar_entry]
    if !(colorbar_entry isa Bool)
        @warn "Non-boolean colorbar_entry ignored."
        colorbar_entry = true
    end

    if !colorbar_entry
        nothing
    elseif isfilledcontour(series)
        cbar_fill
    elseif iscontour(series)
        cbar_lines
    elseif series[:seriestype] ∈ (:heatmap, :surface) ||
           any(series[z] !== nothing for z in [:marker_z, :line_z, :fill_z])
        cbar_gradient
    else
        nothing
    end
end

hascolorbar(series::Series) = colorbar_style(series) !== nothing
hascolorbar(sp::Subplot) =
    sp[:colorbar] != :none && any(hascolorbar(s) for s in series_list(sp))

function get_colorbar_ticks(sp::Subplot; update = true)
    if update || !haskey(sp.attr, :colorbar_optimized_ticks)
        ticks = _transform_ticks(sp[:colorbar_ticks])
        cvals = sp[:colorbar_continuous_values]
        dvals = sp[:colorbar_discrete_values]
        clims = get_clims(sp)
        scale = sp[:colorbar_scale]
        formatter = sp[:colorbar_formatter]
        sp.attr[:colorbar_optimized_ticks] =
            get_ticks(ticks, cvals, dvals, clims, scale, formatter)
    end
    return sp.attr[:colorbar_optimized_ticks]
end

function _update_subplot_colorbars(sp::Subplot)
    # Dynamic callback from the pipeline if needed
    update_clims(sp)
end


text:
# 4x5x3 model

wet3d = falses(4,15,3)
wet3d[2, :, 1] .= true
wet3d[2, :, 2] .= true
wet3d[3:4, 7:9, 1] = true
wet3d[3:4, 12:14, 1] = true
wet3d[3, 7:9, 2] = true
wet3d[3, 12:14, 2] = true
wet3d[2, 2:4, 3] = true
wet3d[2, 7:9, 3] = true
wet3d[3, 12:14, 3] = true
wet3d[3, 7, 3] = true
wet3d[3, 9, 3] = true
wet3d[2, 12, 3] = true

const iwet = (LinearIndices(wet3d))[findall(!iszero, wet3d)] # replaces find(wet3d) :(nwet = length(iwet)


text:
# Linearized AC branch flow model for radial networks.
# Variables: squared voltage magnitude, active power, reactive power.

## Variables

""
function _PM.variable_bus_voltage(pm::BFARadPowerModel; kwargs...)
    _PM.variable_bus_voltage_magnitude_sqr(pm; kwargs...)
    _PM.variable_bus_voltage_angle(pm; kwargs...)
end

"Voltage angle of all buses is that of the reference bus"
function _PM.variable_bus_voltage_angle(pm::BFARadPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    report && _IM.sol_component_fixed(pm, _PM.pm_it_sym, nw, :bus, :va, _PM.ids(pm, nw, :bus), last(first(_PM.ref(pm,nw,:ref_buses)))["va"])
end

# Copied from _PM.variable_branch_power_real(pm::AbstractAPLossLessModels; nw::Int, bounded::Bool, report::Bool)
# Since this model is lossless, active power variables are 1 per branch instead of 2.
""
function _PM.variable_branch_power_real(pm::BFARadPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    p = _PM.var(pm, nw)[:p] = JuMP.@variable(pm.model,
        [(l,i,j) in _PM.ref(pm, nw, :arcs_from)], base_name="$(nw)_p",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branch, l), "p_start")
    )

    if bounded
        flow_lb, flow_ub = _PM.ref_calc_branch_flow_bounds(_PM.ref(pm, nw, :branch), _PM.ref(pm, nw, :bus))

        for arc in _PM.ref(pm, nw, :arcs_from)
            l,i,j = arc
            if !isinf(flow_lb[l])
                JuMP.set_lower_bound(p[arc], flow_lb[l])
            end
            if !isinf(flow_ub[l])
                JuMP.set_upper_bound(p[arc], flow_ub[l])
            end
        end
    end

    for (l,branch) in _PM.ref(pm, nw, :branch)
        if haskey(branch, "pf_start")
            f_idx = (l, branch["f_bus"], branch["t_bus"])
            JuMP.set_start_value(p[f_idx], branch["pf_start"])
        end
    end

    # this explicit type erasure is necessary
    p_expr = Dict{Any,Any}( ((l,i,j), p[(l,i,j)]) for (l,i,j) in _PM.ref(pm, nw, :arcs_from) )
    p_expr = merge(p_expr, Dict( ((l,j,i), -1.0*p[(l,i,j)]) for (l,i,j) in _PM.ref(pm, nw, :arcs_from)))
    _PM.var(pm, nw)[:p] = p_expr

    report && _IM.sol_component_value_edge(pm, _PM.pm_it_sym, nw, :branch, :pf, :pt, _PM.ref(pm, nw, :arcs_from), _PM.ref(pm, nw, :arcs_to), p_expr)
end

# Copied from _PM.variable_ne_branch_power_real(pm::AbstractAPLossLessModels; nw::Int, bounded::Bool, report::Bool)
# and improved by comparing with _PM.variable_branch_power_real(pm::AbstractAPLossLessModels; nw::Int, bounded::Bool, report::Bool).
# Since this model is lossless, active power variables are 1 per branch instead of 2.
""
function _PM.variable_ne_branch_power_real(pm::BFARadPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    p_ne = _PM.var(pm, nw)[:p_ne] = JuMP.@variable(pm.model,
        [(l,i,j) in _PM.ref(pm, nw, :ne_arcs_from)], base_name="$(nw)_p_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_branch, l), "p_start")
    )

    if bounded
        flow_lb, flow_ub = _PM.ref_calc_branch_flow_bounds(_PM.ref(pm, nw, :ne_branch), _PM.ref(pm, nw, :bus))

        for arc in _PM.ref(pm, nw, :ne_arcs_from)
            l,i,j = arc
            if !isinf(flow_lb[l])
                JuMP.set_lower_bound(p_ne[arc], flow_lb[l])
            end
            if !isinf(flow_ub[l])
                JuMP.set_upper_bound(p_ne[arc], flow_ub[l])
            end
        end
    end

    # this explicit type erasure is necessary
    p_ne_expr = Dict{Any,Any}( ((l,i,j), p_ne[(l,i,j)]) for (l,i,j) in _PM.ref(pm, nw, :ne_arcs_from) )
    p_ne_expr = merge(p_ne_expr, Dict(((l,j,i), -1.0*p_ne[(l,i,j)]) for (l,i,j) in _PM.ref(pm, nw, :ne_arcs_from)))
    _PM.var(pm, nw)[:p_ne] = p_ne_expr

    report && _IM.sol_component_value_edge(pm, _PM.pm_it_sym, nw, :ne_branch, :pf, :pt, _PM.ref(pm, nw, :ne_arcs_from), _PM.ref(pm, nw, :ne_arcs_to), p_ne_expr)
end

# Adapted from variable_branch_power_real(pm::BFARadPowerModel; nw::Int, bounded::Bool, report::Bool)
""
function _PM.variable_branch_power_imaginary(pm::BFARadPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    q = _PM.var(pm, nw)[:q] = JuMP.@variable(pm.model,
        [(l,i,j) in _PM.ref(pm, nw, :arcs_from)], base_name="$(nw)_q",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branch, l), "q_start")
    )

    if bounded
        flow_lb, flow_ub = _PM.ref_calc_branch_flow_bounds(_PM.ref(pm, nw, :branch), _PM.ref(pm, nw, :bus))

        for arc in _PM.ref(pm, nw, :arcs_from)
            l,i,j = arc
            if !isinf(flow_lb[l])
                JuMP.set_lower_bound(q[arc], flow_lb[l])
            end
            if !isinf(flow_ub[l])
                JuMP.set_upper_bound(q[arc], flow_ub[l])
            end
        end
    end

    for (l,branch) in _PM.ref(pm, nw, :branch)
        if haskey(branch, "qf_start")
            f_idx = (l, branch["f_bus"], branch["t_bus"])
            JuMP.set_start_value(q[f_idx], branch["qf_start"])
        end
    end

    # this explicit type erasure is necessary
    q_expr = Dict{Any,Any}( ((l,i,j), q[(l,i,j)]) for (l,i,j) in _PM.ref(pm, nw, :arcs_from) )
    q_expr = merge(q_expr, Dict( ((l,j,i), -1.0*q[(l,i,j)]) for (l,i,j) in _PM.ref(pm, nw, :arcs_from)))
    _PM.var(pm, nw)[:q] = q_expr

    report && _IM.sol_component_value_edge(pm, _PM.pm_it_sym, nw, :branch, :qf, :qt, _PM.ref(pm, nw, :arcs_from), _PM.ref(pm, nw, :arcs_to), q_expr)
end

# Adapted from variable_ne_branch_power_real(pm::BFARadPowerModel; nw::Int, bounded::Bool, report::Bool)
""
function _PM.variable_ne_branch_power_imaginary(pm::BFARadPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
    q_ne = _PM.var(pm, nw)[:q_ne] = JuMP.@variable(pm.model,
        [(l,i,j) in _PM.ref(pm, nw, :ne_arcs_from)], base_name="$(nw)_q_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_branch, l), "q_start")
    )

    if bounded
        flow_lb, flow_ub = _PM.ref_calc_branch_flow_bounds(_PM.ref(pm, nw, :ne_branch), _PM.ref(pm, nw, :bus))

        for arc in _PM.ref(pm, nw, :ne_arcs_from)
            l,i,j = arc
            if !isinf(flow_lb[l])
                JuMP.set_lower_bound(q_ne[arc], flow_lb[l])
            end
            if !isinf(flow_ub[l])
                JuMP.set_upper_bound(q_ne[arc], flow_ub[l])
            end
        end
    end

    # this explicit type erasure is necessary
    q_ne_expr = Dict{Any,Any}( ((l,i,j), q_ne[(l,i,j)]) for (l,i,j) in _PM.ref(pm, nw, :ne_arcs_from) )
    q_ne_expr = merge(q_ne_expr, Dict(((l,j,i), -1.0*q_ne[(l,i,j)]) for (l,i,j) in _PM.ref(pm, nw, :ne_arcs_from)))
    _PM.var(pm, nw)[:q_ne] = q_ne_expr

    report && _IM.sol_component_value_edge(pm, _PM.pm_it_sym, nw, :ne_branch, :qf, :qt, _PM.ref(pm, nw, :ne_arcs_from), _PM.ref(pm, nw, :ne_arcs_to), q_ne_expr)
end

## Constraints

"Nothing to do, this model is lossless"
function _PM.constraint_power_losses(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, tm)
end

"Nothing to do, this model is lossless"
function constraint_power_losses_on_off(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, tm, vad_min, vad_max)
end

"Nothing to do, this model is lossless"
function constraint_power_losses_frb(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, tm)
end

"Nothing to do, this model is lossless"
function constraint_power_losses_frb_on_off(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, tm, vad_min, vad_max)
end

"Nothing to do, this model is lossless"
function constraint_power_losses_oltc(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to)
end

"Nothing to do, this model is lossless"
function constraint_power_losses_oltc_on_off(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, vad_min, vad_max)
end

"Nothing to do, this model is lossless"
function constraint_ne_power_losses(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, tm, vad_min, vad_max)
end

"Nothing to do, this model is lossless"
function constraint_ne_power_losses_parallel(pm::BFARadPowerModel, n::Int, br_idx_e, br_idx_c, f_bus, t_bus, f_idx_c, t_idx_c, r_e, x_e, g_sh_fr_e, g_sh_to_e, b_sh_fr_e, b_sh_to_e, r_c, x_c, g_sh_fr_c, g_sh_to_c, b_sh_fr_c, b_sh_to_c, tm, vad_min, vad_max)
end

"Nothing to do, this model is lossless"
function constraint_ne_power_losses_frb(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, tm, vad_min, vad_max)
end

"Nothing to do, this model is lossless"
function constraint_ne_power_losses_frb_parallel(pm::BFARadPowerModel, n::Int, br_idx_e, br_idx_c, f_bus, t_bus, f_idx_c, t_idx_c, r_e, x_e, g_sh_fr_e, g_sh_to_e, b_sh_fr_e, b_sh_to_e, r_c, x_c, g_sh_fr_c, g_sh_to_c, b_sh_fr_c, b_sh_to_c, tm, vad_min, vad_max)
end

"Nothing to do, this model is lossless"
function constraint_ne_power_losses_oltc(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, vad_min, vad_max)
end

"Nothing to do, no voltage angle variables"
function _PM.constraint_voltage_angle_difference(pm::BFARadPowerModel, n::Int, f_idx, angmin, angmax)
end

"Nothing to do, no voltage angle variables"
function _PM.constraint_voltage_angle_difference_on_off(pm::BFARadPowerModel, n::Int, f_idx, angmin, angmax, vad_min, vad_max)
end

"Nothing to do, no voltage angle variables"
function _PM.constraint_ne_voltage_angle_difference(pm::BFARadPowerModel, n::Int, f_idx, angmin, angmax, vad_min, vad_max)
end

"Defines voltage drop over a branch whose `f_bus` is the reference bus"
function constraint_voltage_magnitude_difference_frb(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, b_sh_fr, tm)
    p_fr = _PM.var(pm, n, :p, f_idx)
    q_fr = _PM.var(pm, n, :q, f_idx)
    w_to = _PM.var(pm, n, :w, t_bus)
    # w_fr is assumed equal to 1.0

    JuMP.@constraint(pm.model, (1.0/tm^2) - w_to ==  2*(r*p_fr + x*q_fr))
end

"Defines voltage drop over a transformer branch that has an OLTC"
function constraint_voltage_magnitude_difference_oltc(pm::BFARadPowerModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, b_sh_fr)
    p_fr = _PM.var(pm, n, :p, f_idx)
    q_fr = _PM.var(pm, n, :q, f_idx)
    ttmi = _PM.var(pm, n, :ttmi, i)
    w_to = _PM.var(pm, n, :w, t_bus)
    # w_fr is assumed equal to 1.0 to preserve the linearity of the model

    JuMP.@constraint(pm.model, 1.0*ttmi - w_to == 2*(r*p_fr + x*q_fr))
end

"Complex power is limited by an octagon instead of a circle, so as to keep the model linear"
function _PM.constraint_thermal_limit_from(pm::BFARadPowerModel, n::Int, f_idx, rate_a)
    p_fr = _PM.var(pm, n, :p, f_idx)
    q_fr = _PM.var(pm, n, :q, f_idx)
    c_perp = cos(π/8) # ~0.92
    c_diag = sin(π/8) + cos(π/8) # == cos(π/8) * sqrt(2), ~1.31

    JuMP.@constraint(pm.model, -c_perp*rate_a <= p_fr        <= c_perp*rate_a)
    JuMP.@constraint(pm.model, -c_perp*rate_a <=        q_fr <= c_perp*rate_a)
    JuMP.@constraint(pm.model, -c_diag*rate_a <= p_fr + q_fr <= c_diag*rate_a)
    JuMP.@constraint(pm.model, -c_diag*rate_a <= p_fr - q_fr <= c_diag*rate_a)
end

"Complex power is limited by an octagon instead of a circle, so as to keep the model linear"
function _PM.constraint_thermal_limit_from_on_off(pm::BFARadPowerModel, n::Int, i, f_idx, rate_a)
    p_fr = _PM.var(pm, n, :p, f_idx)
    q_fr = _PM.var(pm, n, :q, f_idx)
    z    = _PM.var(pm, n, :z_branch, i)
    c_perp = cos(π/8) # ~0.92
    c_diag = sin(π/8) + cos(π/8) # == cos(π/8) * sqrt(2), ~1.31

    JuMP.@constraint(pm.model, p_fr        >= -c_perp*rate_a*z)
    JuMP.@constraint(pm.model, p_fr        <=  c_perp*rate_a*z)
    JuMP.@constraint(pm.model,        q_fr >= -c_perp*rate_a*z)
    JuMP.@constraint(pm.model,        q_fr <=  c_perp*rate_a*z)
    JuMP.@constraint(pm.model, p_fr + q_fr >= -c_diag*rate_a*z)
    JuMP.@constraint(pm.model, p_fr + q_fr <=  c_diag*rate_a*z)
    JuMP.@constraint(pm.model, p_fr - q_fr >= -c_diag*rate_a*z)
    JuMP.@constraint(pm.model, p_fr - q_fr <=  c_diag*rate_a*z)
end

"Complex power is limited by an octagon instead of a circle, so as to keep the model linear"
function _PM.constraint_ne_thermal_limit_from(pm::BFARadPowerModel, n::Int, i, f_idx, rate_a)
    p_fr = _PM.var(pm, n, :p_ne, f_idx)
    q_fr = _PM.var(pm, n, :q_ne, f_idx)
    z    = _PM.var(pm, n, :branch_ne, i)
    c_perp = cos(π/8) # ~0.92
    c_diag = sin(π/8) + cos(π/8) # == cos(π/8) * sqrt(2), ~1.31

    JuMP.@constraint(pm.model, p_fr        >= -c_perp*rate_a*z)
    JuMP.@constraint(pm.model, p_fr        <=  c_perp*rate_a*z)
    JuMP.@constraint(pm.model,        q_fr >= -c_perp*rate_a*z)
    JuMP.@constraint(pm.model,        q_fr <=  c_perp*rate_a*z)
    JuMP.@constraint(pm.model, p_fr + q_fr >= -c_diag*rate_a*z)
    JuMP.@constraint(pm.model, p_fr + q_fr <=  c_diag*rate_a*z)
    JuMP.@constraint(pm.model, p_fr - q_fr >= -c_diag*rate_a*z)
    JuMP.@constraint(pm.model, p_fr - q_fr <=  c_diag*rate_a*z)
end

""
function constraint_ne_thermal_limit_from_parallel(pm::BFARadPowerModel, n::Int, br_idx_e, br_idx_c, f_idx_c, rate_a_e, rate_a_c)
    # Suffixes: _e: existing branch; _c: candidate branch; _p: parallel equivalent
    branch_e = _PM.ref(pm, n, :branch, br_idx_e)
    branch_c = _PM.ref(pm, n, :ne_branch, br_idx_c)
    r_e      = branch_e["br_r"]
    r_c      = branch_c["br_r"]
    x_e      = branch_e["br_x"]
    x_c      = branch_c["br_x"]

    r_p      = (r_e*(r_c^2+x_c^2)+r_c*(r_e^2+x_e^2)) / ((r_e+r_c)^2+(x_e+x_c)^2)
    x_p      = (x_e*(r_c^2+x_c^2)+x_c*(r_e^2+x_e^2)) / ((r_e+r_c)^2+(x_e+x_c)^2)
    rate_a_p = min(rate_a_e*sqrt(r_e^2+x_e^2),rate_a_c*sqrt(r_c^2+x_c^2)) / sqrt(r_p^2+x_p^2)

    _PM.constraint_ne_thermal_limit_from(pm, n, br_idx_c, f_idx_c, rate_a_p)
end

"Nothing to do, this model is symmetric"
function _PM.constraint_thermal_limit_to(pm::BFARadPowerModel, n::Int, t_idx, rate_a)
end

"Nothing to do, this model is symmetric"
function _PM.constraint_thermal_limit_to_on_off(pm::BFARadPowerModel, n::Int, i, t_idx, rate_a)
end

"Nothing to do, this model is symmetric"
function _PM.constraint_ne_thermal_limit_to(pm::BFARadPowerModel, n::Int, i, t_idx, rate_a)
end

"Nothing to do, this model is symmetric"
function constraint_ne_thermal_limit_to_parallel(pm::BFARadPowerModel, n::Int, br_idx_e, br_idx_c, f_idx_c, rate_a_e, rate_a_c)
end

## Other functions

"""
Converts the solution data into the data model's standard space, polar voltages and rectangular power.

Bus voltage magnitude `vm` is the square root of `w`.
Voltage magnitude of the reference bus is 1.0 p.u.
Branch OLTC tap ratio `tm` (if applies) is the square root of the inverse of `ttmi`.
"""
function _PM.sol_data_model!(pm::BFARadPowerModel, solution::Dict)
    if haskey(solution["it"]["pm"], "nw")
        nws_sol = solution["it"]["pm"]["nw"]
    else
        nws_sol = Dict("0" => solution)
    end

    for (nw, nw_sol) in nws_sol

        # Find reference bus id
        n = parse(Int, nw)
        ref_buses = _PM.ref(pm, n, :ref_buses)
        if length(ref_buses) == 0
            Memento.error(_LOGGER, "no reference bus found")
        end
        if length(ref_buses) > 1
            Memento.error(_LOGGER, "networks with multiple reference buses are not supported")
        end
        ref_bus_id = first(keys(ref_buses))

        # Bus voltage magnitude `vm` is the square root of `w`
        for (i,bus) in nw_sol["bus"]
            if haskey(bus, "w")
                bus["vm"] = sqrt(bus["w"])
                delete!(bus, "w")
            end
        end

        # The voltage magnitude of the reference bus is 1.0 p.u.
        nw_sol["bus"]["$ref_bus_id"]["vm"] = 1.0

        # OLTC tap ratio `tm` of `branch`es (if applies) is the square root of the inverse of `ttmi`
        for (i,br) in nw_sol["branch"]
            if haskey(br, "ttmi")
                if haskey(nw_sol,"ne_branch") && any([nw_sol["ne_branch"]["$ne_br_id"]["built"] for ne_br_id in ne_branch_ids(pm, parse(Int,i), nw = n)] .== 1.0)
                    # if branch is not built
                    br["tm"] = 0.0
                else
                    br["tm"] = sqrt(1.0/br["ttmi"])
                end
                delete!(br, "ttmi")
            end
        end

        # OLTC tap ratio `tm` of `ne_branch`es (if applies) is the square root of the inverse of `ttmi`
        if haskey(nw_sol,"ne_branch")
            for (i,br) in nw_sol["ne_branch"]
                if haskey(br, "ttmi")
                    if br["built"] == 0.0
                        br["tm"] = 0.0
                    else
                        br["tm"] = sqrt(1.0/br["ttmi"])
                    end
                    delete!(br, "ttmi")
                end
            end
        end

    end
end


text:
include("wannier90/projections.jl")

import Base: ==
"""
Represents an element.
"""
struct Element
    symbol        ::Symbol
    Z             ::Int64
    name          ::String
    atomic_weight ::Float64
    color         ::NTuple{3, Float64}
end

"""
Reads all the elements from the file.
"""
const ELEMENTS = Element[]
open(joinpath(@__DIR__, "..", "assets", "elements.txt"), "r") do f
    while !eof(f)
        line = split(readline(f))
        push!(ELEMENTS, Element(Symbol(line[4]), Meta.parse(line[1]), line[9], Meta.parse(line[10]), (Meta.parse.(line[5:7],)...,)./65535))
    end
end
#undefined element
push!(ELEMENTS, Element(:undef, 0, "undef",0, (0.0, 0.0, 0.0)))
function element(sym::Symbol)
    if tryparse(Int, String(sym)[end:end]) != nothing
        sym = Symbol(String(sym)[1:end-1])
    end
    found = filter(x->x.symbol == sym, ELEMENTS)
    if isempty(found)
       @warn "No element found with symbol '$sym'."
       return ELEMENTS[end]
    end
    return found[1]
end

element(z::Int) = getfirst(x->x.Z == z, ELEMENTS)

abstract type AbstractAtom{T, LT<:Length{T}} end

#Definition of the AbstractAtom interface, each AbstractAtom needs to provide a method `atom(...)` that returns a datastructure with the Atom fields.
elsym(atom::AbstractAtom) = element(atom).symbol
"Extracts all the positions of the atoms and puts them in a vector."
positions(atoms::Vector{<:AbstractAtom}, name::Symbol) = [position_cart(x) for x in filter(y -> name(y) == name, atoms)]
getpseudoset(at::AbstractAtom) = getpseudoset(elsym(at), pseudo(at))

ismagnetic(at::AbstractAtom) = !iszero(sum(magnetization(at)))

"Takes a Vector of atoms and returns a Vector with the atoms having unique symbols."
function Base.unique(atoms::Vector{<:AbstractAtom{T}}) where T <: AbstractFloat
    uni = AbstractAtom{T}[]
    for at in atoms
        if findfirst(x -> isequal_species(x, at), uni) === nothing
            push!(uni, at)
        end
    end
    return uni
end

@with_kw mutable struct DFTU{T}
	l::Int = -1
	U::T   = zero(T)
	J0::T  = zero(T)
	#QE params
	α::T   = zero(T)
	β::T   = zero(T)
	J::Vector{T} = T[zero(T)]
end

function ==(x::DFTU, y::DFTU)
	fnames = fieldnames(DFTU)
	for fn in fnames
		if getfield(x, fn) != getfield(y, fn)
			return false
		end
	end
	return true
end

isdefault(x::DFTU{T}) where {T} =
	x == DFTU{T}()

isdefault(x::Any) = isempty(x)
isdefault(x::AbstractVector) = isempty(x) || all(iszero, x)

Base.string(::Type{Elk}, dftu::DFTU) = "$(dftu.l) $(dftu.U) $(dftu.J0)"

mutable struct Pseudo
	name::String
	dir ::String
	Pseudo() =
		new("", "")
	Pseudo(name::AbstractString, dir::AbstractString) =
		new(name, abspath(dir))
end

Base.isempty(p::Pseudo) =
	isempty(p.name) && isempty(p.dir)

==(p1::Pseudo, p2::Pseudo) =
	p1.name == p2.name && p1.dir == p2.dir

path(p::Pseudo) = joinpath(p.dir, p.name)
# TODO Multiple l per atom in Elk??
#We use angstrom everywhere
@with_kw_noshow mutable struct Atom{T<:AbstractFloat, LT<:Length{T}} <: AbstractAtom{T, LT}
    name          ::Symbol
    element       ::Element
    position_cart ::Point3{LT}
    position_cryst::Point3{T}=zero(Point3{T})
    pseudo        ::Pseudo=Pseudo()
    projections   ::Vector{Projection} = Projection[]
    magnetization ::Vec3{T} = zero(Vec3{T})
    dftu          ::DFTU{T} = DFTU{T}()
end

Atom(name::Symbol, el::Element, pos_cart::Point3{LT}, pos_cryst::Point3{T}; kwargs...) where {T, LT<:Length{T}} =
	Atom{T, LT}(name=name, element=el, position_cart=pos_cart, position_cryst=pos_cryst; kwargs...)
Atom(name::Symbol, el::Symbol, args...; kwargs...) =
	Atom(name, element(el), args...; kwargs...)

#TODO this is a little iffy
Atom(orig_at::Atom, new_pos_cart::Point3, new_pos_cryst::Point3) = Atom(name(orig_at), element(orig_at), new_pos_cart, new_pos_cryst, pseudo(orig_at), projections(orig_at), magnetization(orig_at), dftu(orig_at))
#Easiest way to implement a new abstractatom is to provide a way to access
#the struct holding `name`, `position_cart`, `element`, `pseudo`, `projection` fields
atom(at::Atom) = at

for interface_function in fieldnames(Atom)
	@eval $interface_function(at::AbstractAtom) = atom(at).$interface_function
	@eval export $interface_function
end

function Base.range(at::AbstractAtom)
	projs = projections(at)
	@assert length(projs) != 0 "At $(name(at)) has no defined projections. Please use `setprojections!` first."
	return projs[1].start : projs[end].last
end

setname!(at::AbstractAtom, name::Symbol) =
	atom(at).name = name

length_unit(at::Atom{T, LT}) where {T, LT} = LT

function setpseudo!(at::AbstractAtom, pseudo::Pseudo; print=true)
	print && @info "Pseudo of atom $(name(at)) set to $pseudo."
	!ispath(path(pseudo)) && @warn "Pseudopath $(path(pseudo)) not found."
	atom(at).pseudo = pseudo
end

function setprojections!(at::AbstractAtom, projections::Vector{Projection}; print=true)
    print && @info "Setting projections for atom $(name(at)) to $projections"
    atom(at).projections = projections
end

bondlength(at1::AbstractAtom{T}, at2::AbstractAtom{T}, R=T(0.0)) where T<:AbstractFloat = norm(position_cart(at1) - position_cart(at2) - R)

==(at1::AbstractAtom{T, LT}, at2::AbstractAtom{T, LT}) where {T, LT} =
    name(at1) == name(at2) && norm(position_cart(at1) - position_cart(at2)) < LT(1e-6)

function isequal_species(at1::AbstractAtom, at2::AbstractAtom)
    for f in fieldnames(typeof(at1))
        if f in (:position_cart, :position_cryst, :projections)
            continue
        end
        if getfield(at1, f) != getfield(at2, f)
            return false
        end
    end
    return true
end

#TODO fix documentation
for hub_param in (:U, :J0, :α, :β)
	f = Symbol("set_Hubbard_$(hub_param)!")
	str = "$(hub_param)"
	@eval begin
		"""
			$($(f))(at::AbstractAtom, v::AbstractFloat; print=true)

		Set the Hubbard $($(str)) parameter for the specified atom.

		Example:
			`$($(f))(at, 2.1)'
		"""
		function $f(at::AbstractAtom{T}, v::AbstractFloat; print=true) where {T}
			dftu(at).$(hub_param) = convert(T, v)
			print && @info "Hubbard $($(str)) of atom $(at.name) set to $v"
		end
	end
end

"""
	set_Hubbard_J!(at::AbstractAtom, v::Vector{<:AbstractFloat}; print=true)

Set the Hubbard J parameter for the specified atom.

Example:
	`set_Hubbard_J(at, [2.1])'
"""
function set_Hubbard_J!(at::AbstractAtom{T}, v::Vector{<:AbstractFloat}; print=true) where {T}
	dftu(at).J = convert.(T, v)
	print && @info "Hubbard J of atom $(at.name) set to $v"
end

function position_string(::Type{QE}, at::AbstractAtom; relative=true)
	pos = relative ? round.(position_cryst(at), digits=5) : round.(ustrip.(uconvert.(Ang, position_cart(at))), digits=5)
	return "$(name(at))  $(pos[1]) $(pos[2]) $(pos[3])\n"
end

function set_magnetization!(at::AbstractAtom, mag; print=true)
	at.magnetization = convert(Vec3, mag)
	print && @info "Magnetization of at $(name(at)) was set to $(magnetization(at))"
end

"""
	distance(at1::AbstractAtom, at2::AbstractAtom)

Calculates the distance between the two atoms.
"""
distance(at1::AbstractAtom, at2::AbstractAtom) = norm(at1.position_cart - at2.position_cart)

"""
	set_position!(at::AbstractAtom, pos::AbstractVector{T}, unit_cell::Mat3) where {T<:Real}

Updates the position of the atom to this. The unit cell is used to make sure both `position_cryst` and `position_cart` are correct.
"""
function set_position!(at::AbstractAtom, pos::AbstractVector{T}, unit_cell::Mat3) where {T<:Real}
	atom(at).position_cryst = Point3{T}(pos...)
	atom(at).position_cart  = unit_cell * atom(at).position_cryst
	return at
end

function set_position!(at::AbstractAtom, pos::AbstractVector{T}, unit_cell::Mat3) where {T<:Length}
	atom(at).position_cart  = Point3{T}(pos...)
	atom(at).position_cryst = unit_cell^-1 * atom(at).position_cart
	return at
end

"""
	scale_bondlength!(at1::AbstractAtom, at2::AbstractAtom, scale::Real, cell::Mat3)

Scales the bondlength between two atoms. The center of mass remains the same.
"""
function scale_bondlength!(at1::AbstractAtom, at2::AbstractAtom, scale::Real, cell::Mat3)
	p1 = position_cryst(at1)
	p2 = position_cryst(at2)
	mid = (p1 + p2)/2
	orig_dist = norm(p1 - p2)
	direction = (p1 - mid)/norm(p1-mid)
	new_dist = orig_dist * scale
	new_p1 = mid + new_dist * direction/2
	new_p2 = mid - new_dist * direction/2
	set_position!(at1, new_p1, cell)
	set_position!(at2, new_p2, cell)
end


text:
using Test
import KenyaCoV
include("test_functions.jl");
@testset "Simulation tests" begin
    @test test_no_infecteds()
    @test stays_nonnegative()
end


text:
function Base.cat(a::NamedDimsArray{L}; dims) where {L}
    newL = expand_dimnames(L, dims)
    numerical_dims = dim(newL, dims)
    data = Base.cat(parent(a); dims=numerical_dims) # Base.cat is type unstable
    return NamedDimsArray{newL}(data)
end

# to dispatch on the first _or the second_ argument being the NDA.
for (T, S) in [
    (:NamedDimsArray, :AbstractArray),
    (:AbstractArray, :NamedDimsArray),
    (:NamedDimsArray, :NamedDimsArray),
]
    @eval function Base.cat(a::$T, b::$S, cs::AbstractArray...; dims)
        combL = unify_names_longest(dimnames(a), dimnames(b), dimnames.(cs)...)
        newL = expand_dimnames(combL, dims)
        numerical_dims = dim(newL, dims)
        data = Base.cat(unname(a), unname(b), unname.(cs)...; dims=numerical_dims)
        return NamedDimsArray{newL}(data)
    end
end

function Base.hcat(a::NamedDimsArray{L}) where {L}
    newL = expand_dimnames(L, 2)
    data = Base.hcat(parent(a))
    return NamedDimsArray{newL}(data)
end

Base.vcat(a::NamedDimsArray{L}) where {L} = a

for (T, S) in [
    (:NamedDimsVecOrMat, :NamedDimsVecOrMat),
    (:NamedDimsVecOrMat, :AbstractVecOrMat),
    (:AbstractVecOrMat, :NamedDimsVecOrMat),
]
    for (fun, d) in zip((:vcat, :hcat), (1, 2))
        @eval function Base.$fun(a::$T, b::$S, cs::AbstractVecOrMat...)
            combL = unify_names_longest(dimnames(a), dimnames(b), dimnames.(cs)...)
            newL = expand_dimnames(combL, $d)
            data = Base.$fun(unname(a), unname(b), unname.(cs)...)
            return NamedDimsArray{newL}(data)
        end
    end
end

for (T, S) in [
    (:NamedDimsVector, :NamedDimsVector),
    (:NamedDimsVector, :AbstractVector),
    (:AbstractVector, :NamedDimsVector),
]
    @eval function Base.vcat(a::$T, b::$S, cs::AbstractVector...)
        newL = unify_names(dimnames(a), dimnames(b), dimnames.(cs)...)
        data = vcat(unname(a), unname(b), unname.(cs)...)
        return NamedDimsArray{newL}(data)
    end
end

for (f, nf, tf, tup) in
    [(:vcat, :_named_vcat, :_typed_vcat, ()), (:hcat, :_named_hcat, :_typed_hcat, (:_,))]
    @eval begin
        function Base.reduce(::typeof($f), A::AbstractVector{<:NamedDimsVecOrMat})
            return $nf(mapreduce(dimnames, unify_names_longest, A), A)
        end
        function Base.reduce(::typeof($f), A::NamedDimsVector{<:Any,<:AbstractVecOrMat})
            return $nf(mapreduce(dimnames, unify_names_longest, A), A)
        end
        function Base.reduce(::typeof($f), A::NamedDimsVector{<:Any,<:NamedDimsVecOrMat})
            return $nf(mapreduce(dimnames, unify_names_longest, A), A)
        end

        function $nf(Linner, A)
            Louter = ($tup..., dimnames(A)...)
            Lnew = unify_names_longest(Linner, Louter)
            data = Base.$tf(mapreduce(eltype, promote_type, A), A) # same as Base
            return NamedDimsArray{Lnew}(data)
        end
    end
end


text:
# Arvind Prasadan
# 2019 January
# DMF Package
# Main Module File

"""
Module providing the Dynamic Mode Factorication (cf. Dynamic Mode Decomposition) algorithm.
"""
module DMF

# Packages
# Math
using LinearAlgebra
# Signal Processing
using DSP
# Statistics
using Statistics
using QuantEcon
using StatsBase
using Combinatorics
# SOBI Via R: Dropped
# using RCall

# Functions to export (all)
export 	dmf, # Main Algorithm
	optshrink, # Denoising
	permute_IP, # Best inner product
	eigenvector_error, # Loss function
	autocorrelation, # Autocorrelation function
	autocorr_mat, # Autocorrelation matrix
	gen_arma_sequence, # Generate a realization of an ARMA process
	gen_cos_sequence # Generate a cosine sequence
	# SOBI_Wrapper, # Wrapper for SOBI method

# Code files
include("dmf_function.jl")
include("optshrink.jl")
# include("SOBI_Wrapper.jl")
include("permute_IP.jl")
include("eigenvector_error.jl")
include("autocorrelation.jl")
include("autocorr_mat.jl")
include("gen_arma_sequence.jl")
include("gen_cos_sequence.jl")
end


text:
using DifferentialEquations
using Sundials
using BenchmarkTools
using Plots

using DifferentialEquations, Plots
function orego(du,u,p,t)
  s,q,w = p
  y1,y2,y3 = u
  du[1] = s*(y2+y1*(1-q*y1-y2))
  du[2] = (y3-(1+y1)*y2)/s
  du[3] = w*(y1-y3)
end
p = [77.27,8.375e-6,0.161]
prob = ODEProblem(orego,[1.0,2.0,3.0],(0.0,360.0),p)
sol = solve(prob)
plot(sol)

plot(sol,vars=(1,2,3))

using BenchmarkTools
prob = ODEProblem(orego,[1.0,2.0,3.0],(0.0,50.0),p)
@btime sol = solve(prob,Tsit5())

@btime sol = solve(prob,Rodas5())

function orego(du,u,p,t)
  s,q,w = p
  y1,y2,y3 = u
  du[1] = s*(y2+y1*(1-q*y1-y2))
  du[2] = (y3-(1+y1)*y2)/s
  du[3] = w*(y1-y3)
end
function g(du,u,p,t)
  du[1] = 0.1u[1]
  du[2] = 0.1u[2]
  du[3] = 0.1u[3]
end
p = [77.27,8.375e-6,0.161]
prob = SDEProblem(orego,g,[1.0,2.0,3.0],(0.0,30.0),p)
sol = solve(prob,SOSRI())
plot(sol)

sol = solve(prob,ImplicitRKMil()); plot(sol)

sol = solve(prob,ImplicitRKMil()); plot(sol)

function orego(du,u,p,t)
  s,q,w = p
  y1,y2,y3 = u
  du[1] = s*(y2+y1*(1-q*y1-y2))
  du[2] = (y3-(1+y1)*y2)/s
  du[3] = w*(y1-y3)
end
p = [60.0,1e-5,0.2]
prob = ODEProblem(orego,[1.0,2.0,3.0],(0.0,30.0),p)
sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)

function onecompartment(du,u,p,t)
  Ka,Ke = p
  du[1] = -Ka*u[1]
  du[2] =  Ka*u[1] - Ke*u[2]
end
p = (Ka=2.268,Ke=0.07398)
prob = ODEProblem(onecompartment,[100.0,0.0],(0.0,90.0),p)

tstops = [24,48,72]
condition(u,t,integrator) = t ∈ tstops
affect!(integrator) = (integrator.u[1] += 100)
cb = DiscreteCallback(condition,affect!)
sol = solve(prob,Tsit5(),callback=cb,tstops=tstops)
plot(sol)

function onecompartment_delay(du,u,h,p,t)
  Ka,Ke,τ = p
  delayed_depot = h(p,t-τ)[1]
  du[1] = -Ka*u[1]
  du[2] =  Ka*delayed_depot - Ke*u[2]
end
p = (Ka=2.268,Ke=0.07398,τ=6.0)
h(p,t) = [0.0,0.0]
prob = DDEProblem(onecompartment_delay,[100.0,0.0],h,(0.0,90.0),p)

tstops = [24,48,72]
condition(u,t,integrator) = t ∈ tstops
affect!(integrator) = (integrator.u[1] += 100)
cb = DiscreteCallback(condition,affect!)
sol = solve(prob,MethodOfSteps(Rosenbrock23()),callback=cb,tstops=tstops)
plot(sol)

p = (Ka = 0.5, Ke = 0.1, τ = 4.0)

function f(du, u, p, t)
    du[1] = -p[1]*u[1] + p[2]*u[2]*u[3]
    du[2] = p[1]*u[1] - p[2]*u[2]*u[3] - p[3]*u[2]*u[2]
    du[3] = u[1] + u[2] + u[3] - 1.
end
M = [1 0 0; 0 1 0; 0 0 0.]
p = [0.04, 10^4, 3e7]
u0 = [1.,0.,0.]
tspan = (0., 1e6)
prob = ODEProblem(ODEFunction(f, mass_matrix = M), u0, tspan, p)
sol = solve(prob, Rodas5())
plot(sol, xscale=:log10, tspan=(1e-6, 1e5), layout=(3,1))

# Robertson Equation DAE Implicit form
function h(out, du, u, p, t)
    out[1] = -p[1]*u[1] + p[2]*u[2]*u[3] - du[1]
    out[2] = p[1]*u[1] - p[2]*u[2]*u[3] - p[3]*u[2]*u[2] - du[2]
    out[3] = u[1] + u[2] + u[3] - 1.
end
p = [0.04, 10^4, 3e7]
du0 = [-0.04, 0.04, 0.0]
u0 = [1.,0.,0.]
tspan = (0., 1e6)
differential_vars = [true, true, false]
prob = DAEProblem(h, du0, u0, tspan, p, differential_vars = differential_vars)
sol = solve(prob, IDA())
plot(sol, xscale=:log10, tspan=(1e-6, 1e5), layout=(3,1))

function f(out, da, a, p, t)
   (L, m, g) = p
   u, v, x, y, T = a
   du, dv, dx, dy, dT = da
   out[1] = x*T/(m*L) - du
   out[2] = y*T/(m*L) - g - dv
   out[3] = u - dx
   out[4] = v - dy
   out[5] = u^2 + v^2 - y*g + T/m
   nothing
end

# Release pendulum from top right
u0 = zeros(5)
u0[3] = 1.0
du0 = zeros(5)
du0[2] = 9.81

p = [1,1,9.8]
tspan = (0.,100.)

differential_vars = [true, true, true, true, false]
prob = DAEProblem(f, du0, u0, tspan, p, differential_vars = differential_vars)
sol = solve(prob, IDA())
plot(sol, vars=(3,4))

function f(out, da, a, p, t)
   L1, m1, L2, m2, g = p

   u1, v1, x1, y1, T1,
   u2, v2, x2, y2, T2 = a

   du1, dv1, dx1, dy1, dT1,
   du2, dv2, dx2, dy2, dT2 = da

   out[1]  = x2*T2/(m2*L2) - du2
   out[2]  = y2*T2/(m2*L2) - g - dv2
   out[3]  = u2 - dx2
   out[4]  = v2 - dy2
   out[5]  = u2^2 + v2^2 -y2*g + T2/m2

   out[6]  = x1*T1/(m1*L1) - x2*T2/(m2*L2) - du1
   out[7]  = y1*T1/(m1*L1) - g - y2*T2/(m2*L2) - dv1
   out[8]  = u1 - dx1
   out[9]  = v1 - dy1
   out[10] = u1^2 + v1^2 + T1/m1 +
                (-x1*x2 - y1*y2)/(m1*L2)*T2 - y1*g
   nothing
end

# Release pendulum from top right
u0 = zeros(10)
u0[3] = 1.0
u0[8] = 1.0
du0 = zeros(10)
du0[2] = 9.8
du0[7] = 9.8

p = [1,1,1,1,9.8]
tspan = (0.,100.)

differential_vars = [true, true, true, true, false,
                     true, true, true, true, false]
prob = DAEProblem(f, du0, u0, tspan, p, differential_vars = differential_vars)
sol = solve(prob, IDA())

plot(sol, vars=(3,4))
plot(sol, vars=(8,9))

using DifferentialEquations, Sundials, Plots

# initial condition
function init_brusselator_2d(xyd)
    N = length(xyd)
    u = zeros(N, N, 2)
    for I in CartesianIndices((N, N))
        x = xyd[I[1]]
        y = xyd[I[2]]
        u[I,1] = 22*(y*(1-y))^(3/2)
        u[I,2] = 27*(x*(1-x))^(3/2)
    end
    u
end

N = 32

xyd_brusselator = range(0,stop=1,length=N)

u0 = vec(init_brusselator_2d(xyd_brusselator))

tspan = (0, 22.)

p = (3.4, 1., 10., xyd_brusselator)

brusselator_f(x, y, t) = ifelse((((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) &&
                                (t >= 1.1), 5., 0.)

using LinearAlgebra, SparseArrays
du = ones(N-1)
D2 = spdiagm(-1 => du, 0=>fill(-2.0, N), 1 => du)
D2[1, N] = D2[N, 1] = 1
D2 = 1/step(xyd_brusselator)^2*D2
tmp = Matrix{Float64}(undef, N, N)
function brusselator_2d_op(du, u, (D2, tmp, p), t)
    A, B, α, xyd = p
    dx = step(xyd)
    N = length(xyd)
    α = α/dx^2
    du = reshape(du, N, N, 2)
    u = reshape(u, N, N, 2)
    @views for i in axes(u, 3)
        ui = u[:, :, i]
        dui = du[:, :, i]
        mul!(tmp, D2, ui)
        mul!(dui, ui, D2')
        dui .+= tmp
    end

    @inbounds begin
        for I in CartesianIndices((N, N))
            x = xyd[I[1]]
            y = xyd[I[2]]
            i = I[1]
            j = I[2]

            du[i,j,1] = α*du[i,j,1] + B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
            du[i,j,2] = α*du[i,j,2] + A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
        end
    end
    nothing
end

prob1 = ODEProblem(brusselator_2d_op, u0, tspan, (D2, tmp, p))

sol1 = @time solve(prob1, TRBDF2(autodiff=false));

@gif for t in sol1.t[1]:0.1:sol1.t[end]
    off = N^2
    solt = sol1(t)
    plt1 = surface(reshape(solt[1:off], N, N), zlims=(0, 5), leg=false)
    surface!(plt1, reshape(solt[off+1:end], N, N), zlims=(0, 5), leg=false)
    display(plt1)
end

function brusselator_2d_loop(du, u, p, t)
    A, B, α, xyd = p
    dx = step(xyd)
    N = length(xyd)
    α = α/dx^2
    limit = a -> let N=N
        a == N+1 ? 1 :
        a == 0 ? N :
        a
    end
    II = LinearIndices((N, N, 2))

    @inbounds begin
        for I in CartesianIndices((N, N))
            x = xyd[I[1]]
            y = xyd[I[2]]
            i = I[1]
            j = I[2]
            ip1 = limit(i+1)
            im1 = limit(i-1)
            jp1 = limit(j+1)
            jm1 = limit(j-1)

            ii1 = II[i,j,1]
            ii2 = II[i,j,2]

            du[II[i,j,1]] = α*(u[II[im1,j,1]] + u[II[ip1,j,1]] + u[II[i,jp1,1]] + u[II[i,jm1,1]] - 4u[ii1]) +
            B + u[ii1]^2*u[ii2] - (A + 1)*u[ii1] + brusselator_f(x, y, t)

            du[II[i,j,2]] = α*(u[II[im1,j,2]] + u[II[ip1,j,2]] + u[II[i,jp1,2]] + u[II[i,jm1,2]] - 4u[II[i,j,2]]) +
            A*u[ii1] - u[ii1]^2*u[ii2]
        end
    end
    nothing
end

prob2 = ODEProblem(brusselator_2d_loop, u0, tspan, p)

sol2 = @time solve(prob2, TRBDF2())
sol2_2 = @time solve(prob2, CVODE_BDF())

using SparseDiffTools, SparsityDetection

sparsity_pattern = jacobian_sparsity(brusselator_2d_loop,similar(u0),u0,p,2.0)
jac_sp = sparse(sparsity_pattern)
jac = Float64.(jac_sp)
colors = matrix_colors(jac)
prob3 = ODEProblem(ODEFunction(brusselator_2d_loop, colorvec=colors,jac_prototype=jac_sp), u0, tspan, p)
sol3 = @time solve(prob3, TRBDF2())

using DiffEqOperators
using Sundials
using AlgebraicMultigrid: ruge_stuben, aspreconditioner, smoothed_aggregation
prob6 = ODEProblem(ODEFunction(brusselator_2d_loop, jac_prototype=JacVecOperator{Float64}(brusselator_2d_loop, u0)), u0, tspan, p)
II = Matrix{Float64}(I, N, N)
Op = kron(Matrix{Float64}(I, 2, 2), kron(D2, II) + kron(II, D2))
Wapprox = -I+Op
#ml = ruge_stuben(Wapprox)
ml = smoothed_aggregation(Wapprox)
precond = aspreconditioner(ml)
sol_trbdf2 = @time solve(prob6, TRBDF2(linsolve=LinSolveGMRES())); # no preconditioner
sol_trbdf2 = @time solve(prob6, TRBDF2(linsolve=LinSolveGMRES(Pl=lu(Wapprox)))); # sparse LU
sol_trbdf2 = @time solve(prob6, TRBDF2(linsolve=LinSolveGMRES(Pl=precond))); # AMG
sol_cvodebdf = @time solve(prob2, CVODE_BDF(linear_solver=:GMRES));

function laplacian2d(du, u, p, t)
    A, B, α, xyd = p
    dx = step(xyd)
    N = length(xyd)
    du = reshape(du, N, N, 2)
    u = reshape(u, N, N, 2)
    @inbounds begin
        α = α/dx^2
        limit = a -> let N=N
            a == N+1 ? 1 :
            a == 0 ? N :
            a
        end
        for I in CartesianIndices((N, N))
            x = xyd[I[1]]
            y = xyd[I[2]]
            i = I[1]
            j = I[2]
            ip1 = limit(i+1)
            im1 = limit(i-1)
            jp1 = limit(j+1)
            jm1 = limit(j-1)
            du[i,j,1] = α*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1])
            du[i,j,2] = α*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2])
        end
    end
    nothing
end
function brusselator_reaction(du, u, p, t)
    A, B, α, xyd = p
    dx = step(xyd)
    N = length(xyd)
    du = reshape(du, N, N, 2)
    u = reshape(u, N, N, 2)
    @inbounds begin
        for I in CartesianIndices((N, N))
            x = xyd[I[1]]
            y = xyd[I[2]]
            i = I[1]
            j = I[2]
            du[i,j,1] = B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
            du[i,j,2] = A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
        end
    end
    nothing
end
prob7 = SplitODEProblem(laplacian2d, brusselator_reaction, u0, tspan, p)
sol7 = @time solve(prob7, KenCarp4())
M = MatrixFreeOperator((du,u,p)->laplacian2d(du, u, p, 0), (p,), size=(2*N^2, 2*N^2), opnorm=1000)
prob7_2 = SplitODEProblem(M, brusselator_reaction, u0, tspan, p)
sol7_2 = @time solve(prob7_2, ETDRK4(krylov=true), dt=1)
prob7_3 = SplitODEProblem(DiffEqArrayOperator(Op), brusselator_reaction, u0, tspan, p)
sol7_3 = solve(prob7_3, KenCarp4());

using DiffEqDevTools
abstols = 0.1 .^ (5:8)
reltols = 0.1 .^ (1:4)
sol = solve(prob3,CVODE_BDF(linear_solver=:GMRES),abstol=1/10^7,reltol=1/10^10)
test_sol = TestSolution(sol)
probs = [prob2, prob3, prob6]
setups = [Dict(:alg=>CVODE_BDF(),:prob_choice => 1),
          Dict(:alg=>CVODE_BDF(linear_solver=:GMRES), :prob_choice => 1),
          Dict(:alg=>TRBDF2(), :prob_choice => 1),
          Dict(:alg=>TRBDF2(linsolve=LinSolveGMRES(Pl=precond)), :prob_choice => 3),
          Dict(:alg=>TRBDF2(), :prob_choice => 2)
         ]
labels = ["CVODE_BDF (dense)" "CVODE_BDF (GMRES)" "TRBDF2 (dense)" "TRBDF2 (sparse)" "TRBDF2 (GMRES)"]
wp = WorkPrecisionSet(probs,abstols,reltols,setups;appxsol=[test_sol,test_sol,test_sol],save_everystep=false,numruns=3,
  names=labels, print_names=true, seconds=0.5)
plot(wp)

function henon(dz,z,p,t)
  p₁, p₂, q₁, q₂ = z[1], z[2], z[3], z[4]
  dp₁ = -q₁*(1 + 2q₂)
  dp₂ = -q₂-(q₁^2 - q₂^2)
  dq₁ = p₁
  dq₂ = p₂

  dz .= [dp₁, dp₂, dq₁, dq₂]
  return nothing
end

u₀ = [0.1, 0.0, 0.0, 0.5]
prob = ODEProblem(henon, u₀, (0., 1000.))
sol = solve(prob, Vern9(), abstol=1e-14, reltol=1e-14)

plot(sol, vars=[(3,4,1)], tspan=(0,100))

function henon(ddz,dz,z,p,t)
  p₁, p₂ = dz[1], dz[2]
  q₁, q₂ = z[1], z[2]
  ddq₁ = -q₁*(1 + 2q₂)
  ddq₂ = -q₂-(q₁^2 - q₂^2)

  ddz .= [ddq₁, ddq₂]
end

p₀ = u₀[1:2]
q₀ = u₀[3:4]
prob2 = SecondOrderODEProblem(henon, p₀, q₀, (0., 1000.))
sol = solve(prob2, DPRKN6(), abstol=1e-10, reltol=1e-10)

plot(sol, vars=[(3,4)], tspan=(0,100))

H(p, q, params) = 1/2 * (p[1]^2 + p[2]^2) + 1/2 * (q[1]^2 + q[2]^2 + 2q[1]^2 * q[2] - 2/3*q[2]^3)

prob3 = HamiltonianProblem(H, p₀, q₀, (0., 1000.))
sol = solve(prob3, DPRKN6(), abstol=1e-10, reltol=1e-10)

plot(sol, vars=[(3,4)], tspan=(0,100))

function generate_ics(E,n)
  # The hardcoded values bellow can be estimated by looking at the
  # figures in the Henon-Heiles 1964 article
  qrange = range(-0.4, stop = 1.0, length = n)
  prange = range(-0.5, stop = 0.5, length = n)
  z0 = Vector{Vector{typeof(E)}}()
  for q in qrange
    V = H([0,0],[0,q],nothing)
    V ≥ E && continue
    for p in prange
      T = 1/2*p^2
      T + V ≥ E && continue
      z = [√(2(E-V-T)), p, 0, q]
      push!(z0, z)
    end
  end
  return z0
end

z0 = generate_ics(0.125, 10)

function prob_func(prob,i,repeat)
  @. prob.u0 = z0[i]
  prob
end

ensprob = EnsembleProblem(prob, prob_func=prob_func)
sim = solve(ensprob, Vern9(), EnsembleThreads(), trajectories=length(z0))

plot(sim, vars=(3,4), tspan=(0,10))

using DiffEqGPU

function henon_gpu(dz,z,p,t)
  @inbounds begin
    dz[1] = -z[3]*(1 + 2z[4])
    dz[2] = -z[4]-(z[3]^2 - z[4]^2)
    dz[3] = z[1]
    dz[4] = z[2]
  end
  return nothing
end

z0 = generate_ics(0.125f0, 50)
prob_gpu = ODEProblem(henon_gpu, Float32.(u₀), (0.f0, 1000.f0))
ensprob = EnsembleProblem(prob_gpu, prob_func=prob_func)
sim = solve(ensprob, Tsit5(), EnsembleGPUArray(), trajectories=length(z0))

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
module LassoPath

       using COSMO, JuMP, LinearAlgebra
       using Statistics
       import Plots as Plots
       import ProgressMeter as Pm
       MOI = JuMP.MathOptInterface

       # Write your package code here.

       # The order here matters
       include("utils.jl")
       include("lasso_root.jl")
       include("model_make.jl")
       include("lasso_scop.jl")
       include("prox.jl")
       include("lasso_prox.jl")
       include("lasso_methods.jl")

       export vanderMonde,
              PrintTitle,
              Warn

       export LassoRoot,
              LassoSCOP,
              Lasso2NormProx,
              GetLassoPath,
              VisualizeLassoPath,
              CaptureImportantWeights

       # For testing internal metods
       # export BuildProximalGradientLasso,
       #        OptimizeProximalGradient

end