text:
function checkschema(heap::CachedArrays.CompactHeap, schemas)
    @test length(heap) == length(schemas)
    for (block, schema) in zip(heap, schemas)
        isfree, size = schema
        @test CachedArrays.isfree(block) == isfree
        @test block.size == size
    end
end

@testset "CompactHeap" begin
    @testset "Testing Power of Two" begin
        x = CachedArrays.poweroftwo(10)
        for i in 1:10
            @test CachedArrays.binsize(x, i) == (2^10) * i
        end
        @test CachedArrays.getbin(x, 1024) == 1
        @test CachedArrays.getbin(x, 1025) == 2
        @test CachedArrays.getbin(x, 2 * 1024 - 1) == 2
        @test CachedArrays.getbin(x, 2 * 1024 - 0) == 2
        @test CachedArrays.getbin(x, 2 * 1024 + 1) == 3
    end

    @testset "Heap Fundamentals" begin
        allocator = CachedArrays.AlignedAllocator()
        len = 2^12
        minallocation = 10
        heap = CachedArrays.CompactHeap(allocator, len; minallocation = 10)
        M = 2^minallocation
        @test length(CachedArrays.freelists(heap)) == div(len, M)

        # Perform a min-allocation
        # Heap layout after allocation should be
        #
        # +---------+-----------+-----------+-----------+
        # |   B1    |    B2     |    B3     |    B4     |
        # +---------+-----------+-----------+-----------+
        # |  ptr0   |              Free                 |
        # +---------+-----------+-----------+-----------+
        ptr0 = CachedArrays.alloc(heap, 2^minallocation - CachedArrays.headersize())
        schema = [(false, M), (true, 3 * M)]
        checkschema(heap, schema)

        # Expected State
        # +---------+-----------+-----------+-----------+
        # |   B1    |    B2     |    B3     |    B4     |
        # +---------+-----------+-----------+-----------+
        # |  ptr0   |         ptr1          |   Free    |
        # +---------+-----------+-----------+-----------+
        ptr1 = CachedArrays.alloc(heap, 2 * 2^minallocation - CachedArrays.headersize())
        schema = [(false, M), (false, 2 * M), (true, M)]
        checkschema(heap, schema)

        @test CachedArrays.canalloc(heap, 2 * M) == false
        @test CachedArrays.canalloc(heap, M - CachedArrays.headersize()) == true

        # Expected State
        # +---------+-----------+-----------+-----------+
        # |   B1    |    B2     |    B3     |    B4     |
        # +---------+-----------+-----------+-----------+
        # |  Free   |         ptr1          |   Free    |
        # +---------+-----------+-----------+-----------+
        CachedArrays.free(heap, ptr0)
        schema = [(true, M), (false, 2 * M), (true, M)]
        checkschema(heap, schema)

        # Expected State
        # +---------+-----------+-----------+-----------+
        # |   B1    |    B2     |    B3     |    B4     |
        # +---------+-----------+-----------+-----------+
        # |                   Free                      |
        # +---------+-----------+-----------+-----------+
        CachedArrays.free(heap, ptr1)
        schema = [(true, 4 * M)]
        checkschema(heap, schema)
    end

    @testset "Stress" begin
        allocator = CachedArrays.AlignedAllocator()
        len = 2^27
        heap = CachedArrays.CompactHeap(allocator, len)
        numtests = 10000
        pointers = Set{Ptr{Nothing}}()
        Random.seed!(123)

        timer = TimerOutput()
        num_outer_tests = 3
        @timeit timer "Test Running" for test in Base.OneTo(num_outer_tests)
            for _ in Base.OneTo(numtests)
                # CHange ratios of allocations to frees.
                if test == 1
                    num_allocs = rand(1:20)
                elseif test == 2
                    num_allocs = rand(10:20)
                else
                    num_allocs = rand(1:10)
                end

                for _ in 1:num_allocs
                    # Make a random allocation.
                    sz = rand(0:len - 256)
                    @timeit timer "Allocating" ptr = CachedArrays.alloc(heap, sz)

                    # Throw an error if it doesn't pass.
                    # That way, if it actually doesn't pass, we don't get SPAMMED in the
                    # test window.
                    @timeit timer "Checking 1" passed = CachedArrays.check(heap)
                    !passed && error()

                    if !isnothing(ptr)
                        @test !in(ptr, pointers)
                        push!(pointers, ptr)
                    elseif isnothing(ptr)
                        canalloc = CachedArrays.canalloc(heap, sz)
                        @test sz == 0 || !canalloc
                        canalloc || break
                    end
                end

                if test == 1
                    num_frees = rand(1:10)
                elseif test == 2
                    num_frees = rand(10:20)
                else
                    num_frees = rand(1:11)
                end

                for _ in 1:num_frees
                    isempty(pointers) && break
                    ptr = rand(pointers)
                    delete!(pointers, ptr)
                    @timeit timer "Freeing" CachedArrays.free(heap, ptr)
                    @timeit timer "Checking 2" passed = CachedArrays.check(heap)
                    !passed && error()
                end
            end
        end
        printstyled(stdout, "Timer Results for Compact Heap\n"; color = :cyan, bold = true)
        println(stdout, timer)

        for ptr in pointers
            CachedArrays.free(heap, ptr)
            @test CachedArrays.check(heap)
        end
    end

    @testset "Eviction" begin
        # Make a small-granularity heap so we can wrap our heads around what's going on.
        # This heap will only have 4 blocks, each of size 4096.
        allocator = CachedArrays.AlignedAllocator()
        heap = CachedArrays.CompactHeap(allocator, 4 * 4096; minallocation = 12)
        @test sizeof(heap) == 4 * 4096

        # Get the base pointer for verification purposes.
        base = CachedArrays.basepointer(heap)

        # The first allocation should go into the bottom-most block
        ptrA = CachedArrays.alloc(heap, 10, 1)
        @test ptrA == base + CachedArrays.headersize()

        # Next allocation should go to the next block.
        ptrB = CachedArrays.alloc(heap, 4096 + 10, 2)
        @test ptrB == base + 4096 + CachedArrays.headersize()

        # Last allocation should take the the last available slot.
        ptrC = CachedArrays.alloc(heap, 10, 3)
        @test ptrC == base + 3 * 4096 + CachedArrays.headersize()

        # By now, the cache should be full.
        @test isnothing(CachedArrays.alloc(heap, 10, 4))

        # Now, when we evice ptrA, should get a CB for id 1
        ids = Int[]
        cb = x -> push!(ids, CachedArrays.getid(x))
        CachedArrays.evictfrom!(heap, ptrA, 10; cb = cb)
        @test ids == [1]
        @test CachedArrays.check(heap)

        # Okay, get this pointer back, lets try evicting something larger.
        ptrA = CachedArrays.alloc(heap, 10, 1)
        @test ptrA == base + CachedArrays.headersize()

        empty!(ids)
        CachedArrays.evictfrom!(heap, ptrA, 4096 + 10; cb = cb)
        @test ids == [1,2]

        # Repopulate - lets try to force allocation to go backwards.
        CachedArrays.free(heap, ptrC)

        #####
        ##### New set of tests for easier copy-paste debugging
        #####

        allocator = CachedArrays.AlignedAllocator()
        heap = CachedArrays.CompactHeap(allocator, 4 * 4096; minallocation = 12)
        base = CachedArrays.basepointer(heap)
        ids = Int[]
        cb = x -> push!(ids, CachedArrays.getid(x))

        ptrA = CachedArrays.alloc(heap, 10, 1)
        ptrB = CachedArrays.alloc(heap, 10, 2)
        ptrC = CachedArrays.alloc(heap, 10, 3)
        ptrD = CachedArrays.alloc(heap, 10, 4)

        @test ptrA == base + CachedArrays.headersize()
        @test ptrB == base + 4096 + CachedArrays.headersize()
        @test ptrC == base + 2 * 4096 + CachedArrays.headersize()
        @test ptrD == base + 3 * 4096 + CachedArrays.headersize()

        empty!(ids)
        CachedArrays.evictfrom!(heap, ptrC, 2 * 4096 + 10; cb = cb)
        @test ids == [3, 4, 2]

        ptrD = CachedArrays.alloc(heap, 2 * 4096 + 10)
        @test ptrD == base + 4096 + CachedArrays.headersize()
    end

    @testset "Defragmentation" begin
        # The idea here is to get the heap into a known state with an expected
        # post-defragmentation result, then try defragmentation with an appropriate
        # callback and make sure the right things happen.
        #
        # First, we will do the easy case with no blocks queued for being returned
        # to the heap, then we'll test that case.
        minallocation = CachedArrays.PowerOfTwo(12)
        minallocation_bytes = 2 ^ (minallocation.val)
        allocation_size = 10 * 2 ^ (minallocation.val)

        heap = CachedArrays.CompactHeap(
            CachedArrays.AlignedAllocator(),
            allocation_size;
            minallocation = minallocation,
        )

        oneblock = (2 ^ (minallocation.val) - CachedArrays.headersize())
        twoblock = (2 * 2 ^ (minallocation.val) - CachedArrays.headersize())

        idmap = Dict{UInt64,CachedArrays.Block}()

        # The target state pre-defragmentation is shown below, where each block
        # in the upper region represents a region of size "minallocation"
        #
        # F = Free
        # A = Allocated
        #
        # +---+---+---+---+---+---+---+---+---+---+
        # |   | 1 |   | 2 |   |   | 3 | 4 | 5 |   | Heap Block
        # +---+---+---+---+---+---+---+---+---+---+
        # | F |   A1  | A2|   F   | A3| A4|   A5  | Logical Block
        # +---+-------+---+-------+---+---+-------+
        #         |     |               |     |
        #         |     +---------------+     |
        #         |          siblings         |
        #         +---------------------------+

        ptr_f1 = CachedArrays.alloc(heap, oneblock, UInt(0))
        ptr_a1 = CachedArrays.alloc(heap, twoblock, UInt(1))
        ptr_a2 = CachedArrays.alloc(heap, oneblock, UInt(2))

        ptr_f2 = CachedArrays.alloc(heap, twoblock, UInt(0))
        ptr_a3 = CachedArrays.alloc(heap, oneblock, UInt(3))
        ptr_a4 = CachedArrays.alloc(heap, oneblock, UInt(4))
        ptr_a5 = CachedArrays.alloc(heap, twoblock, UInt(5))

        basepointer = CachedArrays.basepointer(heap)
        headersize = CachedArrays.headersize()
        @test ptr_f1 == basepointer + headersize
        @test ptr_a1 == basepointer + minallocation_bytes + headersize
        @test ptr_a2 == basepointer + 3 * minallocation_bytes + headersize
        @test ptr_f2 == basepointer + 4 * minallocation_bytes + headersize
        @test ptr_a3 == basepointer + 6 * minallocation_bytes + headersize
        @test ptr_a4 == basepointer + 7 * minallocation_bytes + headersize
        @test ptr_a5 == basepointer + 8 * minallocation_bytes + headersize

        CachedArrays.free(heap, ptr_f1)
        CachedArrays.free(heap, ptr_f2)

        # Set siblings
        block_a1 = CachedArrays.unsafe_block(ptr_a1)
        block_a5 = CachedArrays.unsafe_block(ptr_a5)
        @test block_a1.id == 1
        @test block_a5.id == 5
        block_a1.sibling = block_a5
        block_a5.sibling = block_a1
        @test block_a1.id == 1
        @test block_a5.id == 5

        block_a2 = CachedArrays.unsafe_block(ptr_a2)
        block_a4 = CachedArrays.unsafe_block(ptr_a4)
        @test block_a2.id == 2
        @test block_a4.id == 4
        block_a2.sibling = block_a4
        block_a4.sibling = block_a2
        @test block_a2.id == 2
        @test block_a4.id == 4

        block_a3 = CachedArrays.unsafe_block(ptr_a3)
        @test block_a3.id == 3

        # Set block mapping
        for block in (block_a1, block_a2, block_a3, block_a4, block_a5)
            idmap[block.id] = block
        end

        # Time to defrag!
        count = 1
        CachedArrays.defrag!(heap) do id, newblock, _
            @test id == count
            count += 1

            @test haskey(idmap, id)
            idmap[id] = newblock
        end

        # Time to see how we did!
        offset = 0
        for i in 1:length(idmap)
            block = idmap[i]
            @test pointer(block) == CachedArrays.basepointer(heap) + offset
            offset += sizeof(block)
        end

        # Check siblings.
        @test idmap[1].sibling === idmap[5]
        @test idmap[5].sibling === idmap[1]

        @test idmap[2].sibling === idmap[4]
        @test idmap[4].sibling === idmap[2]

        # Number of blocks should be 6 - 5 allocated, one free at the end.
        @test length(heap) == 6
        @test CachedArrays.check(heap)
    end

    @testset "Defragmentation - Queueing" begin
        # The idea here is to get the heap into a known state with an expected
        # post-defragmentation result, then try defragmentation with an appropriate
        # callback and make sure the right things happen.
        #
        # First, we will do the easy case with no blocks queued for being returned
        # to the heap, then we'll test that case.
        minallocation = CachedArrays.PowerOfTwo(12)
        minallocation_bytes = 2 ^ (minallocation.val)
        allocation_size = 10 * 2 ^ (minallocation.val)

        heap = CachedArrays.CompactHeap(
            CachedArrays.AlignedAllocator(),
            allocation_size;
            minallocation = minallocation,
        )

        oneblock = (2 ^ (minallocation.val) - CachedArrays.headersize())
        twoblock = (2 * 2 ^ (minallocation.val) - CachedArrays.headersize())

        idmap = Dict{UInt64,CachedArrays.Block}()

        # The target state pre-defragmentation is shown below, where each block
        # in the upper region represents a region of size "minallocation"
        #
        # F = Free
        # A = Allocated
        #
        #               +- Marked As Queued. Cannot be moved.
        #               |
        #               v
        # +---+---+---+---+---+---+---+---+---+---+
        # |   | 1 |   | 2 |   |   | 3 | 4 | 5 |   | Heap Block
        # +---+---+---+---+---+---+---+---+---+---+
        # | F |   A1  | A2|   F   | A3| A4|   A5  | Logical Block
        # +---+-------+---+-------+---+---+-------+
        #         |     |               |     |
        #         |     +---------------+     |
        #         |          siblings         |
        #         +---------------------------+

        ptr_f1 = CachedArrays.alloc(heap, oneblock, UInt(0))
        ptr_a1 = CachedArrays.alloc(heap, twoblock, UInt(1))
        ptr_a2 = CachedArrays.alloc(heap, oneblock, UInt(2))

        ptr_f2 = CachedArrays.alloc(heap, twoblock, UInt(0))
        ptr_a3 = CachedArrays.alloc(heap, oneblock, UInt(3))
        ptr_a4 = CachedArrays.alloc(heap, oneblock, UInt(4))
        ptr_a5 = CachedArrays.alloc(heap, twoblock, UInt(5))

        basepointer = CachedArrays.basepointer(heap)
        headersize = CachedArrays.headersize()
        @test ptr_f1 == basepointer + headersize
        @test ptr_a1 == basepointer + minallocation_bytes + headersize
        @test ptr_a2 == basepointer + 3 * minallocation_bytes + headersize
        @test ptr_f2 == basepointer + 4 * minallocation_bytes + headersize
        @test ptr_a3 == basepointer + 6 * minallocation_bytes + headersize
        @test ptr_a4 == basepointer + 7 * minallocation_bytes + headersize
        @test ptr_a5 == basepointer + 8 * minallocation_bytes + headersize

        CachedArrays.free(heap, ptr_f1)
        CachedArrays.free(heap, ptr_f2)

        # Set siblings
        block_a1 = CachedArrays.unsafe_block(ptr_a1)
        block_a5 = CachedArrays.unsafe_block(ptr_a5)
        @test block_a1.id == 1
        @test block_a5.id == 5
        block_a1.sibling = block_a5
        block_a5.sibling = block_a1
        @test block_a1.id == 1
        @test block_a5.id == 5

        block_a2 = CachedArrays.unsafe_block(ptr_a2)
        block_a4 = CachedArrays.unsafe_block(ptr_a4)
        @test block_a2.id == 2
        @test block_a4.id == 4
        block_a2.sibling = block_a4
        block_a4.sibling = block_a2
        @test block_a2.id == 2
        @test block_a4.id == 4

        block_a3 = CachedArrays.unsafe_block(ptr_a3)
        @test block_a3.id == 3

        # Set block mapping
        for block in (block_a1, block_a2, block_a3, block_a4, block_a5)
            idmap[block.id] = block
        end

        # Time to defrag!
        block_a2.queued = true
        CachedArrays.defrag!(heap) do id, newblock, _
            @test haskey(idmap, id)
            idmap[id] = newblock
        end

        # Time to see how we did!
        offset = 0
        for i in 1:length(idmap)
            if i == 2
                offset += minallocation_bytes
            end

            block = idmap[i]
            @test pointer(block) == CachedArrays.basepointer(heap) + offset
            offset += sizeof(block)
        end

        # Check siblings.
        @test idmap[1].sibling === idmap[5]
        @test idmap[5].sibling === idmap[1]

        @test idmap[2].sibling === idmap[4]
        @test idmap[4].sibling === idmap[2]

        # Number of blocks should be 7
        # 1 allocated at the beginning.
        # 1 free after it.
        # 4 allocated.
        # 1 free for the remainder.
        @test length(heap) == 7
        @test CachedArrays.check(heap)
    end
end

@testset "Testing Eviction Corner Cases" begin
    # This tests the case where we have the following scenario:
    #
    #                 Evicted      Evicting
    #                    |             |
    #                    V             V
    # +------------+-------------+-----------+
    # |     A      |     B       |     C     |
    # +------------+-------------+-----------+
    #       ∧
    #       |
    # Freed during eviction callback
    #
    # In other words, block B is being evicted and during the eviction process, Block A
    # is freed.
    #
    # Even though block B may be next to block A in the heap, we want to prevent block A and
    # B from merging since that breaks the eviction logic.
    #
    # I.E. Eviction assumes that all blocks actively touched during the eviction process
    # do not belong to any freelist.

    # Small heap - enough for 6 allocations.
    allocator = CachedArrays.AlignedAllocator()
    heap = CachedArrays.CompactHeap(allocator, 6 * 4096; minallocation = 12)

    A = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 10))
    B = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 0))
    C = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 20))

    @test CachedArrays.getid(A) == 10
    @test CachedArrays.getid(B) == 0
    @test CachedArrays.getid(C) == 20

    # Callback function will free "A".
    # This will occur while "B" is being evicted.
    cb = function(block)
        if CachedArrays.getid(block) == CachedArrays.getid(C)
            CachedArrays.free(heap, CachedArrays.datapointer(A))
        end
    end

    CachedArrays.evictfrom!(heap, B, 2 * 4096; cb = cb)
    @test CachedArrays.check(heap)
    @test length(heap) == 1
end

# Lots of corner cases live in eviction.
# Particularly when early aborts are triggered.
# Here, we try to test the corner cases.
@testset "Testing More Eviction Corner Cases" begin
    # Scenarios to try.
    # 1. Block right before the eviction is freed during callback.
    # 2. Block right after current eviction frontier is freed during callback.
    # 3. Abort called during forward eviction frontier expansion. Block is NOT orphaned.
    # 4. Abort called during forward eviction frontier expansion. Block IS orphaned.
    # 5. Abort called during backward eviction frontier expansion. Block is NOT orphaned.
    # 6. Abort called during backward eviction frontier expansion. Block IS orphaned.

    #####
    ##### Case 1.
    #####

    let
        allocator = CachedArrays.AlignedAllocator()
        heap = CachedArrays.CompactHeap(allocator, 6 * 4096; minallocation = 12)

        A = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 10))
        B = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 0))
        C = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 20))

        @test CachedArrays.getid(A) == 10
        @test CachedArrays.getid(B) == 0
        @test CachedArrays.getid(C) == 20

        # Callback function will free "A".
        # This will occur while "B" is being evicted.
        cb = function(block)
            if CachedArrays.getid(block) == CachedArrays.getid(C)
                CachedArrays.free(heap, CachedArrays.datapointer(A))
            end
        end

        CachedArrays.evictfrom!(heap, B, 2 * 4096; cb = cb)
        @test CachedArrays.check(heap)
        @test length(heap) == 1
    end

    #####
    ##### Case 2.
    #####

    let
        allocator = CachedArrays.AlignedAllocator()
        heap = CachedArrays.CompactHeap(allocator, 6 * 4096; minallocation = 12)

        A = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 10))
        B = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 0))
        C = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 20))

        @test CachedArrays.getid(A) == 10
        @test CachedArrays.getid(B) == 0
        @test CachedArrays.getid(C) == 20

        # Callback function will free "C".
        # This will occur while "B" is being evicted.
        cb = function(block)
            if CachedArrays.getid(block) == CachedArrays.getid(B)
                CachedArrays.free(heap, CachedArrays.datapointer(C))
            end
        end

        CachedArrays.evictfrom!(heap, B, 2 * 4096; cb = cb)
        @test CachedArrays.check(heap)
        @test length(heap) == 2
        CachedArrays.free(heap, CachedArrays.datapointer(A))
        @test length(heap) == 1
    end

    #####
    ##### Case 3.
    #####

    let
        allocator = CachedArrays.AlignedAllocator()
        heap = CachedArrays.CompactHeap(allocator, 6 * 4096; minallocation = 12)

        A = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 10))
        B = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 0))
        C = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 20))
        D = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 30))

        @test CachedArrays.getid(A) == 10
        @test CachedArrays.getid(B) == 0
        @test CachedArrays.getid(C) == 20
        @test CachedArrays.getid(D) == 30

        # Callback function will free "C".
        # This will occur while "B" is being evicted.
        cb = function(block)
            return CachedArrays.getid(block) == CachedArrays.getid(C)
        end

        # Evict from B, abort at C.
        # C is still taken and was not freed during the eviction process.
        # Make sure C is still in the heap.
        CachedArrays.evictfrom!(heap, B, 2 * 4096; cb = cb)
        @test CachedArrays.isfree(A) == false
        @test CachedArrays.isfree(B) == true
        @test CachedArrays.isfree(C) == false
        @test CachedArrays.isfree(D) == false
        @test CachedArrays.check(heap)

        @test length(heap) == 5
    end

    #####
    ##### Case 4.
    #####

    let
        allocator = CachedArrays.AlignedAllocator()
        heap = CachedArrays.CompactHeap(allocator, 6 * 4096; minallocation = 12)

        A = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 10))
        B = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 0))
        C = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 20))
        D = CachedArrays.unsafe_block(CachedArrays.alloc(heap, 256, 30))

        @test CachedArrays.getid(A) == 10
        @test CachedArrays.getid(B) == 0
        @test CachedArrays.getid(C) == 20
        @test CachedArrays.getid(D) == 30

        # Callback function will free "C".
        # This will occur while "B" is being evicted.
        cb = function(block)
            if CachedArrays.getid(block) == 20
                CachedArrays.free(heap, CachedArrays.datapointer(C))
                @test C.orphaned == true
                return true
            end
            return false
        end

        # Evict from B, abort at C.
        # C is still freed while "B" is being evicted.
        # We then abort when trying to evict C.
        # Make sure that C gets reclaimed during the cleanup process.
        CachedArrays.evictfrom!(heap, B, 2 * 4096; cb = cb)
        @test CachedArrays.isfree(A) == false
        @test CachedArrays.isfree(B) == true
        # B was merged with C
        @test B.size == 2 * 4096
        @test CachedArrays.isfree(D) == false
        @test CachedArrays.check(heap)

        @test length(heap) == 4 # A, B, D, tail
    end

    #####
    ##### Case 6
    #####
end


text:
"""
A data type to represent a string format. Instances of this type are
generated by the [`@f_str`](@ref) macro.
"""
struct Format{F}
    str::String
    fun::F
end

Base.show(out::IO, fmt::Format) = print(out, "f\"", fmt.str, '"')

# READ ONLY
const DUMMY_BUFFER = UInt8[]

"""
    format(fmt::Fmt.Format, positionals...; keywords...)

Create a formatted string.
"""
format(fmt::Format, positionals...; keywords...) = stringify(fmt.fun(DUMMY_BUFFER, 1, positionals...; keywords...))

"""
    format!(buf::AbstractVector{UInt8}, [pos::Integer = 1,] fmt::Fmt.Format, positionals...; keywords...)

Write a formatted string to `buf` from `pos` and return the number of written bytes.

It is the caller's responsibility to allocate enough room on the buffer to write all data.
"""
format!(buf::AbstractVector{UInt8}, pos::Integer, fmt::Format, positionals...; keywords...) = fmt.fun(buf, pos, positionals...; keywords...)[2]
format!(buf::AbstractVector{UInt8}, fmt::Format, positionals...; keywords...) = fmt.fun(buf, 1, positionals...; keywords...)[2]

"""
    printf([out::IO,] fmt::Fmt.Format, positional...; keywords...)

Output a formatted string to `out` (default: `stdout`).
"""
printf(out::IO, fmt::Format, positionals...; keywords...) = print(out, format(fmt, positionals...; keywords...))
printf(fmt::Format, positionals...; keywords...) = printf(stdout, fmt, positionals...; keywords...)

"""
    @f_str fmt

Create a formatted string or a formatter object from string `fmt`.

If there is a field whose arguments are interpolated by `\$`, it creates a new
string.  Otherwise, it creates a formatting object of the [`Fmt.Format`](@ref)
type, which can be passed to the [`Fmt.format`](@ref) function or the
[`Fmt.printf`](@ref) function to create or output formatted strings.

# Examples
```jldoctest
julia> using Fmt

julia> x, y = 42, -8;

julia> f"({\$x}, {\$y})"  # all variables are interpolated
"(42, -8)"

julia> fmt = f"({x}, {y})"  # x and y are now a placeholder
f"({x}, {y})"

julia> Fmt.format(fmt; x, y)  # substitute variables
"(42, -8)"
```
"""
macro f_str(fstr) compile(unescape_string(fstr)) end

ext: jl


text:
# Simple implementation for tensordot and einsum for Julia,
# giving max compatibility for ForwardDiff and Zygote.

module TensorAnyDiff

using LinearAlgebra
import OMEinsum: EinCode
import Zygote: @adjoint

export @einmm_str, execute, contract

include("tensordot.jl")
include("einengine.jl")

end


text:
@testset "Lennard Jones" begin
    lj = LennardJones()
    @test lj.energy isa Function
    @test lj.force isa Function
end

@testset "Hard Sphere" begin
    struct HardSphere <: Discrete
        energy::Function

        HardSphere() = new((x, y, i, j, d2, u) -> u += Inf)
    end

    hs = HardSphere()
    @test hs.energy isa Function
end


text:
function grpc_result_or_error(result::T,
    status::Task,
    f::Function) where {T<:Union{<:Proto.ProtoType,<:Channel{<:Proto.ProtoType},<:Nothing}}

    if istaskdone(status)
        s = fetch(status)
        if !s.success
            throw(TypeDBClientException(s.message, gRPCServiceCallException(s.grpc_status,fetch(status).message)))
        elseif result === nothing
            throw(TypeDBClientException("something went wrong in gRPC", gRPCServiceCallException(0, "Error not defined by gRPC")))
        end
    end
    f(result)
end


text:
module DotOverload

export getMember
export setMember!
export @dotted

macro dotted(ex::Expr)

  function getExpression(ex, mem)
    if (mem != nothing)
      return Expr(:call, :getMember, ex, mem)
    end
    return ex
  end

  function treewalk(stuff)
    # return whatever is not an expression
    if (!isa(stuff, Expr))
      return stuff, nothing
    end

    # care only about binary expressions ( dot, assignment )
    if (length(stuff.args) == 2)
      h = stuff.head

      left, leftmember = nothing, nothing

      # special case (should be ignored):
      # module.func(p) = 2*p
      if (isa(stuff.args[1], Expr)
        && isa(stuff.args[2], Expr)
        && h == :(=)
        && stuff.args[1].head == :call
        && stuff.args[2].head == :block)

        # ignore for now
        left, leftmember = stuff.args[1], nothing
      else
        left, leftmember = treewalk(stuff.args[1])
      end

      right, rightmember = treewalk(stuff.args[2])

      getleft = getExpression(left, leftmember)
      right = getExpression(right, rightmember)

      # that's what we're looking for
      if ( h == :(.) )
        return getleft, right
      end

      # look for assignment
      if (leftmember != nothing)
        setleft = nothing

        if ( h == :(=) )
          setleft = right
        elseif ( h == :(+=) )
          setleft = :($getleft + $right)
        elseif ( h == :(-=) )
          setleft = :($getleft - $right)
        elseif ( h == :(*=) )
          setleft = :($getleft * $right)
        elseif ( h == :(/=) )
          setleft = :($getleft / $right)
        elseif ( h == :(\=) )
          setleft = :($getleft \ $right)
        elseif ( h == :(÷=) )
          setleft = :($getleft ÷ $right)
        elseif ( h == :(%=) )
          setleft = :($getleft % $right)
        elseif ( h == :(^=) )
          setleft = :($getleft ^ $right)
        elseif ( h == :(&=) )
          setleft = :($getleft & $right)
        elseif ( h == :(|=) )
          setleft = :($getleft | $right)
        elseif ( h == :($=) )
          setleft = :($getleft $ $right)
        elseif ( h == :(>>>=) )
          setleft = :($getleft >>> $right)
        elseif ( h == :(>>=) )
          setleft = :($getleft >> $right)
        elseif ( h == :(<<=) )
          setleft = :($getleft << $right)
        end

        # If it is an assignment expression, return to parent expression
        # and discard the assignment. Instead generate a call expression.
        if (setleft != nothing)
          return Expr(:call, :setMember!, left, leftmember, setleft), nothing
        end
      end

      # if we did not have an assignment or a dot expression
      # we will end here.
      stuff.args[1] = getleft
      stuff.args[2] = right

      return stuff, nothing
    end

    # non binary expressions -> fix all args
    for i = 1:length(stuff.args)
      ex, mem = treewalk(stuff.args[i])
      stuff.args[i] = getExpression(ex, mem)
    end

    return stuff, nothing
  end

  ex, mem = treewalk(ex)
  ex = getExpression(ex, mem)

  if (ex.head == :module)
    return Expr(:escape, Expr(:toplevel, ex))
  end

  return Expr(:escape, ex)
end

function getMember(t, m)
  return Base.getfield(t, m)
end

function setMember!(t, m, v)
  return Base.setfield!(t,m,v)
end

end # module


text:
include("lib_for_robot.jl")
function task18(r::Robot)

    putmarker!(r)

    path = go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    for side in (Nord, Ost, Sud, West)
        while isborder(r,side)==false
            move!(r,side)
            putmarker!(r)
        end
    end

    go_by_path!(r, path)
end


text:
using PowerSystems
using NLsolve
const PSY = PowerSystems

############### Data Network ########################
include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
include(joinpath(dirname(@__FILE__), "data_utils.jl"))
############### Data Network ########################
threebus_file_dir = joinpath(dirname(@__FILE__), "ThreeBusRenewable.raw")

### Case 2 Generators ###

function inv_generic_renewable(static_device, F_Flag)
    if F_Flag == 0
        return PSY.DynamicInverter(
            get_name(static_device),
            1.0, #ω_ref
            converter_regca(), #converter
            outer_control_TypeAB(), #outercontrol
            inner_ctrl_typeB(), #inner_control
            dc_source_lv(),
            no_pll(),
            filt_current(),
        )
    else
        return PSY.DynamicInverter(
            get_name(static_device),
            1.0, #ω_ref
            converter_regca(), #converter
            outer_control_TypeAB_FreqFlag(), #outercontrol
            inner_ctrl_typeB(), #inner_control
            dc_source_lv(),
            no_pll(),
            filt_current(),
        )
    end
end


text:
export LocalOperator, CorrelationFunction, Observable, expect

struct LocalOperator
    name::String
    prefactor::Float64
    op::String
end

struct CorrelationFunction
    name::String
    prefactor::Float64
    op1::String
    op2::String
end

struct Observable
    name::String
    op::MPO
end

function expect(psi::MPS, x::LocalOperator)::Vector{ComplexF64}
    return x.prefactor * ITensors.expect(psi, x.op)
end

function expect(psi::MPS, x::CorrelationFunction)::Matrix{ComplexF64}
    return x.prefactor * correlation_matrix(psi, x.op1, x.op2)
end

function expect(psi::MPS, x::Observable)::Tuple{ComplexF64,ComplexF64,ComplexF64}
    value = inner(psi, x.op, psi)
    squared = inner(x.op, psi, x.op, psi)
    variance = squared - (value^2)
    return value, squared, variance
end


text:
function countsort(a::Vector{<:Integer})
    lo, hi = extrema(a)
    b   = zeros(a)
    cnt = zeros(eltype(a), hi - lo + 1)
    for i in a cnt[i-lo+1] += 1 end
    z = 1
    for i in lo:hi
        while cnt[i-lo+1] > 0
            b[z] = i
            z += 1
            cnt[i-lo+1] -= 1
        end
    end
    return b
end

v = rand(UInt8, 20)
println("# unsorted bytes: $v\n -> sorted bytes: $(countsort(v))")
v = rand(1:2 ^ 10, 20)
println("# unsorted integers: $v\n -> sorted integers: $(countsort(v))")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
###########################################################################
#  Test 1 poligono                                                        #
###########################################################################
#  Input:                                                                 #
#  Mind --> Magnetizzazione indotta   --> A/m                             #
#  Mrem --> Magnetizzazione rimanente --> A/m                             #
#  Iind --> Inclinazione magn. ind.   --> °                               #
#  Irem --> Inclinazione magn. rem.   --> °                               #
#  Dind --> Declinazione magn. ind.   --> °                               #
#  Drem --> Declinazione magn. rem.   --> °                               #
#  Obs  --> Coordinate (x,z) relative ad ogni pto oss.
#  northxax --> orientazione profilo  --> °
###########################################################################
#  Output:                                                                #
#  tmag --> Campo magnetico scalare totale --> nT                         #
###########################################################################

using Test
using Mag2Dpoly

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

#function dotests()

#end

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

function test1()

    Jind = MagnetizVector(mod=4.9,Ideg=90.0,Ddeg=45.0)

    Jrem = MagnetizVector(mod=3.1,Ideg=45.0,Ddeg=0.0)

    northxax = 90.0

    N=101
    xzobs = hcat(LinRange(0.0,100.0,N), -1.0*ones(N))

    vertices  = [35.0 50.0;
                 65.0 50.0;
                 80.0 35.0;
                 65.0 20.0;
	         35.0 20.0;
	         20.0 35.0]

    ind1 = collect(1:6)
    bodyindices = [ind1]
    pbody = MagPolygBodies2D(bodyindices,vertices)

    ## compute total field (superposition of tmag of bodies)
    forwardtype = ["talwani","talwani_red","krav","wonbev"]

    nfwd = length(forwardtype)
    tmag = zeros(size(xzobs,1),nfwd)
    nbody = length(pbody.bo)
    for (f,fwdt) in enumerate(forwardtype)
        tmag[:,f] = tmagpolybodies2Dgen(xzobs,[Jind],[Jrem],northxax,pbody,fwdt)

        # for i=1:nbodya
        #     tmag[:,f] .+= tmagpoly2Dgen(xzobs,Jind,Jrem,northxax,pbody.bo[i],fwdt)
        # end
    end

    ##----------------------------------------------------------

    tmag_ref = [-116.16336912423309;
                -107.60475622349567;
                -98.17723975959588;
                -87.82662842366796;
                -76.49775256995866;
                -64.13484643609536;
                -50.682001858897294;
                -36.08369908036598;
                -20.285419949994587;
                -3.2343484367466506;
                15.119837140202215;
                34.82407638693401;
                55.92094064483901;
                78.44754640867708;
                102.43436783270268;
                127.90394975470824;
                154.8695254168912;
                183.3335469464661;
                213.28614208527696;
                244.70351807578447;
                277.5463434505119;
                311.758151040835;
                347.2638208530303;
                383.9682191188879;
                421.7550886836475;
                460.4863039069605;
                500.0016173506611;
                540.119031663933;
                580.6359235938513;
                621.331023429952;
                661.9673092115332;
                702.2958101588549;
                742.0602315222982;
                781.0022216070385;
                818.8670137044077;
                855.4091061413785;
                890.3976071295564;
                923.62087771655;
                954.890158626987;
                984.0419590104054;
                1010.9391031766961;
                1035.4704568268496;
                1057.5494684492162;
                1077.1117497867892;
                1094.111973515768;
                1108.5203855942902;
                1120.3192191143874;
                1129.499264418154;
                1136.0568061804884;
                1139.991090454424;
                1141.302439274996;
                1139.9910904544242;
                1136.0568061804884;
                1129.499264418154;
                1120.3192191143867;
                1108.52038559429;
                1094.111973515767;
                1077.1117497867892;
                1057.5494684492162;
                1035.4704568268494;
                1010.9391031766954;
                984.0419590104062;
                954.8901586269874;
                923.6208777165492;
                890.3976071295564;
                855.4091061413781;
                818.8670137044072;
                781.0022216070377;
                742.0602315222977;
                702.295810158854;
                661.9673092115327;
                621.3310234299516;
                580.6359235938512;
                540.1190316639322;
                500.0016173506604;
                460.48630390696064;
                421.75508868364693;
                383.9682191188882;
                347.2638208530297;
                311.7581510408341;
                277.5463434505112;
                244.70351807578382;
                213.28614208527665;
                183.33354694646584;
                154.86952541689124;
                127.90394975470798;
                102.43436783270306;
                78.4475464086772;
                55.92094064483902;
                34.82407638693425;
                15.119837140202181;
                -3.234348436746581;
                -20.285419949994342;
                -36.083699080365705;
                -50.68200185889743;
                -64.13484643609462;
                -76.49775256995939;
                -87.82662842366803;
                -98.17723975959623;
                -107.60475622349514;
                -116.16336912423309]

    ##----------------------------------------------
    same = [all(tmag[:,i].≈tmag_ref) for i=1:size(tmag,2)]

    if all(same)
        return true
    else
        for i=1:length(same)
            if same[i]==false
                println("Failed for forward type: $(forwardtype[i])")
            end
        end
        return false
    end
end

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

function test2()

    Jind = MagnetizVector(mod=4.9,Ideg=90.0,Ddeg=45.0)

    Jrem = MagnetizVector(mod=3.1,Ideg=45.0,Ddeg=0.0)

    northxax = 90.0

    N=101
    xzobs = hcat(LinRange(0.0,100.0,N), -1.0*ones(N))

    vertices  = [39.9767752377201 13.278991954030298;
	  43.21168888652307 13.402291449170173;
	  43.317652312473136 14.739290679698811;
	  43.546649812203974 13.203638738574346;
	  43.91560484766585 13.172384962659983;
	  44.14344907863333 13.366977100309308;
	  44.40792260208563 14.433477629153117;
	  44.69197601158024 13.41736322165884;
	  43.051984435135324 12.28972077523347;
	  41.64196146829342 11.82578388310931;
	  40.0 10.0;
	  51.593762919577486 12.985398228772771;
	  52.893379913843816 11.074145035188602;
	  53.03238448156978 10.749216912218833;
	  53.17394928686575 10.785849813929147;
	  53.21679397938986 10.936364994289892;
	  53.31684510201022 10.901007653562868;
	  52.164999382907524 9.566098324208415;
	  51.69760101144281 10.0;
	  60.559336241924306 25.07742274155419;
	  65.61740252955528 25.405022500973224;
	  66.6522482157798 24.330905878795054;
	  66.99370462700534 22.430809598897174;
	  67.95964200423342 20.871100638551827;
	  68.41821282760033 20.46729028637734;
	  67.67999306568265 20.187645182621033;
	  67.42746927845577 19.369846341974842;
	  65.71882888409377 17.507867891105946;
	  65.63379439530597 16.85561419292242;
	  65.52497931177504 14.899927178489119;
	  64.20104498482468 12.937540781866952;
	  62.56031603648637 11.598110231234838;
	  60.78820427692065 10.0]

    ind1 = collect(1:11)
    ind2 = collect(12:19)
    ind3 = collect(20:33)
    bodyindices = [ind1,ind2,ind3]
    pbody = MagPolygBodies2D(bodyindices,vertices)

    ## compute total field (superposition of tmag of bodies)
    forwardtype = ["talwani","talwani_red","krav","wonbev"]

    nfwd = length(forwardtype)
    tmag = zeros(size(xzobs,1),nfwd)
    nbody = length(pbody.bo)
    for (f,fwdt) in enumerate(forwardtype)
        tmag[:,f] = tmagpolybodies2Dgen(xzobs,[Jind,Jind,Jind],[Jrem,Jrem,Jrem],northxax,pbody,fwdt)
        # for i=1:nbody
        #     tmag[:,f] .+= tmagpoly2Dgen(xzobs,Jind,Jrem,northxax,pbody.bo[i],fwdt)
        # end
    end

    ##----------------------------------------------------------

    tmag_ref = [-25.18101203031008;
                -25.814111679386286;
                -26.46591188468871;
                -27.136513821069393;
                -27.82589071680914;
                -28.533854477346786;
                -29.260014742480834;
                -30.003728645158226;
                -30.764039135123586;
                -31.53959923300231;
                -32.32857897143186;
                -33.12855103994676;
                -33.93635025992559;
                -34.74790095972702;
                -35.558005092427436;
                -36.360082555206354;
                -37.14585368518767;
                -37.90495244621838;
                -38.624457626684254;
                -39.28832887474771;
                -39.87673534453047;
                -40.36526832574021;
                -40.724037390920856;
                -40.91666522359578;
                -40.89922361225488;
                -40.61919790412695;
                -40.01463677488;
                -39.013746213150505;
                -37.53532637130842;
                -35.49062356397606;
                -32.78735204013638;
                -29.336765353626436;
                -25.06459486138753;
                -19.926213598986948;
                -13.925267157081024;
                -7.133055532339174;
                0.2966799431282361;
                8.125321786345921;
                16.047468856320766;
                23.73409183413679;
                30.896898575646777;
                37.35947087475763;
                43.11364135233236;
                48.34131017064244;
                53.39281394546982;
                58.7274836881353;
                64.83295011470159;
                72.14348149663033;
                80.97541229042082;
                91.49257447124981;
                103.70790485858134;
                117.51814892303605;
                132.75663244957155;
                149.23869197241677;
                166.77366986321647;
                185.1314145984579;
                203.9747839694488;
                222.7894273624335;
                240.84727102364297;
                257.2295606565289;
                270.9152383168951;
                280.91830730494206;
                286.440275178761;
                286.99698347216184;
                282.4869194996864;
                273.1878553816839;
                259.6912319962327;
                242.79908349000044;
                223.41189102240213;
                202.4297310575936;
                180.67896175091312;
                158.86732337391078;
                137.5640994379945;
                117.19911113444167;
                98.07391651937287;
                80.37955931028789;
                64.21668987444241;
                49.615331309029074;
                36.552740278605995;
                24.96865481639862;
                14.7777648040935;
                5.879555320167113;
                -1.834170357150894;
                -8.473724523855406;
                -14.147563163141468;
                -18.959577242278417;
                -23.00736030015876;
                -26.381212717810055;
                -29.163688024210938;
                -31.429526675334152;
                -33.245857032715975;
                -34.67257173728606;
                -35.76281076807948;
                -36.5635008447628;
                -37.11591518972078;
                -37.45622870164328;
                -37.616051936809455;
                -37.6229334959912;
                -37.50082493720261;
                -37.27050556150594;
                -36.949966659975665]

    ##----------------------------------------------
    same = [all(tmag[:,i].≈tmag_ref) for i=1:size(tmag,2)]

    if all(same)
        return true
    else
        for i=1:length(same)
            if same[i]==false
                println("Failed for forward type: $(forwardtype[i])")
            end
        end
        return false
    end
end

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

#dotests()

@testset "Mag2D polygons" begin

    printstyled("Test 1 \n", bold=true,color=:cyan)
    @test test1()

    printstyled("Test 2 \n", bold=true,color=:cyan)
    @test test2()

end

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


text:
"""
    PatternFold{T, P}
An abstract stype used as an interface for folded patterns such as `VectorFold`.
To implement the interface and inherit from it, a new structure must define three fields:
- `pattern::P`. Note that both `length(::P)` and `rand(:P)` methods must be available
- `gap::T`
- `folds::int`
"""
abstract type PatternFold{T, P} end

"""
    pattern(<:PatternFold)
Return the *pattern* of any `PatternFold`. The pattern defines the occurences of the first fold.
"""
pattern(pf) = pf.pattern

"""
    gap(<:PatternFold)
Return the *gap* between the starts of consecutive folds.
"""
gap(pf) = pf.gap

"""
    folds(<:PatternFold)
Return the number of *folds*. An infinite folded pattern returns `0`.
"""
folds(pf) = pf.folds

# TODO - look if another name is more appropriate
"""
    pattern_length(pf<:PatternFold)
Return the length of the basic pattern of `pf`.
"""
pattern_length(pf) = length(pattern(pf))

"""
    length(pf<:PatternFold)
Return the length of `pf` if unfolded.
"""

Base.length(pf::PatternFold) = pattern_length(pf) * folds(pf)

"""
    rand(pf<:PatternFold)
Returns a random value of `pf` as if it was unfolded.
"""
function rand(pf::PF) where {PF <: PatternFold}
    return rand(pattern(pf)) + rand(0:(folds(pf) - 1)) * gap(pf)
end


text:
export AbstractPhyHermite, AbstractProHermite, CstMaxDegree

struct AbstractPhyHermite <: Hermite
    P::ImmutablePolynomial
    scaled::Bool
end

(P::AbstractPhyHermite)(x::T) where {T<:Real} = P.P(x)*exp(-x^2/2)

struct AbstractProHermite <: Hermite
    P::ImmutablePolynomial
    scaled::Bool
end

(P::AbstractProHermite)(x) = P.P(x)*exp(-x^2/4)

const CstMaxDegree = 40


text:
using RegularizedLeastSquares, Images, PyPlot, Random

# image
N = 256
I = shepp_logan(N)

# sampling operator
idx = sort( shuffle( collect(1:N^2) )[1:div(N^2,2)] )
A = SamplingOp(idx,(N,N))

# generate undersampled data
y = A*vec(I)

# regularizer
reg = Regularization("TV", 0.01; shape=(N,N))

# solver
solver = createLinearSolver("admm",A; reg=reg, ρ=0.1, iterations=20)

Ireco = solve(solver,y)
Ireco = reshape(Ireco,N,N)


text:
# regularize/z-list.jl

include("Aodwt.jl")
include("diffl.jl")
include("diffs.jl")


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

using Random

# main tests

# issue #35800
# tested very early since it can be state-dependent
@test @inferred(mapreduce(x->count(!iszero,x), +, [rand(1)]; init = 0.)) == 1.0

function safe_mapslices(op, A, region)
    newregion = intersect(region, 1:ndims(A))
    return isempty(newregion) ? A : mapslices(op, A, dims = newregion)
end
safe_sum(A::Array{T}, region) where {T} = safe_mapslices(sum, A, region)
safe_prod(A::Array{T}, region) where {T} = safe_mapslices(prod, A, region)
safe_maximum(A::Array{T}, region) where {T} = safe_mapslices(maximum, A, region)
safe_minimum(A::Array{T}, region) where {T} = safe_mapslices(minimum, A, region)
safe_count(A::AbstractArray{T}, region) where {T} = safe_mapslices(count, A, region)
safe_sumabs(A::Array{T}, region) where {T} = safe_mapslices(sum, abs.(A), region)
safe_sumabs2(A::Array{T}, region) where {T} = safe_mapslices(sum, abs2.(A), region)
safe_maxabs(A::Array{T}, region) where {T} = safe_mapslices(maximum, abs.(A), region)
safe_minabs(A::Array{T}, region) where {T} = safe_mapslices(minimum, abs.(A), region)

@testset "test reductions over region: $region" for region in Any[
    1, 2, 3, 4, 5, (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4),
    (1, 2, 3), (1, 3, 4), (2, 3, 4), (1, 2, 3, 4)]
    Areduc = rand(3, 4, 5, 6)
    Breduc = rand(Bool, 3, 4, 5, 6)
    @assert axes(Areduc) == axes(Breduc)

    r = fill(NaN, map(length, Base.reduced_indices(axes(Areduc), region)))
    @test sum!(r, Areduc) ≈ safe_sum(Areduc, region)
    @test prod!(r, Areduc) ≈ safe_prod(Areduc, region)
    @test maximum!(r, Areduc) ≈ safe_maximum(Areduc, region)
    @test minimum!(r, Areduc) ≈ safe_minimum(Areduc, region)
    @test count!(r, Breduc) ≈ safe_count(Breduc, region)

    @test sum!(abs, r, Areduc) ≈ safe_sumabs(Areduc, region)
    @test sum!(abs2, r, Areduc) ≈ safe_sumabs2(Areduc, region)
    @test maximum!(abs, r, Areduc) ≈ safe_maxabs(Areduc, region)
    @test minimum!(abs, r, Areduc) ≈ safe_minabs(Areduc, region)
    @test count!(!, r, Breduc) ≈ safe_count(.!Breduc, region)

    # With init=false
    r2 = similar(r)
    fill!(r, 1)
    @test sum!(r, Areduc, init=false) ≈ safe_sum(Areduc, region) .+ 1
    fill!(r, 2.2)
    @test prod!(r, Areduc, init=false) ≈ safe_prod(Areduc, region)*2.2
    fill!(r, 1.8)
    @test maximum!(r, Areduc, init=false) ≈ fill!(r2, 1.8)
    fill!(r, -0.2)
    @test minimum!(r, Areduc, init=false) ≈ fill!(r2, -0.2)
    fill!(r, 1)
    @test count!(r, Breduc, init=false) ≈ safe_count(Breduc, region) .+ 1

    fill!(r, 8.1)
    @test sum!(abs, r, Areduc, init=false) ≈ safe_sumabs(Areduc, region) .+ 8.1
    fill!(r, 8.1)
    @test sum!(abs2, r, Areduc, init=false) ≈ safe_sumabs2(Areduc, region) .+ 8.1
    fill!(r, 1.5)
    @test maximum!(abs, r, Areduc, init=false) ≈ fill!(r2, 1.5)
    fill!(r, -1.5)
    @test minimum!(abs, r, Areduc, init=false) ≈ fill!(r2, -1.5)
    fill!(r, 1)
    @test count!(!, r, Breduc, init=false) ≈ safe_count(.!Breduc, region) .+ 1

    @test @inferred(sum(Areduc, dims=region)) ≈ safe_sum(Areduc, region)
    @test @inferred(prod(Areduc, dims=region)) ≈ safe_prod(Areduc, region)
    @test @inferred(maximum(Areduc, dims=region)) ≈ safe_maximum(Areduc, region)
    @test @inferred(minimum(Areduc, dims=region)) ≈ safe_minimum(Areduc, region)
    @test @inferred(count(Breduc, dims=region)) ≈ safe_count(Breduc, region)

    @test @inferred(sum(abs, Areduc, dims=region)) ≈ safe_sumabs(Areduc, region)
    @test @inferred(sum(abs2, Areduc, dims=region)) ≈ safe_sumabs2(Areduc, region)
    @test @inferred(maximum(abs, Areduc, dims=region)) ≈ safe_maxabs(Areduc, region)
    @test @inferred(minimum(abs, Areduc, dims=region)) ≈ safe_minabs(Areduc, region)
    @test @inferred(count(!, Breduc, dims=region)) ≈ safe_count(.!Breduc, region)

    @test isequal(
        @inferred(count(Breduc, dims=region, init=0x02)),
        safe_count(Breduc, region) .% UInt8 .+ 0x02,
    )
    @test isequal(
        @inferred(count(!, Breduc, dims=region, init=Int16(0))),
        safe_count(.!Breduc, region) .% Int16,
    )
end

# Combining dims and init
A = Array{Int}(undef, 0, 3)
@test_throws "reducing over an empty collection is not allowed" maximum(A; dims=1)
@test maximum(A; dims=1, init=-1) == reshape([-1,-1,-1], 1, 3)

# Test reduction along first dimension; this is special-cased for
# size(A, 1) >= 16
Breduc = rand(64, 3)
r = fill(NaN, map(length, Base.reduced_indices(axes(Breduc), 1)))
@test sum!(r, Breduc) ≈ safe_sum(Breduc, 1)
@test sum!(abs, r, Breduc) ≈ safe_sumabs(Breduc, 1)
@test sum!(abs2, r, Breduc) ≈ safe_sumabs2(Breduc, 1)
@test sum(Breduc, dims=1) ≈ safe_sum(Breduc, 1)
@test sum(abs, Breduc, dims=1) ≈ safe_sumabs(Breduc, 1)
@test sum(abs2, Breduc, dims=1) ≈ safe_sumabs2(Breduc, 1)

fill!(r, 4.2)
@test sum!(r, Breduc, init=false) ≈ safe_sum(Breduc, 1) .+ 4.2
fill!(r, -6.3)
@test sum!(abs, r, Breduc, init=false) ≈ safe_sumabs(Breduc, 1) .- 6.3
fill!(r, -1.1)
@test sum!(abs2, r, Breduc, init=false) ≈ safe_sumabs2(Breduc, 1) .- 1.1

# Small arrays with init=false
let A = reshape(1:15, 3, 5)
    R = fill(1, 3)
    @test sum!(R, A, init=false) == [36,41,46]
    R = fill(1, 1, 5)
    @test sum!(R, A, init=false) == [7 16 25 34 43]
end
let R = [2]
    A = reshape(1:6, 3, 2)
    @test prod!(R, A, init=false) == [1440]

    # min/max
    @test reduce(max, A, dims=1) == [3 6]
    @test reduce(min, A, dims=2) == reshape([1,2,3], 3, 1)
end

# Small integers
@test @inferred(sum(Int8[1], dims=1)) == [1]
@test @inferred(sum(UInt8[1], dims=1)) == [1]

# Complex types
@test typeof(@inferred(sum([1.0+1.0im], dims=1))) == Vector{ComplexF64}
@test typeof(@inferred(Base.sum(abs, [1.0+1.0im], dims=1))) == Vector{Float64}
@test typeof(@inferred(Base.sum(abs2, [1.0+1.0im], dims=1))) == Vector{Float64}
@test typeof(@inferred(prod([1.0+1.0im], dims=1))) == Vector{ComplexF64}
@test typeof(@inferred(Base.prod(abs, [1.0+1.0im], dims=1))) == Vector{Float64}
@test typeof(@inferred(Base.prod(abs2, [1.0+1.0im], dims=1))) == Vector{Float64}

@testset "heterogeneously typed arrays" begin
    for x in (sum(Union{Float32, Float64}[1.0], dims=1),
              prod(Union{Float32, Float64}[1.0], dims=1))
        @test x == [1.0]
        @test x isa Vector{Float64}
    end

    x = sum(Real[1.0], dims=1)
    @test x == [1.0]
    @test x isa Vector{Real}

    x = mapreduce(cos, +, Union{Int,Missing}[1, 2], dims=1)
    @test x == mapreduce(cos, +, [1, 2], dims=1)
    @test x isa Vector{Float64}
end

@test reduce((a,b) -> a|b, [true false; false false], dims=1, init=false) == [true false]
let R = reduce((a,b) -> a+b, [1 2; 3 4], dims=2, init=0.0)
    @test eltype(R) == Float64
    @test R ≈ [3,7]
end
@test reduce((a,b) -> a+b, [1 2; 3 4], dims=1, init=0) == [4 6]

# inferred return types
@test typeof(@inferred(reduce(+, ones(3,3,3), dims=1, init=0.0))) == Array{Float64, 3}

@testset "empty cases" begin
    A = Matrix{Int}(undef, 0,1)
    @test sum(A) === 0
    @test prod(A) === 1
    @test_throws ["reducing over an empty",
                  "consider supplying `init`"] minimum(A)
    @test_throws "consider supplying `init`" maximum(A)

    @test isequal(sum(A, dims=1), zeros(Int, 1, 1))
    @test isequal(sum(A, dims=2), zeros(Int, 0, 1))
    @test isequal(sum(A, dims=(1, 2)), zeros(Int, 1, 1))
    @test isequal(sum(A, dims=3), zeros(Int, 0, 1))
    @test isequal(prod(A, dims=1), fill(1, 1, 1))
    @test isequal(prod(A, dims=2), fill(1, 0, 1))
    @test isequal(prod(A, dims=(1, 2)), fill(1, 1, 1))
    @test isequal(prod(A, dims=3), fill(1, 0, 1))

    for f in (minimum, maximum)
        @test_throws "reducing over an empty collection is not allowed" f(A, dims=1)
        @test isequal(f(A, dims=2), zeros(Int, 0, 1))
        @test_throws "reducing over an empty collection is not allowed" f(A, dims=(1, 2))
        @test isequal(f(A, dims=3), zeros(Int, 0, 1))
    end
    for f in (findmin, findmax)
        @test_throws ArgumentError f(A, dims=1)
        @test isequal(f(A, dims=2), (zeros(Int, 0, 1), zeros(Int, 0, 1)))
        @test_throws ArgumentError f(A, dims=(1, 2))
        @test isequal(f(A, dims=3), (zeros(Int, 0, 1), zeros(Int, 0, 1)))
    end

end

## findmin/findmax/minimum/maximum

A = [1.0 5.0 6.0;
     5.0 2.0 4.0]
for (tup, rval, rind) in [((1,), [1.0 2.0 4.0], [CartesianIndex(1,1) CartesianIndex(2,2) CartesianIndex(2,3)]),
                          ((2,), reshape([1.0,2.0], 2, 1), reshape([CartesianIndex(1,1),CartesianIndex(2,2)], 2, 1)),
                          ((1,2), fill(1.0,1,1),fill(CartesianIndex(1,1),1,1))]
    @test findmin(A, dims=tup) == (rval, rind)
    @test findmin!(similar(rval), similar(rind), A) == (rval, rind)
    @test isequal(minimum(A, dims=tup), rval)
    @test isequal(minimum!(similar(rval), A), rval)
    @test isequal(minimum!(copy(rval), A, init=false), rval)
end

for (tup, rval, rind) in [((1,), [5.0 5.0 6.0], [CartesianIndex(2,1) CartesianIndex(1,2) CartesianIndex(1,3)]),
                          ((2,), reshape([6.0,5.0], 2, 1), reshape([CartesianIndex(1,3),CartesianIndex(2,1)], 2, 1)),
                          ((1,2), fill(6.0,1,1),fill(CartesianIndex(1,3),1,1))]
    @test findmax(A, dims=tup) == (rval, rind)
    @test findmax!(similar(rval), similar(rind), A) == (rval, rind)
    @test isequal(maximum(A, dims=tup), rval)
    @test isequal(maximum!(similar(rval), A), rval)
    @test isequal(maximum!(copy(rval), A, init=false), rval)
end

@testset "missing in findmin/findmax" begin
    B = [1.0 missing NaN;
         5.0 NaN missing]
    for (tup, rval, rind) in [(1, [5.0 missing missing], [CartesianIndex(2, 1) CartesianIndex(1, 2) CartesianIndex(2, 3)]),
                              (2, [missing; missing],    [CartesianIndex(1, 2) CartesianIndex(2, 3)] |> permutedims)]
        (rval′, rind′) = findmax(B, dims=tup)
        @test all(rval′ .=== rval)
        @test all(rind′ .== rind)
        @test all(maximum(B, dims=tup) .=== rval)
    end

    for (tup, rval, rind) in [(1, [1.0 missing missing], [CartesianIndex(1, 1) CartesianIndex(1, 2) CartesianIndex(2, 3)]),
                              (2, [missing; missing],    [CartesianIndex(1, 2) CartesianIndex(2, 3)] |> permutedims)]
        (rval′, rind′) = findmin(B, dims=tup)
        @test all(rval′ .=== rval)
        @test all(rind′ .== rind)
        @test all(minimum(B, dims=tup) .=== rval)
    end
end

@testset "reducedim_init min/max unorderable handling" begin
    x = Any[1.0, NaN]
    y = [1, missing]
    for (v, rval1, rval2) in [(x, [NaN], x),
                              (y, [missing], y),
                              (Any[1. NaN; 1. 1.], Any[1. NaN], Any[NaN, 1.])]
        for f in (minimum, maximum)
            @test all(f(v, dims=1) .=== rval1)
            @test all(f(v, dims=2) .=== rval2)
        end
    end
end

#issue #23209

A = [1.0 3.0 6.0;
     NaN 2.0 4.0]
for (tup, rval, rind) in [((1,), [NaN 2.0 4.0], [CartesianIndex(2,1) CartesianIndex(2,2) CartesianIndex(2,3)]),
                          ((2,), reshape([1.0, NaN], 2, 1), reshape([CartesianIndex(1,1),CartesianIndex(2,1)], 2, 1)),
                          ((1,2), fill(NaN,1,1),fill(CartesianIndex(2,1),1,1))]
    @test isequal(findmin(A, dims=tup), (rval, rind))
    @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(minimum(A, dims=tup), rval)
    @test isequal(minimum!(similar(rval), A), rval)
    @test isequal(minimum!(copy(rval), A, init=false), rval)
    @test isequal(Base.reducedim!(min, copy(rval), A), rval)
end

for (tup, rval, rind) in [((1,), [NaN 3.0 6.0], [CartesianIndex(2,1) CartesianIndex(1,2) CartesianIndex(1,3)]),
                          ((2,), reshape([6.0, NaN], 2, 1), reshape([CartesianIndex(1,3),CartesianIndex(2,1)], 2, 1)),
                          ((1,2), fill(NaN,1,1),fill(CartesianIndex(2,1),1,1))]
    @test isequal(findmax(A, dims=tup), (rval, rind))
    @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(maximum(A, dims=tup), rval)
    @test isequal(maximum!(similar(rval), A), rval)
    @test isequal(maximum!(copy(rval), A, init=false), rval)
    @test isequal(Base.reducedim!(max, copy(rval), A), rval)
end

# issue #28320
@testset "reducedim issue with abstract complex arrays" begin
let A = Complex[1.5 0.5]
    @test mapreduce(abs2, +, A, dims=2) == reshape([2.5], 1, 1)
    @test sum(abs2, A, dims=2) == reshape([2.5], 1, 1)
    @test prod(abs2, A, dims=2) == reshape([0.5625], 1, 1)
    @test maximum(abs2, A, dims=2) == reshape([2.25], 1, 1)
    @test minimum(abs2, A, dims=2) == reshape([0.25], 1, 1)
end
end

A = [1.0 NaN 6.0;
     NaN 2.0 4.0]
for (tup, rval, rind) in [((1,), [NaN NaN 4.0], [CartesianIndex(2,1) CartesianIndex(1,2) CartesianIndex(2,3)]),
                          ((2,), reshape([NaN, NaN], 2, 1), reshape([CartesianIndex(1,2),CartesianIndex(2,1)], 2, 1)),
                          ((1,2), fill(NaN,1,1),fill(CartesianIndex(2,1),1,1))]
    @test isequal(findmin(A, dims=tup), (rval, rind))
    @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(minimum(A, dims=tup), rval)
    @test isequal(minimum!(similar(rval), A), rval)
    @test isequal(minimum!(copy(rval), A, init=false), rval)
end

for (tup, rval, rind) in [((1,), [NaN NaN 6.0], [CartesianIndex(2,1) CartesianIndex(1,2) CartesianIndex(1,3)]),
                          ((2,), reshape([NaN, NaN], 2, 1), reshape([CartesianIndex(1,2),CartesianIndex(2,1)], 2, 1)),
                          ((1,2), fill(NaN,1,1),fill(CartesianIndex(2,1),1,1))]
    @test isequal(findmax(A, dims=tup), (rval, rind))
    @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(maximum(A, dims=tup), rval)
    @test isequal(maximum!(similar(rval), A), rval)
    @test isequal(maximum!(copy(rval), A, init=false), rval)
end

A = [Inf -Inf Inf  -Inf;
     Inf  Inf -Inf -Inf]
for (tup, rval, rind) in [((1,), [Inf -Inf -Inf -Inf], [CartesianIndex(1,1) CartesianIndex(1,2) CartesianIndex(2,3) CartesianIndex(1,4)]),
                          ((2,), reshape([-Inf -Inf], 2, 1), reshape([CartesianIndex(1,2),CartesianIndex(2,3)], 2, 1)),
                          ((1,2), fill(-Inf,1,1),fill(CartesianIndex(1,2),1,1))]
    @test isequal(findmin(A, dims=tup), (rval, rind))
    @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(minimum(A, dims=tup), rval)
    @test isequal(minimum!(similar(rval), A), rval)
    @test isequal(minimum!(copy(rval), A, init=false), rval)
end

for (tup, rval, rind) in [((1,), [Inf Inf Inf -Inf], [CartesianIndex(1,1) CartesianIndex(2,2) CartesianIndex(1,3) CartesianIndex(1,4)]),
                          ((2,), reshape([Inf Inf], 2, 1), reshape([CartesianIndex(1,1),CartesianIndex(2,1)], 2, 1)),
                          ((1,2), fill(Inf,1,1),fill(CartesianIndex(1,1),1,1))]
    @test isequal(findmax(A, dims=tup), (rval, rind))
    @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(maximum(A, dims=tup), rval)
    @test isequal(maximum!(similar(rval), A), rval)
    @test isequal(maximum!(copy(rval), A, init=false), rval)
end

A = [BigInt(10)]
for (tup, rval, rind) in [((2,), [BigInt(10)], [1])]
    @test isequal(findmin(A, dims=tup), (rval, rind))
    @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(minimum(A, dims=tup), rval)
    @test isequal(minimum!(similar(rval), A), rval)
    @test isequal(minimum!(copy(rval), A, init=false), rval)
end

for (tup, rval, rind) in [((2,), [BigInt(10)], [1])]
    @test isequal(findmax(A, dims=tup), (rval, rind))
    @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(maximum(A, dims=tup), rval)
    @test isequal(maximum!(similar(rval), A), rval)
    @test isequal(maximum!(copy(rval), A, init=false), rval)
end

A = [BigInt(-10)]
for (tup, rval, rind) in [((2,), [BigInt(-10)], [1])]
    @test isequal(findmin(A, dims=tup), (rval, rind))
    @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(minimum(A, dims=tup), rval)
    @test isequal(minimum!(similar(rval), A), rval)
    @test isequal(minimum!(copy(rval), A, init=false), rval)
end

for (tup, rval, rind) in [((2,), [BigInt(-10)], [1])]
    @test isequal(findmax(A, dims=tup), (rval, rind))
    @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(maximum(A, dims=tup), rval)
    @test isequal(maximum!(similar(rval), A), rval)
    @test isequal(maximum!(copy(rval), A, init=false), rval)
end

A = [BigInt(10) BigInt(-10)]
for (tup, rval, rind) in [((2,), reshape([BigInt(-10)], 1, 1), reshape([CartesianIndex(1,2)], 1, 1))]
    @test isequal(findmin(A, dims=tup), (rval, rind))
    @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(minimum(A, dims=tup), rval)
    @test isequal(minimum!(similar(rval), A), rval)
    @test isequal(minimum!(copy(rval), A, init=false), rval)
end

for (tup, rval, rind) in [((2,), reshape([BigInt(10)], 1, 1), reshape([CartesianIndex(1,1)], 1, 1))]
    @test isequal(findmax(A, dims=tup), (rval, rind))
    @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(maximum(A, dims=tup), rval)
    @test isequal(maximum!(similar(rval), A), rval)
    @test isequal(maximum!(copy(rval), A, init=false), rval)
end

A = ["a", "b"]
for (tup, rval, rind) in [((1,), ["a"], [1])]
    @test isequal(findmin(A, dims=tup), (rval, rind))
    @test isequal(findmin!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(minimum(A, dims=tup), rval)
    @test isequal(minimum!(similar(rval), A), rval)
    @test isequal(minimum!(copy(rval), A, init=false), rval)
end

for (tup, rval, rind) in [((1,), ["b"], [2])]
    @test isequal(findmax(A, dims=tup), (rval, rind))
    @test isequal(findmax!(similar(rval), similar(rind), A), (rval, rind))
    @test isequal(maximum(A, dims=tup), rval)
    @test isequal(maximum!(similar(rval), A), rval)
    @test isequal(maximum!(copy(rval), A, init=false), rval)
end

# issue #6672
@test sum(Real[1 2 3; 4 5.3 7.1], dims=2) == reshape([6, 16.4], 2, 1)
@test sum(Any[1 2;3 4], dims=1) == [4 6]
@test sum(Vector{Int}[[1,2],[4,3]], dims=1)[1] == [5,5]

@testset "Issue #10461. region=$region" for region in Any[-1, 0, (-1, 2), [0, 1], (1,-2,3), [0 1;
                                                     2 3], "hello"]
    Areduc = rand(3, 4, 5, 6)

    @test_throws ArgumentError sum(Areduc, dims=region)
    @test_throws ArgumentError prod(Areduc, dims=region)
    @test_throws ArgumentError maximum(Areduc, dims=region)
    @test_throws ArgumentError minimum(Areduc, dims=region)
    @test_throws ArgumentError sum(abs, Areduc, dims=region)
    @test_throws ArgumentError sum(abs2, Areduc, dims=region)
    @test_throws ArgumentError maximum(abs, Areduc, dims=region)
    @test_throws ArgumentError minimum(abs, Areduc, dims=region)
end

# issue #26488
@testset "don't map over initial values not provided" begin
    @test sum(x->x+1, [1], dims=1)[1] === sum(x->x+1, [1]) === 2
    @test prod(x->x+1, [1], dims=1)[1] === prod(x->x+1, [1]) === 2
    @test mapreduce(x->x+1, +, [1], dims=1)[1] === mapreduce(x->x+1, +, [1]) === 2
    @test mapreduce(x->x+1, *, [1], dims=1)[1] === mapreduce(x->x+1, *, [1]) === 2
    @test mapreduce(!, &, [false], dims=1)[1] === mapreduce(!, &, [false]) === true
    @test mapreduce(!, |, [true], dims=1)[1] === mapreduce(!, |, [true]) === false
    @test mapreduce(x->1/x, max, [1], dims=1)[1] === mapreduce(x->1/x, max, [1]) === 1.0
    @test mapreduce(x->-1/x, min, [1], dims=1)[1] === mapreduce(x->-1/x, min, [1]) === -1.0
end

# check type of result
@testset "type of sum(::Array{$T}" for T in [UInt8, Int8, Int32, Int64, BigInt]
    result = sum(T[1 2 3; 4 5 6; 7 8 9], dims=2)
    @test result == hcat([6, 15, 24])
    @test eltype(result) === (T <: Base.SmallSigned ? Int :
                              T <: Base.SmallUnsigned ? UInt :
                              T)
end

@testset "argmin/argmax" begin
    B = reshape(3^3:-1:1, (3, 3, 3))
    @test B[argmax(B, dims=[2, 3])] == maximum(B, dims=[2, 3])
    @test B[argmin(B, dims=[2, 3])] == minimum(B, dims=[2, 3])
end

@testset "in-place reductions with mismatched dimensionalities" begin
    B = reshape(1:24, 4, 3, 2)
    for R in (fill(0, 4), fill(0, 4, 1), fill(0, 4, 1, 1))
        @test @inferred(maximum!(R, B)) == reshape(21:24, size(R))
        @test @inferred(minimum!(R, B)) == reshape(1:4, size(R))
    end
    for R in (fill(0, 1, 3), fill(0, 1, 3, 1))
        @test @inferred(maximum!(R, B)) == reshape(16:4:24, size(R))
        @test @inferred(minimum!(R, B)) == reshape(1:4:9, size(R))
    end
    @test_throws DimensionMismatch maximum!(fill(0, 4, 1, 1, 1), B)
    @test_throws DimensionMismatch minimum!(fill(0, 4, 1, 1, 1), B)
    @test_throws DimensionMismatch maximum!(fill(0, 1, 3, 1, 1), B)
    @test_throws DimensionMismatch minimum!(fill(0, 1, 3, 1, 1), B)
    @test_throws DimensionMismatch maximum!(fill(0, 1, 1, 2, 1), B)
    @test_throws DimensionMismatch minimum!(fill(0, 1, 1, 2, 1), B)
end

# issue #26709
@testset "dimensional reduce with custom non-bitstype types" begin
    struct Variable
        name::Symbol
    end
    struct AffExpr
        vars::Vector{Variable}
    end
    Base.zero(::Union{Variable, Type{Variable}, AffExpr}) = AffExpr(Variable[])
    Base.:+(v::Variable, w::Variable) = AffExpr([v, w])
    Base.:+(aff::AffExpr, v::Variable) = AffExpr([aff.vars; v])
    Base.:+(aff1::AffExpr, aff2::AffExpr) = AffExpr([aff1.vars; aff2.vars])
    Base.:(==)(a::Variable, b::Variable) = a.name == b.name
    Base.:(==)(a::AffExpr, b::AffExpr) = a.vars == b.vars

    @test sum([Variable(:x), Variable(:y)], dims=1) == [AffExpr([Variable(:x), Variable(:y)])]
end

# count
@testset "count: throw on non-bool types" begin
    @test_throws TypeError count([1], dims=1)
    @test_throws TypeError count!([1], [1])
end

@test @inferred(count(false:true, dims=:, init=0x0004)) === 0x0005
@test @inferred(count(isodd, reshape(1:9, 3, 3), dims=:, init=Int128(0))) === Int128(5)

@testset "reduced_index for BigInt (issue #39995)" begin
    for T in [Int8, Int16, Int32, Int64, Int128, BigInt]
        r = T(1):T(2)
        ax = axes(r, 1)
        axred = Base.reduced_index(ax)
        @test axred == Base.OneTo(1)
        @test typeof(axred) === typeof(ax)
        r_red = reduce(+, r, dims = 1)
        @test eltype(r_red) == T
        @test r_red == [3]
    end
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
function rankmodified{T<:Real}(a::Array{T,1})
    indexin(a, a)
end


text:
# Get an unsigned position
function uposition(stream)
    UInt(position(stream))
end

# Function that generates a name based on the name of the cached function
function generate_cache_filename()
    _filename = "_" * string(hash(rand()), base=16) * "_.bin"
    return abspath(_filename)
end

# Function that retrieves one entry from a stream
function load_disk_cache_entry(io::IO,
                               pos_start::UInt,
                               pos_end::UInt,
                               decompressor=Noop)
    seek(io, pos_start)
    cbuf = read(io, pos_end-pos_start)
    buf = transcode(decompressor, cbuf)
    datum = deserialize(IOBuffer(buf))
    return datum
end

# Function that stores one entry to a stream
function store_disk_cache_entry(io::IO,
                                pos_start::UInt,
                                datum,
                                compressor=Noop)
    seek(io, pos_start)
    buf = IOBuffer()
    serialize(buf, datum)
    posbuf = uposition(buf)
    cbuf = transcode(compressor, buf.data[1:posbuf])
    write(io, cbuf)
    return uposition(io)
end

# Function that checks the consistency of the cache and `offset`
# field value of the Cache object
function check_disk_cache(filename::String,
                          offsets::Dict,
                          output_type::Type
                         )::Bool
    open(filename, "r") do fid
        _, decompressor = get_transcoders(filename)
        try
            O = deserialize(fid);
            @assert O == output_type
            for (_, (startpos, endpos)) in offsets
                load_disk_cache_entry(fid, startpos, endpos, decompressor)
            end
            return true
        catch excep
            return false
        end
    end
end

function get_transcoders(filename::T="") where {T<:AbstractString}
    ext = split(filename, ".")[end]
    if ext == "lz4"
        compressor = LZ4FrameCompressor
        decompressor = LZ4FrameDecompressor
    elseif ext == "bz2" || ext == "bzip2"
        compressor = Bzip2Compressor
        decompressor = Bzip2Decompressor
    elseif ext == "gz" || ext == "gzip"
        compressor = GzipCompressor
        decompressor = GzipDecompressor
    else
        compressor = Noop  # no compression
        decompressor = Noop  # no compression
    end
    return compressor, decompressor
end

cache2nt(cache::Cache) = (name=cache.name,
                          filename=cache.filename,
                          func_def=cache.func_def,
                          cache=cache.cache,
                          offsets=cache.offsets,
                          history=cache.history,
                          max_size=cache.max_size)

serialize(stream::IO, cache::Cache) = serialize(stream, cache2nt(cache))

serialize(filename::AbstractString, cache::Cache) = serialize(filename, cache2nt(cache))

function deserialize(stream_or_filename, ::Type{Cache}; func=nothing)
    nt = deserialize(stream_or_filename)
    if nt.func_def === nothing && func === nothing
        throw(ErrorException("Cannot reconstruct cache; use the `func` keyword argument."))
    elseif func === nothing
        func = try
                eval(Meta.parse(nt.func_def))
            catch e
                throw(ErrorException("Cannot reconstruct cache:\n$e"))
            end
    else
        # use keyword func over serialized function definition code
    end
    return Cache(func, nt...)
end


text:
function updateneuron!(n::Neuron)
end
export updateneuron!

function collectmessages!(n::Neuron, c)
end
export collectmessages!

max_stars_repo_name: jbrea/NetSim

max_issues_repo_name: jbrea/NetSim

max_forks_repo_name: jbrea/NetSim


text:
# HellwarthEffectiveFrequency.jl
#   - use Hellwarth et al. 1999 PRB method to reduce multiple phonon modes to a single effective frequency

push!(LOAD_PATH,"../src/") # load module from local directory
using PolaronMobility

# ((freq THz)) ((IR Activity / e^2 amu^-1))
# These data from MAPbI3-Cubic_PeakTable.csv
# https://github.com/WMD-group/Phonons/tree/master/2015_MAPbI3/SimulatedSpectra
# Data published in Brivio2015 (PRB)
# https://doi.org/10.1103/PhysRevB.92.144308
MAPI= [
96.20813558773261 0.4996300522819191
93.13630357703363 1.7139631746083817
92.87834578121567 0.60108592692181
92.4847918585963 0.0058228799414729
92.26701437594754 0.100590086574602
89.43972834606603 0.006278895133832249
46.89209141511332 0.2460894564364346
46.420949316788 0.14174282581124137
44.0380222871706 0.1987196948553428
42.89702947649343 0.011159939465770681
42.67180170168193 0.02557751102757614
41.46971205834201 0.012555230726601503
37.08982543385215 0.00107488277468418
36.53555265689563 0.02126940080871224
30.20608114002676 0.009019481779712388
27.374810898415028 0.03994453721421388
26.363055017011728 0.05011922682554448
9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
#0.022939578929507105 8.355742795827834e-05   # Acoustic modes!
#0.04882611767873102 8.309858592685e-06
#0.07575149723846182 2.778248540373041e-05
]

MAPbBr3 = [
#v / THz IR Activity / e^2 amu^-1
96.25494581101505 0.4856105419754306
93.36827101597564 1.5237045178260684
93.19828072025739 0.54951133126244
92.59570097621418 0.0017497942636829
92.36075758752216 0.054483693922451705
89.49998024367731 0.003950754031148851
47.175408954380885 0.21156964115237242
46.57576693654983 0.11627478349663156
44.3361503231313 0.19163912283684897
42.967352451897554 0.011987727827091217
42.71066583604171 0.027797225985699996
41.490751582679316 0.01234463397562453
37.22033177431264 0.001649578352399309
36.57010484641483 0.019161684306162618
30.438999700507342 0.007365167772350685
27.625465570185764 0.035161502048866464
26.458785163031635 0.038174133532805965
9.947074760143552 0.0012355225072173002
4.493935108724231 0.07714630445769342
4.215763207131703 0.01724042408702847
4.021795301113039 0.06167913225780508
3.2393471899404047 0.040335730459227545
2.914057478007975 0.27361583921530497
2.6223962950449407 0.23493646364819373
2.400751729577614 0.2753271857365242
2.370809478597967 0.04428317547154772
2.0388405406936783 0.0517120158657721
1.4017373572275915 0.0046510016948825
1.3634089519825427 0.0018248126743814707
1.336413930635821 0.009190430969212338
1.2144969737450764 0.013885351107636132
1.0763905048534899 0.0014676450708402083
0.7427782546143306 0.008829646844819504
#0.012830037764855392 6.62960832415985e-06
#0.02156203260482768 4.708045983230229e-05
#0.05156888885824344 2.602755859934495e-05
]

MAPbCl3= [
#v / THz IR Activity / e^2 amu^-1
96.65380905798915 0.4542931751264273
93.30434573145013 1.4046841552487999
93.27126600254645 0.5286008698609284
92.719992888923 0.0015896642905404556
92.45885414527524 0.0866990370170258
89.61967951587117 0.0031392395308660263
47.44029120639781 0.18901278025443116
46.764089990499585 0.09717407330258415
44.60198174803584 0.19436643102794884
43.04971569182655 0.012986533289171532
42.7637228646608 0.02998334671150052
41.49924134521671 0.010714577213688365
37.34275112339727 0.0028788161278648656
36.62941263763097 0.018016387094969555
30.651112964494686 0.006550727272302003
27.89469865012452 0.03515737503204058
26.603746656189923 0.03497660426639104
10.493140658721968 9.69032373922839e-05
5.242459815070028 0.12301715153349158
4.8348539006445845 0.11307729502447525
4.633936235943856 0.009371996258413128
4.1582548301810665 0.08368265471427451
3.9021555976797266 0.44406366947684967
3.3808364269416287 0.3379317982808455
3.099250159114845 0.38405867852912545
2.8466202448045217 0.025463767673713028
2.5900539956188333 0.07754012506771049
2.0931563078880493 0.008533555182635471
1.869999049657922 0.011347668353507898
1.8382471614901006 0.0005371757988
1.829843559303894 0.019270325513936128
1.6503257694697964 0.000950380841742176
1.0962868606064022 0.018510433944796113
#0.032159306112378404 2.774125107895896e-05 #Acoustic
#0.04846756590412458 7.796048392205337e-05
#0.05253089286287104 2.8730156133e-06
]

MAPI_low=MAPI[19:33,:] # Just inorganic components, everything below 10THz; modes 3-18

println("\n\nMAPI: BScheme (athermal)")
println("\t MAPI: (all values)")
HellwarthBScheme(MAPI)
println("\t MAPI: (low-frequency, non molecular IR, only)")
HellwarthBScheme(MAPI_low)

println("\nMAPI: AScheme (thermal)")
println("\t MAPI: (all values)")
HellwarthAScheme(MAPI)
println("\t MAPI: (low-frequency, non molecular IR)")
HellwarthAScheme(MAPI_low)

println("\n\n> MAPI/MAPbBr3/MAPbCl3 - with simple B-scheme...")
for LOs in [MAPI, MAPbBr3, MAPbCl3]
    println("All modes...")
    HellwarthBScheme(LOs)
    println("Low frequency (cage modes only)...")
    HellwarthBScheme(LOs[19:33,:])
    println()
end

println("\n Test summation of Lorentz oscillators to get to static dielectric constant from i.r. modes.")
# Integrate through Lorentz oscillators to get dielectric fn
# Should give 'extra' contribution from these modes, extrapolated to zero omega
function integrate_dielectric(LO,V0)
    summate=sum( (LO[:,2])./(LO[:,1].^2) )
#    summate*4*π/V0
    summate*ε0/V0
end

const Å=1E-10 # angstrom in metres
const r=6.29Å # Sensible cubic cell size
const V0=(r)^3
println("volume: $V0")
const amu=1.66054e-27
const ε0=8.854187817E-12
const eV = const q = const ElectronVolt = 1.602176487e-19;                         # kg m2 / s2

# Change to SI; not actually needed as units cancel everywhere (?)
MAPI_SI = [ MAPI[:,1].*10^12*2*π MAPI[:,2]./(q^2/amu) ]

println(" MAPI: ",integrate_dielectric(MAPI,1.0))
println(" MAPI_low: ",integrate_dielectric(MAPI_low,1.0))
println(" MAPI_SI: ",integrate_dielectric(MAPI_SI,V0))
println(" MAPI_SI: fudged epislon0 ",integrate_dielectric(MAPI_SI,V0)*ε0/(4*π))
println(" MAPI_SI_low: fudged epislon0 ",integrate_dielectric(MAPI_SI[19:33,:],V0)*ε0/(4*π))

println()
#println("From ε_S-ε_Inf, expect this to be: ",ε_S-ε_Inf)


text:
function neg_bin_dist(n, k, p, q)
    return binomial(n-1, k-1) * p^k * q^(n-k)
end

str = readline()
n = parse(Int64, readline())
strvec = split(str, ' ')
intvec = map(x -> parse(Int64, x), strvec)
a = intvec[1]
b = intvec[2]

p = a / b
q = 1 - p

println(round(neg_bin_dist(n, 1, p, q), 3))

max_stars_repo_path: Tutorials/10 Days of Statistics/Day 4/geometric_distribution1.jl

max_stars_repo_name: xuedong/hacker-rank

max_issues_repo_name: xuedong/hacker-rank

max_forks_repo_name: xuedong/hacker-rank


text:
using MLKernels
using JuMP
using MLLabelUtils
using StatsBase
using NearestNeighbors
import MLBase

const learning_pool_enc = LabelEnc.NativeLabels([:U, :Lin, :Lout])
const class_label_enc = LabelEnc.NativeLabels([:inlier, :outlier])
const OPT_PRECISION = 1e-7

@enum ModelState model_created=1 model_initialized=2 model_fitted=3

const Scope = Union{Val{:Subspace}, Val{:Global}}

abstract type ModelException <: Exception end

struct ModelStateException <: ModelException
    actual::ModelState
    expected::ModelState
end

Base.showerror(io::IO, e::ModelStateException) = print(io, "Invalid model state: $(e.actual). Model must at least be in state $(e.expected).")

struct ModelInvariantException <: ModelException
    msg
end

Base.showerror(io::IO, e::ModelInvariantException) = print(io, "Model invariant error. $(e.msg)")

abstract type OCClassifier end

abstract type SVDDClassifier <: OCClassifier end

abstract type SubOCClassifier <: OCClassifier end

abstract type InitializationStrategy end


text:
# Try to open at least all downloaded datasets
@testset "Reading rasters" begin
    for f in remotefiles
        @testset "Reading $f" begin
            ga = GeoArrays.read(f)
        end
    end
end

@testset "Writing rasters" begin

    @testset "Simplest version" begin
        ga = GeoArray(rand(100,200,3))
        GeoArrays.write!("test.tif", ga)
    end
    @testset "Nodata" begin
        ga = GeoArray(Array{Union{Missing, Int32}}(rand(1:10,100,200,3)))
        GeoArrays.write!("test_nodata.tif", ga, 1)
    end
end


text:
import Plots
import Dierckx

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

scale_time_to_days(t) = t / 86400
scale_time_to_hours(t) = t / 3600
scale_z_to_kilometers(t) = t / 10^3
scale_z_to_meters(t) = t
xlabel(::typeof(scale_time_to_hours)) = "Time [hr]"
xlabel(::typeof(scale_time_to_days)) = "Time [days]"
ylabel(::typeof(scale_z_to_kilometers)) = "Altitude (z) [km]"
ylabel(::typeof(scale_z_to_meters)) = "Altitude (z) [m]"

function save_tz_profile(
        profile;
        scale_time = scale_time_to_hours,
        scale_z = scale_z_to_kilometers,
        units = "",
        label = "",
        t_range,
        z_range,
    )
    prof = profile(Float64)

    # We could alternatively use this
    # as it's nicely documenting how
    # to use the profiles, but we want
    # to scale the altitude
    # Plots.contourf(t_range, z_range, (t, z) -> prof(t, z))
    data = prof.(t_range', z_range)
    Plots.contourf(scale_time.(t_range), scale_z.(z_range), data; c = :viridis)
    Plots.xlabel!(xlabel(scale_time))
    Plots.ylabel!(ylabel(scale_z))
    Plots.title!("$(nameof(profile)) $units, ($label)")
    Plots.savefig("tz_$(nameof(profile)).png")
end

function save_z_profile(
        profile;
        xlabel,
        scale_z = scale_z_to_kilometers,
        z_range,
    )
    prof = profile(Float64)
    # data = (z) -> prof(z)
    data = prof.(z_range)
    Plots.plot(data, scale_z.(z_range))
    Plots.xlabel!("$xlabel $(units(xlabel))")
    Plots.ylabel!(ylabel(scale_z))
    Plots.title!("$(nameof(profile))")
    Plots.savefig("z_$(nameof(profile)).png")
end

function save_t_profile(
        profile;
        label,
        scale_time = scale_time_to_hours,
        t_range,
    )
    prof = profile(Float64)
    # data = (t) -> prof(t)
    data = prof.(t_range)
    Plots.plot(scale_time.(t_range), data)
    Plots.xlabel!("Time $(xlabel(scale_time))")
    Plots.ylabel!("$label $(units(label))")
    Plots.title!("$(nameof(profile))")
    Plots.savefig("t_$(nameof(profile)).png")
end


text:
"""
    divergence(f)
"""
divergence(f) = tr(gradient(f))

function divergence(a::AbstractArray{<:Field})
  ag = gradient(a)
  operate_arrays_of_fields(tr,ag)
end

function symmetric_gradient end

"""
    symmetric_gradient(f)
"""
symmetric_gradient(f) = symmetric_part(gradient(f))

function symmetric_gradient(a::AbstractArray{<:Field})
  ag = gradient(a)
  operate_arrays_of_fields(symmetric_part,ag)
end

"""
    const ε = symmetric_gradient

Alias for the symmetric gradient
"""
const ε = symmetric_gradient

"""
    curl(f)
"""
curl(f) = grad2curl(gradient(f))

function curl(a::AbstractArray{<:Field})
  ag = gradient(a)
  operate_arrays_of_fields(grad2curl,ag)
end

"""
    grad2curl(∇f)
"""
function grad2curl(f)
  @abstractmethod
end

grad2curl(a::GridapType) = operate(grad2curl,a)

grad2curl(f::Field) = operate_fields(grad2curl,f)

@inline function grad2curl(∇u::TensorValue{2})
  ∇u[1,2] - ∇u[2,1]
end

@inline function grad2curl(∇u::TensorValue{3})
  c1 = ∇u[2,3] - ∇u[3,2]
  c2 = ∇u[3,1] - ∇u[1,3]
  c3 = ∇u[1,2] - ∇u[2,1]
  VectorValue(c1,c2,c3)
end

function laplacian end

"""
    const Δ = laplacian

Alias for the `laplacian` function
"""
const Δ = laplacian

"""
    laplacian(f)
"""
function laplacian(f)
  g = gradient(f)
  divergence(g)
end

"""
    ∇⋅f

Equivalent to

    divergence(f)
"""
dot(::typeof(∇),f) = divergence(f)
dot(::typeof(∇),f::GridapType) = divergence(f)

function (*)(::typeof(∇),f)
  msg = "Syntax ∇*f has been removed, use ∇⋅f (\\nabla \\cdot f) instead"
  error(msg)
end

function (*)(::typeof(∇),f::GridapType)
  msg = "Syntax ∇*f has been removed, use ∇⋅f (\\nabla \\cdot f) instead"
  error(msg)
end

"""
    outer(∇,f)

Equivalent to

    gradient(f)
"""
outer(::typeof(∇),f) = gradient(f)
outer(::typeof(∇),f::GridapType) = gradient(f)

"""
    outer(f,∇)

Equivalent to

    transpose(gradient(f))
"""
outer(f,::typeof(∇)) = transpose(gradient(f))
outer(f::GridapType,::typeof(∇)) = transpose(gradient(f))

"""
    cross(∇,f)

Equivalent to

    curl(f)
"""
cross(::typeof(∇),f) = curl(f)
cross(::typeof(∇),f::GridapType) = curl(f)

# Automatic differentiation of functions

function gradient(f::Function)
  function grad_f(x)
    _grad_f(f,x,zero(return_type(f,typeof(x))))
  end
end

function _grad_f(f,x,fx)
  VectorValue(ForwardDiff.gradient(f,get_array(x)))
end

function _grad_f(f,x,fx::VectorValue)
  TensorValue(transpose(ForwardDiff.jacobian(y->get_array(f(y)),get_array(x))))
end

function _grad_f(f,x,fx::MultiValue)
  @notimplemented
end

function divergence(f::Function)
  x -> tr(ForwardDiff.jacobian(y->get_array(f(y)),get_array(x)))
end

function curl(f::Function)
  x -> grad2curl(TensorValue(transpose(ForwardDiff.jacobian(y->get_array(f(y)),get_array(x)))))
end

function laplacian(f::Function)
  function lapl_f(x)
    _lapl_f(f,x,zero(return_type(f,typeof(x))))
  end
end

function _lapl_f(f,x,fx)
  tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(f,y), get_array(x)))
end

function _lapl_f(f,x,fx::VectorValue)
  A = length(x)
  B = length(fx)
  a = ForwardDiff.jacobian(y->transpose(ForwardDiff.jacobian(z->get_array(f(z)),y)), get_array(x))
  tr(ThirdOrderTensorValue{A,A,B}(Tuple(transpose(a))))
end

function _lapl_f(f,x,fx::MultiValue)
  @notimplemented
end

function symmetric_gradient(f::Function)
    x -> symmetric_part(_grad_f(f,x,zero(return_type(f,typeof(x)))))
end


text:
# Computes r_hat from Section 3.1.1 of the paper
function rhat(θ, Ψy)
    r = transpose(Ψy) * θ
    r .= exp.(r)
    r ./ mean(r)
end

# Computes the Hessian of the KLIEP loss
KLIEP_Hessian(θ, Ψy) = StatsBase.cov(Ψy, weights(rhat(θ, Ψy)), 2; corrected=false)

# Computes the sufficient statistics for an Ising mode
#
#    Ψ(x, y) = -1. if x == y
#    Ψ(x, y) =  1. if x != y
#
# Observations are columns of X
#
# Output
#    matrix with columns corresponding to observations
#    rows are sufficients statistics Ψ(xa, xb)
#    the order corresponds to the packed order where
#    lower triangular part is mapped row by row
function Ψising(X)
    m, n = size(X)
    p = div(m * (m - 1), 2)
    out = zeros(Float64, p, n)
    for i = 1:n
        ind = 0
        for row = 2:m
            for col = 1:row-1
                ind += 1
                out[ind, i] = xor(X[col, i], X[row, i]) ? -1. : 1.
            end
        end
    end

    # test that all the sufficient statistics
    # have variance different from zero
    for j=1:p
        if iszero(var(view(out, j, :)))
            @warn "a row of Ψ has zero variance"
        end
    end
    out
end

function _find_supp(k::Int, x)
    supp = findall(!iszero, x)
    pos_k = findfirst(isequal(k), supp)
    if pos_k == nothing
        push!(supp, k)
    else
        supp[pos_k], supp[end] = supp[end], supp[pos_k]
    end
    supp
end

function _find_supp(k::Int, x, y)
    supp1 = findall(!iszero, x)
    supp2 = findall(!iszero, y)
    supp3 = union(supp1, supp2)
    pos_k = findfirst(isequal(k), supp3)
    if pos_k == nothing
        push!(supp3, k)
    else
        supp3[pos_k], supp3[end] = supp3[end], supp3[pos_k]
    end
    supp3
end

# maps (i, j) coordinate of a symmetric matrix to a packed format
# diagonal of the matrix is not included
# upper triangular part is mapped col by col
function trimap(i::Integer, j::Integer)
    if i > j
        trimap(j, i)
    else
        i + div((j - 1) * (j - 2), 2)
    end
end

# mapping from packed storage to (i, j)
function itrimap(k::Integer)
    j = convert(Int, ceil((1. + sqrt(1. + 8. * k)) / 2.))
    i = k - div((j - 1) * (j - 2), 2)
    CartesianIndex(i, j)
end

# vector -> symmetric matrix
function unpack(θ)
    m = length(θ)
    p = convert(Int, ceil((1. + sqrt(1. + 8. * m)) / 2.))
    out = zeros(Float64, p, p)
    for i=1:m
        out[itrimap(i)] = θ[i]
    end
    out + out'
end

# symmetric matrix -> vector
function pack(Θ)
    m = size(Θ, 1)
    p = div(m * (m - 1), 2)
    out = zeros(Float64, p)
    for j = 1:(m-1)
        for i = (j+1):m
            out[trimap(i,j)] = Θ[i,j]
        end
    end
    out
end

####################################
#
# loss En[-θ'Ψx(i)] + log( En[exp(θ'Ψy(i))] )
#
####################################
struct CDKLIEPLoss <: CoordinateDifferentiableFunction
    μx::Vector{Float64}   # mean(Ψx, dims = 2)
    Ψy::Matrix{Float64}
    r::Vector{Float64}    # ny dim vector --- stores θ'Ψy(i)
    p::Int64
end

function CDKLIEPLoss(Ψx::Matrix{Float64}, Ψy::Matrix{Float64})
    (p = size(Ψx, 1)) == size(Ψy, 1) || throw(DimensionMismatch())

    CDKLIEPLoss(vec(mean(Ψx, dims=2)), Ψy, zeros(size(Ψy, 2)), p)
end

CoordinateDescent.numCoordinates(f::CDKLIEPLoss) = f.p

function CoordinateDescent.initialize!(f::CDKLIEPLoss, x::SparseIterate)

    mul!(f.r, transpose(f.Ψy), x)

    nothing
end

function CoordinateDescent.gradient(
    f::CDKLIEPLoss,
    x::SparseIterate,
    j::Int64)

    μx, Ψy, r = f.μx, f.Ψy, f.r

    -μx[j] + mean( exp.(r) .* Ψy[j, :] ) / mean(exp, r)
end

function CoordinateDescent.descendCoordinate!(
    f::CDKLIEPLoss,
    g::Union{ProxL1, ProxZero},
    x::SparseIterate,
    j::Int64)

    μx, Ψy, r = f.μx, f.Ψy, f.r

    mean_exp_r = mean(exp, r)
    mean_exp_r_Ψy = mean( exp.(r) .* Ψy[j, :] )
    mean_exp_r_Ψy2 = mean( exp.(r) .* Ψy[j, :].^2. )

    @inbounds grad = -μx[j] + mean_exp_r_Ψy / mean_exp_r
    H = mean_exp_r_Ψy2 / mean_exp_r - mean_exp_r_Ψy^2 / mean_exp_r^2.

    @inbounds oldVal = x[j]
    @inbounds x[j] -= grad / H
    newVal = cdprox!(g, x, j, 1. / H)
    h = newVal - oldVal

    # update internals
    for i=1:length(r)
        @inbounds r[i] += h*Ψy[j, i]
    end
    h
end


text:
using Test
using NQCDynamics

ω = 3.4
β = 11.2
m = 1
n_beads = 400
s = HarmonicRingPolymer(ω, β, m, n_beads)
configs = (reshape(rand(s), (1,1,n_beads)) for _ in 1:5000)

sim = RingPolymerSimulation(Atoms(1), Harmonic(m=m, ω=ω), n_beads, temperature=1/β)

@testset "Energy expectation" begin
    quantum = ω/2 * coth(ω*β/2)
    estimate = Estimators.@estimate total_energy(sim, configs)
    @test quantum ≈ estimate rtol=1e-2
end

@testset "Radius of gyration expectation" begin
    T = 1/β
    quantum = sqrt(T / ω^2 * (ω / (2T) * coth(ω/2T) - 1))
    estimate = Estimators.@estimate radius_of_gyration(sim, configs)
    @test quantum ≈ estimate[1] rtol=1e-2
end

@testset "select_item" begin
    atoms = 10
    beads = 4
    s = HarmonicRingPolymer.(1:atoms, 1, 1, beads)

    r = NonadiabaticDistributions.select_item(s, 0, (3, atoms, beads))
    @test r isa Array
    @test_throws ErrorException("Sample size does not match distribution.") NonadiabaticDistributions.select_item(s, 0, (3, 1, beads))
    @test_throws ErrorException("Sample size does not match distribution.") NonadiabaticDistributions.select_item(s, 0, (3, 1, 1))
    @test_throws ErrorException("`size` (3, 1) must have 3 dimensions.") NonadiabaticDistributions.select_item(s, 0, (3, 1))
end