text:
@testset "strel" begin
    @testset "offset to mask" begin
        se = [CartesianIndex(-1), CartesianIndex(1)]
        se_mask = @inferred strel(Bool, se)
        @test axes(se_mask) == (-1:1,)
        @test se_mask == centered(Bool[1, 1, 1])

        se = [CartesianIndex(-1, -1), CartesianIndex(1, 1)]
        @test strel(se) === se
        se_mask = @inferred strel(Bool, se)
        @test axes(se_mask) == (-1:1, -1:1)
        @test se_mask == centered(Bool[1 0 0; 0 1 0; 0 0 1])
        # internally, Bool is converted to SEMask
        se_mask = @inferred strel(ImageMorphology.SEMask{2}(), se)
        @test se_mask == centered(Bool[1 0 0; 0 1 0; 0 0 1])

        se = CartesianIndices((-1:1, -1:1))
        se_mask = @inferred strel(Bool, se)
        @test se_mask == centered(trues((3, 3)))

        se = CartesianIndices((1:2, -1:1))
        se_mask = @inferred strel(Bool, se)
        @test axes(se_mask) == (-2:2, -1:1)
        @test se_mask == centered(Bool[0 0 0; 0 0 0; 0 1 0; 1 1 1; 1 1 1])

        se = @inferred strel(Bool, CartesianIndex{2}[])
        @test se == centered(reshape(Bool[1], (1, 1))) # v1.7: Bool[1;;]

        # invalid cases
        se = [CartesianIndex(1), CartesianIndex(1, 1)]
        @test_throws ErrorException strel(se)
    end

    @testset "mask to offset" begin
        se = centered(Bool[1, 0, 1])
        se_offsets = @inferred strel(CartesianIndex, se)
        @test se_offsets == [CartesianIndex(-1), CartesianIndex(1)]
        @test se_offsets == strel(CartesianIndex, centered(Bool[1, 1, 1]))

        se = centered(Bool[1 0 0; 0 1 0; 0 0 1])
        se_offsets = @inferred strel(CartesianIndex, se)
        @test se_offsets == [CartesianIndex(-1, -1), CartesianIndex(1, 1)]
        # internally, CartesianIndex is converted to SEOffset
        se_offsets = @inferred strel(ImageMorphology.SEOffset{2}(), se)
        @test se_offsets == [CartesianIndex(-1, -1), CartesianIndex(1, 1)]

        se = centered(trues((3, 3)))
        se_offsets = @inferred strel(CartesianIndex, se)
        @test se_offsets == filter(x -> !iszero(x), vec(CartesianIndices((-1:1, -1:1))))

        se = @inferred strel(CartesianIndex, centered(falses(3, 3)))
        @test isempty(se)

        # test deprecation
        msg = @capture_err strel(CartesianIndex, trues(3, 3))
        @test isempty(msg) || contains(msg, "connectivity mask is expected to be a centered bool array")
        # but only for 1-based array
        se = OffsetArray(trues(3, 3), -1, -1)
        msg = "`connectivity` must be symmetric bool array"
        @test_throws ArgumentError(msg) strel(CartesianIndex, se)
    end

    @testset "center point" begin
        for se in centered.([Bool[1 0 0; 0 1 0; 0 0 1], Bool[1 0 0; 0 0 0; 0 0 1]])
            # center point is always excluded in offsets
            se_offsets = @inferred strel(CartesianIndex, se)
            @test se_offsets == [CartesianIndex(-1, -1), CartesianIndex(1, 1)]
            # but included in mask
            se_mask = @inferred strel(Bool, se_offsets)
            @test se_mask == centered(Bool[1 0 0; 0 1 0; 0 0 1])
        end
    end
end

@testset "strel_chain" begin
    se1 = centered(Bool[0 1 0; 0 1 0; 0 1 0])
    se2 = centered(Bool[0 0 0; 1 1 1; 0 0 0])
    se = @inferred strel_chain(se1, se2)
    @test se isa ImageMorphology.SEChainArray
    @test eltype(se) == Bool
    @test axes(se) == (-1:1, -1:1)
    @test se == strel_box((3, 3))
    @test se == strel_chain([se1, se2]) == strel_chain((se1, se2))
    @test se === strel_chain(se)
    @test strel_size(se) == (3, 3)

    for f in Any[dilate, erode]
        se1 = centered(rand(Bool, 3, 3))
        se2 = centered(rand(Bool, 3, 3))
        img = rand(64, 64)
        out_pipe = f(f(img, se1), se2)
        out_se = f(img, strel_chain(se1, se2))
        R = CartesianIndices(img)[2:(end - 1), 2:(end - 1)] # inner region
        @test out_pipe[R] == out_se[R]
    end
end

@testset "strel_product" begin
    se1 = centered(Bool[1, 1, 1])
    se2 = centered(Bool[1, 1, 1])
    se = strel_product(se1, se2)
    @test_broken @inferred strel_product(se1, se2) # if necessary, fix this inference issue
    @test se isa ImageMorphology.SEChainArray
    @test se == strel_box((3, 3))

    se_list = [centered(rand(Bool, sz...)) for sz in Any[(3,), (3, 3), (3, 3, 3)]]
    se = strel_product(se_list)
    @test ndims(se) == 6
end

@testset "strel_diamond" begin
    @testset "N=1" begin
        img = rand(5)
        se = @inferred strel_diamond(img)
        @test se isa ImageMorphology.SEDiamondArray
        @test eltype(se) == Bool
        @test axes(se) == (-1:1,)
        @test se == centered(Bool[1, 1, 1])

        se = @inferred strel_diamond((5,); r=1)
        @test se == centered(Bool[0, 1, 1, 1, 0])

        se = @inferred strel_diamond((5,))
        @test se == centered(Bool[1, 1, 1, 1, 1])
    end

    @testset "N=2" begin
        img = rand(5, 5)
        se = @inferred strel_diamond(img)
        @test se isa ImageMorphology.SEDiamondArray
        @test eltype(se) == Bool
        @test axes(se) == (-1:1, -1:1)
        @test se == centered(Bool[0 1 0; 1 1 1; 0 1 0])
        @test se == strel_diamond((3, 3), (1, 2); r=1)
        @test se == strel_diamond(img, (1, 2); r=1)
        @test se == strel_diamond(img, :) == strel_diamond(img, (:,))

        se = @inferred strel_diamond(img, (1,))
        @test se == @inferred strel_diamond(img, 1)
        @test se == centered(reshape(Bool[1, 1, 1], (3, 1)))

        se = @inferred strel_diamond((5, 5), (2, 1))
        @test se == strel_diamond((5, 5), 1:2) # not inferable

        se = @inferred strel_diamond((5, 5))
        @test se == centered(Bool[0 0 1 0 0; 0 1 1 1 0; 1 1 1 1 1; 0 1 1 1 0; 0 0 1 0 0])

        se = @inferred strel_diamond((3, 5))
        @test se == centered(Bool[0 1 1 1 0; 1 1 1 1 1; 0 1 1 1 0])

        se = @inferred strel_diamond((3, 5), (1,))
        @test se == @inferred strel_diamond((3, 5), 1)
        @test se == centered(Bool[0 0 1 0 0; 0 0 1 0 0; 0 0 1 0 0])

        se = @inferred strel_diamond((3, 5), (2,))
        @test se == centered(Bool[0 0 0 0 0; 1 1 1 1 1; 0 0 0 0 0])

        se = @inferred strel_diamond((3, 5); r=2)
        @test se == centered(Bool[0 1 1 1 0; 1 1 1 1 1; 0 1 1 1 0])
    end

    @testset "N=3" begin
        img = rand(5, 5, 5)
        se = @inferred strel_diamond(img)
        @test se isa ImageMorphology.SEDiamondArray
        @test eltype(se) == Bool
        @test axes(se) == (-1:1, -1:1, -1:1)
        @test se[:, :, -1] == se[:, :, 1] == centered(Bool[0 0 0; 0 1 0; 0 0 0])
        @test se[:, :, 0] == strel_diamond((3, 3))

        se = @inferred strel_diamond((3, 3, 3), (1, 2))
        @test se[:, :, -1] == se[:, :, 1] == centered(falses((3, 3)))
        @test se[:, :, 0] == strel_diamond((3, 3))
    end

    @testset "strel conversion" begin
        se = strel_diamond((3, 3))
        se_offsets = @inferred strel(CartesianIndex, se)
        @test se_offsets == [CartesianIndex(0, -1), CartesianIndex(-1, 0), CartesianIndex(1, 0), CartesianIndex(0, 1)]
        se_mask = @inferred strel(Bool, se)
        # not BitMatrix as SEDiamondArray provides more information of what the SE is
        @test se_mask isa ImageMorphology.SEDiamondArray
        @test se_mask === se
    end

    # edge cases
    img = rand(5, 5)
    err = ArgumentError("`axes` length should be at least 2")
    @test_throws err strel_diamond((3,), (1, 2))
    err = ArgumentError("size should be odd integers")
    @test_throws err strel_diamond((2, 3))
    err = ArgumentError("dims should be unique")
    @test_throws err strel_diamond((3, 3), (1, 1))
    err = ArgumentError("all `dims` values should be less than or equal to 2")
    @test_throws err strel_diamond((3, 3), (5,))
end

@testset "strel_box" begin
    @testset "N=1" begin
        img = rand(5)
        se = @inferred strel_box(img)
        @test se isa ImageMorphology.SEBoxArray
        @test eltype(se) == Bool
        @test se == centered(Bool[1, 1, 1])

        se = @inferred strel_box((5,))
        @test se == centered(Bool[1, 1, 1, 1, 1])

        se = @inferred strel_box((5,); r=1)
        @test se == centered(Bool[0, 1, 1, 1, 0])
    end

    @testset "N=2" begin
        img = rand(5, 5)
        se = @inferred strel_box(img)
        @test se isa ImageMorphology.SEBoxArray
        @test eltype(se) == Bool
        @test se == centered(Bool[1 1 1; 1 1 1; 1 1 1])
        @test se == strel_box(img; r=1)
        @test se == strel_box(img, :) == strel_box(img, (:,))

        se = @inferred strel_box(img, (1,))
        @test se == @inferred strel_box((3, 1), (1,))
        @test se == @inferred strel_box((3, 1), 1)

        se = @inferred strel_box(img; r=2)
        @test se == centered(trues((5, 5)))

        se = @inferred strel_box(img; r=(1, 2))
        @test se == centered(trues((3, 5)))

        se = @inferred strel_box((5, 5), (2, 1))
        @test se == @inferred strel_box((5, 5), 1:2)

        se = @inferred strel_box((3, 5))
        @test se == centered(Bool[1 1 1 1 1; 1 1 1 1 1; 1 1 1 1 1])

        se = @inferred strel_box((3, 5); r=1)
        @test se == centered(Bool[0 1 1 1 0; 0 1 1 1 0; 0 1 1 1 0])

        se = @inferred strel_box((3, 5); r=(1, 0))
        @test se == @inferred strel_box((3, 5), (1,))
        @test se == @inferred strel_box((3, 5), 1)
        @test se == centered(Bool[0 0 1 0 0; 0 0 1 0 0; 0 0 1 0 0])

        se = @inferred strel_box((3, 5); r=(0, 1))
        @test se == centered(Bool[0 0 0 0 0; 0 1 1 1 0; 0 0 0 0 0])

        se = @inferred strel_box((3, 5), (2,))
        @test se == centered(Bool[0 0 0 0 0; 1 1 1 1 1; 0 0 0 0 0])

        se = @inferred strel_box((3, 5); r=2)
        @test se == centered(Bool[1 1 1 1 1; 1 1 1 1 1; 1 1 1 1 1])
    end

    @testset "N=3" begin
        img = rand(5, 5, 5)
        se = @inferred strel_box(img)
        @test se isa ImageMorphology.SEBoxArray
        @test eltype(se) == Bool
        @test se[:, :, -1] == se[:, :, 1] == centered(trues((3, 3)))
        @test se[:, :, 0] == strel_box((3, 3))

        se = @inferred strel_box((3, 3, 3), (1, 2))
        @test axes(se) == (-1:1, -1:1, -1:1)
        @test se[:, :, 0] == strel_box((3, 3))
    end

    @testset "strel conversion" begin
        se = strel_box((3, 3))
        se_offsets = @inferred strel(CartesianIndex, se)
        @test se_offsets == filter(i -> !iszero(i), vec(CartesianIndices((-1:1, -1:1))))
        se_mask = @inferred strel(Bool, se)
        # not BitMatrix as SEBoxArray provides more information of what the SE is
        @test se_mask isa ImageMorphology.SEBoxArray
        @test se_mask === se
    end
end

@testset "strel_split" begin
    se = strel_diamond((3, 3))
    upper, lower = ImageMorphology.strel_split(se)
    @test upper == centered(Bool[0 1 0; 1 1 0; 0 0 0])
    @test lower == centered(Bool[0 0 0; 0 1 1; 0 1 0])
    se = strel(CartesianIndex, se)
    upper, lower = ImageMorphology.strel_split(se)
    @test upper == [CartesianIndex(0, -1), CartesianIndex(-1, 0)]
    @test lower == [CartesianIndex(1, 0), CartesianIndex(0, 1)]

    # non-symmetric is not allowed
    se = collect(strel_diamond((3, 3)))
    se[1, 1] = true
    msg = "structuring element must be symmetric with respect to its center"
    @test_throws ArgumentError(msg) ImageMorphology.strel_split(se)
end

@testset "strel_type" begin
    @test strel_type(strel_diamond((3, 3))) isa ImageMorphology.SEDiamond
    @test strel_type(strel_box((3, 3))) isa ImageMorphology.SEBox
    @test strel_type([CartesianIndex(1, 2)]) isa ImageMorphology.SEOffset
    @test strel_type(CartesianIndices((-1:1, -1:1))) isa ImageMorphology.SEOffset
    @test strel_type(trues(3, 3)) isa ImageMorphology.SEMask

    msg = "invalid structuring element data type: $(Vector{CartesianIndex})"
    @test_throws ErrorException(msg) strel_type(CartesianIndex[])
end

@testset "strel_size" begin
    se = strel_diamond((5, 5), (1,))
    @test strel_size(se) == strel_size(centered(collect(se))) == (5, 1)

    se = strel_box((5, 5), (2,))
    @test strel_size(se) == strel_size(centered(collect(se))) == (1, 5)

    se = [CartesianIndex(-2, -2), CartesianIndex(1, 1)]
    @test strel_size(se) == (5, 5)
end

@testset "centered" begin
    se = strel_diamond((3, 3))
    @test centered(se) === se
    se = strel_box((3, 3))
    @test centered(se) === se
end


text:
fname = path*"/SampleFiles/fdsn.conf"
hood_reg = Float64[44.8, 46.0, -122.4, -121.0]
rainier_rad = Float64[46.852886, -121.760374, 0.0, 0.1]

printstyled("  FDSNevq\n", color=:light_green)

# FDSNevq
printstyled("    single-server query\n", color=:light_green)
(H,R) = FDSNevq("2011-03-11T05:47:00", mag=[3.0, 9.9], nev=1, src="IRIS", v=0)
(H,R) = FDSNevq("201103110547", mag=[3.0, 9.9], nev=10, src="IRIS", v=0)
@test length(H)==9

printstyled("    single-server query without nev specified\n", color=:light_green)
(H,R) = FDSNevq("2018-06-01",reg=[32.0,38.0,-120.0,-115.0,-50.0,50.0],mag=[2.0,8.0],evw=[0.,375243600.0]);
@test length(H) == length(R)
@test length(H) > 1000

printstyled("    multi-server query\n", color=:light_green)
open("FDSNevq.log", "w") do out
  redirect_stdout(out) do
    ot = replace(split(string(now()),'.')[1], r"[-,:,A-Z,a-z]" => "")
    (H,R) = FDSNevq(ot, mag=[3.0, 9.9], evw=[-86400.0, 0.0], src="all", nev=10, v=2)
  end
end

printstyled("    radius search (rad=)\n", color=:light_green)
(H,R) = FDSNevq("20190101000000", rad=rainier_rad, evw=[31536000.0, 31536000.0], mag=[0.0, 2.9], nev=100, src="IRIS", v=0)

printstyled("    partly-specified region search (reg=)\n", color=:light_green)
(H,R) = FDSNevq("20120601000000", reg=hood_reg, evw=[31536000.0, 31536000.0], mag=[0.0, 2.9], nev=100, src="IRIS", v=0)

# FDSNevt
printstyled("  FDSNevt\n", color=:light_green)
S = FDSNevt("201103110547", "PB.B004..EH?,PB.B004..BS?,PB.B001..BS?,PB.B001..EH?", v=0)

ext: jl


text:
"""
Clears all waveforms and any data added by calls to the
[*wave_add_**] functions.

...
wave_clear(pi, )
...
"""
function wave_clear(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVCLR, 0, 0))
end

"""
Starts a new empty waveform.

You would not normally need to call this function as it is
automatically called after a waveform is created with the
[*wave_create*] function.

...
wave_add_new(pi, )
...
"""
function wave_add_new(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVNEW, 0, 0))
end

"""
Adds a list of pulses to the current waveform.

pulses:= list of pulses to add to the waveform.

Returns the new total number of pulses in the current waveform.

The pulses are interleaved in time order within the existing
waveform (if any).

Merging allows the waveform to be built in parts, that is the
settings for GPIO#1 can be added, and then GPIO#2 etc.

If the added waveform is intended to start after or within
the existing waveform then the first pulse should consist
solely of a delay.

...
G1=4
G2=24

set_mode(pi, G1, pigpio.OUTPUT)
set_mode(pi, G2, pigpio.OUTPUT)

flash_500=[] # flash every 500 ms
flash_100=[] # flash every 100 ms

#                              ON     OFF  DELAY

flash_500.append(pigpio.pulse(1<<G1, 1<<G2, 500000))
flash_500.append(pigpio.pulse(1<<G2, 1<<G1, 500000))

flash_100.append(pigpio.pulse(1<<G1, 1<<G2, 100000))
flash_100.append(pigpio.pulse(1<<G2, 1<<G1, 100000))

wave_clear(pi, ) # clear any existing waveforms

wave_add_generic(pi, flash_500) # 500 ms flashes
f500 = wave_create(pi, ) # create and save id

wave_add_generic(pi, flash_100) # 100 ms flashes
f100 = wave_create(pi, ) # create and save id

wave_send_repeat(pi, f500)

time.sleep(4)

wave_send_repeat(pi, f100)

time.sleep(4)

wave_send_repeat(pi, f500)

time.sleep(4)

wave_tx_stop(pi, ) # stop waveform

wave_clear(pi, ) # clear all waveforms
...
"""
function wave_add_generic(self::Pi, pulses)
    # pigpio message format

    # I p1 0
    # I p2 0
    # I p3 pulses * 12
    ## extension ##
    # III on/off/delay * pulses
    if length(pulses)
        ext = bytearray()
        for p in pulses
        ext.extend(pack("III", p.gpio_on, p.gpio_off, p.delay))
        end
        extents = [ext]
        return _u2i(_pigpio_command_ext(
            self.sl, _PI_CMD_WVAG, 0, 0, length(pulses)*12, extents))
    else
        return 0
    end
end

"""
Adds a waveform representing serial data to the existing
waveform (if any).  The serial data starts [*offset*]
microseconds from the start of the waveform.

user_gpio:= GPIO to transmit data.  You must set the GPIO mode
   to output.
baud:= 50-1000000 bits per second.
data:= the bytes to write.
offset:= number of microseconds from the start of the
   waveform, default 0.
bb_bits:= number of data bits, default 8.
bb_stop:= number of stop half bits, default 2.

Returns the new total number of pulses in the current waveform.

The serial data is formatted as one start bit, [*bb_bits*]
data bits, and [*bb_stop*]/2 stop bits.

It is legal to add serial data streams with different baud
rates to the same waveform.

The bytes required for each character depend upon [*bb_bits*].

For [*bb_bits*] 1-8 there will be one byte per character.
For [*bb_bits*] 9-16 there will be two bytes per character.
For [*bb_bits*] 17-32 there will be four bytes per character.

...
wave_add_serial(pi, 4, 300, 'Hello world')

wave_add_serial(pi, 4, 300, b"Hello world")

wave_add_serial(pi, 4, 300, b'\\x23\\x01\\x00\\x45')

wave_add_serial(pi, 17, 38400, [23, 128, 234], 5000)
...
"""
function wave_add_serial(
    self::Pi, user_gpio, baud, data, offset=0, bb_bits=8, bb_stop=2)
    if length(data)
        extents = [pack("III", bb_bits, bb_stop, offset), data]
        return _u2i(_pigpio_command_ext(
            self.sl, _PI_CMD_WVAS, user_gpio, baud, length(data)+12, extents))
    else
        return 0
    end
end

"""
Creates a waveform from the data provided by the prior calls
to the [*wave_add_**] functions.

Returns a wave id (>=0) if OK,  otherwise PI_EMPTY_WAVEFORM,
PI_TOO_MANY_CBS, PI_TOO_MANY_OOL, or PI_NO_WAVEFORM_ID.

The data provided by the [*wave_add_**] functions is consumed by
this function.

As many waveforms may be created as there is space available.
The wave id is passed to [*wave_send_**] to specify the waveform
to transmit.

Normal usage would be

Step 1. [*wave_clear*] to clear all waveforms and added data.

Step 2. [*wave_add_**] calls to supply the waveform data.

Step 3. [*wave_create*] to create the waveform and get a unique id

Repeat steps 2 and 3 as needed.

Step 4. [*wave_send_**] with the id of the waveform to transmit.

A waveform comprises one or more pulses.

A pulse specifies

1) the GPIO to be switched on at the start of the pulse.
2) the GPIO to be switched off at the start of the pulse.
3) the delay in microseconds before the next pulse.

Any or all the fields can be zero.  It doesn't make any sense
to set all the fields to zero (the pulse will be ignored).

When a waveform is started each pulse is executed in order with
the specified delay between the pulse and the next.

...
wid = wave_create(pi, )
...
"""
function wave_create(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVCRE, 0, 0))
end

"""
This function deletes the waveform with id wave_id.

wave_id:= >=0 (as returned by a prior call to [*wave_create*]).

Wave ids are allocated in order, 0, 1, 2, etc.

...
wave_delete(pi, 6) # delete waveform with id 6

wave_delete(pi, 0) # delete waveform with id 0
...
"""
function wave_delete(self::Pi, wave_id)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVDEL, wave_id, 0))
end

"""
This function is deprecated and has been removed.

Use [*wave_create*]/[*wave_send_**] instead.
"""
function wave_tx_start(self::Pi) # DEPRECATED
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVGO, 0, 0))
end

"""
This function is deprecated and has beeen removed.

Use [*wave_create*]/[*wave_send_**] instead.
"""
function wave_tx_repeat(self::Pi) # DEPRECATED
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVGOR, 0, 0))
end

"""
Transmits the waveform with id wave_id.  The waveform is sent
once.

NOTE: Any hardware PWM started by [*hardware_PWM*] will
be cancelled.

wave_id:= >=0 (as returned by a prior call to [*wave_create*]).

Returns the number of DMA control blocks used in the waveform.

...
cbs = wave_send_once(pi, wid)
...
"""
function wave_send_once(self::Pi, wave_id)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVTX, wave_id, 0))
end

"""
Transmits the waveform with id wave_id.  The waveform repeats
until wave_tx_stop is called or another call to [*wave_send_**]
is made.

NOTE: Any hardware PWM started by [*hardware_PWM*] will
be cancelled.

wave_id:= >=0 (as returned by a prior call to [*wave_create*]).

Returns the number of DMA control blocks used in the waveform.

...
cbs = wave_send_repeat(pi, wid)
...
"""
function wave_send_repeat(self::Pi, wave_id)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVTXR, wave_id, 0))
end

"""
Transmits the waveform with id wave_id using mode mode.

wave_id:= >=0 (as returned by a prior call to [*wave_create*]).
mode:= WAVE_MODE_ONE_SHOT, WAVE_MODE_REPEAT,
    WAVE_MODE_ONE_SHOT_SYNC, or WAVE_MODE_REPEAT_SYNC.

WAVE_MODE_ONE_SHOT: same as [*wave_send_once*].

WAVE_MODE_REPEAT same as [*wave_send_repeat*].

WAVE_MODE_ONE_SHOT_SYNC same as [*wave_send_once*] but tries
to sync with the previous waveform.

WAVE_MODE_REPEAT_SYNC same as [*wave_send_repeat*] but tries
to sync with the previous waveform.

WARNING: bad things may happen if you delete the previous
waveform before it has been synced to the new waveform.

NOTE: Any hardware PWM started by [*hardware_PWM*] will
be cancelled.

wave_id:= >=0 (as returned by a prior call to [*wave_create*]).

Returns the number of DMA control blocks used in the waveform.

...
cbs = wave_send_using_mode(pi, wid, WAVE_MODE_REPEAT_SYNC)
...
"""
function wave_send_using_mode(self::Pi, wave_id, mode)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVTXM, wave_id, mode))
end

"""
Returns the id of the waveform currently being
transmitted.

Returns the waveform id or one of the following special
values

WAVE_NOT_FOUND (9998) - transmitted wave not found.
NO_TX_WAVE (9999) - no wave being transmitted.

...
wid = wave_tx_at(pi, )
...
"""
function wave_tx_at(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVTAT, 0, 0))
end

"""
Returns 1 if a waveform is currently being transmitted,
otherwise 0.

...
wave_send_once(pi, 0) # send first waveform

while wave_tx_busy(pi, ): # wait for waveform to be sent
time.sleep(0.1)

wave_send_once(pi, 1) # send next waveform
...
"""
function wave_tx_busy(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVBSY, 0, 0))
end

"""
Stops the transmission of the current waveform.

This function is intended to stop a waveform started with
wave_send_repeat.

...
wave_send_repeat(pi, 3)

time.sleep(5)

wave_tx_stop(pi, )
...
"""
function wave_tx_stop(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVHLT, 0, 0))
end

"""
This function transmits a chain of waveforms.

NOTE: Any hardware PWM started by [*hardware_PWM*]
will be cancelled.

The waves to be transmitted are specified by the contents
of data which contains an ordered list of [*wave_id*]s
and optional command codes and related data.

Returns 0 if OK, otherwise PI_CHAIN_NESTING,
PI_CHAIN_LOOP_CNT, PI_BAD_CHAIN_LOOP, PI_BAD_CHAIN_CMD,
PI_CHAIN_COUNTER, PI_BAD_CHAIN_DELAY, PI_CHAIN_TOO_BIG,
or PI_BAD_WAVE_ID.

Each wave is transmitted in the order specified.  A wave
may occur multiple times per chain.

A blocks of waves may be transmitted multiple times by
using the loop commands. The block is bracketed by loop
start and end commands.  Loops may be nested.

Delays between waves may be added with the delay command.

The following command codes are supported

Name         @ Cmd & Data @ Meaning
Loop Start   @ 255 0      @ Identify start of a wave block
Loop Repeat  @ 255 1 x y  @ loop x + y*256 times
Delay        @ 255 2 x y  @ delay x + y*256 microseconds
Loop Forever @ 255 3      @ loop forever

If present Loop Forever must be the last entry in the chain.

The code is currently dimensioned to support a chain with
roughly 600 entries and 20 loop counters.

...
#!/usr/bin/env python

import time
import pigpio

WAVES=5
GPIO=4

wid=[0]*WAVES

pi = pigpio.pi() # Connect to local Pi.

set_mode(pi, GPIO, pigpio.OUTPUT);

for i in range(WAVES)
pi.wave_add_generic([
pigpio.pulse(1<<GPIO, 0, 20),
pigpio.pulse(0, 1<<GPIO, (i+1)*200)]);

wid[i] = wave_create(pi, );

pi.wave_chain([
wid[4], wid[3], wid[2],       # transmit waves 4+3+2
255, 0,                       # loop start
wid[0], wid[0], wid[0],    # transmit waves 0+0+0
255, 0,                    # loop start
   wid[0], wid[1],         # transmit waves 0+1
   255, 2, 0x88, 0x13,     # delay 5000us
255, 1, 30, 0,             # loop end (repeat 30 times)
255, 0,                    # loop start
   wid[2], wid[3], wid[0], # transmit waves 2+3+0
   wid[3], wid[1], wid[2], # transmit waves 3+1+2
255, 1, 10, 0,             # loop end (repeat 10 times)
255, 1, 5, 0,                 # loop end (repeat 5 times)
wid[4], wid[4], wid[4],       # transmit waves 4+4+4
255, 2, 0x20, 0x4E,           # delay 20000us
wid[0], wid[0], wid[0],       # transmit waves 0+0+0
])

while wave_tx_busy(pi, )
time.sleep(0.1);

for i in range(WAVES)
wave_delete(pi, wid[i])

stop(pi, )
...
"""
function wave_chain(self::Pi, data)
# I p1 0
# I p2 0
# I p3 len
## extension ##
# s len data bytes
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_WVCHA, 0, 0, length(data), [data]))
end

"""
Returns the length in microseconds of the current waveform.

...
micros = wave_get_micros(pi, )
...
"""
function wave_get_micros(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVSM, 0, 0))
end

"""
Returns the maximum possible size of a waveform in microseconds.

...
micros = wave_get_max_micros(pi, )
...
"""
function wave_get_max_micros(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVSM, 2, 0))
end

"""
Returns the length in pulses of the current waveform.

...
pulses = wave_get_pulses(pi, )
...
"""
function wave_get_pulses(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVSP, 0, 0))
end

"""
Returns the maximum possible size of a waveform in pulses.

...
pulses = wave_get_max_pulses(pi, )
...
"""
function wave_get_max_pulses(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVSP, 2, 0))
end

"""
Returns the length in DMA control blocks of the current
waveform.

...
cbs = wave_get_cbs(pi, )
...
"""
function wave_get_cbs(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVSC, 0, 0))
end

"""
Returns the maximum possible size of a waveform in DMA
control blocks.

...
cbs = wave_get_max_cbs(pi, )
...
"""
function wave_get_max_cbs(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WVSC, 2, 0))
end


text:
using JuLIP
using Base.Test
using JuLIP.Testing

verbose=true

# check whether on CI
isCI = haskey(ENV, "CI")
notCI = !isCI
eam_W4 = nothing

# check whether ASE is available
hasase = true
try
   import ASE
catch
   hasase = false
end

julip_tests = [
   ("testaux.jl", "Miscellaneous"),
   ("test_atoms.jl", "Atoms"),
   ("test_build.jl", "Build"),
   ("testanalyticpotential.jl", "Analytic Potential"),
   ("testpotentials.jl", "Potentials"),
   # ("test_ad.jl", "AD Potentials"),
   ("testvarcell.jl", "Variable Cell"),
   ("testhessian.jl", "Hessian"),
   ("testsolve.jl", "Solve"),
]

# remove testsolve if on Travis
if isCI
   julip_tests = julip_tests[1:end-1]
end

# "testexpvarcell.jl";  # USE THIS TO WORK ON EXPCELL IMPLEMENTATION

# ===== some prototype potentials ======
print("Loading some interatomic potentials . .")
data = joinpath(dirname(@__FILE__), "..", "data") * "/"
eam_Fe = JuLIP.Potentials.EAM(data * "pfe.plt", data * "ffe.plt", data * "F_fe.plt")
print(" .")
eam_W = JuLIP.Potentials.FinnisSinclair(data*"W-pair-Wang-2014.plt", data*"W-e-dens-Wang-2014.plt")
print(" .")
try
   eam_W4 = JuLIP.Potentials.EAM(data * "w_eam4.fs")
catch
   eam_W4 = nothing
end
println(" done.")

println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")
println("  Starting JuLIP Tests")
println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")

@testset "JuLIP" begin
   for (testfile, testid) in julip_tests
      println("=======================")
      println("Testset $(testid)")
      println("=======================")
      @testset "$(testid)" begin include(testfile); end
   end
end


text:
#   Unit "filters" of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019-2022,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements filtering procedures interfacing the DSP.jl package

FilterDesign = Union{ZeroPoleGain, FIRWindow}

#######################################################################
"""
```julia
function filterbank(x         :: Vector{T},
                    sr        :: Int,
                    bandwidth :: IntOrReal    = 2;
                filtkind      :: FilterDesign = Butterworth(2),
                fmin          :: IntOrReal    = bandwidth,
                fmax          :: IntOrReal    = sr÷2,
                ⏩           :: Bool         = true) where T<:Real
```

Pass signal vector `x` throught a bank of band-pass filters,
given sampling rate `sr` and `bandwidth` parameter.

The kind of filters is specified by optional keyword argument `filtkind`,
of type [FilterDesign](@ref), using the DSP package.
By default `filtkind` is a forward-backward (linear phase response)
Butterworth IIR filter of order 2 in each direction
(hence, 4th order total). See [notes on DSP package useful functions](@ref)
for tips on how to design other IIR and FIR filters.

Return 2-tuple `(f, Y)`, where `f` is a vector holding the center
frequencies of the filter bank band-pass regions
and `Y` a matrix holding in the ``i^{th}`` column the signal `x`
band-pass iltered by the ``i^{th}`` band-pass filter.
Hence, ```size(Y, 1)=length(x)``` and ```size(Y, 2)=length(f)```.

The filter bank is designed by means of argument `bandwidth`
and optional keyword arguments `fmin` and `fmax`.
These three arguments can be passed either as integers or real numbers.
All band-pass regions have bandwidth equal to the `bandwidth` argument
and overlap with adjacent band-pass regions.
By default the lower limit of the first band-pass region is set at
`bandwidth` Hz, and successive band-pass regions are defined up to,
but excluding, the Nyquist frequency (``sr/2``).
If `fmin` is specified (in Hz), the center frequency of the first band-pass
region is set as close as possible, but not below, `fmin`.
If `fmax` is specified (in Hz), the center frequency of the last band-pass
region is set as close as possible, but not above, `fmax`.

Here are some examples of filter bank definitions given different arguments
(`sr`=16 in all examples).

| bandwidth | fmin | fmax | center frequencies |       band-pass regions       |
|:---------:|:----:|:----:|:------------------:|:-----------------------------:|
|     4     |   -  |  -   |         4          |              2-6              |
|     2     |   -  |  -   | 2, 3, 4, 5, 6      |    1-3, 2-4, 3-5, 4-6, 5-7    |
|     2     |   3  |  7   | 3, 4, 5, 6         |    2-4, 3-5, 4-6, 5-7         |
|     1     |   3  |  5   | 3, 3.5, 4, 4.5, 5  | 2.5-3.5, 3-4, 3.5-4.5, 4-5, 4.5-5.5 |
|    1.1    |   3  |  5   | 2.75, 3.3, 3.85, 4.4, 4.95 | 2.2-3.3, 2.75-8.85, 3.3-4.4, 3.85-4.95, 4.4-5.5 |
|    1.9    |   3  |  5   | 2.85, 3.8, 4.75    | 1.9-3.8, 2.85-4.75, 3.8-5.7 |

!!! note "Nota Bene"
    At least `sr` samples should be trimmed at the beginning and end of the
    output signal `Y`, as those samples are severely distorted by the filtering
    process.

    If keyword optional argument ⏩ is true (default), the filtering is
    multi-threaded across band-pass filters. See [Threads](@ref).

This function is called by the following functions operating
on time-frequency reprsentations: [`TFanalyticsignal`](@ref),
[`TFamplitude`](@ref), [`TFphase`](@ref), [`meanAmplitude`](@ref),
[`concentration`](@ref), [`meanDirection`](@ref), [`comodulation`](@ref),
[`coherence`](@ref).

**Examples**:
```julia
using FourierAnalysis, DSP, Plots
# generate a sinusoidal + noise
f, sr, t = 8, 128, 512
v=sinusoidal(1., f, sr, t, 0)
x=v+randn(t)
flabels, Y=filterbank(x, 128)
flabels, Y=filterbank(x, 128; fmin=4, fmax=32)
flabels, Y=filterbank(x, 128, 4; fmin=4, fmax=32)
flabels, Y=filterbank(x, 128, 4;
                      filtkind=Chebyshev2(8, 10),
                      fmin=4,
                      fmax=16)
# trick for plotting the signal filtered in the band-pass regions
for i=1:size(Y, 2) Y[:, i].+=convert(eltype(Y), i)*1.5 end
mylabels=Array{String}(undef, 1, length(flabels))
for i=1:length(flabels) mylabels[1, i]=string(flabels[i])*" Hz" end
plot(Y; c=:grey, labels=mylabels)
```
"""
function filterbank(x         :: Vector{T},
                    sr        :: Int,
                    bandwidth :: IntOrReal    = 2;
                filtkind      :: FilterDesign = Butterworth(2),
                fmin          :: IntOrReal    = bandwidth,
                fmax          :: IntOrReal    = sr÷2,
                ⏩           :: Bool         = true) where T<:Real

    df, dm, ff, Bp=digitalfilter, filtkind, filtfilt, Bandpass # aliases for DSP.jl stuff
    fb=_ffilterbank(sr, bandwidth; fmin=fmin, fmax=fmax) # get limits for band-passes
    nf=length(fb)-2 # eliminate last frequency, which may be the Nyquist frequency
    flt=[df(Bp(fb[f], fb[f+2]; fs=sr), dm) for f=1:nf] # f band-pass filters
    x_=[zeros(T, sr); x; zeros(T, sr)]     # Before filtering, `x` is padded with `sr` zeros at both sides in order to symmetrize the distortion of DSP.jl butterworth filters.
    Y=Matrix{T}(undef, length(x_), nf)
    _thread(⏩, nf) ? (@threads for f=1:nf Y[:, f]=ff(flt[f], x_) end) :
                               (for f=1:nf Y[:, f]=ff(flt[f], x_) end)
    return fb[2:end-1], Y[sr+1:end-sr, :] # trim the padded segments in Y
end

##################################################################
# internal function: Create a vector holding necessary information
# to build a filter bank, given sampling rate `rt`, the `bandwidth`
# and optional keyword parameters `fmin` and `fmax`.
# `bandwidth`, `fmax` and `fmin` can be either Integers or Floats.
# Let y be the output vector, a filter bank with length(y)-2 filters
# can be constructed with band-pass regions (in Hz) going from y[i-1] to y[i+1],
# where y[i+1]-y[i-1] = `bandwidth` for all i in [2:end-1].
# All the length(y)-2 filters are centered at y[i].
# By default the filter bank is defined up to, but exclude, the Nyquist frequency.
# If `fmin` is specified (in Hz), the lowest center frequency is >= `fmin`.
# If `fmax` is specified (in Hz), the highest center frequency is <= `fmax`.
#
# Examples:
#
# sr=16
# bw=4
# _ffilterbank(sr, bw) # [2.0, 4.0, 6.0]
#
# bw=2
# _ffilterbank(sr, bw) # [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
# _ffilterbank(sr, bw; fmin=3, fmax=7) # [2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
#
# bw=1
# _ffilterbank(sr, bw; fmin=3, fmax=5) # [2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5]
#
# bw=1.1
# _ffilterbank(sr, bw; fmin=3, fmax=5) # [2.2, 2.75, 3.3, 3.85, 4.4, 4.95, 5.5]
#
# bw=1.9
# _ffilterbank(sr, bw; fmin=3, fmax=5) # [1.9, 2.85, 3.8, 4.75, 5.7]
function _ffilterbank(sr        :: Int,
                      bandwidth :: IntOrReal;
                  fmin          :: IntOrReal = bandwidth,
                  fmax          :: IntOrReal = sr÷2)
    f=collect(range(0, step=bandwidth/2, stop=sr÷2))
    popfirst!(f)
    a=max(1, (findmin([abs(g-fmin) for g∈f])[2])-1)
    b=min(length(f), (findmin([abs(g-fmax)  for g∈f])[2])+1)
    f[b]>=sr÷2 ? (return f[a:b-1]) : (return f[a:b])
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 0b30d8f2-3788-11ec-0155-cd3fa7eec623
begin
	import Pkg
	Pkg.add(url = "https://github.com/frankhuettner/RICoverAlg.jl")
	using RICoverAlg
end

# ╔═╡ 2e419b18-204f-4604-ac62-006112717796
begin
	Pkg.add("Plots")
	using Plots
end

# ╔═╡ 5f1dbaf1-c88e-43d8-8039-2c5eeba35f23
md"To use the package, you have to add it from its github url:"

# ╔═╡ 8432015b-9df7-43cf-afb5-86e4ffac350e
UtilTable = [0   1   0   1;
    		 0   0   1   1;
    		1/2 1/2 1/2 1/2]

# ╔═╡ 69ad0d60-17e7-4c17-9f95-91bacf47cc80
λ = 0.4

# ╔═╡ edbb7ffd-7325-4389-9164-4eebcb446160
md"We will compute the result for different prior beliefs, g = [1+ρ; 1-ρ; 1-ρ; 1+ρ], where ρ iterates in 0.1 steps from -1 to 1."

# ╔═╡ 5ad0d918-9639-42c5-9f21-acbf4bb477cf
ρs = -1:0.1:1

# ╔═╡ b48b5d06-3b0e-47ea-99a6-671fab98743b
results = [risolve(UtilTable, [1+ρ; 1-ρ; 1-ρ; 1+ρ]./4, λ)[1][1] for ρ in ρs]

# ╔═╡ 477f0966-d094-483d-b9ef-14000282f6c0
begin
	plot(-1:.1:1, results, label="λ=0.4", c=:purple, marker=:square, alpha=0.6,
	    xlabel="ρ", ylabel="Unconditional Choice Probability")
	plot!(-1:.1:1, 0.5:-0.0125:0.25, label="λ=0", c=:black, marker=:o, alpha=0.6)
end

# ╔═╡ ba95ef6b-7aae-4c59-88d1-671400d4d7d3
savefig("img/MM_RBT.png")

# ╔═╡ Cell order:
# ╟─5f1dbaf1-c88e-43d8-8039-2c5eeba35f23
# ╠═0b30d8f2-3788-11ec-0155-cd3fa7eec623
# ╠═8432015b-9df7-43cf-afb5-86e4ffac350e
# ╠═69ad0d60-17e7-4c17-9f95-91bacf47cc80
# ╟─edbb7ffd-7325-4389-9164-4eebcb446160
# ╠═5ad0d918-9639-42c5-9f21-acbf4bb477cf
# ╠═b48b5d06-3b0e-47ea-99a6-671fab98743b
# ╠═2e419b18-204f-4604-ac62-006112717796
# ╠═477f0966-d094-483d-b9ef-14000282f6c0
# ╠═ba95ef6b-7aae-4c59-88d1-671400d4d7d3


text:
@testset "DataFrames" begin
    df = DataFrame(x1 = rand(5), x2 = rand(5), y = [:a,:a,:b,:a,:b])

    @test_throws ArgumentError targets(df)

    @test targets(:y, df) == [:a,:a,:b,:a,:b]
    @test targets(row->row[1,:y], df) == [:a,:a,:b,:a,:b]
    @test typeof(targets(:y, df)) == Vector{Symbol}
    @test eltype(targets(x->x, df)) <: SubDataFrame

    @test nobs(@inferred(undersample(:y, df))) === 4

    @test @inferred(getobs(df, 2)) == df[2,:]
    @test @inferred(getobs(datasubset(df, 2))) == df[2,:]

    @test typeof(datasubset(df, 2)) <: SubDataFrame
    @test @inferred(datasubset(df, [1,2,3,5])) == view(df, [1,2,3,5])
    @test @inferred(datasubset(df, 2)) == view(df, 2)
    @test @inferred(datasubset(datasubset(df, 2:3), 2)) == view(df, 3)
    @test @inferred(datasubset(df)) == view(df, 1:5)
end


text:
# SAM Header
# ==========

immutable SAMHeader
    metainfo::Vector{SAMMetaInfo}
end

function SAMHeader()
    return SAMHeader(SAMMetaInfo[])
end

function Base.eltype(::Type{SAMHeader})
    return SAMMetaInfo
end

function Base.length(header::SAMHeader)
    return length(header.metainfo)
end

function Base.start(header::SAMHeader)
    return 1
end

function Base.done(header::SAMHeader, i)
    return i > length(header.metainfo)
end

function Base.next(header::SAMHeader, i)
    return header.metainfo[i], i + 1
end

function Base.find(header::SAMHeader, key::AbstractString)
    return filter(m -> isequalkey(m, key), header.metainfo)
end

function Base.unshift!(header::SAMHeader, metainfo::SAMMetaInfo)
    unshift!(header.metainfo, metainfo)
    return header
end

function Base.push!(header::SAMHeader, metainfo::SAMMetaInfo)
    push!(header.metainfo, metainfo)
    return header
end


text:
##### Beginning of file

Test.@test( Base.VERSION >= VersionNumber("1.0") )

Test.@test( Snapshots.version() > VersionNumber(0) )

Test.@test(
    Snapshots.version() ==
        Snapshots.version(Snapshots)
    )

Test.@test(
    Snapshots.version() ==
        Snapshots.version(first(methods(Snapshots.eval)))
    )

Test.@test(
    Snapshots.version() ==
        Snapshots.version(Snapshots.eval)
    )

Test.@test(
    Snapshots.version() ==
        Snapshots.version(Snapshots.eval, (Any,))
    )

Test.@test( Snapshots.version(TestModuleA) == VersionNumber("1.2.3") )

Test.@test( Snapshots.version(TestModuleB) == VersionNumber("4.5.6") )

Test.@test_throws(
    ErrorException,
    Snapshots._TomlFile(joinpath(mktempdir(),"1","2","3","4")),
    )

##### End of file


text:
module TestGraph

using Test
import TreeParzen: Delayed, Graph, HP

# Delayed objects can be nested
a = HP.Normal(:a, 10.0, 2.0)
b = HP.Uniform(:b, a, 7.0)
ab = Dict(:x => [a, b])
ba = Dict(:x => [b, a])

@test typeof.(Graph.dfs(ab)) == [
    Delayed.Normal,
    HP.Normal,
    Delayed.Normal,
    HP.Normal,
    Delayed.Uniform,
    HP.Uniform,
]
@test typeof.(Graph.dfs(ba)) == [
    Delayed.Normal,
    HP.Normal,
    Delayed.Uniform,
    HP.Uniform,
    Delayed.Normal,
    HP.Normal,
]

# Using the same label more than once is not allowed
space = Dict(
    :x => HP.Normal(:a, 1.0, 2.0),
    :y => HP.Normal(:a, 1.0, 2.0),
)
@test_throws DomainError Graph.checkspace(space)

end
true


text: include("basics.jl")


text:
module UsingOverride

    export @override

    macro override(m)
        esc(_override(m))
    end

    function _override(m)
        me = eval(Main,m)
        ns = names(me)
        ex = Expr[]
        for n in ns
            n == m && continue
            ov = Expr(:(=),n, Expr(:call,:getfield,m, QuoteNode(n) ) )
            check = Expr(:call,:isdefined,:Main,QuoteNode(n))
            push!(ex,:( $check && $ov ))
        end
        use_ex = :( using $m )
        Expr(:block,use_ex,ex...)
    end

end


text:
module TestGetters
using Test
using PhaseSpaceIO
using PhaseSpaceIO.Testing
using PhaseSpaceIO.Getters
using Accessors
using ConstructionBase

@testset "test getters" begin
    p = arbitrary(IAEAParticle{0,1})
    p = @set p.x = 100
    @test x(p) ≈ 100
    p = @set p.E = 10
    p = @set p.weight = 1
    @test energy(p) == 10
    p = @set p.weight = 0.1
    @test energy(p) ≈ 1

    p = arbitrary(EGSParticle{Float32})
    p = @set p.x = 1000f0
    @test x(p) ≈ 1000
    p = @set p.zlast = nothing
    @test p.zlast == nothing
    p = @set p.zlast = 1f0
    @test p.zlast == 1f0
    l = p.latch
    @test getproperties(l) isa NamedTuple
    @test propertynames(l) == propertynames(getproperties(l))
    @test l == setproperties(l, getproperties(l))

    for p in [arbitrary(IAEAParticle{0,1}), arbitrary(EGSParticle{Float32})]
        nt = getproperties(p)
        @test getproperties(p) isa NamedTuple
        @test propertynames(nt) == propertynames(p)
        @test setproperties(p, nt) == p
    end
end

end


text:
A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]

function cost(X::ProductRepr)
    return cost([submanifold_components(X)...])
end
function cost(X::Array{Matrix{Float64},1})
    return -0.5 * norm(transpose(X[1]) * A * X[2])^2
end

function egrad(X::Array{Matrix{Float64},1})
    U = X[1]
    V = X[2]
    AV = A * V
    AtU = transpose(A) * U
    return [-AV * (transpose(AV) * U), -AtU * (transpose(AtU) * V)]
end

function rgrad(M::ProductManifold, X::ProductRepr)
    eG = egrad([submanifold_components(M, X)...])
    x = [submanifold_components(M, X)...]
    return Manifolds.ProductRepr(project.(M.manifolds, x, eG)...)
end

function e2rHess(
    M::Grassmann, x, ξ, eGrad::Matrix{T}, Hess::Matrix{T}
) where {T<:Union{U,Complex{U}}} where {U<:AbstractFloat}
    pxHess = project(M, x, Hess)
    xtGrad = x' * eGrad
    ξxtGrad = ξ * xtGrad
    return pxHess - ξxtGrad
end

function eHess(X::Array{Matrix{Float64},1}, H::Array{Matrix{Float64},1})
    U = X[1]
    V = X[2]
    Udot = H[1]
    Vdot = H[2]
    AV = A * V
    AtU = transpose(A) * U
    AVdot = A * Vdot
    AtUdot = transpose(A) * Udot
    return [
        -(
            AVdot * transpose(AV) * U +
            AV * transpose(AVdot) * U +
            AV * transpose(AV) * Udot
        ),
        -(
            AtUdot * transpose(AtU) * V +
            AtU * transpose(AtUdot) * V +
            AtU * transpose(AtU) * Vdot
        ),
    ]
end

function rhess(M::ProductManifold, X::ProductRepr, H::ProductRepr)
    x = [submanifold_components(M, X)...]
    h = [submanifold_components(M, H)...]
    eG = egrad(x)
    eH = eHess(x, h)
    return Manifolds.ProductRepr(e2rHess.(M.manifolds, x, h, eG, eH)...)
end

@testset "Manopt Trust-Region" begin
    import Random: seed!
    seed!(42)

    N = Grassmann(3, 2)
    M = N × N

    x = random_point(M)

    @test_throws ErrorException trust_regions(
        M, cost, (X) -> rgrad(M, X), x, rhess; ρ_prime=0.3
    )
    @test_throws ErrorException trust_regions(
        M, cost, (X) -> rgrad(M, X), x, rhess; Δ_bar=-0.1
    )
    @test_throws ErrorException trust_regions(M, cost, (X) -> rgrad(M, X), x, rhess; Δ=-0.1)
    @test_throws ErrorException trust_regions(
        M, cost, (X) -> rgrad(M, X), x, rhess; Δ_bar=0.1, Δ=0.11
    )

    X = trust_regions(M, cost, (X) -> rgrad(M, X), x, rhess; Δ_bar=4 * sqrt(2 * 2))
    opt = trust_regions(
        M, cost, (X) -> rgrad(M, X), x, rhess; Δ_bar=4 * sqrt(2 * 2), return_options=true
    )
    @test isapprox(M, X, get_solver_result(opt))

    @test cost(X) + 142.5 ≈ 0 atol = 10.0^(-13)

    XuR = trust_regions(
        M, cost, (X) -> rgrad(M, X), x, rhess; Δ_bar=4 * sqrt(2 * 2), useRandom=true
    )

    @test cost(XuR) + 142.5 ≈ 0 atol = 10.0^(-12)

    XaH = trust_regions(
        M,
        cost,
        (X) -> rgrad(M, X),
        x,
        (p, x, ξ) -> approxHessianFD(
            p,
            x,
            (X) -> rgrad(M, X),
            ξ;
            stepsize=2^(-9),
            vector_transport_method=ProductVectorTransport(
                ProjectionTransport(), ProjectionTransport()
            ),
        );
        stopping_criterion=StopWhenAny(
            StopAfterIteration(2000), StopWhenGradientNormLess(10^(-6))
        ),
        Δ_bar=4 * sqrt(2 * 2),
    )
    @test cost(XaH) + 142.5 ≈ 0 atol = 10.0^(-10)

    ξ = random_tangent(M, x)
    @test_throws ErrorException getHessian(
        SubGradientProblem(M, cost, (X) -> rgrad(M, X)), x, ξ
    )

    # Test the random step trust region
    p = HessianProblem(M, cost, (X) -> rgrad(M, X), rhess, (M, x, ξ) -> ξ)
    o = TrustRegionsOptions(
        x,
        StopAfterIteration(2000),
        10.0^(-8),
        sqrt(manifold_dimension(M)),
        ExponentialRetraction(),
        true,
        0.1,
        1000.0,
    )
    @test step_solver!(p, o, 0) == nothing

    η = truncated_conjugate_gradient_descent(M, cost, (X) -> rgrad(M, X), x, ξ, rhess, 0.5)
    ηOpt = truncated_conjugate_gradient_descent(
        M, cost, (X) -> rgrad(M, X), x, ξ, rhess, 0.5; return_options=true
    )
    @test submanifold_components(get_solver_result(ηOpt)) == submanifold_components(η)
end


text:
int_rules_6_3_11 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 6.3.11*(e*x)^m*(a+b*tanh(c+d*x^n))^p =#
    @apply_utils Antiderivative((~(a') + ~(b') * tanh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (1 / ~n - 1) * (~a + ~b * tanh(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (IGtQ(1 / ~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(a') + ~(b') * coth(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (1 / ~n - 1) * (~a + ~b * coth(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (IGtQ(1 / ~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(a') + ~(b') * tanh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Integral((~a + ~b * tanh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * coth(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Integral((~a + ~b * coth(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * tanh(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * tanh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
    @apply_utils Antiderivative((~(a') + ~(b') * coth(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * coth(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
    @apply_utils Antiderivative((~(a') + ~(b') * tanh(~u)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * tanh(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(a') + ~(b') * coth(~u)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * coth(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * tanh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * tanh(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (IGtQ(Simplify((~m + 1) / ~n), 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * coth(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * coth(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (IGtQ(Simplify((~m + 1) / ~n), 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~x) ^ ~(m') * tanh(~(c') + ~(d') * (~x) ^ ~n) ^ 2, ~x) => (-((~x) ^ ((~m - ~n) + 1)) * tanh(~c + ~d * (~x) ^ ~n)) / (~d * ~n) + Antiderivative((~x) ^ ~m, ~x) + (((~m - ~n) + 1) / (~d * ~n)) * Antiderivative((~x) ^ (~m - ~n) * tanh(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * coth(~(c') + ~(d') * (~x) ^ ~n) ^ 2, ~x) => (-((~x) ^ ((~m - ~n) + 1)) * coth(~c + ~d * (~x) ^ ~n)) / (~d * ~n) + Antiderivative((~x) ^ ~m, ~x) + (((~m - ~n) + 1) / (~d * ~n)) * Antiderivative((~x) ^ (~m - ~n) * coth(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * tanh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Integral((~x) ^ ~m * (~a + ~b * tanh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * coth(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Integral((~x) ^ ~m * (~a + ~b * coth(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * tanh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * tanh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * coth(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * coth(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * tanh(~u)) ^ ~(p'), ~x) => Antiderivative((~e * ~x) ^ ~m * (~a + ~b * tanh(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~m, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * coth(~u)) ^ ~(p'), ~x) => Antiderivative((~e * ~x) ^ ~m * (~a + ~b * coth(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~m, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~x) ^ ~(m') * sech(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * tanh(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(q'), ~x) => (-((~x) ^ ((~m - ~n) + 1)) * sech(~a + ~b * (~x) ^ ~n) ^ ~p) / (~b * ~n * ~p) + (((~m - ~n) + 1) / (~b * ~n * ~p)) * Antiderivative((~x) ^ (~m - ~n) * sech(~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (RationalQ(~m) && (IntegerQ(~n) && (GeQ(~m - ~n, 0) && EqQ(~q, 1))))
    @apply_utils Antiderivative((~x) ^ ~(m') * csch(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * coth(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(q'), ~x) => (-((~x) ^ ((~m - ~n) + 1)) * csch(~a + ~b * (~x) ^ ~n) ^ ~p) / (~b * ~n * ~p) + (((~m - ~n) + 1) / (~b * ~n * ~p)) * Antiderivative((~x) ^ (~m - ~n) * csch(~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (RationalQ(~m) && (IntegerQ(~n) && (GeQ(~m - ~n, 0) && EqQ(~q, 1))))
end

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.3 Hyperbolic tangent/6.3.11 (e x)^m (a+b tanh(c+d x^n))^p.jl


text:
mutable struct VacuumModel{T <: SSDFloat} <: AbstractTemperatureModel{T} end

function VacuumModel(config_file::Dict; T::Type{<:AbstractFloat} = Float32)::VacuumModel
    return VacuumModel{T}(config_file::Dict)
end

function VacuumModel{T}(config_file::Dict; temperature::Union{Missing, T} = missing)::VacuumModel where {T <: SSDFloat}
    m = VacuumModel{T}()
end

function scale_to_given_temperature(m::VacuumModel{T})::NTuple{4,T} where {T <: SSDFloat}
    return T(1), T(1), T(1), T(1)
end


text:
struct PINningSparseSynapseParameter
end

@snn_kw mutable struct PINningSparseSynapse{VIT=Vector{Int32},VFT=Vector{Float32}}
    param::PINningSparseSynapseParameter = PINningSparseSynapseParameter()
    colptr::VIT # column pointer of sparse W
    I::VIT      # postsynaptic index of W
    W::VFT  # synaptic weight
    rI::VFT # postsynaptic rate
    rJ::VFT # presynaptic rate
    g::VFT  # postsynaptic conductance
    P::VFT  # <rᵢrⱼ>⁻¹
    q::VFT  # P * r
    f::VFT  # postsynaptic traget
    records::Dict = Dict()
end

"""
[PINing Sparse Synapse](https://www.ncbi.nlm.nih.gov/pubmed/26971945)
"""
PINningSparseSynapse

function PINningSparseSynapse(pre, post; σ = 1.5, p = 0.0, α = 1, kwargs...)
    w = σ / √(p * pre.N) * sprandn(post.N, pre.N, p)
    rowptr, colptr, I, J, index, W = dsparse(w)
    rI, rJ, g = post.r, pre.r, post.g
    P = α .* (I .== J)
    f, q = zeros(post.N), zeros(post.N)
    PINningSparseSynapse(;@symdict(colptr, I, W, rI, rJ, g, P, q, f)..., kwargs...)
end

function forward!(c::PINningSparseSynapse, param::PINningSparseSynapseParameter)
    @unpack colptr, I, W, rJ, g, P, q = c
    fill!(q, zero(Float32))
    fill!(g, zero(Float32))
    @inbounds for j in 1:(length(colptr) - 1)
        rJj = rJ[j]
        for s = colptr[j]:(colptr[j+1] - 1)
            i = I[s]
            g[i] += W[s] * rJj
            q[i] += P[s] * rJj
        end
    end
end

function plasticity!(c::PINningSparseSynapse, param::PINningSparseSynapseParameter, dt::Float32, t::Float32)
    @unpack colptr, I, W, rI, g, P, q, f = c
    C = 1 / (1 + dot(q, rI))
    @inbounds for j in 1:(length(colptr) - 1)
        for s in colptr[j]:(colptr[j+1] - 1)
            i = I[s]
            P[s] += -C * q[i] * q[j]
            W[s] += C * (f[i] - g[i]) * q[j]
        end
    end
end


text:
## operations on Tape

## Input

struct Input <: AbstractOp
    var::TAny
    argid::Int
end

@inline exec!(tape::Tape, op::Input) = op.var

function record!(tape::Tape, ::Type{Input}, val; argid=-1)
    op = Input(tracked(tape, val), argid)
    push!(tape.ops, op)
    op.var.id = length(tape)
    return op.var
end

## Constant

struct Constant <: AbstractOp
    var::TAny
end
constant(tape::Tape, x::TAny) = x
constant(tape::Tape, x) = record!(tape, Constant, x)

@inline exec!(tape::Tape, op::Constant) = op.var

function record!(tape::Tape, ::Type{Constant}, val)
    op = Constant(tracked(tape, val))
    push!(tape.ops, op)
    op.var.id = length(tape)
    return op.var
end

## Call

"""
Method call
"""
struct Call{Fn, ARGS <: Tuple} <: AbstractOp
    var::TAny                  # tracked var, result of the call
    fn::Fn                     # function to apply to args
    args::ARGS                 # call arguments
    kwargs::Dict{Symbol, Any}  # keyword arguments
    # Call(var::TAny, fn::Fn, args::ARGS; kwargs=Dict()) where {Fn, ARGS} =
    #     new{Fn, ARGS}(var, fn, args, kwargs)
end

function Base.show(io::IO, op::Call)
    args_str = join([var isa TAny ? "%$(var.id)" : var for var in op.args], ", ")
    kwargs_str = isempty(op.kwargs) ? "" : "; " * join(["$k=$v" for (k, v) in op.kwargs], ", ")
    print(io, "Call(%$(op.var.id) = $(op.fn)($(args_str)$kwargs_str))")
end

function record!(tape::Tape, ::Type{Call}, fn::Fn, args::ARGS;
                 kwargs=Dict{Symbol,Any}()) where {Fn, ARGS<:Tuple}
    arg_vals = map(getvalue, args)
    val = to_device(tape.device, fn(arg_vals...; kwargs...))
    var = tracked(tape, val)
    op = Call(var, fn, args, kwargs)
    _record!(tape, op)
    return var
end

"""
Execute operation on a tape, store result to op's var.
"""
function exec!(tape::Tape, op::Call)
    arg_data = map(getvalue, op.args)
    op.var.val = op.fn(arg_data...; op.kwargs...)
    return op.var
end

## Bcast

"""
Broadcasting
"""
struct Bcast{Fn, ARGS <: Tuple} <: AbstractOp
    var::TAny
    fn::Fn
    args::ARGS
end

function Base.show(io::IO, op::Bcast)
    args_str = join([var isa TAny ? "%$(var.id)" : var for var in op.args], ", ")
    print(io, "Bcast(%$(op.var.id) = $(op.fn).($args_str))")
end

function record!(tape::Tape, ::Type{Bcast}, fn::Fn, args::ARGS) where {Fn, ARGS<:Tuple}
    arg_vals = map(getvalue, args)
    # println("**** $(fn). | $(map(typeof, arg_vals)) | $(map(size, arg_vals))")
    val = fn.(arg_vals...)
    var = tracked(tape, val)
    op = Bcast(var, fn, args)
    _record!(tape, op)
    return var
end

"""
Execute operation on a tape, store result to op's var.
"""
function exec!(tape::Tape, op::Bcast)
    arg_data = map(getvalue, op.args)
    op.var.val = op.fn.(arg_data...)
    return op.var
end

## Assign

struct Assign <: AbstractOp
    var::TAny
    src::TAny
end

Base.show(io::IO, op::Assign) = print(io, "Assign(%$(op.var.id) = %$(op.src.id))")

function exec!(tape::Tape, op::Assign)
    op.var.val = op.src.val
    return op.var
end

function record!(tape::Tape, ::Type{Assign}, src::TAny)
    var = tracked(tape, src.val)
    op = Assign(var, src)
    _record!(tape, op)
    exec!(tape, op)
    return var
end

## COMMON UTILS

getvar(op::AbstractOp) = op.var
setvar!(op::AbstractOp, var::TAny) = (op.var = var)

getvalue(op::AbstractOp) = op |> getvar |> getvalue
setvalue!(op::AbstractOp, val) = setvalue!(getvar(op), val)

@inline Base.getproperty(op::AbstractOp, name::Symbol) = _getproperty(op::AbstractOp, Val(name))
@inline _getproperty(x::AbstractOp, ::Val{F}) where F = getfield(x, F)
@inline _getproperty(x::AbstractOp, ::Val{:val}) = getvalue(x)

## mutable structs: writing all trackable fields to the tape

"""
Traverse mutable struct and write all trackable fields to the tape,
keeping mapping from field paths to tracked vars.
"""
function record_struct!(tape::Tape, s, argid::Int; field_path=[])
    # S = typeof(s)
    S = Core.eval(Main, Symbol(typeof(s).name))  # parameter-free type of
    new_fields = []
    for name in fieldnames(S)
        full_field_path = vcat(field_path, name)
        field = getfield(s, name)
        is_trackable = (field isa Real && !isa(field, Bool)) || field isa AbstractArray
        is_concrete = fieldtype(S, name) |> isconcretetype
        if is_trackable && !is_concrete
            var = record!(tape, Input, field; argid=argid)
            push!(new_fields, var)
            # save mapping field_path -> var
            if !haskey(tape.sfields, argid)
                tape.sfields[argid] = Dict{Any,Any}()
            end
            full_field_path_tuple = (full_field_path...,)
            tape.sfields[argid][full_field_path_tuple] = var.id
        elseif isstruct(field)
            var = record_struct!(tape, field, argid; field_path=full_field_path)
            push!(new_fields, var)
        else
            # if field isn't trackable, take it as it is
            push!(new_fields, field)
        end
    end
    # assuming default constuctor
    return S(new_fields...)
end


text:
function Newton_Spectral(nlp:: AbstractNLPModel,
                         x₀ :: AbstractVector;
                         τ₀ :: Float64 = 0.0005,
                         ϵ  :: Float64 = 1e-6,
                         maxiter :: Int = 200)

    x = copy(x₀)
    iter = 0
    f, g = obj(nlp,x), grad(nlp, x)

    while (norm(g, Inf) > ϵ) && (iter <= maxiter)
        H    = Matrix(Symmetric(hess(nlp, x),:L))
        Δ, O = eigen(H)

        # Boost negative values of Δ to 1e-8
        D = Δ .+ max.((1e-8 .- Δ), 0.0)

        d = - O*diagm(1.0 ./ D)*O'*g

        # Simple Armijo backtracking
        hp0 = g'*d
        t   = 1.0
        ft  = obj(nlp, x + t*d)
        while ft > ( f + τ₀*t*hp0)
            t /= 2.0
            ft = obj(nlp, x + t*d)
        end
        x    += t*d
        f, g  = ft, grad(nlp, x)
        iter += 1
    end
    if iter > maxiter @warn "Iteration limit" end

    return x, f, norm(g), iter
end


text:
"""
    abs_val(num)

Program to find the absolute value of a number

# Example
```julia
abs_val(-100) # returns 100
abs_val(0) # returns 0
abs(123.1) # returns 123.1
-1000 == abs_val(-1000) #returns false
1000 == abs_val(1000) #returns true
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function abs_val(num)
    return num < 0 ? -num : num
end

"""
    abs_max(x)

Program to find the max absolute value in a vector

# Example
```julia
abs_max([1,3,4]) # returns 4
abs_max([-3,1,2]) # returns -3
abs_max([-7,-3,6]) #returns -7
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function abs_max(x)
    result = sort(x, by = abs; alg = QuickSort)
    return result[length(result)]
end

"""
    abs_min(num)

Program to find the min absolute value in a vector

# Example
```julia
abs_min([1,3,4]) # returns 1
abs_min([-3,1,2]) # returns 1
abs_min([-7,-3,6]) #returns -3
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function abs_min(x)
    result = sort(x, by = abs; alg = QuickSort)
    return result[1]
end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
abstract type StaticInjection <: Device end
function supports_services(::T) where {T <: Device}
    return true
end

function get_services(::Device)
    return Vector{Service}()
end

get_dynamic_injector(d::StaticInjection) = nothing


text:
using Animations
using GtkReactive
using Gtk: get_gtk_property, visible
using Images
import Interact
import Interact: @map, Widget, Widgets, @layout!, hbox, vbox
using Javis
import Latexify: latexify
using LaTeXStrings
using ReferenceTests
using Test
using VideoIO

function ground(video, action, framenumber)
    background("white")
    sethue("blue")
    return framenumber
end

function ground_black_on_white(video, action, framenumber)
    background("white")
    sethue("black")
    return framenumber
end

function ground_color(color_bg, color_pen, framenumber)
    background(color_bg)
    sethue(color_pen)
    return framenumber
end

function circle_with_color(position, radius, action, color)
    sethue(color)
    circle(position, radius, action)
end

@testset "Unit" begin
    include("unit.jl")
end
@testset "SVG LaTeX tests" begin
    include("svg.jl")
end
@testset "Animations" begin
    include("animations.jl")
end
@testset "Morphing" begin
    include("morphing.jl")
end
@testset "Javis Viewer" begin
    include("viewer.jl")
end
@testset "Shorthands" begin
    include("shorthands.jl")
end
@testset "Layers" begin
    include("layers.jl")
end
@testset "Postprocessing" begin
    include("postprocessing.jl")
end


text:
using Documenter, DocStringExtensions, PowerSimulationsDynamics

makedocs(
    modules = [PowerSimulationsDynamics],
    format = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename = "PowerSimulationsDynamics.jl",
    pages = Any[ # Compat: `Any` for 0.4 compat
        "Home" => "index.md",
        # "User Guide" => "man/guide.md",
        "Tutorials" => Any[
            "Tutorial 0: Data Creation" => "Examples/example_data.md",
            "Tutorial 1: OMIB" => "Examples/example_OMIB.md",
            "Tutorial 2: Dynamic Lines" => "Examples/example_lines.md",
        ],
        "Models" => Any[
            "Network" => "Models/network.md",
            "Reference Frames" => "Models/srf.md",
            "Generator" => "Models/gens.md",
            "Inverter" => "Models/inverters.md",
            "Small Signal" => "Models/small.md",
        ],
    ],
)

deploydocs(
    repo = "github.com/NREL-SIIP/PowerSimulationsDynamics.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#"],
)


text:
module FormulationLattice

import Base: &, |
export Literal, @Literal, @Literals, And, ∧, Or, ∨,
       isliteral, iscnf, isdnf, basicstep, dnf

abstract Clause
Base.show(io::IO, x::Clause) = print(io, x)

# TODO: add negation
immutable Literal <: Clause
    name::UTF8String
end

default_names = ["A","B","C","D","E","F","G","H","I","J","K"]
let counter = 0
    Literal() = (counter += 1; Literal(default_names[counter]))
end

macro Literal(name)
    @assert isa(name, Symbol)
    :($(esc(name)) = Literal($(string(name))))
end

macro Literals(args...)
    names = if isa(args,Symbol)
        [args]
    elseif isa(args,Expr)
        @assert args.head == :tuple
        args.args
    else
        if isa(args,Tuple{Expr})
            args[1].args
        else
            [args...]
        end
    end

    @assert all(x -> (isa(x,Symbol)), names)
    code = Expr(:block)
    lastline = Expr(:tuple)
    for name in names
        push!(code.args, :($(esc(name)) = Literal($(string(name)))))
        push!(lastline.args, :($(esc(name))))
    end
    :($code; $lastline)
end

Base.print(io::IO, l::Literal) = print(io, l.name)

# percolate
percolate{T<:Clause}(::Type{T}, cl) = Clause[cl]
percolate(::Type{Literal}, l::Literal) = Clause[l]
function percolate{T<:Clause}(::Type{T}, cl::T)
    x = map(x -> percolate(T, x), cl.clauses)
    vcat(x...)
end

immutable And <: Clause
    clauses::Vector{Clause}
end

(&)(clauses::Clause...) = And(clauses...)
const ∧ = &

immutable Or <: Clause
    clauses::Vector{Clause}
end
Or( clauses::Clause...) = Or(collect(clauses))
(|)(clauses::Clause...) = Or(clauses...)
const ∨ = |

for (T,connector) in [(:And," ∧ "), (:Or," ∨ ")]; @eval begin
    # helper that percolates like clauses up to toplevel, e.g.
    # A ∧ (B ∧ C) to A ∧ B ∧ C
    function $T(clauses::Clause...)
        length(clauses) > 1 || error("Unary $($T) not allowed")
        terms = Clause[]
        for cl in clauses
            append!(terms, percolate($T, cl))
        end
        $T(terms)
    end

    function Base.print(io::IO, cl::$T)
        numclause = length(cl.clauses)
        if numclause == 0
            print(io,)
        elseif numclause == 1
            print(io,cl.clauses[1])
        else
            inner = join(map(x->sprint(print,x), cl.clauses), $connector)
            print(io, string("(", inner, ")"))
        end
    end
end; end

isliteral(c::Clause) = isa(c, Literal)
isand(c::Clause)     = isa(c,And)
isor(c::Clause)      = isa(c,Or)

Base.copy(x::Literal)      = Literal(copy(x.name))
Base.copy{C<:Clause}(x::C) = C(copy(x.clauses))

# Not sure if this logic holds exactly as I want:
# will probably want to simplify expressions like
# (A ⋁ (B ⋁ C)) ⋀ D
# to
# (A ⋁ B ⋁ C) ⋀ D
# Probably add "sanitizers" to inner constructors to clauses above?
iscnf(c::Literal) = true
iscnf(c::Or)      = all(isliteral, c.clauses)
iscnf(c::And)     = all(c.clauses) do cl
    isliteral(cl) | (isor(cl) && all(isliteral, cl.clauses))
end

isdnf(c::Literal) = true
isdnf(c::Or)      = all(c.clauses) do cl
    isliteral(cl) | (isand(cl) && all(isliteral, cl.clauses))
end
isdnf(c::And)     = all(isliteral, c.clauses)

include("formulations.jl")

function dnf(c::Clause, formtrack=FormulaState[])
    # Compute disjunctive normal form of clause c
    cl = c
    push!(formtrack, FormulaState(cl))
    while !isdnf(cl)
        cl = basicstep(cl)
        push!(formtrack, FormulaState(cl))
    end
    cl
end

function basicstep(c::Or)
    isdnf(c) && return c
    n = length(c.clauses)
    idx = findfirst(x -> !isdnf(x), c.clauses)
    newterm = basicstep(c.clauses[idx])
    before = c.clauses[1:(idx-1)]
    after  = c.clauses[(idx+1):n]
    ret = Or(before..., newterm, after...)
    ret
end

function basicstep(c::And)
    idx = findfirst(x -> isa(x,Or), c.clauses)
    ret = _basicstep(c, idx)
    ret
end

function _basicstep(c::And, idx::Int)
    n = length(c.clauses)
    @assert 1 ≤ idx ≤ n
    cl = c.clauses[idx]
    remaining = c.clauses[setdiff(1:n, idx)]
    @assert isa(cl, Or)
    tip  = cl.clauses[1]
    tail = cl.clauses[2:end]
    before = c.clauses[1:(idx-1)]
    after  = c.clauses[(idx+1):n]
    if length(tail) == 1
        And(before..., tip, after...) ∨ And(before..., tail[1], after...)
    else
        And(before..., tip, after...) ∨ And(before..., Or(tail...), after...)
    end
end

end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using SimilaritySearch
using Test

function test_vectors(create_index, dist::Function, ksearch, nick)
    @testset "indexing vectors with $nick and $dist" begin
        n = 300 # number of items in the dataset
        m = 30  # number of queries
        dim = 3  # vector's dimension

        db = [rand(Float32, dim) for i in 1:n]
        queries = [rand(Float32, dim) for i in 1:m]

        index = create_index(db)
        # testing pushes
        push!(index, dist, rand(Float32, dim))
        @test length(index.db) == n + 1
        perf = Performance(dist, index.db, queries, expected_k=10)
        p = probe(perf, index, dist)
        @show dist, p
        return p
    end
end

function test_sequences(create_index, dist::Function, ksearch, nick)
    @testset "indexing sequences with $nick and $dist" begin
        n = 300 # number of items in the dataset
        m = 30  # number of queries
        dim = 5  # the length of sequences
        V = collect(1:10)  # vocabulary of the sequences

        function create_item()
            s = rand(V, dim)
            if dist == jaccard_distance || dist == dice_distance || dist == intersection_distance
                sort!(s)
                s = unique(s)
            end

            return s
        end
        db = [create_item() for i in 1:n]
        queries = [create_item() for i in 1:m]

        @info "inserting items into the index"
        # index = Laesa(db, dist, σ)
        index = create_index(db)
        # optimize!(index, recall=0.9, use_distances=true)
        @test length(index.db) == n
        perf = Performance(dist, index.db, queries, expected_k=10)
        p = probe(perf, index, dist)
        # @show dist, p
        return p
    end
end

@testset "indexing vectors" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are required
    ksearch = 10
    local index

    for (recall_lower_bound, dist) in [
        (1.0, l2_distance), # 1.0 -> metric, < 1.0 if dist is not a metric
        (1.0, l1_distance),
        (1.0, linf_distance),
        (1.0, lp_distance(3.0)),
        (0.1, lp_distance(0.5))
    ]
        @show recall_lower_bound, dist

        @show Kvp
        p = test_vectors((db) -> fit(Kvp, dist, db, 3, 32), dist, ksearch, "Kvp")
        @test p.recall >= recall_lower_bound * 0.99 # not 1 to allow some "numerical" deviations
    end
end

@testset "indexing sequences" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are required
    ksearch = 10
    local index

    # metric distances should achieve recall=1 (perhaps lesser because of numerical inestability)
    for (recall_lower_bound, dist) in [
        (1.0, jaccard_distance),
        (0.1, dice_distance),
        (0.1, intersection_distance),
        (0.1, common_prefix_distance),
        (1.0, levenshtein_distance),
        (1.0, lcs_distance),
        (1.0, hamming_distance)
    ]
        p = test_sequences((db) -> fit(Kvp, dist, db, 3, 32), dist, ksearch, "Kvp")
        @test p.recall >= recall_lower_bound * 0.99  # not 1 to allow some "numerical" deviations
    end
end


text:
using PkgSRA, Plots, LinearAlgebra, Statistics
using AxisArrays, DataFrames, DataStructures
using BSON: @load
pyplot()

include("paper_plotting.jl")

#####
##### Load LORENZ variables (special top panel)
#####
include("paper_settings.jl")
this_dat_name = DAT_FOLDERNAME*"TMP_dat_flowchart_inkscape_"

# Uncontrolled data
fname = this_dat_name*"uncontrolled_lorenz.bson"
@load fname dat_raw numerical_grad_raw

# Original ODE and controller variables
fname = this_dat_name*"ode_vars.bson";
@load fname dat dat_grad true_grad U_true

# Bayesian variables 1: before control
# fname = this_dat_name*"naive_vars_sindy.bson";
# @load fname sindy_dat_unctr sindy_unctr

# Bayesian variables 1: before control
fname = this_dat_name*"naive_vars_bayes.bson";
@load fname dat_unctr sample_trajectory_noise residual accepted_ind ctr_guess

# SINDY variables 2: after control
fname = this_dat_name*"ctr_vars_sindy.bson";
@load fname sindy_dat_ctr sindy_grad_ctr

# Bayesian variables 2: after control
fname = this_dat_name*"ctr_vars_turing.bson";
@load fname dat_ctr sample_trajectory_noise_ctr sample_gradients_ctr turing_dat_ctr ctr_guess2

#####
##### Plot everything
#####
this_fig_name = FIGURE_FOLDERNAME*"fig_library_of_examples_"

## 1: Attractor view (uncontrolled)
plot_raw = plot3d(dat_raw[1, :], dat_raw[2, :], dat_raw[3, :],
        color=COLOR_DICT["data_uncontrolled"], lw=4; plot_opt...)
        # xlabel!("Lorenz");
        title!("Lorenz")

## 2: 1d time series (controlled)
plot_coordinate = 1;
    ind = 501:1000;
plot_1d = plot(dat[plot_coordinate, ind],
        color=COLOR_DICT["data"], lw=4; plot_opt...)
        ylabel!("X", guidefontsize=24)
        title!("Time Series");
plot_ctr_true = plot(U_true[plot_coordinate,ind],
        color=COLOR_DICT["true"], lw=3;
        plot_opt...);

## 3: Residual with noise envelope (x coordinate)
plot_residual1 = plot(residual[plot_coordinate,ind],
                    ribbon=mean(sample_trajectory_noise),
                    fillalpha=0.5,
                    color=COLOR_DICT["residual"], lw=4;
                    plot_opt...)
    ylabel!("dX/dt", guidefontsize=24)
    title!("Residual");

## 4: Control signal guess (x coordinate)
# Same ind as above
plot_ctr_guess = plot(ctr_guess[plot_coordinate,ind],
                color=COLOR_DICT["control_time"], lw=6;
                plot_opt...)
        title!("Learned Controller");

## 5: "Revealed underlying system"
controlled_ind = 1:5000;
    ind = controlled_ind;
let d = sindy_dat_ctr
    global plot_reconstruction = plot3d(
                d[1, ind], d[2, ind], d[3, ind],
                color=COLOR_DICT["model_controlled"], alpha=0.8, lw=3;
                plot_opt...);
end
    title!("Learned System");

##
## Put it all together
##
lay = @layout [a b c [d; e] f]

plot_final = plot(
        plot_raw, plot_1d, #plot_ctr_true,
        plot_residual1, plot_ctr_guess, plot_ctr_true, plot_reconstruction,
        layout=lay)
plot!(size=(2000, 250))
fname = this_fig_name * "lorenz.pdf";
savefig(plot_final, fname)

#####
##### Second, default-style panel: ROSSLER
#####
this_dat_name = DAT_FOLDERNAME*"TMP_dat_library_of_examples_rossler";
this_fig_name = FIGURE_FOLDERNAME*"fig_library_of_examples_";
plot_final = plot_library(this_dat_name, this_fig_name, "Rossler",
                plot_ind=1001:2000,
                plot_coordinate=3)

fname = this_fig_name * "rossler.pdf";
savefig(plot_final, fname)

#####
##### Van der Pol oscillator
#####
this_dat_name = DAT_FOLDERNAME*"dat_library_of_examples_vdp_";
this_fig_name = FIGURE_FOLDERNAME*"fig_library_of_examples_";
plot_final = plot_library(this_dat_name, this_fig_name, "Van der Pol",
                plot_ind=1001:2000,
                plot_coordinate=1)

fname = this_fig_name * "vanDerPol.pdf";
savefig(plot_final, fname)

#####
##### FitzHugh Nagumo oscillator
#####
this_dat_name = DAT_FOLDERNAME*"dat_library_of_examples_fhn_";
this_fig_name = FIGURE_FOLDERNAME*"fig_library_of_examples_";
plot_final = plot_library(this_dat_name, this_fig_name, "FitzHugh Nagumo",
                plot_ind=1001:2000,
                plot_coordinate=1)

fname = this_fig_name * "fitzHughNagumo.pdf";
savefig(plot_final, fname)

#####
##### Lotka-Volterra predator prey system
#####
this_dat_name = DAT_FOLDERNAME*"dat_library_of_examples_lv_";
this_fig_name = FIGURE_FOLDERNAME*"fig_library_of_examples_";
plot_final = plot_library(this_dat_name, this_fig_name, "Lotka-Volterra",
                plot_ind=2001:3000,
                plot_coordinate=1)

fname = this_fig_name * "lotkaVolterra.pdf";
savefig(plot_final, fname)


text:
function test_phidf(order, α)
    srand(98765)
    dummy_input = rand(64)

    delay = zeros(order + 1)
    for x in dummy_input
        SPTK.phidf!(x, order, α, delay)
        @test all(isfinite.(delay))
    end
end

function test_lspcheck(order)
    srand(98765)
    dummy_lsp = rand(1024)

    # TODO valid check
    lspcheck(dummy_lsp)
end

println("-- test_phidf")
for order in [20, 22, 24]
    for α in [0.35, 0.41, 0.5]
        println(" where order = $order, α = $α")
        test_phidf(order, α)
    end
end

let
    println("-- test_phidf exceptions")
    c = zeros(21)
    delay = zeros(length(c)-1)
    @test_throws ArgumentError SPTK.phidf!(0.0, length(c)-1, 0.41, delay)
end

println("-- test_lspcheck")
for order in [20, 22, 24]
    println(" where order = $order")
    test_lspcheck(order)
end


text:
using CameraUnits, Unitful, FixedPointNumbers
using Base.Test

@testset "CameraUnits" begin
    x = 1cu"pco55"
    @test uconvert(γe, x) == 0.46γe
end

nothing


text:
using Base.Test
# pointwise.jl
import WaveletScattering: Identity, Log1P
# blob.jl
import WaveletScattering: ScatteredBlob
# inputlayer.jl
import WaveletScattering: InputLayer
# inputlayerstate.jl
import WaveletScattering: InputLayerState, kthrange
# node.jl
import WaveletScattering: Node
# path.jl
import WaveletScattering: Path

# Identity
ρ = Identity()
@test ρ([1.0, 2.0, 3.0]) ≈ [1.0, 2.0, 3.0]

# Log1P
threshold = 0.1
ρ = Log1P(threshold)
@test ρ([0.0]) ≈ [0.0]
@test expm1.(ρ([1.0, 2.0, 3.0])) ≈ [1.0, 2.0, 3.0]*threshold
@test_throws DomainError ρ([- 1.001 / threshold])

# Base.map
layer = InputLayer(
        tops = [:signal],
        symbols = [:time, :chunk],
        data = map(Float32, randn(256, 2)))
ranges = ntuple(k -> kthrange(layer, k), ndims(layer.data))
blob = ScatteredBlob(
    DataStructures.SortedDict((Path() => Node(layer.data, ranges),)))
ρ = Identity()
outblob = ρ(blob)


text:
# ---
# title: 99. Recover Binary Search Tree
# id: problem99
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/recover-binary-search-tree/description/>
# hidden: true
# ---
#
# You are given the `root` of a binary search tree (BST), where exactly two
# nodes of the tree were swapped by mistake. _Recover the tree without changing
# its structure_.
#
# **Follow up:** A solution using `O(n)` space is pretty straight forward. Could
# you devise a constant space solution?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/28/recover1.jpg)
#
#
#
#     Input: root = [1,3,null,null,2]
#     Output: [3,1,null,null,2]
#     Explanation: 3 cannot be a left child of 1 because 3 > 1. Swapping 1 and 3 makes the BST valid.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/28/recover2.jpg)
#
#
#
#     Input: root = [3,1,4,null,null,2]
#     Output: [2,1,4,null,null,3]
#     Explanation: 2 cannot be in the right subtree of 3 because 2 < 3. Swapping 2 and 3 makes the BST valid.
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[2, 1000]`.
#   * `-231 <= Node.val <= 231 - 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function gmres(a::AbstractArray{T,2},b::AbstractArray{T,1},x::AbstractArray{T,1}) where {T <: AbstractFloat} # Exact GMRES

    # Create Hessenberg form using Arnoldi iteration
    q, h, c = arnoldi(a,b)

    e1    = zeros(c)
    e1[1] = norm(b)

    g = zeros(T,2,2)
    # Use Givens rotations to transform to upper-triangular form
    @views for i = 2:c
        givens!(g,h[i-1:i,i-1])
        h[i-1:i,:] = g'h[i-1:i,:]
        e1[i-1:i] = g'e1[i-1:i]
    end

    step = q*(h\e1)

    return step

end

function gmres(a::AbstractArray{T,2},b::AbstractArray{T,1},x::AbstractArray{T,1},forcing_term::T) where {T <: AbstractFloat} # Inexact GMRES

    # Create Hessenberg form using Arnoldi iteration
    q, h, c = arnoldi(a,b,forcing_term*norm(b))

    e1    = zeros(c)
    e1[1] = norm(b)

    g = Array{T,2}(undef,2,2)
    # Use Givens rotations to transform to upper-triangular form
    @views for i = 2:c
        givens!(g,h[i-1:i,i-1])
        h[i-1:i,:] = g'h[i-1:i,:]
        e1[i-1:i] = g'e1[i-1:i]
    end

    step = q*(h\e1)

    return step

end

function gmres(a::AbstractArray{T,2},b::AbstractArray{T,1},x::AbstractArray{T,1},forcing_term::T,m::S) where {T <: AbstractFloat, S <: Integer} # Restarted inexact GMRES

    x = copy(x)
    n = length(x)

    m = min(m,n)

    step = Array{T,1}(undef,n)

    g = Array{T,2}(undef,2,2)

    loop_max = ceil(Int,n/m)
    loop_count = 1
    while loop_count <= loop_max

        # Create Hessenberg form using Arnoldi iteration
        q, h, c = arnoldi(a,b,forcing_term*norm(b),m)

        e1    = zeros(c)
        e1[1] = norm(b)

        # Use Givens rotations to transform to upper-triangular form
        @views for i = 2:c
            givens!(g,h[i-1:i,i-1])
            h[i-1:i,:] = g'h[i-1:i,:]
            e1[i-1:i] = g'e1[i-1:i]
        end

        step .= q*(h\e1)

        if maximum(abs,a*step-b) < forcing_term*norm(b)
            return step, true
        else
            x .= x + step
        end
        loop_count += 1
    end

    return step, false

end

function gmres(f::Function,x::Array{T,1}) where {T <: AbstractFloat} # Exact GMRES

    n = length(x)
    a = ForwardDiff.jacobian(f,x)
    b = -f(x)

    # Create Hessenberg form using Arnoldi iteration
    q, h, c = arnoldi(a,b)

    e1    = zeros(c)
    e1[1] = norm(b)

    g = Array{T,2}(undef,2,2)
    # Use Givens rotations to transform to upper-triangular form
    @views for i = 2:c
        givens!(g,h[i-1:i,i-1])
        h[i-1:i,:] = g'h[i-1:i,:]
        e1[i-1:i] = g'e1[i-1:i]
    end

    step = q*(h\e1)

    return step

end

function gmres_sparse(f::Function,x::Array{T,1}) where {T <: AbstractFloat} # Exact GMRES

    n = length(x)
    a = sparse(ForwardDiff.jacobian(f,x))
    b = -f(x)

    # Create Hessenberg form using Arnoldi iteration
    q, h, c = arnoldi(a,b)

    e1    = zeros(c)
    e1[1] = norm(b)

    g = Array{T,2}(undef,2,2)
    # Use Givens rotations to transform to upper-triangular form
    @views for i = 2:c
        givens!(g,h[i-1:i,i-1])
        h[i-1:i,:] = g'h[i-1:i,:]
        e1[i-1:i] = g'e1[i-1:i]
    end

    step = q*(h\e1)

    return step

end

function gmres(f::Function,x::Array{T,1},forcing_term::T) where {T <: AbstractFloat} # Inexact GMRES

    n = length(x)
    a = ForwardDiff.jacobian(f,x)
    b = -f(x)

    # Create Hessenberg form using Arnoldi iteration
    q, h, c = arnoldi(a,b,forcing_term*norm(b))

    e1    = zeros(c)
    e1[1] = norm(b)

    g = Array{T,2}(undef,2,2)
    # Use Givens rotations to transform to upper-triangular form
    @views for i = 2:c
        givens!(g,h[i-1:i,i-1])
        h[i-1:i,:] = g'h[i-1:i,:]
        e1[i-1:i] = g'e1[i-1:i]
    end

    step = q*(h\e1)

    return step

end

function gmres_sparse(f::Function,x::Array{T,1},forcing_term::T) where {T <: AbstractFloat} # Inexact GMRES

    n = length(x)
    a = sparse(ForwardDiff.jacobian(f,x))
    b = -f(x)

    # Create Hessenberg form using Arnoldi iteration
    q, h, c = arnoldi(a,b,forcing_term*norm(b))

    e1    = zeros(c)
    e1[1] = norm(b)

    g = Array{T,2}(undef,2,2)
    # Use Givens rotations to transform to upper-triangular form
    @views for i = 2:c
        givens!(g,h[i-1:i,i-1])
        h[i-1:i,:] = g'h[i-1:i,:]
        e1[i-1:i] = g'e1[i-1:i]
    end

    step = q*(h\e1)

    return step

end

function gmres(f::Function,x::Array{T,1},forcing_term::T,m::S) where {T <: AbstractFloat, S <: Integer} # Restarted inexact GMRES

    x = copy(x)
    n = length(x)

    step = similar(x)

    m = min(m,n)

    j = Array{T,2}(undef,n,n)
    b = Array{T,1}(undef,n)

    g = Array{T,2}(undef,2,2)

    loop_max = ceil(Int,n/m)
    loop_count = 1
    while loop_count <= loop_max

        j .= ForwardDiff.jacobian(f,x)
        b .= -f(x)

        # Create Hessenberg form using Arnoldi iteration
        q, h, c = arnoldi(j,b,forcing_term*norm(b),m)

        e1    = zeros(c)
        e1[1] = norm(b)

        # Use Givens rotations to transform to upper-triangular form
        @views for i = 2:c
            givens!(g,h[i-1:i,i-1])
            h[i-1:i,:] = g'h[i-1:i,:]
            e1[i-1:i] = g'e1[i-1:i]
        end

        step .= q*(h\e1)

        if maximum(abs,j*step-b) < forcing_term*norm(b)
            return step, true
        else
            x .= x + step
        end
        loop_count +=1
    end

    return step, false

end

function gmres_sparse(f::Function,x::Array{T,1},forcing_term::T,m::S) where {T <: AbstractFloat, S <: Integer} # Restarted inexact GMRES

    x = copy(x)
    n = length(x)

    step = similar(x)

    m = min(m,n)

    j = sparse(Array{T,2}(undef,n,n))
    b = Array{T,1}(undef,n)

    g = Array{T,2}(undef,2,2)

    loop_max = ceil(Int,n/m)
    loop_count = 1
    while loop_count <= loop_max

        j .= sparse(ForwardDiff.jacobian(f,x))
        b .= -f(x)

        # Create Hessenberg form using Arnoldi iteration
        q, h, c = arnoldi(j,b,forcing_term*norm(b),m)

        e1    = zeros(c)
        e1[1] = norm(b)

        # Use Givens rotations to transform to upper-triangular form
        @views for i = 2:c
            givens!(g,h[i-1:i,i-1])
            h[i-1:i,:] = g'h[i-1:i,:]
            e1[i-1:i] = g'e1[i-1:i]
        end

        step .= q*(h\e1)

        if maximum(abs,j*step-b) < forcing_term*norm(b)
            return step, true
        else
            x .= x + step
        end
        loop_count +=1
    end

    return step, false

end

function jacobian_free_gmres(f::Function,x::Array{T,1}) where {T <: AbstractFloat} # Jacobian-free exact GMRES

    b = -f(x)

    # Create Hessenberg form using Arnoldi iteration
    q, h, c = arnoldi(f,x,b)

    e1    = zeros(c)
    e1[1] = norm(b)

    g = Array{T,2}(undef,2,2)
    # Use Givens rotations to transform to upper-triangular form
    @views for i = 2:c
        givens!(g,h[i-1:i,i-1])
        h[i-1:i,:] = g'h[i-1:i,:]
        e1[i-1:i] = g'e1[i-1:i]
    end

    step = q*(h\e1)

    return step

end

function jacobian_free_gmres(f::Function,x::Array{T,1},forcing_term::T) where {T <: AbstractFloat} # Jacobian-free inexact GMRES

    b = -f(x)

    # Create Hessenberg form using Arnoldi iteration
    q, h, c = arnoldi(f,x,b,forcing_term*norm(b))

    e1    = zeros(c)
    e1[1] = norm(b)

    g = Array{T,2}(undef,2,2)
    # Use Givens rotations to transform to upper-triangular form
    @views for i = 2:c
        givens!(g,h[i-1:i,i-1])
        h[i-1:i,:] = g'h[i-1:i,:]
        e1[i-1:i] = g'e1[i-1:i]
    end

    step = q*(h\e1)

    return step

end

function jacobian_free_gmres(f::Function,x::Array{T,1},forcing_term::T,m::S) where {T <: AbstractFloat, S <: Integer} # Jacobian-free restarted inexact GMRES

    x = copy(x)
    n = length(x)
    jvec = zeros(n)
    step = similar(x)

    m = min(m,n)

    b = Array{T,1}(undef,n)

    g = Array{T,2}(undef,2,2)

    loop_max = ceil(Int,n/m)
    loop_count = 1
    while loop_count <= loop_max

        b .= -f(x)

        # Create Hessenberg form using Arnoldi iteration
        q, h, c = arnoldi(f,x,b,forcing_term*norm(b),m)

        e1    = zeros(c)
        e1[1] = norm(b)

        # Use Givens rotations to transform to upper-triangular form
        @views for i = 2:c
            givens!(g,h[i-1:i,i-1])
            h[i-1:i,:] = g'h[i-1:i,:]
            e1[i-1:i] = g'*e1[i-1:i]
        end

        step .= q*(h\e1)

        jacvec!(jvec,f,x,step)

        if maximum(abs,jvec-b) < forcing_term*norm(b)
            return step, true
        else
            x .= x + step
        end
        loop_count += 1
    end

    return step, false

end

function jacobian_free_gmres_inplace(f::Function,x::Array{T,1},forcing_term::T,m::S) where {T <: AbstractFloat, S <: Integer} # Jacobian-free restarted inexact GMRES

    x = copy(x)
    n = length(x)
    jvec = zeros(n)
    step = similar(x)

    m = min(m,n)

    b   = Array{T,1}(undef,n)
    ffk = Array{T,1}(undef,n)

    g = Array{T,2}(undef,2,2)

    loop_max = ceil(Int,n/m)
    loop_count = 1
    while loop_count <= loop_max

        f(ffk,x)
        b .= -ffk

        # Create Hessenberg form using Arnoldi iteration
        q, h, c = arnoldi_inplace(f,x,b,forcing_term*norm(b),m)

        e1    = zeros(c)
        e1[1] = norm(b)

        # Use Givens rotations to transform to upper-triangular form
        @views for i = 2:c
            givens!(g,h[i-1:i,i-1])
            h[i-1:i,:] = g'h[i-1:i,:]
            e1[i-1:i] = g'*e1[i-1:i]
        end

        step .= q*(h\e1)

        jacvec_inplace!(jvec,f,x,step)

        if maximum(abs,jvec-b) < forcing_term*norm(b)
            return step, true
        else
            x .= x + step
        end
        loop_count += 1
    end

    return step, false

end


text:
import Base.println, Base.flush

FILE_SIZE = 5000 * 1028

type FileRoller <: IO
    prefix::AbstractString
    folder::AbstractString
    filepath::AbstractString
    file::IO
    byteswritten::Int64
end

function getsuffix(n::Integer)
    str = string(n)
    for i = 1:(4 - length(str))
        str = "0"*str
    end
    str
end

function getfile(folder::AbstractString, prefix::AbstractString)
    i = 1
    getpath(i) = joinpath(folder, "$(prefix).$(getsuffix(i))")
    p = getpath(i)
    while isfile(p)
        p = getpath(i)
        i += 1
    end
    p, open(p, "a")
end

FileRoller(prefix) = FileRoller(prefix, pwd())

FileRoller(prefix, dir) = FileRoller(prefix, dir, (getfile(dir, prefix))..., 0)

function println(f::FileRoller, s::AbstractString)
    if f.byteswritten > FILE_SIZE
        gf = getfile(f.folder, f.prefix)
        f.filepath = gf[1]
        f.file = gf[2]
        f.byteswritten = 0
    end
    f.byteswritten += write(f.file, "$s\n")
end

flush(f::FileRoller) = flush(f.file)


text:
export HEVISolverType

"""
# Description
    HEVISolverType(FT;
        solver_method = ARK2ImplicitExplicitMidpoint,

        linear_max_subspace_size = Int(30)
        linear_atol = FT(-1.0)
        linear_rtol = FT(5e-5)

        nonlinear_max_iterations = Int(10)
        nonlinear_rtol = FT(1e-4)
        nonlinear_ϵ = FT(1.e-10)
        preconditioner_update_freq = Int(50)
    )

This solver type constructs a solver for ODEs with the
additively horizontal explicit vertical explicit~(HEVI) partitioned form.
the equation is assumed to be decomposed as

```math
  \\dot{Q} = [l(Q, t)] + [f(Q, t) - l(Q, t)]
```

where `Q` is the state, `f` is the full tendency and `l` is the vertical implicit
operator. The splitting is done automatically.

# Arguments
- `solver_method` (Function): Function defining the particular additive
    Runge-Kutta method to be used for the HEVI method.
    Default: `ARK2ImplicitExplicitMidpoint`
- `linear_max_subspace_size` (Int): maximal dimension of each (batched)
    Krylov subspace. GEMRES, a iterative linear solver is applied
    Default: `30`
- `linear_atol` (FT): absolute tolerance for linear solver convergence.
    Default: `-1.0`
- `linear_rtol` (FT): relative tolerance for linear solver convergence.
    Default: `5.0e-5`
- `nonlinear_max_iterations` (Int): max number of Newton iterations
    Default: `10`
- `nonlinear_rtol` (FT): relative tolerance for nonlinear solver convergence.
    Default: `1e-4`
- `nonlinear_ϵ` (FT): parameter denoting finite different step size for the
   Jacobian approximation.
   Default: `1e-10`
- `preconditioner_update_freq` (Int): Int denoting how frequent you need
    to update the preconditioner
    -1: no preconditioner;
    positive number, update every freq times.
    Default: `50`
"""
struct HEVISolverType{FT} <: AbstractSolverType
    # Function for the HEVI method
    solver_method::Function

    linear_max_subspace_size::Int
    linear_atol::FT
    linear_rtol::FT

    nonlinear_max_iterations::Int
    nonlinear_rtol::FT
    nonlinear_ϵ::FT

    preconditioner_update_freq::Int

    function HEVISolverType(
        FT;
        solver_method = ARK2ImplicitExplicitMidpoint,
        linear_max_subspace_size = Int(30),
        linear_atol = FT(-1.0),
        linear_rtol = FT(5e-5),
        nonlinear_max_iterations = Int(10),
        nonlinear_rtol = FT(1e-4),
        nonlinear_ϵ = FT(1.e-10),
        preconditioner_update_freq = Int(50),
    )

        return new{FT}(
            solver_method,
            linear_max_subspace_size,
            linear_atol,
            linear_rtol,
            nonlinear_max_iterations,
            nonlinear_rtol,
            nonlinear_ϵ,
            preconditioner_update_freq,
        )
    end
end

"""
    getdtmodel(ode_solver::HEVISolverType, bl)

A function which returns a model representing the dynamics
with the most restrictive time-stepping requirements.
"""
function getdtmodel(ode_solver::HEVISolverType, bl)
    # Most restrictive dynamics are treated implicitly
    return bl
end

"""
# Description
    solversetup(
        ode_solver::HEVISolverType{FT},
        dg,
        Q,
        dt,
        t0,
        diffusion_direction,
    )

Creates an ODE solver using a HEVI-type time-integration
scheme. All horizontal acoustic waves are treated explicitly,
while the 1-D vertical problem is treated implicitly. All
other dynamics (advection, diffusion) is treated explicitly
in the additive Runge-Kutta method.

# Comments:
Currently, the only HEVI-type splitting ClimateMachine can currently
do only involves splitting the acoustic processes; it is not currently
possible to perform more fine-grained separation of tendencies
(for example, including vertical advection or diffusion in the 1-D implicit problem)
"""
function solversetup(
    ode_solver::HEVISolverType,
    dg,
    Q,
    dt,
    t0,
    diffusion_direction,
)
    # All we need to do is create a DGModel for the
    # vertical acoustic waves (determined from the `implicit_model`)
    vdg = DGModel(
        dg.balance_law,
        dg.grid,
        dg.numerical_flux_first_order,
        dg.numerical_flux_second_order,
        dg.numerical_flux_gradient;
        state_auxiliary = dg.state_auxiliary,
        direction = VerticalDirection(),
        diffusion_direction = VerticalDirection(),
        check_for_crashes = dg.check_for_crashes,
    )

    # linear solver relative tolerance rtol which should be slightly smaller than the nonlinear solver tol
    linearsolver = BatchedGeneralizedMinimalResidual(
        dg,
        Q;
        max_subspace_size = ode_solver.linear_max_subspace_size,
        atol = ode_solver.linear_atol,
        rtol = ode_solver.linear_rtol,
    )

    # N(q)(Q) = Qhat  => F(Q) = N(q)(Q) - Qhat
    #
    # F(Q) == 0
    # ||F(Q^i) || / ||F(Q^0) || < tol
    # ϵ is a sensity parameter for this problem, it determines the finite difference Jacobian dF = (F(Q + ϵdQ) - F(Q))/ϵ
    # I have also try larger tol, but tol = 1e-3 does not work
    nonlinearsolver = JacobianFreeNewtonKrylovSolver(
        Q,
        linearsolver;
        tol = ode_solver.nonlinear_rtol,
        ϵ = ode_solver.nonlinear_ϵ,
        M = ode_solver.nonlinear_max_iterations,
    )

    # this is a second order time integrator, to change it to a first order time integrator
    # change it ARK1ForwardBackwardEuler, which can reduce the cost by half at the cost of accuracy
    # and stability
    # preconditioner_update_freq = 50 means updating the preconditioner every 50 Newton solves,
    # update it more freqent will accelerate the convergence of linear solves, but updating it
    # is very expensive
    solver = ode_solver.solver_method(
        dg,
        vdg,
        NonLinearBackwardEulerSolver(
            nonlinearsolver;
            isadjustable = true,
            preconditioner_update_freq = ode_solver.preconditioner_update_freq,
        ),
        Q;
        dt = dt,
        t0 = t0,
        split_explicit_implicit = false,
        variant = NaiveVariant(),
    )

    return solver
end


text:
function get_prob(matrix,str)
	# Grabbing resonator state
  level = eval(parse(str[search(str,',')+1:end-1]))
  # Grabbig qubit state
  if str[search(str,'|')+1] == 'g'
    return(matrix[:, level*2+1, level*2+1])
  else
   return(matrix[:, level*2+2, level*2+2])
  end
end