text:
using Test
using StaticArrays, IntervalSets, LinearAlgebra

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.CartesianVector.(∇ᶜ.(sin.(faces)))
        @test ∂sin ≈ Geometry.Cartesian3Vector.(cos.(centers)) atol = 1e-2

        divᶜ = Operators.DivergenceF2C()
        ∂sin = divᶜ.(Geometry.Cartesian3Vector.(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.CartesianVector.(∇ᶠ.(centers))
        @test ∂z ≈ Geometry.Cartesian3Vector.(ones(FT, face_space)) rtol =
            10 * eps(FT)

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

        ∇ᶠ = Operators.GradientC2F(
            left = Operators.SetGradient(Geometry.Cartesian3Vector(FT(0))),
            right = Operators.SetGradient(Geometry.Cartesian3Vector(FT(0))),
        )
        ∂cos = Geometry.CartesianVector.(∇ᶠ.(cos.(centers)))
        @test ∂cos ≈ Geometry.Cartesian3Vector.(.-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.CartesianVector.(∂.(w .* I.(θ)))
        @test ∂sin ≈ Geometry.Cartesian3Vector.(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.CartesianVector.(∂.(w .* I.(θ)))
        @test ∂sin ≈ Geometry.Cartesian3Vector.(cos.(centers)) atol = 1e-2

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

        ∂sin = Geometry.CartesianVector.(∂.(w .* I.(θ)))
        @test ∂sin ≈ Geometry.Cartesian3Vector.(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.CartesianVector.(∂.(w .* I.(θ)))
        @test ∂sin ≈ Geometry.Cartesian3Vector.(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.Cartesian3Vector.(3π .* cos.(3π .* faces.z)),
                )

            operator = Operators.GradientC2F(
                left = Operators.SetGradient(Geometry.Cartesian3Vector(3π)),
                right = Operators.SetGradient(Geometry.Cartesian3Vector(-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.Cartesian3Vector.(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))
    Δ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.CartesianVector.(∇ᶜ.(sin.(faces)))

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

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

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

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

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

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

        Δh[k] = cs.face_local_geometry.J[1]
        # Errors
        err_grad_sin_c[k] =
            norm(gradsinᶜ .- Geometry.Cartesian3Vector.(cos.(centers)))
        err_div_sin_c[k] = norm(divsinᶜ .- cos.(centers))
        err_grad_z_f[k] = norm(∂zᶠ .- Geometry.Cartesian3Vector.(ones(FT, fs)))
        err_grad_cos_f1[k] =
            norm(gradcosᶠ¹ .- Geometry.Cartesian3Vector.(.-sin.(faces)))
        err_grad_cos_f2[k] =
            norm(gradcosᶠ² .- Geometry.Cartesian3Vector.(.-sin.(faces)))
        err_div_sin_f[k] = norm(
            divsinᶠ .-
            (Geometry.Cartesian3Vector.(cos.(faces))).components.data.:1,
        )
        err_div_cos_f[k] = norm(
            divcosᶠ .-
            (Geometry.Cartesian3Vector.(.-sin.(faces))).components.data.:1,
        )
    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)

    # 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]

end


text:
# This file contains the basic functions of PCG.
import RNG: AbstractRNG

"Return the default multiplier for a certain type."
@inline default_multiplier(::Type{UInt8}) = 0x8d
@inline default_multiplier(::Type{UInt16}) = 0x321d
@inline default_multiplier(::Type{UInt32}) = 0x2c9277b5
@inline default_multiplier(::Type{UInt64}) = 0x5851f42d4c957f2d
@inline default_multiplier(::Type{UInt128}) = 0x2360ed051fc65da44385df649fccf645

"Return the default increment for a certain type."
@inline default_increment(::Type{UInt8}) = 0x4d
@inline default_increment(::Type{UInt16}) = 0xbb75
@inline default_increment(::Type{UInt32}) = 0xac564b05
@inline default_increment(::Type{UInt64}) = 0x14057b7ef767814f
@inline default_increment(::Type{UInt128}) = 0x5851f42d4c957f2d14057b7ef767814f

"Return the default MCG multiplier for a certain type."
@inline mcg_multiplier(::Type{UInt8}) = 0xd9
@inline mcg_multiplier(::Type{UInt16}) = 0xf2d9
@inline mcg_multiplier(::Type{UInt32}) = 0x108ef2d9
@inline mcg_multiplier(::Type{UInt64}) = 0xaef17502108ef2d9
@inline mcg_multiplier(::Type{UInt128}) = 0xf69019274d7f699caef17502108ef2d9

"Return the default MCG unmultiplier for a certain type."
@inline mcg_unmultiplier(::Type{UInt8}) = 0x69
@inline mcg_unmultiplier(::Type{UInt16}) = 0x6d69
@inline mcg_unmultiplier(::Type{UInt32}) = 0xacb86d69
@inline mcg_unmultiplier(::Type{UInt64}) = 0xd04ca582acb86d69
@inline mcg_unmultiplier(::Type{UInt128}) = 0xc827645e182bc965d04ca582acb86d69

@inline uint_index(::Type{UInt8}) = 1
@inline uint_index(::Type{UInt16}) = 2
@inline uint_index(::Type{UInt32}) = 3
@inline uint_index(::Type{UInt64}) = 4
@inline uint_index(::Type{UInt128}) = 5

@inline half_width(::Type{UInt16}) = UInt8
@inline half_width(::Type{UInt32}) = UInt16
@inline half_width(::Type{UInt64}) = UInt32
@inline half_width(::Type{UInt128}) = UInt64

# Shift and Rotate functions

# Rotate functions.
@inline function pcg_rotr{T<:PCGUInt}(value::T, rot::T)
    s = sizeof(T) << 3
    (value >> (rot % s)) | (value << (-rot % s))
end

"General advance functions."
@inline function pcg_advance_lcg{T<:PCGUInt}(state::T, delta::T, cur_mult::T, cur_plus::T)
    acc_mult = 1 % T
    acc_plus = 0 % T
    while delta > 0
        if delta & 1 == 1
            acc_mult *= cur_mult
            acc_plus = acc_plus * cur_mult + cur_plus
        end
        cur_plus = (cur_mult + (1 % T)) * cur_plus
        cur_mult *= cur_mult
        delta = delta >> 1
    end
    acc_mult * state + acc_plus
end

# output_xsh_rs
# Xorshift High, Random Shift.
# return half bits of T.
@inline function pcg_output{T<:Union{pcg_uints[2:end]...}}(state::T, ::Type{PCG_XSH_RS})
    return_bits = sizeof(T) << 2
    bits = return_bits << 1
    spare_bits = bits - return_bits
    op_bits = spare_bits - 5 >= 64 ? 5 :
              spare_bits - 4 >= 32 ? 4 :
              spare_bits - 3 >= 16 ? 3 :
              spare_bits - 2 >= 4  ? 2 :
              spare_bits - 1 >= 1  ? 1 : 0
    mask = (1 << op_bits) - 1
    xshift = op_bits + (return_bits + mask) >> 1
    rshift = op_bits != 0 ? (state >> (bits - op_bits)) & mask : 0 % T
    state = state $ (state >> xshift)
    (state >> (spare_bits - op_bits - mask + rshift)) % half_width(T)
end

# output_xsh_rr
# Xorshift High, Random Rotation.
# return half bits of T.
@inline function pcg_output{T<:Union{pcg_uints[2:end]...}}(state::T, ::Type{PCG_XSH_RR})
    return_bits = sizeof(T) << 2
    bits = return_bits << 1
    spare_bits = bits - return_bits
    op_bits = return_bits >= 128 ? 7 :
              return_bits >=  64 ? 6 :
              return_bits >=  32 ? 5 :
              return_bits >=  16 ? 4 : 3
    mask = (1 << op_bits) - 1
    xshift = (op_bits + return_bits) >> 1
    rot = (state >> (bits - op_bits)) & mask
    state $= (state >> xshift)
    result = (state >> (spare_bits - op_bits)) % half_width(T)
    pcg_rotr(result, rot % half_width(T))
end

# output_rxs_m_xs
# Random Xorshift, Multiplication, Xorshift.
# return the same bits as T.
# Insecure.
@inline function pcg_output{T<:PCGUInt}(state::T, ::Type{PCG_RXS_M_XS})
    bits = return_bits = sizeof(T) << 3
    op_bits = return_bits >= 128 ? 6 :
              return_bits >=  64 ? 5 :
              return_bits >=  32 ? 4 :
              return_bits >=  16 ? 3 : 2
    mask = (1 << op_bits) - 1
    rshift = (state >> (bits - op_bits)) & mask
    state $= state >> (op_bits + rshift)
    state *= mcg_multiplier(T)
    state $ (state >> ((return_bits << 1 + 2) ÷ 3))
end

# output_xsl_rr
# Xorshift Low, Random Rotation.
# return half bits of T.
@inline function pcg_output{T<:Union{UInt64, UInt128}}(state::T, ::Type{PCG_XSL_RR})
    return_bits = sizeof(T) << 2
    bits = return_bits << 1
    spare_bits = bits - return_bits
    op_bits = return_bits >= 128 ? 7 :
              return_bits >=  64 ? 6 :
              return_bits >=  32 ? 5 :
              return_bits >=  16 ? 4 : 3
    mask = (1 << op_bits) - 1
    xshift = (spare_bits + return_bits) >> 1
    rot = (state >> (bits - op_bits)) & mask
    pcg_rotr((state $ (state >> xshift)) % half_width(T), rot % half_width(T))
end

# output_xsl_rr_rr
# Xorshift Low, Random Rotation, Random Rotation.
# return the same bits as T.
# Insecure.
@inline function pcg_output{T<:Union{UInt64, UInt128}}(state::T, ::Type{PCG_XSL_RR_RR})
    half_bits = sizeof(T) << 2
    bits = half_bits << 1
    spare_bits = bits - half_bits
    op_bits = half_bits >= 128 ? 7 :
              half_bits >=  64 ? 6 :
              half_bits >=  32 ? 5 :
              half_bits >=  16 ? 4 : 3
    mask = (1 << op_bits) - 1
    xshift = (spare_bits + half_bits) >> 1
    rot = (state >> (bits - op_bits)) & mask
    state $= state >> xshift
    low_bits = state % half_width(T)
    low_bits = pcg_rotr(low_bits, rot % half_width(T))
    high_bits = (state >> spare_bits) % half_width(T)
    rot2 = low_bits & mask
    high_bits = pcg_rotr(high_bits, rot2 % half_width(T))
    ((high_bits % T) << spare_bits) $ (low_bits % T)
end

# PCGState types, SRandom and Step functions.
"""
```julia
AbstractPCG{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <: AbstractRNG{OutputType}
```

The base abstract type for PCGs.
"""
abstract AbstractPCG{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <:
    AbstractRNG{OutputType}

# pcg_state_XX
# XX is one of the UInt types.
type PCGStateOneseq{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <:
        AbstractPCG{StateType, MethodType, OutputType}
    state::StateType
    PCGStateOneseq() = new()
end

@inline function pcg_srand{T<:PCGUInt}(s::PCGStateOneseq{T}, init_state::T)
    s.state = 0 % T
    pcg_step!(s)
    s.state += init_state
    pcg_step!(s)
    s
end

# pcg_oneseq_XX_step_r
# XX is one of the UInt types.
@inline function pcg_step!{T<:PCGUInt}(s::PCGStateOneseq{T})
    s.state = s.state * default_multiplier(T) + default_increment(T)
end

# pcg_oneseq_XX_advance_r
# XX is one of the UInt types.
@inline function pcg_advance!{T<:PCGUInt}(s::PCGStateOneseq{T}, delta::T)
    s.state = pcg_advance_lcg(s.state, delta, default_multiplier(T), default_increment(T))
end

# pcg_state_XX
# XX is one of the UInt types.
type PCGStateMCG{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <:
        AbstractPCG{StateType, MethodType, OutputType}
    state::StateType
    PCGStateMCG() = new()
end

@inline function pcg_srand{T<:PCGUInt}(s::PCGStateMCG{T}, init_state::T)
    s.state = init_state | 1
    s
end

# pcg_mcg_XX_step_r
# XX is one of the UInt types.
@inline function pcg_step!{T<:PCGUInt}(s::PCGStateMCG{T})
    s.state = s.state * default_multiplier(T)
end

# pcg_mcg_XX_advance_r
# XX is one of the UInt types.
@inline function pcg_advance!{T<:PCGUInt}(s::PCGStateMCG{T}, delta::T)
    s.state = pcg_advance_lcg(s.state, delta, default_multiplier(T), 0 % T)
end

# pcg_state_XX
# XX is one of the UInt types.
type PCGStateUnique{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <:
        AbstractPCG{StateType, MethodType, OutputType}
    state::StateType
    PCGStateUnique() = new()
end

@inline function pcg_srand{T<:PCGUInt}(s::PCGStateUnique{T}, init_state::T)
    s.state = 0 % T
    pcg_step!(s)
    s.state += init_state
    pcg_step!(s)
    s
end

# pcg_unique_XX_step_r
# XX is one of the UInt types.
@inline function pcg_step!{T<:PCGUInt}(s::PCGStateUnique{T})
    s.state = s.state * default_multiplier(T) + (UInt(pointer_from_objref(s)) | 1) % T
end

# pcg_unique_XX_advance_r
# XX is one of the UInt types.
@inline function pcg_advance!{T<:PCGUInt}(s::PCGStateUnique{T}, delta::T)
    s.state = pcg_advance_lcg(s.state, delta, default_multiplier(T),
        (UInt(pointer_from_objref(s)) | 1) % T)
end

# pcg_state_XX
# XX is one of the UInt types.
type PCGStateSetseq{StateType<:PCGUInt, MethodType<:PCGMethod, OutputType<:PCGUInt} <:
        AbstractPCG{StateType, MethodType, OutputType}
    state::StateType
    inc::StateType
    PCGStateSetseq() = new()
end

@inline function pcg_srand{T<:PCGUInt}(s::PCGStateSetseq{T}, init_state::T, init_seq::T)
    s.state = 0
    s.inc = (init_seq << 1) | 1
    pcg_step!(s)
    s.state += init_state
    pcg_step!(s)
    s
end

# pcg_setseq_XX_step_r
# XX is one of the UInt types.
@inline function pcg_step!{T<:PCGUInt}(s::PCGStateSetseq{T})
    s.state = s.state * default_multiplier(T) + s.inc
end

# pcg_setseq_XX_advance_r
# XX is one of the UInt types.
@inline function pcg_advance!{T<:PCGUInt}(s::PCGStateSetseq{T}, delta::T)
    s.state = pcg_advance_lcg(s.state, delta, default_multiplier(T), s.inc)
end

"Return the output of a state for a certain PCG type."
pcg_output

"Initialize a PCG object."
pcg_srand

"Do one iteration step for a PCG object."
pcg_step!

"Advance a PCG object."
pcg_advance!

ext: jl


text:
@testset "1_bfv_basics" begin
  @testset "EncryptionParameters" begin
    @test_nowarn EncryptionParameters(SchemeType.bfv)
  end
  enc_parms = EncryptionParameters(SchemeType.bfv)

  @testset "polynomial modulus degree" begin
    @test_nowarn set_poly_modulus_degree!(enc_parms, 4096)
  end

  @testset "coefficient modulus" begin
    @test_nowarn coeff_modulus_bfv_default(4096)
    @test_nowarn set_coeff_modulus!(enc_parms, coeff_modulus_bfv_default(4096))
  end

  @testset "plain modulus" begin
    @test_nowarn set_plain_modulus!(enc_parms, 1024)
  end

  @testset "SEALContext" begin
    @test_nowarn SEALContext(enc_parms)
  end
  context = SEALContext(enc_parms)

  @testset "parameter_error_message" begin
    @test parameter_error_message(context) == "valid"
  end

  @testset "KeyGenerator" begin
    @test_nowarn KeyGenerator(context)
  end
  keygen = KeyGenerator(context)

  @testset "PublicKey" begin
    @test_nowarn PublicKey()
  end
  public_key_ = PublicKey()

  @testset "create_public_key" begin
    @test_nowarn create_public_key!(public_key_, keygen)
  end

  @testset "SecretKey" begin
    @test_nowarn secret_key(keygen)
  end
  secret_key_ = secret_key(keygen)

  @testset "RelinKeys" begin
    @test_nowarn RelinKeys()
  end
  relin_keys_ = RelinKeys()

  @testset "create_relin_keys" begin
    @test_nowarn create_relin_keys!(relin_keys_, keygen)
  end

  @testset "Encryptor" begin
    @test_nowarn Encryptor(context, public_key_)
  end
  encryptor = Encryptor(context, public_key_)

  @testset "Evaluator" begin
    @test_nowarn Evaluator(context)
  end
  evaluator = Evaluator(context)

  @testset "Decryptor" begin
    @test_nowarn Decryptor(context, secret_key_)
  end
  decryptor = Decryptor(context, secret_key_)

  @testset "Plaintext" begin
    @test_nowarn Plaintext(string(6))
  end
  x_plain = Plaintext(string(6))

  @testset "to_string" begin
    @test to_string(x_plain) == "6"
  end

  x_encrypted = Ciphertext()
  @testset "encrypt!" begin
    @test_nowarn encrypt!(x_encrypted, x_plain, encryptor)
  end

  @testset "length" begin
    @test length(x_encrypted) == 2
  end

  @testset "invariant_noise_budget" begin
    # Next test is fuzzy since the actual noise budget might vary
    @test invariant_noise_budget(x_encrypted, decryptor) in (54, 55, 56)
  end

  x_decrypted = Plaintext()
  @testset "decrypt!" begin
    @test_nowarn decrypt!(x_decrypted, x_encrypted, decryptor)
    @test to_string(x_decrypted) == "6"
  end

  x_sq_plus_one = Ciphertext()
  @testset "square!" begin
    @test_nowarn square!(x_sq_plus_one, x_encrypted, evaluator)
  end

  plain_one = Plaintext("1")
  @testset "add_plain_inplace! and length/noise budget" begin
    @test_nowarn add_plain_inplace!(x_sq_plus_one, plain_one, evaluator)
    @test length(x_sq_plus_one) == 3
    @test invariant_noise_budget(x_sq_plus_one, decryptor) == 33
  end

  decrypted_result = Plaintext()
  @testset "decrypt! and check (x^2 + 1 = 37 = 0x25)" begin
    @test_nowarn decrypt!(decrypted_result, x_sq_plus_one, decryptor)
    @test to_string(decrypted_result) == "25"
  end

  x_plus_one_sq = Ciphertext()
  @testset "add_plain!" begin
    @test_nowarn add_plain!(x_plus_one_sq, x_encrypted, plain_one, evaluator)
  end

  @testset "square_inplace! and length/noise budget" begin
    @test_nowarn square_inplace!(x_plus_one_sq, evaluator)
    @test length(x_plus_one_sq) == 3
    @test invariant_noise_budget(x_plus_one_sq, decryptor) == 33
  end

  @testset "decrypt! and check ((x+1)^2 = 49 = 0x31)" begin
    @test_nowarn decrypt!(decrypted_result, x_plus_one_sq, decryptor)
    @test to_string(decrypted_result) == "31"
  end

  encrypted_result = Ciphertext()
  plain_four = Plaintext("4")
  @testset "compute encrypted_result (4(x^2+1)(x+1)^2)" begin
    @test_nowarn multiply_plain_inplace!(x_sq_plus_one, plain_four, evaluator)
    @test_nowarn multiply!(encrypted_result, x_sq_plus_one, x_plus_one_sq, evaluator)
    @test length(encrypted_result) == 5
    # Next test is fuzzy since the actual noise budget might vary
    @test invariant_noise_budget(encrypted_result, decryptor) in (3, 4, 5)
  end

  x_squared = Ciphertext()
  @testset "compute and relinearize x_squared (x^2)" begin
    @test_nowarn square!(x_squared, x_encrypted, evaluator)
    @test length(x_squared) == 3
    @test_nowarn relinearize_inplace!(x_squared, relin_keys_, evaluator)
    @test length(x_squared) == 2
  end

  @testset "compute x_sq_plus_one (x^2+1) and decrypt" begin
    @test_nowarn add_plain!(x_sq_plus_one, x_squared, plain_one, evaluator)
    @test invariant_noise_budget(x_sq_plus_one, decryptor) == 33
    @test_nowarn decrypt!(decrypted_result, x_sq_plus_one, decryptor)
    @test to_string(decrypted_result) == "25"
  end

  x_plus_one = Ciphertext()
  @testset "compute x_plus_one (x+1)" begin
    @test_nowarn add_plain!(x_plus_one, x_encrypted, plain_one, evaluator)
  end

  @testset "compute and relinearize x_plus_one_sq ((x+1)^2)" begin
    @test_nowarn square!(x_plus_one_sq, x_plus_one, evaluator)
    @test length(x_plus_one_sq) == 3
    @test_nowarn relinearize_inplace!(x_plus_one_sq, relin_keys_, evaluator)
    @test invariant_noise_budget(x_plus_one_sq, decryptor) == 33
  end

  @testset "decrypt (x+1)^2 and check" begin
    @test_nowarn decrypt!(decrypted_result, x_plus_one_sq, decryptor)
    @test to_string(decrypted_result) == "31"
  end

  @testset "compute and relinearize encrypted_result (4(x^2+1)(x+1)^2)" begin
    @test_nowarn multiply_plain_inplace!(x_sq_plus_one, plain_four, evaluator)
    @test_nowarn multiply!(encrypted_result, x_sq_plus_one, x_plus_one_sq, evaluator)
    @test length(encrypted_result) == 3
    @test_nowarn relinearize_inplace!(encrypted_result, relin_keys_, evaluator)
    @test length(encrypted_result) == 2
    # Next test is fuzzy since the actual noise budget might vary
    @test invariant_noise_budget(encrypted_result, decryptor) in (10, 11, 12)
  end

  @testset "decrypt encrypted_result (4(x^2+1)(x+1)^2) and check" begin
    @test_nowarn decrypt!(decrypted_result, encrypted_result, decryptor)
    @test to_string(decrypted_result) == "54"
  end

  @testset "invalid parameters" begin
    @test_nowarn set_poly_modulus_degree!(enc_parms, 2048)
    context = SEALContext(enc_parms)
    @test parameter_error_message(context) == "parameters are not compliant with HomomorphicEncryption.org security standard"
  end
end


text:
"""
    @exception_without_fields(
        exception_name::Symbol,
        docstring::Union{Expr, String},
        error_message_bits::Tuple{Vararg{Union{Expr, String}}},
    ) -> Expr

Create an exception with no fields.

# Arguments
- `exception_name::Symbol`: name of the exception
- `docstring::Union{Expr, String}`: documentation string
- `error_message_bits::Tuple{Vararg{Union{Expr, String}}}`: strings and expressions which
  will be interpolated in the `showerror` output

# Returns
- `Expr`: an exception definition (struct + `showerror` overload)
"""
@exception exception_without_fields context = begin
    docstring = Meta.quot(
        """
            $(exception_name) <: Exception

        $(docstring)

        See also: [`@exception_without_fields`](@ref)
        """
    )
end


text:
"""
    TupleDummy(arrays::Tuple)
Dummy structure for broadcast with multiple outputs.
A simplified SoA implementation.
**This is not part of the interface. Not exported.**
"""
struct TupleDummy{T,N,L,As,AXs} <: AbstractArray{T,N}
    arrays::As
    ax::AXs
    TupleDummy{T,N,L}(arrays::As, ax::AXs) where {T,N,L,As,AXs} =
        new{T,N,L,As,AXs}(arrays, ax)
end

function TupleDummy(arrays::Tuple{AbstractArray,AbstractArray,Vararg{AbstractArray}}) ## at least 2 outputs
    ax = getsame(axes, arrays...)
    Style = typeof(IndexStyle(arrays...))
    ElType = Tuple{eltype.(arrays)...}
    TupleDummy{ElType,length(ax),Style}(arrays, ax)
end
parent(td::TupleDummy) = td.arrays
size(td::TupleDummy, args...) = size(td.arrays[1], args...)
axes(td::TupleDummy) = td.ax
IndexStyle(T::Type{<:TupleDummy}) = T.parameters[3]()
@inline setindex!(td::TupleDummy{T,N,IndexLinear}, value::Tuple, ix::Int) where {T,N} =
    map((a, v) -> unsafe_setindex!(a, v, ix) , td.arrays, value)
@inline setindex!(td::TupleDummy{T,N,IndexCartesian}, value::Tuple, ixs::Vararg{Int,N}) where {T,N} =
    map((a, v) -> unsafe_setindex!(a, v, ixs...) , td.arrays, value)

@inline unalias(dest::TupleDummy, A::AbstractRange) = A
@inline unalias(dest::TupleDummy, A::AbstractArray) =
    mapreduce(x -> mightalias(x, A), |, dest.arrays) ? unaliascopy(A) : A
@inline broadcast_unalias(dest::TupleDummy, src::AbstractArray) =
    mapreduce(x -> x === src, |, dest.arrays) ? src : unalias(dest, src)

## toa_similar
_similar(bc, T) = similar(bc, T)
_similar(bc::Broadcasted{<:DefaultArrayStyle}, ::Type{Bool}) =
    similar(Array{Bool}, axes(bc))
function toa_similar(bc::Broadcasted)
    ElType = combine_eltypes(bc.f, bc.args)
    ElType <: Tuple{Any,Vararg{Any}} && Base.isconcretetype(ElType) ||
        throw(ArgumentError("$ElType is not a legal return type for @tab!"))
    dest = map(T -> _similar(bc, T), tuple(ElType.parameters...))
    TupleDummy{ElType,ndims(bc),IndexLinear}(dest, axes(bc))
end

## materialize
@inline tab_materialize(x) = x
@inline tab_materialize(bc::Broadcasted) = instantiate(bc) |> tab_copy

@inline materialize!(dest::NTuple{N,AbstractArray}, bc::Broadcasted) where {N} =
    (materialize!(TupleDummy(dest), bc); dest)

## copyto
@inline tab_copy(bc::FilledBC) = copy(bc)
@inline tab_copy(bc::Broadcasted) = copyto!(toa_similar(bc), bc) |> parent
@inline _maketuple(x, y) = tuple(x..., y)
@inline tab_copy(bc::Broadcasted{Style{Tuple}}) =
    foldl((x, y::Tuple{Any,Vararg{Any}}) -> _maketuple.(x, y), copy(bc))


text:
###############################################################################
#                       G MATRIX FUNCTIONS
################################################################################
function gmat_switch!(G, θ, covstr, r)
    if covstr.random[r].covtype.s == :SI
        gmat_si!(G, θ[covstr.tr[r]]) # i > r
    elseif covstr.random[r].covtype.s == :DIAG
        gmat_diag!(G, θ[covstr.tr[r]])
    elseif covstr.random[r].covtype.s == :AR
        gmat_ar!(G, θ[covstr.tr[r]])
    elseif covstr.random[r].covtype.s == :ARH
        gmat_arh!(G, θ[covstr.tr[r]])
    elseif covstr.random[r].covtype.s == :CSH
        gmat_csh!(G, θ[covstr.tr[r]])
    elseif covstr.random[r].covtype.s == :CS
        gmat_cs!(G, θ[covstr.tr[r]])
    elseif covstr.random[r].covtype.s == :ARMA
        gmat_arma!(G, θ[covstr.tr[r]])
    elseif covstr.random[r].covtype.s == :TOEP
        gmat_toep!(G, θ[covstr.tr[r]])
    elseif covstr.random[r].covtype.s == :TOEPP
        gmat_toepp!(G, θ[covstr.tr[r]],  covstr.random[r].covtype.p)
    elseif covstr.random[r].covtype.s == :TOEPH
        gmat_toeph!(G, θ[covstr.tr[r]])
    elseif covstr.random[r].covtype.s == :TOEPHP
        gmat_toephp!(G, θ[covstr.tr[r]],  covstr.random[r].covtype.p)
    elseif covstr.random[r].covtype.s == :FUNC
         covstr.random[r].covtype.f.xmat!(G, θ[covstr.tr[r]], covstr.random[r].covtype.p)
    end
    G
end
################################################################################
@inline function zgz_base_inc!(mx::AbstractArray{T}, θ::AbstractArray{T}, covstr, block, sblock) where T
    if covstr.random[1].covtype.s != :ZERO
        #length of random to covstr
        for r = 1:covstr.rn
            G = fill!(Symmetric(Matrix{T}(undef, covstr.q[r], covstr.q[r])), zero(T))
            gmat_switch!(G.data, θ, covstr, r)
            zblock    = view(covstr.z, block, covstr.zrndur[r])
            for i = 1:length(sblock[r])
                @inbounds mulαβαtinc!(view(mx, sblock[r][i], sblock[r][i]), view(zblock, sblock[r][i], :), G)
            end
        end
    end
    #    fill!(mx, zero(T))
    #end
    mx
end
################################################################################
function gmat_si!(mx, θ)
    val = θ[1] ^ 2
    @inbounds @simd for i = 1:size(mx, 1)
        mx[i, i] = val
    end
    nothing
end
function gmat_diag!(mx, θ)
    @inbounds @simd for i = 1:size(mx, 1)
        mx[i, i] = θ[i] ^ 2
    end
    nothing
end
#function gmat_vc!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T
#    nothing
#end
function gmat_ar!(mx, θ)
    de  = θ[1] ^ 2
    s   = size(mx, 1)
    @inbounds @simd for i = 1:s
        mx[i, i] = de
    end
    if s > 1
        for m = 1:s - 1
            @inbounds @simd for n = m + 1:s
                mx[m, n] = de * θ[2] ^ (n - m)
            end
        end
    end
    nothing
end
function gmat_arh!(mx, θ)
    s = size(mx, 1)
    @inbounds @simd for m = 1:s
        mx[m, m] = θ[m]
    end
    if s > 1
        for m = 1:s - 1
            @inbounds @simd for n = m + 1:s
                mx[m, n] = mx[m, m] * mx[n, n] * θ[end] ^ (n - m)
            end
        end
    end
    @inbounds @simd for m = 1:s
        mx[m, m] = mx[m, m] * mx[m, m]
    end
    nothing
end
function gmat_cs!(mx, θ)
    s = size(mx, 1)
    mx .= θ[1]^2
    if s > 1
        for m = 1:s - 1
            @inbounds @simd for n = m + 1:s
                mx[m, n] = mx[m, m] * θ[2]
            end
        end
    end
    nothing
end
function gmat_csh!(mx, θ)
    s = size(mx, 1)
    @inbounds @simd for m = 1:s
        mx[m, m] = θ[m]
    end
    if s > 1
        for m = 1:s - 1
            @inbounds @simd for n = m + 1:s
                mx[m, n] = mx[m, m] * mx[n, n] * θ[end]
            end
        end
    end
    @inbounds @simd for m = 1:s
        mx[m, m] = mx[m, m] * mx[m, m]
    end
    nothing
end
################################################################################
function gmat_arma!(mx, θ)
    de  = θ[1] ^ 2
    s   = size(mx, 1)
    @inbounds @simd for i = 1:s
        mx[i, i] = de
    end
    if s > 1
        for m = 1:s - 1
            @inbounds @simd for n = m + 1:s
                mx[m, n] = de * θ[2] * θ[3] ^ (n - m - 1)
            end
        end
    end
    nothing
end
function gmat_toep!(mx, θ)
    de  = θ[1] ^ 2    #diagonal element
    s   = size(mx, 1) #size
    @inbounds @simd for i = 1:s
        mx[i, i] = de
    end
    if s > 1
        for m = 1:s - 1
            @inbounds @simd for n = m + 1:s
                mx[m, n] = de * θ[n-m+1]
            end
        end
    end
    nothing
end
function gmat_toepp!(mx, θ, p::Int)
    de  = θ[1] ^ 2    #diagonal element
    s   = size(mx, 1) #size
    @inbounds @simd for i = 1:s
        mx[i, i] = de
    end
    if s > 1 && p > 1
        for m = 1:s - 1
            @inbounds @simd for n = m + 1:(m + p - 1 > s ? s : m + p - 1)
                mx[m, n] = de * θ[n - m + 1]
            end
        end
    end
    nothing
end
function gmat_toeph!(mx, θ)
    s = size(mx, 2)
    @inbounds @simd for m = 1:s
        mx[m, m] = θ[m]
    end
    if s > 1
        for m = 1:s - 1
            @inbounds @simd for n = m + 1:s
                mx[m, n] = mx[m, m] * mx[n, n] * θ[n-m+s]
            end
        end
    end
    @inbounds @simd for m = 1:s
        mx[m, m] = mx[m, m] * mx[m, m]
    end
    nothing
end
function gmat_toephp!(mx, θ, p::Int)
    s = size(mx, 2)
    @inbounds @simd for m = 1:s
        mx[m, m] = θ[m]
    end
    if s > 1 && p > 1
        for m = 1:s - 1
            for n = m + 1:(m + p - 1 > s ? s : m + p - 1)
                @inbounds mx[m, n] = mx[m, m] * mx[n, n] * θ[n - m + s]
            end
        end
    end
    @inbounds @simd for m = 1:s
        mx[m, m] = mx[m, m] * mx[m, m]
    end
    nothing
end


text:
function solve()
    n = parse(Int, readline())
    ret = 0
    for x=1:n
        base10 = string(x, base=10)
        base8 = string(x, base=8)
        if '7' ∉ base10*base8
            ret += 1
        end
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Pkg
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
CI && Pkg.activate(@__DIR__)
CI && Pkg.instantiate()
using Documenter, Neighborhood
using DocumenterTools: Themes

for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "style.scss"), String)
    theme = read(joinpath(@__DIR__, "$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "$(w).scss"), header*"\n"*theme)
end
Themes.compile(joinpath(@__DIR__, "light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

# %% actually make the docs
makedocs(
    modules=[Neighborhood],
    sitename= "Neighborhood.jl",
    authors = "George Datseris.",
    format = Documenter.HTML(
        prettyurls = CI,
        assets = [
            "assets/logo.ico",
            asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css),
            ],
        collapselevel = 2,
        ),
    doctest=false,
    pages = [
        "Public API" => "index.md",
        "Dev Docs" => "dev.md",
        ]
)

if CI
    deploydocs(
        repo = "github.com/JuliaNeighbors/Neighborhood.jl.git",
        target = "build",
        push_preview = true
    )
end


text:
using Nemo
using Distributions
using LinearAlgebra
using PolynomialRoots

struct CollocationFunction
    a::Array{Float64}
    b::Array{Float64}
    c::Array{Float64}
    x::Array{Float64}
    leftSlope::Float64
    rightSlope::Float64
end

function evaluate(self::CollocationFunction, z::Float64)
    if z <= self.x[1]
        return self.leftSlope*(z-self.x[1]) + self.a[1]
    elseif z >= self.x[end]
        return self.rightSlope*(z-self.x[end])+self.a[end]
    end
    i = searchsortedfirst(self.x,z)  # x[i-1]<z<=x[i]
    if i > 1
        i -= 1
    end
    h = z-self.x[i]
    return self.a[i] + h*(self.b[i]+h*self.c[i])
end

function solve(self::CollocationFunction, strike::Float64)
    if strike <= self.a[1]
        sn = self.leftSlope
        return (strike-self.a[1])/sn + self.x[1]
    elseif strike >= self.a[end]
        sn = self.rightSlope
        return (strike-self.a[end])/sn + self.x[end]
    end
    i = searchsortedfirst(self.a,strike)  # a[i-1]<strike<=a[i]
    if abs(self.a[i]-strike)< 1e-10
        return self.x[i]
    elseif abs(self.a[i-1]-strike)< 1e-10
        return self.x[i-1]
    elseif i == 1
        i += 1
    end
    x0 = self.x[i-1]
    c = self.c[i-1]
    b = self.b[i-1]
    a = self.a[i-1]
    d = 0
    cc = a + x0*(-b+x0*(c-d*x0)) - strike
    bb = b + x0*(-2*c+x0*3*d)
    aa = -3*d*x0 + c

    allck = roots([cc,bb,aa])
    for ck in allck
        if abs(imag(ck)) < 1e-10 && real(ck) >= self.x[i-1]-1e-10 && real(ck) <= self.x[i]+1e-10
            return real(ck)
        end
    end
    println("ERRROR !!! no roots found in range", allck," ", strike, " ",aa, bb, cc," ", i," ",self.x[i-1]," ",self.x[i])
    return self.x[1]
end

function phi(cx,cv, dx, dv)
    sqx = (dx ./ cx).^2
    sqv = (dv ./ cv).^2
    phi = sqrt.(1 .+ sqx+sqv)
    dphidx = (dx ./ (cx .* cx)) ./ phi
    phi3 = phi .^3
    d2phidx2 = (1.0 ./ (cx .* cx)) ./ phi - (sqx ./ (cx .* cx)) ./ phi3
    dphidv = (dv ./ (cv .* cv)) ./ phi
    d2phidv2 = (1.0 ./ (cv .* cv)) ./ phi - (sqv ./ (cv .* cv)) ./ phi3
    d2phidxdv = -((dx ./ (cx .* cx)) .* (dv ./ (cv .* cv))) ./ phi3
    return phi, dphidx, d2phidx2, dphidv, d2phidv2, d2phidxdv
end

function makeRBFMatrix(PX,PV,cx,cv)
    len = length(PX)
    rx = zeros(typeof(PX[1,1]),len,len)
    rv = zeros(typeof(PX[1,1]),len,len)
    for i = 1:len
        rx[:,i] = PX .- PX[i]
        rv[:,i] = PV .- PV[i]
    end
    I,dIdX, d2IdX2, dIdV, d2IdV2, d2IdXdV = phi(cx,cv,rx,rv)
    return I,dIdX, d2IdX2, dIdV, d2IdV2, d2IdXdV
end

function price(isCall, xFactor, vFactor, spotArray, priceArray, v0, kappa, theta, sigma, rho, r, q, T, cFunc, K, B, N, M, L)
    epsilon = 1e-4
    boundaryCondition = "Dirichlet" #Dirichlet, ZeroGamma, Gamma
    dChi = 4*kappa*theta/(sigma*sigma)
    chiN = 4*kappa*exp(-kappa*T)/(sigma*sigma*(1-exp(-kappa*T)))
    ncx2 = NoncentralChisq(dChi,v0*chiN)
    vmax = quantile(ncx2,1-epsilon)*exp(-kappa*T)/chiN
    vmin = quantile(ncx2,epsilon)*exp(-kappa*T)/chiN
    vmin = max(vmin,1e-4)
    V = collect(range(vmin,stop=vmax,length=L))
    Xspan = 4*sqrt(theta*T)
    logK = log(K)
    Kinv = solve(cFunc,K)
    logKinv = log(Kinv)
    Xmin = logK - Xspan + (r-q)*T - 0.5*v0*T
    Xmax = logK + Xspan + (r-q)*T - 0.5*v0*T
    X = collect(range(Xmin,stop=Xmax,length=M))
    hm = X[2]-X[1]
    S = exp.(X)
    cx = xFactor*(Xmax-Xmin)/(M-1)
    cv = vFactor*(vmax-vmin)/(L-1)
    Sc = zeros(M)
    for i = 1:M
        Sc[i] = evaluate(cFunc,S[i])
    end
    if isCall
        F0 = max.(Sc .- K, 0.0)
    else
        F0 = max.(K .- Sc, 0.0)
    end
    F0smooth = F0
    F = zeros(L*M)
    for j = 1:L
        F[1+(j-1)*M:j*M] = F0smooth
    end
    PX = zeros(Float64,L*M)
    PV = zeros(typeof(PX[1,1]),L*M)
    for j = 1:L
        PX[1+(j-1)*M:j*M]=X
        PV[1+(j-1)*M:j*M] .= V[j]
    end
    Iphi,dIdX, d2IdX2, dIdV, d2IdV2, d2IdXdV = makeRBFMatrix(PX,PV,cx,cv)
    dt = -T/N
    PV2D = zeros(typeof(PX[1,1]),L*M,L*M)
    for j = 1:L*M
        PV2D[:,j] = PV
    end
    dtPV = (0.5*dt) .* PV2D
    A1 = dtPV .* d2IdX2 + (dt*(r-q) .- dtPV) .* dIdX - (r*dt) .* Iphi
    A2 = (sigma^2 .* dtPV) .* d2IdV2 + (dt*kappa .* (theta .- PV2D)) .* dIdV + (dt*rho*sigma .* PV2D) .* d2IdXdV
    A = A1+A2
    etime = @elapsed begin
        #dirichlet
        if boundaryCondition == "Dirichlet"
            for j = 1:L
                A[1+(j-1)*M,:] .= 0
                A[M+(j-1)*M,:] .= 0
            end
        elseif  boundaryCondition == "ZeroGamma"
            #second = 0
            for j = 1:L
                A[1+(j-1)*M,:] .= d2IdX2[1+(j-1)*M,:] .- dIdX[1+(j-1)*M,:]
                A[M+(j-1)*M,:] .= d2IdX2[M+(j-1)*M,:] .- dIdX[M+(j-1)*M,:]
            end
        elseif  boundaryCondition == "Gamma"
            #second = 0
            for j = 1:L
                A[1+(j-1)*M,:] .= (dt*(r-q) ) .* dIdX[1+(j-1)*M,:] - (r*dt) .* Iphi[1+(j-1)*M,:]
                A[M+(j-1)*M,:] .=(dt*(r-q) ) .* dIdX[M+(j-1)*M,:] - (r*dt) .* Iphi[M+(j-1)*M,:]
            end
        end

        reg = 0e-15
        # IphiR = matrix(R, [R(Iphi[i,j]) for i in 1:L*M, j in 1:L*M])
        # AR = matrix(R, [R(A[i,j]) for i in 1:L*M, j in 1:L*M])
        # Iphi = convert(Array{BigFloat},Iphi)
        # A = convert(Array{BigFloat},A)
        lu0 = factorize(Iphi+reg*LinearAlgebra.I)
        a = 1.0 - sqrt(2.0)/2.0
        Li = Iphi + a .* A
        if boundaryCondition == "ZeroGamma"
            for j = 1:L
                Li[1+(j-1)*M,:] .= A[1+(j-1)*M,:]
                Li[M+(j-1)*M,:] .= A[M+(j-1)*M,:]
            end
        end
        lui = factorize(Li)
        # println("Condition Iphi ",cond(Iphi), " A ",cond(Li))
        # d0 = diag(lu0.U)
        # di = diag(lui.U)
        # println("sq ratio ", maximum(d0 .* d0)/minimum(d0 .* d0)," A ", maximum(di .* di)/minimum(di .* di))
        for i = 1:N
            ti = (T*i)/N
            for j = 1:L
                if boundaryCondition == "Dirichlet"
                    if isCall
                        F[1+(j-1)*M] = 0
                        F[M+(j-1)*M] =  max(Sc[M]*exp(-q*ti)-K*exp(-r*ti),0.0)
                    else
                        F[1+(j-1)*M] = max(K*exp(-r*ti)-Sc[1]*exp(-q*ti),0.0)
                        F[M+(j-1)*M] = 0.0
                    end
                elseif  boundaryCondition == "ZeroGamma"
                    F[1+(j-1)*M] = 0
                    F[M+(j-1)*M] = 0
                end
            end
            coeff = lui \ F
            F1 = Iphi*coeff
            for j = 1:L
                if boundaryCondition == "Dirichlet"

                    if isCall
                        F1[1+(j-1)*M] = 0
                        F1[M+(j-1)*M] =  max(Sc[M]*exp(-q*ti)-K*exp(-r*ti),0.0)
                    else
                        F1[1+(j-1)*M] = max(K*exp(-r*ti)-Sc[1]*exp(-q*ti),0.0)
                        F1[M+(j-1)*M] = 0.0
                    end
                elseif  boundaryCondition == "ZeroGamma"
                    F1[1+(j-1)*M] = 0
                    F1[M+(j-1)*M] = 0
                end
            end
            coeff = lui \ F1
            F2 = Iphi*coeff
            for j = 1:L
                if boundaryCondition == "Dirichlet"
                    if isCall
                        F2[1+(j-1)*M] = 0
                        F2[M+(j-1)*M] =  max(Sc[M]*exp(-q*ti)-K*exp(-r*ti),0.0)
                    else
                        F2[1+(j-1)*M] = max(K*exp(-r*ti)-Sc[1]*exp(-q*ti),0.0)
                        F2[M+(j-1)*M] = 0.0
                    end
                elseif  boundaryCondition == "ZeroGamma"
                    F2[1+(j-1)*M] = 0
                    F2[M+(j-1)*M] = 0
                end
            end
            F = (1+sqrt(2.0)) .* F2 - sqrt(2.0) .* F1
        end
    end  #elapsed
    coeff = lu0 \ F

    for i = 1:length(spotArray)
        spot = spotArray[i]
        refPrice = priceArray[i]
        logspot = log(spot)
        phispotv0,dphi,d2phi,dphidv,d2phidv2, dphidxdv = phi(cx,cv,logspot  .- PX,v0 .- PV)
        price = dot(phispotv0,  coeff)
        delta = dot(dphi, coeff)/spot
        gamma = dot(d2phi, coeff)/(spot^2)-delta/spot
        error = price -refPrice
        println(K," ",cx," ",cv," ",xFactor," ",vFactor," Float64 ",N, " ",M," ", L," ", Float64(price)," ", Float64(delta)," ",Float64(gamma)," ",Float64(error), " ",etime)
    end

    for spot in S
        logspot = log(spot)
        phispotv0,dphi,d2phi,dphidv,d2phidv2, dphidxdv = phi(cx,cv,logspot  .- PX,v0 .- PV)
        price = dot(phispotv0,  coeff)
        delta = dot(dphi, coeff)/spot
        gamma = dot(d2phi, coeff)/(spot^2)-delta/spot
        println(spot," LS ",N, " ",M," ", L," ", price," ", delta," ",gamma)
    end
end

function approx_solve!(z::arb_mat, x::arb_mat, y::arb_mat)
    r = ccall((:arb_mat_approx_solve, :libarb), Cint,
    (Ref{arb_mat}, Ref{arb_mat}, Ref{arb_mat}, Int),
    z, x, y, prec(base_ring(x)))
    r == 0 && error("Matrix cannot be inverted numerically")
    nothing
end

function approx_inv(x::arb_mat)
    y = one(parent(x))
    approx_solve!(y, x, y)
    return y
end
function priceInv(isCall, spotArray, priceArray, v0, kappa, theta, sigma, rho, r, q, T, cFunc, K, B, N, M, L)
    epsilon = 1e-3
    dChi = 4*kappa*theta/(sigma*sigma)
    chiN = 4*kappa*exp(-kappa*T)/(sigma*sigma*(1-exp(-kappa*T)))
    ncx2 = NoncentralChisq(dChi,v0*chiN)
    vmax = quantile(ncx2,1-epsilon)*exp(-kappa*T)/chiN
    vmin = quantile(ncx2,epsilon)*exp(-kappa*T)/chiN
    vmin = max(vmin,1e-4)
    V = collect(range(vmin,stop=vmax,length=L))
    Xspan = 8*sqrt(theta*T)
    logK = log(K)
    Kinv = solve(cFunc,K)
    logKinv = log(Kinv)
    Xmin = logK*0 - Xspan + (r-q)*T - 0.5*v0*T
    Xmax = logK*0 + Xspan + (r-q)*T - 0.5*v0*T
    X = collect(range(Xmin,stop=Xmax,length=M))
    hm = X[2]-X[1]
    S = exp.(X)
    cx = 5*(Xmax-Xmin)/(M-1)
    cv = 10*(vmax-vmin)/(L-1)
    Sc = zeros(M)
    for i = 1:M
        Sc[i] = evaluate(cFunc,S[i])
    end
    if isCall
        F0 = max.(Sc .- K,0)
    else
        F0 = max.(K .- Sc,0)
    end
    F0smooth = F0
    F = zeros(L*M)
    for j = 1:L
        F[1+(j-1)*M:j*M] = F0smooth
    end
    PX = zeros(Float64,L*M)
    PV = zeros(typeof(PX[1,1]),L*M)
    for j = 1:L
        PX[1+(j-1)*M:j*M]=X
        PV[1+(j-1)*M:j*M] .= V[j]
    end
    Iphi,dIdX, d2IdX2, dIdV, d2IdV2, d2IdXdV = makeRBFMatrix(PX,PV,cx,cv)
    dt = -T/N
    PV2D = zeros(typeof(PX[1,1]),L*M,L*M)
    for j = 1:L*M
        PV2D[:,j] = PV
    end
    dtPV = (0.5*dt) .* PV2D
    A1 = dtPV .* d2IdX2 + (dt*(r-q) .- dtPV) .* dIdX - (r*dt) .* Iphi
    A2 = (sigma^2 .* dtPV) .* d2IdV2 + (dt*kappa .* (theta .- PV2D)) .* dIdV + (dt*rho*sigma .* PV2D) .* d2IdXdV
    A = A1+A2
    etime = @elapsed begin
        for j = 1:L
            A[1+(j-1)*M,:] .= 0
            A[M+(j-1)*M,:] .= 0
        end
        a = 1.0 - sqrt(2.0)/2.0
        reg = 0e-15
        # Iphi = convert(Array{BigFloat},Iphi)
        # A = convert(Array{BigFloat},A)
        # lu0 = factorize(Iphi+reg*LinearAlgebra.I)
        # IphiInv = inv(Iphi)
        # # Li = Iphi + a .* A
        # # lui = factorize(Li)
        # LiInv = inv(LinearAlgebra.I + a .* (A*IphiInv))

        IphiR = matrix(R, [R(Iphi[i,j]) for i in 1:L*M, j in 1:L*M])
        AR = matrix(R, [R(A[i,j]) for i in 1:L*M, j in 1:L*M])
        FR = zero_matrix(R, L*M, 1)
        for i = 1:L*M
            FR[i,1] = R(F[i])
        end

        IphiInv = approx_inv(IphiR)
        B = AR # Ip * inv(Ip + a*A) = inv (I + a*A*inv(Ip))
        B *= R(a)
        # for i=1:L*M
        #     B[i,i] += 1
        # end
        B += IphiR
        LiInv = approx_inv(B)
        LiInv = IphiR*LiInv
        for i = 1:N
            ti = (T*i)/N
            for j = 1:L
                if isCall
                    FR[1+(j-1)*M,1] = 0
                    FR[M+(j-1)*M,1] =  max(Sc[M]*exp(-q*ti)-K*exp(-r*ti),0.0)
                else
                    FR[1+(j-1)*M,1] = max(K*exp(-r*ti)-Sc[1]*exp(-q*ti),0.0)
                    FR[M+(j-1)*M,1] = 0.0
                end
            end
            #coeff = LiInv * F
            #F1 = Iphi*coeff
            F1 = LiInv * FR
            for j = 1:L
                if isCall
                    F1[1+(j-1)*M,1] = 0
                    F1[M+(j-1)*M,1] =  max(Sc[M]*exp(-q*ti)-K*exp(-r*ti),0.0)
                else
                    F1[1+(j-1)*M,1] = max(K*exp(-r*ti)-Sc[1]*exp(-q*ti),0.0)
                    F1[M+(j-1)*M,1] = 0.0
                end
            end
            # coeff = LiInv * F1
            # F2 = Iphi*coeff
            F2 = LiInv*F1
            for j = 1:L
                if isCall
                    F2[1+(j-1)*M,1] = 0
                    F2[M+(j-1)*M,1] =  max(Sc[M]*exp(-q*ti)-K*exp(-r*ti),0.0)
                else
                    F2[1+(j-1)*M,1] = max(K*exp(-r*ti)-Sc[1]*exp(-q*ti),0.0)
                    F2[M+(j-1)*M,1] = 0.0
                end
            end
            for index=1:L*M
                FR[index,1] = (1+sqrt(2.0)) * F2[index,1] - sqrt(2.0)*F1[index,1]
            end
        end
    end  #elapsed
    coeffR = IphiInv * FR
    coeff=[Float64(coeffR[i,1]) for i in 1:L*M]
    println("coeff",coeff)
    for i = 1:length(spotArray)
        spot = spotArray[i]
        refPrice = priceArray[i]
        logspot = log(spot)
        phispotv0,dphi,d2phi,dphidv,d2phidv2, dphidxdv = phi(cx,cv,logspot  .- PX,v0 .- PV)
        price = dot(phispotv0,  coeff)
        delta = dot(dphi, coeff)/spot
        gamma = dot(phispotv0, coeff)/(spot^2)
        error = price -refPrice
        println(spot," LS ",N, " ",M," ", L," ", price," ", delta," ",gamma," ",error, " in ",etime)
    end
end

A=[0.6266758553145932, 0.8838690008217314 ,0.9511741483703275, 0.9972169412308787 ,1.045230848712316, 1.0932361943842062, 1.1786839882076958, 1.2767419415280061]
B=[0.8329310535215612, 0.5486175716699259, 1.0783076034285555 ,1.1476195823811128 ,1.173600641673776, 1.1472056638621118, 0.918270335988941]
C=[-0.38180731761048253, 3.2009663415588276, 0.8377175268235754, 0.31401193651971954 ,-0.31901463307065175, -1.3834775717464938, -1.9682171790586938]
X=[0.5171192610665245, 0.8894451290344221, 0.972184210805066, 1.013553751690388, 1.05492329257571, 1.0962928334610318, 1.179031915231676, 1.3445100787729636]
leftSlope=0.8329310535215612
rightSlope=0.2668764075068484
cFunc = CollocationFunction(A,B,C,X,leftSlope,rightSlope)
A = X
leftSlope = 1.0
rightSlope = 1.0
C = 0*C
B = [1.0 for i=1:length(B)]
kappa = 0.35
theta = 0.321
sigma = 1.388
rho = -0.63
r = 0.0
q = 0.0
v0=0.133
T=0.4986301369863014
M=30
L=30
N=32
spotArray = [1.0]
K=1.0
priceArray= [0.07278065]
K=0.7
priceArray = [0.00960629]
K=1.4
priceArray = [0.40015225]
K=0.8
priceArray = [0.01942598]
K=1.2
priceArray = [0.20760077]
isCall = false
xFactor=2.0
vFactor=2.0
xFactors = [1+Float64(10*i)/40 for i in 0:40]
xFactors = reverse(xFactors)
for xFactor in xFactors
    vFactor = xFactor
    price(isCall,xFactor,vFactor, spotArray, priceArray, v0, kappa, theta, sigma, rho, r, q, T, cFunc, K, B, N, M, L)
end

max_stars_repo_name: fabienlefloch/pythonjuliapde

max_issues_repo_name: fabienlefloch/pythonjuliapde

max_forks_repo_name: fabienlefloch/pythonjuliapde


text:
# Function to read priors from a csv file. returns a Dict
#
# input :
# -------
# * path to csv file
#
# the csv file should contain at least two columns with:
# * "name" (name of the parameter)
# * "value" (value of the parameter)
# * "upper_bound" (upper bound for the support of parameter "name")
# * "lower_bound" (upper bound for the support of parameter "name")
#
# output :
# -------
# * a dictionary with "name" => [value, lower_bound, upper_bound]
#
# Remark:
#--------
# *
#----------------------------------------------------------------------------
"""
  read_priors(pathToCsv::String)

Function to load priors stored in a csv file. Returns a dictionary. The csv
file should have the following columns "name", "value", "upper_bound", "lower_bound"
"""
 function read_priors(pathToCsv::String)

   # initialize and empty Dictionary:
   #---------------------------------
   dictionary = OrderedDict{String,Array{Float64,1}}()

   # Check the existence of the file
   #----------------------------------
   if isfile(pathToCsv) == false
     error(string(pathToCsv," not found."))
   end

   # If file exists
   #----------------------------------
   dataFrame = CSV.read(pathToCsv)

   # Check whether the csv file has the appropriate columns
   #-------------------------------------------------------
   ListToCheck = [:name, :value, :upper_bound, :lower_bound]

   for colName in ListToCheck
     try
       dataFrame[colName]
     catch
       error("column $(colName) does not exist.")
     end
   end

   # store the number of rows in the dataframe:
   #------------------------------------------
   number_rows = size(dataFrame, 1)
   info(string(number_rows, " prior value(s) found"))

   # Append the dictionary
   #----------------------
   for i=1:number_rows
      dictionary[dataFrame[i,:name]] = [dataFrame[i,:value],  dataFrame[i,:lower_bound], dataFrame[i,:upper_bound]]
   end

   return dictionary

 end

#----------------------------------------------------------------------------
# Function to read moments and weights from a csv file. returns a DataFrame
#
# input :
# -------
# * path to csv file
#
# the csv file should contain at least two columns with:
# * "name" (name of the parameter)
# * "value" (value of the parameter)
# * "weight" (weight associated to moment "name")
#
# output :
# -------
# a dictionary with "name" => [value, weights]
#
# Remark:
#--------
#
#----------------------------------------------------------------------------
"""
  read_empirical_moments(pathToCsv::String)

Function to load empirical moments stored in a csv file. Returns a dictionary.
The csv file should have the following columns "name", "value", "weight"
"""
function read_empirical_moments(pathToCsv::String)

  # initialize and empty Dictionary:
  #---------------------------------
  dictionary = OrderedDict{String,Array{Float64,1}}()

  # Check the existence of the file
  #----------------------------------
  if isfile(pathToCsv) == false
   error(string(pathToCsv," not found"))
  end

  # If file exists:
  #----------------
  dataFrame =  CSV.read(pathToCsv)

  #check whether the csv file has the columns "name" and "value"
  #-------------------------------------------------------------
  ListToCheck = [:name, :value, :weight]

  for colName in ListToCheck
   try
     dataFrame[colName]
   catch
     error("column $(colName) does not exist.")
   end
  end

  # store the number of rows in the dataframe:
  #------------------------------------------
  number_rows = size(dataFrame,1)
  info(string(number_rows, " moment(s) found."))

  # Append the dictionary
  #----------------------
  for i=1:number_rows
     dictionary[dataFrame[i,:name]] = [dataFrame[i,:value],  dataFrame[i,:weight]]
  end

  return dictionary

end

"""
  saveSMMOptim(sMMProblem::SMMProblem; verbose::Bool = false, saveName::String = "")

Function to save a SMMProblem to disk.
"""
function saveSMMOptim(sMMProblem::SMMProblem; verbose::Bool = false, saveName::String = "")

    bestValue = best_candidate(sMMProblem.bbResults)
    bestFitness = best_fitness(sMMProblem.bbResults)

    if verbose == true
        info("Best Value: ", bestValue)
        info("Best Fitness: ", bestFitness)
        info("Saving optimisation to disk")
    end

    # If no name is provided, generate a random name
    #-----------------------------------------------
    if isempty(saveName) == true
        saveName = string(rand(1:Int(1e8)))
    end

    tempfilename = saveName * ".jld2"
    JLD2.@save tempfilename sMMProblem

    if verbose == true
        info("Done.")
    end

end

"""
  loadSMMOptim(saveName::String; verbose::Bool = false)

Function to save a SMMProblem to disk.
"""
function loadSMMOptim(saveName::String; verbose::Bool = false)

    if verbose == true
        info("Loading optimisation from disk")
    end

    #-------------------------------
    # Load BBsetup and BBoptimize
    #--------------------------------
    if isempty(saveName) == true
        error("Please enter a non-empty saveName.")
    end

    tempfilename = saveName * ".jld2"
    JLD2.@load tempfilename sMMProblem

    if verbose == true
        info("Done.")
    end

    return sMMProblem

end


text:
export GramMatrix

import MultivariateMoments: trimat, SymMatrix, getmat
export gram_operate, getmat

"""
    struct GramMatrix{T, MT <: MP.AbstractMonomial, MVT <: AbstractVector{MT}} <: MP.APL{T}
        Q::SymMatrix{T}
        x::MVT
    end

Gram matrix ``x^\\top Q x`` where `Q` is a symmetric matrix indexed by the
vector of monomials `x`.
"""
struct GramMatrix{T, MT <: MP.AbstractMonomial, MVT <: AbstractVector{MT}} <: MP.APL{T} # should be MP.APL{eltype(T)} but it doesn't work
    Q::SymMatrix{T}
    x::MVT
end
# When taking the promotion of a GramMatrix of JuMP.Variable with a Polynomial JuMP.Variable, it should be a Polynomial of AffExpr
MP.coefficienttype(::Type{<:GramMatrix{T}}) where {T} = Base.promote_op(+, T, T)
function MP.constantmonomial(p::GramMatrix)
    if isempty(p.x)
        return MP.constantmonomial(MP.monomialtype(p))
    else
        return MP.constantmonomial(p.x[1])
    end
end
function MP.monomialtype(::Union{GramMatrix{T, MT},
                                 Type{<:GramMatrix{T, MT}}}) where {T, MT}
    return MT
end
MP.polynomialtype(::Type{GramMatrix{T, MT, MVT}}) where {T, MT, MVT} = MP.polynomialtype(MT, MP.coefficienttype(GramMatrix{T, MT, MVT}))
MP.polynomialtype(::Type{GramMatrix{T, MT, MVT}}, ::Type{S}) where {S, T, MT, MVT} = MP.polynomialtype(MT, S)
MP.variables(p::GramMatrix) = MP.variables(p.x)
MP.nvariables(p::GramMatrix) = MP.nvariables(p.x)

Base.zero(::Type{GramMatrix{T, MT, MVT}}) where {T, MT, MVT} = GramMatrix{T, MT, MP.monovectype(MT)}(SymMatrix{T}(T[], 0), MP.emptymonovec(MT))
Base.iszero(p::GramMatrix) = iszero(MP.polynomial(p))

Base.getindex(p::GramMatrix, I...) = getindex(p.Q, I...)
Base.copy(p::GramMatrix) = GramMatrix(copy(p.Q), copy(p.x))

MultivariateMoments.getmat(p::GramMatrix{T}) where {T} = p.Q

function GramMatrix{T}(f::Function, x::AbstractVector{MT}, σ) where {T, MT<:MP.AbstractMonomial}
    GramMatrix{T, MT, MP.monovectype(x)}(trimat(T, f, length(x), σ), x)
end
GramMatrix{T}(f::Function, x::AbstractVector, σ) where T = GramMatrix{T}(f, MP.monomial.(x), σ)
function GramMatrix{T}(f::Function, x::AbstractVector) where T
    σ, X = MP.sortmonovec(x)
    GramMatrix{T}(f, X, σ)
end
GramMatrix(f::Function, x) = GramMatrix{Base.promote_op(f, Int, Int)}(f, x)

function GramMatrix(Q::AbstractMatrix{T}, x, σ) where T
    GramMatrix{T}((i,j) -> Q[σ[i], σ[j]], x)
end
function GramMatrix(Q::AbstractMatrix, x)
    σ, X = MP.sortmonovec(x)
    GramMatrix(Q, X, σ)
end

#function Base.convert{T, PT <: AbstractPolynomial{T}}(::Type{PT}, p::GramMatrix)
#    # coefficienttype(p) may be different than T and MP.polynomial(p) may be different than PT (different module)
#    convert(PT, MP.polynomial(p))
#end
function MP.polynomial(p::GramMatrix)
    MP.polynomial(getmat(p), p.x)
end
function MP.polynomial(p::GramMatrix, ::Type{S}) where {S}
    MP.polynomial(getmat(p), p.x, S)
end

# The `i`th index of output is the index of occurence of `x[i]` in `y`,
# or `0` if it does not occur.
function multi_findsorted(x, y)
    I = zeros(Int, length(x))
    j = 1
    for i in eachindex(x)
        while j ≤ length(y) && x[i] < y[j]
            j += 1
        end
        if j ≤ length(y) && x[i] == y[j]
            I[i] = j
        end
    end
    return I
end

"""
    gram_operate(::typeof(+), p::GramMatrix{S}, q::GramMatrix{T})

Computes the Gram matrix equal to the sum between `p` and `q`. On the opposite,
`p + q` gives a polynomial equal to `p + q`. The polynomial `p + q` can also be
obtained by `polynomial(gram_sum(p, q))`.
"""
function gram_operate(::typeof(+), p::GramMatrix{S}, q::GramMatrix{T}) where {S, T}
    monos = MultivariatePolynomials.mergemonovec([p.x, q.x])
    U = typeof(zero(S) + zero(T))
    n = length(monos)
    Q = SymMatrix(zeros(U, div(n * (n + 1), 2)), n)
    Ip = multi_findsorted(monos, p.x)
    Iq = multi_findsorted(monos, q.x)
    for j in 1:length(monos)
        for i in 1:j
            if !iszero(Ip[j]) && !iszero(Ip[i])
                MultivariateMoments.symmetric_setindex!(
                    Q, Q[i, j] + p.Q[Ip[i], Ip[j]], i, j)
            end
            if !iszero(Iq[j]) && !iszero(Iq[i])
                MultivariateMoments.symmetric_setindex!(
                    Q, Q[i, j] + q.Q[Iq[i], Iq[j]], i, j)
            end
        end
    end
    return GramMatrix(Q, monos)
end

"""
    gram_operate(/, p::GramMatrix{S}, q::GramMatrix{T})

Computes the Gram matrix equal to the sum between `p` and `q`. On the opposite,
`p + q` gives a polynomial equal to `p + q`. The polynomial `p + q` can also be
obtained by `polynomial(gram_sum(p, q))`.
"""
function gram_operate(::typeof(/), q::GramMatrix, α)
    Q = SymMatrix(q.Q.Q / α, q.Q.n)
    return GramMatrix(Q, q.x)
end

Base.:(+)(x::MP.APL, y::GramMatrix) = x + MP.polynomial(y)
Base.:(+)(x::GramMatrix, y::MP.APL) = MP.polynomial(x) + y
Base.:(+)(x::GramMatrix, y::GramMatrix) = MP.polynomial(x) + MP.polynomial(y)
Base.:(-)(x::MP.APL, y::GramMatrix) = x - MP.polynomial(y)
Base.:(-)(x::GramMatrix, y::MP.APL) = MP.polynomial(x) - y
Base.:(-)(x::GramMatrix, y::GramMatrix) = MP.polynomial(x) - MP.polynomial(y)
Base.:(==)(p::MP.APL, q::GramMatrix) = p == MP.polynomial(q)
Base.:(==)(p::GramMatrix, q::GramMatrix) = iszero(p - q)
function Base.isapprox(p::GramMatrix, q::GramMatrix; kwargs...)
    p.x == q.x && isapprox(p.Q, q.Q; kwargs...)
end


text:
"""
    eom_henon4(x, p, n) -> Function

Equations of motion for a 4D Henon system consisting of two identical
Henon maps, `X` and `Y` with unidirectional forcing from `X` to`Y` [1].
The coupling constant `c` controls the strength of the forcing.

The implementation here has adjustable parameters; the default values are the
ones used in [1].

# References
1. Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D.,
Jajcay, N., & Paluš, M. (2018). Comparison of six methods for the detection of
causality in a bivariate time series. Physical Review E, 97(4), 042207.
"""
function eom_henon4(x, p, n)
    a₁, a₂, b₁, b₂, c, k = (p...,)
    x₁, x₂, y₁, y₂ = (x...,)

    dx₁ = a₁ - x₁^2 + k*x₂
    dx₂ = b₁*x₁
    dy₁ = a₂ - (c * x₁ * y₁ + (1 - c)*y₁^2) + k*y₂
    dy₂ = b₂*y₁
    return SVector{4}(dx₁, dx₂, dy₁, dy₂)
end

"""
    henon4(u₀, c) -> DiscreteDynamicalSystem

Initialize an instance of a 4D Henon map system consisting of two identical
Henon systems, ``x`` and ``y`` [1]. The subsystems are unidirectionally
coupled . The coupling constant `c` controls the strength of the forcing.
Synchronization occurs when the values of the coupling constant
``c > 0.7``.

The difference equations are:

```math
\\begin{aligned}
x_1(t+1) &= a_1 - x_1^2(t) + k x_2(t) \\
x_2(t+1) &= x_1(t) \\
y_1(t+1) &= a_1 - [c x_1(t) y_1(t) + (1-c) y_1^2(t)] + k y_2(t) \\
y_2(t+1) &= y_1(t)
\\end{aligned}
```

This system was investigated by Krakovská to study the performance of different
causality detection algorithms.

# References
Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D., Jajcay, N., & Paluš, M. (2018). Comparison of six methods for the detection of causality in a bivariate time series. Physical Review E, 97(4), 042207.
"""
function henon4(u₀, a₁, a₂, b₁, b₂, c, k)
    p = [a₁, a₂, b₁, b₂, c, k]
    logistic_system = DiscreteDynamicalSystem(eom_henon4, u₀, p)
    return logistic_system
end

henon4(;u₀ = rand(4),
        a₁ = 1.4, a₂ = 1.4, b₁ = 1.0, b₂ = 1.0, c = 2.0, k = 0.3) =
    henon4(u₀, a₁, a₂, b₁, b₂, c, k)


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

using Base: llvmcall

#function add1234(x::Tuple{Int32,Int32,Int32,Int32})
#    llvmcall("""%3 = add <4 x i32> %1, %0
#                ret <4 x i32> %3""",
#        Tuple{Int32,Int32,Int32,Int32},
#        Tuple{Tuple{Int32,Int32,Int32,Int32},
#        Tuple{Int32,Int32,Int32,Int32}},
#        (Int32(1),Int32(2),Int32(3),Int32(4)),
#        x)
#end
#
#function add1234(x::NTuple{4,Int64})
#    llvmcall("""%3 = add <4 x i64> %1, %0
#      ret <4 x i64> %3""",NTuple{4,Int64},
#      Tuple{NTuple{4,Int64},NTuple{4,Int64}},
#        (Int64(1),Int64(2),Int64(3),Int64(4)),
#      x)
#end
#
function add1234(x::Tuple{Int32,Int32,Int32,Int32})
    llvmcall("""%3 = extractvalue [4 x i32] %0, 0
      %4 = extractvalue [4 x i32] %0, 1
      %5 = extractvalue [4 x i32] %0, 2
      %6 = extractvalue [4 x i32] %0, 3
      %7 = extractvalue [4 x i32] %1, 0
      %8 = extractvalue [4 x i32] %1, 1
      %9 = extractvalue [4 x i32] %1, 2
      %10 = extractvalue [4 x i32] %1, 3
      %11 = add i32 %3, %7
      %12 = add i32 %4, %8
      %13 = add i32 %5, %9
      %14 = add i32 %6, %10
      %15 = insertvalue [4 x i32] undef, i32 %11, 0
      %16 = insertvalue [4 x i32] %15, i32 %12, 1
      %17 = insertvalue [4 x i32] %16, i32 %13, 2
      %18 = insertvalue [4 x i32] %17, i32 %14, 3
      ret [4 x i32] %18""",Tuple{Int32,Int32,Int32,Int32},
      Tuple{Tuple{Int32,Int32,Int32,Int32},Tuple{Int32,Int32,Int32,Int32}},
        (Int32(1),Int32(2),Int32(3),Int32(4)),
        x)
end

@test add1234(map(Int32,(2,3,4,5))) === map(Int32,(3,5,7,9))
#@test add1234(map(Int64,(2,3,4,5))) === map(Int64,(3,5,7,9))

# Test whether llvmcall escapes the function name correctly
baremodule PlusTest
    using Base: llvmcall
    using Test
    using Base

    function +(x::Int32, y::Int32)
        llvmcall("""%3 = add i32 %1, %0
                    ret i32 %3""",
            Int32,
            Tuple{Int32, Int32},
            x,
            y)
    end
    @test Int32(1) + Int32(2) == Int32(3)
end

# issue #11800
@test_throws ErrorException eval(Expr(:call,Core.Intrinsics.llvmcall,
    """%3 = add i32 %1, %0
       ret i32 %3""", Int32, Tuple{Int32, Int32},
        Int32(1), Int32(2))) # llvmcall must be compiled to be called

# Test whether declarations work properly
function undeclared_ceil(x::Float64)
    llvmcall("""%2 = call double @llvm.ceil.f64(double %0)
        ret double %2""", Float64, Tuple{Float64}, x)
end
@test_throws ErrorException undeclared_ceil(4.2)
@test_throws ErrorException undeclared_ceil(4.2)

function declared_floor(x::Float64)
    llvmcall(
        ("""declare double @llvm.floor.f64(double)""",
         """%2 = call double @llvm.floor.f64(double %0)
            ret double %2"""),
    Float64, Tuple{Float64}, x)
end
@test declared_floor(4.2) ≈ 4.
ir = sprint(code_llvm, declared_floor, Tuple{Float64})
@test contains(ir, "call double @llvm.floor.f64") # should be inlined

function doubly_declared_floor(x::Float64)
    llvmcall(
        ("""declare double @llvm.floor.f64(double)""",
         """%2 = call double @llvm.floor.f64(double %0)
            ret double %2"""),
    Float64, Tuple{Float64}, x+1)-1
end
@test doubly_declared_floor(4.2) ≈ 4.

function doubly_declared2_trunc(x::Float64)
    a = llvmcall(
        ("""declare double @llvm.trunc.f64(double)""",
         """%2 = call double @llvm.trunc.f64(double %0)
            ret double %2"""),
    Float64, Tuple{Float64}, x)
    b = llvmcall(
        ("""declare double @llvm.trunc.f64(double)""",
         """%2 = call double @llvm.trunc.f64(double %0)
            ret double %2"""),
    Float64, Tuple{Float64}, x+1)-1
    a + b
end
@test doubly_declared2_trunc(4.2) ≈ 8.

# Test for single line
function declared_ceil(x::Float64)
    llvmcall(
        ("declare double @llvm.ceil.f64(double)",
         """%2 = call double @llvm.ceil.f64(double %0)
            ret double %2"""),
    Float64, Tuple{Float64}, x)
end
@test declared_ceil(4.2) ≈ 5.0

# Test for multiple lines
function ceilfloor(x::Float64)
    llvmcall(
        ("""declare double @llvm.ceil.f64(double)
            declare double @llvm.floor.f64(double)""",
         """%2 = call double @llvm.ceil.f64(double %0)
            %3 = call double @llvm.floor.f64(double %2)
            ret double %3"""),
    Float64, Tuple{Float64}, x)
end
@test ceilfloor(7.4) ≈ 8.0

# Test for proper declaration extraction
function confuse_declname_parsing()
    llvmcall(
        ("""declare i64 addrspace(0)* @foobar()""",
         """ret void"""),
    Cvoid, Tuple{})
end
confuse_declname_parsing()

# Test for proper mangling of external (C) functions
function call_jl_errno()
    llvmcall(
    (""" declare i32 @jl_errno()""",
    """
    %r = call i32 @jl_errno()
    ret i32 %r
    """),Int32,Tuple{})
end
call_jl_errno()

module ObjLoadTest
    using Base: llvmcall, @ccallable
    using Test
    didcall = false
    @ccallable Cvoid function jl_the_callback()
        global didcall
        didcall = true
        nothing
    end
    # Make sure everything up until here gets compiled
    jl_the_callback(); didcall = false
    function do_the_call()
        llvmcall(
        (""" declare void @jl_the_callback()""",
        """
        call void @jl_the_callback()
        ret void
        """),Cvoid,Tuple{})
    end
    do_the_call()
    @test didcall
end

# Test for proper parenting
if Base.libllvm_version >= v"3.6" # llvm 3.6 changed the syntax for a gep, so just ignore this test on older versions
    local foo
    function foo()
        # this IR snippet triggers an optimization relying
        # on the llvmcall function having a parent module
        Base.llvmcall(
         """%1 = getelementptr i64, i64* null, i64 1
            ret void""",
        Cvoid, Tuple{})
    end
    code_llvm(DevNull, foo, ())
else
    @info "Skipping gep parentage test on llvm < 3.6"
end

module CcallableRetTypeTest
    using Base: llvmcall, @ccallable
    using Test
    @ccallable function jl_test_returns_float()::Float64
        return 42
    end
    function do_the_call()
        llvmcall(
        (""" declare double @jl_test_returns_float()""",
        """
        %1 = call double @jl_test_returns_float()
        ret double %1
        """),Float64,Tuple{})
    end
    @test do_the_call() === 42.0
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
using Test

@testset "/src/PrepareAndMeasure.jl" begin

using BellScenario: PrepareAndMeasure

@testset "num_zero_rows()" begin

    @test 0 == PrepareAndMeasure.num_zero_rows([1 0;0 1])
    @test 1 == PrepareAndMeasure.num_zero_rows([1 1;0 0])
    @test 2 == PrepareAndMeasure.num_zero_rows([0 0;0 0])

    @test 0 == PrepareAndMeasure.num_zero_rows([1 0 1;0 1 0])
    @test 1 == PrepareAndMeasure.num_zero_rows([0 0 0;1 1 1])

    @test 1 == PrepareAndMeasure.num_zero_rows([0 1;0 0;1 0])
    @test 3 == PrepareAndMeasure.num_zero_rows([1 0;0 0;0 1;0 0;0 0])

    @test 2 == PrepareAndMeasure.num_zero_rows([1 0 0]')
end

@testset "vertices_α_interpolable()" begin
    @test PrepareAndMeasure.vertices_α_interpolable([1 0 0;0 1 1;0 0 0], [1 0 0;0 1 1;0 0 0])
    @test PrepareAndMeasure.vertices_α_interpolable([1 0 0;0 1 1;0 0 0], [0 1 0;1 0 1;0 0 0])
    @test PrepareAndMeasure.vertices_α_interpolable([1 1 0;0 0 1;0 0 0], [0 0 0;1 1 1;0 0 0])
    @test PrepareAndMeasure.vertices_α_interpolable([0 0 0;1 1 1;0 0 0], [1 1 0;0 0 1;0 0 0])

    @test !PrepareAndMeasure.vertices_α_interpolable([1 0 0;0 1 1;0 0 0], [0 0 0;1 0 0;0 1 1])
    @test !PrepareAndMeasure.vertices_α_interpolable([1 0 1;0 1 0;0 0 0], [0 0 0;0 1 0;1 0 1])
end

@testset "vertices_β_interpolable()" begin
    @test PrepareAndMeasure.vertices_β_interpolable([1 0 0;0 1 1;0 0 0], [1 0 0;0 1 1;0 0 0])
    @test PrepareAndMeasure.vertices_β_interpolable([1 1 0;0 0 1;0 0 0], [0 0 0;1 1 1;0 0 0])
    @test PrepareAndMeasure.vertices_β_interpolable([0 0 0;1 1 1;0 0 0], [1 1 0;0 0 1;0 0 0])
    @test PrepareAndMeasure.vertices_β_interpolable([1 0 0;0 1 1;0 0 0], [0 0 0;1 0 0;0 1 1])
    @test PrepareAndMeasure.vertices_β_interpolable([1 0 1;0 1 0;0 0 0], [0 0 0;0 1 0;1 0 1])
    @test PrepareAndMeasure.vertices_β_interpolable([1 0 0;0 1 1;0 0 0], [0 0 0;0 0 0;1 1 1])

    @test !PrepareAndMeasure.vertices_β_interpolable([1 0 0;0 1 1;0 0 0], [0 1 0;1 0 1;0 0 0])
    @test !PrepareAndMeasure.vertices_β_interpolable([1 0 0;0 1 1;0 0 0], [0 1 0;0 0 0;1 0 1])
end

@testset "vertices_interpolable()" begin
    @test PrepareAndMeasure.vertices_interpolable([1 0 0;0 1 1;0 0 0], [0 1 0;1 0 1;0 0 0])
    @test PrepareAndMeasure.vertices_interpolable([1 0 0;0 1 1;0 0 0], [1 0 0;0 1 1;0 0 0])
    @test PrepareAndMeasure.vertices_interpolable([1 1 0;0 0 1;0 0 0], [0 0 0;1 1 1;0 0 0])

    @test !PrepareAndMeasure.vertices_interpolable([1 0 1;0 1 0;0 0 0], [1 0 0;0 0 0;0 1 1])
end

end


text:
using Random

function makeBernoulliModel(θ::Float64, ys::Vector{Int64} = Vector{Int64}(undef, 0))
  n::Int64 = length(ys)
  @inline function G(p::Int64, px::Int64Particle)
    @inbounds return px.x == ys[p] ? 1.0 : 0.0
  end
  @inline function μ(px::Int64Particle, rng::RNG) where RNG <: AbstractRNG
    px.x = rand(rng) < θ ? 1 : 0
  end
  @inline function simulate(px::Int64Particle, ::RNG) where RNG <: AbstractRNG
    return px.x
  end
  return Model(μ, G, n, simulate, Int64Particle, Int64)
end