text:
mutable struct Chip
    gfx::Vector{UInt8}      # Total amount of pixels: 2048
    key::Vector{UInt8}

    pc::UInt16              # Program counter
    opcode::UInt16          # Current opcode
    I::UInt16               # Index register
    sp::UInt16              # Stack pointer

    V::Vector{UInt8}        # V-regs (V0-VF)
    stack::Vector{UInt16}   # Stack (16 levels)
    memory::Vector{UInt8}   # Memory (size=4k)

    delay_timer::Char       # Delay timer
    sound_timer::Char       # Sound timer

    drawFlag::Bool
end

"Construct Chip-8"
function Chip()
    pc = 0x201      # Program counter stars at 0x200 + 1 (1-based index)
    opcode = 0x00   # Reset current opcode
    I = 0x00        # Reset index register
    sp = 0x00       # Reset stack pointer

    # clear display
    gfx = zeros(UInt8, 2048)

    # clear stack
    stack = zeros(UInt16, 16)
    V = zeros(UInt8, 16)
    key = V

    # clear memory
    memory = zeros(UInt8, 4096)

    # load fontset
    memory[1:80] = [
        0xF0, 0x90, 0x90, 0x90, 0xF0, #0
        0x20, 0x60, 0x20, 0x20, 0x70, #1
        0xF0, 0x10, 0xF0, 0x80, 0xF0, #2
        0xF0, 0x10, 0xF0, 0x10, 0xF0, #3
        0x90, 0x90, 0xF0, 0x10, 0x10, #4
        0xF0, 0x80, 0xF0, 0x10, 0xF0, #5
        0xF0, 0x80, 0xF0, 0x90, 0xF0, #6
        0xF0, 0x10, 0x20, 0x40, 0x40, #7
        0xF0, 0x90, 0xF0, 0x90, 0xF0, #8
        0xF0, 0x90, 0xF0, 0x10, 0xF0, #9
        0xF0, 0x90, 0xF0, 0x90, 0x90, #A
        0xE0, 0x90, 0xE0, 0x90, 0xE0, #B
        0xF0, 0x80, 0x80, 0x80, 0xF0, #C
        0xE0, 0x90, 0x90, 0x90, 0xE0, #D
        0xF0, 0x80, 0xF0, 0x80, 0xF0, #E
        0xF0, 0x80, 0xF0, 0x80, 0x80  #F
    ]

    # reset timers
    delay_timer = 0x00
    sound_timer = 0x00

    Chip(gfx, key,
         pc, opcode, I, sp,
         V, stack, memory,
         delay_timer, sound_timer,
         true)
end

"Read in c8 file"
function loadApplication(c8::Chip, fname::String)
    i = 512 + 1
    open(fname) do f
        while !eof(f)
            buffer = read(f, UInt8)
            c8.memory[i] = buffer

            i += 1
        end
    end
end

"Start the emulation"
function emulateCycle(c8::Chip)
    # fetch the opcode
    c8.opcode = UInt16(c8.memory[c8.pc]) << 8 |
                c8.memory[c8.pc + 1]

    # get the first 4 bits
    first4 = c8.opcode & 0xF000

    # get X and Y
    X = c8.opcode & 0x0f00 >> 8 + 1
    Y = c8.opcode & 0x00f0 >> 4 + 1

    # process opcode
    if first4 == 0x0000
        # Clears the screen.
        if c8.opcode == 0x00e0
            fill!(c8.gfx, 0x00)
        # Returns from a subroutine.
        elseif c8.opcode == 0x00ee
            c8.sp -= 1
            c8.pc = c8.stack[c8.sp + 1]
            c8.pc += 2
        else
            warn("Unknown opcode")
        end
    # 0x1NNN Jumps to address NNN.
    elseif first4 == 0x1000
        c8.pc = c8.opcode & 0x0FFF + 0x0001
    # 0x2NNN Calls subroutine at NNN.
    elseif first4 == 0x2000
        c8.stack[c8.sp+1] = c8.pc
        c8.sp += 1
        c8.pc = c8.opcode & 0x0fff
    # 0x3XNN Skips the next instruction if VX equals NN.
    elseif first4 == 0x3000
        if c8.V[X] == c8.opcode & 0x00ff
            c8.pc += 4
        else
            c8.pc += 2
        end
    # 0x4XNN Skips the next instruction if VX doesn't equal NN.
    elseif first4 == 0x4000
        if c8.V[X] != c8.opcode & 0x00ff
            c8.pc += 4
        else
            c8.pc += 2
        end
    # 0x5XY0 Skips the next instruction if VX equals VY.
    elseif first4 == 0x5000
        if c8.V[X] == c8.V[Y]
            c8.pc += 4
        else
            c8.pc += 2
        end
    # 0x6XNN Sets VX to NN
    elseif first4 == 0x6000
        c8.V[X] = c8.opcode & 0x00ff
        c8.pc += 2
    # 0x7XNN Adds NN to VX
    elseif first4 == 0x7000
        c8.V[X] += c8.opcode & 0x00ff
        c8.pc += 2
    # 0x8XY* Manipulation of VX and VY
    elseif first4 == 0x8000
        last4 = c8.opcode & 0x000f

        # Sets VX to the value of VY.
        if last4 == 0x000
            c8.V[X] = c8.V[Y]
        # Sets VX to VX or VY. (Bitwise OR operation)
        elseif last4 == 0x0001
            c8.V[X] |= c8.V[Y]
        # Sets VX to VX and VY. (Bitwise AND operation)
        elseif last4 == 0x0002
            c8.V[X] &= c8.V[Y]
        # Sets VX to VX xor VY.
        elseif last4 == 0x0003
            c8.V[X] ⊻= c8.V[Y]
        # Adds VY to VX.
        # VF is set to 1 when there's a carry, and to 0 when there isn't.
        elseif last4 == 0x0004
            if c8.V[Y] > (0xff - c8.V[X])
                c8.V[16] = 1
            else
                c8.V[16] = 0
            end
            c8.V[X] += c8.V[Y]
        # VY is subtracted from VX.
        # VF is set to 0 when there's a borrow, and 1 when there isn't.
        elseif last4 == 0x0005
            if c8.V[Y] > c8.V[X]
                c8.V[16] = 0
            else
                c8.V[16] = 1
            end
            c8.V[X] -= c8.V[Y]
        # Shifts VX right by one
        # VF is set to the value of the least significant bit of VX
        # before the shift (modern implementation of 0x8XY6)
        elseif last4 == 0x0006
            c8.V[16] = c8.V[X] & 0x1
            c8.V[X] >>= 1
        # Sets VX to VY minus VX.
        # VF is set to 0 when there's a borrow, and 1 when there isn't.
        elseif last4 == 0x0007
            if c8.V[X] > c8.V[Y]
                c8.V[16] = 0
            else
                c8.V[16] = 1
            end
            c8.V[X] = c8.V[Y] - c8.V[X]
        # Shifts VX left by one
        # VF is set to the value of the most significant bit of VX
        # before the shift (modern implementation of 0x8XYE)
        elseif last4 == 0x000E
            c8.V[16] = c8.V[X] >> 7
            c8.V[X] <<= 1
        else
            warn("Unknown opcode [0x8000]")
        end

        c8.pc += 2
    # 0x9XY0 Skips the next instruction if VX doesn't equal VY.
    elseif first4 == 0x9000
        if c8.V[X] != c8.V[Y]
            c8.pc += 4
        else
            c8.pc += 2
        end
    # Sets I to the address NNN.
    elseif first4 == 0xa000
        I = c8.opcode & 0x0fff
        c8.pc += 2
    # Jumps to the address NNN plus V0.
    elseif first4 == 0xb000
        c8.pc = (c8.opcode & 0x0fff) + c8.V[1]
    # Sets VX to the result of a bitwise and operation on a
    # random number (Typically: 0 to 255) and NN.
    elseif first4 == 0xc000
        c8.V[X] = rand(0x00:0xff) & (c8.opcode * 0x00ff)
        c8.pc += 2
    # Draws a sprite at coordinate (VX, VY) that has a width of 8 pixels
    # and a height of N pixels. Each row of 8 pixels is read as
    # bit-coded starting from memory location I; I value doesn’t
    # change after the execution of this instruction.
    # As described above, VF is set to 1 if any screen pixels are
    # flipped from set to unset when the sprite is drawn, and to 0
    # if that doesn’t happen
    elseif first4 == 0xd000
        height = c8.opcode & 0x000f
        pixel = 0x0000

        V[16] = 0x00

        for yline in 0:(height-1)
            pixel = memory[c8.I + yline + 1]
            for xline in 0:7
                if (pixel & (0x80 >> (xline))) != 0
                    if(gfx[(V[X] + xline + ((V[Y] + yline) * 64)) + 1] == 1)
                        V[16] = 1
                        gfx[V[X] + xline + ((V[Y] + yline) * 64) + 1] ⊻= 1
                    end
                end
            end
        end

        c8.drawFlag = true
        c8.pc += 2
    elseif first4 == 0xe000
        # Skips the next instruction if the key stored in VX is pressed.
        if (c8.opcode & 0x00ff) == 0x009e
            if c8.key[c8.V[X]] != 0
                c8.pc += 4
            else
                c8.pc += 2
            end
        # Skips the next instruction if the key stored in VX isn't pressed
        elseif (c8.opcode & 0x00ff) == 0x00a1
            if c8.key[c8.V[X]] == 0
                c8.pc += 4
            else
                c8.pc += 2
            end
        else
            warn("Unknown opcode")
        end
    elseif first4 == 0xf000
        last4 = c8.opcode * 0x00ff
        # Sets VX to the value of the delay timer.
        if last4 == 0x0007
           c8.V[X] = c8.delay_timer
           c8.pc += 2
        # A key press is awaited, and then stored in VX.
        # (Blocking Operation. All instruction halted until next key event)
        elseif last4 == 0x000a
            keyPresss = false

            for i in 1:16
                if key[i] != 0
                    c8.V[X] = i
                    keyPress = true
                end
            end

            if !keyPress
                return
            end

            c8.pc += 2
        # Sets the delay timer to VX.
        elseif last4 == 0x0015
            c8.delay_timer = c8.V[X]
            c8.pc += 2
        # Sets the sound timer to VX.
        elseif last4 == 0x0018
            c8.sound_timer = c8.V[X]
            c8.pc += 2
        # Adds VX to I
        elseif last4 == 0x001e
            if c8.I + V[X] > 0xfff
                c8.V[16] = 1
            else
                c8.V[16] = 0
            end
            c8.I += c8.V[X]
            c8.pc += 2
        # Sets I to the location of the sprite for the character in VX.
        # Characters 0-F (in hexadecimal) are represented by a 4x5 font.
        elseif last4 == 0x0029
            c8.I = c8.V[X] * 0x5
            c8.pc += 2
        # Stores the Binary-coded decimal representation of VX at the addresses I, I plus 1, and I plus 2
        elseif last4 == 0x0033
            c8.memory[c8.I]     = V[X] / 100
            c8.memory[c8.I + 1] = (V[X] / 10) % 10
            c8.memory[c8.I + 2] = (V[X] % 100) % 10
            c8.pc += 2
        # Stores V0 to VX (including VX) in memory starting at address I.
        # I is increased by 1 for each value written.
        elseif last4 == 0x0055
            for i in 1:X
                c8.memory[I + i] = c8.V[i]
            end

            c8.I += X + 1
            c8.pc += 2
        # Fills V0 to VX (including VX) with values from memory starting at address I.
        # I is increased by 1 for each value written.
        elseif last4 == 0x0065
            for i in 1:X
                c8.V[i] = c8.memory[I + i]
            end

            c8.I += X + 1
            c8.pc += 2
        end
    # temporary for testing
    else
        warn("Unrecognized opcode")
        c8.pc += 2
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3c84ab1c-f504-11ea-25f0-6f2b3cbf472b
using Pkg, DrWatson

# ╔═╡ 3c84dde2-f504-11ea-1b8c-95f180a567f3
begin
	@quickactivate "StatisticalRethinkingTuring"
	using StatisticalRethinking
end

# ╔═╡ f89aef6a-f503-11ea-306e-e983e6719b50
md"## Clip-03-20-24t.jl"

# ╔═╡ 3c85558c-f504-11ea-3eb4-afb3815d4382
md"### snippet 3.20"

# ╔═╡ 3c8ceeee-f504-11ea-1413-eb6c84b6c00e
pdf.(Binomial(2, 0.7), 0:2)

# ╔═╡ 3c906b5a-f504-11ea-3892-835bb1c47540
md"### snippet 3.21"

# ╔═╡ 3ca25e0a-f504-11ea-3c6b-11d90448ccf0
rand(Binomial(2, 0.7))

# ╔═╡ 3ca30f4e-f504-11ea-3c4e-33ee204e71de
md"### snippet 3.22"

# ╔═╡ 3caea962-f504-11ea-2856-3da15f13e8c1
rand(Binomial(2, 0.7), 10)

# ╔═╡ 3caf7ab8-f504-11ea-0d8b-fb1bff09a842
begin
	N = 100000
	d = rand(Binomial(2, 0.7), N);
	[length(filter(e -> e == i, d)) for i in 0:2] / N |> display

	d = rand(Binomial(9, 0.7), N);
	h1 = histogram([filter(e -> e == i, d) for i in 0:9];
	 bins=-0.5:1:9.5, color=:lightblue, leg=false, xticks=0:9, bar_width=0.2)
end

# ╔═╡ 3cbca094-f504-11ea-3dbb-c7d342341fb2
md"## End of clip-03-20-24t.jl"

# ╔═╡ Cell order:
# ╟─f89aef6a-f503-11ea-306e-e983e6719b50
# ╠═3c84ab1c-f504-11ea-25f0-6f2b3cbf472b
# ╠═3c84dde2-f504-11ea-1b8c-95f180a567f3
# ╟─3c85558c-f504-11ea-3eb4-afb3815d4382
# ╠═3c8ceeee-f504-11ea-1413-eb6c84b6c00e
# ╟─3c906b5a-f504-11ea-3892-835bb1c47540
# ╠═3ca25e0a-f504-11ea-3c6b-11d90448ccf0
# ╟─3ca30f4e-f504-11ea-3c4e-33ee204e71de
# ╠═3caea962-f504-11ea-2856-3da15f13e8c1
# ╠═3caf7ab8-f504-11ea-0d8b-fb1bff09a842
# ╟─3cbca094-f504-11ea-3dbb-c7d342341fb2

ext: jl


text:
module ZStackAnalysis

include("./Preprocessing.jl")
include("./Transforms.jl")
using Interpolations
using .Preprocessing
using .Transforms

function interpolate_zstack(zstack, scale)
    scale_x, scale_y, scale_z = scale/minimum(scale)
    nz, nx, ny = size(zstack)
    grids = ((0:nz-1)*scale_z, (0:nx-1)*scale_x, (0:ny-1)*scale_y)
    algorithm = BSpline(Linear())
    itp = interpolate(zstack, algorithm)
    sitp = Interpolations.scale(itp, grids...)
    return sitp(0:(nz-1)*scale_z, 0:(nx-1)*scale_x, 0:(ny-1)*scale_y)
end

export Preprocessing
export interpolate_zstack
export LineFilterTransform, OrientationFilterTransform, line_segment
end


text:
@testset "join" begin
    t1 = NDSparse(Columns(([1,1,2,2], [1,2,1,2])), [1,2,3,4])
    t2 = NDSparse(Columns(([0,2,2,3], [1,1,2,2])), [1,2,3,4])

    j1 = innerjoin(t1,t2)
    j2 = innerjoin(+, t1,t2)

    lj1 = leftjoin(t1,t2)
    lj2 = leftjoin(+, t1,t2)

    mj1 = merge(t1,t2)
    mj2 = merge(t2,t1)
    nn = Ref{Int}()
    nn2 = Ref{Int}()
    for n=1:5
        nn[] = n
        d1 = distribute(t1, n)
        for n2 = 1:5
            nn2[] = n2
            d2 = distribute(t2, n2)
            @test isequal(collect(innerjoin(d1, d2)), j1)
            @test isequal(collect(innerjoin(+, d1, d2)), j2)

            @test isequal(collect(leftjoin(d1, d2)), lj1)
            @test isequal(collect(leftjoin(+, d1, d2)), lj2)

            @test isequal(collect(merge(d1, d2)), mj1)
            @test isequal(collect(merge(d2, d1)), mj2)
        end
    end

    t1 = NDSparse([:msft,:ibm,:ge], [1,3,4], [100,200,150])
    t2 = NDSparse([:ibm,:msft,:aapl,:ibm], [0,0,0,2], [100,99,101,98])

    aj = asofjoin(t1,t2)

    for n=1:5
        d1 = distribute(t1, n)
        for n2 = 1:5
            d2 = distribute(t2, n2)
            @test collect(asofjoin(d1, d2)) == aj
        end
    end
end

@testset "broadcast" begin
    t1 = NDSparse(Columns(([1,2,3,4],[1,1,2,2])), [5,6,7,8])
    t2 = NDSparse(Columns(([0,1,2,4],)),[0,10,100,10000])

    for n=1:4
        for m=1:4
            d1 = distribute(t1, n)
            d2 = distribute(t2, m)

            @test collect(d1 .+ d2) == t1 .+ t2
            @test collect(broadcast(+, d1, d2, dimmap=(0,1))) == broadcast(+, t1, t2, dimmap=(0,1))
        end
    end
end

@testset "join with missing" begin
    y = rand(10)
    z = rand(10)

    # DIndexedTable
    t  = table((x=1:10,   y=y), pkey=:x, chunks=2)
    t2 = table((x=1:2:20, z=z), pkey=:x, chunks=2)

    # DNDSparse
    nd = ndsparse((x=1:10,), (y=y,), chunks=2)
    nd2 = ndsparse((x=1:2:20,), (z=z,), chunks=2)

    @testset "how = :left" begin
        # Missing
        z_left = Union{Float64,Missing}[missing for i in 1:10]
        z_left[1:2:9] = z[1:5]
        t_left = table((x = 1:10, y = y, z = z_left))
        nd_left = ndsparse((x=1:10,), (y=y, z=z_left))
        @test isequal(join(t, t2; how=:left), t_left)
        @test isequal(join(nd, nd2; how=:left), nd_left)

        # DataValue
        z_left2 = [DataValue{Float64}() for i in 1:10]
        z_left2[1:2:9] = z[1:5]
        t_left2 = table((x=1:10, y = y, z = z_left2))
        nd_left2 = ndsparse((x=1:10,), (y=y, z=z_left2))
        @test isequal(join(t, t2, how=:left, missingtype=DataValue), t_left2)
        @test isequal(join(nd, nd2, how=:left, missingtype=DataValue), nd_left2)
    end
end


text:
function assignterminals!(system::CVSystem,old=Dict("a"=>0),restart="no")
    # terminal properties adapted from Danielsen (1998)
    # R2Total = 0.3*mmHgToPa/cm3Tom3;
    R2Total = 0.3*mmHgToPa/cm3Tom3;
    R3Total = 0.21*mmHgToPa/cm3Tom3;
    # R3Total = 0.22*mmHgToPa/cm3Tom3;
    # R4Total = 0.003*mmHgToPa/cm3Tom3;
    # R5Total = 0.01*mmHgToPa/cm3Tom3;
    R4Total = 0.003*mmHgToPa/cm3Tom3;
    R5Total = 0.01*mmHgToPa/cm3Tom3;
    C1Total = 0.01*cm3Tom3/mmHgToPa;
    C2Total = 1.64*cm3Tom3/mmHgToPa;
    C3Total = 1.81*cm3Tom3/mmHgToPa;
    C4Total = 13.24*cm3Tom3/mmHgToPa;
    C5Total = 73.88*cm3Tom3/mmHgToPa;
    L5Total = 5e-5*mmHgToPa/cm3Tom3;
    V01 = 370*cm3Tom3;
    V02 = 370*cm3Tom3;
    V03 = 400*cm3Tom3;
    # V04 = 596*cm3Tom3;
    V04 = 500*cm3Tom3;
    # V05 = 1938*cm3Tom3;
    V05 = 1400*cm3Tom3;

    numterminals = 0;
    numlower = 0;
    numupper = 0;

    for i = 1:length(system.branches.ID)
        if isempty(system.branches.children[i])
            numterminals += 1
            if system.branches.group[i] == "lower"
                numlower += 1
            else
                numupper += 1
            end
        end
    end

    # divide flow between upper and lower body
    if numupper > 0
        lowerflowfraction = 0.7;
    else
        lowerflowfraction = 1;
    end

    # divide flow between systemic veins and vena cava
    venousfractionofR = 0.1;
    venousfractionofL = 0.05;
    venousfractionofC = 0.95;
    venousfractionofV0 = 0.9;

    # divide total peripheral impedance according to desired flow divisions
    R2Lower = numterminals*R2Total;
    R2Upper = numterminals*R2Total;
    R3Lower = numterminals*R3Total;
    R3Upper = numterminals*R3Total;
    R4Lower = numterminals*R4Total;
    R4Upper = numterminals*R4Total;
    C1Lower = C1Total/numterminals;
    C1Upper = C1Total/numterminals;
    C2Lower = C2Total/numterminals;
    C2Upper = C2Total/numterminals;
    C3Lower = C3Total/numterminals;
    C3Upper = C3Total/numterminals;
    C4Lower = C4Total/numterminals;
    C4Upper = C4Total/numterminals;
    V01Lower = V01/numterminals;
    V01Upper = V01/numterminals;
    V02Lower = V02/numterminals;
    V02Upper = V02/numterminals;
    V03Lower = V03/numterminals;
    V03Upper = V03/numterminals;
    V04Lower = V04/numterminals;
    V04Upper = V04/numterminals;

    C5Lower = venousfractionofC*C5Total/numterminals;
    C5Upper = venousfractionofC*C5Total/numterminals;
    V05Lower = venousfractionofV0*V05/numterminals;
    V05Upper = venousfractionofV0*V05/numterminals;

    if numupper > 0
        system.svc.C = 0.25*(1-venousfractionofC)*C5Total;
        system.ivc.C = 0.75*(1-venousfractionofC)*C5Total;
        system.svc.V0 = 0.5*(1-venousfractionofV0)*V05;
        system.ivc.V0 = 0.5*(1-venousfractionofV0)*V05;

        RL = R5Total/lowerflowfraction;
        RU = R5Total/(1-lowerflowfraction);

        system.ivc.R = RL*(venousfractionofR/(1-venousfractionofR) + 1)^-1;
        system.svc.R = RU*(venousfractionofR/(1-venousfractionofR) + 1)^-1;

        Rl = RL - system.ivc.R;
        Ru = RU - system.svc.R;

        R5Lower = numlower*Rl;
        R5Upper = numupper*Ru;

        LL = L5Total/lowerflowfraction;
        LU = L5Total/(1-lowerflowfraction);

        system.ivc.L = LL*(venousfractionofL/(1-venousfractionofL) + 1)^-1;
        system.svc.L = LU*(venousfractionofL/(1-venousfractionofL) + 1)^-1;

        Ll = LL - system.ivc.L;
        Lu = LU - system.svc.L;

        L5Lower = numlower*Ll;
        L5Upper = numupper*Lu;
    else
        system.ivc.C = (1-venousfractionofC)*C5Total;
        system.ivc.V0 = (1-venousfractionofV0)*V05;

        RL = R5Total/lowerflowfraction;
        system.ivc.R = RL*(venousfractionofR/(1-venousfractionofR) + 1)^-1;

        Rl = RL - system.ivc.R;
        R5Lower = numlower*Rl;

        LL = L5Total/lowerflowfraction;
        system.ivc.L = LL*(venousfractionofL/(1-venousfractionofL) + 1)^-1;

        Ll = LL - system.ivc.L;
        L5Lower = numlower*Ll;
    end

    # construct and define terminals
    for i in 1:length(system.branches.ID)
        if isempty(system.branches.children[i])
            system.branches.term[i] = CVModule.ArterialTerminal();
            group = system.branches.group[i];
            if group == "lower"
                if !isempty(system.branches.term[i].C)
                    system.branches.term[i].C[1] += C1Lower;
                else
                    push!(system.branches.term[i].C,C1Lower);
                end
                push!(system.branches.term[i].C,C2Lower);
                push!(system.branches.term[i].C,C3Lower);
                if restart == "no"
                    push!(system.branches.term[i].C,C4Lower);
                    push!(system.branches.term[i].C,C5Lower);
                elseif restart == "yes"
                    temp = old[i]["C"];
                    push!(system.branches.term[i].C,temp[4]);
                    push!(system.branches.term[i].C,temp[5]);
                end
                if !isempty(system.branches.term[i].R)
                    system.branches.term[i].R[1] += (system.solverparams.rho*
                        system.branches.c0[i][end]/system.branches.A0[i][end]);
                else
                    push!(system.branches.term[i].R,system.solverparams.rho*
                        system.branches.c0[i][end]/system.branches.A0[i][end]);
                end
                if restart == "no"
                    push!(system.branches.term[i].R,R2Lower);
                    push!(system.branches.term[i].R,R3Lower);
                elseif restart == "yes"
                    temp = old[i]["R"];
                    push!(system.branches.term[i].R,temp[2]);
                    push!(system.branches.term[i].R,temp[3]);
                end
                push!(system.branches.term[i].R,R4Lower);
                push!(system.branches.term[i].R,R5Lower);
                if !isempty(system.branches.term[i].V0)
                    system.branches.term[i].V0[1] += V01Lower;
                else
                    push!(system.branches.term[i].V0,V01Lower);
                end
                push!(system.branches.term[i].V0,V02Lower);
                push!(system.branches.term[i].V0,V03Lower);
                if restart == "no"
                    push!(system.branches.term[i].V0,V04Lower);
                    push!(system.branches.term[i].V0,V05Lower);
                elseif restart == "yes"
                    temp = old[i]["V0"];
                    push!(system.branches.term[i].V0,temp[4]);
                    push!(system.branches.term[i].V0,temp[5]);
                end
                append!(system.branches.term[i].L,[zeros(4),L5Lower;])
            else
                if !isempty(system.branches.term[i].C)
                    system.branches.term[i].C[1] += C1Upper;
                else
                    push!(system.branches.term[i].C,C1Upper);
                end
                push!(system.branches.term[i].C,C2Upper);
                push!(system.branches.term[i].C,C3Upper);
                if restart == "no"
                    push!(system.branches.term[i].C,C4Upper);
                    push!(system.branches.term[i].C,C5Upper);
                elseif restart == "yes"
                    temp = old[i]["C"];
                    push!(system.branches.term[i].C,temp[4]);
                    push!(system.branches.term[i].C,temp[5]);
                end
                if !isempty(system.branches.term[i].R)
                    system.branches.term[i].R[1] += (system.solverparams.rho*
                        system.branches.c0[i][end]/system.branches.A0[i][end]);
                else
                    push!(system.branches.term[i].R,system.solverparams.rho*
                        system.branches.c0[i][end]/system.branches.A0[i][end]);
                end
                if restart == "no"
                    push!(system.branches.term[i].R,R2Upper);
                    push!(system.branches.term[i].R,R3Upper);
                elseif restart == "yes"
                    temp = old[i]["R"];
                    push!(system.branches.term[i].R,temp[2]);
                    push!(system.branches.term[i].R,temp[3]);
                end
                push!(system.branches.term[i].R,R4Upper);
                push!(system.branches.term[i].R,R5Upper);
                if !isempty(system.branches.term[i].V0)
                    system.branches.term[i].V0[1] += V01Upper;
                else
                    push!(system.branches.term[i].V0,V01Upper);
                end
                push!(system.branches.term[i].V0,V02Upper);
                push!(system.branches.term[i].V0,V03Upper);
                if restart == "no"
                    push!(system.branches.term[i].V0,V04Upper);
                    push!(system.branches.term[i].V0,V05Upper);
                elseif restart == "yes"
                    temp = old[i]["V0"];
                    push!(system.branches.term[i].V0,temp[4]);
                    push!(system.branches.term[i].V0,temp[5]);
                end
                append!(system.branches.term[i].L,[zeros(4),L5Upper;])
            end
        else
            system.branches.term[i] = CVModule.ArterialTerminal();
            if !isempty(system.branches.term[i].C)
                system.branches.term[i].C[1] += NaN;
            else
                push!(system.branches.term[i].C,NaN);
            end
            if !isempty(system.branches.term[i].R)
                system.branches.term[i].R[1] += NaN;
            else
                push!(system.branches.term[i].R,NaN);
            end
            if !isempty(system.branches.term[i].V0)
                system.branches.term[i].V0[1] += NaN;
            else
                push!(system.branches.term[i].V0,NaN);
            end
            push!(system.branches.term[i].L,NaN);
        end
    end

end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using Documenter
using Literate

const literate_dir = joinpath(@__DIR__, "literate")
const generated_dir = joinpath(@__DIR__, "src", "generated")

@info "Loading Catlab.jl"
using Catlab,
  # Must explicitly load all submodules for which docs are needed.
  Catlab.Theories, Catlab.CategoricalAlgebra, Catlab.WiringDiagrams,
  Catlab.Graphics, Catlab.Programs, Catlab.LinearAlgebra

@info "Building Literate.jl docs"

# XXX: Work around old LaTeX distribution in GitHub CI.
if haskey(ENV, "GITHUB_ACTIONS")
    import TikzPictures
    TikzPictures.standaloneWorkaround(true)
end

# Set Literate.jl config if not being compiled on recognized service.
config = Dict{String,String}()
if !(haskey(ENV, "GITHUB_ACTIONS") || haskey(ENV, "GITLAB_CI"))
    config["nbviewer_root_url"] = "https://nbviewer.jupyter.org/github/mu-ziq/cat-ml/blob/gh-pages/dev"
    config["repo_root_url"] = "https://github.com/mu-ziq/cat-ml/blob/master/docs"
end

for (root, dirs, files) in walkdir(literate_dir)
    out_dir = joinpath(generated_dir, relpath(root, literate_dir))
    for file in files
        if last(splitext(file)) == ".jl"
            Literate.markdown(joinpath(root, file), out_dir;
        config=config, documenter=true, credit=false)
            Literate.notebook(joinpath(root, file), out_dir;
        execute=true, documenter=true, credit=false)
        end
    end
end

@info "Building Documenter.jl docs"
makedocs(
  modules=[Catlab],
  format=Documenter.HTML(),
  sitename="CatML",
  doctest=false,
  checkdocs=:none,
  pages=Any[
    "Catlab.jl" => "index.md",
    "Vignettes" => Any[
      "generated/composejl_wiring_diagrams.md",
    ],
  ]
)

@info "Deploying docs"
deploydocs(
  target="build",
  repo="github.com/mu-ziq/CatML.git",
  branch="gh-pages"
)


text:
# Retrieves politicians from Google and inserts/updates them accordingly

function get_politicians(api_key, conn, zipcode)
    url = "https://www.googleapis.com/civicinfo/v2/representatives"
    query = Dict("address" => zipcode, "key" => api_key)
    response = get(url; query = query)
    response_data = Requests.json(response)

    positions = Dict()

    for office in response_data["offices"]
        for index in office["officialIndices"]
            positions[index + 1] = office["name"]
        end
    end

    for (i, politician) in enumerate(response_data["officials"])
        data = Dict("name" => politician["name"], "zipcodes" => string(zipcode))

        if haskey(positions, i)
            if haskey(POSITIONS, positions[i])
                data["position"] = POSITIONS[positions[i]]
            else
                data["position"] = positions[i]
            end
        end

        if haskey(politician, "party")
            party = politician["party"]

            if party == "Democratic"
                party = "Democrat"
            end

            if lowercase(party) != "unknown"
                data["party"] = party
            end
        end

        if haskey(politician, "photoUrl")
            data["photo_url"] = politician["photoUrl"]
        end

        if haskey(politician, "emails")
            data["email"] = politician["emails"][1]
        end

        if haskey(politician, "phones")
            data["phone"] = politician["phones"][1]
        end

        if haskey(politician, "urls")
            data["website"] = politician["urls"][1]
        end

        if haskey(politician, "channels")
            for channel in politician["channels"]
                data[lowercase(channel["type"])] = channel["id"]
            end
        end

        mysql_stmt_prepare(conn, "SELECT * FROM politicians WHERE name=?")
        mysql_execute(conn, [MYSQL_TYPE_STRING], [politician["name"]])

        existing_data = 0

        for row in MySQLRowIterator(conn, [MYSQL_TYPE_STRING], [politician["name"]])
            existing_data = row # there should only ever be one
        end

        if existing_data == 0
            values_str = ""
            keys = []
            types = []
            values = []

            for (j, (key, val)) in enumerate(data)
                values_str = string(values_str, j == length(data) ? "?" : "?, ")
                push!(keys, key)
                push!(types, MYSQL_TYPE_STRING)
                push!(values, val)
            end

            command = string("INSERT INTO politicians (", join(keys, ", "), ") VALUES (", values_str, ")")

            mysql_stmt_prepare(conn, command)
            mysql_execute(conn, types, values)
        else
            new_data = Dict()

            zipcodes = map(y -> parse(Int64, y), split(existing_data[3], ","))

            if !in(zipcode, zipcodes)
                push!(zipcodes, zipcode)
                new_data["zipcodes"] = join(zipcodes, ",")
            end

            check(pos, key) = begin
                if haskey(data, key)
                    if isnull(existing_data[pos])
                        new_data[key] = data[key]
                    elseif data[key] != get(existing_data[pos])
                        new_data[key] = data[key]
                    end
                end
            end

            check(4, "position")
            check(5, "party")
            check(6, "photo_url")
            check(7, "email")
            check(8, "phone")
            check(9, "website")
            check(10, "facebook")
            check(11, "googleplus")
            check(12, "twitter")
            check(13, "youtube")

            if length(new_data) > 0
                keys_values_str = ""
                types = []
                values = []

                for (j, (key, val)) in enumerate(new_data)
                    keys_values_str = string(keys_values_str, key, "=?", j == length(new_data) ? "" : ",")
                    push!(types, MYSQL_TYPE_STRING)
                    push!(values, val)
                end

                command = string("UPDATE politicians SET ", keys_values_str, " WHERE id=", existing_data[1])

                mysql_stmt_prepare(conn, command)
                mysql_execute(conn, types, values)
            end
        end
    end
end


text:
gather(input::AbstractArray{T}, index::CuArray{Int}) where T = gather(cu(Array(input)), index)

function gather(input::CuMatrix{T}, index::CuArray{Int}) where T
    out = CUDA.zeros(T, size(input,1), size(index)...)
    @inbounds for ind = CartesianIndices(index)
        view(out, :, ind) .= view(input, :, index[ind])
    end
    return out
end

function gather_indices(X::CuArray{T}) where T
    Y = gather_indices(Array(X))
    cuY = Dict{T,CuVector}(k => cu(Tuple.(v)) for (k, v) in Y)
    cuY
end

function numerical_cmp(X::CuArray{T}, Y::CuArray) where T
    Z = map((x,y) -> sign(x - y)^2, X, Y)
    Z = map(x -> (one(T) - x)^2, Z)
    Z
end


text:
struct FilterState{T1, T2}
    μ::T1
    Σ::T2
end

abstract type SF end # supertype of all synaptic filters

struct BlockSF{T1, T2} <: SF
    numblocks::Int
    blocksize::Int
    SModel::T1
    OModel::T2
end

struct DiagSF{T1, T2} <: SF
    dim::Int
    SModel::T1
    OModel::T2
end

struct FullSF{T1, T2} <: SF
    dim::Int
    SModel::T1
    OModel::T2
end

struct GradientRule{T1, T2, T3} <: SF
    dim::Int
    SModel::T1
    OModel::T2
    η::T3 # learning rate
end

function FilterState(filter::DiagSF)
    dim = filter.dim
    μ = zeros(dim)
    Σ = zeros(dim)
    return FilterState(μ, Σ)
end

function FilterState(filter::BlockSF)
    numblocks = filter.numblocks
    blocksize = filter.blocksize
    μ = zeros(numblocks * blocksize)
    Σ = zeros(blocksize, blocksize, numblocks)
    return FilterState(μ, Σ)
end

function FilterState(filter::FullSF)
    dim = filter.dim
    μ = zeros(dim)
    Σ = zeros(dim, dim)
    return FilterState(μ, Σ)
end

function FilterState(filter::GradientRule)
    dim = filter.dim
    μ = zeros(dim)
    Σ = filter.η
    return FilterState(μ, Σ)
end

function update!(state::FilterState, filter::SF, obs::NeuronObs, dt)
    SModel = filter.SModel
    OModel = filter.OModel

    τ = SModel.τ
    σs = SModel.σs

    μ = state.μ
    Σ = state.Σ

    g0 = OModel.g0
    β = OModel.β

    x = obs.x
    y = obs.y

    _filter_update!(μ, Σ, τ, σs, g0, β, x, y, dt)
    return state
end

function expterm(μ::AbstractVector{T1}, x::AbstractVector{T2}, v::AbstractVector{T1}) where {T1, T2}
    sum = zero(T1)
    for i in eachindex(μ)
        @inbounds sum += x[i] * (μ[i] + v[i] / 2)
    end
    return sum
end

# update for diagonal filter
function _filter_update!(μ, Σ::AbstractVector, τ, σs, g0, β, x, y, dt)
    α = dt / τ
    α2 = 2*α

    v = β .* Σ .* x
    γ = g0 * dt * exp(β * expterm(μ, x, v))

    μ .+= -μ .* α .+ v .* (y - γ)
    Σ .-= γ .* v .* v .+ (Σ .- σs) .* α2
    return nothing
end

# update for block filter
function _filter_update!(μ, Σ::AbstractArray{T, 3}, τ, σs, g0, β, x, y, dt) where T
    numblocks = size(Σ, 3)
    blocksize = size(Σ, 2)

    α = dt / τ
    α2 = 2*α

    # temporary array for storing the Σ * x for each block
    V = zeros(blocksize, numblocks)
    u = 0.

    @inbounds for i in 1:numblocks
        μ1 = view(μ, (i-1)*blocksize+1:i*blocksize)
        Σ1 = view(Σ, :, :, i)

        x1 = view(x, (i-1)*blocksize+1:i*blocksize)

        v = view(V, :, i)

        v .= β .* (Σ1 * x1)
        u += β * expterm(μ1, x1, v)
    end

    @inbounds for i in 1:numblocks
        μ1 = view(μ, (i-1)*blocksize+1:i*blocksize)
        Σ1 = view(Σ, :, :, i)

        v = view(V, :, i)
        γ = g0 * dt * exp(u)

        for j in 1:size(Σ1, 2)
            @inbounds μ1[j] += -μ1[j] * α + v[j] * (y - γ)
            for i in 1:size(Σ1, 1)
                @inbounds Σ1[i, j] -= γ * v[i] * v[j] + (Σ1[i, j] - σs * (i==j)) * α2
            end
        end
    end
    return nothing
end

# update for full filter
function _filter_update!(μ, Σ::AbstractMatrix, τ, σs, g0, β, x, y, dt)
    α = dt / τ
    α2 = 2*α
    v = β * (Σ * x)
    γ = g0 * dt * exp(β * expterm(μ, x, v))

    # explicit loops for 2.5x speedup and less allocations
    for j in 1:size(Σ, 2)
        @inbounds μ[j] += -μ[j] * α + v[j] * (y - γ)
        for i in 1:size(Σ, 1)
            @inbounds Σ[i, j] -= γ * v[i] * v[j] + (Σ[i, j] - σs * (i==j)) * α2
        end
    end
    return nothing
end

function _filter_update!(μ, Σ::AnyCuMatrix, τ, σs, g0, β, x, y, dt)
    α = dt / τ
    α2 = 2*α
    v = β * (Σ * x)
    γ = g0 * dt * exp(β * (dot(μ, x) + dot(x, v) / 2))

    μ .+= -μ .* α .+ v .* (y - γ)
    view(Σ, diagind(Σ)) .-= σs
    Σ .-= γ .* v * transpose(v) .+ Σ .* α2
    view(Σ, diagind(Σ)) .+= σs
    return nothing
end

# update for gradient rule
function _filter_update!(μ, Σ::Number, τ, σs, g0, β, x, y, dt)
    v = (β * Σ) .* x
    γ = g0 * dt * exp(β * expterm(μ, x, v))

    μ .+= v .* (y - γ)
    return nothing
end


text:
#################### Safe searches ####################

#"""
# proper Markdown comments are not possible
#
#This part creates functions which enable the search for different nodes in the
#tree. It is possible to look for a node via its name, its binary representation
#or to find the root.
#This functionality can be extended by adding more fields to the nodes and the
#meta programmming part here.

# These functions are safe, in the sense that they properly exit if the node is
# not found.
#"""
for (sym, my_type) in [(:binary, :String), (:name, :String), (:root ,:Bool), (:num, :Int64)]
    # extend the list to look for more fields in the node
    @eval function $(Symbol(string("find_by_$sym")))(tree::T, identifier::$my_type)::T  where T<:GeneralNode
        # create each function and make it so it only accepts the correct type
        local all_nodes = post_order(tree) # make sure all_nodes only belongs to this function
        for node in all_nodes
            if node.$sym == identifier
                # return the node if it is found
                return node
            end # if
        end # for
        # the node is not found. Therefore throw an error!
        throw("The node identified by $identifier is not in the tree.")
    end # function
end

#################### Unsafe searches ####################

"""
    find_num(root::T, num::Int64)  where T<:GeneralNode

Find a node by its number. The function assumes that the node is present in the
tree.

Do not use this function if you are unsure whether the node is in the tree at all.

Returns reference to Node.

* `root` : root Node of tree to be searched.

* `num` : number of desired Node.
"""
function find_num(root::T, num::I)::T  where {T<:GeneralNode, I<:Integer}
    po = post_order(root)
    store = Node[]
    found = find_num(root, num, store)
    if length(store) == 0
        throw(ArgumentError("Node not found"))
    else
        return store[1]
    end
end

"""
    find_num(root::T, num::Int64, rn::Vector{T})::Bool  where T<:GeneralNode

Do a post order traversal to find the node corresponding to the `num`.

Returns true if node is found, false otherwise. Desired Node is pushed to rn.

* `root` : root Node of tree to be searched.

* `num` : number of desired Node.

* `rn` : Vector of Nodes; desired Node is pushed to this vector when found.
"""
function find_num(root::T, num::I, rn::Vector{T})::Bool  where {T<:GeneralNode, I<:Integer}
    # if the current node is the correct one store it in rn
    if root.num === num
        push!(rn, root)
        found = true
    else
        found = false
    end

    if !found
        # if the node is not yet found continue
        for child in root.children
            found = find_num(child, num,  rn)
        end
    end # if
    return found
end

"""
    find_binary(root::T, bin::String)::T where T<:GeneralNode

Find a node by its binary representation. The function assumes that the node is
present in the tree.

Do not use this function if you are unsure whether the node is in the tree at all.

Returns a reference to the desired Node.

* `root` : root Node of tree to search.

* `bin` : binary representation of desired Node as a String.
"""
function find_binary(root::T, bin::String)::T where T<:GeneralNode
    rv = root
    for i in split(bin[1:end-1],",")[2:end]
        rv = rv.children[parse(Int64, i)+1]
    end
    rv
end

"""
    find_root(node::Node)::Node

Finds the root of tree indicated by Node.

Returns reference to root Node of the tree.

* `node` : Node in Tree of interest.
"""
function find_root(node::Node)::Node
    while node.root == false
        node = node.mother
    end # while
    return node
end


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings(
    description="This program computes and/or verifies an approximate sum of squares decomposition for element
    `Adj_n + k·Op_n - l·Δ_n ~ Σξ_i*ξ_i ∈ Σ²₂ ISAut(F_n)`.
    Parameters `-n`, `-k` and `-l` must be provided when executing the script. \n",
    preformatted_description=true)
    @add_arg_table! s begin
        "-n"
            help = "Perform computations in automorphism group of the free group on `n` generators."
            arg_type = Int
            required = true
        "-k"
            help = "Compute sum of squares for `Adj_n + k·Op_n ∈ ISAut(F_n)`."
            arg_type = Float64
            required = true
        "--lambda", "-l"
            help = "upper bound on the amount of `Δ_n` contained in the element."
            arg_type = Float64
            required = true
        "--tol"
            help = "solvers numerical tolerance"
            required = false
            default = 1e-12
    end

    return parse_args(s)
end


text:
include("Hiker.jl")
import .Hiker
using Test

@test Hiker.answer() == 42

max_stars_repo_name: cyber-dojo-start-points/julia-test

max_issues_repo_name: cyber-dojo-start-points/julia-test

max_forks_repo_name: cyber-dojo-start-points/julia-test


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_polydata

"""
    vtk_polydata(;kwargs...)
    vtk_polydata(children::Any;kwargs...)
    vtk_polydata(children_maker::Function;kwargs...)

A PolyData component.
PolyData is exposing a vtkPolyData to a downstream filter
It takes the following set of properties:
  - points: [x, y, z, x, y, z, ...],
  - verts: [cellSize, pointId0, pointId1, ..., cellSize, pointId0, ...]
  - lines: [cellSize, pointId0, pointId1, ..., cellSize, pointId0, ...]
  - polys: [cellSize, pointId0, pointId1, ..., cellSize, pointId0, ...]
  - strips: [cellSize, pointId0, pointId1, ..., cellSize, pointId0, ...]
Cell connectivity helper property:
  - connectivity: 'manual', // [manual, points, triangles, strips]
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
- `id` (String; optional): The ID used to identify this component.
- `port` (Real; optional): downstream connection port
- `points` (Array of Reals; optional): xyz coordinates
- `verts` (Array of Reals; optional): verts cells
- `lines` (Array of Reals; optional): lines cells
- `polys` (Array of Reals; optional): polys cells
- `strips` (Array of Reals; optional): strips cells
- `connectivity` (String; optional): Type of connectivity `manual` or implicit such as `points`, `triangles`, `strips`
"""
function vtk_polydata(; kwargs...)
        available_props = Symbol[:children, :id, :port, :points, :verts, :lines, :polys, :strips, :connectivity]
        wild_props = Symbol[]
        return Component("vtk_polydata", "PolyData", "dash_vtk", available_props, wild_props; kwargs...)
end

vtk_polydata(children::Any; kwargs...) = vtk_polydata(;kwargs..., children = children)
vtk_polydata(children_maker::Function; kwargs...) = vtk_polydata(children_maker(); kwargs...)

max_stars_repo_name: ericbdaniels/dash-vtk

max_issues_repo_name: ericbdaniels/dash-vtk

max_forks_repo_name: ericbdaniels/dash-vtk


text:
module ImportKeysightBin

export importkeysightbin

@enum WaveformType begin
	unknownwaveform = 0
	normal = 1
	peakdetect = 2
	average = 3
	logic = 6
end

@enum Unit begin
	unknownunit = 0
	volt = 1
	second = 2
	constant = 3
	ampere = 4
	decibel = 5
	hertz = 6
end

@enum BufferType begin
	unknowbuf = 0
	float32 = 1
	maxfloat = 2
	minfloat = 3
	char = 6
end

"""
    importkeysightbin(io)
    importkeysightbin(filename)

Import waveform data from the binary file format used by _Keysight_
(formerly _Agilent_) oscilloscopes.

The return value is a tuple containing a tuple of the `x` and `y` values
for each channel found in the file, and a dictionary containing metadata
read from the file header.

# Examples
```julia
data = importkeysightbin(io)
ch1, ch2, ch3, ch4, metadata = importkeysightbin(io)
(x1, y1), (x2, y2), (x3, y3), (x4, y4), metadata = importkeysightbin(io)
```
"""
function importkeysightbin(f::IO)
	# File header
	cookie = String(read(f, 2))
	if cookie != "AG"
		error("Unsupported file format")
	end
	meta = Dict{Symbol,Any}()
	meta[:version] = String(read(f, 2))
	meta[:size] = read(f, Int32)
	meta[:nwaveforms] = read(f, Int32)
	meta[:waveforms] = Vector{Dict{Symbol,Any}}(undef, meta[:nwaveforms])
	channels = Vector{Tuple}()

	# Waveforms
	for wf_n in 1:meta[:nwaveforms]
		wf = Dict{Symbol,Any}()
		meta[:waveforms][wf_n] = wf
		wf[:headersize] = read(f, Int32)
		wf[:type] = WaveformType(read(f, Int32))
		wf[:buffers] = read(f, Int32)
		wf[:points] = read(f, Int32)
		wf[:count] = read(f, Int32)
		wf[:xdisplayrange] = read(f, Float32)
		wf[:xdisplayorigin] = read(f, Float64)
		wf[:xincrement] = read(f, Float64)
		wf[:xorigin] = read(f, Float64)
		wf[:xunits] = Unit(read(f, Int32))
		wf[:yunits] = Unit(read(f, Int32))
		wf[:date] = rstrip(String(read(f, 16)), '\0')
		wf[:time] = rstrip(String(read(f, 16)), '\0')
		wf[:frame] = rstrip(String(read(f, 24)), '\0')
		wf[:label] = rstrip(String(read(f, 16)), '\0')
		wf[:timetags] = read(f, Float64)
		wf[:segmentindex] = read(f, UInt32)
		wf[:datasets] = Vector{Dict{Symbol,Any}}(undef, wf[:buffers])
		x = range(start=wf[:xorigin], step=wf[:xincrement], length=wf[:points])
		y = zeros(wf[:points], wf[:buffers])
		push!(channels, (x, y))

		# Data sets
		for ds_n in 1:wf[:buffers]
			ds = Dict{Symbol,Any}()
			wf[:datasets][ds_n] = ds
			ds[:headersize] = read(f, Int32)
			ds[:buffertype] = BufferType(read(f, Int16))
			ds[:bytesperpoint] = read(f, Int16)
			ds[:buffersize] = read(f, Int32)
			y[:,ds_n] = reinterpret(Float32, read(f, ds[:buffersize]))
		end
	end
	(channels..., meta)
end

function importkeysightbin(filename::AbstractString)
	f = open(filename)
	data = importkeysightbin(f)
	close(f)
	return data
end

end


text:
using OffsetArrays
using OffsetArrays: IdOffsetRange, IdentityUnitRange, no_offset_view

@testset "Offset Axes" begin
    oneto = Base.OneTo(2)
    a = Axis(oneto)
    b = OffsetAxis(1, oneto)
    #b = Axis(oneto, IdOffsetRange(oneto, 1))
    c = Axis(OffsetAxis(1, oneto), OffsetAxis(1, oneto))
    #c = Axis(IdOffsetRange(oneto, 1), IdOffsetRange(oneto, 1))
    d = Axis(OffsetAxis(2, oneto), OffsetAxis(1, oneto))
    #d = Axis(IdOffsetRange(oneto, 2), IdOffsetRange(oneto, 1))
    e = SimpleAxis(values(b))
    f = SimpleAxis(2:3)

    @test b[2] == c[2] == d[2]
    @test_throws BoundsError b[1]

    # `keys(v) == OneTo(10)` so `firstindex(keys(v)) == 1` but `firstindex(values(v)) == 2`
    @test b[2] == 2
    # `keys(v) == IdOffsetRange(oneto, 1)` so `firstindex(keys(v)) == 2` and `firstindex(values(v)) == 2`
    @test c[2] == 2
    # `keys(v) == IdOffsetRange(oneto, 2)` so `firstindex(keys(v)) == 3` and `firstindex(values(v)) == 2`
    @test d[isequal(3)] == 2
    @test e[isequal(2)] == 2

    @test values(b[in(1:2)]) == 2:3
    @test values(c[in(2:3)]) == 2:3
    @test values(d[in(3:4)]) == 2:3

    oneto = Base.OneTo(2)
    a = Axis(oneto)
    b = Axis(["a", "b"], IdOffsetRange(oneto, 1))
    @test b["a"] == 2
    @test b[["a", "b"]] == [2, 3]

end

#= TODO
@testset "Single-entry arrays in dims 0:5" begin
    for n = 0:5
        for z in (OffsetArray(ones(Int,ntuple(d->1,n)), ntuple(x->x-1,n)),
                  fill!(OffsetArray{Float64}(undef, ntuple(x->x:x, n)), 1),
                  fill!(OffsetArray{Float64}(undef, ntuple(x->x:x, n)...), 1),
                  fill!(OffsetArray{Float64,n}(undef, ntuple(x->x:x, n)), 1),
                  fill!(OffsetArray{Float64,n}(undef, ntuple(x->x:x, n)...), 1))
            @test length(LinearIndices(AxisArray(z))) == 1
            @test axes(z) == ntuple(x->x:x, n)
            @test z[1] == 1
        end
    end
    a0 = reshape([3])
    a = OffsetArray(a0)
    @test axes(a) == ()
    @test ndims(a) == 0
    @test a[] == 3
end
=#

@testset "Offset range construction" begin
    r = -2:5
    y = AxisArray(OffsetArray(r, r))
    @test axes(y) == (r,)
    y = AxisArray(OffsetArray(r, (r,)))
    @test axes(y) == (r,)
end

@testset "Traits" begin
    A0 = [1 3; 2 4]
    A = AxisArray(OffsetArray(A0, (-1,2)))                   # IndexLinear
    S = AxisArray(OffsetArray(view(A0, 1:2, 1:2), (-1,2)))   # IndexCartesian
    @test axes(A) == axes(S) == (0:1, 3:4)
    @test size(A) == size(A0)
    @test size(A, 1) == size(A0, 1)
    @test length(A) == length(A0)
    @test A == OffsetArray(A0, 0:1, 3:4)
end

@testset "Scalar indexing" begin
    A0 = [1 3; 2 4]
    A = AxisArray(OffsetArray(A0, (-1,2)));
    S = AxisArray(OffsetArray(view(A0, 1:2, 1:2), (-1,2)));

    @test @inferred(A[0,3]) == @inferred(A[0,3,1]) == @inferred(A[1]) == @inferred(S[0,3]) == @inferred(S[0,3,1]) == @inferred(S[1]) == 1
    @test A[1,3] == A[1,3,1] == A[2] == S[1,3] == S[1,3,1] == S[2] == 2
    @test A[0,4] == A[0,4,1] == A[3] == S[0,4] == S[0,4,1] == S[3] == 3
    @test A[1,4] == A[1,4,1] == A[4] == S[1,4] == S[1,4,1] == S[4] == 4
    @test @inbounds(A[0,3]) == @inbounds(A[0,3,1]) == @inbounds(A[1]) == @inbounds(S[0,3]) == @inbounds(S[0,3,1]) == @inbounds(S[1]) == 1
    @test @inbounds(A[1,3]) == @inbounds(A[1,3,1]) == @inbounds(A[2]) == @inbounds(S[1,3]) == @inbounds(S[1,3,1]) == @inbounds(S[2]) == 2
    @test @inbounds(A[0,4]) == @inbounds(A[0,4,1]) == @inbounds(A[3]) == @inbounds(S[0,4]) == @inbounds(S[0,4,1]) == @inbounds(S[3]) == 3
    @test @inbounds(A[1,4]) == @inbounds(A[1,4,1]) == @inbounds(A[4]) == @inbounds(S[1,4]) == @inbounds(S[1,4,1]) == @inbounds(S[4]) == 4
    @test_throws BoundsError A[1,1]
    @test_throws BoundsError S[1,1]
    @test_throws BoundsError A[0,3,2]
    @test_throws BoundsError A[0,3,0]
    Ac = copy(A)
    Ac[0,3] = 10
    @test Ac[0,3] == 10
    Ac[0,3,1] = 11
    @test Ac[0,3] == 11
    @inbounds Ac[0,3,1] = 12
    @test Ac[0,3] == 12

    y = AxisArray(OffsetArray{Float64}(undef, -1:1, -7:7, -128:512, -5:5, -1:1, -3:3, -2:2, -1:1))
    y[-1,-7,-128,-5,-1,-3,-2,-1] = 14
    y[-1,-7,-128,-5,-1,-3,-2,-1] += 5
    @test y[-1,-7,-128,-5,-1,-3,-2,-1] == 19
end

@testset "CartesianAxisCore" begin
    A0 = [1 3; 2 4]
    A = AxisArray(OffsetArray(A0, (-1,2)));
    S = AxisArray(OffsetArray(view(A0, 1:2, 1:2), (-1,2)));

    @test A[CartesianIndex((0,3))] == S[CartesianIndex((0,3))] == 1
    @test A[CartesianIndex((0,3)),1] == S[CartesianIndex((0,3)),1] == 1
    @test @inbounds(A[CartesianIndex((0,3))]) == @inbounds(S[CartesianIndex((0,3))]) == 1
    @test @inbounds(A[CartesianIndex((0,3)),1]) == @inbounds(S[CartesianIndex((0,3)),1]) == 1
    @test_throws BoundsError A[CartesianIndex(1,1)]
    @test_throws BoundsError A[CartesianIndex(1,1),0]
    @test_throws BoundsError A[CartesianIndex(1,1),2]
    @test_throws BoundsError S[CartesianIndex(1,1)]
    @test_throws BoundsError S[CartesianIndex(1,1),0]
    @test_throws BoundsError S[CartesianIndex(1,1),2]
    @test eachindex(A) == 1:4
    @test eachindex(S) == CartesianIndices(IdentityUnitRange.((0:1,3:4)))
end

@testset "Vector indexing" begin
    A0 = [1 3; 2 4]
    A_1 = OffsetArray(A0, (-1,2))
    S_1 = OffsetArray(view(A0, 1:2, 1:2), (-1,2))

    A = AxisArray(OffsetArray(A0, (-1,2)))
    S = AxisArray(OffsetArray(view(A0, 1:2, 1:2), (-1,2)))

    s = S[:, 3]
    s == S_1[:,3]
    o = OffsetArray([1,2], (parent(A).offsets[1],))
    @test A[:, 3] == S[:, 3] == OffsetArray([1,2], (parent(A).offsets[1],))
    @test A[:, 4] == S[:, 4] == OffsetArray([3,4], (parent(A).offsets[1],))
    @test_throws BoundsError A[:, 1]
    @test_throws BoundsError S[:, 1]
    @test A[0, :] == S[0, :] == OffsetArray([1,3], (parent(A).offsets[2],))
    @test A[1, :] == S[1, :] == OffsetArray([2,4], (parent(A).offsets[2],))
    @test_throws BoundsError A[2, :]
    @test_throws BoundsError S[2, :]
    @test A[0:1, 3] == S[0:1, 3] == [1,2]
    @test A[[1,0], 3] == S[[1,0], 3] == [2,1]
    @test A[0, 3:4] == S[0, 3:4] == [1,3]
    @test A[1, [4,3]] == S[1, [4,3]] == [4,2]
    @test A[:, :] == S[:, :] == A
end

#= TODO
@testset "Vector indexing with offset ranges" begin
    r = AxisArray(OffsetArray(8:10, -1:1))
    r1 = r[0:1]
    @test r1 == 9:10
    r1 = (8:10)[AxisArray(OffsetArray(1:2, -5:-4))]
    @test axes(r1) == (IdentityUnitRange(-5:-4),)
    @test parent(r1) == 8:9
    r1 = AxisArray(OffsetArray(8:10, -1:1))[OffsetArray(0:1, -5:-4)]
    @test axes(r1) == (IdentityUnitRange(-5:-4),)
    @test parent(r1) == 9:10
end

@testset "view" begin
    A0 = [1 3; 2 4]
    A = AxisArray(OffsetArray(A0, (-1,2)))

    S = view(A, :, 3)
    @test S == AxisArray(OffsetArray([1,2], (A.offsets[1],)))
    @test S[0] == 1
    @test S[1] == 2
    @test_throws BoundsError S[2]
    @test axes(S) == (IdentityUnitRange(0:1),)
    S = view(A, 0, :)
    @test S == AxisArray(OffsetArray([1,3], (A.offsets[2],)))
    @test S[3] == 1
    @test S[4] == 3
    @test_throws BoundsError S[1]
    @test axes(S) == (IdentityUnitRange(3:4),)
    S = view(A, 0:0, 4)
    @test S == [3]
    @test S[1] == 3
    @test_throws BoundsError S[0]
    @test axes(S) == (Base.OneTo(1),)
    S = view(A, 1, 3:4)
    @test S == [2,4]
    @test S[1] == 2
    @test S[2] == 4
    @test_throws BoundsError S[3]
    @test axes(S) == (Base.OneTo(2),)
    S = view(A, :, :)
    @test S == A
    @test S[0,3] == S[1] == 1
    @test S[1,3] == S[2] == 2
    @test S[0,4] == S[3] == 3
    @test S[1,4] == S[4] == 4
    @test_throws BoundsError S[1,1]
    @test axes(S) == IdentityUnitRange.((0:1, 3:4))
    S = view(A, axes(A)...)
    @test S == A
    @test S[0,3] == S[1] == 1
    @test S[1,3] == S[2] == 2
    @test S[0,4] == S[3] == 3
    @test S[1,4] == S[4] == 4
    @test_throws BoundsError S[1,1]
    @test axes(S) == OffsetArrays.IdOffsetRange.((0:1, 3:4))
    # issue 100
    S = view(A, axes(A, 1), 3)
    @test S == A[:, 3]
    @test S[0] == 1
    @test S[1] == 2
    @test_throws BoundsError S[length(S)]
    @test axes(S) == (OffsetArrays.IdOffsetRange(0:1), )
    # issue 100
    S = view(A, 1, axes(A, 2))
    @test S == A[1, :]
    @test S[3] == 2
    @test S[4] == 4
    @test_throws BoundsError S[1]
    @test axes(S) == (OffsetArrays.IdOffsetRange(3:4), )

    A0 = collect(reshape(1:24, 2, 3, 4))
    A = OffsetArray(A0, (-1,2,1))
    S = view(A, axes(A, 1), 3:4, axes(A, 3))
    @test S == A[:, 3:4, :]
    @test S[0, 1, 2] == A[0, 3, 2]
    @test S[0, 2, 2] == A[0, 4, 2]
    @test S[1, 1, 2] == A[1, 3, 2]
    @test axes(S) == (OffsetArrays.IdOffsetRange(0:1), Base.OneTo(2), OffsetArrays.IdOffsetRange(2:5))
end

@testset "similar" begin
    A0 = [1 3; 2 4]
    A = AxisArray(OffsetArray(A0, (-1, 2)))

    B = similar(A, Float32)
    @test isa(B, AxisArray(OffsetArray{Float32,2}))
    @test axes(B) == axes(A)
    B = similar(A, (3,4))
    @test isa(B, Array{Int,2})
    @test size(B) == (3,4)
    @test axes(B) == (Base.OneTo(3), Base.OneTo(4))
    B = similar(A, (-3:3,1:4))
    @test isa(B, OffsetArray{Int,2})
    @test axes(B) == IdentityUnitRange.((-3:3, 1:4))
    B = similar(parent(A), (-3:3,1:4))
    @test isa(B, OffsetArray{Int,2})
end

@testset "reshape" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))

    B = reshape(A0, -10:-9, 9:10)
    @test isa(B, OffsetArray{Int,2})
    @test parent(B) === A0
    @test axes(B) == IdentityUnitRange.((-10:-9, 9:10))
    B = reshape(A, -10:-9, 9:10)
    @test isa(B, OffsetArray{Int,2})
    @test pointer(parent(B)) === pointer(A0)
    @test axes(B) == IdentityUnitRange.((-10:-9, 9:10))
    b = reshape(A, -7:-4)
    @test axes(b) == (IdentityUnitRange(-7:-4),)
    @test isa(parent(b), Vector{Int})
    @test pointer(parent(b)) === pointer(parent(A))
    @test parent(b) == A0[:]
    a = OffsetArray(rand(3,3,3), -1:1, 0:2, 3:5)
    # Offset axes are required for reshape(::OffsetArray, ::Val) support
    b = reshape(a, Val(2))
    @test isa(b, OffsetArray{Float64,2})
    @test pointer(parent(b)) === pointer(parent(a))
    @test axes(b) == IdentityUnitRange.((-1:1, 1:9))
    b = reshape(a, Val(4))
    @test isa(b, OffsetArray{Float64,4})
    @test pointer(parent(b)) === pointer(parent(a))
    @test axes(b) == (axes(a)..., IdentityUnitRange(1:1))

    @test reshape(OffsetArray(-1:0, -1:0), :, 1) == reshape(-1:0, 2, 1)
    @test reshape(OffsetArray(-1:2, -1:2), -2:-1, :) == reshape(-1:2, -2:-1, 2)
end

@testset "AxisCore with OffsetArray axes" begin
    A0 = [1 3; 2 4]

    i1 = OffsetArray([2,1], (-5,))
    i1 = OffsetArray([2,1], -5)
    b = A0[i1, 1]
    @test axes(b) == (IdentityUnitRange(-4:-3),)
    @test b[-4] == 2
    @test b[-3] == 1
    b = A0[1,i1]
    @test axes(b) == (IdentityUnitRange(-4:-3),)
    @test b[-4] == 3
    @test b[-3] == 1
    v = view(A0, i1, 1)
    @test axes(v) == (IdentityUnitRange(-4:-3),)
    v = view(A0, 1:1, i1)
    @test axes(v) == (Base.OneTo(1), IdentityUnitRange(-4:-3))

    for r in (1:10, 1:1:10, StepRangeLen(1, 1, 10), LinRange(1, 10, 10))
        for s in (IdentityUnitRange(2:3), OffsetArray(2:3, 2:3))
            @test axes(r[s]) == axes(s)
        end
    end
end

@testset "logical indexing" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))

    @test A[A .> 2] == [3,4]
end

@testset "copyto!" begin
    a = AxisArray(OffsetArray{Int}(undef, (-3:-1,)))
    fill!(a, -1)
    copyto!(a, (1,2))   # non-array iterables
    @test a[-3] == 1
    @test a[-2] == 2
    @test a[-1] == -1
    fill!(a, -1)
    copyto!(a, -2, (1,2))
    @test a[-3] == -1
    @test a[-2] == 1
    @test a[-1] == 2
    @test_throws BoundsError copyto!(a, 1, (1,2))
    fill!(a, -1)
    copyto!(a, -2, (1,2,3), 2)
    @test a[-3] == -1
    @test a[-2] == 2
    @test a[-1] == 3
    @test_throws BoundsError copyto!(a, -2, (1,2,3), 1)
    fill!(a, -1)
    copyto!(a, -2, (1,2,3), 1, 2)
    @test a[-3] == -1
    @test a[-2] == 1
    @test a[-1] == 2

    b = 1:2    # copy between AbstractArrays
    bo = AxisArray(OffsetArray(1:2, (-3,)))
    if VERSION < v"1.5-"
        @test_throws BoundsError copyto!(a, b)
        fill!(a, -1)
        copyto!(a, bo)
        @test a[-3] == -1
        @test a[-2] == 1
        @test a[-1] == 2
    else
        # the behavior of copyto! is corrected as the documentation says "first n element"
        # https://github.com/JuliaLang/julia/pull/34049
        fill!(a, -1)
        copyto!(a, bo)
        @test a[-3] == 1
        @test a[-2] == 2
        @test a[-1] == -1
    end
    fill!(a, -1)
    copyto!(a, -2, bo)
    @test a[-3] == -1
    @test a[-2] == 1
    @test a[-1] == 2
    @test_throws BoundsError copyto!(a, -4, bo)
    @test_throws BoundsError copyto!(a, -1, bo)
    fill!(a, -1)
    copyto!(a, -3, b, 2)
    @test a[-3] == 2
    @test a[-2] == a[-1] == -1
    @test_throws BoundsError copyto!(a, -3, b, 1, 4)
    am = OffsetArray{Int}(undef, (1:1, 7:9))  # for testing linear indexing
    fill!(am, -1)
    copyto!(am, b)
    @test am[1] == 1
    @test am[2] == 2
    @test am[3] == -1
    @test am[1,7] == 1
    @test am[1,8] == 2
    @test am[1,9] == -1
end

@testset "map" begin
    am = AxisArray(OffsetArray{Int}(undef, (1:1, 7:9)))  # for testing linear indexing
    fill!(am, -1)
    copyto!(am, 1:2)

    dest = similar(am)
    map!(+, dest, am, am)
    @test dest[1,7] == 2
    @test dest[1,8] == 4
    @test dest[1,9] == -2
end

@testset "reductions" begin
    A = AxisArray(OffsetArray(rand(4,4), (-3,5)))
    @test maximum(A) == maximum(parent(A))
    @test minimum(A) == minimum(parent(A))
    @test extrema(A) == extrema(parent(A))
    C = similar(A)
    cumsum!(C, A, dims = 1)
    @test parent(C) == cumsum(parent(A), dims = 1)
    @test parent(cumsum(A, dims = 1)) == cumsum(parent(A), dims = 1)
    cumsum!(C, A, dims = 2)
    @test parent(C) == cumsum(parent(A), dims = 2)
    R = similar(A, (1:1, 6:9))
    maximum!(R, A)
    @test parent(R) == maximum(parent(A), dims = 1)
    R = similar(A, (-2:1, 1:1))
    maximum!(R, A)
    @test parent(R) == maximum(parent(A), dims = 2)
    amin, iamin = findmin(A)
    pmin, ipmin = findmin(parent(A))
    @test amin == pmin
    @test A[iamin] == amin
    @test amin == parent(A)[ipmin]
    amax, iamax = findmax(A)
    pmax, ipmax = findmax(parent(A))
    @test amax == pmax
    @test A[iamax] == amax
    @test amax == parent(A)[ipmax]

    amin, amax = extrema(parent(A))
    @test clamp.(A, (amax+amin)/2, amax) == OffsetArray(clamp.(parent(A), (amax+amin)/2, amax), axes(A))
end

# v  = OffsetArray([1,1e100,1,-1e100], (-3,))*1000
# v2 = OffsetArray([1,-1e100,1,1e100], (5,))*1000
# @test isa(v, OffsetArray)
# cv  = OffsetArray([1,1e100,1e100,2], (-3,))*1000
# cv2 = OffsetArray([1,-1e100,-1e100,2], (5,))*1000
# @test isequal(cumsum_kbn(v), cv)
# @test isequal(cumsum_kbn(v2), cv2)
# @test isequal(sum_kbn(v), sum_kbn(parent(v)))

@testset "Collections" begin
    A = OffsetArray(rand(4,4), (-3,5))

    @test unique(A, dims=1) == OffsetArray(parent(A), 0, first(axes(A, 2)) - 1)
    @test unique(A, dims=2) == OffsetArray(parent(A), first(axes(A, 1)) - 1, 0)
    v = OffsetArray(rand(8), (-2,))
    @test sort(v) == OffsetArray(sort(parent(v)), v.offsets)
    @test sortslices(A; dims=1) == OffsetArray(sortslices(parent(A); dims=1), A.offsets)
    @test sortslices(A; dims=2) == OffsetArray(sortslices(parent(A); dims=2), A.offsets)
    @test sort(A, dims = 1) == OffsetArray(sort(parent(A), dims = 1), A.offsets)
    @test sort(A, dims = 2) == OffsetArray(sort(parent(A), dims = 2), A.offsets)

    @test mapslices(v->sort(v), A, dims = 1) == OffsetArray(mapslices(v->sort(v), parent(A), dims = 1), A.offsets)
    @test mapslices(v->sort(v), A, dims = 2) == OffsetArray(mapslices(v->sort(v), parent(A), dims = 2), A.offsets)
end

@testset "rot/reverse" begin
    A = OffsetArray(rand(4,4), (-3,5))

    @test rotl90(A) == OffsetArray(rotl90(parent(A)), A.offsets[[2,1]])
    @test rotr90(A) == OffsetArray(rotr90(parent(A)), A.offsets[[2,1]])
    @test reverse(A, dims = 1) == OffsetArray(reverse(parent(A), dims = 1), A.offsets)
    @test reverse(A, dims = 2) == OffsetArray(reverse(parent(A), dims = 2), A.offsets)
end

@testset "fill" begin
    B = fill(5, 1:3, -1:1)
    @test axes(B) == (1:3,-1:1)
    @test all(B.==5)
end

@testset "broadcasting" begin
    A = OffsetArray(rand(4,4), (-3,5))

    @test A.+1 == OffsetArray(parent(A).+1, A.offsets)
    @test 2*A == OffsetArray(2*parent(A), A.offsets)
    @test A+A == OffsetArray(parent(A)+parent(A), A.offsets)
    @test A.*A == OffsetArray(parent(A).*parent(A), A.offsets)
end

@testset "@inbounds" begin
    a = OffsetArray(zeros(7), -3:3)
    unsafe_fill!(x) = @inbounds(for i in axes(x,1); x[i] = i; end)
    function unsafe_sum(x)
        s = zero(eltype(x))
        @inbounds for i in axes(x,1)
            s += x[i]
        end
        s
    end
    unsafe_fill!(a)
    for i = -3:3
        @test a[i] == i
    end
    @test unsafe_sum(a) == 0
end

@testset "Resizing OffsetVectors" begin
    local a = OffsetVector(rand(5),-3)
    axes(a,1) == -2:2
    length(a) == 5
    resize!(a,3)
    length(a) == 3
    axes(a,1) == -2:0
    @test_throws ArgumentError resize!(a,-3)
end

####
#### type defined for testing no_offset_view
####

struct NegativeArray{T,N,S <: AbstractArray{T,N}} <: AbstractArray{T,N}
    parent::S
end

# Note: this defines the axes-of-the-axes to be OneTo.
# In general this isn't recommended, because
#    positionof(A, i, j, ...) == map(getindex, axes(A), (i, j, ...))
# is quite desirable, and this requires that the axes be "identity" ranges, i.e.,
# `r[i] == i`.
# Nevertheless it's useful to test this on a "broken" implementation
# to make sure we still get the right answer.
Base.axes(A::NegativeArray) = map(n -> (-n):(-1), size(A.parent))

Base.size(A::NegativeArray) = size(A.parent)

function Base.getindex(A::NegativeArray{T,N}, I::Vararg{Int,N}) where {T,N}
    getindex(A.parent, (I .+ size(A.parent) .+ 1)...)
end

@testset "no offset view" begin
    # OffsetArray fallback
    A = randn(3, 3)
    O1 = OffsetArray(A, -1:1, 0:2)
    O2 = OffsetArray(O1, -2:0, -3:(-1))
    @test no_offset_view(O2) ≡ A

    # generic fallback
    A = collect(reshape(1:12, 3, 4))
    N = NegativeArray(A)
    @test N[-3, -4] == 1
    V = no_offset_view(N)
    @test collect(V) == A

    # bidirectional
    B = BidirectionalVector([1, 2, 3])
    pushfirst!(B, 0)
    OB = OffsetArrays.no_offset_view(B)
    @test axes(OB, 1) == 1:4
    @test collect(OB) == 0:3
end

@testset "no nesting" begin
    A = randn(2, 3)
    x = A[2, 2]
    O1 = AxisArray(OffsetArray(A, -1:0, -1:1))
    O2 = AxisArray(OffsetArray(O1, 0:1, 0:2))
    @test parent(O1) ≡ parent(O2)
    @test eltype(O1) ≡ eltype(O2)
    O2[1, 1] = x + 1            # just a sanity check
    @test A[2, 2] == x + 1
end

@testset "mutating functions for OffsetVector" begin
    # push!
    o = AxisArray(OffsetVector(Int[], -1))
    @test push!(o) === o
    @test axes(o, 1) == 0:-1
    @test push!(o, 1) === o
    @test axes(o, 1) == 0:0
    @test o[end] == 1
    @test push!(o, 2, 3) === o
    @test axes(o, 1) == 0:2
    @test o[end-1:end] == [2, 3]
    # pop!
    o = AxisArray(OffsetVector([1, 2, 3], -1))
    @test pop!(o) == 3
    @test axes(o, 1) == 0:1
    # empty!
    o = AxisArray(OffsetVector([1, 2, 3], -1))
    @test empty!(o) === o
    @test axes(o, 1) == 0:-1
end

=#

@testset "iteration" begin
    A0 = [1 3; 2 4]
    A = AxisArray(OffsetArray(A0, (-1, 2)))

    let a
        for (a,d) in zip(A, A0)
            @test a == d
        end
    end
end


text:
using StochasticPrograms
using JuMP
using GLPK
using Ipopt
using Distributed
using Distributions
using LinearAlgebra
using Logging
using Test

import StochasticPrograms: probability, expected

include("decisions/decisions.jl")
TestDecisionVariable.runtests()
TestDecisionConstraint.runtests()
TestDecisionObjective.runtests()
TestSolve.runtests()
include("problems/problem_load.jl")
include("functional_tests.jl")
include("solver_tests.jl")
include("distributed/distributed_tests.jl")


text:
s = """
function f(s,t,m)
s.x = m[1,1] * t.x + m[1,2] * t.y + m[1,3] * t.z
s.y = m[2,1] * t.x + m[2,2] * t.y + m[2,3] * t.z
s.z = m[3,1] * t.x + m[3,2] * t.y + m[3,3] * t.z
return(s)
end
"""

msgs = lintstr(s)
@test isempty(msgs)

# can you spot the error?
s = """
function f(s,t,m)
s.x = m[1,1] * t.x + m[1,2] * t.y + m[1,3] * t.z
s.y = m[2,1] * t.x + m[2,2] * t.y + m[2,3] * t.z
s.z = m[3,1] * t.x + m[3,2] * t.x + m[3,3] * t.z
return(s)
end
"""

ctx = LintContext()
filter!(i->!(i==Lint.LintIgnore(:W651, "")), ctx.ignore)
msgs = lintstr(s, ctx)
@test length(msgs) == 1
@test msgs[1].code == :W651
@test occursin("looks different", msgs[1].message)


text:
module VisualGrowthDynamics

import Base.Iterators: product
import StaticArrays: SVector
import GrowthDynamics: TumorConfigurations, Lattices
import .TumorConfigurations: TumorConfiguration
import .Lattices: HexagonalLattice, CubicLattice, neighbors
import LightGraphs: neighborhood

using GeometryBasics
using Colors

import AbstractPlotting: Plot, plot!, default_theme, to_value
using Makie

include("ColorFunctions.jl")

function default_theme(scene, ::Type{<: Plot(TumorConfiguration{HexagonalLattice{Int64}})})
	Attributes(
		markersize=1,
		size=(800,800)
	)
end

function plot!(p::Plot(TumorConfiguration{HexagonalLattice{Int64}}))
	state = to_value(p[1])
	L = state.lattice.Na

    hex_vertices = map([
			# Vertices of a hex with flat edge down
			Point2f0(1,0),
		    Point2f0(cos(pi/3),-sin(pi/3)),
		    Point2f0(-cos(pi/3),-sin(pi/3)),
		    Point2f0(-1,0),
		    Point2f0(-cos(pi/3),sin(pi/3)),
		    Point2f0(cos(pi/3),sin(pi/3))
			]/L/Float32(cos(pi/6))) do x
				# Rotation by π/6 --> pointy edge down
				[Point2f0(cos(pi/6),sin(pi/6)) Point2f0(-sin(pi/6),cos(pi/6))]*x
			end

	hex_mesh=GeometryBasics.mesh(hex_vertices)
	off = (sqrt(5)-sqrt(3))/2

	hex_positions = reshape(Point2f0[Point2f0(
		2*(q+ifelse(r%2==0,0.,1/2))/L-1.0,
		2*(r+1/2+off/2*ifelse(r==0,0.,-r))/L-1.0)
		 for (r,q) in Iterators.product(0:L-1,0:L-1)], L^2)

	meshscatter!(p, hex_positions, marker=hex_mesh, markersize=1,
	 	color=lift(to_value(p[:color]), Node(state)),
		shading=false, raw=false)

	p
end

## Cubic lattice
function default_theme(scene::SceneLike, ::Type{<: Plot(TumorConfiguration{CubicLattice{Int64}})})
	Attributes(
		markersize=1,
		resolution=(800,800),
		genotype=1
	)
end

function AbstractPlotting.plot!(p::Plot(TumorConfiguration{CubicLattice{Int64}}))
	global myp = p
	@show typeof(p)
	state = to_value(p[1])
	L = state.lattice.Na

	geno_filter = state.meta.genotypes[
		neighborhood(state.Phylogeny, to_value(p[:genotype]), typemax(Int64); dir=:in)
		]

	flt = BitArray( x in geno_filter for x in state.lattice.data)

	if haskey(p.attributes, :plane)
		P = to_value(p[:plane])
		#B = Lattices.intersect_lattice_with_plane(state.lattice, P)
		indices = product(1:L, 1:L, 1:L)
		B = BitArray(Lattices.euclidean_dist(SVector{3}(I), P) <= 1/2 for I in indices)
		#@show size(B)
		flt = flt .& B

		arrows!(fill(Point3f0(P.p), 3), 10.0*[P.u, P.v, P.w], arrowsize=0.1, arrowlength=10, linewidth=3)
	end

	if all(x->x==false, flt)
		@info "No cells to draw."
	end

	positions = Point3f0[Point3f0(l,m,n)
				 for (l,m,n) in Iterators.product(0:L-1,0:L-1,0:L-1)][flt]
    positions = reshape(positions, length(positions))

	color = (to_value(p[:color])(state))[reshape(flt, length(flt))]

    cube = GeometryBasics.HyperRectangle(0,0,0, 1,1,1)
	meshscatter!(p, positions, marker=cube, markersize=1,
	 	 color=color, raw=false)

	p
end

end # module


text:
const ObjectOccupancyArray{O, N} = ObjectIndexableArray{O, BitArray{N}, Bool, N}

function ObjectOccupancyArray(objects::Tuple{Vararg{AbstractObject}}, dims::Integer...)
    grid = falses(length(objects), dims...)
    return ObjectOccupancyArray{typeof(objects), length(dims) + 1}(grid)
end


text:
function compute_nlines(meth)
    lines = collect(map(x->x.line,
        filter(x->isa(x,LineNumberNode),
        Base.uncompressed_ast(meth.lambda_template))))
    isempty(lines) ? 1 : maximum(lines) - meth.line
end

const filemap = Dict{Symbol, Vector{Method}}()
function register_meth(meth)
    !haskey(filemap, meth.file) && (filemap[meth.file] = Vector{Method}())
    push!(filemap[meth.file], meth)
end
function initial_sweep()
    const visited_modules = Set{Core.Module}()
    const workqueue = Set{Core.Module}()
    push!(workqueue, Main)
    push!(visited_modules, Main)
    while !isempty(workqueue)
        mod = pop!(workqueue)
        for sym in names(mod, true)
            !isdefined(mod, sym) && continue
            Base.isdeprecated(mod, sym) && continue
            x = getfield(mod, sym)
            if isa(x, Core.Module)
                !(x in visited_modules) && push!(workqueue, x)
                push!(visited_modules, x)
                continue
            else
                mt = methods(x)
                for m in mt
                    register_meth(m)
                end
            end
        end
    end
    for (k,v) in filemap
        sort!(v, by = m->m.line)
    end
end

type MethSource <: LocationSource
    bp::Breakpoint
    fT::Type
    meth_predicate
    spec_predicate
    function MethSource(bp::Breakpoint,fT::Type, meth_predicate=meth->true, spec_predicate=spec->true)
        !haskey(TracedTypes, fT) && (TracedTypes[fT] = Set{MethSource}())
        this = new(bp,fT,meth_predicate,spec_predicate)
        push!(TracedTypes[fT], this)
        finalizer(this,function (this)
            pop!(TracedTypes[this.fT], this)
            if isempty(TracedTypes[this.fT])
                delete!(TracedMethods, this.fT)
            end
        end)
        this
    end
end

function Base.show(io::IO, source::MethSource)
    print(io,"Any matching method added to ",source.fT)
end
const TracedTypes = Dict{Type,Set{MethSource}}()

type FileLineSource <: LocationSource
    bp::Breakpoint
    fname::Symbol
    line::Int
    function FileLineSource(bp, fname, line)
        this = new(bp, fname, line)
        push!(FLBPs, this)
        this
    end
end
function Base.show(io::IO, source::FileLineSource)
    print(io,"Any method reaching ",source.fname,":",source.line)
end
const FLBPs = Vector{FileLineSource}()

function fire(s::FileLineSource, meth)
    add_meth_to_bp!(s.bp, meth)
end

function fire(s::MethSource, meth)
    s.meth_predicate(meth) || return
    add_meth_to_bp!(s.bp, meth, s.spec_predicate)
end

function newmeth_tracer(x::Ptr{Void})
    meth = unsafe_pointer_to_objref(x)::Method
    fT = meth.lambda_template.specTypes.parameters[1]
    if haskey(TracedTypes, fT)
        for source in TracedTypes[fT]
            fire(source, meth)
        end
    end
    for source in FLBPs
        contains(string(meth.file), string(source.fname)) || continue
        if meth.line <= source.line <= meth.line + compute_nlines(meth)
            fire(source, meth)
        end
    end
    register_meth(meth)
    nothing
end

global did_arm_breakfile = false
function arm_breakfile()
    global did_arm_breakfile
    did_arm_breakfile && return
    did_arm_breakfile = true
    initial_sweep()
end

function methods_for_line(meths, line)
    ret = Vector{Method}()
    for meth in meths
        meth.line > line && break
        if line < meth.line + compute_nlines(meth)
            push!(ret, meth)
        end
    end
    ret
end

function breakpoint(file::AbstractString, line::Int)
    arm_breakfile()
    bp = Breakpoint()
    found = false
    for (fname, meths) in filemap
        contains(string(fname), file) || continue
        found = true
        for meth in methods_for_line(meths, line)
            add_meth_to_bp!(bp, meth)
        end
    end
    if !found
        warn("No file $file found in loaded packages or included files.")
    end
    unshift!(bp.sources, FileLineSource(bp, file, line))
    bp
end


text:
using Batched, BenchmarkTools

A = BatchedMatrix(rand(10, 10, 30, 5))
@benchmark 0.1 * A

@profiler for i in 1:4000
    0.1 * A
end

B = rand(10, 10, 30, 5)

@benchmark 0.1 * B

@profiler for i in 1:4000
    0.1 * B
end

t = Broadcast.broadcasted(*, 0.1, B)

@which copy(t)

@which copyto!(similar(B), t)

convert(Broadcast.Broadcast)

Base.isconcretetype(Float64)


text:
"""
    boundaries(domain::AbstractDomain)

Return an index set giving the boundaries of the domain.
"""
boundaries(domain::AbstractDomain)

"""
    area(domain::AbstractDomain)

Return the area of the domain.
"""
area(domain::AbstractDomain)

"""
    boundary_parameterization(t, domain::AbstractDomain, i::Integer)

Compute a parameterization of boundary number `i`.

The parameterization goes from t = 0 to t = 1. The enumeration of
the boundaries depends on the type of domain used.
"""
boundary_parameterization(t, domain::AbstractDomain, i::Integer)

"""
    boundary_points(domain::AbstractDomain, i::Integer, n::Integer; distribution)

Return `n` points from boundary number `i`, together with information
about which boundary they are from.

The enumeration of the boundary depends on the type of domain used.
The distribution of the points depends on `distribution`, the default
is `:chebyshev` for planar domains and `:linear` for spherical ones.
possible options depends on the domain.
"""
boundary_points(domain::AbstractDomain, i::Integer, n::Integer; distribution = :chebyshev)

"""
    boundary_points(domain::AbstractDomain, n::Integer)

Return `n` points taken from all boundaries of the domain, together
with information about which boundary they are from.
"""
function boundary_points(
    domain::AbstractDomain,
    n::Integer;
    distribution = ifelse(domain isa AbstractPlanarDomain, :chebyshev, :linear),
)
    m = length(boundaries(domain))
    res = [
        boundary_points(domain, i, div(n, m) + ifelse(n % m >= i, 1, 0); distribution)
        for i in boundaries(domain)
    ]
    return vcat(getindex.(res, 1)...), vcat(getindex.(res, 2)...)
end

"""
    boundary_points(
        domain::AbstractDomain,
        u::AbstractEigenfunction,
        N::Integer,
        n::Integer;
        distribution,
)

Return `n` points from boundaries of the domain on which the
eigenfunction is not identically zero, together with information about
which boundary they are from.

The boundaries where the egienfunction is not identically equal to
zero on is determined by `active_boundaries(domain, u)`.

`N` should be set to the number of expansion terms used for `u`. It's
not used in the default implementation but some type of eigenfunctions
need access to this to compute good boundary points.
"""
function boundary_points(
    domain::AbstractDomain,
    u::AbstractEigenfunction,
    ::Integer,
    n::Integer;
    distribution = ifelse(domain isa AbstractPlanarDomain, :chebyshev, :linear),
)
    active = active_boundaries(domain, u)
    m = length(active)
    res = [
        boundary_points(
            domain,
            active[i],
            div(n, m) + ifelse(n % m >= i, 1, 0);
            distribution,
        ) for i in eachindex(active)
    ]
    return vcat(getindex.(res, 1)...), vcat(getindex.(res, 2)...)
end

function boundary_points(
    domain::AbstractDomain,
    u::CombinedEigenfunction,
    ::Integer,
    n::Integer;
    distribution = ifelse(domain isa AbstractPlanarDomain, :root_exponential, :linear),
)
    active = active_boundaries(domain, u)
    m = length(active)
    per_boundary = fill(div(n, m), m) + [n % m >= i for i = 1:m]
    res = [
        getindex.(
            boundary_points(
                domain,
                active[i],
                ifelse(
                    active[i] ∈ u.even_boundaries,
                    2per_boundary[i] - 1,
                    per_boundary[i],
                );
                distribution,
            ),
            Ref(1:per_boundary[i]),
        ) for i in eachindex(active)
    ]
    return vcat(getindex.(res, 1)...), vcat(getindex.(res, 2)...)
end

"""
    interior_points(domain::AbstractDomain, n::Integer; rng = MersenneTwister(42))

Return n points taken uniformly at random from the interior of the
domain.
"""
interior_points(domain::AbstractDomain, n::Integer; rng = MersenneTwister(42))


text:
#!/usr/bin/env julia

using Gadfly
using RDatasets
iris = dataset("datasets", "iris")

p = plot(iris, x=:SepalLength, y=:SepalWidth, Geom.point);

img = SVG("iris_plot.svg", 14cm, 8cm)
draw(img, p)

# Look at this, coloring by species.
#p2 = plot(iris, x=:SepalLength, y=:SepalWidth, color=:Species, Geom.point);
#img2 = SVG("iris_plot2.svg", 14cm, 8cm)
#draw(img2, p2)

# Egads.
# ./hello.jl  42.91s

# Same time for 2 plots, though.
# ./hello.jl  49.89s

# Load iris data alone.
# ./hello.jl  15.18s
# With using Gadfly
# hello.jl  20.15s
# entire first plot.
# 48.74
# So plotting is about 28 seconds.

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
"""Validates the contents of generators."""
function validate(gen::GenClasses)::Bool
    is_valid = true

    for field in fieldnames(typeof(gen))
        generators = getfield(gen, field)
        if !isnothing(generators)
            if !validate_devices(generators)
                is_valid = false
            end
        end
    end

    @debug "GenClasses validation" is_valid
    return is_valid
end

"""Validates the contents of a Generator."""
function validate(generator::Generator)::Bool
    is_valid = true

    @debug "Generator validation" generator.name is_valid
    return is_valid
end


text:
using Distributions
#Using the NonCentral T Distribution
#Use an F-test to calculate power for a two-sided test
function CRT_clustersize(n1, rho, d, alpha, power_level)
    n2 = 3
    while n2 < 100
        ncp = d/(sqrt(4*(1+(n1-1)*rho)/(n1*n2)))
        #println(ncp)
        crit = quantile(FDist(1,n2-2),1-alpha)
        power_calc = 1 - cdf(NoncentralF(1,n2-2,ncp^2),crit)
        #println("p",power_calc)
        if power_calc > power_level
            cluster_num = n2
            return cluster_num
         end
        n2 +=1
    end
end

max_stars_repo_path: Cluster Randomized Trial/CRTPower.jl

max_stars_repo_name: dshehtanian/julia-learn

max_issues_repo_path: Cluster Randomized Trial/CRTPower.jl

max_issues_repo_name: dshehtanian/julia-learn

max_forks_repo_path: Cluster Randomized Trial/CRTPower.jl

max_forks_repo_name: dshehtanian/julia-learn


text:
###############################################################################
#
# Soil optical parameters
#
###############################################################################
"""
    mutable struct SoilOpticals{FT}

A struct of soil optical parameters

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef mutable struct SoilOpticals{FT}
    # these arrays should match the WL in wls
    "Shortwave albedo"
    albedo_SW    ::Array{FT,1}
    "Shortwave albedo for SIF"
    albedo_SW_SIF::Array{FT,1}
    "Shortwave Emissivity"
    emsvty_SW    ::Array{FT,1}
    "Longwave albedo"
    albedo_LW    ::Array{FT,1}
    "Soil surface temperature `[K]`"
    soil_skinT   ::FT
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
twoj(j) = twice(Integer,j)
function twojp1(j)
	t = twoj(j)
	t + one(t)
end

function filledelements(j)
	N = twojp1(j)
	Np = (N+1) >> 1
	return (1 + N - Np)*Np
end

struct WignerdMatrixContainer{T} <: AbstractMatrix{T}
	j :: HalfInt
	v :: Vector{T}

	function WignerdMatrixContainer{T}(j::HalfInt, v::Vector{T}) where {T}

		j >= zero(j) || throw(ArgumentError("j must be ≥ 0"))

		if length(v) < filledelements(j)
			throw(ArgumentError("array provided is not large enough"*
			" to store all values. Need an array of "*
			"length $(filledelements(j)) whereas the array provided "*
			"has length = $(length(v))"))
		end

		new{T}(j, v)
	end
end

function WignerdMatrixContainer{T}(::UndefInitializer, j) where {T}
	jh = HalfInt(j)
	N = filledelements(jh)
	v = Vector{T}(undef, N)
	WignerdMatrixContainer{T}(jh, v)
end

Base.size(d::WignerdMatrixContainer) = (twojp1(d.j), twojp1(d.j))
Base.fill!(d::WignerdMatrixContainer, val) = fill!(d.v, val)
Base.IndexStyle(::Type{<:WignerdMatrixContainer}) = IndexCartesian()

function Base.:(==)(d1::WignerdMatrixContainer, d2::WignerdMatrixContainer)
	d1.j == d2.j || return false
	N = filledelements(d1.j)
	for i = 1:N
		d1.v[i] == d2.v[i] || return false
	end
	return true
end

function flatind(j, m_ind, n_ind)
	# This is m-major, so m increases faster than n
	# Store only the left triangular quadrant
	N = twojp1(j)
	indskip = (2 + N - n_ind)*(n_ind - 1)
	indskip + (m_ind - n_ind) + 1
end

function flatind_phase(j, m_ind, n_ind)
	N = twojp1(j)
	Np = (N + 1) >> 1

	neg = (-1)^(m_ind - n_ind)
	phase = one(neg)

	# The left quadrant (m,n <= 0, n <= m <= -n) is stored
	# We evaluate the others using symmetry relations
	if n_ind > Np && N - n_ind < m_ind <= n_ind
		# right
		# (m, n) → (-m, -n)
		n_ind, m_ind = N - n_ind + 1, N - m_ind + 1
		phase = neg
	elseif m_ind < n_ind
		# top
		# (m, n) → (n, m)
		m_ind, n_ind = n_ind, m_ind
		phase = neg
	elseif n_ind < m_ind && m_ind + n_ind > N + 1
		# bottom
		# (m, n) → (-n ,-m)
		m_ind, n_ind = N - n_ind + 1, N - m_ind + 1
	end

	ind = flatind(j, m_ind, n_ind)

	return ind, phase
end

@propagate_inbounds function Base.getindex(d::WignerdMatrixContainer, m_ind::Int, n_ind::Int, I...)
	@boundscheck checkbounds(d, m_ind, n_ind, I...)
	ind, phase = flatind_phase(d.j, m_ind, n_ind)
	d.v[ind]*phase
end

@propagate_inbounds function Base.setindex!(d::WignerdMatrixContainer, val, m_ind::Int, n_ind::Int, I...)
	@boundscheck checkbounds(d, m_ind, n_ind, I...)
	ind, phase = flatind_phase(d.j, m_ind, n_ind)
	d.v[ind] = val*phase
	d
end

abstract type AbstractWignerMatrix{T} <: AbstractHalfIntegerMatrix{T} end

Base.size(w::AbstractWignerMatrix) = size(parent(w))
Base.axes(w::AbstractWignerMatrix) = axes(parent(w))
Base.axes(w::AbstractWignerMatrix, d) = axes(parent(w), d)

"""
	WignerdMatrix(j, beta::Real)

Compute the Wigner d-matrix in the z-y-z convention for the
angular momentum `j` and the angle `beta`.
The resulting matrix has a size of `(2j+1, 2j+1)` and axes `(-j:j, -j:j)`.
The angular momentum `j` needs to be either an `Integer` or a half-integer.
The matrix elements are real in this convention.

# Examples

```jldoctest
julia> d = WignerdMatrix(0.5, π)
2×2 WignerdMatrix{Float64} for j = 1/2 and beta = π with indices -1/2:1/2×-1/2:1/2:
  0.0  1.0
 -1.0  0.0

julia> d[0.5, 0.5]
0.0
```
"""
struct WignerdMatrix{T, B<:Real} <: AbstractWignerMatrix{T}
	beta :: B
	dj :: SpinMatrix{T,WignerdMatrixContainer{T}}
end

Base.parent(d::WignerdMatrix) = d.dj
Base.collect(d::WignerdMatrix) = collect(parent(d))

@propagate_inbounds function Base.getindex(d::WignerdMatrix{T}, m::Real, n::Real) where {T}
	parent(d)[m,n]
end

@propagate_inbounds function Base.getindex(d::WignerdMatrix{T, <:Union{ScaledPi,Irrational{:π}}}, m::Real, n::Real) where {T}
	@boundscheck checkbounds(d, m, n)
	wignerdmatrixelement(T, sphericaldegree(d), m, n, d.beta)
end

#= Linear indexing has to be done in a roundabout way through Cartesian indexing
to avoid undefined references.
=#
@inline function Base.getindex(d::WignerdMatrix, i::Int)
	@boundscheck checkbounds(d, i)
	mn = Tuple(@inbounds eachindex(d)[i])
	@inbounds d[mn...]
end

function Base.:(==)(d1::WignerdMatrix, d2::WignerdMatrix)
	d1.beta == d2.beta && parent(d1.dj) == parent(d2.dj)
end

function Base.:(==)(d1::WignerdMatrix{<:Any,<:ScaledPi}, d2::WignerdMatrix{<:Any,<:ScaledPi})
	(d1.beta == d2.beta && sphericaldegree(d1) == sphericaldegree(d2)) || return false

	j = sphericaldegree(d1)
	for n = -j:zero(j)
		m = nonzerofirstindex(d1.beta, n)
		d1[m,n] == d2[m,n] || return false
	end
	return true
end

"""
	WignerDMatrix(j, α::Real, β::Real, γ::Real)

Construct the Wigner D-matrix in the z-y-z convention for the angular momentum `j`
and the Euler angles `(α, β, γ)`. The angular momentum `j` needs to be an `Integer`
or a half-integer. The resulting matrix has a size of `(2j+1, 2j+1)` and axes
`(-j:j, -j:j)`.

# Examples

```jldoctest
julia> D = WignerDMatrix(0.5, 0, 0, 0)
2×2 WignerDMatrix{Complex{Float64}} for j = 1/2, alpha = 0, beta = 0 and gamma = 0 with indices -1/2:1/2×-1/2:1/2:
 1.0+0.0im  0.0-0.0im
 0.0+0.0im  1.0+0.0im

julia> D[-0.5, 0.5]
-0.0 + 0.0im
```
"""
struct WignerDMatrix{T<:Complex,A<:Real,W<:WignerdMatrix,G<:Real} <: AbstractWignerMatrix{T}
	alpha :: A
	dj :: W
	gamma :: G
end

function WignerDMatrix{T}(α::A, dj::W, γ::G) where {W<:WignerdMatrix,A<:Real,G<:Real,T<:Complex}
	WignerDMatrix{T,A,W,G}(α, dj, γ)
end

function WignerDMatrix(α::A, dj::WignerdMatrix{R}, γ::G) where {R,A<:Real,G<:Real}
	Tphase = Complex{promote_type(HalfInt,promote_type(A,G))}
	T = promote_type(Tphase, R)
	WignerDMatrix{T}(α, dj, γ)
end

Base.parent(D::WignerDMatrix) = D.dj

WignerDphase(m, α::Real) = cis(-m * α)
WignerDphase(m, α::SpecialAngle) = cis_special(-m, α)
WignerDphase(m, α::ZeroRadians) = 1 # return an integer to avoid complex multiplication

@inline function Base.getindex(D::WignerDMatrix{T}, m::Real, n::Real) where {T}
	@boundscheck checkbounds(D, m, n)
	dmn = @inbounds D.dj[m,n]
	val = dmn * WignerDphase(m, D.alpha) * WignerDphase(n, D.gamma)
	convert(T, val)
end

# Linear indexing forces Cartesian indexing
@inline function Base.getindex(D::WignerDMatrix, i::Int)
	@boundscheck checkbounds(D, i)
	mn = Tuple(@inbounds eachindex(D)[i])
	@inbounds D[mn...]
end

function Base.:(==)(d1::WignerDMatrix, d2::WignerDMatrix)
	d1.alpha == d2.alpha && d1.gamma == d2.gamma && d1.dj == d2.dj
end

eulerangles(d::WignerdMatrix) = (zero(d.beta), d.beta, zero(d.beta))
eulerangles(D::WignerDMatrix) = (D.alpha, D.dj.beta, D.gamma)

sphericaldegree(d::WignerdMatrixContainer) = d.j
sphericaldegree(d::SpinMatrix) = d.j
sphericaldegree(d::Union{WignerdMatrix, WignerDMatrix}) = sphericaldegree(d.dj)

Base.IndexStyle(d::Type{<:AbstractWignerMatrix}) = IndexCartesian()

function Base.similar(d::WignerdMatrix, ::Type{T}) where {T}
	jh = HalfInt(sphericaldegree(d))
	w = WignerdMatrixContainer{T}(undef, jh)
	s = SpinMatrix(w, jh)
	WignerdMatrix(d.beta, s)
end

function Base.similar(D::WignerDMatrix, ::Type{T}) where {T}
	d = similar(D.dj)
	WignerDMatrix{T}(D.alpha, d, D.gamma)
end

HalfIntegerArrays.SpinMatrix(d::AbstractWignerMatrix) = SpinMatrix(collect(d))

Base.dataids(d::WignerdMatrixContainer) = Base.dataids(d.v)
Base.dataids(d::AbstractWignerMatrix) = Base.dataids(parent(d))
function Broadcast.broadcast_unalias(dest::AbstractWignerMatrix, src::AbstractWignerMatrix)
	parent(dest) === parent(src) ? src : Broadcast.unalias(dest, src)
end

# Display

function Base.showarg(io::IO, d::WignerdMatrix, toplevel)
	j = sphericaldegree(d)
	if toplevel
		print(io,"WignerdMatrix{$(eltype(d))} for j = ", j, " and beta = ", d.beta)
	else
		print(io,"WignerdMatrix{$(eltype(d))}(",j,", ",d.beta,")")
	end
end
function Base.showarg(io::IO, D::WignerDMatrix, toplevel)
	α, β, γ = eulerangles(D)
	j = sphericaldegree(D)
	if toplevel
		print(io,"WignerDMatrix{$(eltype(D))} for j = ", j,
			", alpha = ", α,", beta = ", β, " and gamma = ", γ)
	else
		print(io,"WignerDMatrix{$(eltype(D))}(", j,
			", ", α, ", ", β, ", ", γ,")")
	end
end

function Base.replace_in_print_matrix(d::WignerdMatrix{<:Any, <:ScaledPi}, m::HalfInteger, n::HalfInteger, s::AbstractString)
    nonzerocheck(d.beta, m, n, sphericaldegree(d)) ? s : Base.replace_with_centered_mark(s)
end

function Base.replace_in_print_matrix(D::WignerDMatrix{<:Any, <:Any, <:WignerdMatrix{<:Any, <:ScaledPi}},
	m::HalfInteger, n::HalfInteger, s::AbstractString)

	_,β,_ = eulerangles(D)
	nonzerocheck(β, m, n, sphericaldegree(D)) ? s : Base.replace_with_centered_mark(s)
end


text:
using Revise
using LinearAlgebra
using BenchmarkTools
using Strided
using TensorOperations
using Tullio
using LoopVectorization

sizes = ceil.(Int, 2 .^(2:1.5:20))

function benchmark_sum(sizes)
    times = zeros(length(sizes), 3)
    for (i, s) in enumerate(sizes)
        A = randn(Float64, s)
        times[i, 1] = @belapsed sum($A)
        Strided.disable_threads()
        times[i, 2] = @belapsed @strided sum($A)
        Strided.enable_threads()
        times[i, 3] = @belapsed @strided sum($A)
        println("step $i: size $s => times = $(times[i, :])")
    end
    return times
end

function benchmark_permute(sizes, p = (4,3,2,1))
    times = zeros(length(sizes), 4)
    for (i, s) in enumerate(sizes)
        A = randn(Float64, s .* one.(p))
        B = similar(A)
        times[i, 1] = @belapsed copy!($B, $A)
        times[i, 2] = @belapsed permutedims!($B, $A, $p)
        Strided.disable_threads()
        times[i, 3] = @belapsed @strided permutedims!($B, $A, $p)
        Strided.enable_threads()
        times[i, 4] = @belapsed @strided permutedims!($B, $A, $p)
        println("step $i: size $s => times = $(times[i, :])")
    end
    return times
end
permute_times1 = benchmark_permute(sizes, (4,3,2,1))
permute_times2 = benchmark_permute(sizes, (2,3,4,1))
permute_times3 = benchmark_permute(sizes, (3,4,1,2))

function benchmark_mul(sizesm, sizesk = sizesm, sizesn = sizesm)
    N = Threads.nthreads()
    @assert length(sizesm) == length(sizesk) == length(sizesn)
    times = zeros(length(sizesm), 4)
    @inbounds for i = 1:length(sizesm)
        m = sizesm[i]
        k = sizesk[i]
        n = sizesn[i]
        A = randn(Float64, (m,k))
        B = randn(Float64, (k,n))
        C = randn(Float64, (m,n))

        BLAS.set_num_threads(1) # base case: single-threaded blas
        times[i, 1] = @belapsed mul!($C,$A,$B)
        BLAS.set_num_threads(N) # multithreaded blas
        times[i, 2] = @belapsed mul!($C,$A,$B)
        Strided.disable_threaded_mul() # same, except for small overhead from strided
        times[i, 3] = @belapsed @strided mul!($C,$A,$B)
        BLAS.set_num_threads(1) # single-threaded blas with strided multithreading
        Strided.enable_threaded_mul()
        times[i, 4] = @belapsed @strided mul!($C,$A,$B)
        println("step $i: sizes $((m,k,n)) => times = $(times[i, :])")
    end
    return times
end

function tensorcontraction!(wEnv, hamAB,hamBA,rhoBA,rhoAB,w,v,u)
    @tensor wEnv[-1,-2,-3] =
    	hamAB[7,8,-1,9]*rhoBA[4,3,-3,2]*conj(w[7,5,4])*u[9,10,-2,11]*conj(u[8,10,5,6])*v[1,11,2]*conj(v[1,6,3]) +
    	hamBA[1,2,3,4]*rhoBA[10,7,-3,6]*conj(w[-1,11,10])*u[3,4,-2,8]*conj(u[1,2,11,9])*v[5,8,6]*conj(v[5,9,7]) +
    	hamAB[5,7,3,1]*rhoBA[10,9,-3,8]*conj(w[-1,11,10])*u[4,3,-2,2]*conj(u[4,5,11,6])*v[1,2,8]*conj(v[7,6,9]) +
    	hamBA[3,7,2,-1]*rhoAB[5,6,4,-3]*v[2,1,4]*conj(v[3,1,5])*conj(w[7,-2,6])
    return wEnv
end

function benchmark_tensorcontraction(sizes)
    N = Threads.nthreads()
    times = zeros(length(sizes), 5)
    @inbounds for (i, s) in enumerate(sizes)
        hamAB = randn(Float64, (s,s,s,s))
        hamBA = randn(Float64, (s,s,s,s))
        rhoAB = randn(Float64, (s,s,s,s))
        rhoBA = randn(Float64, (s,s,s,s))
        v = randn(Float64, (s,s,s))
        w = randn(Float64, (s,s,s))
        u = randn(Float64, (s,s,s,s))
        wEnv = randn(Float64, (s,s,s))

        tensorcontraction!(wEnv, hamAB,hamBA,rhoBA,rhoAB,w,v,u)

        BLAS.set_num_threads(1)
        Strided.disable_threads()
        Strided.disable_threaded_mul()
        times[i,1] = @belapsed tensorcontraction!($wEnv,$hamAB,$hamBA,$rhoBA,$rhoAB,$w,$v,$u)

        BLAS.set_num_threads(1)
        Strided.enable_threads()
        Strided.disable_threaded_mul()
        times[i,2] = @belapsed tensorcontraction!($wEnv,$hamAB,$hamBA,$rhoBA,$rhoAB,$w,$v,$u)

        BLAS.set_num_threads(N)
        Strided.disable_threads()
        Strided.disable_threaded_mul()
        times[i,3] = @belapsed tensorcontraction!($wEnv,$hamAB,$hamBA,$rhoBA,$rhoAB,$w,$v,$u)

        BLAS.set_num_threads(N)
        Strided.enable_threads()
        Strided.disable_threaded_mul()
        times[i,4] = @belapsed tensorcontraction!($wEnv,$hamAB,$hamBA,$rhoBA,$rhoAB,$w,$v,$u)

        BLAS.set_num_threads(1)
        Strided.enable_threads()
        Strided.enable_threaded_mul()
        times[i,5] = @belapsed tensorcontraction!($wEnv,$hamAB,$hamBA,$rhoBA,$rhoAB,$w,$v,$u)

        println("step $i: size $s => times = $(times[i, :])")
    end
    return times
end