text:
using Colors, FixedPointNumbers, Test
using InteractiveUtils # for `subtypes`

@testset "Utilities" begin
    @testset "Mat3x3" begin
        M3x3 = Colors.Mat3x3([1 2 3; 4 5 6; 7 8 9])
        @test @inferred(M3x3[4]) === 2
        @test @inferred(M3x3[2,3]) === 6
        @test @inferred(M3x3[CartesianIndex(3, 2)]) === 8
    end

    @test Colors.cbrt01(0.6) ≈ cbrt(big"0.6") atol=eps(1.0)
    @test Colors.cbrt01(0.6f0) === cbrt(0.6f0)

    # issue #351
    xs = max.(rand(1000), 1e-4)
    @noinline l_pow_x_y() = for x in xs; x^2.4 end
    @noinline l_pow12_5() = for x in xs; Colors.pow12_5(x) end
    l_pow_x_y(); t_pow_x_y = @elapsed l_pow_x_y()
    l_pow12_5(); t_pow12_5 = @elapsed l_pow12_5()
    if t_pow12_5 > t_pow_x_y
        @warn "Optimization technique in `pow12_5` may have the opposite effect."
    end
    @noinline l_exp_y_log_x() = for x in xs; exp(1/2.4 * log(x)) end
    @noinline l_pow5_12() = for x in xs; Colors.pow5_12(x) end
    l_exp_y_log_x(); t_exp_y_log_x = @elapsed l_exp_y_log_x()
    l_pow5_12(); t_pow5_12 = @elapsed l_pow5_12()
    if t_pow5_12 > t_exp_y_log_x
        @warn "Optimization technique in `pow5_12` may have the opposite effect."
    end

    @test Colors.pow5_12(0.6) ≈ Colors.pow5_12(big"0.6") atol=1e-6
    @test Colors.pow5_12(0.6f0) ≈ Colors.pow5_12(big"0.6") atol=1e-6
    @test Colors.pow5_12(0.6N0f16) ≈ Colors.pow5_12(big"0.6") atol=1e-6
    @test Colors.pow12_5(0.6) ≈ Colors.pow12_5(big"0.6") atol=1e-6
    @test Colors.pow12_5(0.6N0f16) ≈ Colors.pow12_5(big"0.6") atol=1e-6

    # TODO: migrate to ColorTypes.jl
    @testset "atan360" begin
        for y in (-Inf, -0.0, +0.0, Inf, NaN), x in (-Inf, -0.0, +0.0, Inf, NaN)
            z0 = atand(y, x)
            z = signbit(z0) ? 360 + z0 : z0
            @test Colors.atan360(y, x) === z
            @test Colors.atan360(Float32(y), Float32(x)) === Float32(z)
        end
        for theta in 0.0:15.0:360.0
            yf64, xf64 = (sind(theta), cosd(theta)) .* 65.4321
            ybf, xbf = big(yf64), big(xf64)
            @test Colors.atan360(yf64, xf64) ≈ Colors.atan360(ybf, xbf) rtol=3eps(Float64)
            yf32, xf32 = Float32(yf64), Float32(xf64)
            ybf, xbf = big(yf32), big(xf32)
            @test Colors.atan360(yf32, xf32) ≈ Colors.atan360(ybf, xbf) rtol=3eps(Float32)
        end
    end

    @testset "polar_to_cartesian" begin
        @test Colors.polar_to_cartesian(50.0, NaN) === (NaN, NaN)
        @test Colors.polar_to_cartesian(0.0f0, Inf32) === (NaN32, NaN32)
        @test Float64.(Colors.polar_to_cartesian(big"2.0", big"120.0")) === (-1.0, sqrt(3.0))

        theta = -360.0:15.0:540.0
        yxbf = [sincos(deg2rad(big(t))) .* 64.5 for t in theta]
        xyf64 = Colors.polar_to_cartesian.(64.5, theta)
        xyf32 = Colors.polar_to_cartesian.(64.5f0, Float32.(theta))

        @test all(((a, b),) -> isapprox(a[1], b[2], atol=1e-14) &&
                               isapprox(a[2], b[1], atol=1e-14), zip(xyf64, yxbf))
        @test all(((a, b),) -> isapprox(a[1], b[2], atol=1e-6) &&
                               isapprox(a[2], b[1], atol=1e-6), zip(xyf32, yxbf))
    end

    @testset "hex" begin
        @test hex(RGB(1,0.5,0)) == "FF8000"
        @test hex(RGBA(1,0.5,0,0.25)) == "FF800040"
        @test hex(ARGB(1,0.5,0,0.25)) == "40FF8000"
        @test hex(HSV(30,1.0,1.0)) == "FF8000"

        @test hex(RGB(1,0.5,0), :AUTO) == "FF8000"
        @test hex(RGBA(1,0.5,0,0.25), :AUTO) == "FF800040"
        @test hex(ARGB(1,0.5,0,0.25), :AUTO) == "40FF8000"
        @test hex(HSV(30,1.0,1.0), :AUTO) == "FF8000"
        @test hex(RGB(1,0.5,0), :SomethingUnknown) == "FF8000"

        @test hex(RGB(1,0.5,0), :S) == "FF8000"
        @test hex(RGBA(1,0.5,0,0.25), :S) == "FF800040"
        @test hex(ARGB(1,0.5,0,0.25), :S) == "40FF8000"
        @test hex(RGB(1,0.533,0), :S) == "F80"
        @test hex(RGBA(1,0.533,0,0.267), :S) == "F804"
        @test hex(ARGB(1,0.533,0,0.267), :S) == "4F80"

        @test hex(RGB(1,0.5,0), :s) == "ff8000"
        @test hex(RGBA(1,0.5,0,0.25), :s) == "ff800040"
        @test hex(ARGB(1,0.5,0,0.25), :s) == "40ff8000"
        @test hex(RGB(1,0.533,0), :s) == "f80"
        @test hex(RGBA(1,0.533,0,0.267), :s) == "f804"
        @test hex(ARGB(1,0.533,0,0.267), :s) == "4f80"

        @test hex(RGB(1,0.5,0), :RGB) == "F80"
        @test hex(RGBA(1,0.5,0,0.25), :RGB) == "F80"
        @test hex(ARGB(1,0.5,0,0.25), :RGB) == "F80"
        @test hex(HSV(30,1.0,1.0), :RGB) == "F80"

        @test hex(RGB(1,0.5,0), :rgb) == "f80"
        @test hex(RGBA(1,0.5,0,0.25), :rgb) == "f80"
        @test hex(ARGB(1,0.5,0,0.25), :rgb) == "f80"
        @test hex(HSV(30,1.0,1.0), :rgb) == "f80"

        @test hex(RGB(1,0.5,0), :ARGB) == "FF80"
        @test hex(RGBA(1,0.5,0,0.25), :ARGB) == "4F80"
        @test hex(ARGB(1,0.5,0,0.25), :ARGB) == "4F80"
        @test hex(HSV(30,1.0,1.0), :ARGB) == "FF80"

        @test hex(RGB(1,0.5,0), :argb) == "ff80"
        @test hex(RGBA(1,0.5,0,0.25), :argb) == "4f80"
        @test hex(ARGB(1,0.5,0,0.25), :argb) == "4f80"
        @test hex(HSV(30,1.0,1.0), :argb) == "ff80"

        @test hex(RGB(1,0.5,0), :RGBA) == "F80F"
        @test hex(RGBA(1,0.5,0,0.25), :RGBA) == "F804"
        @test hex(ARGB(1,0.5,0,0.25), :RGBA) == "F804"
        @test hex(HSV(30,1.0,1.0), :RGBA) == "F80F"

        @test hex(RGB(1,0.5,0), :rgba) == "f80f"
        @test hex(RGBA(1,0.5,0,0.25), :rgba) == "f804"
        @test hex(ARGB(1,0.5,0,0.25), :rgba) == "f804"
        @test hex(HSV(30,1.0,1.0), :rgba) == "f80f"

        @test hex(RGB(1,0.5,0), :rrggbb) == "ff8000"
        @test hex(RGBA(1,0.5,0,0.25), :rrggbb) == "ff8000"
        @test hex(ARGB(1,0.5,0,0.25), :rrggbb) == "ff8000"
        @test hex(HSV(30,1.0,1.0), :rrggbb) == "ff8000"

        @test hex(RGB(1,0.5,0), :RRGGBB) == "FF8000"
        @test hex(RGBA(1,0.5,0,0.25), :RRGGBB) == "FF8000"
        @test hex(ARGB(1,0.5,0,0.25), :RRGGBB) == "FF8000"
        @test hex(HSV(30,1.0,1.0), :RRGGBB) == "FF8000"

        @test hex(RGB(1,0.5,0), :rrggbb) == "ff8000"
        @test hex(RGBA(1,0.5,0,0.25), :rrggbb) == "ff8000"
        @test hex(ARGB(1,0.5,0,0.25), :rrggbb) == "ff8000"
        @test hex(HSV(30,1.0,1.0), :rrggbb) == "ff8000"

        @test hex(RGB(1,0.5,0), :AARRGGBB) == "FFFF8000"
        @test hex(RGBA(1,0.5,0,0.25), :AARRGGBB) == "40FF8000"
        @test hex(ARGB(1,0.5,0,0.25), :AARRGGBB) == "40FF8000"
        @test hex(HSV(30,1.0,1.0), :AARRGGBB) == "FFFF8000"

        @test hex(RGB(1,0.5,0), :aarrggbb) == "ffff8000"
        @test hex(RGBA(1,0.5,0,0.25), :aarrggbb) == "40ff8000"
        @test hex(ARGB(1,0.5,0,0.25), :aarrggbb) == "40ff8000"
        @test hex(HSV(30,1.0,1.0), :aarrggbb) == "ffff8000"

        @test hex(RGB(1,0.5,0), :RRGGBBAA) == "FF8000FF"
        @test hex(RGBA(1,0.5,0,0.25), :RRGGBBAA) == "FF800040"
        @test hex(ARGB(1,0.5,0,0.25), :RRGGBBAA) == "FF800040"
        @test hex(HSV(30,1.0,1.0), :RRGGBBAA) == "FF8000FF"

        @test hex(RGB(1,0.5,0), :rrggbbaa) == "ff8000ff"
        @test hex(RGBA(1,0.5,0,0.25), :rrggbbaa) == "ff800040"
        @test hex(ARGB(1,0.5,0,0.25), :rrggbbaa) == "ff800040"
        @test hex(HSV(30,1.0,1.0), :rrggbbaa) == "ff8000ff"

        @test hex(Gray(0.5)) == "808080"
        @test hex(AGray(1.0, 0.5), :aarrggbb) == "80ffffff"
        @test hex(Gray{Bool}(1)) == "FFFFFF"

        # clamping
        @test hex(RGB(2.0,-1.0,0.5)) == "FF0080"
        @test hex(ARGB(2.0,-1.0,0.5,-10), :rrggbbaa) == "ff008000"
        @test hex(AHSV(30,1.0,1.0,-10), :RRGGBBAA) == "FF800000"
        @test hex(Gray(2.0)) == "FFFFFF"
        @test hex(AGray(1.0, -0.5), :rrggbbaa) == "ffffff00"
    end

    @testset "normalize_hue" begin
        function test_normalize_hue(h::T) where T<:AbstractFloat
            hn, he = normalize_hue(h), T(mod(h, BigFloat(360)))
            !signbit(hn) & (hn <= 360) || return false
            hn ≈ he || (hn == 360 && he == 0) || (hn == 0 && he == 360)
        end

        hues = -1500:1:1500
        @test all(h -> test_normalize_hue(Float64(h)), hues)
        @test all(h -> test_normalize_hue(prevfloat(Float64(h))), hues)
        @test all(h -> test_normalize_hue(nextfloat(Float64(h))), hues)

        @test all(h -> test_normalize_hue(Float32(h)), hues)
        @test all(h -> test_normalize_hue(prevfloat(Float32(h))), hues)
        @test all(h -> test_normalize_hue(nextfloat(Float32(h))), hues)

        @test all(h -> test_normalize_hue(Float16(h)), hues)
        @test all(h -> test_normalize_hue(prevfloat(Float16(h))), hues)
        @test all(h -> test_normalize_hue(nextfloat(Float16(h))), hues)

        @test normalize_hue(HSV{Float64}( 876.5, 0.4, 0.3)) === HSV{Float64}(156.5, 0.4, 0.3)
        @test normalize_hue(HSL{Float32}( 87.65, 0.4, 0.3)) === HSL{Float32}(87.65, 0.4, 0.3)
        @test normalize_hue(HSI{Float16}(-876.5, 0.4, 0.3)) === HSI{Float16}(203.5, 0.4, 0.3)

        @test normalize_hue(AHSV{Float32}(360.5, 0.4, 0.3, 0.2)) === AHSV{Float32}(0.5, 0.4, 0.3, 0.2)
        @test normalize_hue(HSLA{Float64}(360.5, 0.4, 0.3, 0.2)) === HSLA{Float64}(0.5, 0.4, 0.3, 0.2)

        @test normalize_hue(LCHab{Float64}(30, 40,  5.6e7)) === LCHab{Float64}(30, 40, 200)
        @test normalize_hue(LCHuv{Float32}(30, 40, -5.6e7)) === LCHuv{Float32}(30, 40, 160)

        @test normalize_hue(ALCHab{Float32}(30, 40, -0.5, 0.6)) === ALCHab{Float32}(30, 40, 359.5, 0.6)
        @test normalize_hue(LCHuvA{Float64}(30, 40, -0.5, 0.6)) === LCHuvA{Float64}(30, 40, 359.5, 0.6)
    end

    @testset "mean_hue" begin
        @test @inferred(mean_hue(0.0, 90.0)) === 45.0
        @test @inferred(mean_hue(90.0f0, 270.0f0)) === 180.0f0
        @test @inferred(mean_hue(90, 271)) === 0.5

        @test @inferred(mean_hue(HSV(50, 1, 1), HSV(100, 1, 0.5))) === 75.0
        @test @inferred(mean_hue(HSL(50, 1, 1), HSL(100, 0, 0.5))) === 50.0
        @test @inferred(mean_hue(Lab(50, 0, 10), Lab(50, 10, -10))) === 22.5f0
        @test @inferred(mean_hue(ALuv(50, 0, 0), ALuv(50, 0, 10))) === 90.0f0

        @test @inferred(mean_hue(LCHab(50, 60, 270), LCHab(90, 80, 80))) === 355.0f0
        @test @inferred(mean_hue(LCHuv(50, 60, 80), LCHuv(90, 0, 70))) === 80.0f0
        @test @inferred(mean_hue(LCHabA(90.0, 0, 70), LCHabA(50, 60, 80))) === 80.0

        @test_throws Exception mean_hue(RGB(0.1, 0.2, 0.3), RGB(0.4, 0.5, 0.6))
        @test_throws Exception mean_hue(LChab(10, 20, 30), LCHuv(10, 20, 30))
    end

    @testset "delta_h" begin
        @test @inferred(Colors.delta_h(LCHab(50, 60, 135.123f0), LCHab(90, 80, 135.123f0))) === 0.0f0
        @test @inferred(Colors.delta_h(LCHuv(50, 60, 80), LCHuv(90,  0, 70))) === 0.0f0
        @test @inferred(Colors.delta_h(LCHab(50, 60, 70), LCHab(90, 80, 190))) ≈ -120.0f0
        @test @inferred(Colors.delta_h(LCHab(90, 80, 70), LCHab(50, 60, 310.0))) ≈ 120.0
        @test_throws DomainError Colors.delta_h(LCHab(50, -1f-9, 70), LCHab(50, 60, 70))
        @test @inferred(Colors.delta_h(Lab(50, -10, 10), Lab(50, 10, -10))) ≈ 28.284271f0
        @test @inferred(Colors.delta_h(ALuv(50, 0, 0), ALuv(50, 0, 10))) === 0.0f0

        a, b = Lab(50, 4, 80), Lab(60, 5, 101)
        @test Colors.delta_h(a, b) ≈ -0.04444444f0 atol = 1f-4
        dhb = Colors.delta_h(Lab{BigFloat}(a), Lab{BigFloat}(b))
        @test Colors.delta_h(Lab{Float64}(a), Lab{Float64}(b)) ≈ dhb atol = 1e-13
    end

    @testset "delta_c" begin
        a, b = Lab(50, 40, -31), Lab(90, -30, 41)
        @test @inferred(Colors.delta_c(a, b)) ≈ -0.19721946f0 atol = 1f-7
        @test @inferred(Colors.delta_c(LCHab(a), LCHab(b))) ≈ -0.19721946f0 atol = 1f-6
    end

    # test utility function weighted_color_mean
    parametric2 = [GrayA,AGray32,AGray]
    parametric3 = ColorTypes.parametric3
    parametric4A = setdiff(subtypes(ColorAlpha),[GrayA])
    parametricA4 = setdiff(subtypes(AlphaColor),[AGray32,AGray])
    colortypes = vcat(parametric3,parametric4A,parametricA4)
    colorElementTypes = [Float16,Float32,Float64,BigFloat,N0f8,N6f10,N4f12,N2f14,N0f16]

    @testset "weighted_color_mean" begin
        for T in colorElementTypes
            c1 = Gray{T}(1)
            c2 = Gray{T}(0)
            @test weighted_color_mean(0.5,c1,c2) == Gray{T}(0.5)
        end
        gray_b1 = Gray{Bool}(1)
        gray_b0 = Gray{Bool}(0)
        @test @inferred(weighted_color_mean(0, gray_b1, gray_b0)) === gray_b0
        @test @inferred(weighted_color_mean(1, gray_b1, gray_b0)) === gray_b1
        @test @inferred(weighted_color_mean(0.5, gray_b1, gray_b1)) === gray_b1
        @test_throws InexactError weighted_color_mean(0.5, gray_b1, gray_b0)
        @test_throws DomainError weighted_color_mean(-1, gray_b1, gray_b0)

        for C in parametric2
            for T in colorElementTypes
                c1 = C(T(0), T(1))
                c2 = C(T(1), T(0))
                cx = C(T(0.5), T(0.5))
                @test @inferred(weighted_color_mean(0.5, c1, c2)) == cx
            end
        end

        for C in colortypes
            for T in (Float16,Float32,Float64,BigFloat)
                if C<:Color
                    c1 = C(T(1), T(1), T(0))
                    c2 = C(T(0), T(1), T(1))
                    cx = C(T(0.5), T(1.0), T(0.5))
                    @test @inferred(weighted_color_mean(0.5, c1, c2)) == cx
                else
                    c1 = C(T(1), T(1), T(0), T(0))
                    c2 = C(T(0), T(1), T(0), T(1))
                    cx = C(T(0.5), T(1.0), T(0.0), T(0.5))
                    @test @inferred(weighted_color_mean(0.5, c1, c2)) == cx
                end
            end
        end

        @test_throws DomainError weighted_color_mean(-0.01, RGB(1,1,0), RGB(0,1,1))
        @test_throws DomainError weighted_color_mean( 1.01, RGB(1,1,0), RGB(0,1,1))

        # promotion
        @test weighted_color_mean(0.5, RGB{N0f8}(1,1,0), RGB{Float32}(0,1,1)) === RGB{Float32}(0.5,1,0.5)
        @test_throws ArgumentError weighted_color_mean(0.5, RGB24(1,1,0), RGB{Float32}(0,1,1))

        @test weighted_color_mean(0.5, HSV(-360.0,1,0), HSV(180.0,0,1)) === HSV{Float64}(270,0.5,0.5)
        alchab1 = ALCHab{Float32}(0,100,-360,1)
        alchab2 = ALCHab{Float32}(100,0,-180,0)
        @test weighted_color_mean(0.5, alchab1, alchab2) === ALCHab{Float32}(50,50,90,0.5)
        lchuva1 = LCHuvA{Float16}(0,100, 90,1)
        lchuva2 = LCHuvA{Float16}(100,0,810,0)
        @test weighted_color_mean(0.5, lchuva1, lchuva2) === LCHuvA{Float16}(50,50,90,0.5)

        @testset "with collections" begin
            colors = [RGBA(1, 0, 0, 1) RGBA(0, 0, 1, 1)
                      RGBA(0, 1, 0, 1) RGBA(0, 0, 0, 0)]
            c = RGBA{N0f8}(0.1, 0.2, 0.3, 0.6)
            @test weighted_color_mean([0.1 0.3; 0.2 0.4],   colors) === c # matrix
            @test weighted_color_mean((0.1, 0.2, 0.3, 0.4), colors) === c # tuple
            @test weighted_color_mean(0.1:0.1:0.4,          colors) === c # range

            @test weighted_color_mean((0, 1), (Gray(1), Gray(0.0))) === Gray{N0f8}(0)
            @test weighted_color_mean((1, 0), (Gray(1.0), Gray(0))) === Gray{Float64}(1)

            @test_throws DimensionMismatch weighted_color_mean((0.1, 0.2, 0.3), colors)
            @test_throws Exception weighted_color_mean((0.1, 0.2), (RGB(1, 1, 1), HSV(0, 1, 1)))
        end
    end

    @testset "range" begin
        # `range` uses `weighted_color_mean` which is extensively tested.
        # Therefore it suffices to test the function using gray colors.
        for T in colorElementTypes
            c1 = Gray(T(1))
            c2 = Gray(T(0))
            linc1c2 = range(c1,stop=c2,length=43)
            @test length(linc1c2) == 43
            @test linc1c2[1] == c1
            @test linc1c2[end] == c2
            @test linc1c2[22] == Gray(T(0.5))
            @test typeof(linc1c2) == Array{Gray{T},1}
            if VERSION >= v"1.1"
                @test range(c1,c2,length=43) == range(c1,stop=c2,length=43)
            end
        end
    end
end


text:
using Match

# Functions for performing PPU cycles.

function ppuread(console::Console, address::UInt16)::UInt8
  ppu = console.ppu
  address = address & 0x3FFF
  if address < 0x2000
    read(console.mapper, console.cartridge, address)
  elseif address < 0x3F00
    mode = console.cartridge.mirror
    @inbounds ppu.nameTableData[(mirroraddress(mode, address) & 0x07FF) + 0x01]
  else
    readpalette(ppu, address & 0x1F)
  end
end

function ppuwrite!(console::Console, address::UInt16, val::UInt8)
  ppu = console.ppu
  address = address & 0x3FFF
  if address < 0x2000
    write!(console.mapper, console.cartridge, address, val)
  elseif address < 0x3F00
    mode = console.cartridge.mirror
    @inbounds ppu.nameTableData[(mirroraddress(mode, address) & 0x07FF) + 0x01] = val
  else
    writepalette!(ppu, address & 0x1F, val)
  end
end

function ppureset!(ppu::PPU)
  ppu.cycle = 340
  ppu.scanline = 240
  ppu.frame = 0
  writecontrol!(ppu, 0x00)
  writemask!(ppu, 0x00)
  writeoamaddress!(ppu, 0x00)
end

function readpalette(ppu::PPU, address::UInt16)::UInt8
  if address >= 16 && (address & 3) == 0
    address -= 0x10
  end
  @inbounds return ppu.paletteData[address + 1]
end

function writepalette!(ppu::PPU, address::UInt16, val::UInt8)
  if address >= 16 && (address & 3) == 0
    address -= 0x10
  end
  @inbounds ppu.paletteData[address + 1] = val
end

function ppureadregister(console::Console, address::UInt16)::UInt8
  ppu = console.ppu
  @match address begin
    0x2002 => readstatus!(ppu)
    0x2004 => readoamdata(ppu)
    0x2007 => readdata!(console)
    _ => 0
  end
end

function ppuwriteregister!(console::Console, address::UInt16, val::UInt8)
  ppu = console.ppu
  ppu.reg = val
  @match address begin
    0x2000 => begin
      writecontrol!(ppu, val)
    end
    0x2001 => begin
      writemask!(ppu, val)
    end
    0x2003 => begin
      writeoamaddress!(ppu, val)
    end
    0x2004 => begin
      writeoamdata!(ppu, val)
    end
    0x2005 => begin
      writescroll!(ppu, val)
    end
    0x2006 => begin
      writeaddress!(ppu, val)
    end
    0x2007 => begin
      writedata!(console, val)
    end
    0x4014 => begin
      writedma!(console, val)
    end
  end
end

# $2000: PPUCTRL
function writecontrol!(ppu::PPU, val::UInt8)
  ppu.flagNameTable = (val >> 0) & 3
  ppu.flagIncrement = (val >> 2) & 1
  ppu.flagSpriteTable = (val >> 3) & 1
  ppu.flagBackgroundTable = (val >> 4) & 1
  ppu.flagSpriteSize = (val >> 5) & 1
  ppu.flagMasterSlave = (val >> 6) & 1
  ppu.nmiOutput = ((val >> 7) & 1) == 1
  nmichange!(ppu)
  # t: ....BA.. ........ = d: ......BA
  ppu.t = (ppu.t & 0xF3FF) | ((val & 0x0003) << 10)
end

# $2001: PPUMASK
function writemask!(ppu::PPU, val::UInt8)
  ppu.flagGrayscale = (val >> 0) & 1
  ppu.flagShowLeftBackground = (val >> 1) & 1
  ppu.flagShowLeftSprites = (val >> 2) & 1
  ppu.flagShowBackground = (val >> 3) & 1
  ppu.flagShowSprites = (val >> 4) & 1
  ppu.flagRedTint = (val >> 5) & 1
  ppu.flagGreenTint = (val >> 6) & 1
  ppu.flagBlueTint = (val >> 7) & 1
end

# $2002: PPUSTATUS
function readstatus!(ppu::PPU)::UInt8
  result = ppu.reg & 0x1F
  result |= ppu.flagSpriteOverflow << 5
  result |= ppu.flagSpriteZeroHit << 6
  if ppu.nmiOccurred
    result |= 0x01 << 7
  end
  ppu.nmiOccurred = false
  nmichange!(ppu)
  # w:                   = 0
  ppu.w = 0
  result
end

# $2003: OAMADDR
function writeoamaddress!(ppu::PPU, val::UInt8)
  ppu.oamAddress = val
end

# $2004: OAMDATA (read)
function readoamdata(ppu::PPU)::UInt8
  @inbounds ppu.oamData[ppu.oamAddress + 1]
end

# $2004: OAMDATA (write)
function writeoamdata!(ppu::PPU, val::UInt8)
  @inbounds ppu.oamData[ppu.oamAddress + 1] = val
  ppu.oamAddress += 0x01
end

# $2005: PPUSCROLL
function writescroll!(ppu::PPU, val::UInt8)
  if ppu.w == 0
    # t: ........ ...HGFED = d: HGFED...
    # x:               CBA = d: .....CBA
    # w:                   = 1
    ppu.t = (ppu.t & 0xFFE0) | (UInt16(val) >> 3)
    ppu.x = val & 0x07
    ppu.w = 1
  else
    # t: .CBA..HG FED..... = d: HGFEDCBA
    # w:                   = 0
    ppu.t = (ppu.t & 0x8FFF) | ((val & 0x0007) << 12)
    ppu.t = (ppu.t & 0xFC1F) | ((val & 0x00F8) << 2)
    ppu.w = 0
  end
end

# $2006: PPUADDR
function writeaddress!(ppu::PPU, val::UInt8)
  if ppu.w == 0
    # t: ..FEDCBA ........ = d: ..FEDCBA
    # t: .X...... ........ = 0
    # w:                   = 1
    ppu.t = (ppu.t & 0x80FF) | ((val & 0x003F) << 8)
    ppu.w = 1
  else
    # t: ........ HGFEDCBA = d: HGFEDCBA
    # v                    = t
    # w:                   = 0
    ppu.t = (ppu.t & 0xFF00) | val
    ppu.v = ppu.t
    ppu.w = 0
  end
end

# $2007: PPUDATA (read)
function readdata!(console::Console)::UInt8
  ppu = console.ppu
  val = ppuread(console, ppu.v)
  # emulate buffered reads
  if ppu.v & 0x3FFF < 0x3F00
    buffered = ppu.bufferedData
    ppu.bufferedData = val
    val = buffered
  else
    ppu.bufferedData = ppuread(console, ppu.v - 0x1000)
  end
  # increment address
  if ppu.flagIncrement == 0
    ppu.v += 0x01
  else
    ppu.v += 0x20
  end
  return val
end

# $2007: PPUDATA (write)
function writedata!(console::Console, val::UInt8)
  ppu = console.ppu
  ppuwrite!(console, ppu.v, val)
  if ppu.flagIncrement == 0
    ppu.v += 0x01
  else
    ppu.v += 0x20
  end
end

# $4014: OAMDMA
function writedma!(console::Console, val::UInt8)
  ppu = console.ppu
  cpu = console.cpu
  address = UInt16(val) << 8
  for i in 0:255
    @inbounds ppu.oamData[ppu.oamAddress + 1] = cpuread(console, address)
    ppu.oamAddress += 0x01
    address += 0x0001
  end
  cpu.stall += Int32(513)
  if cpu.cycles & 1 == 1
    cpu.stall += 0x01
  end
end

# NTSC Timing Helper Functions

function incrementx!(ppu::PPU)
  # increment hori(v)
  # if coarse X == 31
  if (ppu.v & 0x001F) == 31
    # coarse X = 0
    ppu.v &= 0xFFE0
    # switch horizontal nametable
    ppu.v ⊻= 0x0400
  else
    # increment coarse X
    ppu.v += 0x01
  end
end

function incrementy!(ppu::PPU)
  # increment vert(v)
  # if fine Y < 7
  if (ppu.v&0x7000) != 0x7000
    # increment fine Y
    ppu.v += 0x1000
  else
    # fine Y = 0
    ppu.v &= 0x8FFF
    # let y = coarse Y
    y::UInt16 = (ppu.v & 0x03E0) >> 5
    if y == 29
      # coarse Y = 0
      y = 0x0000
      # switch vertical nametable
      ppu.v ⊻= 0x0800
    elseif y == 31
      # coarse Y = 0, nametable not switched
      y = 0x0000
    else
      # increment coarse Y
      y += 0x0001
    end
    # put coarse Y back into v
    ppu.v = (ppu.v & 0xFC1F) | (y << 5)
  end
end

function copyx!(ppu::PPU)
  # hori(v) = hori(t)
  # v: .....F.. ...EDCBA = t: .....F.. ...EDCBA
  ppu.v = (ppu.v & 0xFBE0) | (ppu.t & 0x041F)
end

function copyy!(ppu::PPU)
  # vert(v) = vert(t)
  # v: .IHGF.ED CBA..... = t: .IHGF.ED CBA.....
  ppu.v = (ppu.v & 0x841F) | (ppu.t & 0x7BE0)
end

function nmichange!(ppu::PPU)
  nmi = ppu.nmiOutput && ppu.nmiOccurred
  if nmi && !ppu.nmiPrevious
    # TODO: this fixes some games but the delay shouldn't have to be so
    # long, so the timings are off somewhere
    ppu.nmiDelay = 15
  end
  ppu.nmiPrevious = nmi
end

function setverticalblank!(ppu::PPU)
  temp = ppu.front
  ppu.front = ppu.back
  ppu.back = temp
  ppu.nmiOccurred = true
  nmichange!(ppu)
end

function clearverticalblank!(ppu::PPU)
  ppu.nmiOccurred = false
  nmichange!(ppu)
end

function fetchnametablebyte!(console::Console, ppu::PPU)
  address = 0x2000 | (ppu.v & 0x0FFF)
  ppu.nameTableByte = ppuread(console, address)
end

function fetchattributetablebyte!(console::Console, ppu::PPU)
  address = 0x23C0 | (ppu.v & 0x0C00) | ((ppu.v >> 4) & 0x38) | ((ppu.v >> 2) & 0x07)
  shift = ((ppu.v >> 4) & 4) | (ppu.v & 2)
  ppu.attributeTableByte = ((ppuread(console, address) >> shift) & 0x03) << 2
end

function fetchlowtilebyte!(console::Console, ppu::PPU)
  fineY = (ppu.v >> 12) & 0x07
  table = ppu.flagBackgroundTable
  tile = ppu.nameTableByte
  address = 0x1000 * table + tile * 0x0010 + fineY
  ppu.lowTileByte = ppuread(console, address)
end

function fetchhightilebyte!(console::Console, ppu::PPU)
  fineY = (ppu.v >> 0x0C) & 0x07
  table = ppu.flagBackgroundTable
  tile = ppu.nameTableByte
  address = 0x1000 * table + tile * 0x0010 + fineY
  ppu.highTileByte = ppuread(console, address + 0x08)
end

function storetiledata!(ppu::PPU)
  data = 0x00000000
  for i = 1:8
    a = ppu.attributeTableByte
    p1 = (ppu.lowTileByte & 0x80) >> 7
    p2 = (ppu.highTileByte & 0x80) >> 6
    ppu.lowTileByte <<= 1
    ppu.highTileByte <<= 1
    data <<= 4
    data |= a | p1 | p2
  end
  ppu.tileData |= data
end

function fetchtiledata(ppu::PPU)::UInt32
  ppu.tileData >> 32
end

function backgroundpixel(ppu::PPU)::UInt8
  if ppu.flagShowBackground == 0x00
    0x00
  else
    data = fetchtiledata(ppu) >> ((0x07 - ppu.x) << 2)
    data & 0x0F
  end
end

function spritepixel(ppu::PPU)::Tuple{UInt8, UInt8}
  if ppu.flagShowSprites != 0
    for i = 0:(ppu.spriteCount-1)
      @inbounds offset = (ppu.cycle - Int32(1)) - Int32(ppu.spritePositions[i + 1])
      if offset < 0 || offset > 7
        continue
      end
      offset = Int32(7) - offset
      @inbounds color = (ppu.spritePatterns[i + 1] >> (offset * 4)) & 0x0F
      if color & 3 == 0
        continue
      end
      return (UInt8(i), UInt8(color))
    end
  end
  (0x00, 0x00)
end

function renderpixel!(ppu::PPU, x::Int32, y::Int32)
  background = backgroundpixel(ppu)
  i, sprite = spritepixel(ppu)
  ppux = ppu.x
  if ppux < 8 && ppu.flagShowLeftBackground == 0
    background = 0x00
  end
  if ppux < 8 && ppu.flagShowLeftSprites == 0
    sprite = 0x00
  end
  b = (background & 3) != 0
  s = (sprite & 3) != 0
  color = 0x00
  if !b && !s
    color = 0x00
  elseif !b && s
    color = sprite | 0x10
  elseif b && !s
    color = background
  else
    @inbounds if ppu.spriteIndexes[i + 1] == 0 && ppux < 255
      ppu.flagSpriteZeroHit = 1
    end
    @inbounds if ppu.spritePriorities[i + 1] == 0
      color = sprite | 0x10
    else
      color = background
    end
  end
  pixel = readpalette(ppu, UInt16(color)) & 0x3F
  @inbounds ppu.back[y + 1, x] = pixel
end

function fetchspritepattern(console::Console, i::Int32, row::Int32)::UInt32
  ppu = console.ppu
  @inbounds tile = ppu.oamData[4 * i + 2]
  @inbounds attributes = ppu.oamData[4 * i + 3]
  address = 0x0000
  if ppu.flagSpriteSize == 0
    if (attributes & 0x80) == 0x80
      row = Int32(7) - row
    end
    table = ppu.flagSpriteTable
    address = 0x1000 * table + 0x0010 * tile + UInt16(row)
  else
    if (attributes & 0x80) == 0x80
      row = Int32(15) - row
    end
    table = tile & 0x01
    tile &= 0xFE
    if row > 7
      tile += 0x01
      row -= Int32(8)
    end
    address = 0x1000 * table + 0x0010 * tile + UInt16(row)
  end
  a = (attributes & 0x03) << 0x02
  lowTileByte = ppuread(console, address)
  highTileByte = ppuread(console, address + 0x08)
  data = 0x00000000
  for i = 1:8
    p1 = 0x00
    p2 = 0x00
    if (attributes & 0x40) == 0x40
      p1 = (lowTileByte & 0x01) << 0x00
      p2 = (highTileByte & 0x01) << 0x01
      lowTileByte >>= 1
      highTileByte >>= 1
    else
      p1 = (lowTileByte & 0x80) >> 0x07
      p2 = (highTileByte & 0x80) >> 0x06
      lowTileByte <<= 1
      highTileByte <<= 1
    end
    data <<= 4
    data |= a | p1 | p2
  end
  data
end

function evaluatesprites!(console::Console)
  ppu = console.ppu
  h::Int32 = ppu.flagSpriteSize == 0 ? 8 : 16
  count::Int32 = 0
  for i = Int32(0):Int32(63)
    @inbounds y = ppu.oamData[4 * i + 1]
    @inbounds a = ppu.oamData[4 * i + 3]
    @inbounds x = ppu.oamData[4 * i + 4]
    row = ppu.scanline - y
    if row < 0 || row >= h
      continue
    end
    if count < 8
      @inbounds ppu.spritePatterns[count + 1] = fetchspritepattern(console, i, row)
      @inbounds ppu.spritePositions[count + 1] = x
      @inbounds ppu.spritePriorities[count + 1] = (a >> 5) & 1
      @inbounds ppu.spriteIndexes[count + 1] = i
    end
    count += Int32(1)
  end
  if count > 8
    count = Int32(8)
    ppu.flagSpriteOverflow = 1
  end
  ppu.spriteCount = count
end

# tick updates cycle, scanline and frame counters
function tick!(console::Console, ppu::PPU)
  if ppu.nmiDelay > 0
    ppu.nmiDelay -= 0x01
    if ppu.nmiDelay == 0 && ppu.nmiOutput && ppu.nmiOccurred
      triggernmi!(console.cpu)
    end
  end

  if ppu.flagShowBackground != 0 || ppu.flagShowSprites != 0
    if ppu.f == 1 && ppu.scanline == 261 && ppu.cycle == 339
      ppu.cycle = 0
      ppu.scanline = 0
      ppu.frame += 1
      ppu.f ⊻= 0x01
      return
    end
  end
  ppu.cycle += Int32(1)
  if ppu.cycle > 340
    ppu.cycle = 0
    ppu.scanline += Int32(1)
    if ppu.scanline > 261
      ppu.scanline = 0
      ppu.frame += 1
      ppu.f ⊻= 1
    end
  end
end

# Step executes a single PPU cycle
function ppustep!(console::Console, ppu::PPU)
  tick!(console, ppu)

  renderingEnabled = ppu.flagShowBackground != 0x00 || ppu.flagShowSprites != 0x00
  scanline = ppu.scanline
  preLine = (scanline == Int32(261))
  visibleLine = (scanline < Int32(240))
  # postLine := scanline == 240
  renderLine = (preLine || visibleLine)
  cycle = ppu.cycle
  preFetchCycle = (cycle >= Int32(321) && cycle <= Int32(336))
  visibleCycle = (cycle >= Int32(1) && cycle <= Int32(256))
  fetchCycle = (preFetchCycle || visibleCycle)

  # background logic
  if renderingEnabled
    if visibleLine && visibleCycle
      renderpixel!(ppu, cycle, scanline)
    end
    if renderLine && fetchCycle
      tdb = ppu.tileData
      ppu.tileData <<= 4
      b = ppu.cycle & 7
      @match b begin
        1 => begin
          fetchnametablebyte!(console, ppu)
        end
        3 => begin
          fetchattributetablebyte!(console, ppu)
        end
        5 => begin
          fetchlowtilebyte!(console, ppu)
        end
        7 => begin
          fetchhightilebyte!(console, ppu)
        end
        0 => begin
          storetiledata!(ppu)
        end
      end
    end
    if preLine && ppu.cycle >= 280 && ppu.cycle <= 304
      copyy!(ppu)
    end
    if renderLine
      if fetchCycle && (ppu.cycle & 7 == 0)
        incrementx!(ppu)
      end
      if ppu.cycle == 256
        incrementy!(ppu)
      end
      if ppu.cycle == 257
        copyx!(ppu)
      end
    end
  end

  # sprite logic
  if renderingEnabled
    if ppu.cycle == 257
      if visibleLine
        evaluatesprites!(console)
      else
        ppu.spriteCount = 0
      end
    end
  end

  # vblank logic
  if ppu.scanline == 241 && ppu.cycle == 1
    setverticalblank!(ppu)
  end
  if preLine && ppu.cycle == 1
    clearverticalblank!(ppu)
    ppu.flagSpriteZeroHit = 0
    ppu.flagSpriteOverflow = 0
  end
  0x00
end

ext: jl


text:
function mapper_parserecordacomp(x::AbstractString)
  line = chomp(x)
  separator = ","
  isempty(line) && throw(MapperException("line is empty"))
  s = split(line, separator)
  length(s) < 6 && throw(MapperException("malformed line has less than six elements when split by ,"))

  departuretime = UDFDepartureTime(s[5])
  totalflighttime = UDFTotalFlightTime(s[6])
  arrivaltime = departuretime + Dates.Minute(totalflighttime)
  UDFRAcomp(
    UDFPassengerId(s[1]),
    UDFFlightId(s[2]),
    UDFAirportCode(s[3]),
    UDFAirportCode(s[4]),
    departuretime,
    arrivaltime,
    totalflighttime
  )
end

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
#
# CART-related types and functions
#

struct Cart <: Wrapper
    pyo::PyObject

    function Cart(x::DataSet, y::Vector; threshold=nothing, threshold_type=">",
                         include=nothing, exclude=nothing, kwargs...)

        pandasDF = pandas_dataframe(x; include=include, exclude=exclude)

        cart = rhodium.Cart(pandasDF, y; threshold=threshold, threshold_type=threshold_type,
                            kwargs...)
        new(cart)
    end
end

function show_tree(c::Cart; kwargs...)
    fig = c.pyo.show_tree(kwargs...)
    return fig
end

function print_tree(c::Cart; coi=nothing, all=true, kwargs...)
    c.pyo.print_tree(;coi=coi, all=all, kwargs...)
    nothing
end

function save(c::Cart, file, format="png"; kwargs...)
    c.pyo[:save](file, format; kwargs...)
end

function save_pdf(c::Cart, file; feature_names=nothing, kwargs...)
    save(c, file, "pdf"; feature_names=feature_names, kwargs...)
end

function save_png(c::Cart, file; feature_names=nothing, kwargs...)
    save(c, file, "png"; feature_names=feature_names, kwargs...)
end


text:
using DelimitedFiles

"""
    modify the values of vp,vs,rho in Crust1.0 with the upper layer's value, so that we won't have smoothing error.
"""
function modify_sediments(basedir::String, tag::String)
    filename = joinpath(basedir, "crust1." * tag)
    data = readdlm(filename)
    for i in 1:size(data)[1]
        if (data[i,3] == 0. ) && (data[i,4] == 0.) && (data[i,5] == 0.)
            continue
        elseif (data[i,3] != 0. ) && (data[i,4] == 0.) && (data[i,5] == 0.)
            data[i,4] = data[i,3]
            data[i,5] = data[i,3]
        elseif (data[i,3] != 0. ) && (data[i,4] != 0.) && (data[i,5] == 0.)
            data[i,5] = data[i,4]
        elseif (data[i,3] != 0. ) && (data[i,4] != 0.) && (data[i,5] != 0.)
            continue
        else
            @error "abnormal 0 points"
        end
    end

    filename_old = filename * ".old"
    run(`mv $(filename) $(filename_old)`)
    writedlm(filename, data)
end

function modify_sediments_all(basedir::String)
    tags = ["vp","vs","rho"]
    for tag in tags
        modify_sediments(basedir, tag)
    end
end

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts


text:
using PyCall
# using Conda
# Conda.add("ply")
# Conda.add("yaml")
runpy = pyimport("runpy")
push!(pyimport("sys")."path", pwd())
mymodule = pyimport("itf1788")
runpy.run_module("itf1788", run_name="__main__")


text:
using TimeSeriesClustering
using Test
using Random
using JLD2

@load normpath(joinpath(dirname(@__FILE__),"reference_generation","run_clust.jld2")) reference_results
data = "CEP_GER1"
ts_input_data = load_timeseries_data(Symbol(data))

include("test_utils.jl")

@testset "ClustResultAll" begin
    method = "hierarchical"
    repr = "medoid"
    Random.seed!(1111)
    ref_all = reference_results["$data-$method-$repr-ClustResultAll"]
    t_all = run_clust(ts_input_data;method=method,representation=repr,n_clust=5,n_init=10,get_all_clust_results=true)
    test_ClustResult(t_all,ref_all)
    for i = 1:length(t_all.centers_all)
        @test all(t_all.centers_all[i] .≈ ref_all.centers_all[i])
        @test all(t_all.weights_all[i] .≈ ref_all.weights_all[i])
        @test all(t_all.clustids_all[i] .≈ ref_all.clustids_all[i])
        @test all(t_all.cost_all[i] .≈ ref_all.cost_all[i])
        @test all(t_all.iter_all[i] .≈ ref_all.iter_all[i])
    end
end


text:
export @evalpoly_derivative, @goertzel_deriv, @horner_deriv

"""
    evalpoly(::Type{T}, degrees, coefficients, var::Union{Symbol, Expr})

Evaluate the polynomial defined by the degrees and coefficients.
"""
function evalpoly(::Type{T}, degrees::AbstractVector, coefficients::AbstractVector, var) where T
    normalized_coeffs = normalized_coefficients(T, degrees, coefficients)

    if length(normalized_coeffs) == 1
        return normalized_coeffs[1]
    elseif length(normalized_coeffs) == 2
        a, b = normalized_coeffs
        :(muladd($b, $var, $a))
    else
        :(@evalpoly($var, $(normalized_coeffs...)))
    end

end

"""
    evalpoly_derivative!(exprs, ::Type{T}, degrees, coefficients, var::Union{Symbol, Expr})

Evaluate the polynomial and its derivative defined by the degrees and coefficients.
"""
function evalpoly_derivative!(exprs, ::Type{T}, degrees::AbstractVector, coefficients::AbstractVector, var) where T
    normalized_coeffs = normalized_coefficients(T, degrees, coefficients)

    if length(normalized_coeffs) == 1
        if isa(normalized_coeffs[1], Symbol)
            val = normalized_coeffs[1]
        else
            @gensym val
            push!(exprs, :($val = $(normalized_coeffs[1])))
        end
        # return zero directly
        dval = :(zero($T))
    elseif length(normalized_coeffs) == 2
        @gensym val
        a, b = normalized_coeffs
        push!(exprs, :($val = muladd($b, $var, $a)))
        if isa(b, Symbol)
            dval = b
        else
            @gensym dval
            push!(exprs, :($dval = $b))
        end
    else
        @gensym dval val
        push!(exprs, :(($val, $dval) = @evalpoly_derivative($var, $(normalized_coeffs...))))
    end

    return val, dval
end

"""
    @evalpoly_derivative(z, c...)

Evaluate the polynomial ``\\sum_k c[k] z^{k-1}`` for the coefficients `c[1]`, `c[2]`, ... and its derivative;
that is, the coefficients are given in ascending order by power of `z`.  This macro expands
to efficient inline code that uses either Horner's method or, for complex `z`, a more
efficient Goertzel-like algorithm.

```jldoctest
julia> @evalpoly(3, 1, 0, 1)
(10, 6)

julia> @evalpoly(2, 1, 0, 1)
(5, 4)

julia> @evalpoly(2, 1, 1, 1)
(7, 5)
```
"""
macro evalpoly_derivative(z, p...)
    R = Expr(:macrocall, Symbol("@horner_deriv"), (), :tt, map(esc, p)...)
    C = Expr(:macrocall, Symbol("@goertzel_deriv"), (), :tt, map(esc, p)...)
    :(let tt = $(esc(z))
          isa(tt, Complex) ? $C : $R
      end)
end

macro horner_deriv(x, p...)
    escaped_x = :($(esc(x)))
    escaped_p = map(esc, p)
    horner_deriv_impl(escaped_x, escaped_p)
end

function horner_deriv_impl(x, p)
    if length(p) == 1
        return :($(p[1]), zero($(p[1])))
    end

    exprs = Expr[]
    push!(exprs, :(val = $(p[end])))
    push!(exprs, :(dval = val))
    for k = length(p)-1:-1:1
        if k < length(p) - 1
            push!(exprs, :(dval = muladd(dval, $x, val)))
        end
        push!(exprs, :(val = muladd(val, $x, $(p[k]))))
    end

    Expr(:block,
        exprs...,
        :(val, dval))
end

macro goertzel_deriv(z, p...)
    escaped_z = :($(esc(z)))
    escaped_p = map(esc, p)
    goertzel_deriv_impl(escaped_z, escaped_p)
end

function goertzel_deriv_impl(z, p)
    if length(p) == 1
        return :($(p[end]), zero($(p[end])))
    end

    if length(p) == 2
        return :(
            muladd($(p[end]), $z, $(p[end-1])),
            $(p[end])
        )
    end

    exprs = Expr[]
    goertzel_setup!(exprs, z)
    goertzel_deriv_main!(exprs, z, p)

    Expr(:block, exprs...)
end

function goertzel_setup!(exprs, z)
    push!(exprs,
        :(x = real($z)),
        :(y = imag($z)),
        :(r = x + x),
        :(s = -muladd(x, x, y*y)))
end

function goertzel_deriv_main!(exprs, z, p)
    @assert length(p) > 2

    a = :($(p[end]))
    b = :($(p[end-1]))
    local c, d
    # d = b
    for i = length(p)-2:-1:1
        ei = Symbol("e", i)
        push!(exprs, :($ei = $a)) #e_i
        a = :(muladd(r, $ei, $b))
        b = :(muladd(s, $ei, $(p[i])))

        if i == length(p) - 2
            c = ei
        elseif i == length(p) - 3
            d = ei
        else
            fi = Symbol("f", i)
            push!(exprs, :($fi = $c)) #e_i
            c = :(muladd(r, $fi, $d))
            d = :(muladd(s, $fi, $ei))
        end
    end
    r1 = :r1 # r1
    push!(exprs, :($r1 = $a))
    r2 = b
    push!(exprs, :(val = muladd($r1, $z, $r2)))

    # we handle the cases that length(p) < 3 before this routines is called
    if length(p) == 3
        lm = :($(p[end]))
        push!(exprs, :(t = y + y))
        dval = quote
            if isa(r1, Complex)
                complex(real(r1) - t * imag($lm), imag(r1) + t * real($lm))
            else
                complex(r1 - t * imag($lm), t * real($lm))
            end
        end
        push!(exprs, :((val, $dval)))
    else
        dr1 = :dr1
        push!(exprs, :($dr1 = $c))
        dr0 = d
        push!(exprs, :(dr0 = muladd($dr1, $z, $dr0)))

        push!(exprs, :(t = y + y))
        dval = quote
            if isa(r1, Complex)
                complex(real(r1) - t * imag(dr0), imag(r1) + t * real(dr0))
            else
                complex(r1 - t * imag(dr0), t * real(dr0))
            end
        end
        push!(exprs, :((val, $dval)))
    end
end


text:
import Base: size, convert

import SoftConfidenceWeighted: init, fit!, predict, SCW1, SCW2

function split_dataset(X, y, training_ratio=0.8)
    assert(0.0 <= training_ratio <= 1.0)

    split_point = convert(Int64, size(X, 2)*training_ratio)
    training = X[:, 1:split_point-1], y[1:split_point-1]
    test = X[:, split_point:end], y[split_point:end]
    return training, test
end

function calc_accuracy(y_pred, y_true)
    n_correct = 0
    for (a, b) in zip(y_pred, y_true)
        if a == b
            n_correct += 1
        end
    end

    return n_correct / length(y_pred)
end

function test_batch(X, y, type_; training_ratio = 0.8, C = 1.0, ETA = 1.0)
    model = init(C = C, ETA = ETA, type_ = type_)

    training, test = split_dataset(X, y, training_ratio)

    X, labels = training
    model = fit!(model, X, labels)

    X, y_true = test
    y_pred = predict(model, X)

    accuracy = calc_accuracy(y_pred, y_true)
    assert(accuracy == 1.0)

    println("BATCH")
    println("\ttype: $type_")
    println("\taccuracy: $accuracy")
    println("")
end

function test_online(X, y, type_; training_ratio=0.8, C=1.0, ETA=1.0)
    model = init(C = C, ETA = ETA, type_ = type_)

    training, test = split_dataset(X, y, training_ratio)

    X, labels = training
    for i in 1:size(X, 2)
        model = fit!(model, view(X, :, i), labels[i])
    end

    X, y_true = test

    y_pred = Int64[]
    for i in 1:size(X, 2)
        r = predict(model, view(X, :, i))
        append!(y_pred, r)
    end

    accuracy = calc_accuracy(y_pred, y_true)
    assert(accuracy == 1.0)

    println("ONLINE")
    println("\ttype: $type_")
    println("\taccuracy: $accuracy")
    println("")
end

X = readdlm("data/julia_array/digitsX.txt")
y = readdlm("data/julia_array/digitsy.txt")

println("TEST DIGITS\n")

# Dense matrix
test_batch(X, y, SCW1, training_ratio=0.8)
test_batch(X, y, SCW2, training_ratio=0.8)

test_online(X, y, SCW1, training_ratio=0.8)
test_online(X, y, SCW2, training_ratio=0.8)

X = sparse(X)

# Sparse matrix
test_batch(X, y, SCW1, training_ratio=0.8)
test_batch(X, y, SCW2, training_ratio=0.8)

test_online(X, y, SCW1, training_ratio=0.8)
test_online(X, y, SCW2, training_ratio=0.8)


text:
# NONLINEAR 1-D LONG WAVE MODEL FOR PROPAGATION OVER CONSTANT DEPTH
#CEMAN: COASTAL ECOSYSTEM MANAGEMENT NETWORK
# BY: GERMAN RIVILLAS-OSPINA
# DATE: JUNIO  DE 2020
#LUGAR: BARRANQUILLA-COLOMBIA

using Plots

# CONSTANT VALUES
const g=9.81
const T=100.0
const dt=0.5
const dx=100.0
const Hₒ=200.0

function LW()

    # INITIAL CONDITIONS
    C=(g*Hₒ)^0.5
    L=C*T

    # DEFINING THE DOMAIN
    xmax=0:20.0:L
    xin=length(xmax)
    tin=round.(Int,L)

    #DEFINING VALUES TO PHYSICAL VARIABLES
    η=zeros(1,xin)
    HₒA=fill(Hₒ, (1,xin))
    H=fill(Hₒ, (1,xin))
    Hₙ=zeros(1,xin)
    H₁=zeros(1,xin)
    H₂=zeros(1,xin)
    Hᵣ=zeros(1,xin)
    U=zeros(1,xin)
    U₁=zeros(1,xin)
    U₂=zeros(1,xin)
    Uₙ=zeros(1,xin)

    η₀=1.5
    dtmax=round.(Int,tin/2)
    N=0
    η₁=0.0
    η₂=0.0
    pi=4*atan(1.)
println(xmax)
    while N < dtmax

        #TEMPORAL INCREMENTS
        N=N+1

        #COURANT NUMBER
        TERM1=(N-1)*dt
        TERM2=dx/C

        # FREE SURFACE ELEVATION CALCULATION
        if TERM1<TERM2

            D=sin(2*pi*(N)*(dt/T))
            η[1]=η₀*D+η₁

        else

            A=sin(2*pi*(N-1)*dt/T)
            B=sin((2*pi*(N-1)*dt)/T-(dx/L))
            η₁=η[1]-η₀*A
            η₂=η[2]-η₀*B

            η₁=η₁+(dt/dx)*C*(η₂-η₁)
            D=sin(2*pi*(N)*(dt/T))
            η[1]=η₀*D+η₁

        end

        for i = 2:xin-1
                # DETERMINING WATER DEPTH VARIATION
                Hₙ[1]=HₒA[1]+η[1]

                H₁[i]=(H[i+1]+H[i])*U[i+1]

                H₂[i]=(H[i]+H[i-1])*U[i]

                Hᵣ[i]=H₁[i]-H₂[i]

                Hₙ[i]=H[i]-(dt/(2*dx))*Hᵣ[i]

                # WATER VELOCITIES FIELD CALCULATION CENTRAL SCHEME
                η[i]=Hₙ[i]-HₒA[i]

                U₁[i]=(U[i+1]+U[i])^2-(U[i]+U[i-1])^2
                U₂[i]=g*(η[i]-η[i-1])

                Uₙ[i]=U[i]-(dt/(8dx))U₁[i]-(dt/dx)U₂[i]
                #println(Uₙ[i],"   ", η[i])
        end

            Uₙ[1]=Uₙ[2]
            Uₙ[xin]=Uₙ[xin-1]*(g/Hₙ[xin-1])^0.5

        for i = 1:xin
            H[i]=Hₙ[i]
            U[i]=Uₙ[i]

        end
    end
    #CONVERT RANGE TO VECTORS
    # RETURN U AND η AS VECTORS, NOT 2D ARRAYS
    return collect(xmax), U[1,:], η[1,:]

end
# DEFINING THE RESULTS IN VECTORS
x, u, η = LW()

# PLOTTING RESULTS

    plot(
        x,
        [η u],
         title = ["Surface Elevation" "Velocity field"],
         layout = (2, 1),
         color = ["blue" "green"],
         legend = false
         )

# SAVING PLOT RESULTS IN LOCAL MACHINE
    savefig("G:\\My Drive\\Uninorte\\Modelacion\\Julia\\LongWaves\\result.png")

##_____________________

max_stars_repo_name: cemanetwork/LongWaves

max_issues_repo_name: cemanetwork/LongWaves

max_forks_repo_name: cemanetwork/LongWaves


text:
function main(args)
    return args
end
# Load inline arguments for worker unit
meta_key_file, num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this = parse(Int64, this)

println("Dataset: " * meta_key_file)
println("Worker " * string(this) * " of " * string(num_workers))
#
# Load packages and functions
include("../src/init.jl")
include("../src/fct.jl")
include("../src/worker_fct.jl")
try using UpROOT
catch
    Pkg.add("UpROOT")
    using UpROOT
end
try using LegendUpROOTIO
catch
    Pkg.add(url="https://github.com/legend-exp/LegendUpROOTIO.jl")
    using LegendUpROOTIO
end
using LegendUpROOTIO: read_mgtevent, raw2mgtevent, mgdo2legend
using DataFrames
data_type = "cal" # phy
dataset = "v07.01"
current_dir = pwd()

#
# Gather filepaths to Tier1 and Tier4 data
meta_keys = CSV.read(meta_key_file, Table);

channels = 0:1:36
event_step = Int(1e4)

current_dir = pwd()
ddir = "/remote/ceph/group/gerda/data/phase2/blind/" * dataset * "/gen/"
cd(ddir)
filenames1 = []
filenames4 = []
for meta_key in meta_keys
    filename1 = ddir * glob("tier1/ged/" * data_type * "/" * split(meta_key[1], "-")[2] * "/" * meta_key[1] * "*.root")[1]
    filename4 = ddir * glob("tier4/all/" * data_type * "/" * split(meta_key[1], "-")[2] * "/" * meta_key[1] * "*.root")[1]
    push!(filenames1, filename1)
    push!(filenames4, filename4)
end
cd(current_dir)

#
# Set base path depending on data type
if data_type == "phy"
    base_path = "../../waveforms/test/data/raw_" * dataset * "/"
elseif data_type == "cal"
    base_path = "../../waveforms/test/calib/raw_" * dataset * "/"
end

#
# Load workers log file
log_file = base_path * "log-" * string(this) * ".json"
if !isdir(base_path)
    mkpath(base_path)
end
if !isfile(log_file)
    open(log_file, "w") do f
        JSON.print(f, Dict(), 4)
    end
    global log = Dict()
else
    global log = JSON.parsefile(log_file)
end;

#
# Select range for worker depending on number of workers
filenames1 = get_share_for_worker(filenames1, num_workers, this)
filenames4 = get_share_for_worker(filenames4, num_workers, this)

for i in eachindex(filenames4)
    start_t = now()
    run = parse(Int64, split(split(filenames4[i], "gerda-run")[2], "-")[1])
    filename1 = filenames1[i]
    filename4 = filenames4[i]

    @info("Run " * string(run) * " | file " * string(i) * " of total " * string(length(filenames4)))
    log = JSON.parsefile(log_file)
    if length(keys(log)) >= 1
        estimation = 0
        for k in keys(log)
            estimation += log[k]
        end
        estimation /= length(keys(log))
        fileinfo = stat(filename4)
        est_t = round(fileinfo.size * estimation / 60, digits=1)
        @info("Estimated time for this file: " * string(est_t) * " min")
        total_filesize = 0
        for f in i:length(filenames4)
            total_filesize += stat(filenames4[f]).size
        end
        est_t = total_filesize * estimation / 60
        unit = " s"
        if est_t/60/60/24 > 1
            est_t /= 60*60*24
            unit = " d"
        elseif est_t/60/60 > 1
            est_t /= 60*60
            unit = " h"
        elseif est_t/60 > 1
            est_t /= 60
            unit = " min"
        end
        @info("Estimated time till completion: " * string(round(est_t, digits=1)) * unit)
    end
    @info("------------------------------")

    if !(filename4 in keys(JSON.parsefile(log_file)))
        @info("Load Tier4")
        @time tier4 = Table(TFile(filename4)["tier4"]);
        temp_E = sum.(tier4.energy[:])

        index_1 = findall(x->x == 1, tier4.multiplicity[:])
        index_2 = findall(y->y > 300, temp_E[index_1])
        index_3 = findall(x->x == 0, tier4.isBL[index_1[index_2]])
        index_4 = findall(x->x == 0, tier4.isTP[index_1[index_2[index_3]]])
        filtered_index = index_1[index_2[index_3[index_4]]]
        if length(filtered_index) > 0
            @info("Number of events: " * string(length(temp_E)))
            @info("Number of events after filter: " * string(length(filtered_index)))
            steps = []
            s = 1
            while s <= length(filtered_index)
                a = s
                b = s + event_step - 1 <= length(filtered_index) ? s+event_step-1 : length(filtered_index)
                push!(steps, [a,b])
                s += event_step
            end
            for s in steps
                result = Table( energy       = [],
                        run          = [],
                        channel      = [],
                        AEvetoed     = [],
                        datasetID    = [],
                        AEclassifier = [],
                        auxWaveform  = [],
                        waveform     = [])
                Base.run(`clear`)
                run_str = split(basename(filename4), "-")[2]
                @info("Run " * run_str * " | file " * string(i) * " of total " * string(length(filenames4)))
                @info("Step " * string(findfirst(x->x == s, steps)) * " of " * string(length(steps)))
                file = base_path * run_str * "/"
                if !isdir(file)
                    mkpath(file)
                end
                file *= basename(filename4) * lpad(findfirst(x->x == s, steps), 4, "0") * ".h5"
                if !isfile(file)

                    tmp = filtered_index[s[1]:s[2]]
                    @info("Event " * string(tmp[1]) * " until " * string(tmp[end]))
                    tier4 = Table(TFile(filename4)["tier4"])[tmp];
                    @info("Load Tier1")
                    @time treeTier1 = TFile(filename1)["MGTree"].event;
                    @info("Apply first filter to Tier1 & Tier4")
                    treeTier1 = TypedTables.Table(raw2mgtevent.(treeTier1[tmp]))
                    @time auxWaveforms = treeTier1.fAuxWaveforms
                    @time waveforms = treeTier1.fWaveforms
                    if length(tmp) > 0
                        run_str = split(basename(filename4), "-")[2]
                        run = parse(Int64, split(run_str, "run")[2])

                        @showprogress for ch in channels

                            tmp = findall(y->y[ch + 1] != 0, tier4.energy)
                            if length(tmp) > 0
                                auxWaveform = []
                                waveform = []
                                for i in tmp
                                    if data_type == "phy"
                                        push!(auxWaveform, auxWaveforms[i][ch + 1].wf)
                                        push!(waveform, waveforms[i][ch + 1].wf)
                                    else
                                        push!(auxWaveform, auxWaveforms[i][1].wf)
                                        push!(waveform, waveforms[i][1].wf)
                                    end
                                end

                                AEvetoed     = []
                                datasetID    = []
                                AEclassifier = []
                                for i in tmp
                                    push!(AEvetoed,     tier4[i].isAoEvetoed[ch+1])
                                    push!(datasetID,    tier4[i].datasetID[ch+1])
                                    push!(AEclassifier, tier4[i].AoEclassifier[ch+1])
                                end

                                append!(result, Table(  energy       = sum.(tier4.energy[tmp]),
                                                        run          = Int.(zeros(length(tmp)) .+ run),
                                                        channel      = Int.(zeros(length(tmp)) .+ ch),
                                                        AEvetoed     = AEvetoed,
                                                        datasetID    = datasetID,
                                                        AEclassifier = AEclassifier,
                                                        auxWaveform  = StructArray{RDWaveform}(Array{RDWaveform,1}(auxWaveform)),
                                                        waveform     = StructArray{RDWaveform}(Array{RDWaveform,1}(waveform))
                                        )
                                    )
                            end
                        end
                    end
                    if size(result,1) > 0
                        file = base_path * run_str * "/"
                        file *= basename(filename4) * lpad(findfirst(x->x == s, steps), 4, "0") * ".h5"
                        HDF5.h5open(file, "w") do h5f
                            LegendHDF5IO.writedata( h5f, "data", Table( energy       = float.(result.energy),
                                                                        run          = Int.(result.run),
                                                                        channel      = Int.(result.channel),
                                                                        AEvetoed     = Int.(result.AEvetoed),
                                                                        datasetID    = Int.(result.datasetID),
                                                                        AEclassifier = float.(result.AEclassifier),
                                                                        auxWaveform  = StructArray{RDWaveform}(Array{RDWaveform,1}(result.auxWaveform)),
                                                                        waveform     = StructArray{RDWaveform}(Array{RDWaveform,1}(result.waveform))))
                        end
                    end
                end
            end
            log = JSON.parsefile(log_file)
            dt = now() - start_t
            fileinfo = stat(filenames4[1])
            fileinfo.size
            log[filename4] = (dt.value / 1000) / fileinfo.size
            open(log_file, "w") do f
                JSON.print(f, log, 4)
            end
        end
    end
#     IJulia.clear_output(true)
    Base.run(`clear`)
    log = JSON.parsefile(log_file)
    @info(string(length(log)) * " of " * string(length(filenames4)) * " done!")
end

max_stars_repo_path: calib/01 - Data Calibration root2h5.jl

max_issues_repo_path: calib/01 - Data Calibration root2h5.jl

max_forks_repo_path: calib/01 - Data Calibration root2h5.jl


text:
# NegativeBinomial distribution

export NegativeBinomial
import Base
using StatsFuns
using SpecialFunctions

@parameterized NegativeBinomial(r,p) ≪ CountingMeasure(ℤ[0:∞])

(d::NegativeBinomial ≪ ::CountingMeasure{IntegerRange{a,b}}) where {a,b} = a ≤ 0 && b ≥ d.n

(::CountingMeasure{IntegerRange{a,b}} ≪ ::NegativeBinomial) where {a,b} = a ≥ 0 && b ≤ d.n

###############################################################################
@kwstruct NegativeBinomial(r, p)

function logdensity(d::NegativeBinomial{(:r, :p)}, y)
    (r, p) = (d.r, d.p)
    return -log(y + r) - logbeta(r, y+1) + y * log(p) + r * log1p(-p)
end

###############################################################################
@kwstruct NegativeBinomial(r, logitp)

function logdensity(d::NegativeBinomial{(:r, :logitp)}, y)
    (r, logitp) = (d.r, d.logitp)
    return -log(y + r) - logbeta(r, y+1) - y * log1pexp(-logitp) - r * log1pexp(logitp)
end

###############################################################################
@kwstruct NegativeBinomial(r,λ)
# mean λ, as in Poisson
# Converges to Poisson as r→∞

function logdensity(d::NegativeBinomial{(:r, :λ)}, y)
    (r, λ) = (d.r, d.λ)
    return -log(y + r) - logbeta(r, y+1) + y * log(λ) + r * log(r) - (y + r) * log(r + λ)
end

distproxy(d::NegativeBinomial{(:r, :p)}) = Dists.NegativeBinomial(d.r, d.p)
distproxy(d::NegativeBinomial{(:n,:logitp)}) = Dists.NegativeBinomial(d.n, logistic(d.logitp))

function distproxy(d::NegativeBinomial{(:r,:λ)})
    p = d.λ / (d.r + d.λ)
    return Dists.NegativeBinomial(d.r, p)
end

# https://en.wikipedia.org/wiki/Negative_binomial_distribution#Gamma%E2%80%93Poisson_mixture
function Base.rand(rng::AbstractRNG, d::NegativeBinomial{(:r,:λ)})
    r = d.r
    λ = d.λ
    μ = rand(rng, Dists.Gamma(r, λ/r))
    return rand(rng, Dists.Poisson(μ))
end

asparams(::Type{<:NegativeBinomial}, ::Val{:p}) = as𝕀
asparams(::Type{<:NegativeBinomial}, ::Val{:logitp}) = asℝ
asparams(::Type{<:NegativeBinomial}, ::Val{:r}) = asℝ₊
asparams(::Type{<:NegativeBinomial}, ::Val{:λ}) = asℝ₊


text:
module TestSetup

using Intan

include("board_test.jl")
include("gui_test.jl")
include("sorting_test.jl")
include("registers_test.jl")
include("tasks_test.jl")
include("filter_test.jl")

end


text:
for (T1, T2) in ((:T, :(float(T))), (:(Union{Missing,T}), :(Union{Missing,float(T)})))
  @eval begin

    # unweighted windowed function application

    function rolling(fun::Function, data::AbstractVector{$T1}, windowspan::Int) where {T}
        nvals  = nrolled(length(data), windowspan)
        offset = windowspan - 1
        result = zeros($T2, nvals)

        @inbounds for idx in eachindex(result)
            result[idx] = fun( view(data, idx:idx+offset) )
        end

        return result
    end

    # weighted windowed function application
    function rolling(fun::Function, data::AbstractVector{$T1}, windowspan::Int, weighting::F) where
                        {T, N<:Number, F<:Vector{N}}

        length(weighting) != windowspan &&
           throw(WeightsError(length(weighting), windowspan))

        nvals  = nrolled(length(data), windowspan)
        offset = windowspan - 1
        result = zeros($T2, nvals)
        curwin = zeros($T2, windowspan)
        v = view(weighting,1:length(weighting))

        @inbounds for idx in eachindex(result)
           for i=1:windowspan
               j = i + idx - 1
               curwin[i] = data[j] * v[i]
           end
           result[idx] = fun( curwin )
        end

        return result
    end

    rolling(fun::Function, data::AbstractVector{$T1}, windowspan::Int, weighting::W) where
                    {T, N<:Number, W<:AbstractWeights} =
        rolling(fun, data, windowspan, weighting.values)

    # unweighted windowed function tapering

    function tapers(fun::Function, data::AbstractVector{$T1}) where {T}
        nvals  = length(data)
        result = zeros($T2, nvals)

        @inbounds for idx in nvals:-1:1
            result[idx] = fun( view(data, 1:idx) )
        end

        return result
    end

    function tapers(fun::Function, data::AbstractVector{$T1}, ntocopy::Int) where {T}
        nvals  = length(data)
        ntocopy = min(nvals, max(0, ntocopy))

        result = zeros($T2, nvals)
        if ntocopy > 0
           result[1:ntocopy] = data[1:ntocopy]
        end
        ntocopy += 1

        @inbounds for idx in nvals:-1:ntocopy
            result[idx] = fun( view(data, 1:idx) )
        end

        return result
    end

    function tapers(fun::Function, data::AbstractVector{$T1},
                    trailing_data::AbstractVector{$T1}) where {T}

        ntrailing = axes(trailing_data)[1].stop
        nvals  = length(data) + ntrailing
        result = zeros($T2, nvals)

        result[1:ntrailing] = trailing_data[1:ntrailing]
        ntrailing += 1

        @inbounds for idx in nvals:-1:ntrailing
            result[idx] = fun( view(data, 1:idx) )
        end

        return result
    end

    # weighted windowed function tapering

    function tapers(fun::Function, data::AbstractVector{$T1}, weighting::F) where
                     {T, N<:Number, F<:Vector{N}}

        nvals  = length(data)
        nweighting = length(weighting)

        nweighting == nvals || throw(WeightsError(nweighting, nvals))

        result = zeros($T2, nvals)

        @inbounds for idx in nvals:-1:1
            wts = normalize(view(weighting, (nweighting-idx+1):nweighting))
            result[idx] = fun( view(data, 1:idx) .* wts  )
        end

        return result
    end

    tapers(fun::Function, data::AbstractVector{$T1}, weighting::W) where
                    {T,W<:AbstractWeights} =
        tapers(fun, data, weighting.values)

  end
end


text:
using Test
import Roots.newton, Roots.halley

@testset "Test Newton, Halley, Schroder methods" begin

    @test abs(newton(sin, cos, 0.5) - 0.0) <= 100*eps(1.0)
    @test newton(cos, x -> -sin(x), 1.0)  ≈ pi/2
    @test newton(x -> x^2 - 2x - 1, x -> 2x - 2, 3.0)  ≈ 2.414213562373095
    @test abs(newton(x -> exp(x) - cos(x), x -> exp(x) + sin(x), 3.0) - 0.0) <= 1e-14
    @test halley(x -> x^2 - 2x - 1,x -> 2x - 2,x -> 2, 3.0)  ≈ 2.414213562373095
    a = halley(x -> exp(x) - cos(x),
               x -> exp(x) + sin(x),
               x -> exp(x) + cos(x), 3.0)
    @test abs(a - 0.0) <= 1e-14

    # find_zero calls
    @test find_zero((x -> x^2 - 2x - 1,x -> 2x - 2), 3.0, Roots.Newton())  ≈ 2.414213562373095
    @test find_zero((x -> x^2 - 2x - 1,x -> 2x - 2,x -> 2), 3.0, Roots.Halley())  ≈ 2.414213562373095
    @test find_zero((x -> x^2 - 2x - 1,x -> 2x - 2,x -> 2), 3.0, Roots.Schroder())  ≈ 2.414213562373095

    ## test with Complex input

    @test real(Roots.newton(x ->  x^3 - 1, x ->  3x^2, 1+im)) ≈ 1.0
    @test real(Roots.newton(x ->  x^3 - 1, x ->  3x^2, 1+10im)) ≈ (-1/2)

    ## Issue #143 test with new interface
    Roots.newton(sin, cos, 3.0) ≈ π # uses find_zero
    Roots.newton((sin,cos), 3.0) ≈ π # uses simple

    fdf = x -> (sin(x), sin(x)/cos(x))  # (f, f/f')
    @test Roots.find_zero(fdf, 3.0, Roots.Newton())  ≈ π # uses find_zero
    Roots.newton(fdf, 3.0)  ≈ π # uses simple

    fdfdf = x -> (sin(x), sin(x)/cos(x), -cos(x)/sin(x), sin(x)/cos(x))   # (f, f/f', f'/f'', f''/f''')
    @test Roots.find_zero(fdfdf, 3.0, Roots.Halley()) ≈ π

    # check that functions with multiple return values can work with other
    # methods
    for M in [Roots.Schroder(), Roots.Halley(), Roots.Newton(), Roots.Order1(), Roots.Order0()]
        @test Roots.find_zero(fdfdf, 3.0, M) ≈ π  # can pass function to others
    end

    for M in [Roots.Bisection(), Roots.A42(), Roots.AlefeldPotraShi()]
        @test Roots.find_zero(fdfdf, (3.0, 4.0), M) ≈ π  # can pass function to others
    end

    @test find_zero(x -> (x^2 -2, (x^2-2)/2x), 1.0, Roots.Newton(), Roots.Bisection()) ≈ sqrt(2)

end


text:
export ProjectionVI

function ProjectionVI(model::AbstractVIModel, x0::AbstractVector; rho0::Float64 = 0.5, kwargs...)
  stp = GenericStopping(model, x0; kwargs...)
  return ProjectionVI(stp, rho0 = rho0)
end

function abresidual!(model, xk, rho, Fx) # xk + rho * F(xk)
  residual!(model, xk, Fx)
  Fx .*= rho
  Fx .+= xk
  return Fx
end

function ProjectionVI(stp::AbstractStopping; rho0::Float64 = 0.5)
  xk = stp.current_state.x
  xkp = similar(xk)
  rho = rho0
  Fx = similar(xk)

  OK = update_and_start!(stp)
  while !OK
    abresidual!(stp.pb, xk, rho, Fx)
    project!(stp.pb, Fx, xk) # possible failure here

    if norm(xk - xkp, Inf) < stp.meta.atol * rho
      stp.meta.optimal = true
    end
    OK = update_and_stop!(stp, x = xk)
  end

  return xk
end


text:
# uhdOverNetwork.jl
# This module provides bind to monitor uhdOverNetwork based device from a Host PC
# As the uhdOverNetwork is a SoC based device with ARM, the Host <-> uhdOverNetwork link is assumed to be done with the use of sockets.
# With this module we construct utils to use the uhdOverNetwork device with same functions as it is a classic UHD device.

# This is quite similar with the uhd_network_mode with some important difference
# - We use a specific communication system socket
# - uhdOverNetwork side: a Julia script has to be launched. This script can also be modified to embed additional processing

module SDROverNetworks

# --- Module dependency
using Sockets
using ZMQ
using Printf

# --- Print radio config
include("../Printing.jl");
using .Printing

# --- Method extension
import Sockets:send;
import Sockets:recv;
import ZMQ:recv;
import ZMQ:send;
import ZMQ:close;
import Sockets:close;

# --- Symbol exportation
export openUhdOverNetwork;
export updateCarrierFreq!;
export updateSamplingRate!;
export updateGain!;
export recv;
export recv!;
export print;

struct SocketsuhdOverNetwork
    ip::IPAddr;
    rtcSocket::Socket;
    rttSocket::Socket;
    brSocket::Socket;
end
mutable struct Configuration
	carrierFreq::Float64;
	samplingRate::Float64;
	gain::Union{Float64,Int};
	Antenna::String;
    packetSize::Int;
end
mutable struct MD
    intPart::Int32;
    fracPart::Cdouble;
	error::Int32;
end
mutable struct SDROverNetworkRx
    sockets::SocketsuhdOverNetwork;
	carrierFreq::Float64;
	samplingRate::Float64;
	gain::Union{Int,Float64};
	antenna::String;
    packetSize::Csize_t;
    released::Int;
end
mutable struct SDROverNetworkTx
    sockets::SocketsuhdOverNetwork;
	carrierFreq::Float64;
	samplingRate::Float64;
	gain::Union{Int,Float64};
	antenna::String;
    packetSize::Csize_t;
    released::Int;
end
mutable struct SDROverNetwork
    radio::String;
    rx::SDROverNetworkRx;
    tx::SDROverNetworkTx;
end

export SDROverNetwork;

function initSockets(ip::String)
    # --- Format e310 IP address
    e310Address  =  IPv4(ip);
    # --- Configuration Socket
    # Socket used to send configuration and get configuration back
    rtcSocket   = ZMQ.Socket(REQ);
    ZMQ.connect(rtcSocket,"tcp://$e310Address:5555");

    # --- Tx socket
    # Socket used to transmit data from Host to e310 (Tx link)
    rttSocket = ZMQ.Socket(REP);
    ZMQ.connect(rttSocket,"tcp://$e310Address:9999");

    # --- Rx socket
    # Socket used for e310 to broacast Rx stream
    brSocket    =  Socket(SUB);   # Define IPv4 adress
    tcpSys		= string("tcp://$e310Address:1111");
    ZMQ.subscribe(brSocket);
    ZMQ.connect(brSocket,tcpSys);

    # --- Global socket packet
    sockets     = SocketsuhdOverNetwork(e310Address,rtcSocket,rttSocket,brSocket);
    return sockets
end

"""
Open a uhdOverNetwork remote device and initialize the sockets
--- Syntax

openUhdOverNetwork(carrierFreq,samplingRate,gain,antenna="RX2";ip="192.168.10.11")
# --- Input parameters
- carrierFreq	: Desired Carrier frequency [Union{Int,Float64}]
- samplingRate	: Desired bandwidth [Union{Int,Float64}]
- gain		: Desired Rx Gain [Union{Int,Float64}]
Keywords
- ip	  : uhdOverNetwork IP address
- antenna		: Desired Antenna alias  (default "TX-RX") [String]
# --- Output parameters
- structuhdOverNetwork    : Structure with uhdOverNetwork parameters [SDROverNetwork]
"""
function openUhdOverNetwork(carrierFreq,samplingRate,gain;antenna="RX2",addr="192.168.10.11")
    # --- Create the Sockets
    sockets = initSockets(addr);
    # --- Create the initial configuration based on input parameters
    rx    = SDROverNetworkRx(
        sockets,
        carrierFreq,
        samplingRate,
        gain,
        antenna,
        0,
        0
    );
    tx    = SDROverNetworkTx(
        sockets,
        carrierFreq,
        samplingRate,
        gain,
        antenna,
        0,
        0
    );
    # --- Instantiate the complete radio
    radio = SDROverNetwork(
                           addr,
                           rx,
                           tx
                          );
    # --- Update the radio based on input parameters
    updateCarrierFreq!(radio,carrierFreq);
    updateSamplingRate!(radio,samplingRate);
    updateGain!(radio,gain);
    # Get socket size
    requestConfig!(radio);
    # --- Print the configuration
    print(radio);
    # --- Return the final object
    return radio;
end

function Base.close(radio::SDROverNetwork)
    # @info "coucou"
    # --- We close here all the related sockets
    close(radio.rx.sockets.rtcSocket);
    close(radio.rx.sockets.rttSocket);
    close(radio.rx.sockets.brSocket);
end

sendConfig(uhdOverNetwork::SDROverNetwork,mess) = send(uhdOverNetwork.rx.sockets.rtcSocket,mess)

function updateCarrierFreq!(uhdOverNetwork::SDROverNetwork,carrierFreq)
    # --- Create char with command to be transmitted
    # strF        = "global carrierFreq = $carrierFreq";
    strF        = "Dict(:updateCarrierFreq=>$carrierFreq);";
    # --- Send the command
    sendConfig(uhdOverNetwork,strF);
    # --- Get the effective radio configuration
    config = getuhdOverNetworkConfig(uhdOverNetwork);
    # --- Update the uhdOverNetwork object based on real radio config
    uhdOverNetwork.rx.carrierFreq = config.carrierFreq;
    uhdOverNetwork.tx.carrierFreq = config.carrierFreq;
    return uhdOverNetwork.rx.carrierFreq;
end

function updateSamplingRate!(uhdOverNetwork::SDROverNetwork,samplingRate)
    # --- Create char with command to be transmitted
    strF        = "Dict(:updateSamplingRate=>$samplingRate);";
    # --- Send the command
    sendConfig(uhdOverNetwork,strF);
    # --- Get the effective radio configuration
    config = getuhdOverNetworkConfig(uhdOverNetwork);
    # --- Update the uhdOverNetwork object based on real radio config
    uhdOverNetwork.rx.samplingRate = config.samplingRate;
    uhdOverNetwork.tx.samplingRate = config.samplingRate;
    return uhdOverNetwork.rx.samplingRate;
end

function updateGain!(uhdOverNetwork::SDROverNetwork,gain)
    # --- Create char with command to be transmitted
    strF        = "Dict(:updateGain=>$gain);";
    # --- Send the command
    sendConfig(uhdOverNetwork,strF);
    # --- Get the effective radio configuration
    config = getuhdOverNetworkConfig(uhdOverNetwork);
    # --- Update the uhdOverNetwork object based on real radio config
    uhdOverNetwork.rx.gain = config.gain;
    uhdOverNetwork.tx.gain = config.gain;
    return uhdOverNetwork.rx.gain;
end

function requestConfig!(uhdOverNetwork::SDROverNetwork);
    # --- Create char with command to be transmitted
    strF        = "Dict(:requestConfig=>1);";
    # --- Send the command
    sendConfig(uhdOverNetwork,strF);
    # --- Get the effective radio configuration
    config = getuhdOverNetworkConfig(uhdOverNetwork);
    uhdOverNetwork.rx.packetSize = config.packetSize;
    uhdOverNetwork.tx.packetSize = config.packetSize;
end

function setRxMode(uhdOverNetwork::SDROverNetwork)
    # --- Create char with command to be transmitted
    strF        = "Dict(:mode=>:rx);";
    # --- Send the command
    sendConfig(uhdOverNetwork,strF);
    receiver = recv(uhdOverNetwork.rx.sockets.rtcSocket);
end
function recv(uhdOverNetwork::SDROverNetwork,packetSize)
    # --- Create container
    sig = Vector{Complex{Cfloat}}(undef,packetSize);
    # --- fill the stuff
    recv!(sig,uhdOverNetwork);
    return sig;
end
function recv!(sig::Vector{Complex{Cfloat}},uhdOverNetwork::SDROverNetwork;packetSize=0,offset=0)
    # setRxMode(uhdOverNetwork);
	# --- Defined parameters for multiple buffer reception
	filled		= false;
	# --- Fill the input buffer @ a specific offset
	if offset == 0
		posT		= 0;
	else
		posT 		= offset;
	end
	# --- Managing desired size and buffer size
	if packetSize == 0
		# --- Fill all the buffer
		packetSize	= length(sig);
	else
		packetSize 	= packetSize;
		# --- Ensure that the allocation is possible
		@assert packetSize < (length(sig)+posT) "Impossible to fill the buffer (number of samples > residual size";
	end
	while !filled
		# --- Get a buffer: We should have radio.packetSize or less
		# radio.packetSize is the complex size, so x2
		(posT+uhdOverNetwork.rx.packetSize> packetSize) ? n = packetSize - posT : n = uhdOverNetwork.rx.packetSize;
        # --- UDP recv. This allocs. This is bad. No idea how to use prealloc pointer without rewriting the stack.
        tmp = reinterpret(Complex{Cfloat},recv(uhdOverNetwork.rx.sockets.brSocket));
        sig[posT .+ (1:n)] .= @view tmp[1:n];
		# --- Update counters
		posT += n;
		# --- Breaking flag
		(posT == packetSize) ? filled = true : filled = false;
	end
	return posT
end

#FIXME: We have setTxMode call before each Tx. Shall we do a setRxMode before each rx frame ?
function setTxMode(uhdOverNetwork::SDROverNetwork)
    # --- Create char with command to be transmitted
    strF        = "Dict(:mode=>:tx);";
    # --- Send the command
    sendConfig(uhdOverNetwork,strF);
    receiver = recv(uhdOverNetwork.rx.sockets.rtcSocket);
end
function setTxBufferMode(uhdOverNetwork::SDROverNetwork)
    # --- Create char with command to be transmitted
    strF        = "Dict(:mode=>:txbuffer);";
    # --- Send the command
    sendConfig(uhdOverNetwork,strF);
    receiver = recv(uhdOverNetwork.rx.sockets.rtcSocket);
end

function sendBuffer(buffer::Vector{Complex{Cfloat}},uhdOverNetwork::SDROverNetwork)
  # --- Create char with command to be transmitted
    strF        = "Dict(:buffer:=>$buffer);";
    # --- Send the command
    sendConfig(uhdOverNetwork,strF);
    config = getuhdOverNetworkConfig(uhdOverNetwork);
    return config;
end

function send(sig::Vector{Complex{Cfloat}},uhdOverNetwork::SDROverNetwork,cyclic=false;maxNumSamp=nothing)
    # --- Setting radio in Tx mode
    # setTxMode(uhdOverNetwork);
    nS  = 0;
    it  = length(sig);
    if cyclic == true
        # ----------------------------------------------------
        # --- We handle Tx through metadata
        # ----------------------------------------------------
        # setTxBufferMode(uhdOverNetwork):
        config = sendBuffer(sig,uhdOverNetwork);
    else
        # ----------------------------------------------------
        # --- Using RTT socket to handle data exchange
        # ----------------------------------------------------
        try
            # --- First while loop is to handle cyclic transmission
            # It turns to false in case of interruption or cyclic to false
            while (true)
                # --- Wait for RTT
                rtt = ZMQ.recv(uhdOverNetwork.tx.sockets.rttSocket);
                # --- Sending data to Host
                ZMQ.send(uhdOverNetwork.tx.sockets.rttSocket,sig);
                # --- Update counter
                nS += it;
                # --- Detection of cyclic mode
                (maxNumSamp !== nothing && nS > maxNumSamp) && break
                (cyclic == false ) && break
                # --- Forcing refresh
                yield();
            end
        catch e;
            # --- Interruption handling
            print(e);
            print("\n");
            @info "Interruption detected";
            return 0;
        end
    end
    return nS;
end

function getuhdOverNetworkConfig(uhdOverNetwork::SDROverNetwork)
    receiver = recv(uhdOverNetwork.rx.sockets.rtcSocket);
    res =  Meta.parse(String(receiver))
    config = eval(res);
    return Configuration(config...);
end

function Base.print(uhdOverNetwork::SDROverNetwork)
    strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n Rx Gain: %2.2f dB\n",uhdOverNetwork.rx.carrierFreq/1e6,uhdOverNetwork.rx.samplingRate/1e6,uhdOverNetwork.rx.gain);
    @inforx "Current uhdOverNetwork Configuration in Rx mode\n$strF";
    strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n Rx Gain: %2.2f dB\n",uhdOverNetwork.tx.carrierFreq/1e6,uhdOverNetwork.tx.samplingRate/1e6,uhdOverNetwork.tx.gain);
    @infotx "Current uhdOverNetwork Configuration in Tx mode\n$strF";
end

function getMD(uhdOverNetwork::SDROverNetwork)
    # --- Create char with command to be transmitted
    strF        = "Dict(:requestMD=>1);";
    # --- Send the command
    sendConfig(uhdOverNetwork,strF);
    # --- Get the MD back
    receiver = recv(uhdOverNetwork.rx.sockets.rtcSocket);
    res =  Meta.parse(String(receiver))
    # --- Convert to a MD structure
    md  = eval(res)
    return md;
end

end


text:
using Base: @deprecate

function Base.:|>(f::Transducer, g::Transducer)
    if Base.JLOptions().depwarn == 1
        @warn """
        **Quick upgrade guide**: `f::Transducer |> g::Transducer` is deprecated.
        Instead of

            collect(Filter(f) |> Map(g), xs)
            foldl(+, Filter(f) |> Map(g), xs)

        it's now recommended to use `|>` with input collection

            xs |> Filter(f) |> Map(g) |> collect
            foldl(+, xs |> Filter(f) |> Map(g))

        If there is no input collection; e.g.,

            foldl(right, GroupBy(key, Filter(f) |> Map(g), push!!), xs)

        use `opcompose` instead:

            foldl(right, GroupBy(key, opcompose(Map(f), Filter(g)), push!!), xs)

        For more information, see:
        https://juliafolds.github.io/Transducers.jl/dev/howto/upgrade-to-ixf/
        """ maxlog = 1
    end
    Base.depwarn(
        "`f::Transducer |> g::Transducer` is deprecated. " *
        "Use `xs |> f |> g |> collect` instead of `collect(f |> g, xs)`. " *
        "For more information, see " *
        "https://juliafolds.github.io/Transducers.jl/dev/howto/upgrade-to-ixf/",
        :|>,
    )
    return g ∘ f
end

@deprecate Distinct() Unique()
@deprecate TeeZip(xf) ZipSource(xf) false
@deprecate Keep(f) opcompose(Map(f), NotA(Nothing))

"""
    Keep(f)

Pass non-`nothing` output of `f` to the inner reducing step.

!!! warning

    `Keep(f)` is a deprecated. Use `... |> Map(f) |> NotA(Nothing)`.
    If `f` does not return a `Some`, [`KeepSomething`](@ref) can also
    be used.

# Examples
```julia
julia> using Transducers

julia> xf = Keep() do x
           if x < 3
               x + 1
           end
       end;

julia> collect(xf, 1:5)
2-element Vector{Int64}:
 2
 3
```
"""
Keep

# It wasn't a public API but to play on the safe side:
@deprecate induction(itr) extract_transducer(itr) false

function Base.mapfoldl(xform::Transducer, step, itr; kwargs...)
    Base.depwarn(_MAPFOLDL_DEPWARN, (:mapfoldl, Symbol("mapfoldl##kw")))
    return foldl(whencomplete(step, step), xform, itr; kwargs...)
end

# disambiguation
function Base.mapfoldl(xform::Transducer, step, itr::Foldable; kwargs...)
    Base.depwarn(_MAPFOLDL_DEPWARN, (:mapfoldl, Symbol("mapfoldl##kw")))
    return foldl(whencomplete(step, step), xform, itr; kwargs...)
end

function _mapreduce(xform, step, itr, kwargs)
    Base.depwarn(_MAPREDUCE_DEPWARN, (:mapreduce, Symbol("mapreduce##kw")))
    return reduce(whencomplete(step, step), xform, itr; kwargs...)
end

Base.mapreduce(xform::Transducer, step, itr; kwargs...) =
    _mapreduce(xform, step, itr, kwargs)

# Disambiguation:
Base.mapreduce(xform::Transducer, step, itr::AbstractArray; kwargs...) =
    _mapreduce(xform, step, itr, kwargs)
Base.mapreduce(xform::Transducer, step, itr::AbstractArrayOrBroadcasted; kwargs...) =
    _mapreduce(xform, step, itr, kwargs)
Base.mapreduce(xform::Transducer, step, itr::Base.SkipMissing{<:AbstractArray}; kwargs...) =
    _mapreduce(xform, step, itr, kwargs)
Base.mapreduce(xform::Transducer, step, itr::Number; kwargs...) =
    _mapreduce(xform, step, itr, kwargs)

import Base: reduce
@deprecate reduce(rf, xf::Transducer, itr; kw...) foldxt(rf, xf, itr; kw...) false
@deprecate reduce(rf, itr::Foldable; kw...) foldxt(rf, itr; kw...) false
using Base: reduce

"""
    reduce(rf, [xf,] itr)

Multi-threaded reduce.

!!! warning
    `reduce` is deprecated. Use [`foldxt`](@ref) instead.
"""
reduce

@deprecate dreduce(rf, xf, itr; kw...) foldxd(rf, xf, itr; kw...)
@deprecate dreduce(rf, itr; kw...) foldxd(rf, itr; kw...)

"""
    dreduce(rf, [xf,] itr)

Distributed reduce.

!!! warning
    `dreduce` is deprecated. Use [`foldxd`](@ref) instead.
"""
dreduce


text:
using PowerModels
using NetDecOPF
using DualDecomposition
using Gurobi, Ipopt
const DD = DualDecomposition

sub_optimizer = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
optimizer = optimizer_with_attributes(Gurobi.Optimizer, "OutputFlag" => 0)

file = "./examples/case5.m"
data = parse_file(file)
partitions = [[2, 3, 4], [1, 5]]
dn_model = decompose(data, partitions, ACRPowerModel, NetDecOPF.build_acopf_with_free_lines)
algo = init_DD_algo(dn_model)
set_subnet_optimizer!(dn_model, sub_optimizer)
LM = DD.BundleMaster(BM.ProximalMethod, optimizer)
DD.run!(algo, LM)


text:
module ExponentialFamilies

using LinearAlgebra
using Distributions

include("types.jl")
include("normal.jl")

end # module


text:
# ------------------------------------------------------------------- #
# Stochastic Volatility Exogenous Process
# ------------------------------------------------------------------- #
"""
This is the type that holds the parameters for a constant volatility
exogenous process from the BCFL paper. Namely the process follows:

```
\log(\hat{z}_{t+1}) = (1 - 2γ) \log(\hat{z}_{t}) + (σz1/2) ϵ1 - (σz2/2) ϵ2
\log(g_{t+1}) = (σz1/2) ϵ1 + (σz2/2) ϵ2
v_{t+1} = (1-ϕv) vbar + ϕv v_{t} + τ ϵ3
```

The number of epsilons we use for quadrature in each i is given by nϵi
"""
immutable ExogSV <: AbstractExog
    γ::Float64
    σz2::Float64
    vbar::Float64
    φv::Float64
    τ::Float64
    nϵ1::Int
    nϵ2::Int
    nϵ3::Int
end

# Small self-explanatory helper functions
init(e::ExogSV) = [0.0, e.vbar, 0.0]
n_exog(e::ExogSV) = 3
_unpack(e::ExogSV) = (e.γ, e.σz2, e.vbar, e.φv, e.τ, e.nϵ1, e.nϵ2, e.nϵ3)

const _sv_exog_ordering_docstr_extra ="""
Because lzhp and vp are time t+1 state variables, but lgp is not, the ordering
of variables here is significant. Specifically, at points in the code it is
assumed that lzhp comes first and lgp comes last. This is also true in the
constant volatiltiy version of the model where there is
"""

"""
Takes one step for exogenous process, given current levels and innovations

$_sv_exog_ordering_docstr_extra
"""
function exog_step(exog::ExogSV, lzh, v, ϵ1, ϵ2, ϵ3)
    γ, σz2, vbar, φv, τ, nϵ1, nϵ2, nϵ3 = _unpack(exog)
    lzhp = (1-2γ)*lzh .+ 0.5*((sqrt(v).*ϵ1) .- σz2*ϵ2)
    lgp = 0.0.*lzh .+ 0.5*sqrt(v).*ϵ1 .+ 0.5*σz2*ϵ2
    vp = (1-φv)*vbar + φv*v .+ τ*ϵ3

    lzhp, vp, lgp
end

# API methods on above function
exog_step(exog::ExogSV, current::Vector, ϵ::Vector) =
    exog_step(exog, current[1], current[2], ϵ[1], ϵ[2], ϵ[3])

exog_step(exog::ExogSV, current::Vector, ϵ::Matrix) =
    exog_step(exog, current[1], current[2], ϵ[:, 1], ϵ[:, 2], ϵ[:, 3])

# ------------------------------------------------------------------- #
# Stochastic Volatility Model
# ------------------------------------------------------------------- #
type ModelSV <: AbstractModel
    agent1::Agent
    agent2::Agent
    exog::ExogSV
    m_ls::NTuple{3,LinSpace{Float64}}
    grid::Matrix{Float64}
    grid_transpose::Matrix{Float64}
    lλ☆p_pol::Matrix{Float64}
    ϵ::Matrix{Float64}
    Π::Vector{Float64}
    lzhp::Matrix{Float64}
    vp::Matrix{Float64}
    lgp::Matrix{Float64}
    a1_itp::ScaledInterpolation
    b2_itp::ScaledInterpolation
    other::OtherItems

    """
    Given the paramters of the model as well as information on the coarseness
    of the grids, do initial setup work to define the model. This involves
    the following steps:

    - Construct instances of `Agent` to hold preference paramters for
      both agents
    - Construct an object to hold the paramters of the exogenous process
    - Define grids over state variables (lλ☆, lzh, v)
    - Define an initial guess for the policy of lλ☆p
    - Build quadrature weights and nodes to approximate the expectations
    - Precompute solutions the allocation in terms of lλ☆ and lzh (these
      decisions are given by the associated FOC and are independent of v)

    Note that all parameters match the notation from the paper and are settable
    one at a time using keyword argument syntax (`ModelSV(β2=0.99)`) changes
    only β2 relative to the default values. The only notational exception is
    that the standard deviation of the shock to `z2` is `σz2` here and
    `sqrt(v)` in the paper.

    Other arguments to this constructor are:

    - `lλ☆bds`: bound on grid for lλ☆. The grid is assumed to be symmetric,
    extending from `-lλ☆bds` to `lλ☆bds`.
    - `nX` for various `X`: the number of grid points for variable `X`
    - `σz2`: standard deviation of shock to `lz2` (see above)
    """
    function ModelSV(;ρ1::Real=-1.0, α1::Real=-9.0, β1::Real=0.98,
                      ω1::Real=0.1, σ1::Real=0.0,
                      ρ2::Real=-1.0, α2::Real=-9.0, β2::Real=0.98,
                      ω2::Real=0.1, σ2::Real=0.0,
                      γ::Real=0.10, vbar::Real=0.015^2, σz2::Real=0.015,
                      φv::Real=0.9, τ::Real=0.74e-5, lλ☆bds=4.6,
                      nλ☆::Int=451, nlzh::Int=13, nv::Int=13,
                      nϵ1::Int=5, nϵ2::Int=5, nϵ3::Int=5)

        # pack up agents and exog
        agent1 = Agent(ρ1, α1, β1, ω1, σ1)
        agent2 = Agent(ρ2, α2, β2, ω2, σ2)
        exog = ExogSV(γ, σz2, vbar, φv, τ, nϵ1, nϵ2, nϵ3)

        # Get linspaces that define grids  for lλ☆, zh, vh
        lλ☆_grid = linspace(-lλ☆bds, lλ☆bds, nλ☆)
        lλ☆_bds = (-lλ☆bds, lλ☆bds)

        # Get exogenous grid bounds from simulation of 5,000,000 periods
        lzh_bds, v_bds, _ = get_exog_bounds(exog)
        lzh_grid = linspace(lzh_bds[1], lzh_bds[2], nlzh)
        v_grid = linspace(v_bds[1], v_bds[2], nv)
        m_ls = (lλ☆_grid, lzh_grid, v_grid)

        # Create a grid
        grid = gridmake(collect(lλ☆_grid), collect(lzh_grid), collect(v_grid))
        grid_transpose = grid'

        # Initial guess at policy for lλ☆_{t+1} at each exog state in t+1 is
        # the time t level on the grid
        lλ☆_pol = repeat(grid_transpose[1, :], inner=[nϵ1*nϵ2*nϵ3, 1])

        # build exog grid for next period.  Both lgp, ξp are Neps × Ngrid
        ϵ, Π = qnwnorm([nϵ1, nϵ2, nϵ3], [0.0, 0.0, 0.0], eye(3))
        lzhp, vp, lgp = exog_step(exog, grid[:, 2]', grid[:, 3]',
                                  ϵ[:, 1], ϵ[:, 2], ϵ[:, 3])

        # incompletely initialize the ModelSV so we can pass it to routines that
        # construct the `OtherItems` for us
        m = new(agent1, agent2, exog, m_ls, grid, grid_transpose, lλ☆_pol,
                ϵ, Π, lzhp, vp, lgp)

        # Precompute solutions for a1 and b2 on a very fine grid for lλ☆ and
        # lzh. Then store coefficients of corresponding interpoland
        a1_itp, b2_itp = solve_for_a1b2(m, lλ☆_bds, lzh_bds)
        m.a1_itp = a1_itp
        m.b2_itp = b2_itp

        # Use interpolands for a1, b2 to find these variables on the grid we
        # will be solving the rest of the model on
        a1 = Float64[a1_itp[grid_transpose[1:2, i]...] for i=1:nλ☆*nlzh*nv]
        b2 = Float64[b2_itp[grid_transpose[1:2, i]...] for i=1:nλ☆*nlzh*nv]

        # now extract the rest of the allocation and package into the
        # `OtherItems` object, attching it to the model
        a2 = exp(grid[:, 2]) - a1
        b1 = 1.0 ./ exp(grid[:, 2]) - b2
        c1 = _agg_c(m.agent1, a1, b1)
        c2 = _agg_c(m.agent2, b2, a2)
        m.other = OtherItems(a1, b2, c1, c2)

        m
    end
end

"""
Extracts time t+1 levels of lzh, v, and lgp on the grid in state i at time t for
all time t+1 states

$_sv_exog_ordering_docstr_extra
"""
function get_exog_tomorrow(m::ModelSV, i::Int)
    lzhp_i = m.lzhp[:, i]
    lgp_i = m.lgp[:, i]
    vp_i = m.vp[:, i]
    lzhp_i, vp_i, lgp_i
end

"""
Extracts time t+1 levels of lzh, v, and lgp on the grid in state j at time t+1
and i at time t

$_sv_exog_ordering_docstr_extra
"""
function get_exog_tomorrow(m::ModelSV, i::Int, j::Int)
    lzhp_ij = m.lzhp[j, i]
    lgp_ij = m.lgp[j, i]
    vp_ij = m.vp[j, i]
    lzhp_ij, vp_ij, lgp_ij
end

jld_prefix(::ModelSV) = "curr_soln_cv"
jld_prefix(::Type{ModelSV}) = "curr_soln_cv"


text:
module PropertyFiles

using Dates: now

import Base: haskey, delete!, keys, values, pairs, sizehint!
import Base: keytype, valtype

export  Properties,
        setprop,
        getprop,
        getprop!

export  haskey,
        delete!,
        keys,
        values,
        pairs,
        sizehint!,
        keytype,
        valtype

export  store,
        load

struct Properties
    properties::Dict{String,String}
end

"""
    Properties()

Constructor for the `Properties` type.
"""
Properties() = Properties(Dict{String,String}())

include("overloads.jl")
include("store.jl")

"""
    setprop(prop::Properties, key::String, value::String)

Puts a property into the property-dictionary.
Only `String` keys and values can be used.
"""
function setprop(prop::Properties, key::String, value::String)
    prop.properties[key] = value
    nothing
end

"""
    setprop(prop::Properties, key::String, value)

Puts a property into the property-dictionary.
The `string()` constructor is used to
convert non-string values to string.
The property-dictionary is not returned.
"""
function setprop(prop::Properties, key::String, value)
    prop.properties[key] = string(value)
    nothing
end

"""
    getprop(prop::Properties, key::String, default::Union{String,Nothing}=nothing)

Return the value stored for the given key or the given default value, if there's no such key.
If the key does not exists and no default value is given, `nothing` will be returned.
"""
function getprop(prop::Properties, key::String, default::Union{String,Nothing}=nothing)
    get(prop.properties, key, default)
end

"""
    getprop!(prop::Properties, key::String, default::String)

Return the value stored for the given key.
If there's no such key, add the `key=>default` pair to the property-dictionary.
See `Base.get!`.
"""
function getprop!(prop::Properties, key::String, default::String)
    get!(prop.properties, key, default)
end

end # module


text:
abstract type Outbound end

function JSON.Writer.CompositeTypeWrapper(t::Outbound)
    fns = collect(fieldnames(typeof(t)))
    dels = Int[]
    for i = 1:length(fns)
        f = fns[i]
        if getfield(t, f) isa Missing
            push!(dels, i)
        end
    end
    deleteat!(fns, dels)
    JSON.Writer.CompositeTypeWrapper(t, Tuple(fns))
end

function JSON.lower(a::Outbound)
    if nfields(a) > 0
        JSON.Writer.CompositeTypeWrapper(a)
    else
        nothing
    end
end

function field_allows_missing(field::Expr)
    field.head == :(::) && field.args[2] isa Expr &&
    field.args[2].head == :curly && field.args[2].args[1] == :Union &&
    any(i -> i == :Missing, field.args[2].args)
end

function field_type(field::Expr, typename::String)
    if field.args[2] isa Expr && field.args[2].head == :curly && field.args[2].args[1] == :Union
        if length(field.args[2].args) == 3 && (field.args[2].args[2] == :Missing || field.args[2].args[3] == :Missing)
            return field.args[2].args[2] == :Missing ? field.args[2].args[3] : field.args[2].args[2]
        else
            # We return Any for now, which will lead to no type conversion
            return :Any
        end
    else
        return field.args[2]
    end
end

function get_kwsignature_for_field(field::Expr)
    fieldname = field.args[1]
    fieldtype = field.args[2]
    default_value = field_allows_missing(field) ? missing : :(error("You must provide a value for the $fieldname field."))

    return Expr(:kw, Expr(Symbol("::"), fieldname, fieldtype), default_value)
end

macro dict_readable(arg)
    tname = arg.args[2] isa Expr ? arg.args[2].args[1] : arg.args[2]
    count_real_fields = count(field -> !(field isa LineNumberNode), arg.args[3].args)
    ex = quote
        $((arg))

        $(count_real_fields > 0 ? :(
        function $tname(; $((get_kwsignature_for_field(field) for field in arg.args[3].args if !(field isa LineNumberNode))...))
            $tname($((field.args[1] for field in arg.args[3].args if !(field isa LineNumberNode))...))
        end
        ) : nothing)

        function $tname(dict::Dict)
        end
    end

    fex = :($((tname))())
    for field in arg.args[3].args
        if !(field isa LineNumberNode)
            fieldname = string(field.args[1])
            fieldtype = field_type(field, string(tname))
            if fieldtype isa Expr && fieldtype.head == :curly && fieldtype.args[2] != :Any
                f = :($(fieldtype.args[2]).(dict[$fieldname]))
            elseif fieldtype != :Any
                f = :($(fieldtype)(dict[$fieldname]))
            else
                f = :(dict[$fieldname])
            end
            if field_allows_missing(field)
                f = :(haskey(dict, $fieldname) ? $f : missing)
            end
            push!(fex.args, f)
        end
    end
    push!(ex.args[end].args[2].args, fex)
    return esc(ex)
end


text:
function formatinfo(s::Spec, ::Missing)
    size = 7  # ncodeunits("missing")
    return size + paddingsize(s, size), nothing
end

function formatfield(data::Vector{UInt8}, p::Int, s::Spec, ::Missing, ::Nothing)
    align = default(s.align, ALIGN_RIGHT)
    pw = paddingwidth(s, 7)
    p = padleft(data, p, s.fill, align, pw)
    p = @copy data p "missing"
    p = padright(data, p, s.fill, align, pw)
    return p
end