text:
using Test
using StaticArrays, IntervalSets, LinearAlgebra
using JET

import ClimaCore: slab, Domains, Meshes, Topologies, Spaces, Fields, Operators
import ClimaCore.Domains: Geometry

import ClimaCore.Operators: half, PlusHalf

@testset "PlusHalf" begin
    @test half + 0 == half
    @test half < half + 1
    @test half <= half + 1
    @test !(half > half + 1)
    @test !(half >= half + 1)
    @test half != half + 1
    @test half + half == 1
    @test half - half == 0
    @test half + 3 == 3 + half == PlusHalf(3)
    @test min(half, half + 3) == half
    @test max(half, half + 3) == half + 3

    @test collect(half:(2 + half)) == [half, 1 + half, 2 + half]

    @test_throws InexactError convert(Int, half)
    @test_throws InexactError convert(PlusHalf, 1)
    @test_throws InexactError convert(PlusHalf{Int}, 1)
end

@testset "Scalar Field FiniteDifferenceSpaces" begin
    for FT in (Float32, Float64)
        domain = Domains.IntervalDomain(
            Geometry.ZPoint{FT}(0.0),
            Geometry.ZPoint{FT}(pi);
            boundary_tags = (:left, :right),
        )
        @test eltype(domain) === Geometry.ZPoint{FT}

        mesh = Meshes.IntervalMesh(domain; nelems = 16)

        center_space = Spaces.CenterFiniteDifferenceSpace(mesh)
        face_space = Spaces.FaceFiniteDifferenceSpace(center_space)

        @test sum(ones(FT, center_space)) ≈ pi
        @test sum(ones(FT, face_space)) ≈ pi

        centers = getproperty(Fields.coordinate_field(center_space), :z)
        @test sum(sin.(centers)) ≈ FT(2.0) atol = 1e-2

        faces = getproperty(Fields.coordinate_field(face_space), :z)
        @test sum(sin.(faces)) ≈ FT(2.0) atol = 1e-2

        ∇ᶜ = Operators.GradientF2C()
        ∂sin = Geometry.WVector.(∇ᶜ.(sin.(faces)))
        @test ∂sin ≈ Geometry.WVector.(cos.(centers)) atol = 1e-2

        divᶜ = Operators.DivergenceF2C()
        ∂sin = divᶜ.(Geometry.WVector.(sin.(faces)))
        @test ∂sin ≈ cos.(centers) atol = 1e-2

        # Center -> Face operator
        # first order convergence at boundaries
        ∇ᶠ = Operators.GradientC2F(
            left = Operators.SetValue(FT(0)),
            right = Operators.SetValue(FT(pi)),
        )
        ∂z = Geometry.WVector.(∇ᶠ.(centers))
        @test ∂z ≈ Geometry.WVector.(ones(FT, face_space)) rtol = 10 * eps(FT)

        ∇ᶠ = Operators.GradientC2F(
            left = Operators.SetValue(FT(1)),
            right = Operators.SetValue(FT(-1)),
        )
        ∂cos = Geometry.WVector.(∇ᶠ.(cos.(centers)))
        @test ∂cos ≈ Geometry.WVector.(.-sin.(faces)) atol = 1e-1

        ∇ᶠ = Operators.GradientC2F(
            left = Operators.SetGradient(Geometry.WVector(FT(0))),
            right = Operators.SetGradient(Geometry.WVector(FT(0))),
        )
        ∂cos = Geometry.WVector.(∇ᶠ.(cos.(centers)))
        @test ∂cos ≈ Geometry.WVector.(.-sin.(faces)) atol = 1e-2

        # test that broadcasting into incorrect field space throws an error
        empty_centers = zeros(FT, center_space)
        @test_throws Exception empty_centers .= ∇ᶠ.(cos.(centers))
    end
end

@testset "Test composed stencils" begin
    for FT in (Float32, Float64)
        domain = Domains.IntervalDomain(
            Geometry.ZPoint{FT}(0.0),
            Geometry.ZPoint{FT}(pi);
            boundary_tags = (:left, :right),
        )
        @test eltype(domain) === Geometry.ZPoint{FT}

        mesh = Meshes.IntervalMesh(domain; nelems = 16)

        center_space = Spaces.CenterFiniteDifferenceSpace(mesh)
        face_space = Spaces.FaceFiniteDifferenceSpace(center_space)

        centers = getproperty(Fields.coordinate_field(center_space), :z)
        w = ones(FT, face_space)
        θ = sin.(centers)

        # 1) we set boundaries on the 2nd operator
        I = Operators.InterpolateC2F()
        ∂ = Operators.GradientF2C(
            left = Operators.SetValue(FT(0)),
            right = Operators.SetValue(FT(0)),
        )

        ∂sin = Geometry.WVector.(∂.(w .* I.(θ)))
        @test ∂sin ≈ Geometry.WVector.(cos.(centers)) atol = 1e-2

        # can't define Neumann conditions on GradientF2C
        ∂ = Operators.GradientF2C(
            left = Operators.SetGradient(FT(1)),
            right = Operators.SetGradient(FT(-1)),
        )

        @test_throws Exception ∂.(w .* I.(θ))

        # 2) we set boundaries on the 1st operator
        I = Operators.InterpolateC2F(
            left = Operators.SetValue(FT(0)),
            right = Operators.SetValue(FT(0)),
        )
        ∂ = Operators.GradientF2C()

        ∂sin = Geometry.WVector.(∂.(w .* I.(θ)))
        @test ∂sin ≈ Geometry.WVector.(cos.(centers)) atol = 1e-2

        I = Operators.InterpolateC2F(
            left = Operators.SetGradient(Geometry.WVector(FT(1))),
            right = Operators.SetGradient(Geometry.WVector(FT(-1))),
        )
        ∂ = Operators.GradientF2C()

        ∂sin = Geometry.WVector.(∂.(w .* I.(θ)))
        @test ∂sin ≈ Geometry.WVector.(cos.(centers)) atol = 1e-2

        # 3) we set boundaries on both: 2nd should take precedence
        I = Operators.InterpolateC2F(
            left = Operators.SetValue(FT(NaN)),
            right = Operators.SetValue(FT(NaN)),
        )
        ∂ = Operators.GradientF2C(
            left = Operators.SetValue(FT(0)),
            right = Operators.SetValue(FT(0)),
        )

        ∂sin = Geometry.WVector.(∂.(w .* I.(θ)))
        @test ∂sin ≈ Geometry.WVector.(cos.(centers)) atol = 1e-2

        # test that broadcasting into incorrect field space throws an error
        empty_faces = zeros(FT, face_space)
        @test_throws Exception empty_faces .= ∂.(w .* I.(θ))

        # 5) we set boundaries on neither
        I = Operators.InterpolateC2F()
        ∂ = Operators.GradientF2C()

        # TODO: should we throw something else?
        @test_throws BoundsError ∂.(w .* I.(θ))
    end
end

@testset "Composite Field FiniteDifferenceSpaces" begin
    for FT in (Float32, Float64)
        domain = Domains.IntervalDomain(
            Geometry.ZPoint{FT}(0.0),
            Geometry.ZPoint{FT}(pi);
            boundary_tags = (:left, :right),
        )

        @test eltype(domain) === Geometry.ZPoint{FT}
        mesh = Meshes.IntervalMesh(domain; nelems = 16)

        center_space = Spaces.CenterFiniteDifferenceSpace(mesh)
        face_space = Spaces.FaceFiniteDifferenceSpace(center_space)

        FieldType = NamedTuple{(:a, :b), Tuple{FT, FT}}

        center_sum = sum(ones(FieldType, center_space))
        @test center_sum isa FieldType
        @test center_sum.a ≈ FT(pi)
        @test center_sum.b ≈ FT(pi)

        face_sum = sum(ones(FieldType, face_space))
        @test face_sum isa FieldType
        @test face_sum.a ≈ FT(pi)
        @test face_sum.b ≈ FT(pi)
    end
end

"""
    convergence_rate(err, Δh)

Estimate convergence rate given vectors `err` and `Δh`

    err = C Δh^p+ H.O.T
    err_k ≈ C Δh_k^p
    err_k/err_m ≈ Δh_k^p/Δh_m^p
    log(err_k/err_m) ≈ log((Δh_k/Δh_m)^p)
    log(err_k/err_m) ≈ p*log(Δh_k/Δh_m)
    log(err_k/err_m)/log(Δh_k/Δh_m) ≈ p

"""
convergence_rate(err, Δh) =
    [log(err[i] / err[i - 1]) / log(Δh[i] / Δh[i - 1]) for i in 2:length(Δh)]

@testset "Face -> Center interpolation (uniform and stretched)" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    n_elems_seq = 2 .^ (5, 6, 7, 8)
    stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(0.5))
    for (i, stretch_fn) in enumerate(stretch_fns)
        err, Δh = zeros(FT, length(n_elems_seq)), zeros(FT, length(n_elems_seq))

        for (k, n) in enumerate(n_elems_seq)
            interval = Geometry.ZPoint(a) .. Geometry.ZPoint(b)
            domain = Domains.IntervalDomain(
                interval;
                boundary_tags = (:left, :right),
            )
            mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)

            cs = Spaces.CenterFiniteDifferenceSpace(mesh)
            fs = Spaces.FaceFiniteDifferenceSpace(cs)

            cent_field_exact = zeros(FT, cs)
            cent_field = zeros(FT, cs)
            face_field = zeros(FT, fs)

            centers = Fields.coordinate_field(cs)
            faces = Fields.coordinate_field(fs)

            face_field .= sin.(3π .* faces.z)
            cent_field_exact .= sin.(3π .* centers.z)
            operator = Operators.InterpolateF2C()
            cent_field .= operator.(face_field)

            Δh[k] = cs.face_local_geometry.J[1]
            err[k] =
                norm(parent(cent_field) .- parent(cent_field_exact)) /
                length(parent(cent_field_exact))
        end
        conv = convergence_rate(err, Δh)
        # conv should be approximately 2 for second order-accurate stencil.
        @test 1.5 ≤ conv[1] ≤ 3
        @test 1.5 ≤ conv[3] ≤ 3
        if i == 1
            @test conv[1] ≤ conv[2] ≤ conv[3]
        end
        @test err[3] ≤ err[2] ≤ err[1] ≤ 1e-2
    end
end

@testset "Center -> Face interpolation (uniform and stretched)" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    n_elems_seq = 2 .^ (5, 6, 7, 8)
    stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(0.5))
    for (i, stretch_fn) in enumerate(stretch_fns)
        err, Δh = zeros(FT, length(n_elems_seq)), zeros(FT, length(n_elems_seq))

        for (k, n) in enumerate(n_elems_seq)
            interval = Geometry.ZPoint(a) .. Geometry.ZPoint(b)
            domain = Domains.IntervalDomain(
                interval;
                boundary_tags = (:left, :right),
            )
            mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)

            cs = Spaces.CenterFiniteDifferenceSpace(mesh)
            fs = Spaces.FaceFiniteDifferenceSpace(cs)

            face_field_exact = zeros(FT, fs)
            cent_field = zeros(FT, cs)
            face_field = zeros(FT, fs)

            centers = Fields.coordinate_field(cs)
            faces = Fields.coordinate_field(fs)

            cent_field .= sin.(3π .* centers.z)
            face_field_exact .= sin.(3π .* faces.z)

            operator = Operators.InterpolateC2F(
                left = Operators.SetValue(0.0),
                right = Operators.SetValue(0.0),
            )
            face_field .= operator.(cent_field)

            Δh[k] = cs.face_local_geometry.J[1]
            err[k] =
                norm(parent(face_field) .- parent(face_field_exact)) /
                length(parent(face_field_exact))
        end
        conv = convergence_rate(err, Δh)
        # conv should be approximately 2 for second order-accurate stencil.
        @test 1.5 ≤ conv[1] ≤ 3
        @test 1.5 ≤ conv[3] ≤ 3
        if i == 1
            @test conv[1] ≤ conv[2] ≤ conv[3]
        end
        @test err[3] ≤ err[2] ≤ err[1] ≤ 1e-2
    end
end

@testset "∂ Center -> Face interpolation (uniform and stretched)" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    n_elems_seq = 2 .^ (5, 6, 7, 8)
    stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(0.5))
    for (i, stretch_fn) in enumerate(stretch_fns)
        err, Δh = zeros(FT, length(n_elems_seq)), zeros(FT, length(n_elems_seq))
        for (k, n) in enumerate(n_elems_seq)
            interval = Geometry.ZPoint(a) .. Geometry.ZPoint(b)
            domain = Domains.IntervalDomain(
                interval;
                boundary_tags = (:left, :right),
            )
            mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)

            cs = Spaces.CenterFiniteDifferenceSpace(mesh)
            fs = Spaces.FaceFiniteDifferenceSpace(cs)

            face_field_exact = Geometry.Covariant3Vector.(zeros(FT, fs))
            cent_field = zeros(FT, cs)
            face_field = Geometry.Covariant3Vector.(zeros(FT, fs))

            centers = Fields.coordinate_field(cs)
            faces = Fields.coordinate_field(fs)

            cent_field .= sin.(3π .* centers.z)
            face_field_exact .=
                Geometry.CovariantVector.(
                    Geometry.WVector.(3π .* cos.(3π .* faces.z)),
                )

            operator = Operators.GradientC2F(
                left = Operators.SetGradient(Geometry.WVector(3π)),
                right = Operators.SetGradient(Geometry.WVector(-3π)),
            )

            face_field .= operator.(cent_field)

            Δh[k] = cs.face_local_geometry.J[1]
            err[k] = norm(face_field .- face_field_exact)
        end
        conv = convergence_rate(err, Δh)
        # conv should be approximately 2 for second order-accurate stencil.
        @test err[3] ≤ err[2] ≤ err[1] ≤ 0.1
        @test conv[1] ≈ 2 atol = 0.1
        @test conv[2] ≈ 2 atol = 0.1
        @test conv[3] ≈ 2 atol = 0.1
        @test conv[1] ≤ conv[2] ≤ conv[3]
    end
end

@testset "∂ Face -> Center interpolation (uniform and stretched)" begin
    FT = Float64
    a, b = FT(0.0), FT(1.0)
    n_elems_seq = 2 .^ (5, 6, 7, 8)
    stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(0.5))
    for (i, stretch_fn) in enumerate(stretch_fns)
        err, Δh = zeros(FT, length(n_elems_seq)), zeros(FT, length(n_elems_seq))
        for (k, n) in enumerate(n_elems_seq)
            interval = Geometry.ZPoint(a) .. Geometry.ZPoint(b)
            domain = Domains.IntervalDomain(
                interval;
                boundary_tags = (:left, :right),
            )
            mesh = Meshes.IntervalMesh(domain, stretch_fn, nelems = n)

            cs = Spaces.CenterFiniteDifferenceSpace(mesh)
            fs = Spaces.FaceFiniteDifferenceSpace(cs)

            cent_field_exact = Geometry.Covariant3Vector.(zeros(FT, cs))
            cent_field = Geometry.Covariant3Vector.(zeros(FT, cs))
            face_field = zeros(FT, fs)

            centers = Fields.coordinate_field(cs)
            faces = Fields.coordinate_field(fs)

            face_field .= sin.(3π .* faces.z)
            cent_field_exact .=
                Geometry.CovariantVector.(
                    Geometry.WVector.(3π .* cos.(3π .* centers.z)),
                )

            operator = Operators.GradientF2C()

            cent_field .= operator.(face_field)

            Δh[k] = cs.face_local_geometry.J[1]
            err[k] = norm(cent_field .- cent_field_exact)
        end
        conv = convergence_rate(err, Δh)
        # conv should be approximately 2 for second order-accurate stencil.
        @test err[3] ≤ err[2] ≤ err[1] ≤ 0.1
        @test conv[1] ≈ 2 atol = 0.1
        @test conv[2] ≈ 2 atol = 0.1
        @test conv[3] ≈ 2 atol = 0.1
        @test conv[1] ≤ conv[2] ≤ conv[3]
    end
end

@testset "∂ Center -> Face and ∂ Face-> Center (uniform)" begin
    FT = Float64
    n_elems_seq = 2 .^ (5, 6, 7, 8)

    err_grad_sin_c = zeros(FT, length(n_elems_seq))
    err_div_sin_c = zeros(FT, length(n_elems_seq))
    err_grad_z_f = zeros(FT, length(n_elems_seq))
    err_grad_cos_f1 = zeros(FT, length(n_elems_seq))
    err_grad_cos_f2 = zeros(FT, length(n_elems_seq))
    err_div_sin_f = zeros(FT, length(n_elems_seq))
    err_div_cos_f = zeros(FT, length(n_elems_seq))
    err_curl_sin_f = zeros(FT, length(n_elems_seq))
    Δh = zeros(FT, length(n_elems_seq))

    for (k, n) in enumerate(n_elems_seq)
        domain = Domains.IntervalDomain(
            Geometry.ZPoint{FT}(0.0),
            Geometry.ZPoint{FT}(pi);
            boundary_tags = (:left, :right),
        )
        mesh = Meshes.IntervalMesh(domain; nelems = n)

        cs = Spaces.CenterFiniteDifferenceSpace(mesh)
        fs = Spaces.FaceFiniteDifferenceSpace(cs)

        centers = getproperty(Fields.coordinate_field(cs), :z)
        faces = getproperty(Fields.coordinate_field(fs), :z)

        # Face -> Center operators:
        # GradientF2C
        # f(z) = sin(z)
        ∇ᶜ = Operators.GradientF2C()
        gradsinᶜ = Geometry.WVector.(∇ᶜ.(sin.(faces)))

        # DivergenceF2C
        # f(z) = sin(z)
        divᶜ = Operators.DivergenceF2C()
        divsinᶜ = divᶜ.(Geometry.WVector.(sin.(faces)))

        # Center -> Face operators:
        # GradientC2F, SetValue
        # f(z) = z
        ∇ᶠ⁰ = Operators.GradientC2F(
            left = Operators.SetValue(FT(0)),
            right = Operators.SetValue(FT(pi)),
        )
        ∂zᶠ = Geometry.WVector.(∇ᶠ⁰.(centers))

        # GradientC2F, SetValue
        # f(z) = cos(z)
        ∇ᶠ¹ = Operators.GradientC2F(
            left = Operators.SetValue(FT(1)),
            right = Operators.SetValue(FT(-1)),
        )
        gradcosᶠ¹ = Geometry.WVector.(∇ᶠ¹.(cos.(centers)))

        # GradientC2F, SetGradient
        # f(z) = cos(z)
        ∇ᶠ² = Operators.GradientC2F(
            left = Operators.SetGradient(Geometry.WVector(FT(0))),
            right = Operators.SetGradient(Geometry.WVector(FT(0))),
        )
        gradcosᶠ² = Geometry.WVector.(∇ᶠ².(cos.(centers)))

        # DivergenceC2F, SetValue
        # f(z) = sin(z)
        divᶠ⁰ = Operators.DivergenceC2F(
            left = Operators.SetValue(Geometry.WVector(zero(FT))),
            right = Operators.SetValue(Geometry.WVector(zero(FT))),
        )
        divsinᶠ = divᶠ⁰.(Geometry.WVector.(sin.(centers)))

        # DivergenceC2F, SetDivergence
        # f(z) = cos(z)
        divᶠ¹ = Operators.DivergenceC2F(
            left = Operators.SetDivergence(FT(0)),
            right = Operators.SetDivergence(FT(0)),
        )
        divcosᶠ = divᶠ¹.(Geometry.WVector.(cos.(centers)))

        curlᶠ = Operators.CurlC2F(
            left = Operators.SetValue(Geometry.Covariant1Vector(zero(FT))),
            right = Operators.SetValue(Geometry.Covariant1Vector(zero(FT))),
        )
        curlsinᶠ = curlᶠ.(Geometry.Covariant1Vector.(sin.(centers)))

        Δh[k] = cs.face_local_geometry.J[1]
        # Errors
        err_grad_sin_c[k] = norm(gradsinᶜ .- Geometry.WVector.(cos.(centers)))
        err_div_sin_c[k] = norm(divsinᶜ .- cos.(centers))
        err_grad_z_f[k] = norm(∂zᶠ .- Geometry.WVector.(ones(FT, fs)))
        err_grad_cos_f1[k] = norm(gradcosᶠ¹ .- Geometry.WVector.(.-sin.(faces)))
        err_grad_cos_f2[k] = norm(gradcosᶠ² .- Geometry.WVector.(.-sin.(faces)))
        err_div_sin_f[k] =
            norm(divsinᶠ .- (Geometry.WVector.(cos.(faces))).components.data.:1)
        err_div_cos_f[k] = norm(
            divcosᶠ .- (Geometry.WVector.(.-sin.(faces))).components.data.:1,
        )
        err_curl_sin_f[k] =
            norm(curlsinᶠ .- Geometry.Contravariant2Vector.(cos.(faces)))
    end

    # GradientF2C conv, with f(z) = sin(z)
    conv_grad_sin_c = convergence_rate(err_grad_sin_c, Δh)
    # DivergenceF2C conv, with f(z) = sin(z)
    conv_div_sin_c = convergence_rate(err_div_sin_c, Δh)
    # GradientC2F conv, with f(z) = z, SetValue
    conv_grad_z = convergence_rate(err_grad_z_f, Δh)
    # GradientC2F conv, with f(z) = cos(z), SetValue
    conv_grad_cos_f1 = convergence_rate(err_grad_cos_f1, Δh)
    # GradientC2F conv, with f(z) = cos(z), SetGradient
    conv_grad_cos_f2 = convergence_rate(err_grad_cos_f2, Δh)
    # DivergenceC2F conv, with f(z) = sin(z), SetValue
    conv_div_sin_f = convergence_rate(err_div_sin_f, Δh)
    # DivergenceC2F conv, with f(z) = cos(z), SetDivergence
    conv_div_cos_f = convergence_rate(err_div_cos_f, Δh)
    # CurlC2F with f(z) = sin(z), SetValue
    conv_curl_sin_f = convergence_rate(err_curl_sin_f, Δh)

    # GradientF2C conv, with f(z) = sin(z)
    @test err_grad_sin_c[3] ≤ err_grad_sin_c[2] ≤ err_grad_sin_c[1] ≤ 0.1
    @test conv_grad_sin_c[1] ≈ 2 atol = 0.1
    @test conv_grad_sin_c[2] ≈ 2 atol = 0.1
    @test conv_grad_sin_c[3] ≈ 2 atol = 0.1
    @test conv_grad_sin_c[1] ≤ conv_grad_sin_c[2] ≤ conv_grad_sin_c[3]

    # DivergenceF2C conv, with f(z) = sin(z)
    @test err_div_sin_c[3] ≤ err_div_sin_c[2] ≤ err_div_sin_c[1] ≤ 0.1
    @test conv_div_sin_c[1] ≈ 2 atol = 0.1
    @test conv_div_sin_c[2] ≈ 2 atol = 0.1
    @test conv_div_sin_c[3] ≈ 2 atol = 0.1
    @test conv_div_sin_c[1] ≤ conv_div_sin_c[2] ≤ conv_div_sin_c[3]

    # GradientC2F conv, with f(z) = z, SetValue
    @test norm(err_grad_z_f) ≤ 200 * eps(FT)
    # Convergence rate for this case is noisy because error very small

    # GradientC2F conv, with f(z) = cos(z), SetValue
    @test err_grad_cos_f1[3] ≤ err_grad_cos_f1[2] ≤ err_grad_cos_f1[1] ≤ 0.1
    @test conv_grad_cos_f1[1] ≈ 1.5 atol = 0.1
    @test conv_grad_cos_f1[2] ≈ 1.5 atol = 0.1
    @test conv_grad_cos_f1[3] ≈ 1.5 atol = 0.1
    # @test conv_grad_cos_f1[1] ≤ conv_grad_cos_f1[2] ≤ conv_grad_cos_f1[3]

    # GradientC2F conv, with f(z) = cos(z), SetGradient
    @test err_grad_cos_f2[3] ≤ err_grad_cos_f2[2] ≤ err_grad_cos_f2[1] ≤ 0.1
    @test conv_grad_cos_f2[1] ≈ 2 atol = 0.1
    @test conv_grad_cos_f2[2] ≈ 2 atol = 0.1
    @test conv_grad_cos_f2[3] ≈ 2 atol = 0.1
    @test conv_grad_cos_f2[1] ≤ conv_grad_cos_f2[2] ≤ conv_grad_cos_f2[3]

    # DivergenceC2F conv, with f(z) = sin(z), SetValue
    @test err_div_sin_f[3] ≤ err_div_sin_f[2] ≤ err_div_sin_f[1] ≤ 0.1
    @test conv_div_sin_f[1] ≈ 2 atol = 0.1
    @test conv_div_sin_f[2] ≈ 2 atol = 0.1
    @test conv_div_sin_f[3] ≈ 2 atol = 0.1
    @test conv_div_sin_f[1] ≤ conv_div_sin_f[2] ≤ conv_div_sin_f[3]

    # DivergenceC2F conv, with f(z) = cos(z), SetDivergence
    @test err_div_cos_f[3] ≤ err_div_cos_f[2] ≤ err_div_cos_f[1] ≤ 0.1
    @test conv_div_cos_f[1] ≈ 2 atol = 0.1
    @test conv_div_cos_f[2] ≈ 2 atol = 0.1
    @test conv_div_cos_f[3] ≈ 2 atol = 0.1
    @test conv_div_cos_f[1] ≤ conv_div_cos_f[2] ≤ conv_div_cos_f[3]

    # CurlC2F with f(z) = sin(z), SetValue
    @test err_curl_sin_f[3] ≤ err_curl_sin_f[2] ≤ err_curl_sin_f[1] ≤ 0.1
    @test conv_curl_sin_f[1] ≈ 2 atol = 0.1
    @test conv_curl_sin_f[2] ≈ 2 atol = 0.1
    @test conv_curl_sin_f[3] ≈ 2 atol = 0.1
    @test conv_curl_sin_f[1] ≤ conv_curl_sin_f[2] ≤ conv_curl_sin_f[3]
end

@testset "Biased interpolation" begin
    FT = Float64
    n_elems = 10

    domain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(0.0),
        Geometry.ZPoint{FT}(pi);
        boundary_tags = (:bottom, :top),
    )
    mesh = Meshes.IntervalMesh(domain; nelems = n_elems)

    cs = Spaces.CenterFiniteDifferenceSpace(mesh)
    fs = Spaces.FaceFiniteDifferenceSpace(cs)

    zc = getproperty(Fields.coordinate_field(cs), :z)
    zf = getproperty(Fields.coordinate_field(fs), :z)

    function field_wrapper(space, nt::NamedTuple)
        cmv(z) = nt
        return cmv.(Fields.coordinate_field(space))
    end

    field_vars() = (; y = FT(0))

    cfield = field_wrapper(cs, field_vars())
    ffield = field_wrapper(fs, field_vars())

    cy = cfield.y
    fy = ffield.y

    cyp = parent(cy)
    fyp = parent(fy)

    # C2F biased operators
    LBC2F = Operators.LeftBiasedC2F(; bottom = Operators.SetValue(10))
    @. cy = cos(zc)
    @. fy = LBC2F(cy)
    fy_ref = [FT(10), [cyp[i] for i in 1:length(cyp)]...]
    @test all(fy_ref .== fyp)

    RBC2F = Operators.RightBiasedC2F(; top = Operators.SetValue(10))
    @. cy = cos(zc)
    @. fy = RBC2F(cy)
    fy_ref = [[cyp[i] for i in 1:length(cyp)]..., FT(10)]
    @test all(fy_ref .== fyp)

    # F2C biased operators
    LBF2C = Operators.LeftBiasedF2C(; bottom = Operators.SetValue(10))
    @. cy = cos(zc)
    @. cy = LBF2C(fy)
    cy_ref = [i == 1 ? FT(10) : fyp[i] for i in 1:length(cyp)]
    @test all(cy_ref .== cyp)

    RBF2C = Operators.RightBiasedF2C(; top = Operators.SetValue(10))
    @. cy = cos(zc)
    @. cy = RBF2C(fy)
    cy_ref = [i == length(cyp) ? FT(10) : fyp[i + 1] for i in 1:length(cyp)]
    @test all(cy_ref .== cyp)
end


text:
using Base.Test
using BCBI_base

@testset "Test Adding and Precompiling Packages" begin
    println("--------------------------------")
    println("Testin Base Pkgs")
    println("--------------------------------")
    @test length(BCBI_base.add(BCBI_base.base_pkgs)) == 0
    println("--------------------------------")
    println("Testing Plotting Pkgs")
    println("--------------------------------")
    @test length(BCBI_base.add(BCBI_base.plotting_pkgs)) == 0
    println("--------------------------------")
    println("Testing Datasets Pkgs")
    println("--------------------------------")
    @test length(BCBI_base.add(BCBI_base.datasets_pkg)) == 0
    println("--------------------------------")
    println("Testing Python/R Pkgs")
    println("--------------------------------")
    @test length(BCBI_base.add(BCBI_base.external_pkgs)) == 0
    println("--------------------------------")
    println("Testing Clone Pkgs")
    println("--------------------------------")
    @test length(BCBI_base.clone(BCBI_base.clone_pkgs)) == 0
    println("--------------------------------")
    println("Testing Checkout Pkgs")
    println("--------------------------------")
    @test length(BCBI_base.checkout(BCBI_base.checkout_pkgs)) == 0
    println("--------------------------------")
    println("Precompiling Installed Pkgs")
    println("--------------------------------")
    @test length(using_all()) == 0
    BCBI_base.clean_up()
    @test check_installed() == 0
end

ext: jl


text:
#=Copyright (c) 2017 Gabriel Goh

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.=#

# f(x) = ((1 + ((x[1]+x[2]+1)^2)*(19 - 14*x[1]+3*x[1]^2-14*x[2]+6*x[1]*x[2]+3*x[2]^2))*
#         (30 + ((2*x[1] - 3*x[2])^2)*(18 - 32*x[1] +
#             12*x[1]^2 + 48*x[2] - 36*x[1]*x[2] + 27*x[2]^2))  )
# ∇f(x) = ForwardDiff.gradient(f,x)

# Beale
# f(x) = ( (1.5 - x[1] + x[1]*x[2])^2
#      + (2.25 - x[1] + x[1]*x[2]^2)^2
#      + (2.625 - x[1] + x[1]*x[2]^3)^2 )
# ∇f(x) = ForwardDiff.gradient(f,x)

# Rosenbrock
# f(x) = 100*(x[2] - x[1]^2)^2 + (x[1] - 1)^2
# ∇f(x) = ForwardDiff.gradient(f,x)

# Monkey Saddle
# f(x) = x[1]^3 - 3*x[1]*x[2]^2
# ∇f(x) = [3*x[1] - 3*x[2]^2; -6*x[1]*x[2]]

A = diagm([1; 1])
f(x) = 0.5*vecdot(x,A*x)
∇f(x) = A*x

α = 1.5; β = 0.51

w = [1;1]; z = zeros(2)

N = 2000
fW = zeros(N)
wp = w
W = zeros(2,N)
for i = 1:N
    if i > 2
        #if vecdot(∇f(w), w - wp) > 0; println(i); z = zeros(2); end
    end
    W[:,i] = vecdot(∇f(w), w - wp)
    wp = w
    z = β*z + ∇f(w)
    w = w - α*(∇f(w) + β*z)
    println(f(w))
    if abs(f(w)) < 100; fW[i] = f(w); else; fW[i] = NaN; end
    #W[:,i] = w
end

plot(fW)

max_stars_repo_name: sbeleidy/post--momentum

max_issues_repo_name: sbeleidy/post--momentum

max_forks_repo_name: sbeleidy/post--momentum


text:
module Simulation

using Statistics, Random, Distributions

using ..Graphs
using ..Graphs.Utility
using ..Graphs.Intervals

# ------------------------------------------------------------------------
# types

"""
	struct Rates
		snp :: Float64
		hgt :: Float64
		del :: Float64
		inv :: Float64
	end

Store the rates of evolution of mutation `snp`, recombination `hgt`, deletion `del`, and inversion `inv`.
"""
struct Rates
    snp :: Float64
    hgt :: Float64
    del :: Float64
    inv :: Float64
end
Rates(;snp=0, hgt=0, del=0, inv=0) = Rates(snp, hgt, del, inv)

"""
	struct Params
		N    :: Int
		L    :: Int
		σ    :: Int
		rate :: Rates
	end

Store all parameters of a single recombinative Wright-Fisher model.
`N` is the population size.
`L` is the expected genome size of all descendants.
`σ` is the variance of genome size of all descendants.
`rate` is the various rate of evolutionary processes.
"""
struct Params
    N    :: Int
    L    :: Int
    σ    :: Int
    rate :: Rates
end

Params(; N=100, L=Int(1e6), σ=Int(1e5), snp=0, hgt=0, del=0, inv=0) = Params(N, L, σ, Rates(snp, hgt, del, inv))

# bitpacked: 30 bytes(ancestor) | 30 bytes (location) | 3 bytes (mutation) | 1 byte strand
"""
	Sequence = Array{UInt64,1}

A bitpacked array of sequence state. Each UInt64 bits are interpreted as
  > 30 bytes(ancestor) | 30 bytes (location) | 3 bytes (mutation) | 1 byte strand
"""
Sequence = Array{UInt64,1}
# states of mutuation
# 0** -> no mutation, ** ignored
# 1** -> mutation given by **
const shift = (
    ancestor = 34,
    location = 4,
    mutation = 1,
    strand   = 0,
)

mask(n::Int) = (64 ≤ n || n ≤ 0) ? 0 : UInt64((1 << n)-1)

const unpack = (
    ancestor = (x::UInt64) -> (x>>shift.ancestor) & mask(30),
    location = (x::UInt64) -> (x>>shift.location) & mask(30),
    mutation = (x::UInt64) -> (x>>shift.mutation) & mask(2),
    strand   = (x::UInt64) -> x & mask(1),

    # helpers
    snp = (x::UInt64) -> (x >> 3) & 1
)

ancestor(n::Int,L::Int)   = [(UInt64(n) << shift.ancestor) | (UInt64(l) << shift.location) for l in 1:L]
population(L::Array{Int}) = [ ancestor(n,len) for (n,len) in enumerate(L)]

# ------------------------------------------------------------------------
# evolution operators

"""
	mutate!(s::Sequence, at::Int)

Apply a random mutation to sequence `s` at locus `at`.
"""
mutate!(s::Sequence, at::Int) = s[at] |= ((rand(0:3) << 1) | 8)
"""
	delete!(s::Sequence, from::Int, to::Int)

Delete the interval `from`:`to` from sequence `s`.
"""
delete!(s::Sequence, from::Int, to::Int) = splice!(s, from:to, [])
"""
	insert!(acceptor::Sequence, donor::Sequence, at::Int)

Insert sequence `donor` into `acceptor` at locus `at`.
"""
insert!(acceptor::Sequence, donor::Sequence, at::Int) = splice!(acceptor, at:at-1, donor)
"""
	invert!(s::Sequence, from::Int, to::Int)

Replace the interval `from`:`to` of sequence `s` with its reverse complement.
"""
function invert!(s::Sequence, from::Int, to::Int)
    s[from:to] .⊻= 1        # complement
    reverse!(s, from, to)   # reverse
end

"""
	model(param::Params)

Return an evolution function based upon parameters `param`.
"""
function model(param::Params)
    parent = Array{Int}(undef, param.N)
    indel = Normal(param.L, param.σ)

    int(x) = Int(round(x))

    return function(population::Array{Sequence}, offspring::Array{Sequence})
        rand!(parent, 1:param.N)

        for n in 1:param.N
            copy!(offspring[n], population[parent[n]])

            # test to mutate
            L = length(offspring[n])
            μ = rand(Poisson(L*param.rate.snp))
            if μ > 0
                for locus in rand(1:L, μ)
                    mutate!(offspring[n], locus)
                end
            end

            # test to pick up donation
            if rand() ≤ param.rate.hgt
                donor = population[rand(1:param.N)]

                from = rand(1:length(donor))
                area = int(rand(indel))
                if area > length(offspring[n])
                    d  = area - length(offspring[n])
                    to = (from + d - 1) % length(donor) + 1

                    if to > from
                        donor = donor[from:to]
                    else
                        donor = vcat(donor[from:length(donor)], donor[1:to])
                    end

                    site = rand(1:length(offspring[n]))
                    insert!(offspring[n], donor, site)
                end
            end

            # test to delete
            if rand() ≤ param.rate.del
                from = rand(1:length(offspring[n]))
                area = int(rand(indel))
                if area < length(offspring[n])
                    d  = length(offspring[n]) - area
                    to = (from + d - 1) % length(offspring[n]) + 1

                    if to > from
                        delete!(offspring[n], from, to)
                    else
                        delete!(offspring[n], from, length(offspring[n]))
                        delete!(offspring[n], 1, to)
                    end
                end
            end

            # test to invert
            if rand() ≤ param.rate.inv
                from = rand(1:length(offspring[n]))
                area = int(rand(indel))
                d  = abs(area - length(offspring[n]))
                to = (from + d - 1) % length(offspring[n]) + 1

                if to > from
                    invert!(offspring[n], from, to)
                else
                    invert!(offspring[n], from, length(offspring[n]))
                    invert!(offspring[n], 1, to)
                end
            end
        end
    end
end

"""
	pancontig!(s::Sequence, ancestor::Dict{Int,Array{Interval}})

Return the ancestral tiling imprinted upon Sequence `s`. Modifies ancestor in place.
"""
function pancontig!(s::Sequence, ancestor::Dict{Int,Array{Interval}})
	isolate = NamedTuple{(:loci, :strand, :ancestor), Tuple{Interval,Bool,Int}}[]

	# helper functions
	function push_ancestor!(a, loci)
		if a ∈ keys(ancestor)
			atoms = ancestor[a]
			disjoint = IntervalSet(loci) \ IntervalSet(atoms[1].lo, atoms[end].hi, atoms)
			overlaps = [ x for atom in ancestor[a] for x in partition(atom,loci) ]
			if length(disjoint) == 0
				ancestor[a] = overlaps
			else
				ancestor[a] = sort(vcat(overlaps, disjoint.Is))
			end
        else
			ancestor[a] = [loci]
        end
	end

	function push_isolate!(a, loci, strand)
		push!(isolate, (
			loci=loci,
			strand=strand,
			ancestor=a
		))
	end

	function interval(left, right)
		left < right && return Interval(left, right+1)
		return Interval(right, left+1)
	end

	function increments(prev, curr, strand)
		strand && return prev+1 == curr # forward strand
		return prev-1 == curr # reverse strand
	end

	l  = 1
    id = unpack.ancestor(s[1])
    start = last = unpack.location(s[1])
	strand = unpack.strand(s[1]) == 0

	for (r,state) in enumerate(s[2:end])
        locus = unpack.location(state)
        ident = unpack.ancestor(state)
		polar = unpack.strand(state) == 0

        # still in tile
		if (polar == strand) && (ident == id) && increments(last, locus, strand)
            last = locus
            continue
        end

		push_ancestor!(id, interval(start, last))
		push_isolate!(id, interval(start, last), strand)

		l  = r+1
        id = ident
        start = last = locus
		strand = polar
    end

	# hanging block
	if start != last
		push_ancestor!(id, interval(start, last))
		push_isolate!(id, interval(start, last), strand)
	end

	return isolate
end

"""
	pancontigs(s::Sequence)

Return the ancestral tiling imprinted upon a set of Sequences `isolates`.
"""
function pancontigs(isolates::Array{Sequence})
	ancestors = Dict{Int, Array{Interval}}()
	isolates  = [ pancontig!(isolate, ancestors) for isolate in isolates ]

	return isolates, ancestors
end

function nodes(isolate, ancestors, blocks)
	ancestor = ancestors[isolate.ancestor]
	subpath  = Graphs.Node[]
	for interval in ancestor
		if isolate.loci ⊇ interval
			block = blocks[(isolate.ancestor,interval)]
			node  = Graphs.Node(block, isolate.strand)

			block.mutate[node] = Graphs.SNPMap()
			block.insert[node] = Graphs.InsMap()
			block.delete[node] = Graphs.DelMap()

			push!(subpath, node)
		end
	end
	return subpath
end

const ntab = UInt8['A', 'C', 'G', 'T' ]
"""
	nucleotide(sequence::Array{Sequence}, ancestor::Array{Array{UInt8,1},1})

Generate the set of extant sequences from the ancestral mosiacs `sequence` and the original sequences `ancestor`.
"""
function nucleotide(sequence::Array{Sequence}, ancestor::Array{Array{UInt8,1},1})
	function transform(x::UInt64)
		if unpack.snp(x) == 1
			return ntab[unpack.mutation(x)+1]
		else
			ident  = unpack.ancestor(x)
			locus  = unpack.location(x)
			return ancestor[ident][locus]
		end
	end
	transform(s::Sequence) = transform.(s)

	return transform.(sequence)
end

"""
	run(evolve!::Function, time::Int, initial::Array{Array{UInt8,1},1}; graph=false)

The high level API of the module.
Evolves a set of initial sequences `initial` for `time` generations using the one-step evolution function `evolve!`
If graph is true, the function will return the pangraph associated to the ancestral tiling.
"""
function run(evolve!::Function, time::Int, initial::Array{Array{UInt8,1},1}; graph=false)
	sequence  = population(length.(initial))
	offspring = [ Array{UInt64,1}(undef,maximum(length.(initial))) for _ in 1:length(initial) ]

    for _ in 1:time
        evolve!(sequence, offspring)
        sequence, offspring = offspring, sequence
    end

	if graph
		isolates, ancestors = pancontigs(sequence)

		# build the evolved pangraph from ancestral intervals
		blocks = Dict(
			(ancestor,interval) => Graphs.Block(initial[ancestor][interval])
			for (ancestor,intervals) in ancestors for interval in intervals
		)
		paths = [
			let
				node = [ n for interval in isolate for n in nodes(interval, ancestors, blocks) ]
				Graphs.Path("isolate_$(i)", node, 0, true, [])
			end for (i,isolate) in enumerate(isolates)
		]

		G = Graphs.Graph(Dict(pair(b) for b in values(blocks)), Dict(pair(p) for p in paths))
		Graphs.detransitive!(G)
		Graphs.finalize!(G)

		return nucleotide(sequence,initial), G
	end

	return nucleotide(sequence,initial), nothing
end

randseq(len::Int) = rand(UInt8.(['A','C','G','T']), len)

function test()
    N = 10
    L = Int(1e5)
	evolve! = model(Params(;N=N,L=L,snp=1e-3,hgt=0,inv=0))
	ancestors = [randseq(L) for _ in 1:N]

	sequences, G = run(evolve!, 20, ancestors; graph=true)
    open("data/synthetic/test.fa", "w") do io
        for (i,sequence) in enumerate(sequences)
            write_fasta(io, "isolate_$(i)", sequence)
        end
    end

    open("data/synthetic/test.json", "w") do io
        marshal(io, G; fmt=:json)
    end
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
module apo

include("apoCode.jl")

greet() = println("Hi there.")

export po, greet

end


text:
mutable struct AllMask
    featureMask
    valueMask
end

function getMask(mask::Dict)
    #получение маски из прочитанного словаря (из аттрибутов) в новый словарь с featureMask и valueMask
    newMask = Dict{String, Any}()
    for (key, val) in mask
        if !isa(val,Dict)
            newMask[key] = AllMask(val,val)
            newMask["Not"*key] = AllMask(val,0)
        else
            for (key2, val2) in val
                if val2 != 0 && !isequal(key2,"Mask")
                    newMask[key*"."*key2] = AllMask(val["Mask"],val2)
                end
            end
        end
    end
    return(newMask)
end

function parseType(atr)
    #преобразование аттрибутов в Dict маски типов
    mask = Dict{String, Any}()
    for str in atr
        if match(r"Type",str[1]) !== nothing
            sbstr = split(str[1],".")
            if length(sbstr[2:end]) > 1
                if !haskey(mask, sbstr[2])
                    mask[sbstr[2]] = Dict{String, Any}()
                end
                mask[sbstr[2]][sbstr[3]] = str[2]
            else
                mask[sbstr[2]] = str[2]
            end
        end
    end
    return(mask)
end

#получение типа int по набору признаков
function detIntType(mask::Dict,types::String)
    type = mask[types].featureMask & mask[types].valueMask
    typeint = Int(type)
end

#возвращает по типу имена входящих в него признаков
function detStringType(mask::Dict,type::Int)
    fields = mask.keys[mask.slots .!= 0]
    r = [checkType(type,mask[f]) for f in fields]
    fields[r]
end

#устанавливает в types биты, соответствующие bitName
function setbit(mask::Dict,types::Vector{T},bitName::Vector{String}) where T
    m = [mask[btN] for btN in bitName if haskey(mask, btN)]
    L = length(types)
    newTypes = Array{Int,1}()
    for i = 1:L
        type = types[i]
        for m1 in m
            type = type-type&m1.featureMask + m1.featureMask & m1.valueMask
        end
        push!(newTypes,type)
    end
    return newTypes
end

#построение дерева классов на основе класса и подкласса, с сохранением данных индексации попавших в категорию элементов
function buildQRStree(QPoint,Class, SubClass)

    unCl = unique(Class)
    QPoint = QPoint[1:length(Class)] #обрезаем по кол-ву классифицированных комплексов

    unSubCl = unique(SubClass)

    treeClasses = Dict{String,Any}()
    treeClasses["nodes"] = []
    treeClasses["name"] = "allClasses"
    treeClasses["size"] = length(Class)

    for i = 1:length(unCl)
        isCl = Class.==unCl[i]
        node =  Dict{String,Any}()
        node["nodes"] = []
        node["name"] = unCl[i]
        node["size"] = sum(isCl)
        node["QPoint"] = QPoint[isCl]
        for j = 1:length(unSubCl) #теперь добавляем подклассы
            isSubCl= (SubClass.== unSubCl[j]) .& isCl

            subnode =  Dict{String,Any}()
            subnode["nodes"] = []
            subnode["name"] = unSubCl[j]
            subnode["size"] = sum(isSubCl)
            subnode["QPoint"] = QPoint[isSubCl]
            push!(node["nodes"],subnode)
        end
        push!(treeClasses["nodes"],node)
    end

    result = JSON.json(treeClasses)
    return result
end


text:
function lx_reversed(com, _)
  brace_content = Franklin.content(com.braces[1])
  html = Franklin.fd2html(brace_content, internal=true)
  return "~~~" * replace(html, "<ol>" => "<ol reversed>") * "~~~"
end

max_stars_repo_name: jbmuir/JackBMuir

max_issues_repo_name: jbmuir/JackBMuir

max_forks_repo_name: jbmuir/JackBMuir


text:
# Author: Yuhang(Steven) Wang
# Date: 02/28/2017
# License: MIT/X11
include("main.jl")
import JSON

main(JSON.parse(ARGS[1]))

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
module FractionalCalculus

using QuadGK
using SpecialFunctions: gamma
using LinearAlgebra, InvertedIndices, SpecialMatrices
using SymbolicUtils, SymbolicUtils.Rewriters
using Symbolics

include("./Derivative/Derivative.jl")
include("./Integral/Integral.jl")

export FracDiffAlg, FracIntAlg
export fracdiff, fracint

# Export symbolics computing APIs
export semidiff, semiint
export SEMIDIFFRULES, SEMIINTRULES
export Incomplete_beta, AuxiliaryFresnelSin, AuxiliaryFresnelCos, Struve, MStruve, Legendre, SinIntegral, HyperSinIntegral, Hypergeometric1F1

# Export fractional derivative releating API
export Caputo, GL, RLDiff, Hadamard

# Caputo sense fractional derivative
export CaputoDirect, CaputoPiecewise, CaputoDiethelm, CaputoHighPrecision, CaputoHighOrder

# Grünwald Letnikov sense fractional derivative
export GLDirect, GLMultiplicativeAdditive, GLLagrangeThreePointInterp, GLFiniteDifference, GLHighPrecision

# Riemann Liouville sense fractional derivative
export RLDiffApprox, RLDiffMatrix, RLLinearSplineInterp, RLG1, RLD

# Hadamard sense fractional derivative
export HadamardLRect, HadamardRRect, HadamardTrap

# Riesz sense fractional derivative
export RieszSymmetric, RieszOrtigueira

# Macros for convenient computing
export @fracdiff, @semifracdiff

#Export fractional integral relating API
export RLInt, HadamardInt

export RLDirect, RLPiecewise
export RLIntApprox, RLLinearInterp, RLIntMatrix
export RLIntSimpson, RLIntTrapezoidal, RLIntRectangular, RLIntCubicSplineInterp

export HadamardMat

# Macros for convenient computing
export @fracint, @semifracint

# Auxiliary functions
export RieszMatrix, omega, B, genfun, first_order

end


text:
push!(LOAD_PATH, "/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/main")
using DHC_2DUtils
using DHC_tests
using Profile
using BenchmarkTools
using FFTW
using HDF5
using Test
using SparseArrays
using Plots
theme(:dark)

function fink_filter_bank_csafe(c, L; nx=256, wd=1, pc=1, shift=false, Omega=false, safety_on=true)
    #c     - sets the scale sampling rate (1 is dyadic, 2 is half dyadic)
    #L     - number of angular bins (usually 8*pc or 16*pc)
    #wd    - width of the wavelets (default 1, wd=2 for a double covering)
    #pc    - plane coverage (default 1, full 2pi 2)
    #shift - shift in θ by 1/2 of the θ spacing
    #Omega - true= append Omega filter (all power beyond Nyquist) so the sum of filters is 1.0

    # -------- assertion errors to make sure arguments are reasonable
    @test wd <= L/2

    # -------- set parameters
    dθ   = pc*π/L
    θ_sh = shift ? dθ/2 : 0.0
    dx   = nx/2-1

    im_scale = convert(Int8,log2(nx))
    # -------- number of bins in radial direction (size scales)
    J = (im_scale-3)*c + 1

    # -------- allocate output array of zeros
    filt      = zeros(nx, nx, J*L+(Omega ? 2 : 1))
    psi_index = zeros(Int32, J, L)
    psi_ind_in= zeros(Int32, J*L+(Omega ? 2 : 1), 2)
    theta     = zeros(Float64, L)
    j_value   = zeros(Float64, J)

    # -------- compute the required wd
    j_rad_exp = zeros(J)
    for j_ind = 1:J
        j = (j_ind-1)/c
        jrad  = im_scale-j-2
        j_rad_exp[j_ind] = 2^(jrad)
    end

    wd_j = max.(ceil.(L./(pc.*π.*j_rad_exp)),wd)

    if !safety_on
        wd_j.=wd
    end

    # loop over wd from small to large
    ## there is some uneeded redundancy added by doing this esp in l loop
    for wd in sort(unique(wd_j))
        # -------- allocate theta and logr arrays
        θ    = zeros(nx, nx)
        logr = zeros(nx, nx)

        wdθ  = wd*dθ
        norm = 1.0/sqrt(wd)
        # -------- loop over l
        for l = 0:L-1
            θ_l        = dθ*l+θ_sh
            theta[l+1] = θ_l

        # -------- allocate anggood BitArray
            anggood = falses(nx, nx)

        # -------- loop over pixels
            for x = 1:nx
                sx = mod(x+dx,nx)-dx -1    # define sx,sy so that no fftshift() needed
                for y = 1:nx
                    sy = mod(y+dx,nx)-dx -1
                    θ_pix  = mod(atan(sy, sx)+π -θ_l, 2*π)
                    θ_good = abs(θ_pix-π) <= wdθ

                    # If this is a pixel we might use, calculate log2(r)
                    if θ_good
                        anggood[y, x] = θ_good
                        θ[y, x]       = θ_pix
                        r2            = sx^2 + sy^2
                        logr[y, x]    = 0.5*log2(max(1,r2))
                    end
                end
            end
            angmask = findall(anggood)
        # -------- compute the wavelet in the Fourier domain
        # -------- the angular factor is the same for all j
            F_angular = norm .* cos.((θ[angmask].-π).*(L/(2*wd*pc)))

        # -------- loop over j for the radial part
        #    for (j_ind, j) in enumerate(1/c:1/c:im_scale-2)
            j_ind_w_wd = findall(wd_j.==wd)
            for j_ind in j_ind_w_wd
                j = (j_ind-1)/c
                j_value[j_ind] = 1+j  # store for later
                jrad  = im_scale-j-2
                Δj    = abs.(logr[angmask].-jrad)
                rmask = (Δj .<= 1) #deprecating the 1/c to 1, constant width

        # -------- radial part
                F_radial = cos.(Δj[rmask] .* (π/2)) #deprecating c*π/2 to π/2
                ind      = angmask[rmask]
        #      Let's have these be (J,L) if you reshape...
        #        f_ind    = (j_ind-1)*L+l+1
                f_ind    = j_ind + l*J
                filt[ind, f_ind] = F_radial .* F_angular[rmask]
                psi_index[j_ind,l+1] = f_ind
                psi_ind_in[f_ind,:] = [j_ind-1,l]
            end
        end
    end

    # -------- phi contains power near k=0 not yet accounted for
    filter_power = (sum(filt.*filt, dims=3))[:,:,1]

    # -------- for plane half-covered (pc=1), add other half-plane
    if pc == 1
        filter_power .+= circshift(filter_power[end:-1:1,end:-1:1],(1,1))
    end

    # -------- compute power required to sum to 1.0
    i0 = round(Int16,nx/2-2)
    i1 = round(Int16,nx/2+4)
    center_power = 1.0 .- fftshift(filter_power)[i0:i1,i0:i1]
    zind = findall(center_power .< 1E-15)
    center_power[zind] .= 0.0  # set small numbers to zero
    phi_cen = zeros(nx, nx)
    phi_cen[i0:i1,i0:i1] = sqrt.(center_power)

    # -------- before adding ϕ to filter bank, renormalize ψ if pc=1
    if pc==1 filt .*= sqrt(2.0) end  # double power for half coverage

    # -------- add result to filter array
    phi_index  = J*L+1
    filt[:,:,phi_index] .= fftshift(phi_cen)
    psi_ind_in[phi_index,:] = [J,0]

    # -------- metadata dictionary
    info=Dict{String,Any}()
    info["npix"]         = nx
    info["j_value"]      = j_value
    info["theta_value"]  = theta
    info["psi_index"]    = psi_index
    info["phi_index"]    = phi_index
    info["J_L"]          = psi_ind_in
    info["pc"]           = pc
    info["wd"]           = wd_j
    info["fs_center_r"]  = j_rad_exp

    if Omega     # append a filter containing the rest (outside Nyquist)
        filter_power += filt[:,:,phi_index].^2
        edge_power    = 1.0 .- filter_power
        zind          = findall(edge_power .< 1E-15)
        edge_power[zind]     .= 0.0  # set small numbers to zero
        Omega_index           = J*L+2
        info["Omega_index"]   = Omega_index
        filt[:,:,Omega_index] = sqrt.(edge_power)
    end

    return filt, info
end

function fink_filter_hash_csafe(c, L; nx=256, wd=1, pc=1, shift=false, Omega=false, safety_on=true)
    # -------- compute the filter bank
    filt, hash = fink_filter_bank_csafe(c, L; nx=nx, wd=wd, pc=pc, shift=shift, Omega=Omega, safety_on=safety_on)

    # -------- list of non-zero pixels
    flist = fink_filter_list(filt)

    # -------- pack everything you need into the info structure
    hash["filt_index"] = flist[1]
    hash["filt_value"] = flist[2]

    # -------- compute matrix that projects iso coeffs, add to hash
    S1_iso_mat = S1_iso_matrix(hash)
    hash["S1_iso_mat"] = S1_iso_mat
    S2_iso_mat = S2_iso_matrix(hash)
    hash["S2_iso_mat"] = S2_iso_mat
    hash["num_iso_coeff"] = size(S1_iso_mat)[1] + size(S2_iso_mat)[1] + 2
    hash["num_coeff"] = size(S1_iso_mat)[2] + size(S2_iso_mat)[2] + 2

    return hash
end

fhash = fink_filter_hash_csafe(1,8,wd=2,pc=1)
filt_test = fink_filter_bank_csafe(1,8,wd=2,pc=1)[1]

plot_filter_bank_QA(filt_test, fhash; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_csafe.png")

fhash = fink_filter_hash_csafe(2,8,wd=2,pc=1)
filt_test = fink_filter_bank_csafe(2,8,wd=2,pc=1)[1]

plot_filter_bank_QA(filt_test, fhash; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_csafe_C2.png")

fhash = fink_filter_hash_csafe(2,8,wd=1,pc=1)

fhash = fink_filter_hash(2,8,wd=1,pc=1)

#In both formulations, only the smallest FS wavelet gets widened by our original criteria, even with c=2

fhash = fink_filter_hash_csafe(2,8,wd=1,pc=1)

function fink_filter_bank_csafe(c, L; nx=256, wd=1, pc=1, shift=false, Omega=false, safety_on=true, wd_cutoff=2)
    #c     - sets the scale sampling rate (1 is dyadic, 2 is half dyadic)
    #L     - number of angular bins (usually 8*pc or 16*pc)
    #wd    - width of the wavelets (default 1, wd=2 for a double covering)
    #pc    - plane coverage (default 1, full 2pi 2)
    #shift - shift in θ by 1/2 of the θ spacing
    #Omega - true= append Omega filter (all power beyond Nyquist) so the sum of filters is 1.0

    # -------- assertion errors to make sure arguments are reasonable
    @test wd <= L/2

    # -------- set parameters
    dθ   = pc*π/L
    θ_sh = shift ? dθ/2 : 0.0
    dx   = nx/2-1

    im_scale = convert(Int8,log2(nx))
    # -------- number of bins in radial direction (size scales)
    J = (im_scale-3)*c + 1
    normj = 1/sqrt(c)

    # -------- allocate output array of zeros
    filt      = zeros(nx, nx, J*L+(Omega ? 2 : 1))
    psi_index = zeros(Int32, J, L)
    psi_ind_in= zeros(Int32, J*L+(Omega ? 2 : 1), 2)
    theta     = zeros(Float64, L)
    j_value   = zeros(Float64, J)

    # -------- compute the required wd
    j_rad_exp = zeros(J)
    for j_ind = 1:J
        j = (j_ind-1)/c
        jrad  = im_scale-j-2
        j_rad_exp[j_ind] = 2^(jrad)
    end

    wd_j = max.(ceil.(wd_cutoff.*L./(pc.*π.*j_rad_exp)),wd)

    if !safety_on
        wd_j.=wd
    end

    # loop over wd from small to large
    ## there is some uneeded redundancy added by doing this esp in l loop
    for wd in sort(unique(wd_j))
        # -------- allocate theta and logr arrays
        θ    = zeros(nx, nx)
        logr = zeros(nx, nx)

        wdθ  = wd*dθ
        norm = 1.0/(sqrt(wd))
        # -------- loop over l
        for l = 0:L-1
            θ_l        = dθ*l+θ_sh
            theta[l+1] = θ_l

        # -------- allocate anggood BitArray
            anggood = falses(nx, nx)

        # -------- loop over pixels
            for x = 1:nx
                sx = mod(x+dx,nx)-dx -1    # define sx,sy so that no fftshift() needed
                for y = 1:nx
                    sy = mod(y+dx,nx)-dx -1
                    θ_pix  = mod(atan(sy, sx)+π -θ_l, 2*π)
                    θ_good = abs(θ_pix-π) <= wdθ

                    # If this is a pixel we might use, calculate log2(r)
                    if θ_good
                        anggood[y, x] = θ_good
                        θ[y, x]       = θ_pix
                        r2            = sx^2 + sy^2
                        logr[y, x]    = 0.5*log2(max(1,r2))
                    end
                end
            end
            angmask = findall(anggood)
        # -------- compute the wavelet in the Fourier domain
        # -------- the angular factor is the same for all j
            F_angular = norm .* cos.((θ[angmask].-π).*(L/(2*wd*pc)))

        # -------- loop over j for the radial part
        #    for (j_ind, j) in enumerate(1/c:1/c:im_scale-2)
            j_ind_w_wd = findall(wd_j.==wd)
            for j_ind in j_ind_w_wd
                j = (j_ind-1)/c
                j_value[j_ind] = 1+j  # store for later
                jrad  = im_scale-j-2
                Δj    = abs.(logr[angmask].-jrad)
                rmask = (Δj .<= 1) #deprecating the 1/c to 1, constant width

        # -------- radial part
                F_radial = normj .* cos.(Δj[rmask] .* (π/2)) #deprecating c*π/2 to π/2
                ind      = angmask[rmask]
        #      Let's have these be (J,L) if you reshape...
        #        f_ind    = (j_ind-1)*L+l+1
                f_ind    = j_ind + l*J
                filt[ind, f_ind] = F_radial .* F_angular[rmask]
                psi_index[j_ind,l+1] = f_ind
                psi_ind_in[f_ind,:] = [j_ind-1,l]
            end
        end
    end

    # -------- phi contains power near k=0 not yet accounted for
    filter_power = (sum(filt.*filt, dims=3))[:,:,1]

    # -------- for plane half-covered (pc=1), add other half-plane
    if pc == 1
        filter_power .+= circshift(filter_power[end:-1:1,end:-1:1],(1,1))
    end

    # -------- compute power required to sum to 1.0
    i0 = round(Int16,nx/2-2)
    i1 = round(Int16,nx/2+4)
    center_power = 1.0 .- fftshift(filter_power)[i0:i1,i0:i1]
    zind = findall(center_power .< 1E-15)
    center_power[zind] .= 0.0  # set small numbers to zero
    phi_cen = zeros(nx, nx)
    phi_cen[i0:i1,i0:i1] = sqrt.(center_power)

    # -------- before adding ϕ to filter bank, renormalize ψ if pc=1
    if pc==1 filt .*= sqrt(2.0) end  # double power for half coverage

    # -------- add result to filter array
    phi_index  = J*L+1
    filt[:,:,phi_index] .= fftshift(phi_cen)
    psi_ind_in[phi_index,:] = [J,0]

    # -------- metadata dictionary
    info=Dict{String,Any}()
    info["npix"]         = nx
    info["j_value"]      = j_value
    info["theta_value"]  = theta
    info["psi_index"]    = psi_index
    info["phi_index"]    = phi_index
    info["J_L"]          = psi_ind_in
    info["pc"]           = pc
    info["wd"]           = wd_j
    info["wd_cutoff"]    = wd_cutoff
    info["fs_center_r"]  = j_rad_exp

    if Omega     # append a filter containing the rest (outside Nyquist)
        filter_power += filt[:,:,phi_index].^2
        edge_power    = 1.0 .- filter_power
        zind          = findall(edge_power .< 1E-15)
        edge_power[zind]     .= 0.0  # set small numbers to zero
        Omega_index           = J*L+2
        info["Omega_index"]   = Omega_index
        filt[:,:,Omega_index] = sqrt.(edge_power)
    end

    return filt, info
end

fhash = fink_filter_hash_csafe(2,8,wd=2,pc=1)
filt_test = fink_filter_bank_csafe(2,8,wd=2,pc=1)[1]

plot_filter_bank_QA(filt_test, fhash; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_csafe_C2_wide.png")

#Now, only the smallest FS wavelet with wd=2 is bumped up to 3... visually seems a good choice.

fhash = fink_filter_hash_csafe(1,8,wd=2,pc=1)
filt_test = fink_filter_bank_csafe(1,8,wd=2,pc=1)[1]

plot_filter_bank_QA(filt_test, fhash; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_csafe_C1_wide.png")

fhash = fink_filter_hash_csafe(1,8,wd=1,pc=1)
filt_test = fink_filter_bank_csafe(1,8,wd=1,pc=1)[1]

plot_filter_bank_QA(filt_test, fhash; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_csafe_wd1_C1_wide.png")

# with wd=1, smallest is widened to 3, next to 2, rest to 1

heatmap(fftshift(dropdims(sum(filt_test[:,:,:].^2,dims=3),dims=3)))

fhash = fink_filter_bank(3,16,wd=2,pc=1)

heatmap(fftshift(dropdims(sum(fhash[1].^2,dims=3),dims=3)))

fhash = fink_filter_bank_csafe(4,16,wd=2,pc=1)

heatmap(fftshift(dropdims(sum(fhash[1].^2,dims=3),dims=3)))

fhash = fink_filter_hash_csafe(1,8,wd=1,pc=1,safety_on=false)
filt_test = fink_filter_bank_csafe(1,8,wd=1,pc=1,safety_on=false)[1]

plot_filter_bank_QA(filt_test, fhash; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_csafe_wd1_C1_wide_unsafe.png")

filt_test = fink_filter_bank_csafe(1,8,wd=1,pc=2,safety_on=false)[1]
heatmap(fftshift(dropdims(sum(filt_test[:,:,:].^2,dims=3),dims=3)))

filt_test = fink_filter_bank_csafe(2,8,wd=4,pc=2)[1]
heatmap(fftshift(dropdims(sum(filt_test[:,:,:].^2,dims=3),dims=3)))

filt_test = fink_filter_bank_csafe(3,8,wd=4,pc=2)[1]
heatmap(fftshift(dropdims(sum(filt_test[:,:,:].^2,dims=3),dims=3)))

function fink_filter_bank_csafe(c, L; nx=256, wd=1, pc=1, shift=false, Omega=false, safety_on=true, wd_cutoff=2)
    #c     - sets the scale sampling rate (1 is dyadic, 2 is half dyadic)
    #L     - number of angular bins (usually 8*pc or 16*pc)
    #wd    - width of the wavelets (default 1, wd=2 for a double covering)
    #pc    - plane coverage (default 1, full 2pi 2)
    #shift - shift in θ by 1/2 of the θ spacing
    #Omega - true= append Omega filter (all power beyond Nyquist) so the sum of filters is 1.0

    # -------- assertion errors to make sure arguments are reasonable
    @test wd <= L/2

    # -------- set parameters
    dθ   = pc*π/L
    θ_sh = shift ? dθ/2 : 0.0
    dx   = nx/2-1

    im_scale = convert(Int8,log2(nx))
    # -------- number of bins in radial direction (size scales)
    J = (im_scale-3)*c + 1
    normj = 1/sqrt(c)

    # -------- allocate output array of zeros
    filt      = zeros(nx, nx, J*L+(Omega ? 2 : 1))
    psi_index = zeros(Int32, J, L)
    psi_ind_in= zeros(Int32, J*L+(Omega ? 2 : 1), 2)
    theta     = zeros(Float64, L)
    j_value   = zeros(Float64, J)

    # -------- compute the required wd
    j_rad_exp = zeros(J)
    for j_ind = 1:J
        j = (j_ind-1)/c
        jrad  = im_scale-j-2
        j_rad_exp[j_ind] = 2^(jrad)
    end

    wd_j = max.(ceil.(wd_cutoff.*L./(pc.*π.*j_rad_exp)),wd)

    if !safety_on
        wd_j.=wd
    end

    # loop over wd from small to large
    ## there is some uneeded redundancy added by doing this esp in l loop
    for wd in sort(unique(wd_j))
        # -------- allocate theta and logr arrays
        θ    = zeros(nx, nx)
        logr = zeros(nx, nx)

        wdθ  = wd*dθ
        norm = 1.0/(sqrt(wd))
        # -------- loop over l
        for l = 0:L-1
            θ_l        = dθ*l+θ_sh
            theta[l+1] = θ_l

        # -------- allocate anggood BitArray
            anggood = falses(nx, nx)

        # -------- loop over pixels
            for x = 1:nx
                sx = mod(x+dx,nx)-dx -1    # define sx,sy so that no fftshift() needed
                for y = 1:nx
                    sy = mod(y+dx,nx)-dx -1
                    θ_pix  = mod(atan(sy, sx)+π -θ_l, 2*π)
                    θ_good = abs(θ_pix-π) <= wdθ

                    # If this is a pixel we might use, calculate log2(r)
                    if θ_good
                        anggood[y, x] = θ_good
                        θ[y, x]       = θ_pix
                        r2            = sx^2 + sy^2
                        logr[y, x]    = 0.5*log2(max(1,r2))
                    end
                end
            end
            angmask = findall(anggood)
        # -------- compute the wavelet in the Fourier domain
        # -------- the angular factor is the same for all j
            F_angular = norm .* cos.((θ[angmask].-π).*(L/(2*wd*pc)))

        # -------- loop over j for the radial part
        #    for (j_ind, j) in enumerate(1/c:1/c:im_scale-2)
            j_ind_w_wd = findall(wd_j.==wd)
            for j_ind in j_ind_w_wd
                j = (j_ind-1)/c
                j_value[j_ind] = 1+j  # store for later
                jrad  = im_scale-j-2
                Δj    = abs.(logr[angmask].-jrad)
                rmask = (Δj .<= 1) #deprecating the 1/c to 1, constant width

        # -------- radial part
                F_radial = normj .* cos.(Δj[rmask] .* (π/2)) #deprecating c*π/2 to π/2
                ind      = angmask[rmask]
        #      Let's have these be (J,L) if you reshape...
        #        f_ind    = (j_ind-1)*L+l+1
                f_ind    = j_ind + l*J
                filt[ind, f_ind] = F_radial .* F_angular[rmask]
                psi_index[j_ind,l+1] = f_ind
                psi_ind_in[f_ind,:] = [j_ind-1,l]
            end
        end
    end

    # -------- phi contains power near k=0 not yet accounted for
    filter_power = (sum(filt.*filt, dims=3))[:,:,1]

    # -------- for plane half-covered (pc=1), add other half-plane
    if pc == 1
        filter_power .+= circshift(filter_power[end:-1:1,end:-1:1],(1,1))
    end

    # -------- compute power required to sum to 1.0
    i0 = round(Int16,nx/2-2)
    i1 = round(Int16,nx/2+4)
    center_power = 1.0 .- fftshift(filter_power)[i0:i1,i0:i1]
    zind = findall(center_power .< 1E-15)
    center_power[zind] .= 0.0  # set small numbers to zero
    phi_cen = zeros(nx, nx)
    phi_cen[i0:i1,i0:i1] = sqrt.(center_power)

    # -------- before adding ϕ to filter bank, renormalize ψ if pc=1
    if pc==1 filt .*= sqrt(2.0) end  # double power for half coverage

    # -------- add result to filter array
    phi_index  = J*L+1
    filt[:,:,phi_index] .= fftshift(phi_cen)
    psi_ind_in[phi_index,:] = [J,0]

    # -------- metadata dictionary
    info=Dict{String,Any}()
    info["npix"]         = nx
    info["j_value"]      = j_value
    info["theta_value"]  = theta
    info["psi_index"]    = psi_index
    info["phi_index"]    = phi_index
    info["J_L"]          = psi_ind_in
    info["pc"]           = pc
    info["wd"]           = wd_j
    info["wd_cutoff"]    = wd_cutoff
    info["fs_center_r"]  = j_rad_exp

    if Omega     # append a filter containing the rest (outside Nyquist)
        filter_power += filt[:,:,phi_index].^2
        edge_power    = 1.0 .- filter_power
        zind          = findall(edge_power .< 1E-15)
        edge_power[zind]     .= 0.0  # set small numbers to zero
        Omega_index           = J*L+2
        info["Omega_index"]   = Omega_index
        filt[:,:,Omega_index] = sqrt.(edge_power)
    end

    return filt, info
end

filt_test = fink_filter_bank_csafe(1,8,wd=1,pc=1)[1]
heatmap(fftshift(dropdims(sum(filt_test[:,:,:].^2,dims=3),dims=3)))

filter_power = (sum(filt.*filt, dims=3))[:,:,1]

# -------- for plane half-covered (pc=1), add other half-plane
if pc == 1
    filter_power .+= circshift(filter_power[end:-1:1,end:-1:1],(1,1))
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
# for models that are "exported" learning networks (return a Node as
# their fit-result; see MLJ docs:
abstract type ProbabilisticNetwork <: Probabilistic end
abstract type DeterministicNetwork <: Deterministic end
abstract type  UnsupervisedNetwork <: Unsupervised end

## ABSTRACT NODES AND SOURCE NODES

abstract type AbstractNode <: MLJType end

# K is :target, :input, :weight or :unknown
mutable struct Source{K} <: AbstractNode
    data  # training data
end

"""
    Xs = source(X)
    ys = source(y, kind=:target)
    ws = source(w, kind=:weight)

Defines, respectively, learning network `Source` objects for wrapping
some input data `X` (`kind=:input`), some target data `y`, or some
sample weights `w`.  The values of each variable `X, y, w` can be
anything, even `nothing`, if the network is for exporting as a
stand-alone model only. For training and testing the unexported network,
appropriate vectors, tables, or other data containers are expected.

    Xs = source()
    ys = source(kind=:target)
    ws = source(kind=:weight)

Define source nodes wrapping `nothing` instead of concrete data. Such
definitions suffice if a learning network is to be exported without
testing.

The calling behaviour of a `Source` object is this:

    Xs() = X
    Xs(rows=r) = selectrows(X, r)  # eg, X[r,:] for a DataFrame
    Xs(Xnew) = Xnew

See also: [`@from_network`](@ref], [`sources`](@ref),
[`origins`](@ref), [`node`](@ref).

"""
function source(X; kind=:input)
    kind in [:input, :target, :weights] ||
        @warn "`Source` kind is neither :input, :target, or :weights. "
    return Source{kind}(X)
end

source(X::Source; args...) = X
source(; args...) = source(nothing; args...)

kind(::Source{k}) where k = k
Base.isempty(X::Source) = X.data == nothing

is_stale(s::Source) = false

# make source nodes callable:
function (s::Source)(; rows=:)
    rows == (:) && return s.data
    return selectrows(s.data, rows)
end
(s::Source)(Xnew) = Xnew

"""
    rebind!(s)

Attach new data `X` to an existing source node `s`.
"""
function rebind!(s::Source, X)
    s.data = X
    return s
end

"""
    origins(N)

Return a list of all origins of a node `N` accessed by a call `N()`.
These are the source nodes of the acyclic directed graph (DAG)
associated with the learning network terminating at `N`, if edges
corresponding to training arguments are excluded. A `Node` object
cannot be called on new data unless it has a unique origin.

Not to be confused with `sources(N)` which refers to the same graph
but without the training edge deletions.

See also: [`node`](@ref), [`source`](@ref).
"""
origins(s::Source) = [s,]

# Function to incrementally append to `nodes1` all elements in
# `nodes2`, excluding any elements previously added (or any element of
# `nodes1` in its initial state):
function _merge!(nodes1, nodes2)
    for x in nodes2
        if !(x in nodes1)
            push!(nodes1, x)
        end
    end
    return nodes1
end

# Note that `fit!` has already been defined for any AbstractMachine in
# machines.jl

mutable struct NodalMachine{M<:Model} <: AbstractMachine{M}
    model::M
    previous_model::M # for remembering the model used in last call to `fit!`
    fitresult
    cache
    args::Tuple{Vararg{AbstractNode}}
    report
    frozen::Bool
    previous_rows   # for remembering the rows used in last call to `fit!`
    state::Int      # number of times fit! has been called on machine
    upstream_state  # for remembering the upstream state in last call to `fit!`

    function NodalMachine{M}(model::M, args::AbstractNode...) where M<:Model

        # check number of arguments for model subtypes:
        !(M <: Supervised) || length(args) > 1 ||
            throw(ArgumentError("Wrong number of arguments. " *
                        "Use `machine(model, X, y)` or "*
                        "`machine(model, X, y, w)` for supervised models."))

        if M <: Unsupervised && !(M <: Static)
            length(args) == 1 ||
                throw(ArgumentError("Wrong number of arguments. Use " *
                        "`machine(model, X)` for an unsupervised model "*
                        "(or `machine(model)` if there are no training "*
                        "arguments (\"static\" tranformers"))
        end

        machine = new{M}(model)
        machine.frozen = false
        machine.state = 0
        machine.args = args

        machine.upstream_state = Tuple([state(arg) for arg in args])

        return machine
    end
end

NodalMachine(model::M, args...) where M<:Model = NodalMachine{M}(model, args...)

# Note: freeze! and thaw! are possibly not used within MLJ itself.

"""
    freeze!(mach)

Freeze the machine `mach` so that it will never be retrained (unless
thawed).

See also [`thaw!`](@ref).
"""
function freeze!(machine::NodalMachine)
    machine.frozen = true
end

"""
    thaw!(mach)

Unfreeze the machine `mach` so that it can be retrained.

See also [`freeze!`](@ref).
"""
function thaw!(machine::NodalMachine)
    machine.frozen = false
end

"""
    is_stale(mach)

Check if a machine `mach` is stale.

See also [`fit!`](@ref)
"""
function is_stale(machine::NodalMachine)
    !isdefined(machine, :fitresult) ||
        machine.model != machine.previous_model ||
        reduce(|,[is_stale(arg) for arg in machine.args])
end

"""
    state(mach)

Return the state of a machine, `mach`.
"""
state(machine::NodalMachine) = machine.state

## NODES

struct Node{T<:Union{NodalMachine, Nothing}} <: AbstractNode
    operation   # that can be dispatched on a fit-result (eg, `predict`) or a static operation
    machine::T  # is `nothing` for static operations
    args::Tuple{Vararg{AbstractNode}}  # nodes where `operation` looks for its arguments
    origins::Vector{Source}
    nodes::Vector{AbstractNode}  # all upstream nodes, order consistent with DAG order (the node "tape")

    function Node{T}(operation,
                     machine::T,
                     args::AbstractNode...) where {M<:Model, T<:Union{NodalMachine{M},Nothing}}

        # check the number of arguments:
        if machine === nothing && isempty(args)
            throw(error("`args` in `Node(::Function, args...)` must be non-empty. "))
        end

        origins_ = unique(vcat([origins(arg) for arg in args]...))
        # length(origins_) == 1 ||
        #     @warn "A node referencing multiple origins when called " *
        #           "has been defined:\n$(origins_). "

        # initialize the list of upstream nodes:
        nodes_ = AbstractNode[]

        # merge the lists from arguments:
        for arg in args
            _merge!(nodes_, nodes(arg))
        end

        # merge the lists from training arguments:
        if machine !== nothing
            for arg in machine.args
                _merge!(nodes_, nodes(arg))
            end
        end

        return new{T}(operation, machine, args, origins_, nodes_)
    end
end

origins(X::Node) = X.origins

"""
    nodes(N)

Return all nodes upstream of a node `N`, including `N` itself, in an
order consistent with the extended directed acyclic graph of the
network. Here "extended" means edges corresponding to training
arguments are included.
"""
nodes(X::Node) = AbstractNode[X.nodes..., X]
nodes(S::Source) = AbstractNode[S, ]

"""
    is_stale(N)

Check if a node `N` is stale.
"""
function is_stale(X::Node)
    (X.machine !== nothing && is_stale(X.machine)) ||
        reduce(|, [is_stale(arg) for arg in X.args])
end

state(s::Source) = (state = 0, )

"""
    state(N)

Return the state of a node `N`
"""
function state(W::Node)
    mach = W.machine
    state_ = W.machine === nothing ? 0 : state(W.machine)
    if mach === nothing
        trainkeys   = []
        trainvalues = []
    else
        trainkeys   = (Symbol("train_arg", i) for i in eachindex(mach.args))
        trainvalues = (state(arg) for arg in mach.args)
    end
    keys = tuple(:state,
                 (Symbol("arg", i) for i in eachindex(W.args))...,
                 trainkeys...)
    values = tuple(state_,
                   (state(arg) for arg in W.args)...,
                   trainvalues...)
    return NamedTuple{keys}(values)
end

# autodetect type parameter:
Node(operation, machine::M, args...) where M <: Union{NodalMachine,Nothing} =
    Node{M}(operation, machine, args...)

# constructor for static operations:
Node(operation, args::AbstractNode...) = Node(operation, nothing, args...)

# make nodes callable:
(y::Node)(; rows=:) =
    (y.operation)(y.machine, [arg(rows=rows) for arg in y.args]...)
function (y::Node)(Xnew)
    length(y.origins) == 1 ||
        error("Node $y has multiple origins and cannot be called "*
              "on new data. ")
    return (y.operation)(y.machine, [arg(Xnew) for arg in y.args]...)
end

# and for the special case of static operations:
(y::Node{Nothing})(; rows=:) =
    (y.operation)([arg(rows=rows) for arg in y.args]...)
function (y::Node{Nothing})(Xnew)
    length(y.origins) == 1 ||
        error("Node $y has multiple origins and cannot be called "*
              "on new data. ")
    return (y.operation)([arg(Xnew) for arg in y.args]...)
end

"""
    fit!(N::Node; rows=nothing, verbosity::Int=1, force::Bool=false)

Train all machines in the learning network terminating at node `N`, in an
appropriate order. These machines are those returned by `machines(N)`.

"""
function fit!(y::Node; rows=nothing, verbosity::Int=1,
              force::Bool=false)

    if rows === nothing rows = (:) end

    # get non-source nodes:
    nodes_ = filter(nodes(y)) do n
        n isa Node
    end

    # get machines to fit:
    machines = map(n -> n.machine, nodes_)
    machines = filter(unique(machines)) do mach
        mach !== nothing
    end

    for mach in machines
        fit!(mach; rows=rows, verbosity=verbosity, force=force)
    end

    return y
end

# allow arguments of `Nodes` and `NodalMachine`s to appear
# at REPL:
istoobig(d::Tuple{AbstractNode}) = length(d) > 10

# overload show method

function _recursive_show(stream::IO, X::AbstractNode)
    if X isa Source
        printstyled(IOContext(stream, :color=>SHOW_COLOR), handle(X), color=:blue)
    else
        detail = (X.machine === nothing ? "(" : "($(handle(X.machine)), ")
        operation_name = typeof(X.operation).name.mt.name
        print(stream, operation_name, "(")
        if X.machine !== nothing
            color = (X.machine.frozen ? :red : :green)
            printstyled(IOContext(stream, :color=>SHOW_COLOR), handle(X.machine),
                        bold=SHOW_COLOR)
            print(stream, ", ")
        end
        n_args = length(X.args)
        counter = 1
        for arg in X.args
            _recursive_show(stream, arg)
            counter >= n_args || print(stream, ", ")
            counter += 1
        end
        print(stream, ")")
    end
end

function Base.show(stream::IO, ::MIME"text/plain", X::Node)
    id = objectid(X)
    description = string(typeof(X).name.name)
    str = "$description @ $(handle(X))"
    printstyled(IOContext(stream, :color=>SHOW_COLOR), str, color=:blue)
#    if !(X isa Source)
        print(stream, " = ")
        _recursive_show(stream, X)
#    end
end

function Base.show(stream::IO, ::MIME"text/plain", machine::NodalMachine)
    id = objectid(machine)
    description = string(typeof(machine).name.name)
    str = "$description @ $(handle(machine))"
    printstyled(IOContext(stream, :color=>SHOW_COLOR), str, bold=SHOW_COLOR)
    print(stream, " = ")
    print(stream, "machine($(machine.model), ")
    n_args = length(machine.args)
    counter = 1
    for arg in machine.args
        printstyled(IOContext(stream, :color=>SHOW_COLOR), handle(arg), bold=SHOW_COLOR)
        counter >= n_args || print(stream, ", ")
        counter += 1
    end
    print(stream, ")")
end

## SYNTACTIC SUGAR FOR LEARNING NETWORKS

"""
    N = node(f::Function, args...)

Defines a `Node` object `N` wrapping a static operation `f` and arguments
`args`. Each of the `n` elements of `args` must be a `Node` or `Source`
object. The node `N` has the following calling behaviour:

    N() = f(args[1](), args[2](), ..., args[n]())
    N(rows=r) = f(args[1](rows=r), args[2](rows=r), ..., args[n](rows=r))
    N(X) = f(args[1](X), args[2](X), ..., args[n](X))

    J = node(f, mach::NodalMachine, args...)

Defines a dynamic `Node` object `J` wrapping a dynamic operation `f`
(`predict`, `predict_mean`, `transform`, etc), a nodal machine `mach` and
arguments `args`. Its calling behaviour, which depends on the outcome of
training `mach` (and, implicitly, on training outcomes affecting its
arguments) is this:

    J() = f(mach, args[1](), args[2](), ..., args[n]())
    J(rows=r) = f(mach, args[1](rows=r), args[2](rows=r), ..., args[n](rows=r))
    J(X) = f(mach, args[1](X), args[2](X), ..., args[n](X))

Generally `n=1` or `n=2` in this latter case.

    predict(mach, X::AbsractNode, y::AbstractNode)
    predict_mean(mach, X::AbstractNode, y::AbstractNode)
    predict_median(mach, X::AbstractNode, y::AbstractNode)
    predict_mode(mach, X::AbstractNode, y::AbstractNode)
    transform(mach, X::AbstractNode)
    inverse_transform(mach, X::AbstractNode)

Shortcuts for `J = node(predict, mach, X, y)`, etc.

Calling a node is a recursive operation which terminates in the call
to a source node (or nodes). Calling nodes on *new* data `X` fails unless the
number of such nodes is one.

See also: [`source`](@ref), [`origins`](@ref).

"""
node = Node

# if `args` in `machine(model, args...)` is empty or includes an
# `AbstractNode`, a `NodalMachine` is to be created.
machine(m::Unsupervised, a::AbstractNode...) =
    NodalMachine(m, a...)
# machine(m::Static, a1::AbstractNode, a::AbstractNode...) =
#     NodalMachine(m, a...)
machine(m::Supervised, a1::AbstractNode, a::AbstractNode...) =
    NodalMachine(m, a1, a...)
machine(model::Model, X, y::AbstractNode) = NodalMachine(model, source(X), y)
machine(model::Model, X::AbstractNode, y) = NodalMachine(model, X, source(y))

MMI.matrix(X::AbstractNode)      = node(matrix, X)
MMI.table(X::AbstractNode)       = node(table, X)
Base.vcat(args::AbstractNode...) = node(vcat, args...)
Base.hcat(args::AbstractNode...) = node(hcat, args...)

Statistics.mean(X::AbstractNode)   = node(v->mean.(v), X)
Statistics.median(X::AbstractNode) = node(v->median.(v), X)
StatsBase.mode(X::AbstractNode)    = node(v->mode.(v), X)

Base.log(X::AbstractNode) = node(v->log.(v), X)
Base.exp(X::AbstractNode) = node(v->exp.(v), X)

+(y1::AbstractNode, y2::AbstractNode) = node(+, y1, y2)
+(y1, y2::AbstractNode) = node(+, y1, y2)
+(y1::AbstractNode, y2) = node(+, y1, y2)

*(lambda::Real, y::AbstractNode) = node(y->lambda*y, y)

"""
    selectcols(X::AbstractNode, c)

Returns `Node` object `N` such that `N() = selectcols(X(), c)`.
"""
MMI.selectcols(X::AbstractNode, r) = node(XX->selectcols(XX, r), X)

"""
    selectrows(X::AbstractNode, r)

Returns a `Node` object `N` such that `N() = selectrows(X(), r)` (and
`N(rows=s) = selectrows(X(rows=s), r)`).

"""
MMI.selectrows(X::AbstractNode, r) = node(XX->selectrows(XX, r), X)

# for accessing and setting model hyperparameters at node:
getindex(n::Node{<:NodalMachine{<:Model}}, s::Symbol) =
    getproperty(n.machine.model, s)
setindex!(n::Node{<:NodalMachine{<:Model}}, v, s::Symbol) =
    setproperty!(n.machine.model, s, v)

## INSPECTING LEARNING NETWORKS

"""
    tree(N)

Return a description of the tree `N` defined by the learning network
terminating at a given node.

"""
function tree(W::Node)
    mach = W.machine
    if mach === nothing
        value2 = nothing
        endkeys = []
        endvalues = []
    else
        value2 = mach.model
        endkeys = (Symbol("train_arg", i) for i in eachindex(mach.args))
        endvalues = (tree(arg) for arg in mach.args)
    end
    keys = tuple(:operation,  :model,
                 (Symbol("arg", i) for i in eachindex(W.args))...,
                 endkeys...)
    values = tuple(W.operation, value2,
                   (tree(arg) for arg in W.args)...,
                   endvalues...)
    return NamedTuple{keys}(values)
end
tree(s::Source) = (source = s,)

"""
   args(tree; train=false)

Return a vector of the top level args of the tree associated with a node.
If `train=true`, return the `train_args`.
"""
function args(tree; train=false)
    keys_ = filter(keys(tree) |> collect) do key
        match(Regex("^$("train_"^train)arg[0-9]*"), string(key)) !== nothing
    end
    return [getproperty(tree, key) for key in keys_]
end

"""
    models(N::AbstractNode)

A vector of all models referenced by a node `N`, each model appearing
exactly once.

"""
function models(W::AbstractNode)
    models_ = filter(flat_values(tree(W)) |> collect) do model
        model isa Model
    end
    return unique(models_)
end

"""
    sources(N::AbstractNode; kind=:any)

A vector of all sources referenced by calls `N()` and `fit!(N)`. These
are the sources of the directed acyclic graph associated with the
learning network terminating at `N`, including training edges. The
return value can be restricted further by specifying `kind=:input`,
`kind=:target`, `kind=:weight`, etc.

Not to be confused with `origins(N)` which refers to the same graph with edges
corresponding to training arguments deleted.

See also: [`origins`](@ref), [`source`](@ref).
"""
function sources(W::AbstractNode; kind=:any)
    if kind == :any
        sources_ = filter(flat_values(tree(W)) |> collect) do value
            value isa Source
        end
    else
        sources_ = filter(flat_values(tree(W)) |> collect) do value
            value isa Source{kind}
        end
    end
    return unique(sources_)
end

"""
    machines(N)

List all machines in the learning network terminating at node `N`.

"""
function machines(W::Node)
    if W.machine === nothing
        return vcat((machines(arg) for arg in W.args)...) |> unique
    else
        return vcat(Any[W.machine, ],
                    (machines(arg) for arg in W.args)...,
                    (machines(arg) for arg in W.machine.args)...) |> unique
    end
end
machines(W::Source) = Any[]

## MANIPULATING LEARNING NETWORKS

"""
    reset!(N::Node)

Place the learning network terminating at node `N` into a state in
which `fit!(N)` will retrain from scratch all machines in its
dependency tape. Does not actually train any machine or alter
fit-results. (The method simply resets `m.state` to zero, for every
machine `m` in the network.)

"""
function reset!(W::Node)
    for mach in machines(W)
        mach.state = 0 # to do: replace with dagger object
    end
end


text:
using Luxor

origin = Point(0, 0)
scale = 1
r = 90 * scale
c1 = Point(0, -r) + origin
c2 = Point(-r * 3^0.5 / 2, r / 2) + origin
c3 = Point(0, r) + origin
c4 = Point(r * 3^0.5 / 2, r / 2) + origin
box_size = 15 * scale
circle_size = 30 * scale

box_node1 = Point(0, -box_size) + origin
box_node2 = Point(-box_size, box_size) + origin
box_node3 = Point(0, box_size) + origin
box_node4 = Point(box_size, box_size) + origin

c_node1 = c1 + Point(0, circle_size)
c_node2 = c2 + Point(circle_size * 3^0.5 / 2, -circle_size / 2)
c_node3 = c3 + Point(0, -circle_size)
c_node4 = c4 + Point(-circle_size * 3^0.5 / 2, -circle_size / 2)

arc_node1 = box_node2 + (c_node2 - box_node2) * 0.2
arc_node2 = box_node4 + (c_node4 - box_node4) * 0.2

arrow_width = 5 * scale
aha = pi / 4

Drawing(320, 320, "docs/src/assets/logo.svg")
background(1, 1, 1, 0)
Luxor.origin()

setcolor(0.0, 0.0, 0.0) # black
box(Point(-box_size, -box_size) + origin, Point(box_size, box_size) + origin, :fill)
arc2r(origin, arc_node2, arc_node1, :stroke)
setcolor(0.251, 0.388, 0.847)  # dark blue
circle(c1, circle_size, :fill)
arrow(c_node1, box_node1, linewidth = arrow_width, arrowheadangle = aha)
setcolor(0.796, 0.235, 0.2)  # dark red
circle(c2, circle_size, :fill)
arrow(box_node2, c_node2, linewidth = arrow_width, arrowheadangle = aha)
setcolor(0.22, 0.596, 0.149) # dark green
circle(c3, circle_size, :fill)
arrow(box_node3, c_node3, linewidth = arrow_width, arrowheadangle = aha)
setcolor(0.584, 0.345, 0.698) # dark purple
circle(c4, circle_size, :fill)
arrow(box_node4, c_node4, linewidth = arrow_width, arrowheadangle = aha)

finish()


text:
"""
Class for Dispatching Forward Payoff

		forward=Forward(T::num) where {num <: Number}

Where:\n
		T	=	Time to maturity of the Forward.
"""
mutable struct Forward{num <: Number} <: EuropeanPayoff{num}
    T::num
    function Forward(T::num) where {num <: Number}
        if T <= 0.0
            error("Time to Maturity must be positive")
        else
            return new{num}(T)
        end
    end
end

export Forward;

function payoff(S::AbstractMatrix{num}, optionData::Forward, rfCurve::AbstractZeroRateCurve, T1::num2 = maturity(optionData)) where {num <: Number, num2 <: Number}
    r = rfCurve.r
    T = optionData.T
    NStep = size(S, 2) - 1
    index1 = round(Int, T / T1 * NStep) + 1
    @views ST = S[:, index1]

    return ST * exp(-integral(r, T))
end


text:
##############################################################################
##
## Functions related to the interface with PrettyTables.jl.
##
##############################################################################

# Default DataFrames highlighter for text backend.
#
# This highlighter changes the text color to gray in cells with `nothing`,
# `missing`, `#undef`, and types related to DataFrames.jl.
function _pretty_tables_highlighter_func(data, i::Integer, j::Integer)
    try
        cell = data[i, j]
        return ismissing(cell) ||
            cell === nothing ||
            cell isa Union{AbstractDataFrame, GroupedDataFrame,
                           DataFrameRow, DataFrameRows,
                           DataFrameColumns}
    catch e
        if isa(e, UndefRefError)
            return true
        else
            rethrow(e)
        end
    end
end

const _PRETTY_TABLES_HIGHLIGHTER = Highlighter(_pretty_tables_highlighter_func,
                                               Crayon(foreground = :dark_gray))

# Default DataFrames formatter for text backend.
#
# This formatter changes how the following types are presented when rendering
# the data frame:
#     - missing;
#     - nothing;
#     - Cells with types related to DataFrames.jl.

function _pretty_tables_general_formatter(v, i::Integer, j::Integer)
    if typeof(v) <: Union{AbstractDataFrame, GroupedDataFrame, DataFrameRow,
                          DataFrameRows, DataFrameColumns}

        # Here, we must not use `print` or `show`. Otherwise, we will call
        # `_pretty_table` to render the current table leading to a stack
        # overflow.
        return sprint(summary, v)
    elseif ismissing(v)
        return "missing"
    elseif v === nothing
        return ""
    else
        return v
    end
end

# Formatter to align the floating points as in Julia array printing.
#
# - `float_cols` contains the IDs of the columns that must be formatted.
# - `indices` is a vector of vectors containing the indices of each elements
#   in the data frame.
# - `padding` is a vector of vectors containing the padding of each element for
#   each row.
# - `compact_printing` must be a boolean indicating if we should enable the
#   `:compact` option of `io` when converting the number to string.

function _pretty_tables_float_formatter(v, i::Integer, j::Integer,
                                        float_cols::Vector{Int},
                                        align_col::Vector{Int},
                                        compact_printing::Bool)
    isempty(float_cols) && return v

    # We apply this formatting only to the columns that contains only floats.
    ind_col = findfirst(==(j), float_cols)

    if ind_col !== nothing
        align_col_i = align_col[ind_col]

        # Convert the value to text.
        str = sprint(print, v, context = :compact => compact_printing)

        if v == "missing"
            pad = align_col_i - 8
        else
            # We want to align everything at '.'.
            id_dp = findfirst('.', str)

            # Compute the require padding to align the cell.
            #
            # If a decimal point is not found, then assume that the
            # entire text should be aligned before the alignment column. This
            # can happen with a custom `AbstractFloat` type.
            pad = id_dp !== nothing ? align_col_i - id_dp :
                                      align_col_i - (textwidth(str) + 1)
        end

        # Pad cannot be negative.
        # NOTE: This is just a failsafe, `pad` will be negative only if we have
        # a bug. However, in this case, we can avoid an error when printing.
        pad < 0 && (pad = 0)

        return " "^pad * str
    end

    # The formatter is applied to all tables' cells. Hence, we must return the
    # input value `v` unchanged if this cell is not part of a column that has
    # floating point numbers.
    return v
end


text:
module EMLflash
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...stage_separators.EMLtank
	import EMLtypes.length
	include("flash/flash.jl")
	include("flash/flash_steady.jl")
	include("flash/FlashPHSteady.jl")
end


text:
struct PolicyAcc
    total::Int64
    failed::Int64
end

getfail(x::PolicyAcc) = x.failed
getpass(x::PolicyAcc) = x.total - x.failed

function getpct(x::PolicyAcc)
    round(getpass(x) / x.total, digits=6)
end


text:
using WebSockets, Test, Base64, Random
import Base: BufferStream, convert
import WebSockets.HTTP, Sockets, Dates
import WebSockets.MbedTLS
import Sockets.@ip_str

@testset "WebSockets" begin
    include("logformat.jl")
    include("client_server_functions.jl")

    printstyled(color=:blue, "\nBase.Show\n")
    @testset "Base.show" begin
       include("show_test.jl");sleep(1)
    end

    printstyled(color=:blue, "\nWebsocketLogger\n")
    @testset "WebSocketLogger" begin
       include("test_websocketlogger.jl");sleep(1)
    end

    printstyled(color=:blue, "\nFragment and unit\n")
    @testset "Fragment and unit" begin
       include("frametest.jl");sleep(1)
    end

    printstyled(color=:blue, "\nHandshake\n")
    @testset "Handshake" begin
        include("handshaketest.jl");sleep(1)
    end

    printstyled(color=:blue, "\nClient test\n")
    @testset "Client" begin
        include("client_test.jl");sleep(1)
    end

    printstyled(color=:blue, "\nClient_listen\n")
    @testset "Client_listen" begin
        include("client_listen_test.jl");sleep(1)
    end

    printstyled(color=:blue, "\nClient_serverWS\n")
    @testset "Client_serverWS" begin
        include("client_serverWS_test.jl");sleep(1)
    end

    printstyled(color=:blue, "\nAbrupt close & error handling\n")
    @testset "Abrupt close & error handling" begin
       include("error_test.jl");sleep(1)
    end
    if !@isdefined(OLDLOGGER)
        Logging.global_logger(OLDLOGGER)
    end
end


text:
using KitBase, KitBase.FastGaussQuadrature
using ProgressMeter, Plots

begin
    # case
    matter = "photon"
    case = "linesource"
    space = "1d1f1v"
    nSpecies = 1
    flux = "kfvs"
    collision = "bgk"
    interpOrder = 2
    limiter = "vanleer"
    boundary = "fix"
    cfl = 0.5
    maxTime = 0.2

    # physical space
    x0 = 0
    x1 = 1
    nx = 100
    pMeshType = "uniform"
    nxg = 0

    # velocity space
    vMeshType = "rectangle"
    umin = -1
    umax = 1
    nu = 28
    nug = 0

    # gas
    knudsen = 0.01
    mach = 0.0
    prandtl = 1
    inK = 0
    omega = 0.5
    alphaRef = 1.0
    omegaRef = 0.5
end

set = Setup(
    matter,
    case,
    space,
    flux,
    collision,
    nSpecies,
    interpOrder,
    limiter,
    boundary,
    cfl,
    maxTime,
)

pSpace = PSpace1D(x0, x1, nx, nxg)

points, weights = gausslegendre(nu)
vSpace = VSpace1D(points[1], points[end], nu, points, ones(nu).*(points[end]-points[1])/(nu-1), weights)

μᵣ = ref_vhs_vis(knudsen, alphaRef, omegaRef)
gas = Gas(knudsen, mach, prandtl, inK, 3.0, omega, alphaRef, omegaRef, μᵣ)

fL = 0.5 * ones(nu)
wL = moments_conserve(fL, points, weights)
primL = conserve_prim(wL, 3.0)
fR = 0.0001 * ones(nu)
wR = moments_conserve(fR, points, weights)
primR = conserve_prim(wR, 3.0)
ib = IB1F(wL, primL, fL, primL, wR, primR, fR, primR)

ks = SolverSet(set, pSpace, vSpace, gas, ib, @__DIR__)

ctr = Array{ControlVolume1D1F}(undef, ks.pSpace.nx)
face = Array{Interface1D1F}(undef, ks.pSpace.nx + 1)
for i in eachindex(ctr)
    ctr[i] = ControlVolume1D1F(
        ks.pSpace.x[i],
        ks.pSpace.dx[i],
        ks.ib.wR,
        ks.ib.primR,
        ks.ib.fR,
    )
end

for i = 1:ks.pSpace.nx+1
    face[i] = Interface1D1F(ks.ib.wL, ks.ib.fL)
end

function fboundary(
    ff::T2,
    f::T4,
    u::T5,
    dt,
    rot = 1,
) where {
    T2<:AbstractArray{<:AbstractFloat,1},
    T4<:AbstractArray{<:AbstractFloat,1},
    T5<:AbstractArray{<:AbstractFloat,1},
}

    δ = heaviside.(u .* rot)

    fWall = 0.5 .* δ .+ f .* (1.0 .- δ)

    @. ff = u * fWall * dt

    return nothing

end

function step(
    ffL::T2,
    w::T3,
    f::T4,
    ffR::T2,
    u::T5,
    weights::T5,
    τ,
    dx,
    dt,
) where {
    T2<:AbstractArray{<:AbstractFloat,1},
    T3<:AbstractArray{<:AbstractFloat,1},
    T4<:AbstractArray{<:AbstractFloat,1},
    T5<:AbstractArray{<:AbstractFloat,1},
}
    M = sum(weights .* f)
    for i in eachindex(u)
        f[i] += (ffL[i] - ffR[i]) / dx + (M - f[i]) / τ * dt
    end
    w[1] = sum(weights .* f)
end

dt = cfl * pSpace.dx[1]
nt = maxTime / dt |> floor |> Int

anim = @animate for iter = 1:150#nt
    #reconstruct!(ks, ctr)

    fboundary(
        face[1].ff,
        ctr[1].f,
        vSpace.u,
        dt,
    )
    @inbounds for i = 2:nx
        flux_kfvs!(
            face[i].ff,
            ctr[i-1].f,
            ctr[i].f,
            vSpace.u,
            dt,
            ctr[i-1].sf,
            ctr[i].sf,
        )
    end

    @inbounds for i = 1:nx-1
        step(
            face[i].ff,
            ctr[i].w,
            ctr[i].f,
            face[i+1].ff,
            vSpace.u,
            vSpace.weights,
            1.0,
            pSpace.dx[i],
            dt,
        )
    end

    pltx = ks.pSpace.x[1:ks.pSpace.nx]
    plty = zeros(ks.pSpace.nx, 6)
    for i in eachindex(pltx)
        plty[i, 1] = ctr[i].w[1]
    end
    plot(pltx, plty[:, 1], label = "Density", lw = 2, xlabel = "x")
end

gif(anim, "propagation.gif")

begin
    pltx = ks.pSpace.x[1:ks.pSpace.nx]
    plty = zeros(ks.pSpace.nx, 6)
    for i in eachindex(pltx)
        plty[i, 1] = ctr[i].w[1]
    end
    plot(pltx, plty[:, 1], label = "Density", lw = 2, xlabel = "x")
end


text:
function field_getter(
    obuf,
    prefix :: Symbol,
    struct_type,
    map :: Dict{Symbol, Any})

    for key in keys(map)
        jl_name = symbol(prefix, key)
        field_type, field_offset = map[key]
        field_offset = convert(UInt64, field_offset)

        getter_code = quote
            function $jl_name(ptr::Ptr{$struct_type})
                field_ptr = convert(Ptr{$field_type}, ptr + $field_offset)
                return unsafe_load(field_ptr)
               end
        end

        # Debug output of generated code
        if false && endswith(string(jl_name), "size")
            println(getter_code)
        end

        push!(obuf, getter_code)
    end
end

function gen_fieldgetters(obuf)

    field_getter(
        obuf,
        :queue_info_,
        :hsa_queue_t,
        Dict(
            :type => (:hsa_queue_type_t, 0),
            :features => (:hsa_queue_feature_t, 4),
            :base_address => (UInt64, 8), # TODO handle ifdefs for machine_model/endianness
            :doorbell_signal => (UInt64, 16),
            :size => (UInt32, 24),
            # UInt32 reserved
            :id => (UInt64, 32),
        )
    )
end


text:
type Tokenizer
    dict::IdDict
    tagset::Tagset
    model
end

function Tokenizer()
    dict = IdDict(String["UNKNOWN", " ", "\n"])
    T = Float32
    x = GraphNode()
    x = window(x, (7,), pad=3)
    x = GraphNode(constant, x)
    x = Embedding(T,300,10)(x)
    x = Linear(T,70,70)(x)
    x = relu(x)
    x = Linear(T,70,3)(x)
    f = compile(x)
    Tokenizer(dict, IOE(), f)
end

function (t::Tokenizer)(chars::Vector{Char})
    unk = t.dict["UNKNOWN"]
    x = map(chars) do c
        get(t.dict, string(c), unk)
    end
    y = t.model(x).data
    tags = argmax(y,1)
    decode(t.tagset, tags)
end
(t::Tokenizer)(str::String) = t(Vector{Char}(str))

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
using Distributed
using ArgParse
using RoME, IncrementalInference, DistributedFactorGraphs
@everywhere using RoME, IncrementalInference, DistributedFactorGraphs

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "files"
            help = "Txt describing Which .tar.gz to resolve"
            required = true
    end

    return parse_args(s)
end

pargs = parse_commandline()

fid = open(pargs["files"], "r")
loadFiles = readlines(fid)
close(fid)

for thefile in loadFiles

fg = initfg()
loadDFG(thefile, Main, fg)
# loadDFG(joinpath(@__DIR__,"fg_x37.tar.gz"), Main, fg)

dontMarginalizeVariablesAll!(fg)

Base.mv(thefile, thefile*"_moved", force=true)
try
  @show imgfile = split(splitpath(thefile)[end], '.')[1]
  Base.mv(joinpath("background", imgfile*".png"), joinpath("background", "old", imgfile*".png"), force=true)
catch
  #
end

# solve only last segment
vars = ls(fg, r"x\d", solvable=1) |> sortDFG
nrSolved = floor(Int, length(vars)/10)*10+1
nrSt = maximum([nrSolved-50, 1])

vars = vars[nrSt:nrSolved]
union!(vars, [:l1])

dontSolve = setdiff(ls(fg), vars)

(x->setSolvable!(fg, x, 0)).(dontSolve)

## Solve and look again

getSolverParams(fg).drawtree=true

tree = solveTree!(fg)

##

## resave with fixed lag?
# defaultFixedLagOnTree!(fg,30,limitfixeddown=false)

saveDFG(fg, thefile)

end

#
#
# plotVariableBeliefs(sfg, r"x\d")[1]
#

#