text:
# Reduction per Harris #3
function reduce_3(
    accum,
    input,
    num_samples,
)
    # define needed incides
    threads_per_block = blockDim().x
    block_idx = blockIdx().x
    thread_idx = threadIdx().x
    sample_idx = (block_idx - 1) * threads_per_block + thread_idx

    # allocate the shared memory for the partial sum
    shmem = @cuDynamicSharedMem(Float32, threads_per_block)

    # each thread loads one element from global to shared memory
    @inbounds if sample_idx <= num_samples
        shmem[thread_idx] = input[sample_idx]
    end

    # wait until all finished
    sync_threads()

    # do (partial) reduction in shared memory
    s::UInt32 = threads_per_block ÷ 2
    @inbounds while s != 0
        sync_threads()
        if thread_idx - 1 < s
            shmem[thread_idx] += shmem[thread_idx + s]
        end

        s ÷= 2
    end

    # first thread returns the result of reduction to global memory
    @inbounds if thread_idx == 1
        accum[blockIdx().x] = shmem[1]
    end

    return nothing
end

# Complex reduction per Harris #3
function reduce_cplx_3(
    accum_re,
    accum_im,
    input_re,
    input_im,
    num_samples,
)
    # define needed incides
    threads_per_block = iq_offset = blockDim().x
    block_idx = blockIdx().x
    thread_idx = threadIdx().x
    sample_idx = (block_idx - 1) * threads_per_block + thread_idx

    # allocate the shared memory for the partial sum
    # double the memory for complex values, accessed via
    # iq_offset
    shmem = @cuDynamicSharedMem(Float32, (2 * threads_per_block))

    # each thread loads one element from global to shared memory
    @inbounds if sample_idx <= num_samples
        shmem[thread_idx + 0 * iq_offset] = input_re[sample_idx]
        shmem[thread_idx + 1 * iq_offset] = input_im[sample_idx]
    end

    # wait until all finished
    sync_threads()

    # do (partial) reduction in shared memory
    s::UInt32 = threads_per_block ÷ 2
    @inbounds while s != 0
        sync_threads()
        if thread_idx - 1 < s
            shmem[thread_idx + 0 * iq_offset] += shmem[thread_idx + 0 * iq_offset + s]
            shmem[thread_idx + 1 * iq_offset] += shmem[thread_idx + 1 * iq_offset + s]
        end

        s ÷= 2
    end

    # first thread returns the result of reduction to global memory
    @inbounds if thread_idx == 1
        accum_re[blockIdx().x] = shmem[1 + 0 * iq_offset]
        accum_im[blockIdx().x] = shmem[1 + 1 * iq_offset]
    end

    return nothing
end

# Complex reduction per Harris #3, multi correlator, multi antenna
function reduce_cplx_multi_3(
    accum_re,
    accum_im,
    input_re,
    input_im,
    num_samples,
    num_ants::NumAnts{NANT},
    correlator_sample_shifts::SVector{NCOR, Int64}
) where {NANT, NCOR}
    # define needed incides
    threads_per_block = iq_offset = blockDim().x
    block_idx = blockIdx().x
    thread_idx = threadIdx().x
    sample_idx = (block_idx - 1) * threads_per_block + thread_idx

    # allocate the shared memory for the partial sum
    # double the memory for complex values, accessed via
    # iq_offset
    shmem = @cuDynamicSharedMem(Float32, (2 * threads_per_block, NANT, NCOR))
    @inbounds for antenna_idx = 1:NANT
        for corr_idx = 1:NCOR
            shmem[thread_idx + 0 * iq_offset, antenna_idx, corr_idx]
            shmem[thread_idx + 1 * iq_offset, antenna_idx, corr_idx]
        end
    end

    # each thread loads one element from global to shared memory
    @inbounds if sample_idx <= num_samples
        for antenna_idx = 1:NANT
            for corr_idx = 1:NCOR
                shmem[thread_idx + 0 * iq_offset, antenna_idx, corr_idx] = input_re[sample_idx, antenna_idx, corr_idx]
            	shmem[thread_idx + 1 * iq_offset, antenna_idx, corr_idx] = input_im[sample_idx, antenna_idx, corr_idx]
            end
        end
    end

    # wait until all finished
    sync_threads()

    # do (partial) reduction in shared memory
    s::UInt32 = threads_per_block ÷ 2
    @inbounds while s != 0
        sync_threads()
        if thread_idx - 1 < s
            for antenna_idx = 1:NANT
                for corr_idx = 1:NCOR
                    shmem[thread_idx + 0 * iq_offset, antenna_idx, corr_idx] += shmem[thread_idx + 0 * iq_offset + s, antenna_idx, corr_idx]
                    shmem[thread_idx + 1 * iq_offset, antenna_idx, corr_idx] += shmem[thread_idx + 1 * iq_offset + s, antenna_idx, corr_idx]
                end
            end
        end

        s ÷= 2
    end

    # first thread returns the result of reduction to global memory
    @inbounds if thread_idx == 1
        for antenna_idx = 1:NANT
            for corr_idx = 1:NCOR
                accum_re[blockIdx().x, antenna_idx, corr_idx] = shmem[1 + 0 * iq_offset, antenna_idx, corr_idx]
                accum_im[blockIdx().x, antenna_idx, corr_idx] = shmem[1 + 1 * iq_offset, antenna_idx, corr_idx]
            end
        end
    end

    return nothing
end

# Complex reduction per Harris #3, multi correlator, multi antenna
function reduce_cplx_multi_31(
    accum_re,
    accum_im,
    input_re,
    input_im,
    num_samples,
    num_ants::NumAnts{NANT},
    correlator_sample_shifts::SVector{NCOR, Int64}
) where {NANT, NCOR}
    # define needed incides
    threads_per_block = iq_offset = blockDim().x
    block_idx = blockIdx().x
    thread_idx = threadIdx().x
    sample_idx = (block_idx - 1) * threads_per_block + thread_idx

    # allocate the shared memory for the partial sum
    # double the memory for complex values, accessed via
    # iq_offset
    shmem = @cuDynamicSharedMem(Float32, (2 * threads_per_block, NANT, NCOR))
    @inbounds for antenna_idx = 1:NANT
        for corr_idx = 1:NCOR

            shmem[thread_idx + 0 * iq_offset, antenna_idx, corr_idx]
            shmem[thread_idx + 1 * iq_offset, antenna_idx, corr_idx]

            # each thread loads one element from global to shared memory
            @inbounds if sample_idx <= num_samples
                shmem[thread_idx + 0 * iq_offset, antenna_idx, corr_idx] = input_re[sample_idx, antenna_idx, corr_idx]
            	shmem[thread_idx + 1 * iq_offset, antenna_idx, corr_idx] = input_im[sample_idx, antenna_idx, corr_idx]
            end

            # wait until all finished
            sync_threads()

            # do (partial) reduction in shared memory
            s::UInt32 = threads_per_block ÷ 2
            @inbounds while s != 0
                sync_threads()
                if thread_idx - 1 < s
                    shmem[thread_idx + 0 * iq_offset, antenna_idx, corr_idx] += shmem[thread_idx + 0 * iq_offset + s, antenna_idx, corr_idx]
                    shmem[thread_idx + 1 * iq_offset, antenna_idx, corr_idx] += shmem[thread_idx + 1 * iq_offset + s, antenna_idx, corr_idx]
                end
                s ÷= 2
            end

            # first thread returns the result of reduction to global memory
            @inbounds if thread_idx == 1
                accum_re[blockIdx().x, antenna_idx, corr_idx] = shmem[1 + 0 * iq_offset, antenna_idx, corr_idx]
                accum_im[blockIdx().x, antenna_idx, corr_idx] = shmem[1 + 1 * iq_offset, antenna_idx, corr_idx]
            end
        end
    end

    return nothing
end

# Reduction per Harris #4
function reduce_4(
    accum,
    input,
    num_samples
)
    ## define needed incides
    # local thread index
    tid = threadIdx().x
    # global thread index
    i = (blockIdx().x - 1) * (blockDim().x * 2) + threadIdx().x

    # allocate the shared memory for the partial sum
    shmem = @cuDynamicSharedMem(Float32, blockDim().x)

    # local sum variable
    mysum = 0.0f0

    # each thread loads one element from global to shared memory
    # AND
    # does the first level of reduction
    @inbounds if i <= length(input)
        mysum = input[i]
        if i + blockDim().x <= length(input)
            mysum += input[i + blockDim().x]
        end

        shmem[tid] = mysum
    end

    # wait until all finished
    sync_threads()

    # do (partial) reduction in shared memory
    s::UInt32 = blockDim().x ÷ 2
    @inbounds while s != 0
        sync_threads()
        if tid - 1 < s
            shmem[tid] += shmem[tid + s]
        end

        s ÷= 2
    end

    # first thread returns the result of reduction to global memory
    @inbounds if tid == 1
        accum[blockIdx().x] = shmem[1]
    end

    return nothing
end

# Complex reduction per Harris #4
function reduce_cplx_4(
    accum_re,
    accum_im,
    input_re,
    input_im
)
    ## define needed incides
    # local thread index
    tid = threadIdx().x
    # global thread index
    i = (blockIdx().x - 1) * (blockDim().x * 2) + threadIdx().x

    # allocate the shared memory for the partial sum
    shmem = @cuDynamicSharedMem(Float32, 2 * blockDim().x)
    # wipe values
    shmem[tid + 0 * blockDim().x] = 0.0f0
    shmem[tid + 1 * blockDim().x] = 0.0f0

    # local sum variable
    mysum_re = mysum_im = 0.0f0

    # each thread loads one element from global to shared memory
    # AND
    # does the first level of reduction
    @inbounds if i <= length(input_re)
        shmem[tid + 0 * blockDim().x] = input_re[i]
        shmem[tid + 1 * blockDim().x] = input_im[i]

        if i + blockDim().x <= length(input_re)
            shmem[tid + 0 * blockDim().x] += input_re[i + blockDim().x]
            shmem[tid + 1 * blockDim().x] += input_im[i + blockDim().x]
        end
    end

    # wait until all finished
    sync_threads()

    # do (partial) reduction in shared memory
    s::UInt32 = blockDim().x ÷ 2
    @inbounds while s != 0
        sync_threads()
        if tid - 1 < s
            shmem[tid + 0 * blockDim().x] += shmem[tid + s + 0 * blockDim().x]
            shmem[tid + 1 * blockDim().x] += shmem[tid + s + 1 * blockDim().x]
        end

        s ÷= 2
    end

    # first thread returns the result of reduction to global memory
    @inbounds if tid == 1
        accum_re[blockIdx().x] = shmem[1 + 0 * blockDim().x]
        accum_im[blockIdx().x] = shmem[1 + 1 * blockDim().x]
    end

    return nothing
end

# Complex reduction per Harris #4, multi correlator, multi antenna
function reduce_cplx_multi_4(
    accum_re,
    accum_im,
    input_re,
    input_im,
    num_samples,
    num_ants::NumAnts{NANT},
    correlator_sample_shifts::SVector{NCOR, Int64}
) where {NANT, NCOR}
    ## define needed incides
    # local thread index
    tid = threadIdx().x
    # global thread index
    i = (blockIdx().x - 1) * (blockDim().x * 2) + threadIdx().x

    # allocate the shared memory for the partial sum
    shmem = @cuDynamicSharedMem(Float32, (2 * blockDim().x, NANT, NCOR))
    # wipe values
    @inbounds for antenna_idx = 1:NANT
        for corr_idx = 1:NCOR
            shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] = 0.0f0
            shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] = 0.0f0
        end
    end

    # each thread loads one element from global to shared memory
    # AND
    # does the first level of reduction
    @inbounds if i <= num_samples
        for antenna_idx = 1:NANT
            for corr_idx = 1:NCOR
                shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] = input_re[i, antenna_idx, corr_idx]
                shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] = input_im[i, antenna_idx, corr_idx]

                if i + blockDim().x <= num_samples
                    shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += input_re[i + blockDim().x, antenna_idx, corr_idx]
                    shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += input_im[i + blockDim().x, antenna_idx, corr_idx]
                end
            end
        end
    end

    # wait until all finished
    sync_threads()

    # do (partial) reduction in shared memory
    s::UInt32 = blockDim().x ÷ 2
    @inbounds while s != 0
        sync_threads()
        if tid - 1 < s
            for antenna_idx = 1:NANT
                for corr_idx = 1:NCOR
                    shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + s + 0 * blockDim().x, antenna_idx, corr_idx]
                    shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + s + 1 * blockDim().x, antenna_idx, corr_idx]
                end
            end
        end

        s ÷= 2
    end

    # first thread returns the result of reduction to global memory
    @inbounds if tid == 1
        for antenna_idx = 1:NANT
            for corr_idx = 1:NCOR
                accum_re[blockIdx().x, antenna_idx, corr_idx] = shmem[1 + 0 * blockDim().x, antenna_idx, corr_idx]
                accum_im[blockIdx().x, antenna_idx, corr_idx] = shmem[1 + 1 * blockDim().x, antenna_idx, corr_idx]
            end
        end
    end

    return nothing
end

@inline function warp_reduce(shmem, tid)
    shmem[tid] += shmem[tid + 32]
    shmem[tid] += shmem[tid + 16]
    shmem[tid] += shmem[tid + 8]
    shmem[tid] += shmem[tid + 4]
    shmem[tid] += shmem[tid + 2]
    shmem[tid] += shmem[tid + 1]
end

function reduce_5(accum, input, num_samples)
    # define indices
    tid = threadIdx().x
    idx = (2 * blockDim().x) * (blockIdx().x - 1) + threadIdx().x

    # define shared memory
    shmem = CuDynamicSharedArray(Float32, blockDim().x)

    # each thread loads one element into the shared memory
    # and performs the first level of reduction
    if idx <= num_samples
        @inbounds shmem[tid] = input[idx]
        if idx + blockDim().x <= num_samples
            @inbounds shmem[tid] += input[idx + blockDim().x]
        end
    end

    # wait untill all threads have finished
    sync_threads()

    # do tree-like (partial) reduction in shared memory
    s = blockDim().x ÷ 2
    while s > 32
        if tid - 1 < s
            @inbounds shmem[tid] += shmem[tid + s]
        end
        sync_threads()
        s ÷= 2
    end

    # once three size = warp size, do warp reduction
    # assumes block size >= 64
    if tid <= 32
        warp_reduce(shmem, tid)
    end

    # first thread holds the result of the (partial) reduction
    if tid == 1
        @inbounds accum[blockIdx().x] = shmem[1]
    end

    return nothing
end

@inline function warp_reduce_cplx(shmem, tid)
    shmem[tid + 0 * blockDim().x] += shmem[tid + 0 * blockDim().x + 32]
    shmem[tid + 0 * blockDim().x] += shmem[tid + 0 * blockDim().x + 16]
    shmem[tid + 0 * blockDim().x] += shmem[tid + 0 * blockDim().x + 8]
    shmem[tid + 0 * blockDim().x] += shmem[tid + 0 * blockDim().x + 4]
    shmem[tid + 0 * blockDim().x] += shmem[tid + 0 * blockDim().x + 2]
    shmem[tid + 0 * blockDim().x] += shmem[tid + 0 * blockDim().x + 1]
    shmem[tid + 1 * blockDim().x] += shmem[tid + 1 * blockDim().x + 32]
    shmem[tid + 1 * blockDim().x] += shmem[tid + 1 * blockDim().x + 16]
    shmem[tid + 1 * blockDim().x] += shmem[tid + 1 * blockDim().x + 8]
    shmem[tid + 1 * blockDim().x] += shmem[tid + 1 * blockDim().x + 4]
    shmem[tid + 1 * blockDim().x] += shmem[tid + 1 * blockDim().x + 2]
    shmem[tid + 1 * blockDim().x] += shmem[tid + 1 * blockDim().x + 1]
end

# Complex reduction per Harris #5
function reduce_cplx_5(
    accum_re,
    accum_im,
    input_re,
    input_im,
    num_samples
)
    ## define needed incides
    # local thread index
    tid = threadIdx().x
    # global thread index
    idx = (2 * blockDim().x) * (blockIdx().x - 1) + threadIdx().x

    # define shared memory
    shmem = CuDynamicSharedArray(Float32, 2 * blockDim().x)

    # each thread loads one element into the shared memory
    # and performs the first level of reduction
    if idx <= num_samples
        @inbounds shmem[tid + 0 * blockDim().x] = input_re[idx]
        @inbounds shmem[tid + 1 * blockDim().x] = input_im[idx]

        if idx + blockDim().x <= num_samples
            @inbounds shmem[tid + 0 * blockDim().x] += input_re[idx + blockDim().x]
            @inbounds shmem[tid + 1 * blockDim().x] += input_im[idx + blockDim().x]
        end
    end

    # wait untill all threads have finished
    sync_threads()

    # do tree-like (partial) reduction in shared memory
    s::UInt32 = blockDim().x ÷ 2
    while s > 32
        sync_threads()
        if tid - 1 < s
            @inbounds shmem[tid + 0 * blockDim().x] += shmem[tid + s + 0 * blockDim().x]
            @inbounds shmem[tid + 1 * blockDim().x] += shmem[tid + s + 1 * blockDim().x]
        end

        s ÷= 2
    end

    # once three size = warp size, do warp reduction
    # assumes block size >= 64
    if tid <= 32
        warp_reduce_cplx(shmem, tid)
    end

    # first thread holds the result of the (partial) reduction
    if tid == 1
        @inbounds accum_re[blockIdx().x] = shmem[1 + 0 * blockDim().x]
        @inbounds accum_im[blockIdx().x] = shmem[1 + 1 * blockDim().x]
    end

    return nothing
end

@inline function warp_reduce_cplx_multi(shmem, tid, antenna_idx, corr_idx)
    shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 0 * blockDim().x + 32, antenna_idx, corr_idx]
    shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 0 * blockDim().x + 16, antenna_idx, corr_idx]
    shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 0 * blockDim().x + 8, antenna_idx, corr_idx]
    shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 0 * blockDim().x + 4, antenna_idx, corr_idx]
    shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 0 * blockDim().x + 2, antenna_idx, corr_idx]
    shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 0 * blockDim().x + 1, antenna_idx, corr_idx]
    shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 1 * blockDim().x + 32, antenna_idx, corr_idx]
    shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 1 * blockDim().x + 16, antenna_idx, corr_idx]
    shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 1 * blockDim().x + 8, antenna_idx, corr_idx]
    shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 1 * blockDim().x + 4, antenna_idx, corr_idx]
    shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 1 * blockDim().x + 2, antenna_idx, corr_idx]
    shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + 1 * blockDim().x + 1, antenna_idx, corr_idx]
end

# Complex reduction per Harris #5, multi correlator, multi antenna
function reduce_cplx_multi_5(
    accum_re,
    accum_im,
    input_re,
    input_im,
    num_samples,
    num_ants::NumAnts{NANT},
    correlator_sample_shifts::SVector{NCOR, Int64}
) where {NANT, NCOR}
    ## define needed incides
    # local thread index
    tid = threadIdx().x
    # global thread index
    idx = (2 * blockDim().x) * (blockIdx().x - 1) + threadIdx().x

    # allocate the shared memory for the partial sum
    shmem = CuDynamicSharedArray(Float32, 2 * blockDim().x)

    # each thread loads one element from global to shared memory
    # AND
    # does the first level of reduction
    if idx <= num_samples
        for antenna_idx = 1:NANT
            for corr_idx = 1:NCOR
                @inbounds shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] = input_re[idx, antenna_idx, corr_idx]
                @inbounds shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] = input_im[idx, antenna_idx, corr_idx]

                if idx + blockDim().x <= num_samples
                    @inbounds shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += input_re[idx + blockDim().x, antenna_idx, corr_idx]
                    @inbounds shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += input_im[idx + blockDim().x, antenna_idx, corr_idx]
                end
            end
        end
    end

    # wait until all finished
    sync_threads()

    # do (partial) reduction in shared memory
    s::UInt32 = blockDim().x ÷ 2
    while s > 32
        sync_threads()
        if tid - 1 < s
            for antenna_idx = 1:NANT
                for corr_idx = 1:NCOR
                    @inbounds shmem[tid + 0 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + s + 0 * blockDim().x, antenna_idx, corr_idx]
                    @inbounds shmem[tid + 1 * blockDim().x, antenna_idx, corr_idx] += shmem[tid + s + 1 * blockDim().x, antenna_idx, corr_idx]
                end
            end
        end

        s ÷= 2
    end

    # once three size = warp size, do warp reduction
    # assumes block size >= 64
    if tid <= 32
        for antenna_idx = 1:NANT
            for corr_idx = 1:NCOR
                warp_reduce_cplx_multi(shmem, tid, antenna_idx, corr_idx)
                sync_threads()
            end
        end
    end

    # first thread returns the result of reduction to global memory
    if tid == 1
        for antenna_idx = 1:NANT
            for corr_idx = 1:NCOR
                @inbounds accum_re[blockIdx().x, antenna_idx, corr_idx] = shmem[1 + 0 * blockDim().x, antenna_idx, corr_idx]
                @inbounds accum_im[blockIdx().x, antenna_idx, corr_idx] = shmem[1 + 1 * blockDim().x, antenna_idx, corr_idx]
            end
        end
    end

    return nothing
end

# Complex reduction per Harris #5, multi correlator, multi antenna
function reduce_cplx_multi_nant_5(
    accum_re,
    accum_im,
    input_re,
    input_im,
    num_samples,
    num_ants::NumAnts{NANT},
) where {NANT}
    ## define needed incides
    # local thread index
    tid = threadIdx().x
    # local antenna index
    antenna_idx = blockIdx().y
    # local sample index
    sample_idx = (blockIdx().x - 1) * (2 * blockDim().x) + threadIdx().x;

    # allocate the shared memory for the partial sum
    shmem = CuDynamicSharedArray(Float32, 2 * blockDim().x)

    # each thread loads one element from global to shared memory
    if sample_idx <= num_samples
        @inbounds shmem[tid + 0 * blockDim().x] = input_re[sample_idx, antenna_idx]
        @inbounds shmem[tid + 1 * blockDim().x] = input_im[sample_idx, antenna_idx]
        if sample_idx + blockDim().x <= num_samples
            @inbounds shmem[tid + 0 * blockDim().x] += input_re[sample_idx + blockDim().x, antenna_idx]
            @inbounds shmem[tid + 1 * blockDim().x] += input_im[sample_idx + blockDim().x, antenna_idx]
        end
    end

    # wait until all finished
    sync_threads()

    # do (partial) reduction in shared memory
    s::UInt32 = blockDim().x ÷ 2
    while s > 32
        sync_threads()
        if tid - 1 < s
            @inbounds shmem[tid + 0 * blockDim().x] += shmem[tid + s + 0 * blockDim().x]
            @inbounds shmem[tid + 1 * blockDim().x] += shmem[tid + s + 1 * blockDim().x]
        end

        s ÷= 2
    end

    if tid <= 32
        warp_reduce_cplx(shmem, tid)
    end

    # first thread returns the result of reduction to global memory
    if tid == 1
        @inbounds accum_re[blockIdx().x, antenna_idx] = shmem[1 + 0 * blockDim().x]
        @inbounds accum_im[blockIdx().x, antenna_idx] = shmem[1 + 1 * blockDim().x]
    end

    return nothing
end

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
#!/bin/env julia

# adapted from https://github.com/JuliaEditorSupport/julia-vim/blob/master/autoload/generate_latex_symbols_table.jl

# needs postprocessing due to composing chars (removing whitespace around symbols)

using REPL

open("julia.json","w") do f
    println(f, "[")
    println(f, join([string("  [ \" ", unicode, " \", \"\\", latex, "\"],\n") for (latex,unicode) in sort!(vcat(collect(REPL.REPLCompletions.latex_symbols),collect(REPL.REPLCompletions.emoji_symbols)), by=x->x[2])]))
    println(f, "]")
end

ext: jl


text:
function NevilleInterpolator()
    return NevilleInterpolator(())
end

function interpolate(obj::NevilleInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "interpolate", PolynomialFunctionLagrangeForm, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end


text:
struct Quintic{T} <: SPHKernel
    n_neighbours::Int64
    norm_1D::T
    norm_2D::T
    norm_3D::T
end

"""
    Quintic(T::DataType=Float64, n_neighbours::Integer=64)

Set up a `Quintic` kernel for a given DataType `T`.
"""
Quintic(T::DataType=Float64, n_neighbours::Integer=216) = Quintic{T}(n_neighbours, 243.0/40.0, 15309.0/(478.0*π), 2187.0/(40.0*pi))

"""
    kernel_value_1D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

Evaluate quintic spline at position ``u = \\frac{x}{h}``.
"""
function kernel_value_1D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

    n = kernel.norm_1D * h_inv

    if u < 1/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        u_m13 = T(1/3) - u
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 -
                6u_m23*u_m23*u_m23*u_m23*u_m23 +
                15u_m13*u_m13*u_m13*u_m13*u_m13 ) * n |> T
    elseif u < 2/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 -
                6u_m23*u_m23*u_m23*u_m23*u_m23 ) * n |> T
    elseif u < 1
        u_m1  = 1 - u
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_1D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

Evaluate the derivative of the Quintic spline at position ``u = \\frac{x}{h}``.
"""
function kernel_deriv_1D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

    n = kernel.norm_1D * h_inv^2

    if u < 1/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        u_m13 = T(1/3) - u
        return ( -5u_m1*u_m1*u_m1*u_m1 +
                 30u_m23*u_m23*u_m23*u_m23  -
                 75u_m13*u_m13*u_m13*u_m13 ) * n |> T
    elseif u < 2/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        return ( -5u_m1*u_m1*u_m1*u_m1 +
                 30u_m23*u_m23*u_m23*u_m23 - 75 ) * n |> T
    elseif u < 1
        u_m1 =  1 - u
        return -5u_m1*u_m1*u_m1*u_m1 * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_1D(kernel::Quintic{T}, density::Real, m::Real, h_inv::Real) where T

Does not do anything for the BSplines. Implemented for stability.
"""
function bias_correction_1D(kernel::Quintic{T}, density::Real, m::Real, h_inv::Real) where T
    return density
end

"""
    kernel_value_2D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

Evaluate quintic spline at position ``u = \\frac{x}{h}``.
"""
function kernel_value_2D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

    n = kernel.norm_2D * h_inv^2

    if u < 1/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        u_m13 = T(1/3) - u
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 -
                6u_m23*u_m23*u_m23*u_m23*u_m23 +
               15u_m13*u_m13*u_m13*u_m13*u_m13 ) * n |> T
    elseif u < 2/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 -
                6u_m23*u_m23*u_m23*u_m23*u_m23 ) * n |> T
    elseif u < 1
        u_m1  = 1 - u
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_2D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

Evaluate the derivative of the Quintic spline at position ``u = \\frac{x}{h}``.
"""
function kernel_deriv_2D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

    n = kernel.norm_2D * h_inv^3

    if u < 1/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        u_m13 = T(1/3) - u
        return ( -5u_m1*u_m1*u_m1*u_m1 +
                 30u_m23*u_m23*u_m23*u_m23 -
                 75u_m13*u_m13*u_m13*u_m13 ) * n |> T
    elseif u < 2/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        return ( -5u_m1*u_m1*u_m1*u_m1 +
                 30u_m23*u_m23*u_m23*u_m23 - 75 ) * n |> T
    elseif u < 1
        u_m1  = 1 - u
        return -5u_m1*u_m1*u_m1*u_m1 * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_2D(kernel::Quintic{T}, density::Real, m::Real, h_inv::Real) where T

Does not do anything for the BSplines. Implemented for stability.
"""
function bias_correction_2D(kernel::Quintic{T}, density::Real, m::Real, h_inv::Real) where T
    return density
end

"""
    kernel_value_3D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

Evaluate quintic spline at position ``u = \\frac{x}{h}``.
"""
function kernel_value_3D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

    n = kernel.norm_3D * h_inv^3

    if u < 1/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        u_m13 = T(1/3) - u
        return (  u_m1*u_m1*u_m1*u_m1*u_m1 -
                 6u_m23*u_m23*u_m23*u_m23*u_m23 +
                15u_m13*u_m13*u_m13*u_m13*u_m13 ) * n |> T
    elseif u < 2/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 -
                6u_m23*u_m23*u_m23*u_m23*u_m23 ) * n |> T
    elseif u < 1
        u_m1 = 1 - u
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_3D(kernel::Quintic, u::Real, h_inv::Real)

Evaluate the derivative of the Quintic spline at position ``u = \\frac{x}{h}``.
"""
function kernel_deriv_3D(kernel::Quintic{T}, u::Real, h_inv::Real) where T

    n = kernel.norm_3D * h_inv^4

    if u < 1/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        u_m13 = T(1/3) - u
        return ( -5u_m1*u_m1*u_m1*u_m1 +
                30u_m23*u_m23*u_m23*u_m23 -
                75u_m13*u_m13*u_m13*u_m13 ) * n |> T
    elseif u < 2/3
        u_m1  = 1 - u
        u_m23 = T(2/3) - u
        return ( -5u_m1*u_m1*u_m1*u_m1 +
                 30u_m23*u_m23*u_m23*u_m23 - 75 ) * n |> T
    elseif u < 1
        u_m1  = 1 - u
        return ( -5u_m1*u_m1*u_m1*u_m1 ) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_3D(kernel::Quintic{T}, density::Real, m::Real, h_inv::Real) where T

Does not do anything for the BSplines. Implemented for stability.
"""
function bias_correction_3D(kernel::Quintic{T}, density::Real, m::Real, h_inv::Real) where T
    return density |> T
end


text:
module PkgTest

include("test-fit.jl")
include("test-blb.jl")
# include("test-memory-profiling.jl")

end


text:
# Code to get PiecewiseDeterministicMarkovProcesses working with ForwardDiff
# Makes use of the trick in http://docs.juliadiffeq.org/latest/basics/faq.html#Are-the-native-Julia-solvers-compatible-with-autodifferentiation?-1
using ForwardDiff

struct DiffCache{T<:AbstractArray, S<:AbstractArray}
	rate::T
	dual_rate::S
end

function DiffCache(u::AbstractArray{T}, siz, ::Type{Val{chunk_size}}) where {T, chunk_size}
	DiffCache(u, zeros(ForwardDiff.Dual{nothing,T,chunk_size}, siz...))
end

dualcache(u::AbstractArray, N=Val{ForwardDiff.pickchunksize(length(u))}) = DiffCache(u, size(u), N)

get_rate(dc::DiffCache, u::AbstractArray{T}) where T<:ForwardDiff.Dual = reinterpret(T, dc.dual_rate)
get_rate(dc::DiffCache, u::AbstractArray) = dc.rate


text:
#InspectDR: Base functionnality and types
#-------------------------------------------------------------------------------

#AnnotationArea
#TODO: Themes: Define StyleInfo/style?
#TODO: Themes: Vector/Dict of StyleInfo/Layout?

#==Abstracts
===============================================================================#
abstract Plot
abstract PlotCanvas

#==Constants
===============================================================================#
const DInf = convert(DReal, Inf)
const DNaN = convert(DReal, NaN)

const COLOR_TRANSPARENT = ARGB32(0, 0, 0, 0)

const COLOR_BLACK = RGB24(0, 0, 0)
const COLOR_WHITE = RGB24(1, 1, 1)

const COLOR_RED = RGB24(1, 0, 0)
const COLOR_GREEN = RGB24(0, 1, 0)
const COLOR_BLUE = RGB24(0, 0, 1)

#==Plot-level structures
===============================================================================#
#=Input structures are high-level constructs describing the plot.  They will
be broken down to lower-level structures later on.
=#

type LineAttributes <: AttributeList
	style
	width #[0, 10]
	color
end
line(;style=:solid, width=1, color=COLOR_BLACK) =
	LineAttributes(style, width, color)

#"line" constructor:
#TODO: Inadequate
#eval(genexpr_attriblistbuilder(:line, LineAttributes, reqfieldcnt=0))

type GlyphAttributes <: AttributeList #Don't use "Symbol" - name used by Julia
#==IMPORTANT:
Edge width & color taken from LineAttributes
==#
	shape #because "type" is reserved
	size #of glyph.  edge width taken from LineAttributes
	color #glyph linecolor. = nothing to match line.color
	fillcolor
end
glyph(;shape=:none, size=3, color=nothing, fillcolor=nothing) =
	GlyphAttributes(shape, size, color, fillcolor)

#"glyph" constructor:
#TODO: Inadequate
#eval(genexpr_attriblistbuilder(:glyph, GlyphAttributes, reqfieldcnt=0))

type GridAttributes <: AttributeList
	#Bool values
	vmajor
	vminor
	hmajor
	hminor
end
grid(;vmajor=false, vminor=false, hmajor=false, hminor=false) =
	GridAttributes(vmajor, vminor, hmajor, hminor)

#Dispatchable scale:
immutable AxisScale{T}
	call{T}(t::Type{AxisScale{T}}) = error("$t not supported")
	call(::Type{AxisScale{:lin}}) = new{:lin}()
	call(::Type{AxisScale{:log2}}) = new{:log2}()
	call(::Type{AxisScale{:log10}}) = new{:log10}()
	call(::Type{AxisScale{:dB10}}) = new{:dB10}()
	call(::Type{AxisScale{:dB20}}) = new{:dB20}()

	#Aliases:
	call(::Type{AxisScale{:log}}) = new{:log10}()
end
AxisScale(t::Symbol) = AxisScale{t}()

#Specifies configuration of axes:
abstract Axes

#Rectilinear axis (ex: normal cartesian +logarithmic, ...):
immutable AxesRect <: Axes
	xscale::AxisScale
	yscale::AxisScale
end
AxesRect(xscale::Symbol, yscale::Symbol) = AxesRect(AxisScale{xscale}(), AxisScale{yscale}())
AxesRect() = AxesRect(:lin, :lin)

#Curvilinear axes (ex: polar plots):
immutable AxesCurv <: Axes
	rscale::AxisScale #Radial scale could be
end
AxesCurv(rscale::Symbol=:lin) = AxesCurv(AxisScale{rscale}())

immutable AxesSmith{T} <: Axes
	ref::Float64 #Y/Zref
	call{T}(t::Type{AxesSmith{T}}, ref::Real) = error("$t not supported")
	call(::Type{AxesSmith{:Z}}, ref::Real) = new{:Z}(ref)
	call(::Type{AxesSmith{:Y}}, ref::Real) = new{:Y}(ref)
end
AxesSmith(t::Symbol, ref::Real=1) = AxesSmith{t}(ref)
AxesSmith(t::Symbol, ref::Void) = AxesSmith(t)

type Waveform{T}
	id::DisplayString
	ds::T
	line::LineAttributes
	glyph::GlyphAttributes
	ext::PExtents2D
end

#Input waveform:
typealias IWaveform Waveform{IDataset}

#Display waveform (concrete (x, y) pair):
typealias DWaveform Waveform{Vector{Point2D}}

#=TODO:
typealias DWaveformCplx Waveform{PointComplex}
to track x-values of complex data??
=#

type Annotation
	title::DisplayString
	xlabel::DisplayString
	ylabel::DisplayString
	timestamp::DisplayString
end
Annotation(;title="") = Annotation(title, "", "", Libc.strftime(time()))

type Font
	name::AbstractString
	_size::Float64
	bold::Bool
end
Font(name::AbstractString, _size::Real; bold::Bool=false) = Font(name, _size, bold)
Font(_size::Real; bold::Bool=false) = Font(defaults.fontname, _size, bold)

#Legend layout/style
type LegendLStyle
	enabled::Bool
	#autosize::Bool Auto-compute width from text.
	font::Font
	#position: left/right/...
	width::Float64
	#border::Bool
	#wborder
	vgap::Float64 #Vertical spacing between (% of text height).
	linegap::Float64 #Spacing between line and label (% of M character).
	linelength::Float64 #length of line segment to display.
end

#Plot layout
#TODO: Split Layout into "StyleInfo" - which includes Layout??
type Layout
	htitle::Float64 #Title allocation
	waxlabel::Float64 #Vertical axis label allocation (width)
	haxlabel::Float64 #Horizontal axis label allocation (height)
	wnolabels::Float64 #Width to use where no labels are displayed

	wticklabel::Float64 #y-axis values allocation (width)
	hticklabel::Float64 #x-axis values allocation (height)
	hlabeloffset::Float64
	vlabeloffset::Float64

	tframe::Float64 #Frame thickness

	wdata::Float64 #Suggested width of data (graph) area
	hdata::Float64 #Suggested height of data (graph) area

	fnttitle::Font
	fntaxlabel::Font
	fntticklabel::Font
	fnttime::Font

	grid::GridAttributes
	legend::LegendLStyle
	xlabelformat::TickLabelStyle
	ylabelformat::TickLabelStyle
	showtimestamp::Bool
end
Layout(;fontname::AbstractString=defaults.fontname, fontscale=defaults.fontscale) =
	Layout(
	20*fontscale, 20*fontscale, 20*fontscale, 45, #Title/main labels
	60*fontscale, 15*fontscale, 3, 7, 2, #Ticks/frame
	defaults.wdata, defaults.hdata,
	Font(fontname, fontscale*14, bold=true), #Title
	Font(fontname, fontscale*14), Font(fontname, fontscale*12), #Axis/Tick labels
	Font(fontname, fontscale*8), #Time stamp
	GridAttributes(true, false, true, false),
	LegendLStyle(false, Font(fontname, fontscale*12), 100.0, .25, 0.5, 20),
	TickLabelStyle(), TickLabelStyle(), defaults.showtimestamp
)

#Tag data as being part of a given coordinate system:
immutable CoordSystem{ID}; end
typealias DeviceCoord CoordSystem{:dev}
typealias NormCoord CoordSystem{:norm}
typealias DataCoord CoordSystem{:data}

immutable TypedCoord{CT<:CoordSystem}
	v::DReal
end
#Annotation coordinates can match data or be normalized to plot bounds (0 -> 1):
typealias AnnotationCoord Union{TypedCoord{NormCoord}, TypedCoord{DataCoord}}
coord(s::Symbol, v::DReal) = TypedCoord{CoordSystem{s}}(v)

type TextAnnotation
	text::DisplayString
	pt::Point2D #Data coordinates - set NaN to use offsets only
	xoffset::DReal #Normalized to [0,1] plot bounds
	yoffset::DReal #Normalized to [0,1] plot bounds
	font::Font
	color::Colorant
	angle::DReal #Degrees
	align::Symbol #tl, tc, tr, cl, cc, cr, bl, bc, br
end
#Don't use "text"... high probability of collisions when exported...
atext(text::AbstractString; x::Real=DNaN, y::Real=DNaN, xoffset=0, yoffset=0,
	font=Font(10), color=COLOR_BLACK, angle=0, align=:bl) =
	TextAnnotation(text, Point2D(x,y), xoffset, yoffset, font, color, angle, align)

type HVMarker
	vmarker::Bool #false: hmarker
	pos::DReal
	line::LineAttributes
end
vmarker(pos, line=InspectDR.line()) = HVMarker(true, pos, line)
hmarker(pos, line=InspectDR.line()) = HVMarker(false, pos, line)

#2D plot.
type Plot2D <: Plot
	layout::Layout
	axes::Axes
	annotation::Annotation

	#Plot extents (access using getextents):
	ext_data::PExtents2D #Maximum extents of data
	ext_full::PExtents2D #Defines "full" zoom when combined with ext_data
	ext::PExtents2D #Current/active extents

	data::Vector{IWaveform}
	markers::Vector{HVMarker}
	atext::Vector{TextAnnotation}

	#Display data cache:
	invalid_ddata::Bool #Is cache of display data invalid?
	display_data::Vector{DWaveform} #Clipped to current extents

	#Maximum # of x-pts in display:
	#TODO: move to layout?
	xres::Int
end

Plot2D(;title="") = Plot2D(Layout(), AxesRect(), Annotation(title=title),
	PExtents2D(), PExtents2D(), PExtents2D(), [], [], [], true, [], 1000
)

type Multiplot
	title::DisplayString

	subplots::Vector{Plot}
	ncolumns::Int

	#Default width/height of plots
	wplot::Float64
	hplot::Float64

	htitle::Float64 #Title height allocation
	fnttitle::Font
end
Multiplot(;title="", ncolumns=1, titlefont=Font(defaults.fontname, 20),
		wplot=defaults.wplot, hplot=defaults.hplot) =
	Multiplot(title, [], ncolumns, wplot, hplot, 30, titlefont)

#==Constructor-like functions
===============================================================================#

#axes function:
#Interface is a bit irregular, but should be easy to use...
#-------------------------------------------------------------------------------
function axes(a1::Symbol)
	if :smith == a1
		return AxesSmith(:Z)
	elseif :polar == a1
		return AxesCurv()
	else
		throw(MethodError(axes, (a1,)))
	end
end

function axes(a1::Symbol, a2::Symbol; ref = nothing)
	if :smith == a1
		return AxesSmith(a2, ref)
	elseif ref != nothing
		error("cannot set ref for axes(:$a1, :$a2)")
	end

	if :polar == a1
		return AxesCurv(a2)
	else
		return AxesRect(a1, a2)
	end
end

function axes(a1::Symbol, a2::Symbol, a3::Symbol)
	if :polar == a1
		return AxesCurv(a2, a3)
	else
		throw(MethodError(axes, (a1,a2,a3)))
	end
end

#==Accessors
===============================================================================#
_width(style::LegendLStyle) = style.width

getextents(d::IWaveform) = getextents(d.ds)
function getextents(dlist::Vector{IWaveform})
	result = PExtents2D(DNaN, DNaN, DNaN, DNaN)
	for d in dlist
		result = union(result, d.ext)
	end
	return result
end

#==Mutators
===============================================================================#
settitle(mplot::Plot2D, title::DisplayStringArg) =
	(mplot.annotation.title = DisplayString(title))

settitle(mplot::Multiplot, title::DisplayStringArg) =
	(mplot.title = DisplayString(title))

#=="add" interface
===============================================================================#

function _add(mp::Multiplot, plot::Plot2D)
	push!(mp.subplots, plot)
	return plot
end
_add{T<:Plot}(mp::Multiplot, ::Type{T}) = _add(mp, T())

#Set dataf1=false to overwrite optimizations for functions of 1 argument.
function _add(plot::Plot2D, x::Vector, y::Vector; id::AbstractString="", dataf1=true)
	if dataf1
		dataf1 = isincreasing(x) #Can we use optimizations?
	end
	ext = PExtents2D() #Don't care at the moment
	ds = IWaveform(id, IDataset{dataf1}(x, y), line(), glyph(), ext)
	push!(plot.data, ds)
	return ds
end

function _add(plot::Plot2D, marker::HVMarker)
	push!(plot.markers, marker)
	return marker
end

function _add(plot::Plot2D, a::TextAnnotation)
	push!(plot.atext, a)
	return a
end

#==Mapping/interpolation functions
===============================================================================#

#Mapping functions, depending on axis type:
#-------------------------------------------------------------------------------
datamap{T<:Number}(::Type{T}, ::AxisScale) = (x::T)->DReal(x)
datamap{T<:Number}(::Type{T}, ::AxisScale{:dB10}) = (x::T)->(5*log10(DReal(abs2(x))))
datamap{T<:Number}(::Type{T}, ::AxisScale{:dB20}) = (x::T)->(10*log10(DReal(abs2(x))))
datamap{T<:Number}(::Type{T}, ::AxisScale{:log2}) =
	(x::T)->(x<0? DNaN: log2(DReal(x)))
datamap{T<:Number}(::Type{T}, ::AxisScale{:log10}) =
	(x::T)->(x<0? DNaN: log10(DReal(x)))
#TODO: find a way to show negative values for log10?

datamap_rev{T<:Number}(::Type{T}, ::AxisScale) = (x::T)->DReal(x)
#NOTE: Values dBs remain in dBs for readability
datamap_rev{T<:Number}(::Type{T}, ::AxisScale{:log2}) = (x::T)->(2.0^x)
datamap_rev{T<:Number}(::Type{T}, ::AxisScale{:log10}) = (x::T)->(10.0^x)

datamap_rev{T<:Number}(v::T, s::AxisScale) = datamap_rev(T,s)(v) #One-off conversion

#Extents mapping functions, depending on axis type:
#-------------------------------------------------------------------------------
extentsmap(::AxisScale) = (x::DReal)->x #Most axis types don't need to re-map extents.
#NOTE: Extents in dBs remain extents in dBs
extentsmap(t::AxisScale{:log2}) = datamap(DReal, t)
extentsmap(t::AxisScale{:log10}) = datamap(DReal, t)

extentsmap_rev(::AxisScale) = (x::DReal)->x #Most axis types don't need to re-map extents.
#NOTE: Extents in dBs remain extents in dBs
extentsmap_rev(t::AxisScale{:log2}) = datamap_rev(DReal, t)
extentsmap_rev(t::AxisScale{:log10}) = datamap_rev(DReal, t)

#==Plot extents
===============================================================================#

function invalidate_extents(plot::Plot2D)
	#If extents are no longer valid, neither is display cache:
	plot.invalid_ddata = true
end

function invalidate_datalist(plot::Plot2D)
	plot.invalid_ddata = true
end

rescale(ext::PExtents2D, axes::Axes) = ext #Default
function rescale(ext::PExtents2D, axes::AxesRect)
	xmap = extentsmap(axes.xscale)
	ymap = extentsmap(axes.yscale)
	return PExtents2D(
		xmap(ext.xmin), xmap(ext.xmax),
		ymap(ext.ymin), ymap(ext.ymax)
	)
end
rescale_rev(ext::PExtents2D, axes::Axes) = ext #Default
function rescale_rev(ext::PExtents2D, axes::AxesRect)
	xmap = extentsmap_rev(axes.xscale)
	ymap = extentsmap_rev(axes.yscale)
	return PExtents2D(
		xmap(ext.xmin), xmap(ext.xmax),
		ymap(ext.ymin), ymap(ext.ymax)
	)
end

#Accessor:
getextents(plot::Plot2D) = plot.ext
getextents_xfrm(plot::Plot2D) = rescale(plot.ext, plot.axes)

#Full extents are always merged (ext_full expected to be incomplete):
getextents_full(plot::Plot2D) = merge(plot.ext_data, plot.ext_full)

#Set active plot extents using data coordinates:
function setextents(plot::Plot2D, ext::PExtents2D, hallowed::Bool=true, vallowed::Bool=true)
	xmin = ext.xmin; xmax = ext.xmax; ymin = ext.ymin; ymax = ext.ymax
	if !hallowed
		xmin = plot.ext.xmin
		xmax = plot.ext.xmax
	end
	if !vallowed
		ymin = plot.ext.ymin
		ymax = plot.ext.ymax
	end
	#Automatically fill-in any NaN fields, if possible:
	plot.ext = merge(getextents_full(plot), PExtents2D(xmin, xmax, ymin, ymax))
	invalidate_extents(plot)
end

#Set active plot extents using xfrm display coordinates:
setextents_xfrm(plot::Plot2D, ext::PExtents2D, hallowed::Bool=true, vallowed::Bool=true) =
	setextents(plot, rescale_rev(ext, plot.axes), hallowed, vallowed)

function setextents_full(plot::Plot2D, ext::PExtents2D)
	plot.ext_full = ext
end

#==Plot/graph bounding boxes
===============================================================================#

aspect_square(::Axes) = false
aspect_square(::AxesSmith) = true

#Returns a centered bounding box with square aspect ratio.
function squarebounds(bb::BoundingBox)
	w = width(bb); h = height(bb)
	xmin = bb.xmin; xmax = bb.xmax
	ymin = bb.ymin; ymax = bb.ymax
	if w < h
		hdim = w / 2
		c = (ymin + ymax) / 2
		ymin = c - hdim; ymax = c + hdim
	else
		hdim = h / 2
		c = (xmin + xmax) / 2
		xmin = c - hdim; xmax = c + hdim
	end
	return BoundingBox(xmin, xmax, ymin, ymax)
end

#Get bounding box of graph (plot data area):
function graphbounds(plotb::BoundingBox, lyt::Layout)
	xmin = plotb.xmin + lyt.waxlabel + lyt.wticklabel
	xmax = plotb.xmax
	xmax -= lyt.legend.enabled? _width(lyt.legend): lyt.wnolabels
	ymin = plotb.ymin + lyt.htitle
	ymax = plotb.ymax - lyt.haxlabel - lyt.hticklabel

	#Avoid division by zero, inversions, ...
	if xmin >= xmax
		c = (xmin + xmax)/2
		xmin = c - 0.5
		xmax = c + 0.5
	end
	if ymin >= ymax
		c = (ymin + ymax)/2
		ymin = c - 0.5
		ymax = c + 0.5
	end

	return BoundingBox(xmin, xmax, ymin, ymax)
end

function graphbounds(plotb::BoundingBox, lyt::Layout, axes::Axes)
	graphbb = graphbounds(plotb, lyt)
	if aspect_square(axes)
		graphbb = squarebounds(graphbb)
	end
	return graphbb
end

#Get bounding box of entire plot:
function plotbounds(lyt::Layout, graphw::Float64, graphh::Float64)
	xmax = graphw + lyt.waxlabel + lyt.wticklabel
	xmax += lyt.legend.enabled? _width(lyt.legend): lyt.wnolabels
	ymax = graphh + lyt.htitle + lyt.haxlabel + lyt.hticklabel
	return BoundingBox(0, xmax, 0, ymax)
end

#Get suggested plot bounds:
function plotbounds(lyt::Layout, axes::Axes)
	wdata = lyt.wdata; hdata = lyt.hdata
	if aspect_square(axes)
		wdata = hdata = min(wdata, hdata)
	end
	return plotbounds(lyt, wdata, hdata)
end

#==Pre-processing display data
===============================================================================#

function _reduce(input::IWaveform, ext::PExtents2D, xres_max::Integer)
	return DWaveform(input.id, _reduce(input.ds, ext, xres_max), input.line, input.glyph, input.ext)
end

_reduce(inputlist::Vector{IWaveform}, ext::PExtents2D, xres_max::Integer) =
	map((input)->_reduce(input, ext, xres_max::Integer), inputlist)

#Rescale input dataset:
#-------------------------------------------------------------------------------
function _rescale{T<:Number}(d::Vector{T}, scale::AxisScale)
	#Apparently, passing functions as arguments is not efficient in Julia.
	#-->Specializing on AxisScale, hoping to improve efficiency on dmap:
	dmap = datamap(T, scale)

	result = Array(DReal, length(d))
	for i in 1:length(d)
		result[i] = dmap(d[i])
	end
	return result
end
_rescale{T<:Number}(d::Vector{T}, scale::AxisScale{:lin}) = d #Optimization: Linear scale does not need re-scaling
_rescale{T<:IDataset}(input::T, xscale::AxisScale, yscale::AxisScale) =
	T(_rescale(input.x, xscale), _rescale(input.y, yscale))

function _rescale(input::IWaveform, xscale::AxisScale, yscale::AxisScale)
	ds = _rescale(input.ds, xscale, yscale)
	return IWaveform(input.id, ds, input.line, input.glyph, getextents(ds))
end

#Specialized on xscale/yscale for efficiency:
_rescale(inputlist::Vector{IWaveform}, xscale::AxisScale, yscale::AxisScale) =
	map((input)->_rescale(input, xscale, yscale), inputlist)

_rescale(inputlist::Vector{IWaveform}, axes::AxesRect) = _rescale(inputlist, axes.xscale, axes.yscale)

_rescale(pt::Point2D, xscale::AxisScale, yscale::AxisScale) =
	Point2D(datamap(DReal, xscale)(pt.x), datamap(DReal, yscale)(pt.y))
_rescale(pt::Point2D, axes::AxesRect) = _rescale(pt, axes.xscale, axes.yscale)

#Preprocess input dataset (rescale/reduce quantity of data/...):
#-------------------------------------------------------------------------------
#   (Updates display_data)
function preprocess_data(plot::Plot2D)
	#TODO: Find a way to preprocess x-vectors referencing same data only once?

	#Rescale data
	wfrmlist = _rescale(plot.data, plot.axes)

	#Update extents:
	plot.ext_data = rescale_rev(getextents(wfrmlist), plot.axes) #Update max extents
	setextents(plot, plot.ext) #Update extents, resolving any NaN fields.
	ext = getextents_xfrm(plot) #Read back extents, in transformed coordinates

	#Reduce data:
	plot.display_data = _reduce(wfrmlist, ext, plot.xres)
	plot.invalid_ddata = false

	#NOTE: Rescaling before data reduction is somewhat inefficient, but makes
	#      it easier to interpolate data in _reduce step.
	#TODO: Find a way to efficiently reduce before re-scaling???
end

#==High-level display functions
===============================================================================#
function update_ddata(plot::Plot2D)
	invalidate_extents(plot) #Always compute below:
	#TODO: Conditionnaly compute (only when data changed/added)?

	if plot.invalid_ddata
		preprocess_data(plot)
	end
end

#Last line


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdtimeline

"""
    'feffery'_antdtimeline(;kwargs...)

An AntdTimeline component.

Keyword arguments:
- `id` (String; optional)
- `className` (String; optional)
- `items` (optional): . items has the following type: Array of lists containing elements 'content', 'color', 'icon', 'iconStyle', 'contentStyle', 'label'.
Those elements have the following types:
  - `content` (String; optional)
  - `color` (String; optional)
  - `icon` (String; optional)
  - `iconStyle` (Dict; optional)
  - `contentStyle` (Dict; optional)
  - `label` (String; optional)s
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `mode` (a value equal to: 'left', 'alternate', 'right'; optional)
- `pending` (String; optional)
- `reverse` (Bool; optional)
- `style` (Dict; optional)
"""
function 'feffery'_antdtimeline(; kwargs...)
        available_props = Symbol[:id, :className, :items, :loading_state, :mode, :pending, :reverse, :style]
        wild_props = Symbol[]
        return Component("'feffery'_antdtimeline", "AntdTimeline", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
@testset "Test @instate" begin

  function algo(stp :: AbstractStopping, n :: Int)

    OK = start!(stp)
    x  = stp.current_state.x

    while !OK
        x = x .+ 1

        OK = update_and_stop!(stp, x = x)
    end

    return stp
  end

  function instatealgo(stp :: AbstractStopping, n :: Int)

    state = stp.current_state
    x  = state.x
    OK = start!(stp)

    @instate state while !OK
        x = x .+ 1

        OK = stop!(stp)
    end

    return stp
  end

  n = 10
  stp1 = GenericStopping(x->0., zeros(n), max_iter = n, rtol = 0.)
  stp2 = GenericStopping(x->0., zeros(n), max_iter = n, rtol = 0.)

  stp1 = algo(stp1, n)
  stp2 = algo(stp2, n)

  @test stp1.current_state.x == stp2.current_state.x

  #=
    Suggestions:
     - It would be better to say stp.current_state instead of state.
  =#
end


text:
# Different types of site that can be counted
# ===========================================
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

abstract Site
abstract Mutation <: Site

"""
A `Certain` site describes a site where both of two aligned sites are not an
ambiguity symbol or a gap.
"""
immutable Certain <: Site end

"""
An `Ambiguous` site describes a site where either of two aligned sites are an
ambiguity symbol.
"""
immutable Ambiguous <: Site end

"""
An `Gap` site describes a site where either of two aligned sites are a
gap symbol '-'.
"""
immutable Gap <: Site end

"""
A `Match` site describes a site where two aligned nucleotides are the
same biological symbol.
"""
immutable Match <: Site end

"""
A `Mismatch` site describes a site where two aligned nucleotides are not the
same biological symbol.
"""
immutable Mismatch <: Site end

"""
A `Mismatch` site describes a site where two aligned nucleotides are definately
conserved. By definately conserved this means that the symbols of the site are
non-ambiguity symbols, and they are the same symbol.
"""
immutable Conserved <: Mutation end

"""
A `Mutated` site describes a site where two aligned nucleotides are definately
mutated. By definately mutated this means that the symbols of the site are
non-ambiguity symbols, and they are not the same symbol.
"""
immutable Mutated <: Mutation end

"""
A `Transition` site describes a site where two aligned nucleotides are definately
mutated, and the type of mutation is a transition mutation.
In other words, the symbols must not be ambiguity symbols, and they must
be different such that they constitute a transition mutation: i.e. A<->G, or C<->T.
"""
immutable Transition <: Mutation end

"""
A `Transversion` site describes a site where two aligned nucleotides are definately
mutated, and the type of mutation is a transversion mutation.
In other words, the symbols must not be ambiguity symbols, and they must
be different such that they constitute a transversion mutation: i.e. A<->C,
A<->T, G<->T, G<->C.
"""
immutable Transversion <: Mutation end


text:
function act_1_deck(cur_T_bar::AbstractSymbol=symbols(:T_bar); cur_kwargs...)
  cur_dict = Dict(
    :deck => :act_1,
    :T_bar => cur_T_bar,
    :is_pulsed => false,
    :H => 1.65,
    :Q => 42.5,
    :epsilon => 0.25,
    :kappa_95 => 2.1,
    :delta_95 => 0.4,
    :nu_n => 0.27,
    :nu_T => 1.15,
    :Z_eff => 2.11,
    :f_D => 0.75,

    :A => 2.5,
    :l_i => 0.35906,
    # :rho_m => .82,

    :N_G => 1.0,
    :eta_T => 0.575,
    :eta_RF => 0.4,

    :max_beta_N => 0.0427,
    :max_q_95 => 2.5,

    :max_P_E => 1500.0,

    :max_P_W => 2.45,
    :max_q_DV => 13.7,

    :B_CS => 12.77,
    :k => 0.333,
    :sigma_CS => 660,
    :sigma_TF => 660,
    :J_CS => 50,
    :J_TF => 230,
    :C_ejima => 0.3,
    :tau_FT => 1.6e9,

    :eta_CD => 0.188,
    :wave_theta => 135.0,
    :eta_LH => 0.6
  )

  merge!(cur_dict, Dict(cur_kwargs))

  BaseReactor(cur_dict)
end


text:
function _gkp_global(g::AbstractGraph{T}) where T <: Integer
    n = nv(g)
    sinks = get_sinks(g)
    sources = get_sources(g)
    source_idx = falses(n)
    source_idx[sources] .= true
    visited = falses(n)
    visited[sinks] .= true

    w = 1 ./ indegree(g)
    od = outdegree(g)
    kp = zeros(n)
    kp[sinks] .= 1.0

    layer = sinks
    while length(layer) > 0
        new = T[]

        for v in layer
            inb = inneighbors(g, v)

            for nb in inb
                kp[nb] += w[v] * kp[v]
                od[nb] -= 1

                if !source_idx[nb] && od[nb] == 0
                    if !visited[nb]
                        push!(new, nb)
                        visited[nb] = true
                    end
                end
            end
        end
        layer = new
    end

    kp[sinks] .= 0
    return kp
end

function _gkp_local(g::AbstractGraph{T}) where T <: Integer
    n = nv(g)
    sinks = get_sinks(g)
    sources = get_sources(g)
    source_idx = falses(n)
    source_idx[sources] .= true
    visited = falses(n)
    visited[sinks] .= true
    l = zeros(Int, n)
    l[sinks] .= 1

    w = 1 ./ indegree(g)
    od = outdegree(g)
    kp = zeros(n)
    kp[sinks] .= 1.0

    layer = sinks; lnum = 1
    while length(layer) > 0
        new = T[]

        for v in layer
            inb = inneighbors(g, v)
            l[v] = lnum

            for nb in inb
                kp[nb] += w[v] * kp[v]
                od[nb] -= 1

                if !source_idx[nb] && od[nb] == 0
                    if !visited[nb]
                        push!(new, nb)
                        visited[nb] = true
                    end
                end
            end
        end
        layer = new; lnum += 1
    end

    l[sources] .= lnum
    kp[sinks] .= 0.0
    return kp, l

end

"""
    GKP(normalize=:none)(g)

Struct representing Genetic Knowledge Persistence (GKP) vertex weights.
This is a modified version of the method proposed in Martinelli & Nomaler (2014)
where citations from future patents to patents before the curent layer are not ignored
in the persistence computation.

The struct can be called to compute weights directly or
can be passed to the `mainpath` function for dispatch.
"""
Base.@kwdef struct GKP <: MainPathVertexWeight
    normalize::Symbol = :none
end

function (x::GKP)(g)
    if x.normalize == :global
        res = _gkp_global(g)
        return res ./ maximum(res)
    elseif x.normalize == :local
        res, layers = _gkp_local(g)
        for l in unique(layers)
            res[layers .== l] .= res[layers .== l] ./ maximum(res[layers .== l])
        end
        res[isnan.(res)] .= 0.0
        return res
    elseif x.normalize == :none
        return _gkp_global(g)
    else
        throw(ArgumentError("Viable options for normalize are :global, :local or :none."))
    end
end


text:
using
    Documenter,
    AutocorrelationShell

DocMeta.setdocmeta!(AutocorrelationShell, :DocTestSetup, :(using AutocorrelationShell); recursive=true)

makedocs(
    sitename="AutocorrelationShell.jl",
    modules = [AutocorrelationShell],
    authors = "Shozen Dan, Rishi Subramanian, Christina Chang, and Naoki Saito",
    repo="https://github.com/ShozenD/AutocorrelationShell.jl/blob/{commit}{path}#{line}",

    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ShozenD.github.io/AutocorrelationShell.jl",
        assets=String[],
    ),

    pages = [
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "AC Wavelet 1D" => "acw1d.md",
        "AC Wavelet 2D" => "acw2d.md",
        "AC Wavelet Packets" => "acwpt.md",
    ]
)

deploydocs(
    repo = "https://github.com/ShozenD/AutocorrelationShell.jl"
)

max_stars_repo_name: ShozenD/AutocorrelationShell

max_issues_repo_name: ShozenD/AutocorrelationShell

max_forks_repo_name: ShozenD/AutocorrelationShell


text:
## First order perturbation approach
## Local solution for log-linear model around a non-stochastic steady state
## Blanchard and Kahn (1980)
## Based on Exercises for Practical DSGE Modelling
module DsgeSolves
using Distributions

type DSGE
  b::Int64
  f::Int64
  s::Int64
  A::Array
  B::Array
  C::Array
  v::Distribution
end

function DSGE(b,f,s)
  A = zeros(b+f,b+f)
  B = zeros(b+f,b+f)
  C = zeros(b+f,s)
  v = Normal(0,1)
  return DSGE(b,f,s,A,B,C,v)
end

function create_matrices(m)
  A = inv(m.A)*m.B
  B = inv(m.A)*m.C
  (λ,P) = real.(eig(A))
  λ = real(λ)
  t = sortrows([λ P'])
  Λ = diagm(t[:,1])
  P = t[:,2:(m.b+m.f+1)]'
  pstar = inv(P)
  R = pstar*B
 return A, B, Λ, P, pstar, R
end

function create_partition(m,pstar,Λ,R)
  Λ1 = Λ[1:m.b,1:m.b]
  Λ2 = Λ[(m.b+1):(m.b+m.f),(m.b+1):(m.b+m.f)]

  P11 = pstar[1:m.b,1:m.b]
  P12 = pstar[1:m.b,(m.b+1):(m.b+m.f)]
  P21 = pstar[(m.b+1):(m.b+m.f),1:m.b]
  P22 = pstar[(m.b+1):(m.b+m.f),(m.b+1):(m.b+m.f)]

  R1 = R[1:m.s,1:m.s]
  R2 = R[(m.b+1:m.s),1:m.s]
  return Λ1, Λ2, P11, P12, P21, P22, R1, R2
end

function simulating_dsge(m,pstar,Λ,R,Σ,rep)
  (Λ1, Λ2, P11, P12, P21, P22, R1, R2) = create_partition(m,pstar,Λ,R)
  w = zeros(m.b,rep)
  y = zeros(m.f,rep)
    for i in 2:rep
      w[:,i] = real(inv(P11-P12*inv(P22)*P21)*Λ1*(P11-P12*inv(P22)*P21))*w[:,i-1]+real((P11-P12*inv(P22)*P21)*R1)*Σ*rand(m.v,m.s)
      y[:,i] = -real(inv(P22)*P21)*w[:,i]
    end
return w, y
end

function irf_dsge(m,pstar,Λ,R,shk,rep)
  (Λ1, Λ2, P11, P12, P21, P22, R1, R2) = create_partition(m,pstar,Λ,R)
  shock = zeros(m.s,rep)
  shock[shk,2] = 1
  w = zeros(m.b,rep)
  y = zeros(m.f,rep)
    for i in 2:rep
      w[:,i] = real(inv(P11-P12*inv(P22)*P21)*Λ1*(P11-P12*inv(P22)*P21))*w[:,i-1]+real((P11-P12*inv(P22)*P21)*R1)*shock[:,i]
      y[:,i] = -real(inv(P22)*P21)*w[:,i]
    end
return w, y
end

function fevd(m,pstar,Λ,R,rep) ##restart from this function using a for loop
                               ## with the function irf_dsge
  (Λ1, Λ2, P11, P12, P21, P22, R1, R2) = create_partition(m,pstar,Λ,R)
  irf_w = zeros(m.b,rep,m.b*m.s)
  irf_y = zeros(m.f,rep,m.f*m.s)
  for j in 1:m.s
    shock = zeros(m.s,rep)
    shock[j,2] = 1
    for i in 2:rep
      w[:,i] = real(inv(P11-P12*inv(P22)*P21)*Λ1*(P11-P12*inv(P22)*P21))*w[:,i-1]+real((P11-P12*inv(P22)*P21)*R1)*shock[:,i]
      y[:,i] = -real(inv(P22)*P21)*w[:,i]
    end
    irf_w[:,:,j] = w
    irf_y[:,:,j] = y
  end
  sum(irf_w[1,:,1].^2)Σ[j,j]/(sum(irf_w[1,:,1].^2) + sum(irf_w[1,:,2].^2))
    irf_w[1,i,1]
    den =

return irf_w, irf_y
end

export DSGE, create_matrices, simulating_dsge, irf_dsge

end # end module


text:
include("style.jl")
using InteractiveChaos, Random
using DynamicalBilliards
import Makie
using Makie: to_color, RGBf0, RGBAf0
using DynamicalSystems

# Set style to book colors
InteractiveChaos.obcolor(::Antidot) = to_color(COLORS[1])
InteractiveChaos.obcolor(::Obstacle) = to_color(COLORS[1])
InteractiveChaos.obfill(o::Antidot) = RGBAf0(0,0,0,0)
InteractiveChaos.obls(::Antidot) = nothing

# %% Circle billiard plot
bd = Billiard(Antidot(Float32[0, 0], 1.0f0, false))

# make three particles
js = [(0.2, sind(30)), (1.5, sind(30 + sqrt(2)/10))]

colors = [to_color(COLORS[i+1]) for i in 1:length(js)]

ps = Particle{Float32}[]
for (ξ, s) in js
    x, y = from_bcoords(Float32(ξ), Float32(s), bd)
    push!(ps, Particle(x..., atan(y[2], y[1])))
end

scene, bmapax = billiard_bmap_plot(bd, ps;
colors = colors, tail = 100000, steps = 100003, backgroundcolor = RGBf0(1,1,1),
ms = 10, vr = 0.1)
Makie.ylims!(bmapax, 0.2, 0.9)
# Makie.save(joinpath(figdir, "circlebilliard.png"), scene)

# %% Chaotic billiard (sinai)
InteractiveChaos.obcolor(::Antidot) = to_color(COLORS[1])
InteractiveChaos.obfill(o::Antidot) = RGBAf0(0,0,0,0)
InteractiveChaos.obls(::Antidot) = nothing

colors = [to_color(COLORS[i+1]) for i in 1:2]

bd = billiard_sinai(0.25f0, 1f0, 1f0)

ps = particlebeam(0.2f0, 0.75, π + π/4 + 0.414235, 100, 0.002)
scene, bmapax = billiard_bmap_plot(bd, ps; colors = colors,
tail = 3100, steps = 3400, backgroundcolor = RGBf0(1,1,1),
ms = 10, vr = 0.05)

# %% Mushroom
w = 0.2f0
bd = billiard_mushroom(1f0, w, 1f0, 0f0)

pc = MushroomTools.randomchaotic(1f0, w, 1f0)
pc = Particle(-0.01f0, 0.2f0, sqrt(3f0)/2)

pr = Particle(0.0f0, 1.2f0, 0.0f0)
pr2 = Particle(0.0f0, 1.9f0, 0.0f0)

ps = [pc, pr, pr2]
colors = [to_color(COLORS[i+1]) for i in 1:length(ps)]

scene, bmapax = billiard_bmap_plot(bd, ps;
colors = colors, tail = 100000, steps = 1000003, backgroundcolor = RGBf0(1,1,1),
ms = 5)

# interactive_billiard_bmap(bd)

# bmapax.xticklabelsize = 32
# bmapax.yticklabelsize = 32
bmapax.xlabelsize = 40
bmapax.ylabelsize = 40
bmapax.ylabelpadding = 15
# Makie.save(joinpath(figdir, "mushroom.png"), scene)

# %% Natural measure of henon map
using DynamicalSystems
# using PyPlot; using3D()

ε = 0.0001
x, y, p = begin
    ds = Systems.henon()
    X = trajectory(ds, 100000000; Ttr = 100)
    p, b = ChaosTools.binhist(ε, X)
    x = [a[1] for a in b]
    y = [a[2] for a in b]
    x, y, p ./ maximum(p)
end

x ./= 2

# z = zeros(length(x))
# cmap = get_cmap("viridis")
# cs = cmap.(sqrt.(p ./ maximum(p)))
# PyPlot.figure()
# PyPlot.bar3D(x, y, z, ε/2, ε/2, p, cs)
# PyPlot.axis("off")

# Makie alternative:
using Makie
using Makie: Point3f0, Vec3f0, Rect3D

scene = Makie.meshscatter(vec(Point3f0.(x, y, 0.0)),
    markersize=Vec3f0.(2ε, 2ε, p), marker=Rect3D(Vec3f0(0), Vec3f0(1)),
    limits=Rect3D(Vec3f0(0), Vec3f0(1)),
    color = clamp.(p, 0, 0.2))
ylims!(scene, 1.5 .* extrema(y)...)
xlims!(scene, extrema(x)...)
xticks!(scene; )
zlabel!(scene, "ρ")
display(scene)

scene
# %%

# Same but for standard map
ε = 0.01
us = [
    SVector(0.000767806, 0.00054896),
    SVector(π, 0.5),
    SVector(π, 1.0),
]
sm = Systems.standardmap(; k = 0.9)

scene = Scene()
for u in us
    x, y, p = begin
        X = trajectory(sm, 10^7, u; Ttr = 100)
        p, b = ChaosTools.binhist(ε, X)
        x = [a[1] for a in b]
        y = [a[2] for a in b]
        x, y, p ./ maximum(p)
    end
    using Makie

    Makie.meshscatter!(scene, vec(Point3f0.(x, y, 0.0)),
    markersize=Vec3f0.(ε, ε, p), marker=Rect3D(Vec3f0(0), Vec3f0(1)),
    limits=Rect3D(Vec3f0(0), Vec3f0(1)),
    color = p, colormap = :viridis)
end

xlims!(scene, 0, 2π)
ylims!(scene, 0, 2π)
xlabel!(scene, "θ")
ylabel!(scene, "p")
zlabel!(scene, "ρ")
display(scene)

# Same but for mushroom billiard
w = 0.2
bd = billiard_mushroom(1, w, 1, 0)

pc = Particle(-0.01, 0.2, sqrt(3)/2)
pr = Particle(0.0, 1.2, 0.0)
pr2 = Particle(0.0, 1.9, 0.0)

cmaps = [:blues, :reds, :greens]
particles = [pc, pr, pr2]

ε = 0.01

scene = Scene()
for i in 1:3
    x, y, p = begin
        bmap, = boundarymap(particles[i], bd, 10^6)
        X = Dataset(bmap)
        p, b = ChaosTools.binhist(ε, X)
        x = [a[1]/2 for a in b]
        y = [a[2] for a in b]
        x, y, p ./ maximum(p)
    end

    Makie.meshscatter!(scene, vec(Point3f0.(x, y, 0.0)),
    markersize=Vec3f0.(ε, ε, p), marker=Rect3D(Vec3f0(0), Vec3f0(1)),
    limits=Rect3D(Vec3f0(0), Vec3f0(1)),
    color = p, colormap = cmaps[i])
end

display(scene)

# %% Standard map island shrinking
using DynamicalSystems
grid = 0.0:0.5:2π
sm = Systems.standardmap()
ks = [0.1, 1.0, 4.0]
fig, axs = subplots(1,3; sharey = true)

for i in 1:3
    set_parameter!(sm, 1, ks[i])
    ax = axs[i]
    for θ in grid, p in grid
        u0 = SVector(θ, p)
        tr = trajectory(sm, 1000, u0)
        c = lyapunovs(sm, 4000, 1; u0 = u0)[1] > 0.01 ? "C0" : "C2"
        ax.plot(columns(tr)...; ls = "None", marker = "o", ms = 0.2, c = c)
    end
    ax.set_xlim(0, 2π)
    ax.set_ylim(0, 2π)
    ax.text(4.5, 5.5, "\$k=$(ks[i])\$"; bbox = bbox)
    ax.set_xticks(0:2:6)
end

axs[1].set_ylabel("\$p\$", labelpad = -5)
axs[2].set_xlabel("\$\\theta\$", labelpad = -15)
fig.subplots_adjust(left = 0.05, bottom = 0.14, right = 0.95, top = 0.97, wspace = 0.06)

# fsave(joinpath(figdir, "kam_sm"), fig)

# %% Recurrences
using DynamicalSystems, DynamicalBilliards
w = 0.2
pc = Particle(-0.01, 0.2, sqrt(3)/2)
pr = Particle(0.0, 1.9, 0.0)

ps = [pc, pr]
cs = ["C1", "C3"]
bd = billiard_mushroom(1, w, 1)

zbox = ((420, 130), (460, 170))

fig, axs = subplots(1, 3)
εs = [0.1, 0.1]
for i in 1:2
    ax = axs[i==1 ? 1 : 3]
    p = ps[i]
    N = 500
    bmap, = boundarymap(p, bd, N)
    tr = regularize(Dataset(bmap))
    R = RecurrenceMatrix(tr, εs[i])
    x, y = coordinates(R)
    ax.scatter(x, y, s = 1, color = cs[i])
    ax.set_aspect("equal")
    ℓ = dl_entropy(R)
    r = rt_entropy(R)
    ax.set_xlim(0, N)
    ax.set_ylim(0, N)
    ax.set_xticks(0:250:500)
    ax.set_yticks(0:250:500)
    ax.grid()
    ax.set_title("\$H_\\ell = $(round(ℓ;digits= 1)),\\, H_r =$(round(r;digits=1))\$")
end
axs[3].set_yticklabels([])

# add zoomin
axis_zoomin!(axs[2], axs[1], zbox, zbox, "C0")
axs[2].axis("off")
axs[2].set_aspect("equal")
bmap, = boundarymap(pc, bd, 500)
tr = regularize(Dataset(bmap))
R = RecurrenceMatrix(tr, εs[1])
x, y = coordinates(R)
axs[2].scatter(x, y, s = 10, color = cs[1])
axs[2].set_xlim(zbox[1][1], zbox[2][1])
axs[2].set_ylim(zbox[1][2], zbox[2][2])

fig.subplots_adjust(left = 0.08, bottom = 0.01, right = 0.95, top = 0.97, hspace = 0.1)
# fsave(joinpath(figdir, "recurrence"), fig)

max_stars_repo_name: moritz-g/ExercisesRepo

max_issues_repo_name: moritz-g/ExercisesRepo

max_forks_repo_name: moritz-g/ExercisesRepo


text:
# Configuration for running simple neutral copy error simulation
#=
Recommended command line to run:
> julia -p 4 -L NeutralEvolution.jl nrun.jl examples/example1
> julia -p 4 -L NeutralEvolution.jl nrun.jl examples/example1 1  # include 1 to set random number seed
=#
export simtype
@everywhere simtype = 3
@everywhere const N = 3        # population size
const mutstddev = 0.04
const ngens = 2
const initial_value = 1.0
const num_trials=1
const record_interval = 1
const log_error = false
const wright_fisher_copy = false
const conformist_probability=0.0
const neutral = true


text:
abstract type AbstractParams{T<:Union{Int64,Nothing}} end
abstract type AbstractParamsB1{T<:Union{Int64,Nothing}} <: AbstractParams{T} end
abstract type AbstractParamsNoB1{T<:Union{Int64,Nothing}} <: AbstractParams{T} end

abstract type AbstractDispersion end
abstract type AbstractSharedCash end

struct ParamsB1{T<:Union{Int64, Nothing}} <: AbstractParamsB1{T}
    N::Int64
    m::T
    W::Float64
    u::Float64
    u₁::Float64
    Δτ::Float64
	ParamsB1(N, m::T, W, u, u₁) where T<:Union{Int64, Nothing} = new{T}(N, m, W, u, u₁, W/N)
end

struct ParamsB1_pinned{T<:Union{Int64, Nothing}} <: AbstractParamsB1{T}
    N::Int64
    m::T
    W::Float64
    u::Float64
    u₁::Float64
    Δτ::Float64
	i_pinned::Int64
	ParamsB1_pinned(N::Int64, m::T, W::Float64, u::Float64, u₁::Float64) where T<:Union{Int64, Nothing} = new{T}(N, m, W, u, u₁, W/N, div(N,2)+1)
end

struct ParamsB1_pinned²{T<:Union{Int64, Nothing}} <: AbstractParamsB1{T}
    N::Int64
    m::T
    W::Float64
    u::Float64
    u₁::Float64
    Δτ::Float64
	i_pinned::Int64
	ParamsB1_pinned²(N::Int64, m::T, W::Float64, u::Float64, u₁::Float64) where T<:Union{Int64, Nothing} = new{T}(N, m, W, u, u₁, W/N, div(N,2)+1)
end

mutable struct ParamsNoB1{T<:Union{Int64, Nothing}} <: AbstractParamsNoB1{T}
    N::Int64
    m::T
    W::Float64
    u::Float64
    Δτ::Float64
	ParamsNoB1(N::Int64, m::T, W::Float64, u::Float64) where T<:Union{Int64, Nothing} = new{T}(N, m, W, u, W/N)
end

function set_period!(params::ParamsNoB1, W)
	params.W = W
	params.Δτ = W/params.N
end

struct ParamsNoB1_pinned{T<:Union{Int64, Nothing}} <: AbstractParamsNoB1{T}
    N::Int64
    m::T
    W::Float64
    u::Float64
    Δτ::Float64
	i_pinned::Int64
	ParamsNoB1_pinned(N::Int64, m::T, W::Float64, u::Float64) where T<:Union{Int64, Nothing} = new{T}(N, m, W, u, W/N, div(N,2)+1)
end

function get_u₀(β::AbstractFloat, psamples_raw::Vector{NTuple{3, Float64}})
    s = 0.0
	for (εₚ⁺, εₚ⁻, wₚ) in psamples_raw
        κₚ = sqrt(εₚ⁻^2 + 1.0)
        s += wₚ*(tanh(β*(εₚ⁺+κₚ)/2) + tanh(β*(κₚ-εₚ⁺)/2))/κₚ
    end
    return s/4.0
end
function get_u₀(psamples_raw)
    s = 0.0
	for (εₚ⁺, εₚ⁻, wₚ) in psamples_raw
        κₚ = sqrt(εₚ⁻^2 + 1.0)
        s += wₚ/κₚ
    end
    return s/2.0
end
function get_u₀(psamples_raw::Vector{NTuple{3, Float64}}, γ::AbstractFloat)
    s = 0.0
	for (εₚ⁺, εₚ⁻, wₚ) in psamples_raw
        κₚ = sqrt(εₚ⁻^2 + γ^2)
        s += wₚ/κₚ
    end
    return s/2.0
end
get_u₀(W, m::Int64, psamples_raw) = get_u₀(W*m, psamples_raw)
get_u₀(W, m::Nothing, psamples_raw) = get_u₀(psamples_raw)

for Params_constr in (:ParamsB1, :ParamsB1_pinned, :ParamsB1_pinned²)
	@eval begin
		function $Params_constr(N::Int64, m::Union{Int64, Nothing}, W::Float64, a::Float64, psamples_raw::Vector{NTuple{3, Float64}})
			u₀ = get_u₀(W, nothing, psamples_raw)
			u₁ = u₀*a
			return $Params_constr(N + 1 - rem(N, 2), m, W, u₀, u₁)
		end
	end
end

for Params_constr in (:ParamsNoB1, :ParamsNoB1_pinned)
	@eval begin
		function $Params_constr(N::Int64, m::Union{Int64, Nothing}, W::Float64, psamples_raw::Vector{NTuple{3, Float64}})
			u₀ = get_u₀(W, nothing, psamples_raw)
			return $Params_constr(N, m, W, u₀)
			#return $Params_constr(N + 1 - rem(N, 2), m, W, u₀)
		end
	end
end

########

kw(W) = find_zero(k -> (W-4K(k^2)), (0.0, 1.0))

get_length(params::AbstractParamsB1) = 2params.N
get_length(params::AbstractParamsNoB1) = params.N

get_τs(params::AbstractParams) = collect(0.5:params.N-0.5)*params.Δτ

bs_init_func(τs, k::AbstractFloat) = k*Jacobi.sn.(τs, k^2)
bs_init_func(τs, k::Nothing) = tanh.(τs)
function get_bs(params::AbstractParams)
    W = params.W
    if W < 40.0
        return get_bs(params, kw(W))
    end
    τs = get_τs(params)
    try
        k = kw(W)
        return _fill_bs(τs, W,  k)
    catch e
        return _fill_bs(τs, W,  nothing)
    end
end
function get_bs(params::AbstractParams, k)
    bs = zeros(get_length(params))
    bs[1:params.N] .= bs_init_func(get_τs(params), k)
    pin_bs!(bs, params)
    return bs
end
function _fill_bs(τs, W, k)
    N = length(τs)
    @assert N>10
    Nhalf = div(N, 2) + 1
    Nquarter = div(Nhalf, 2) + 1
    bs = similar(τs)
    bs[1:Nquarter] = bs_init_func(τs[1:Nquarter], k)
    bs[Nquarter+1:Nhalf] .= bs_init_func(W/2 .- τs[Nquarter+1:Nhalf], k)
    bs[end:-1:Nhalf+1] .= -bs[2:N-Nhalf+1]
    return bs
end

########

get_pinned_idxs(params::Union{ParamsNoB1, ParamsB1}) = Int64[]
get_pinned_idxs(params::Union{ParamsB1_pinned, ParamsNoB1_pinned}) = [1, params.i_pinned]
function get_pinned_idxs(params::Union{ParamsB1_pinned²})
    N, i_pinned = params.N, params.i_pinned
    return [1, i_pinned, N+div(N,4)+1, N+div(N,4)+i_pinned]
end

get_free_idxs(params::AbstractParams) = [i for i=Base.OneTo(get_length(params)) if !(i in get_pinned_idxs(params))]

@inline pin_bs!(bs::AbstractVector, params::AbstractParams) = pin_bs!(bs, get_pinned_idxs(params))
@inline pin_bs!(bs::AbstractVector, params::Union{ParamsB1, ParamsNoB1}) = bs
function pin_bs!(bs::AbstractVector{T}, pinned_idxs::Vector{Int64}) where {T}
    for i in pinned_idxs
        bs[i] = zero(T)
    end
    return bs
end

@inline project_onto(a, params::AbstractParams) = project_onto(a, get_free_idxs(params))
@inline project_onto(a, params::Union{ParamsB1, ParamsNoB1}) = a
@inline project_onto(a::AbstractVector, free_idxs::Vector{Int64}) = view(a, free_idxs)
@inline project_onto(a::AbstractMatrix, free_idxs::Vector{Int64}) = view(a, free_idxs)

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

abstract type AbstractOptAlg end

struct NewtonOptAlg <: AbstractOptAlg end
struct BFGSOptAlg <: AbstractOptAlg end

generate_cash(::NewtonOptAlg, params::AbstractParams) = GH_Cash(params)
generate_cash(::BFGSOptAlg, params::AbstractParams) = G_Cash(params)

struct GH_Cash <: AbstractSharedCash
	∂Fs::Vector{SharedVector{Float64}}
	∂²Fs::Vector{SharedMatrix{Float64}}
	function GH_Cash(params::AbstractParams)
		np = nprocs()
		N′ = get_length(params)
		∂Fs = Vector{SharedVector{Float64}}()
		∂²Fs = Vector{SharedMatrix{Float64}}()
        ∂Fs = [@fetchfrom pid SharedArrays.SharedVector{Float64}(N′; pids = procs()) for pid=procs()]
        ∂²Fs = [@fetchfrom pid SharedArrays.SharedMatrix{Float64}(N′, N′; pids = procs()) for pid=procs()]
		new(∂Fs, ∂²Fs)
	end
end

mutable struct G_Cash <: AbstractSharedCash
	∂Fs::Vector{SharedVector{Float64}}
	∂Us::DArray{SMatrix{2,2,Complex{Float128}}, 1, Vector{SMatrix{2,2,Complex{Float128}}}}
	U_cashes::DArray{SMatrix{2,2,Complex{Float128}}, 1, Vector{SMatrix{2,2,Complex{Float128}}}}
	function G_Cash(params::AbstractParams)
		np = nprocs()
		N = params.N
		N′ = get_length(params)
        ∂Fs = [SharedVector{Float64}(N′) for pid=procs()]
		∂Us = DArray([@spawnat i LinearAlgebra.ones(StaticArrays.SMatrix{2,2, Complex{Float128}}, N′) for i=procs()])
		U_cashes = DArray([@spawnat i LinearAlgebra.ones(StaticArrays.SMatrix{2,2, Complex{Float128}}, N+1) for i=procs()])
		g_cash = new(∂Fs, ∂Us, U_cashes)
		finalizer(free_darrays, g_cash)
		return g_cash
	end
end

function free_darrays(g_cash::G_Cash)
	close(g_cash.∂Us)
	close(g_cash.U_cashes)
end

mutable struct G_Cash2 <: AbstractSharedCash
    ΔFss::DArray{Float64, 1, Vector{Float64}}
    ∂Fss::DArray{Vector{Float64}, 1, Vector{Vector{Float64}}}
	∂Us::DArray{SMatrix{2,2,Complex{Float128}}, 1, Vector{SMatrix{2,2,Complex{Float128}}}}
	U_cashes::DArray{SMatrix{2,2,Complex{Float128}}, 1, Vector{SMatrix{2,2,Complex{Float128}}}}
    function G_Cash2(params::AbstractParams, psamples)
		np = nprocs()
		N = params.N
		N′ = get_length(params)
        ΔFss = DArray([@spawnat i LinearAlgebra.zeros(Float64, length(localpart(psamples))) for i=procs()])
        ∂Fss = DArray([@spawnat i [LinearAlgebra.zeros(Float64, N′) for j=1:length(localpart(psamples))] for i=procs()])
		∂Us = DArray([@spawnat i LinearAlgebra.ones(StaticArrays.SMatrix{2,2, Complex{Float128}}, N′) for i=procs()])
		U_cashes = DArray([@spawnat i LinearAlgebra.ones(StaticArrays.SMatrix{2,2, Complex{Float128}}, N+1) for i=procs()])
		g_cash = new(ΔFss, ∂Fss, ∂Us, U_cashes)
		finalizer(free_darrays, g_cash)
		return g_cash
    end
end

function free_darrays(g_cash::G_Cash2)
	close(g_cash.ΔFss)
	close(g_cash.∂Fss)
	close(g_cash.∂Us)
	close(g_cash.U_cashes)
end

export G_Cash2

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

struct ReducedDispersion <: AbstractDispersion
    α::Float64
    P::Float64
	μ::Float64
    Λ::Float64
end

wfunc(y::Float64) = y==0.0 ? 0.0 : y^2*log((1+sqrt(1-y^6))/abs(y)^3)
@inline wfunc_gauss(y::Real) = y==0.0 ? 0.0 : log((1+sqrt(1-y^2))/abs(y))
@inline wfunc_supplementary(y::Real) = -1/y*wfunc_gauss(y)

function get_quadrature(rdisp::ReducedDispersion, rtol::Float64 = 1e-6, limit::Int64 = 200)
	α, P, μ, Λ = rdisp.α, rdisp.P, rdisp.μ, rdisp.Λ
	f(x) = 1/x*wfunc_gauss((x-P)/Λ)
	if abs(P)<Λ
		I, E, segs = quadgk_cauchy(f, P-Λ, P+Λ, 0.0, P; rtol=rtol, limit=limit)
	else
		I, E, segs = quadgk_custom(f, P-Λ, P+Λ, P; rtol = rtol, limit=limit)
	end
	x₁, w₁, x₂, w₂ = construct_quadrature(segs)
    w₂ .*= wfunc_gauss.((x₂.-P)./Λ)./(α*2*π^2)
	return x₂, w₂
end

function get_psamples(rdisp::ReducedDispersion, rtol::Float64 = 1e-6, limit::Int64 = 200)
	ε⁻s, ws = get_quadrature(rdisp, rtol, limit)
	μ = rdisp.μ
	return [(μ, ε⁻s[i], ws[i]) for i in eachindex(ws)]
end

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

(rdisp::ReducedDispersion)(y::Float64, dy::Float64) = (-rdisp.μ, rdisp.P + rdisp.Λ*y^3, wfunc(y)*6dy*rdisp.Λ/(rdisp.α*(2π)^2))::NTuple{3,Float64}
apply_rdisp(rdisp::ReducedDispersion, y::Float64, w::Float64) = (-rdisp.μ, rdisp.P + rdisp.Λ*y^3, 6w*rdisp.Λ/(rdisp.α*(2π)^2))::NTuple{3,Float64}

function get_psamples_new(rdisp::ReducedDispersion, Nₚ)
	ys, ws = gauss(wfunc, Nₚ, -1.0, 1.0)
	return [apply_rdisp(rdisp, ys[i], ws[i]) for i = 1:Nₚ]
end

function get_psamples_old(rdisp::ReducedDispersion, Nₚ)
    psamples_raw = rdisp.(range(-1.0,1.0; length = Nₚ), 2.0/(Nₚ-1))
	psamples_raw[1] = (psamples_raw[1][1], psamples_raw[1][2], 0.5*psamples_raw[1][3])
	psamples_raw[end] = (psamples_raw[end][1], psamples_raw[end][2], 0.5*psamples_raw[end][3])
	return psamples_raw
end

function get_psamples(rdisp::ReducedDispersion, Nₚ)
	temp_quad(f, _a, _b; kwargs...) = quadgk(f, -1.0, 0.0, 1.0; kwargs...)
	ys, ws = gauss(wfunc_gauss, Nₚ, -1.0, 1.0; quad=temp_quad)
	return [(-rdisp.μ, rdisp.P + rdisp.Λ*ys[i], 2*rdisp.Λ*ws[i]/(rdisp.α*(2pi)^2)) for i in eachindex(ws)]
end

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

widen(psamples_raw, W, m) = widen(psamples_raw, W, m, one(Float128))
function widen(psamples_raw::Vector{NTuple{3, Float64}}, W::Float64, m::Int64, γ::Float128)
	β′ = Float128(W)*m
	psamples_widened = Vector{Tuple{Float64, Float64, Float64, Float128, Float128}}()
	for psample in psamples_raw
		εₚ⁺′ = Float128(psample[1])
		κₚ⁰ = sqrt(γ^2 + psample[2]^2)
        a = exp(-β′*(κₚ⁰ - εₚ⁺′)) + exp(-β′*(κₚ⁰ + εₚ⁺′))
        b = log(one(Float128) + exp(-2β′*κₚ⁰) + a)
		push!(psamples_widened, (psample..., a, b))
	end
	return psamples_widened
end

function widen(psamples_raw::Vector{NTuple{3, Float64}}, W::Float64, m::Nothing, γ::Float128)
	psamples_widened = Vector{Tuple{Float64, Float64, Float64, Float128, Float128}}()
	for psample in psamples_raw
		εₚ⁺′ = Float128(psample[1])
		κₚ⁰ = sqrt(γ^2 + psample[2]^2)
		a = W*(εₚ⁺′ - κₚ⁰)
		b = max(a, zero(Float128))
		push!(psamples_widened, (psample..., a, b))
	end
	return psamples_widened
end

function separate_psamples_old(psamples_raw::Vector{Tuple{Float64,Float64,Float64, Float128, Float128}})
	np = nprocs()
	Nₚ = length(psamples_raw)
	chunk_size, num_extra = divrem(Nₚ, np)
	psamples = Vector{Vector{eltype(psamples_raw)}}(undef, np)
	for cid = 0:(num_extra-1)
		psamples[cid+1] = psamples_raw[((chunk_size+1)*cid+1):(chunk_size+1)*(cid+1)]
	end
	n₀ = (chunk_size+1)*num_extra
	for cid = 0:(np-num_extra-1)
		psamples[num_extra+cid+1] = psamples_raw[(n₀+chunk_size*cid+1):min(Nₚ, n₀+chunk_size*(cid+1))]
	end
	return psamples
end

separate_psamples(psamples_raw::AbstractVector) = distribute(psamples_raw; procs = procs())

function generate_psamples(rdisp::ReducedDispersion, Nₚ)
	return separate_psamples(get_psamples(rdisp, Nₚ))
end

struct Dispersion <:AbstractDispersion
	α::Float64
	α′::Float64
	μ::Float64
	Λ::Float64
	P::Float64
end
Dispersion(α, α′, μ, Λ) = Dispersion(α, α′, μ, Λ, 0.0)

wfunc_φ(φ) = 1/sqrt(1-φ^2)

macro unpack(var, args...)
	exs = [:($arg = $(esc(var)).$arg) for arg in args]
    return Expr(:block, exs...)
end

function get_psamples_bad(disp::Dispersion, Nₜ, Nᵩ)
	#@unpack disp α α′ μ Λ P
	α, α′, μ, Λ, P = disp.α, disp.α′, disp.μ, disp.Λ, disp.P
	φs, wᵩs = gauss(wfunc_φ, Nᵩ, -1.0, 1.0)
	ts, wₜs = gauss(Nₜ, 0.0, 1.0)
	c1 = 4*Λ/((α + α′)*(2π)^2)
	c2 = (α - α′)/(α + α′)*Λ
	psamples_raw = Vector{NTuple{3, Float64}}()
	for iₜ = 1:length(ts)
		t = ts[iₜ]
		wₜ = wₜs[iₜ]
		εₚ⁺ = c2*t - μ
		for iᵩ = 1:length(φs)
			wₚ = c1*wₜ*wᵩs[iᵩ]
			εₚ⁻ = Λ*t*φs[iᵩ] + P
			push!(psamples_raw, (εₚ⁺, εₚ⁻, wₚ))
		end
	end
	return psamples_raw
end

function get_psamples_simple(disp::Dispersion, Nₜ, Nᵩ)
	#@unpack disp α α′ μ Λ P
	α, α′, μ, Λ, P = disp.α, disp.α′, disp.μ, disp.Λ, disp.P
	φs, wᵩs = [(i-0.5)*π/Nᵩ for i=1:Nᵩ], fill!(zeros(Nᵩ), π/Nᵩ)
	ts, wₜs = gauss(Nₜ, 0.0, 1.0)
	c1 = 4*Λ/((α + α′)*(2π)^2)
	c2 = (α - α′)/(α + α′)*Λ
	psamples_raw = Vector{NTuple{3, Float64}}()
	for iₜ = 1:Nₜ
		t = ts[iₜ]
		wₜ = wₜs[iₜ]
		εₚ⁺ = c2*t - μ
		for iᵩ = 1:Nᵩ
			wₚ = c1*wₜ*wᵩs[iᵩ]
			εₚ⁻ = Λ*t*cos(φs[iᵩ]) + P
			push!(psamples_raw, (εₚ⁺, εₚ⁻, wₚ))
		end
	end
	return psamples_raw
end

get_psamples(disp::Dispersion, Nₜ, Nᵩ) = get_psamples_simple(disp, Nₜ, Nᵩ)

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

function get_static_energy(a₂, rdisp, int_rtol::Real=1e-6, limits::Int64 = 200)
	psamples_raw = get_psamples(rdisp, int_rtol, limits)
    u′ = get_u₀(psamples_raw)
    u′′ = u′*(1+a₂)
    g = find_zero(x -> (u′′ - get_u₀(psamples_raw, x)), (1e-20, 1e20))
    s = 0.0
	for (εₚ⁺, εₚ⁻, wₚ) in psamples_raw
        s += wₚ*(abs(εₚ⁻) - sqrt(εₚ⁻^2 + g^2))
    end
    s += g^2*u′′
end

export get_static_energy

θₕ(x) = x>0.0 ? 1.0 : 0.0

function get_u₀_quad(rdisp, γ²)
    α, P, μₐ, Λ = rdisp.α, rdisp.P, abs(rdisp.μ), rdisp.Λ
    f(x) = θₕ(sqrt(x^2 + γ²) - μₐ)*wfunc_gauss((x-P)/Λ)/(sqrt(x^2 + γ²)*4α*π^2)
    return quadgk_custom(f, P-Λ, P, P+Λ; rtol = 1e-6)[1]
end
function get_u₀_quad(rdisp, γ², β)
    α, P, μₐ, Λ = rdisp.α, rdisp.P, abs(rdisp.μ), rdisp.Λ
	temp_func(x, μₐ, γ², β) = 0.5*(tanh(0.5β*(μₐ + sqrt(x^2+γ²))) + tanh(0.5β*(sqrt(x^2 + γ²) - μₐ)))
	f(x) = θₕ(sqrt(x^2 + γ²) - μₐ)*wfunc_gauss((x-P)/Λ)*temp_func(x, μₐ, γ², β)/(sqrt(x^2 + γ²)*4α*π^2)
    return quadgk_custom(f, P-Λ, P, P+Λ; rtol = 1e-6)[1]
end
function calc_se_quad(rdisp, γ², u)
    α, P, μₐ, Λ = rdisp.α, rdisp.P, abs(rdisp.μ), rdisp.Λ
    f(x) = wfunc_gauss((x-P)/Λ)*(max(μₐ, abs(x)) - max(μₐ, sqrt(x^2 + γ²)))/(2α*π^2)
    return quadgk(f, P-Λ, P, P+Λ)[1] + u*γ²
end
function calc_se_quad(rdisp, γ², u, β)
	β′ = BigFloat(β)
    α, P, μₐ, Λ = rdisp.α, rdisp.P, abs(rdisp.μ), rdisp.Λ
	energy_func(x, μₐ, γ², β′) = Float64(log((cosh(β′*μₐ) + cosh(β′*x))/(cosh(β′*μₐ) + cosh(β′*sqrt(x^2 + γ²))))/β′)
    f(x) = wfunc_gauss((x-P)/Λ)*energy_func(x, μₐ, γ², β′)/(2α*π^2)
    return quadgk(f, P-Λ, 0.0, P, P+Λ)[1] + u*γ²
end
function get_static_energy_quad(a₂, rdisp, interval::NTuple{2,Float64} = (1e-30, 1e30))
    u = get_u₀_quad(rdisp, 1.0)
    u′ = u*(1+a₂)
    #print(u, "\n",u′, "\n")
    f(x) = (u′ - get_u₀_quad(rdisp, x))
    #print(f(interval[1]),"\n", f(interval[2]), "\n\n")
    γ² = find_zero(f, interval)
    return calc_se_quad(rdisp, γ², u′), γ²
end

get_static_energy_quad(a₂, rdisp, β::Nothing, interval::NTuple{2,Float64} = (1e-30,1e30)) = get_static_energy_quad(a₂, rdisp, interval)
function get_static_energy_quad(a₂, rdisp, β::AbstractFloat, interval::NTuple{2,Float64} = (1e-40,1e40))
	u = get_u₀_quad(rdisp, 1.0)
	u′ = u*(1+a₂)
	γ² = find_zero(x -> (u′ - get_u₀_quad(rdisp, x, β)), interval)
	return calc_se_quad(rdisp, γ², u′, β), γ²
end

function get_u₀_quad0(rdisp, β)
    α, P, μₐ, Λ = rdisp.α, rdisp.P, abs(rdisp.μ), rdisp.Λ
	temp_func(x, μₐ, β) = 0.5*(tanh(0.5β*(μₐ + abs(x))) + tanh(0.5β*(abs(x) - μₐ)))
	f(x) = θₕ(abs(x) - μₐ)*wfunc_gauss((x-P)/Λ)*temp_func(x, μₐ, β)/(abs(x)*4α*π^2)
    return quadgk_custom(f, P-Λ, P, P+Λ; rtol = 1e-6)[1]
end
function get_cutoff_β(a₂, rdisp)
	u = get_u₀_quad(rdisp, 1.0)
	u′ = u*(1+a₂)
	β = find_zero(x -> (u′ - get_u₀_quad0(rdisp, x)), (1e-20, 1000.0))
	return β
end
export get_static_energy_quad, get_cutoff_β


text:
"""
This functions takes two strings and calculates de edit distance, which is the number
of changes of two stringsn when both display different lengths:

```@example
    s1 = "ATCTCGT"
    s2 = "ACTCCTC"
    edit_distance(s1,s2)[1] = 3.0
    edit_distance(s1,s2)[2] =  [0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0;
                               1.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0;
                               2.0 1.0 1.0 1.0 2.0 3.0 4.0 5.0;
                               3.0 2.0 1.0 2.0 1.0 2.0 3.0 4.0;
                               4.0 3.0 2.0 1.0 2.0 2.0 2.0 3.0;
                               5.0 4.0 3.0 2.0 1.0 2.0 3.0 2.0;
                               6.0 5.0 4.0 3.0 2.0 2.0 3.0 3.0;
                               7.0 6.0 5.0 4.0 3.0 3.0 2.0 3.0]
```
This functions outputs the edit distance and the alignment matrix as well.
"""
function edit_distance(s1::String, s2::String)
    D::Array = zeros(length(s1)+1,length(s2)+1)
    D[:,1]= 0:length(s1)
    D[1,:] = 0:length(s2)
    for i in 2:length(s1)+1
        for j in 2:length(s2)+1
            (s1[i-1] != s2[j-1]) ? δ = 1 : δ = 0
            D[i,j] = min(D[i-1,j-1]+δ, D[i-1,j]+1, D[i,j-1]+1)
        end
    end
    return D[length(s1)+1, length(s2)+1],D
end


text:
@testset "neighbor features" begin
    scene=Scene(Entity{VehicleState, BicycleModel, Int}, 100)
    roadway=gen_straight_roadway(3, 200.0, lane_width=3.0)
    push!(scene,Entity(VehicleState(VecSE2(30.0,3.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),1))
    push!(scene,Entity(VehicleState(VecSE2(20.0,3.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),2))
    push!(scene,Entity(VehicleState(VecSE2(40.0,3.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),3))
    push!(scene,Entity(VehicleState(VecSE2(20.0,0.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),4))
    push!(scene,Entity(VehicleState(VecSE2(40.0,0.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),5))
    push!(scene,Entity(VehicleState(VecSE2(20.0,6.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),6))
    push!(scene,Entity(VehicleState(VecSE2(40.0,6.0,0.0), roadway, 0.0),
        BicycleModel(VehicleDef(AgentClass.CAR, 4.826, 1.81)),7))
    @test find_neighbor(scene, roadway, scene[1]) == NeighborLongitudinalResult(3,10.0)
    @test find_neighbor(scene, roadway, scene[1], rear=true) == NeighborLongitudinalResult(2,10.0)
    @test find_neighbor(scene, roadway, scene[1], lane=leftlane(roadway, scene[1])) == NeighborLongitudinalResult(7,10.0)
    @test find_neighbor(scene, roadway, scene[1], lane=leftlane(roadway, scene[1]), rear=true) == NeighborLongitudinalResult(6,10.0)
    @test find_neighbor(scene, roadway, scene[1], lane=rightlane(roadway, scene[1])) == NeighborLongitudinalResult(5,10.0)
    @test find_neighbor(scene, roadway, scene[1], lane=rightlane(roadway, scene[1]), rear=true) == NeighborLongitudinalResult(4,10.0)

    trajdata = get_test_trajdata(roadway)
    scene = trajdata[1]
    @test find_neighbor(scene, roadway, scene[1]) == NeighborLongitudinalResult(2, 3.0)
    @test find_neighbor(scene, roadway, scene[2]) == NeighborLongitudinalResult(nothing, 250.0)
    scene = trajdata[2]
    @test find_neighbor(scene, roadway, scene[1]) == NeighborLongitudinalResult(2, 4.0)
    @test find_neighbor(scene, roadway, scene[2]) == NeighborLongitudinalResult(nothing, 250.0)

    roadway = gen_stadium_roadway(1)
    scene = Scene(Entity{VehicleState, VehicleDef, Int64}, 2)
    scene.n = 2
    def = VehicleDef(AgentClass.CAR, 2.0, 1.0)

    function place_at!(i, s)
        roadproj = proj(VecSE2(0.0,0.0,0.0), roadway)
        roadind = RoadIndex(roadproj.curveproj.ind, roadproj.tag)
        roadind = move_along(roadind, roadway, s)
        frenet = Frenet(roadind, roadway[roadind].s, 0.0, 0.0)
        state = VehicleState(frenet, roadway, 0.0)
        scene[i] = Entity(state, def, i)
    end

    place_at!(1, 0.0)
    place_at!(2, 0.0)
    foreinfo = find_neighbor(scene, roadway, scene[1], max_distance=Inf)
    @test isapprox(foreinfo.Δs, 0.0)

    place_at!(2, 100.0)
    foreinfo = find_neighbor(scene, roadway, scene[1], max_distance=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 100.0)

    foreinfo = find_neighbor(scene, roadway, scene[2], max_distance=Inf)
    @test foreinfo.ind == 1
    @test isapprox(foreinfo.Δs, 277.07, atol=1e-2)

    place_at!(2, 145.0)
    foreinfo = find_neighbor(scene, roadway, scene[1], max_distance=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 145.0, atol=1e-5)

    place_at!(2, 240.0)
    foreinfo = find_neighbor(scene, roadway, scene[1], max_distance=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 240.0, atol=1e-5)

    place_at!(1, 240.0)
    foreinfo = find_neighbor(scene, roadway, scene[1], max_distance=Inf)
    @test foreinfo.ind == 2
    @test isapprox(foreinfo.Δs, 0.0, atol=1e-5)

    ################################################
    # get_frenet_relative_position

    place_at!(1, 0.0)
    place_at!(2, 10.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.origin == scene[1].state.posF.roadind
    @test frp.target.ind.i == 1
    @test isapprox(frp.target.ind.t, 0.1, atol=0.01)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 10.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    place_at!(1, 10.0)
    place_at!(2, 20.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.origin == scene[1].state.posF.roadind
    @test frp.target.ind.i == 1
    @test isapprox(frp.target.ind.t, 0.2, atol=0.01)
    @test frp.target.tag == scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 10.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    place_at!(1, 0.0)
    place_at!(2, 120.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.target.tag == LaneTag(1, 1)
    @test scene[1].state.posF.roadind.tag == LaneTag(6, 1)
    @test isapprox(frp.Δs, 120.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)

    place_at!(1, 0.0)
    place_at!(2, 250.0)
    frp = get_frenet_relative_position(scene[2].state.posG, scene[1].state.posF.roadind, roadway, max_distance_fore=Inf)
    @test frp.target.tag != scene[1].state.posF.roadind.tag
    @test isapprox(frp.Δs, 250.0, atol=1e-2)
    @test isapprox(frp.t, 0.0, atol=1e-5)
    @test isapprox(frp.ϕ, 0.0, atol=1e-7)
end

@testset "feature extraction" begin
    roadway = gen_straight_roadway(4, 100.0)

    scene = Scene([Entity(VehicleState(VecSE2( 0.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
        Entity(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
            ])

    # test each feature individually
    pos1 = extract_feature(featuretype(posgx), posgx, roadway, [scene], 1)
    pos2 = extract_feature(featuretype(posgx), posgx, roadway, [scene], 2)
    poss = extract_feature(featuretype(posgx), posgx, roadway, [scene], [1,2])
    @test poss[1][1] == pos1[1]
    @test poss[2][1] == pos2[2]
    @test pos1[1] == 0.0
    @test pos2[2] == 10.0

    scene = Scene([Entity(VehicleState(VecSE2(1.1,1.2,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1)])
    posy = extract_feature(featuretype(posgy), posgy, roadway, [scene], 1)
    @test posy[1] == 1.2
    posθ = extract_feature(featuretype(posgθ), posgθ, roadway, [scene], 1)
    @test posθ[1] == 0.0
    poss = extract_feature(featuretype(posfs), posfs, roadway, [scene], 1)
    @test poss[1] == 1.1
    post = extract_feature(featuretype(posft), posft, roadway, [scene], 1)
    @test post[1] == 1.2
    posϕ = extract_feature(featuretype(posfϕ), posfϕ, roadway, [scene], 1)
    @test posϕ[1] == 0.0

    scene = Scene([Entity(VehicleState(VecSE2(1.1,1.2,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
                Entity(VehicleState(VecSE2(1.5,1.2,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2)])
    coll = extract_feature(featuretype(iscolliding), iscolliding, roadway, [scene], 1)
    @test coll[1]

    d = extract_feature(featuretype(distance_to(1)), distance_to(1), roadway, [scene], 1)
    @test d[1] == 0.0
    d = extract_feature(featuretype(distance_to(2)), distance_to(2), roadway, [scene], 1)
    @test d[1] ≈ 0.4

    df = extract_feature(featuretype(turn_rate_g), turn_rate_g, roadway, [scene, scene], [1])
    @test df[1][1] === missing
    @test df[1][2] == 0.0

    # extract multiple features
    roadway = gen_straight_roadway(3, 1000.0, lane_width=1.0)
    scene = Scene([
                Entity(VehicleState(VecSE2( 0.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
                Entity(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
            ])

    dfs = extract_features((iscolliding, markerdist_left, markerdist_right), roadway, [scene], [1,2])
    @test isapprox(dfs[1][1,3], 0.5)
    @test isapprox(dfs[2][1,3], 0.5)
    @test isapprox(dfs[1][1,2], 0.5)
    @test isapprox(dfs[2][1,2], 0.5)

    # integration testing, all the features
    feature_list = (posgx, posgy, posgθ, posfs, posft, posfϕ, vel, velfs, velft, velgx, velgy,
                    time_to_crossing_right, time_to_crossing_left,
                    estimated_time_to_lane_crossing, iswaiting, acc, accfs, accft, jerk,
                    jerkft, turn_rate_g, turn_rate_f, isbraking, isaccelerating,
                    lane_width, lane_offset_left, lane_offset_right, has_lane_left,
                    has_lane_right, lane_curvature)
    dfs = extract_features(feature_list, roadway, [scene, scene], [1,2])
    for id=[1,2]
        @test ncol(dfs[id]) == length(feature_list)
        @test nrow(dfs[id]) == 2
    end

    scene = Scene([
            Entity(VehicleState(VecSE2( 1.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Entity(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
            Entity(VehicleState(VecSE2(12.0,1.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 3),
            Entity(VehicleState(VecSE2( 0.0,1.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 4),
        ])
    dfs= extract_features((dist_to_front_neighbor, front_neighbor_speed, time_to_collision), roadway, [scene], [1,2,3,4])
    @test isapprox(dfs[1][1,1], 9.0)

    # type inferrence
    for fun in feature_list
        @inferred featuretype(fun)
    end

end # features

@testset "lidar sensor" begin
    #=
    Create a scene with 7 vehicles, an ego in the middle and 6 surrounding
    vehicles. Use this setup to check that lidar sensors work for various
    range and angle settings.
    =#
    num_veh = 7
    ego_index = 1
    roadway = gen_straight_roadway(4, 400.)
    scene = Scene(Entity{VehicleState,VehicleDef,Int64}, num_veh)
    # order: ego, fore, rear, left, right, fore_fore, fore_fore_fore
    speeds = [10., 15., 15., 0., -5, 20., 20.]
    positions = [200., 220., 150., 200., 200., 240., 350.]
    lanes = [2,2,2,4,1,2,2]
    for i in 1:num_veh
        lane = roadway.segments[1].lanes[lanes[i]]
        road_idx = RoadIndex(proj(VecSE2(0.0, 0.0, 0.0), lane, roadway))
        veh_state = VehicleState(Frenet(road_idx, roadway), roadway, speeds[i])
        veh_state = move_along(veh_state, roadway, positions[i])
        veh_def = VehicleDef(AgentClass.CAR, 2., 2.)
        push!(scene, Entity(veh_state, veh_def, i))
    end

    # basic lidar with sufficient range for all vehicles
    nbeams = 4
    lidar = LidarSensor(nbeams, max_range=200., angle_offset=0.)
    observe!(lidar, scene, roadway, ego_index)
    # order: right, fore, left, back
    @test isapprox(lidar.ranges[1], 2.0, atol=4)
    @test isapprox(lidar.ranges[2], 19.0, atol=4)
    @test isapprox(lidar.ranges[3], 5.0, atol=4)
    @test isapprox(lidar.ranges[4], 49.0, atol=4)
    @test isapprox(lidar.range_rates[1], 0.0, atol=4)
    @test isapprox(lidar.range_rates[2], 5.0, atol=4)
    @test isapprox(lidar.range_rates[3], 0.0, atol=4)
    @test isapprox(lidar.range_rates[4], -5.0, atol=4)

    # angles are such that all the vehicles are missed
    nbeams = 4
    lidar = LidarSensor(nbeams, max_range=200., angle_offset=π/4)
    observe!(lidar, scene, roadway, ego_index)
    # order: right, fore, left, back
    @test isapprox(lidar.ranges[1], 200.0, atol=4)
    @test isapprox(lidar.ranges[2], 200.0, atol=4)
    @test isapprox(lidar.ranges[3], 200.0, atol=4)
    @test isapprox(lidar.ranges[4], 200.0, atol=4)
    @test isapprox(lidar.range_rates[1], 0.0, atol=4)
    @test isapprox(lidar.range_rates[2], 0.0, atol=4)
    @test isapprox(lidar.range_rates[3], 0.0, atol=4)
    @test isapprox(lidar.range_rates[4], 0.0, atol=4)

    # range short enough that it misses the rear vehicle
    nbeams = 4
    lidar = LidarSensor(nbeams, max_range=20., angle_offset=0.)
    observe!(lidar, scene, roadway, ego_index)
    # order: right, fore, left, back
    @test isapprox(lidar.ranges[1], 2.0, atol=4)
    @test isapprox(lidar.ranges[2], 19.0, atol=4)
    @test isapprox(lidar.ranges[3], 5.0, atol=4)
    @test isapprox(lidar.ranges[4], 20.0, atol=4)
    @test isapprox(lidar.range_rates[1], 0.0, atol=4)
    @test isapprox(lidar.range_rates[2], 5.0, atol=4)
    @test isapprox(lidar.range_rates[3], 0.0, atol=4)
    @test isapprox(lidar.range_rates[4], 0.0, atol=4)

end


text:
module ComradeNested

using Comrade
using AbstractMCMC
using TupleVectors
using Reexport

@reexport using NestedSamplers

Comrade.samplertype(::Type{<:Nested}) = Comrade.IsCube()

function AbstractMCMC.sample(post::Comrade.TransformedPosterior, sampler::Nested, args...; kwargs...)
    ℓ = logdensityof(post)
    model = NestedModel(ℓ, identity)

    samples, stats = sample(model, sampler, args...; chain_type=Array, kwargs...)
    weights = samples[:, end]
    chain = transform.(Ref(post), eachrow(samples[:,1:end-1]))
    return TupleVector(chain), merge((;weights,), stats)
end

end