text:
using FiniteDifferences, ForwardDiff
using LinearAlgebra: I
using StatsBase: AbstractWeights, pweights
import Manifolds: mean!, median!

include("utils.jl")

struct TestEuclidean{N} <: Manifold end
struct TestEuclideanMetric <: Metric end

Manifolds.manifold_dimension(::TestEuclidean{N}) where {N} = N
function Manifolds.local_metric(M::MetricManifold{<:TestEuclidean,<:TestEuclideanMetric}, x)
    return Diagonal(1.0:manifold_dimension(M))
end

struct TestSphere{N,T} <: Manifold
    r::T
end

struct TestSphericalMetric <: Metric end

Manifolds.manifold_dimension(::TestSphere{N}) where {N} = N
function Manifolds.local_metric(M::MetricManifold{<:TestSphere,<:TestSphericalMetric}, x)
    r = base_manifold(M).r
    d = allocate(x)
    d[1] = r^2
    d[2] = d[1] * sin(x[1])^2
    return Diagonal(d)
end
sph_to_cart(θ, ϕ) = [cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)]

struct BaseManifold{N} <: Manifold end
struct BaseManifoldMetric{M} <: Metric end
struct DefaultBaseManifoldMetric <: Metric end
struct NotImplementedMetric <: Metric end

Manifolds.manifold_dimension(::BaseManifold{N}) where {N} = N
Manifolds.inner(::BaseManifold, x, v, w) = 2 * dot(v, w)
Manifolds.exp!(::BaseManifold, y, x, v) = y .= x + 2 * v
Manifolds.log!(::BaseManifold, v, x, y) = v .= (y - x) / 2
Manifolds.project!(::BaseManifold, w, x, v) = w .= 2 .* v
Manifolds.project!(::BaseManifold, y, x) = (y .= x)
Manifolds.injectivity_radius(::BaseManifold) = Inf
Manifolds.injectivity_radius(::BaseManifold, ::Any) = Inf
Manifolds.injectivity_radius(::BaseManifold, ::AbstractRetractionMethod) = Inf
Manifolds.injectivity_radius(::BaseManifold, ::ExponentialRetraction) = Inf
Manifolds.injectivity_radius(::BaseManifold, ::Any, ::AbstractRetractionMethod) = Inf
Manifolds.injectivity_radius(::BaseManifold, ::Any, ::ExponentialRetraction) = Inf
Manifolds.local_metric(::MetricManifold{BaseManifold{N},BaseManifoldMetric{N}},x) where N = 2*one(x*x')
Manifolds.exp!(M::MetricManifold{BaseManifold{N},BaseManifoldMetric{N}}, y, x, v) where N = exp!(base_manifold(M), y, x, v)
Manifolds.vector_transport_to!(::BaseManifold, vto, x, v, y, ::ParallelTransport) = (vto .= v)
Manifolds.get_basis(::BaseManifold{N},x,B::DefaultOrthonormalBasis) where {N} = CachedBasis(B, [(Matrix{eltype(x)}(I, N, N)[:,i]) for i in 1:N])
Manifolds.get_coordinates!(::BaseManifold, Y, p, X, ::DefaultOrthonormalBasis) = (Y .= X)
Manifolds.get_vector!(::BaseManifold, Y, p, X, ::DefaultOrthonormalBasis) = (Y .= X)
Manifolds.default_metric_dispatch(::BaseManifold, ::DefaultBaseManifoldMetric) = Val(true)
Manifolds.projected_distribution(M::BaseManifold, d) = ProjectedPointDistribution(M, d, project!, rand(d))
Manifolds.projected_distribution(M::BaseManifold, d, x) = ProjectedPointDistribution(M, d, project!, x)
Manifolds.mean!(::BaseManifold, y, x::AbstractVector, w::AbstractVector; kwargs...) = fill!(y, 1)
Manifolds.median!(::BaseManifold, y, x::AbstractVector, w::AbstractVector; kwargs...) = fill!(y, 2)
Manifolds.mean!(::MetricManifold{BaseManifold{N},BaseManifoldMetric{N}}, y, x::AbstractVector, w::AbstractVector; kwargs...) where {N} = fill!(y, 3)
Manifolds.median!(::MetricManifold{BaseManifold{N},BaseManifoldMetric{N}}, y, x::AbstractVector, w::AbstractVector; kwargs...) where {N} = fill!(y, 4)

function Manifolds.flat!(::BaseManifold, v::FVector{Manifolds.CotangentSpaceType}, x, w::FVector{Manifolds.TangentSpaceType})
    v.data .= 2 .* w.data
    return v
end
function Manifolds.sharp!(
    ::BaseManifold,
    v::FVector{Manifolds.TangentSpaceType},
    x,
    w::FVector{Manifolds.CotangentSpaceType},
)
    v.data .= w.data ./ 2
    return v
end

@testset "Metrics" begin
    # some tests failed due to insufficient accuracy for a particularly bad RNG state
    Random.seed!(42)
    @testset "Metric Basics" begin
        #one for MetricManifold, one for Manifold & Metric
        @test length(methods(is_default_metric)) == 2
    end

    @testset "solve_exp_ode error message" begin
        E = TestEuclidean{3}()
        g = TestEuclideanMetric()
        M = MetricManifold(E, g)

        x = [1.0, 2.0, 3.0]
        v = [2.0, 3.0, 4.0]
        @test_throws ErrorException exp(M, x, v)
        if VERSION ≥ v"1.1"
            using OrdinaryDiffEq
            exp(M, x, v)
        end
    end
    @testset "Local Metric Error message" begin
        M = MetricManifold(BaseManifold{2}(), NotImplementedMetric())
        @test_throws ErrorException local_metric(M, [3, 4])
    end
    @testset "scaled Euclidean metric" begin
        n = 3
        E = TestEuclidean{n}()
        g = TestEuclideanMetric()
        M = MetricManifold(E, g)
        @test repr(M) == "MetricManifold(TestEuclidean{3}(), TestEuclideanMetric())"

        if VERSION ≥ v"1.3"
            @test E |> TestEuclideanMetric() === M
            @test E |> TestEuclideanMetric === M
        end

        G = Diagonal(1.0:n)
        invG = inv(G)
        @test manifold_dimension(M) == n
        @test base_manifold(M) === E
        @test metric(M) === g

        @test_throws ErrorException local_metric_jacobian(E, zeros(3))
        @test_throws ErrorException christoffel_symbols_second_jacobian(E, zeros(3))

        for vtype in (Vector, MVector{n})
            x, v, w = vtype(randn(n)), vtype(randn(n)), vtype(randn(n))

            @test check_manifold_point(M, x) == check_manifold_point(E, x)
            @test check_tangent_vector(M, x, v) == check_tangent_vector(E, x, v)

            @test local_metric(M, x) ≈ G
            @test inverse_local_metric(M, x) ≈ invG
            @test det_local_metric(M, x) ≈ *(1.0:n...)
            @test log_local_metric_density(M, x) ≈ sum(log.(1.0:n)) / 2
            @test inner(M, x, v, w) ≈ dot(v, G * w) atol = 1e-6
            @test norm(M, x, v) ≈ sqrt(dot(v, G * v)) atol = 1e-6

            if VERSION ≥ v"1.1"
                T = 0:0.5:10
                @test geodesic(M, x, v, T) ≈ [x + t * v for t in T] atol = 1e-6
            end

            @test christoffel_symbols_first(M, x) ≈ zeros(n, n, n) atol = 1e-6
            @test christoffel_symbols_second(M, x) ≈ zeros(n, n, n) atol = 1e-6
            @test riemann_tensor(M, x) ≈ zeros(n, n, n, n) atol = 1e-6
            @test ricci_tensor(M, x) ≈ zeros(n, n) atol = 1e-6
            @test ricci_curvature(M, x) ≈ 0 atol = 1e-6
            @test gaussian_curvature(M, x) ≈ 0 atol = 1e-6
            @test einstein_tensor(M, x) ≈ zeros(n, n) atol = 1e-6

            fdm = forward_fdm(2, 1)
            @test christoffel_symbols_first(M, x; backend=fdm) ≈ zeros(n, n, n) atol=1e-6
            @test christoffel_symbols_second(M, x; backend=fdm) ≈ zeros(n, n, n) atol=1e-6
            @test riemann_tensor(M, x; backend=fdm) ≈ zeros(n, n, n, n) atol=1e-6
            @test ricci_tensor(M, x; backend=fdm) ≈ zeros(n, n) atol=1e-6
            @test ricci_curvature(M, x; backend=fdm) ≈ 0 atol=1e-6
            @test gaussian_curvature(M, x; backend=fdm) ≈ 0 atol=1e-6
            @test einstein_tensor(M, x; backend=fdm) ≈ zeros(n, n) atol=1e-6

            @test christoffel_symbols_first(M, x; backend=:forwarddiff) ≈ zeros(n, n, n) atol=1e-6
            @test christoffel_symbols_second(M, x; backend=:forwarddiff) ≈ zeros(n, n, n) atol=1e-6
            @test riemann_tensor(M, x; backend=:forwarddiff) ≈ zeros(n, n, n, n) atol=1e-6
            @test ricci_tensor(M, x; backend=:forwarddiff) ≈ zeros(n, n) atol=1e-6
            @test ricci_curvature(M, x; backend=:forwarddiff) ≈ 0 atol=1e-6
            @test gaussian_curvature(M, x; backend=:forwarddiff) ≈ 0 atol=1e-6
            @test einstein_tensor(M, x; backend=:forwarddiff) ≈ zeros(n, n) atol=1e-6
        end
    end

    @testset "scaled Sphere metric" begin
        n = 2
        r = 10 * rand()
        θ, ϕ = π * rand(), 2π * rand()
        Sr = TestSphere{n,Float64}(r)
        S = Manifolds.Sphere(n)
        g = TestSphericalMetric()
        M = MetricManifold(Sr, g)

        @test manifold_dimension(M) == n
        @test base_manifold(M) === Sr
        @test metric(M) === g

        for vtype in (Vector, MVector{n})
            x = vtype([θ, ϕ])
            G = Diagonal(vtype([1, sin(θ)^2])) .* r^2
            invG = Diagonal(vtype([1, 1 / sin(θ)^2])) ./ r^2
            v, w = normalize(randn(n)), normalize(randn(n))

            @test local_metric(M, x) ≈ G atol = 1e-6
            @test inverse_local_metric(M, x) ≈ invG atol = 1e-6
            @test det_local_metric(M, x) ≈ r^4 * sin(θ)^2 atol = 1e-6
            @test log_local_metric_density(M, x) ≈ 2 * log(r) + log(sin(θ)) atol = 1e-6
            @test inner(M, x, v, w) ≈ dot(v, G * w) atol = 1e-6
            @test norm(M, x, v) ≈ sqrt(dot(v, G * v)) atol = 1e-6

            xcart = sph_to_cart(θ, ϕ)
            vcart =
                [
                    cos(ϕ) * cos(θ) -sin(ϕ) * sin(θ)
                    sin(ϕ) * cos(θ) cos(ϕ) * sin(θ)
                    -sin(θ) 0
                ] * v

            if VERSION ≥ v"1.1" && (!Sys.iswindows() || Sys.ARCH == :x86_64)
                @testset "numerically integrated geodesics for $vtype" begin
                    T = 0:0.1:1
                    @test isapprox(
                        [sph_to_cart(yi...) for yi in geodesic(M, x, v, T)],
                        geodesic(S, xcart, vcart, T);
                        atol = 1e-3,
                        rtol = 1e-3,
                    )
                end
            end

            Γ₁ = christoffel_symbols_first(M, x)
            for i = 1:n, j = 1:n, k = 1:n
                if (i, j, k) == (1, 2, 2) || (i, j, k) == (2, 1, 2)
                    @test Γ₁[i, j, k] ≈ r^2 * cos(θ) * sin(θ) atol = 1e-6
                elseif (i, j, k) == (2, 2, 1)
                    @test Γ₁[i, j, k] ≈ -r^2 * cos(θ) * sin(θ) atol = 1e-6
                else
                    @test Γ₁[i, j, k] ≈ 0 atol = 1e-6
                end
            end

            Γ₂ = christoffel_symbols_second(M, x)
            for l = 1:n, i = 1:n, j = 1:n
                if (l, i, j) == (1, 2, 2)
                    @test Γ₂[l, i, j] ≈ -cos(θ) * sin(θ) atol = 1e-6
                elseif (l, i, j) == (2, 1, 2) || (l, i, j) == (2, 2, 1)
                    @test Γ₂[l, i, j] ≈ cot(θ) atol = 1e-6
                else
                    @test Γ₂[l, i, j] ≈ 0 atol = 1e-6
                end
            end

            R = riemann_tensor(M, x)
            for l = 1:n, i = 1:n, j = 1:n, k = 1:n
                if (l, i, j, k) == (2, 1, 1, 2)
                    @test R[l, i, j, k] ≈ -1 atol = 2e-6
                elseif (l, i, j, k) == (2, 1, 2, 1)
                    @test R[l, i, j, k] ≈ 1 atol = 2e-6
                elseif (l, i, j, k) == (1, 2, 1, 2)
                    @test R[l, i, j, k] ≈ sin(θ)^2 atol = 1e-6
                elseif (l, i, j, k) == (1, 2, 2, 1)
                    @test R[l, i, j, k] ≈ -sin(θ)^2 atol = 1e-6
                else
                    @test R[l, i, j, k] ≈ 0 atol = 1e-6
                end
            end

            @test ricci_tensor(M, x) ≈ G ./ r^2 atol = 2e-6
            @test ricci_curvature(M, x) ≈ 2 / r^2 atol = 2e-6
            @test gaussian_curvature(M, x) ≈ 1 / r^2 atol = 2e-6
            @test einstein_tensor(M, x) ≈ ricci_tensor(M, x) - gaussian_curvature(M, x) .* G atol =
                1e-6
        end
    end

    @testset "Metric decorator" begin
        M = BaseManifold{3}()
        g = BaseManifoldMetric{3}()
        MM = MetricManifold(M, g)
        if VERSION ≥ v"1.3"
            @test DefaultBaseManifoldMetric(BaseManifold{3}()) === MetricManifold(BaseManifold{3}(),DefaultBaseManifoldMetric())
            MT = DefaultBaseManifoldMetric()
            @test MT(BaseManifold{3}()) === MetricManifold(BaseManifold{3}(),DefaultBaseManifoldMetric())
        end
        g2 = DefaultBaseManifoldMetric()
        MM2 = MetricManifold(M, g2)

        @test (@inferred Manifolds.default_metric_dispatch(MM)) ===
              (@inferred Manifolds.default_metric_dispatch(base_manifold(MM), metric(MM)))
        @test (@inferred Manifolds.default_metric_dispatch(MM2)) ===
              (@inferred Manifolds.default_metric_dispatch(base_manifold(MM2), metric(MM2)))
        @test (@inferred Manifolds.default_metric_dispatch(MM2)) === Val(true)
        @test is_default_metric(MM) == is_default_metric(base_manifold(MM), metric(MM))
        @test is_default_metric(MM2) == is_default_metric(base_manifold(MM2), metric(MM2))
        @test is_default_metric(MM2)
        @test Manifolds.default_decorator_dispatch(MM) === Manifolds.default_metric_dispatch(MM)
        @test Manifolds.default_decorator_dispatch(MM2) === Manifolds.default_metric_dispatch(MM2)

        @test convert(typeof(MM2), M) == MM2
        @test_throws ErrorException convert(typeof(MM), M)
        x = [0.1, 0.2, 0.4]
        v = [0.5, 0.7, 0.11]
        w = [0.13, 0.17, 0.19]
        y = allocate(x)

        @test inner(M, x, v, w) == 2 * dot(v, w)
        @test inner(MM, x, v, w) === inner(M, x, v, w)
        @test norm(MM, x, v) === norm(M, x, v)
        @test exp(M, x, v) == x + 2 * v
        @test exp(MM2, x, v) == exp(M, x, v)
        @test exp!(MM, y, x, v) === exp!(M, y, x, v)
        @test retract!(MM, y, x, v) === retract!(M, y, x, v)
        @test retract!(MM, y, x, v, 1) === retract!(M, y, x, v, 1)
        # without a definition for the metric from the embedding, no projection possible
        @test_throws ErrorException log!(MM, w, x, y) === project!(M, w, x, y)
        @test_throws ErrorException project!(MM, w, x, v) ===
                                    project!(M, w, x, v)
        @test_throws ErrorException project!(MM, y, x) === project!(M, y, x)
        @test_throws ErrorException vector_transport_to!(MM, w, x, v, y) ===
                                    vector_transport_to!(M, w, x, v, y)
        # without DiffEq, these error
        # @test_throws ErrorException exp(MM,x, v, 1:3)
        # @test_throws ErrorException exp!(MM, y, x, v)
        # these always fall back anyways.
        @test zero_tangent_vector!(MM, v, x) === zero_tangent_vector!(M, v, x)

        @test injectivity_radius(MM, x) === injectivity_radius(M, x)
        @test injectivity_radius(MM) === injectivity_radius(M)
        @test injectivity_radius(MM, ProjectionRetraction()) === injectivity_radius(M, ProjectionRetraction())
        @test injectivity_radius(MM, ExponentialRetraction()) === injectivity_radius(M, ExponentialRetraction())
        @test injectivity_radius(MM) === injectivity_radius(M)

        @test is_manifold_point(MM, x) === is_manifold_point(M, x)
        @test is_tangent_vector(MM, x, v) === is_tangent_vector(M, x, v)

        @test_throws ErrorException local_metric(MM2, x)
        @test_throws ErrorException local_metric_jacobian(MM2, x)
        @test_throws ErrorException christoffel_symbols_second_jacobian(MM2, x)
        # MM falls back to nondefault error
        @test_throws ErrorException projected_distribution(MM, 1, x)
        @test_throws ErrorException projected_distribution(MM, 1)
        # normal_tvector is not implemented
        @test_throws ErrorException normal_tvector_distribution(MM, x, 0.2)

        @test inner(MM2, x, v, w) === inner(M, x, v, w)
        @test norm(MM2, x, v) === norm(M, x, v)
        @test distance(MM2, x, y) === distance(M, x, y)
        @test exp!(MM2, y, x, v) === exp!(M, y, x, v)
        @test exp(MM2, x, v) == exp(M, x, v)
        @test log!(MM2, v, x, y) === log!(M, v, x, y)
        @test log(MM2, x, y) == log(M, x, y)
        @test retract!(MM2, y, x, v) === retract!(M, y, x, v)
        @test retract!(MM2, y, x, v, 1) === retract!(M, y, x, v, 1)

        @test project!(MM2, y, x) === project!(M, y, x)
        @test project!(MM2, w, x, v) === project!(M, w, x, v)
        @test vector_transport_to!(MM2, w, x, v, y) == vector_transport_to!(M, w, x, v, y)
        @test zero_tangent_vector!(MM2, v, x) === zero_tangent_vector!(M, v, x)
        @test injectivity_radius(MM2, x) === injectivity_radius(M, x)
        @test injectivity_radius(MM2) === injectivity_radius(M)
        @test injectivity_radius(MM2, x, ExponentialRetraction()) === injectivity_radius(M, x, ExponentialRetraction())
        @test injectivity_radius(MM2, ExponentialRetraction()) === injectivity_radius(M, ExponentialRetraction())
        @test injectivity_radius(MM2, x, ProjectionRetraction()) === injectivity_radius(M, x, ProjectionRetraction())
        @test injectivity_radius(MM2, ProjectionRetraction()) === injectivity_radius(M, ProjectionRetraction())
        @test is_manifold_point(MM2, x) === is_manifold_point(M, x)
        @test is_tangent_vector(MM2, x, v) === is_tangent_vector(M, x, v)

        a = Manifolds.projected_distribution(M, Distributions.MvNormal(zero(zeros(3)), 1.0))
        b = Manifolds.projected_distribution(
            MM2,
            Distributions.MvNormal(zero(zeros(3)), 1.0),
        )
        @test isapprox(Matrix(a.distribution.Σ), Matrix(b.distribution.Σ))
        @test isapprox(a.distribution.μ, b.distribution.μ)
        @test get_basis(M, x, DefaultOrthonormalBasis()).data ==
              get_basis(MM2, x, DefaultOrthonormalBasis()).data
        @test_throws ErrorException get_basis(MM, x, DefaultOrthonormalBasis())
        cov = flat(M, x, FVector(TangentSpace, v))
        cow = flat(M, x, FVector(TangentSpace, w))
        @test cov.data ≈ flat(MM, x, FVector(TangentSpace, v)).data
        cotspace = CotangentBundleFibers(M)
        cotspace2 = CotangentBundleFibers(MM)
        @test cov.data ≈ 2 * v
        @test inner(M, x, v, w) ≈ inner(cotspace, x, cov.data, cow.data)
        @test inner(MM, x, v, w) ≈ inner(cotspace, x, cov.data, cow.data)
        @test inner(MM, x, v, w) ≈ inner(cotspace2, x, cov.data, cow.data)
        @test sharp(M, x, cov).data ≈ v

        xsample = [[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]
        w = pweights([0.5, 0.5])
        # test despatch with results from above
        @test mean(M, xsample, w) ≈ ones(3)
        @test mean(MM2, xsample, w) ≈ ones(3)
        @test mean(MM, xsample, w) ≈ 3 .* ones(3)

        @test median(M, xsample, w) ≈ 2 .* ones(3)
        @test median(MM2, xsample, w) ≈ 2 * ones(3)
        @test median(MM, xsample, w) ≈ 4 .* ones(3)
    end

    @testset "Metric decorator dispatches" begin
        M = BaseManifold{3}()
        g = BaseManifoldMetric{3}()
        MM = MetricManifold(M, g)
        x = [1,2,3]
        # nonmutating always go to parent for allocation
        for f in [exp, flat, inverse_retract, log, mean, median, project]
            @test Manifolds.decorator_transparent_dispatch(f, MM) === Val{:parent}()
        end
        for f in [sharp, retract, get_vector, get_coordinates]
            @test Manifolds.decorator_transparent_dispatch(f, MM) === Val{:parent}()
        end
        for f in [vector_transport_along, vector_transport_direction, vector_transport_to]
            @test Manifolds.decorator_transparent_dispatch(f,MM) === Val{:parent}()
        end
        for f in [get_basis, inner, normal_tvector_distribution, projected_distribution]
            @test Manifolds.decorator_transparent_dispatch(f,MM) === Val{:intransparent}()
        end
        for f in [get_coordinates!, get_vector!]
            @test Manifolds.decorator_transparent_dispatch(f,MM) === Val{:intransparent}()
        end

        # mirroring ones are mostly intransparent despite for a few cases - e.g. dispatch/default last variables
        for f in [exp!, flat!, inverse_retract!, log!, mean!, median!]
            @test Manifolds.decorator_transparent_dispatch(f,MM) === Val{:intransparent}()
        end
        for f in [norm, project!, sharp!, retract!]
            @test Manifolds.decorator_transparent_dispatch(f,MM) === Val{:intransparent}()
        end
        for f in [vector_transport_along!, vector_transport_direction!, vector_transport_to!]
            @test Manifolds.decorator_transparent_dispatch(f,MM) === Val{:intransparent}()
        end

        @test Manifolds.decorator_transparent_dispatch(exp!,MM,x,x,x,x) === Val{:parent}()
        @test Manifolds.decorator_transparent_dispatch(inverse_retract!, MM, x, x, x, LogarithmicInverseRetraction()) === Val{:parent}()
        @test Manifolds.decorator_transparent_dispatch(retract!, MM, x, x, x, ExponentialRetraction()) === Val{:parent}()

    end
end


text:
export write_conll

function write_conll(dir::String)
    lines = String[]
    for file in readdir(dir)
        endswith(file, ".pdfanno") || continue
        pdffile = file[1:end-4]
        println(pdffile)
        pdffile = "$dir/$pdffile"
        try
            l = format_conll(pdffile)
            append!(lines, l)
        catch
        end
        push!(lines, "")
    end
    open("$dir.conll", "w") do io
        for line in lines
            println(io, line)
        end
    end
end

function format_conll(pdffile::String)
    chars = readpdf(pdffile)

    pdfanno = readpdfanno(pdffile * "anno")
    spans = pdfanno.spans
    label_s = "abstract_start"
    label_e = "abstract_end"
    i = findfirst(s -> s.label == label_s, spans)
    ind_s = start(spans[i].range)
    i = findfirst(s -> s.label == label_e, spans)
    ind_e = start(spans[i].range)

    texts = filter(c -> c._type == "TEXT", chars[ind_s:ind_e])
    tokens = tokenize_space(texts)
    char2tokenid = Dict()
    for i = 1:length(tokens)
        for c in tokens[i].chars
            char2tokenid[c] = i
        end
    end

    labels = ["O" for i=1:length(tokens)]
    for span in spans
        (span.label == label_s || span.label == label_e) && continue
        s = start(span.range)
        l = last(span.range)
        label = span.label
        startswith(label,"Material_") && (label = "Material")
        bid = char2tokenid[chars[s]]
        eid = char2tokenid[chars[l]]
        if bid == eid
            labels[bid] = "S-$label"
        else
            labels[bid] = "B-$label"
            labels[eid] = "E-$label"
            for i = bid+1:eid-1
                labels[i] = "I-$label"
            end
        end
    end

    lines = String[]
    buffer = []
    b = false
    for i = 1:length(tokens)
        t = tokens[i]
        l = labels[i]
        l == "O" || (b = true)
        push!(buffer, (t,l))
        if t.str[end] == '.' || i == length(tokens)
            if b
                for (_t,_l) in buffer
                    str = _t.str
                    str[end] == '.' && length(str) > 1 && (str = str[1:end-1])
                    push!(lines, "$str\t$_l")
                end
                push!(lines, "")
            end
            b = false
            empty!(buffer)
        end
    end
    lines
end

ext: jl


text:
using Koma, MAT, JLD2
using MRIReco
#Phantom
phantom = brain_phantom2D()
phantom.Δw .*= 0 #Removing off-resonance
#Scanner
sys = Scanner()
#Seq: Excitation and radial acquisition
TB1 = 1e-3
EX = PulseDesigner.RF_hard(1/(2π*γ*TB1),TB1,sys) #α = γ B1 T
FOV, Nr = 23e-2, 100
RAD = PulseDesigner.radial_base(FOV, Nr, sys)
#MRF parameters
path = @__DIR__
data = matread(path*"/MRFinput.mat")
α, ϕ, TE, TR = data["FA"]/180*π, data["PH"]/180*π, data["TE"]*1e-3, data["TR"]*1e-3
A = α.*exp.(1im*ϕ)
#Inversion pulse
TINV = 50e-3
INV = (π+0im)*EX + Delay(TINV)
#Delays
Ta = dur(RAD)
delayTE = [Delay(TE[n]-Ta) for n=1:length(α)]
delayTR = [Delay(TR[n]-TE[n]-Ta) for n=1:length(α)]
#MRF with rotated spokes
φ, Nφ = (√5 + 1)/2, 7;
Δθ = π/(φ+Nφ-1) # Tiny golden angle with Nφ = 7
NTRs = 200 #Number of TRs
seq = INV + sum([A[n]*EX + delayTE[n] + rotz((n-1)*Δθ)*RAD + delayTR[n] for n=1:NTRs])
seq.DEF["Nz"] = NTRs
jldsave("./examples/3.mrf_ismrm/mrf.seqk"; seq=seq)
# plot_seq(seq)
## Simulation
fingerprint = simulate(phantom, seq, sys); #This takes like 5 min for NTRs = 500.
## Output ISMRMRD
raw_ismrmrd = Koma.rawSignalToISMRMRD([fingerprint;;],seq;phantom,sys)
fname = "MRF_signal_Δθ_$(floor(Int64, Δθ/π*180))_NTRs_$NTRs"
fout = ISMRMRDFile("./examples/3.mrf_ismrm/$fname.h5")
save(fout, raw_ismrmrd)
## Output MAT
# matwrite("./examples/3.mrf_ismrm/$fname.mat", Dict("signal"=>fingerprint))


text:
using SafeTestsets

@safetestset "Utility Tests" begin

    @safetestset "find_index" begin
        using ACTRModels, Test
        include("Utility_Functions.jl")

        chunks = populateMemory()
        declarative = Declarative(memory=chunks)
        actr = ACTR(declarative = declarative)
        idx = find_index(actr; attribute = :wings)
        @test chunks[idx].slots.attribute == :wings
        @test chunks[idx].slots.object == :bird
        @test chunks[idx].slots.value == :True

        idx = find_index(actr, ==, !=; attribute = :locomotion, value = :swimming)
        @test chunks[idx].slots.value == :flying
        @test chunks[idx].slots.object == :bird
        @test chunks[idx].slots.attribute == :locomotion
    end

    @safetestset "find_indices" begin
        using ACTRModels, Test
        include("Utility_Functions.jl")

        chunks = populateMemory()
        declarative = Declarative(memory=chunks)
        actr = ACTR(declarative = declarative)
        idx = find_indices(actr; attribute = :locomotion, value = :swimming)
        @test length(idx) == 3
        chunk_set = chunks[idx]
        @test all(x->x.slots.attribute == :locomotion && x.slots.value == :swimming, chunk_set)

        idx = find_indices(actr, ==, !=; attribute = :category, value = :bird)
        @test length(idx) == 4
        chunk_set = chunks[idx]
        @test all(x->x.slots.attribute == :category && x.slots.value != :bird, chunk_set)
    end
end


text:
module steady_state
function logit(x, fixedParam)
    return 1 ./ (1.0 .+ exp.(-fixedParam.lambda .* x))
end
#HC evolution function, g for gender, a for age
function geth_I(x,y,geg,hga)
    #gea is the probability of hc growth for gender g
    #hga is the human capital distribution of gender g in stage a
    #get the job productivity level from y, an integer from 1 to 7
    jobLevel = y%7 == 0 ? 7 : y%7

    hga_I_xy = hga[x,y]*(1-geg[jobLevel])
    #check human capital level, if 1, human capital measure decreases
    if(x%7==1)
        return hga_I_xy
    elseif(x%7==0)#check human capital level, if 0, human capital measure increases
        hga_I_xy = hga[x,y] + hga[x-1,y]*geg[jobLevel]
    else #check human capital level, there's flowin and flow out
        hga_I_xy = hga_I_xy + hga[x-1,y]*geg[jobLevel]
    end
    return hga_I_xy
end
#Searching and Matching Stage
#Component of change that is attributable to change in employment
#A for age
#g for gender
function getEmpChange(A,i,j,ugA_I,SgA,hgA_I,fxp,v)
    dgA = -SgA[i,:] .+ SgA[i,j];
    IgA = Int.(dgA .!= 0);
    #determine the age and the associated parameter
    if (cmp(A,"NC")==0)
        lambdau_a = fxp.lambdau_NC;
    else
        lambdau_a = fxp.lambdau_YC;
    end
    part1 = lambdau_a*ugA_I[i]*v[j];
    part2 = fxp.lambdae*v[j]*sum(hgA_I[i,:].*IgA.*logit(dgA, fxp).*(SgA[i,:].>0));
    part3 = -fxp.lambdae*hgA_I[i,j]*sum(v.*IgA.*logit(-dgA, fxp));
    return part1+part2+part3;
end
#get the change in distribution for human capital at NC
function get_hgNC_II_ij(g,i,j,ugNC_I,SgNC,hgNC_I,fxp,v)
    empChange = getEmpChange("NC",i,j,ugNC_I,SgNC,hgNC_I,fxp,v);
     #determine the gender and the associated parameters
    deltag_NC =  g == "m" ? fxp.deltam_NC : fxp.deltaf_NC;
    return hgNC_I[i,j]+(SgNC[i,j]>0)*(empChange)- hgNC_I[i,j]*(fxp.gamma + fxp.chi + deltag_NC);
end
#get the change in distribution for human capital at YC
function get_hgYC_II_ij(g,i,j,ugYC_I,SgYC,hgYC_I,hgPL_I,fxp,v)
    empChange = getEmpChange("YC",i,j,ugYC_I,SgYC,hgYC_I,fxp,v);
    #determine the gender and the associated parameters
    deltag_YC =  g == "m" ? fxp.deltam_YC : fxp.deltaf_YC;
    etag =  g == "m" ? fxp.etam : fxp.etaf;
    return hgYC_I[i,j]+etag*hgPL_I[i,j]+(SgYC[i,j]>0)*(empChange)- hgYC_I[i,j]*(fxp.gamma + deltag_YC + fxp.chi);
end
#get the change in distribution for human capital at D
function get_hgD_II_ij(i,j,SgD,ugD_I,hgNC_I,hgYC_I,hgPL_I,hgD_I,fxp,v)
    empChange = getEmpChange("D",i,j,ugD_I,SgD,hgD_I,fxp,v);
    return hgD_I[i,j]+ fxp.gamma*(hgNC_I[i,j]+ hgYC_I[i,j]+hgPL_I[i,j]) + (SgD[i,j]>0)*(empChange) - hgD_I[i,j]*(fxp.phi + fxp.delta);
end
export get_hgD_II_ij,get_hgYC_II_ij,get_hgNC_II_ij,geth_I
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module PowerModelsACDC

# import Compat
import JuMP
import Memento
import PowerModels
const _PM = PowerModels
import InfrastructureModels
# import InfrastructureModels: ids, ref, var, con, sol, nw_ids, nws, optimize_model!, @im_fields
const _IM = InfrastructureModels

import JuMP: with_optimizer
export with_optimizer

# Create our module level logger (this will get precompiled)
const _LOGGER = Memento.getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via `getlogger(PowerModels)`
# NOTE: If this line is not included then the precompiled `_PM._LOGGER` won't be registered at runtime.
__init__() = Memento.register(_LOGGER)

include("prob/acdcopf.jl")
include("prob/acdcpf.jl")
include("prob/acdcopf_bf.jl")
include("prob/tnepopf.jl")
include("prob/acdctnepopf.jl")
include("prob/tnepopf_bf.jl")
include("prob/acdctnepopf_bf.jl")
include("prob/mp_tnepopf.jl")
include("prob/mp_acdctnepopf.jl")
include("prob/mp_tnepopf_bf.jl")
include("prob/mp_acdctnepopf_bf.jl")

include("core/solution.jl")
include("core/data.jl")
include("core/variabledcgrid.jl")
include("core/variableconv.jl")
include("core/base.jl")
include("core/constraint.jl")
include("core/objective.jl")
include("core/relaxation_scheme.jl")

include("formdcgrid/acp.jl")
include("formdcgrid/dcp.jl")
include("formdcgrid/wr.jl")
include("formdcgrid/wrm.jl")
include("formdcgrid/bf.jl")
include("formdcgrid/lpac.jl")
include("formdcgrid/shared.jl")

include("formconv/acp.jl")
include("formconv/dcp.jl")
include("formconv/wr.jl")
include("formconv/wrm.jl")
include("formconv/bf.jl")
include("formconv/lpac.jl")
include("formconv/shared.jl")

include("core/constraint_template.jl")
include("io/multinetwork.jl")
include("io/results.jl")
end


text:
using CUDA

# Detect whether there is a display available
display_mode = true
try
   run(pipeline(`xset -q`, `tail -n1`, `grep "Monitor is On"`));
catch e
    global display_mode = false
end
if display_mode
    using ImageView, Images
end

function visualise(image::Array{Float32,3})
    display_mode ? imshow(image) : println("No display available")
end

function visualise(image::Array{Float32,1}, DIMX::Integer, DIMY::Integer, DIMZ::Integer)
    output = zeros(Float32, DIMX, DIMY, DIMZ);
    for i = 1:DIMX
        for j = 1:DIMY
            for k = 1:DIMZ
                index = calculate_index(i, j, k, DIMX, DIMY, DIMZ, 4)
                output[i, j, k] = image[index]
            end
        end
    end
    display_mode ? imshow(output) : println("No display available")
end

function visualise(c_image::CuArray{Float32,3})
    image = Array(c_image)
    display_mode ? imshow(image) : println("No display available")
end

function visualise(c_image::CuArray{Float32,1}, DIMX::Integer, DIMY::Integer, DIMZ::Integer)
    image = Array(c_image);
    visualise(image, DIMX, DIMY, DIMZ)
end

function save_image(image::Array{Float32,1})
    io = open("recon.raw","w")
    for i in eachindex(image)
        write(io, image[i])
    end
    close(io)
end

function save_image(image::Array{Float32,3})
    save_image(image[:])
end

function save_image(image::Array{Float32,1}, DIMX::Integer, DIMY::Integer, DIMZ::Integer)
    output = zeros(Float32, DIMX, DIMY, DIMZ);
    for i = 1:DIMX
        for j = 1:DIMY
            for k = 1:DIMZ
                index = calculate_index(i, j, k, DIMX, DIMY, DIMZ, 4)
                output[i, j, k] = image[index]
            end
        end
    end
    save_image(output)
end

function save_image(c_image::CuArray{Float32,1}, DIMX::Integer, DIMY::Integer, DIMZ::Integer)
    image = Array(c_image);
    save_image(image, DIMX, DIMY, DIMZ)
end

max_stars_repo_name: michielvangendt/PET-Julia

max_issues_repo_name: michielvangendt/PET-Julia

max_forks_repo_name: michielvangendt/PET-Julia


text:
using SymbolicRegression, BenchmarkTools, Random
using SymbolicRegression: eval_tree_array

const SUITE = BenchmarkGroup()

options = Options(; binary_operators=(+, -, /, *), unary_operators=(cos, exp))
tree = Node(
    2,
    Node(
        1,
        Node(3, Node(1, Node(1.0f0), Node(2)), Node(2, Node(-1.0f0))),
        Node(1, Node(3), Node(4)),
    ),
    Node(
        4,
        Node(3, Node(1, Node(1.0f0), Node(2)), Node(2, Node(-1.0f0))),
        Node(1, Node(3), Node(4)),
    ),
)

_X = randn(MersenneTwister(0), Float32, 5, 1000)

SUITE["evaluation"] = BenchmarkGroup()
for T in [Float32, Float64, BigFloat]
    X = map(x -> convert(T, x), _X)
    f = eval_tree_array
    SUITE["evaluation"][string(T)] = @benchmarkable ($f)(tree, $X, options)
end


text:
module Irr
    export sum_fracts

    function sum_fracts(l::Array{Array{Int,1},1})
        result = sum( x -> x[1]//x[2], l )
        if denominator(result) == 1 || numerator(result) == 0
            numerator(result)
        else
            [numerator(result),denominator(result)]
        end
    end

    sum_fracts(l::Any) = nothing

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
export learningUpdate!

function learningUpdate!(evs::dcdp_Emax, t::dcdp_tmpvars, p::dcdp_primitives, σ::Real, dograd::Bool)

    Πψtmp  = t.Πψtmp
    ψspace = _ψspace(p)

    d2plus = 2:_dmax(p.wp)+1  # TODO: would be good to make this more general.
    lrn2inf = inf_fm_lrn(p.wp)
    exp2lrn = exploratory_learning(p.wp)

    @views EV1   =  evs.EV[:,:,  lrn2inf]
    @views dEV1  = evs.dEV[:,:,:,lrn2inf]

    @views EV2   =  evs.EV[ :,:,  exp2lrn]
    @views dEV2  = evs.dEV[ :,:,:,exp2lrn]
    @views dEVσ2 = evs.dEVσ[:,:  ,exp2lrn]

    # ubVtilde = u[:,:,2:dmaxp1] + β * Πψ ⊗ I * EV[:,:,2:dmaxp1]
    _βΠψ!(Πψtmp, ψspace, σ, p.β)
    A_mul_B_md!(EV2, Πψtmp, EV1, 2)

    if dograd
        # dubVtilde/dθ = du/dθ[:,:,:,2:dmaxp1] + β * Πψ ⊗ I * dEV/dθ[:,:,:,2:dmaxp1]
        A_mul_B_md!(dEV2, Πψtmp, dEV1, 2)

        # ∂EVtilde/∂σ[:,:,2:dmaxp1] = ∂u/∂σ[:,:,2:dmaxp1] + β * dΠψ/dσ ⊗ I * EV[:,:,2:dmaxp1]
        _βΠψdθρ!(Πψtmp, ψspace, σ, p.β)
        A_mul_B_md!(dEVσ2, Πψtmp, EV1, 2)
    end
end


text:
"""
Implementation of the Kronecker Graphical Lasso algorithm for tensor graphical model estimation

Author: Wayne Wang
Last modified: 07/29/2021
"""

function kglasso(X::AbstractArray{<:Real}; Ψ0_list::AbstractArray = Array{Array{Float64, 2}, 1}(),
    λ_list::AbstractArray = Array{Float64, 1}(), N_iter::Int = 2, ϵ::Real=1e-5)
    """
    The Kronecker graphical lasso algorithm, which applies the Glasso algorithm alternatively
    Input:
        - X: multi-dimensional array (tensor) of dimension d_1 × … × d_K × N
    Output:
        - Ψ_list: {Ψ_1,…,Ψ_k}
    """
    # initial setup
    dimensions = size(X)
    K = length(dimensions) - 1
    d_list = dimensions[1:K]
    d = prod(d_list)
    N = dimensions[end]

    # initilize the precision matrices if not provided
    if length(Ψ0_list) == 0
        for k = 1:K
            push!(Ψ0_list, 1.0 * I(d_list[k]))
        end
    end

    # set penalty parameters if not provided
    if length(λ_list) == 0
        for k = 1:K
            push!(λ_list, 5 * sqrt(d_list[k] * log(d) / N))
        end
    end

    S_list = [zeros((dk, dk)) for dk in d_list]
    Ψ_list = copy(Ψ0_list)

    for _ = 1:N_iter
        for k = 1:K
            S_list[k] .= Tilde_S_k(k, N, X, Ψ_list)
            Ψ_list[k] .= glasso(S_list[k], λ_list[k])
        end
    end

    return Ψ_list
end

function Tilde_S_k(k::Int, N::Int, X::AbstractArray, Ψ_list::AbstractArray)
    """
    Function to compute:
        (d_k / N*d) * ∑_{i=1}^N V_i^k (V_i^k)^T
    where V_i^k = [X_i × {Ψ_1^(1/2),…,Ψ_(k-1)^(1/2),I,Ψ_(k+1)^(1/2),…,Ψ_K^(1/2)}]_(k);
    here × is the tensor product and []_(k) the mode-k matricication.
    """
    Ψ_list = copy(Ψ_list)
    d_list = size.(Ψ_list, 1)
    dk = d_list[k]
    d = prod(d_list)
    K = length(d_list)
    S = zeros((dk, dk))
    X_mode_k = zeros(dk, Int(d / dk))
    X_mode_k_trans = zeros(Int(d / dk), dk)

    deleteat!(Ψ_list, k)
    reverse!(Ψ_list)
    KP_Ψ_list = (length(Ψ_list) == 1) ? Ψ_list[1] : kron(Ψ_list...) # handles case K=2

    @inbounds for i = 1:N
        ## TODO: modify to work with various K ##
        # str_to_eval = string("X[", join([":," for _ = 1:K]), "i]") # adapt to different K
        # X_mode_k .= tenmat(eval(Meta.parse(str_to_eval)), k)
        X_mode_k .= tenmat(X[:, :, i], k)
        mul!(X_mode_k_trans, KP_Ψ_list, copy(X_mode_k'))
        mul!(S, X_mode_k, X_mode_k_trans, dk / (d * N), 1.0)
    end

    return S
end

# # Testing
# include("utils.jl")
# include("utils_sim.jl")
# include("glasso.jl")

# X = gen_kronecker_data("kp", "sb", 2, 1000, [5,10,15])
# Ψ_hat_list = kglasso(X)


text:
using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
      Test, BandedMatrices, FillArrays

function second_derivative_stencil(N)
  A = zeros(N,N+2)
  A[1,1:5] = [0.916667,   -1.66667,    0.5,   0.333333,  -0.0833333]
  A[2,1:5] = [-0.0833333,   1.33333,   -2.5,   1.33333,   -0.0833333]
  A[3,1:5] = [-0.0833333,   0.333333,   0.5,  -1.66667,    0.916667]
  A[4,1:5] = [0.916667,   -4.66667,    9.5,  -8.66667,    2.91667]
  for i in 5:N-4, j in 5:N-2
      (j-i==0 || j-i==2) && (A[i,j]=1)
      j-i==1 && (A[i,j]=-2)
  end
  A[end,end-4:end] = reverse(A[1,1:5])
  A[end-1,end-4:end] = reverse(A[2,1:5])
  A[end-2,end-4:end] = reverse(A[3,1:5])
  A[end-3,end-4:end] = reverse(A[4,1:5])

  return A
end

# Analytic solutions to higher order operators.
# Do not modify unless you are completely certain of the changes.
function fourth_deriv_approx_stencil(N)
    A = zeros(N,N+2)
    A[1,1:8] = [3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3]
    A[2,1:8] = [2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6]
    A[N-1,N-5:end] = reverse([2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6], dims=2)
    A[N,N-5:end] = reverse([3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3], dims=2)
    for i in 3:N-2
        A[i,i-2:i+4] = [-1/6 2.0 -13/2 28/3 -13/2 2.0 -1/6]
    end
    return A
end

function second_deriv_fourth_approx_stencil(N)
    A = zeros(N,N+2)
    A[1,1:6] = [5/6 -1.25 -1/3 7/6 -0.5 5/60]
    A[N,N-3:end] = reverse([5/6 -1.25 -1/3 7/6 -0.5 5/60], dims=2)
    for i in 2:N-1
        A[i,i-1:i+3] = [-1/12 4/3 -5/2 4/3 -1/12]
    end
    return A
end

function convert_by_multiplication(::Type{Array}, A::AbstractDerivativeOperator{T}, N::Int=A.dimension) where T
    @assert N >= A.stencil_length # stencil must be able to fit in the matrix
    mat = zeros(T, (N, N+2))
    v = zeros(T, N+2)
    for i=1:N+2
        v[i] = one(T)
        #=
            calculating the effect on a unit vector to get the matrix of transformation
            to get the vector in the new vector space.
        =#
        mul!(view(mat,:,i), A, v)
        v[i] = zero(T)
    end
    return mat
end

# Tests the corrrectness of stencils, along with concretization.
# Do not modify the following test-set unless you are completely certain of your changes.
@testset "Correctness of Stencils" begin
    N = 20
    L1 = UpwindDifference(1,2, 1.0, N, t->1.0)
    correct = [-3/2, 2.0, -1/2]
    @test L1.stencil_coefs ≈ correct

    L1 = UpwindDifference(1,3, 1.0, N, t->1.0)
    correct = [-2/6, -3/6, 6/6, -1/6]
    @test L1.stencil_coefs ≈ correct
end

@testset "Taking derivatives" begin
    N = 20
    x = 0:1/(N-1):1
    y = 2x.^2 .- 3x .+2
    y_ = y[2:end-1]
    dx = x[2] - x[1]
    # y = x.^2

    # Dirichlet BC with fixed end points
    Q = RobinBC((1.0, 0.0, y[1]), (1.0, 0.0, y[end]), 1.0)
    U = UpwindDifference(1,2, dx, N-2, t->1.0)
    A = CenteredDifference(1,2, dx, N-2)
    D1 = CenteredDifference(1,2, dx, N-4) # For testing whether the array is constant

    res1 = U*Q*y_
    res2 = A*Q*y_

    @test res1[3:end-2] ≈ res2[1:end-4] # shifted due to upwind operators
    # It is shifted by a constant value so its first derivative has to be 0
    @test D1*(res1[3:end-2] - res2[3:end-2]) ≈ zeros(12) atol=10.0^(-6)

    # @test res3[3:end-2] ≈ res4[1:end-4] # shifted due to upwind operators
    # It is shifted by a constant value so its first derivative has to be 0
    # @test D1*(res4[3:end-2] - res4[3:end-2]) ≈ zeros(12) atol=10.0^(-6)

    y = 3x.^3 .- 4x.^2 .+ 2x .+ 1
    y_ = y[2:end-1]
    Q = RobinBC((1.0, 0.0, y[1]), (1.0, 0.0, y[end]), 1.0)
    U = UpwindDifference(2,2, 1.0, N-2, t->1.0)
    A = CenteredDifference(2,2, 1.0, N-2)
    res1 = U*Q*y_
    res2 = A*Q*y_

    @test res1 ≈ res2 # shifted due to upwind operators
end

@testset "Full and Sparse functions:" begin
    N = 10
    d_order = 2
    approx_order = 2
    correct = second_derivative_stencil(N)
    A = UpwindDifference(d_order,approx_order,1.0,N,t->1.0)

    @test convert_by_multiplication(Array,A,N) ≈ correct atol=10.0^(-4)
    @test Array(A) ≈ second_derivative_stencil(N) atol=10.0^(-4)
    @test sparse(A) ≈ second_derivative_stencil(N) atol=10.0^(-4)
    @test BandedMatrix(A) ≈ second_derivative_stencil(N) atol=10.0^(-4)
    @test_broken opnorm(A, Inf) ≈ opnorm(correct, Inf) atol=10.0^(-4)

    # testing higher derivative and approximation concretization
    N = 20
    d_order = 4
    approx_order = 4
    A = UpwindDifference(d_order,approx_order,1.0,N,t->1.0)
    correct = convert_by_multiplication(Array,A,N)

    @test Array(A) ≈ correct
    @test sparse(A) ≈ correct
    @test BandedMatrix(A) ≈ correct

    N = 100
    d_order = 8
    approx_order = 8
    A = UpwindDifference(d_order,approx_order,1.0,N,t->1.0)
    correct = convert_by_multiplication(Array,A,N)

    @test Array(A) ≈ correct
    @test sparse(A) ≈ correct
    @test BandedMatrix(A) ≈ correct

    # testing correctness of multiplication
    N = 1000
    d_order = 4
    approx_order = 10
    y = collect(1:1.0:N+2).^4 - 2*collect(1:1.0:N+2).^3 + collect(1:1.0:N+2).^2;
    y = convert(Array{BigFloat, 1}, y)

    A = UpwindDifference(d_order,approx_order,one(BigFloat),N,t->1.0)
    correct = convert_by_multiplication(Array,A,N)
    @test Array(A) ≈ correct
    @test sparse(correct) ≈ correct
    @test BandedMatrix(A) ≈ correct
    @test A*y ≈ Array(A)*y
end

@testset "Indexing tests" begin
    N = 1000
    d_order = 4
    approx_order = 10

    A = UpwindDifference(d_order,approx_order,1.0,N,t->1.0)
    @test A[1,1] == Array(A)[1,1]
    @test A[10,20] == 0.0

    correct = Array(A)
    for i in 1:N
        @test A[i,i] ≈ correct[i,i] atol=10^-9
    end

    # Indexing Tests
    N = 1000
    d_order = 2
    approx_order = 2

    A = UpwindDifference(d_order,approx_order,1.0,N,t->1.0)
    M = Array(A,N)
    @test A[1,1] == M[1,1]
    @test A[1:4,1] == M[1:4,1]
    @test A[5,2:10] == M[5,2:10]
    @test A[60:100,500:600] == M[60:100,500:600]

    d_order = 4
    approx_order = 10

    A = UpwindDifference(d_order,approx_order,1.0,N,t->1.0)
    M = Array(A,N)
    @test A[1,1] == M[1,1]
    @test A[1:4,1] == M[1:4,1]
    @test A[5,2:10] == M[5,2:10]
    @test A[524,:] == M[524,:]
    @test A[60:100,500:600] == M[60:100,500:600]
end

@testset begin "Operations on matrices"
    N = 51
    M = 101
    d_order = 2
    approx_order = 2

    xarr = range(0,stop=1,length=N)
    yarr = range(0,stop=1,length=M)
    dx = xarr[2]-xarr[1]
    dy = yarr[2]-yarr[1]
    F = [x^2+y for x = xarr, y = yarr]

    A = UpwindDifference(d_order,approx_order,dx,length(xarr)-2,t->1.0)
    B = UpwindDifference(d_order,approx_order,dy,length(yarr),t->1.0)

    @test A*F ≈ 2*ones(N-2,M)
    F*B
    # A*F*B

    G = [x^2+y^2 for x = xarr, y = yarr]

    @test A*G ≈ 2*ones(N-2,M) atol=1e-2
    G*B
    # A*G*B
end


text:
using ChainRulesCore
using LinearAlgebra
using KrylovKit

@doc raw"
    num_grad(f, K::Real; [δ = 1e-5])

return the numerical gradient of `f` at `K` calculated with
`(f(K+δ/2) - f(K-δ/2))/δ`

# example

```jldoctest; setup = :(using TensorNetworkAD)
julia> TensorNetworkAD.num_grad(x -> x * x, 3) ≈ 6
true
```
"
function num_grad(f, K; δ::Real=1e-5)
    if eltype(K) == ComplexF64
        (f(K + δ / 2) - f(K - δ / 2)) / δ +
            (f(K + δ / 2 * 1.0im) - f(K - δ / 2 * 1.0im)) / δ * 1.0im
    else
        (f(K + δ / 2) - f(K - δ / 2)) / δ
    end
end

@doc raw"
    num_grad(f, K::AbstractArray; [δ = 1e-5])

return the numerical gradient of `f` for each element of `K`.

# example

```jldoctest; setup = :(using TensorNetworkAD, LinearAlgebra)
julia> TensorNetworkAD.num_grad(tr, rand(2,2)) ≈ I
true
```
"
function num_grad(f, a::AbstractArray; δ::Real=1e-5)
    b = Array(copy(a))
    df = map(CartesianIndices(b)) do i
        foo = x -> (ac = copy(b); ac[i] = x; f(_arraytype(a)(ac)))
        num_grad(foo, b[i], δ=δ)
    end
    return _arraytype(a)(df)
    # map(CartesianIndices(a)) do i
    #     foo = x -> (ac = copy(a); ac[i] = x; f(ac))
    #     num_grad(foo, a[i], δ=δ)
    # end
end

# patch since it's currently broken otherwise
function ChainRulesCore.rrule(::typeof(Base.typed_hvcat), ::Type{T}, rows::Tuple{Vararg{Int}}, xs::S...) where {T,S}
    y = Base.typed_hvcat(T, rows, xs...)
    function back(ȳ)
        return NoTangent(), NoTangent(), NoTangent(), permutedims(ȳ)...
    end
    return y, back
end

# improves performance compared to default implementation, also avoids errors
# with some complex arrays
function ChainRulesCore.rrule(::typeof(LinearAlgebra.norm), A::AbstractArray)
    n = norm(A)
    function back(Δ)
        return NoTangent(), Δ .* A ./ (n + eps(0f0)), NoTangent()
    end
    return n, back
end

"""
    ChainRulesCore.rrule(::typeof(leftenv), Au::AbstractArray{T}, Ad::AbstractArray{T}, M::AbstractArray{T}, FL::AbstractArray{T}; kwargs...) where {T}

```
           ┌──   Au  ──┐
           │     │     │
dM    = - FL ──    ──  ξl
           │     │     │
           └──   Ad  ──┘

           ┌──       ──┐
           │     │     │
dAu  = -  FL ──  M ──  ξl
           │     │     │
           └──   Ad  ──┘

           ┌──   Au  ──┐
           │     │     │
dAd  = -  FL  ── M ──  ξl
           │     │     │
           └──       ──┘
```
"""

function ChainRulesCore.rrule(::typeof(leftenv), Au::AbstractArray{T}, Ad::AbstractArray{T}, M::AbstractArray{T}, FL::AbstractArray{T}; kwargs...) where {T}
    λl, FL = leftenv(Au, Ad, M, FL)
    # @show λl
    function back((dλ, dFL))
        dFL -= Array(ein"abc,abc ->"(conj(FL), dFL))[] * FL
        ξl, info = linsolve(FR -> ein"((abc,ceh),dgeb),fgh -> adf"(Au, FR, M, Ad), conj(dFL), -λl, 1; maxiter = 1)
        @assert info.converged==1
        # errL = ein"abc,cba ->"(FL, ξl)[]
        # abs(errL) > 1e-1 && throw("FL and ξl aren't orthometric. err = $(errL)")
        dAu = -ein"((adf,fgh),dgeb),ceh -> abc"(FL, Ad, M, ξl)
        dAd = -ein"((adf,abc),dgeb),ceh -> fgh"(FL, Au, M, ξl)
        dM = -ein"(adf,abc),(fgh,ceh) -> dgeb"(FL, Au, Ad, ξl)
        return NoTangent(), conj(dAu), conj(dAd), conj(dM), NoTangent()...
    end
    return (λl, FL), back
end

"""
    function ChainRulesCore.rrule(::typeof(norm_FL), Au::AbstractArray{T}, Ad::AbstractArray{T}, FL::AbstractArray{T}; kwargs...) where {T}

```
           ┌──       ──┐
           │     │     │
dAu  = -  FL ────┼──── ξl
           │     │     │
           └──   Ad  ──┘

           ┌──   Au  ──┐
           │     │     │
dAd  = -  FL  ───┼───  ξl
           │     │     │
           └──       ──┘
```
"""

function ChainRulesCore.rrule(::typeof(norm_FL), Au::AbstractArray{T}, Ad::AbstractArray{T}, FL::AbstractArray{T}; kwargs...) where {T}
    λl, FL = norm_FL(Au, Ad, FL)
    function back((dλ, dFL))
        dFL -= Array(ein"ab,ab ->"(conj(FL), dFL))[] * FL
        ξl, info = linsolve(FR -> ein"(be,acb), dce -> ad"(FR,Au,Ad), conj(dFL), -λl, 1; maxiter = 1)
        # @assert info.converged==1
        # errL = ein"abc,cba ->"(FL, ξl)[]
        # abs(errL) > 1e-1 && throw("FL and ξl aren't orthometric. err = $(errL)")
        dAu = -ein"(ad,dce), be -> acb"(FL, Ad, ξl)
        dAd = -ein"(ad,acb), be -> dce"(FL, Au, ξl)
        return NoTangent(), conj(dAu), conj(dAd), NoTangent()...
    end
    return (λl, FL), back
end


text:
using Symmetry
using Test

rt2 = sqrt(2.0)
rt3 = sqrt(3.0)
twocos72 = 2*cosd(72.0)
twocos144 = 2*cosd(144.0)

include("pgs/Cn.jl")
include("pgs/Cnv.jl")
include("pgs/Cnh.jl")
include("pgs/Sn.jl")
include("pgs/Dn.jl")
include("pgs/Dnd.jl")
include("pgs/Dnh.jl")
allofem = ["C1", "C2",  "C3",  "C4",  "C5",  "C6",
                 "C2v", "C3v", "C4v", "C5v", "C6v",
                 "C2h", "C3h", "C4h", "C5h", "C6h",
                 "S4",  "S6",  "S8",  "S10",
                 "D2",  "D3",  "D4",  "D5",  "D6",
                 "D2d", "D3d", "D4d", "D5d", "D6d",
                 "D2h", "D3h", "D4h", "D5h", "D6h"]

@testset "Symmetry" begin
    include("test_symel_gen.jl")
    include("test_ctab_gen.jl")
end


text:
using Rankers.Storages
using LOLTools.SummonerV4: SummonerV4, http_action
using LOLTools.SpectatorV4
using LOLToolsDataDragon.DataDragon.en_US

api_key = get(ENV, "RIOT_TOKEN", "")
region = "kr"
summonerName = "DRX Deft"
read_from_file = true # false
summonerdto = SummonerV4.summoner_by_name(api_key, region, summonerName; action=store_with(region, http_action, debug=true))
encryptedSummonerId = summonerdto.id
current_game = SpectatorV4.active_games(api_key, region, encryptedSummonerId; action=update_with(region, http_action, read_from_file=read_from_file, debug=true))
for participant in current_game.participants
    print("*", participant.summonerName)
    print(" ---- ")
    println(en_US.SummonerSpells[participant.spell1Id].name, ' ', en_US.SummonerSpells[participant.spell2Id].name)
end


text:
using Test, Logging

const curdir = @__DIR__
const scripts_dir = normpath(joinpath(curdir, "..", "_literate"))

function finish(bko, bke, wro, rdo, wre, rde)
    redirect_stdout(bko)
    redirect_stderr(bke)
    close(wro)
    close(rdo)
    close(wre)
    close(rde)
end

function try_run(tf)
    bko = stdout
    bke = stderr
    (rdo, wro) = redirect_stdout()
    (rde, wre) = redirect_stderr()
    @sync begin
        Logging.disable_logging(Logging.LogLevel(3_000))
        try
            include(tf)
        catch e
            Logging.disable_logging(Logging.Debug)
            finish(bko, bke, wro, rdo, wre, rde)
            throw(e)
        end
        Logging.disable_logging(Logging.Debug)
    end
    finish(bko, bke, wro, rdo, wre, rde)
end

"""
    strip_code

Gets the julia code out, unless if the block starts with `# notest` in which
case it's ignored (useful for plots and ScikitLearn things which will fail on
travis).
"""
function strip_code(fpath)
    lines = readlines(fpath)
    jc = IOBuffer()
    record = true
    is_ci = get(ENV, "CI", "false") == "true"
    for line in lines
        if is_ci && startswith(line, "# notest")
            record = false
            continue
        end
        if record
            # is it code?
            startswith(line, "# ") && continue
            # is it a savefig ?
            startswith(line, "savefig") && continue
            # it is
            write(jc, line * "\n")
        else
            # toggle it back on if new block
            if startswith(line, "# ")
                record = true
            end
        end
    end
    tp = tempname()
    write(tp, take!(jc))
    return tp
end

for (root, _, files) in walkdir(scripts_dir), file in files

    # NOTE: if want to run a single file  in isolation, uncomment line below

#    splitdir(file)[2] ∉ ("ISL-lab-8.jl",) && continue
#    startswith(splitdir(file)[2], "ISL-lab-4") || continue

    file == ".DS_Store" && continue

    @testset "testing $file" begin
        println("\n\n>> looking at $file ...")
        path = joinpath(root, file)
        tf = tempname()
        write(tf, """
        module Tester
            import ..strip_code
            include(strip_code("$path"))
        end
        """)
        @test begin
            try_run(tf); true
        end
    end
end
println("All done 🍻")

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
using MNIST
sigmoid(z) = 1 ./ (1 + exp(-z))

# load MNIST data
println("Loading MNIST dataset...")
X_train, y_train = traindata()
X_test, y_test = testdata()
y_train, y_test = y_train', y_test'

# mean normalization
X_train /= 255.0
X_test /= 255.0

# make data more convenient for further operations
y_train = y_train .+ 1
y_test = y_test .+ 1
y_train = convert(Array{Int64}, y_train)
y_test = convert(Array{Int64}, y_test)

n, m = size(X_train)
K = 10
h = 256
W1 = 0.01 * rand(h, n)
b1 = zeros(h, 1)
W2 = 0.01 * rand(K, h)
b2 = zeros(K,1)

alpha = 0.2
reg = 0.0

X = X_train
y = y_train

for i = 1:20
    # forward prop 1 hidden layer
    z2 = W1 * X .+ b1 # h x m
    a2 = sigmoid(z2) # h x m
    z3 = W2 * a2 .+ b2 # K x m

    # softmax forw prop
    mat = exp(z3)
    probs = mat ./ sum(mat, 1)

    I = sub2ind((K, m), collect(y), 1:m)
    data_loss = -sum(log(probs[I]))/m
    reg_loss = 0.5 * reg * (sum(W1 .* W1) + sum(W2 .* W2))
    loss = data_loss + reg_loss

    println(loss)

    er3 = probs # K x m
    er3[I] -= 1
    er2 = (W2' * er3) .* (a2 .* (1-a2)) # h x m

    dW2 = er3 * a2' # (K x m) * (m x h) = (K x h)
    db2 = sum(er3, 2)
    dW1 = er2 * X' # (h x m) * (m x n) = (h x n)
    db1 = sum(er2, 2)

    W1 -= alpha * ((dW1 / m) + reg * W1)
    b1 -= alpha * (db1 / m)
    W2 -= alpha * ((dW2 / m) + reg * W2)
    b2 -= alpha * (db2 / m)
end

max_stars_repo_name: ilkerkesen/ufldl-tutorial

max_issues_repo_name: ilkerkesen/ufldl-tutorial

max_forks_repo_name: ilkerkesen/ufldl-tutorial


text:
######### StanSample Bernoulli example  ###########

using AxisKeys
using StanSample

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

sm = SampleModel("bernoulli", bernoulli_model);
rc = stan_sample(sm; data=bernoulli_data,
  num_threads=6, num_cpp_chains=1, num_chains=6, seed=123);

if success(rc)
  chns = read_samples(sm, :keyedarray)
  chns |> display
end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
# Never return Int and also never return abstract types.
# E.g. eltype((1.,1)) == Real
"""
    suff(x) -> Type
If `x` stores BigFloats, `suff` returns BigFloat, else `suff` returns `Float64`.
"""
suff(x::T) where T<:AbstractFloat = T
suff(x::Real) = Float64
suff(x::Complex) = suff(real(x))
suff(x::AbstractArray) = suff(x[1])
suff(x::Tuple) = suff(promote(x...)[1])
suff(x::Union{Nothing, Missing}) = Float64
suff(x) = eltype(x)
suff(args...) = try suff(promote(args...)[1]) catch;  suff(args[1]) end

suff(x::Num) = Num
# Allow for differentiation through suff arrays.
suff(x::T) where T<:ForwardDiff.Dual = T

"""
    MaximalNumberOfArguments(F::Function) -> Int
Infers argument structure of given function, i.e. whether it is of the form `F(x)` or `F(x,y)` or `F(x,y,z)` etc. and returns maximal number of accepted arguments of all overloads of `F` as integer.
"""
MaximalNumberOfArguments(F::Function) = maximum([length(Base.unwrap_unionall(m.sig).parameters)-1 for m in methods(F)])

"""
    GetArgLength(F::Function; max::Int=100) -> Int
Attempts to determine input structure of `F`, i.e. whether it accepts `Number`s or `AbstractVector`s and of what length.
This is achieved by successively evaluating the function on `rand(i)` until the evaluation no longer throws errors.
As a result, `GetArgLength` will be unable to determine the correct input structure if `F` errors on `rand(i)`.
!!! note
    Does NOT discriminate between `Real` and `Vector{Real}` of length one, i.e. `Real`↦`+1`.
    To disciminate between these two options, use `_GetArgLength()` instead.
"""
GetArgLength(F::Function; max::Int=100) = _GetArgLength(F; max=max) |> abs

"""
    _GetArgLength(F::Function; max::Int=100) -> Int
!!! note
    Discriminates between `Real` and `Vector{Real}` of length one, i.e.:
    `Real`↦`-1` and `x::AbstractVector{<:Real}`↦`length(x)`.
"""
function _GetArgLength(F::Function; max::Int=100)
    num = MaximalNumberOfArguments(F)
    if num == 1
         _GetArgLengthOutOfPlace(F; max=max)
    elseif num == 2
        _GetArgLengthInPlace(F; max=max)[2]
    else
        throw("Given function $F appears to take $num number of arguments. It should either accept only 1 argument or 2 for functions which mutate their first argument.")
    end
end
function _GetArgLengthOutOfPlace(F::Function; max::Int=100)
    @assert max > 1
    try     F(rand());  return -1    catch; end
    for i in 1:(max+1)
        try
            res = F(rand(i))
            isnothing(res) && throw("Function returned Nothing for i=$i.")
        catch y
            (isa(y, BoundsError) || isa(y, MethodError) || isa(y, DimensionMismatch) || isa(y, ArgumentError) || isa(y, AssertionError)) && continue
            @warn "Encountered apparent error in specification of function for i=$i: $y."
            continue
        end
        i < (max + 1) ? (return i) : throw(ArgumentError("Function input appears to have >$max components, aborting. Either increase keyword max or ensure function errors on rand(i) for i larger than true component length."))
    end
end
function _GetArgLengthInPlace(F::Function; max::Int=100)
    @assert max > 1
    res = 0.;     Res = zeros(max);     RES = zeros(max,max);   RESS = zeros(max,max,max)
    function _TryOn(output, input)
        try
            F(output, input)
        catch y
            if !(isa(y, BoundsError) || isa(y, MethodError) || isa(y, DimensionMismatch) || isa(y, ArgumentError) || isa(y, AssertionError))
                @warn "Encountered apparent error in specification of function."
                rethrow()
            end
            return nothing
        end;    (FindSubHyperCube(output), (input isa Number ? -1 : length(input)))
    end
    function TryAll(input)
        for AnOut in (res, Res, RES, RESS)
            Z = _TryOn(AnOut, input)
            !isnothing(Z) && return Z
        end
        nothing
    end
    X = TryAll(rand());    !isnothing(X) && return X
    for i in 1:max
        X = TryAll(rand(i))
        !isnothing(X) && return X
    end
    throw(ArgumentError("Function input appears to have >$max components, aborting. Either increase keyword max or ensure function errors on rand(i) for i larger than true component length."))
end

FindSubHyperCube(N::Number, tester::Function=(x->x!=0.0)) = -1
FindSubHyperCube(Vec::AbstractVector{<:Number}, tester::Function=(x->x!=0.0)) = findlast(tester, Vec)

# Still highly inefficient but apparently works without bugs
function FindSubHyperCube(Mat::AbstractArray{<:Number, dim}, tester::Function=(x->x!=0.0)) where dim
    Size = [size(Mat)...];    Biggest = copy(Size)
    for ind in 1:dim
        for val in 1:Size[ind]
            if all(!tester, view(Mat, [Colon() for j in 1:ind-1]..., val, [Colon() for j in ind+1:dim]...))
                if all(!tester, view(Mat, [1:Biggest[i] for i in 1:ind-1]..., val+1:Size[ind], [Colon() for j in ind+1:dim]...))
                    Biggest[ind] = val-1;   break
                end
            end
        end
    end;    Tuple(Biggest)
end

"""
    KillAfter(F::Function, args...; timeout::Real=5, verbose::Bool=false, kwargs...)
Tries to evaluate a given function `F` before a set `timeout` limit is reached and interrupts the evaluation and returns `nothing` if necessary.
NOTE: The given function is evaluated via F(args...; kwargs...).
"""
function KillAfter(F::Function, args...; timeout::Real=5, verbose::Bool=false, kwargs...)
    Res = nothing
    G() = try F(args...; kwargs...) catch Err
        if verbose
            if Err isa DivideError
                @warn "KillAfter: Could not evaluate given Function $(nameof(F)) before timeout limit of $timeout seconds was reached."
            else
                @warn "KillAfter: Could not evaluate given Function $(nameof(F)) because error was thrown: $Err."
            end
        end
    end
    task = @async(G())
    if timedwait(()->istaskdone(task), timeout) == :timed_out
        @async(Base.throwto(task, DivideError())) # kill task
    else
        Res = fetch(task)
    end;    Res
end

"""
    Builder(Fexpr::Union{<:AbstractVector{<:Num},<:Num}, args...; inplace::Bool=false, parallel::Bool=false, kwargs...)
Builds `RuntimeGeneratedFunctions` from expressions via build_function().
"""
function Builder(Fexpr::Union{<:AbstractArray{<:Num},<:Num}, args...; inplace::Bool=false, parallel::Bool=false, kwargs...)
    parallelization = parallel ? Symbolics.MultithreadedForm() : Symbolics.SerialForm()
    Res = if (Fexpr isa Num && args[1] isa Num)
        # build_function throws error when using parallel keyword for R⟶R functions
        Symbolics.build_function(Fexpr, args...; expression=Val{false}, kwargs...)
    else
        Symbolics.build_function(Fexpr, args...; expression=Val{false}, parallel=parallelization, kwargs...)
    end
    try
        Res[inplace ? 2 : 1]
    catch;
        Res
    end
end


text:
module Remarkable
using HTTP
using JSON
using ZipFile
using Dates
using Crayons
using Parameters: @with_kw, type2dict
using UUIDs: uuid4
using AbstractTrees: AbstractTrees, print_tree
export RemarkableClient,
    list_items, get_item, delete_item!, create_folder!, download_pdf, upload_pdf!
export print_tree

include("const.jl")
include("client.jl")
include("documents.jl")
include("api.jl")

end


text:
using SequentialLOB
using StylizedFacts
using Plots

# Configuration Arguments
num_paths = 1
M = 400
T = 2299
p₀ = 238.75
L = 200
Δx = L / M

# Free-Parameters
λ = 1.0
D = 1.0
Δt = (Δx^2) / (2.0 * D)
σ = 1.5
μ = 1.0
nu = 0.5
α_slob = 100.0
α_lob = 0.0
println(Δx)
println(Δt)
# LOB Model

lob_model = SLOB(num_paths,
T, p₀, M, L, D, σ, nu, α_lob, SourceTerm(λ, μ))

sample_price_paths = lob_model(7136)
plt = plot(
    1:T + 1, sample_price_paths[:,1], legend=false,
    xlab="Time", ylab="Mid-Price",
)
for i in 2:num_paths
    plot!(plt, 1:T + 1, sample_price_paths[:,i], legend=false)
end
display(plt)
savefig(plt, "$(dirname(@__FILE__))/lob_model_init_pp.pdf")

# SLOB Model

slob_model = SLOB(num_paths,
T, p₀, M, L, D, σ, nu, α_slob, SourceTerm(λ, μ))

sample_price_paths = slob_model(19822)
plt = plot(
    1:T + 1, sample_price_paths[:,1], legend=false,
    xlab="Time", ylab="Mid-Price",
)
for i in 2:num_paths
    plot!(plt, 1:T + 1, sample_price_paths[:,i], legend=false)
end
display(plt)
savefig(plt, "$(dirname(@__FILE__))/slob_model_init_pp.pdf")

# Compare LOB to SLOB

lob_mid_prices = lob_model(7136)
lob_stylized_facts = StylizedFactsPlot(lob_mid_prices[:,1])
plt = StylizedFacts.plot_log_returns(lob_stylized_facts, "")
display(plt)
savefig(plt, "$(dirname(@__FILE__))/lob_model_init_log_rets.pdf")

slob_mid_prices = slob_model(19822)
slob_stylized_facts = StylizedFactsPlot(slob_mid_prices[:,1])
plt = StylizedFacts.plot_log_returns(slob_stylized_facts, "")
display(plt)
savefig(plt, "$(dirname(@__FILE__))/slob_model_init_log_rets.pdf")

max_stars_repo_name: GantZA/CALOBMTMD

max_issues_repo_name: GantZA/CALOBMTMD

max_forks_repo_name: GantZA/CALOBMTMD


text:
function AggregateBoundedPropagator(arg0::Collection)
    return AggregateBoundedPropagator((Collection,), arg0)
end

function get_initial_state(obj::AggregateBoundedPropagator)
    return jcall(obj, "getInitialState", SpacecraftState, ())
end

function get_max_date(obj::AggregateBoundedPropagator)
    return jcall(obj, "getMaxDate", AbsoluteDate, ())
end

function get_min_date(obj::AggregateBoundedPropagator)
    return jcall(obj, "getMinDate", AbsoluteDate, ())
end

function get_pv_coordinates(obj::AggregateBoundedPropagator, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function reset_initial_state(obj::AggregateBoundedPropagator, arg0::SpacecraftState)
    return jcall(obj, "resetInitialState", void, (SpacecraftState,), arg0)
end


text:
@enum(WindowStyle,
window_none = 0,
window_titlebar = 1 << 0,
window_resize = 1 << 1,
window_close = 1 << 2,
window_fullscreen = 1 << 3,
window_defaultstyle = (1 << 0) | (1 << 1) | (1 << 2))

type ContextSettings
    depth_bits::UInt32
    stencil_bits::UInt32
    antialiasing_level::UInt32
    major_version::UInt32
    minor_version::UInt32
    attribute_flags::UInt32

    function ContextSettings(depth_bits::Integer, stencil_bits::Integer, antialiasing_level::Integer, major_version::Integer, minor_version::Integer)
        new(depth_bits, stencil_bits, antialiasing_level, major_version, minor_version, 0)
    end
    function ContextSettings()
        new(0, 0, 0, 2, 2, 0)
    end
end

function RenderWindow(mode::VideoMode, title::AbstractString, settings::ContextSettings, style::WindowStyle...)
    style_int = 0
    for i = 1:length(style)
        style_int |= Int(style[i])
    end
    settings_ptr = Ref(settings)
    icon = Image(joinpath(dirname(@__FILE__),"..","..","..","assets","sfmljl_icon.png"))
    window = RenderWindow(ccall((:sfRenderWindow_create, libcsfml_graphics), Ptr{Void}, (VideoMode, Ptr{Cchar}, UInt32, Ref{ContextSettings},), mode, title, style_int, settings_ptr))
    # set_icon(window, icon, 64, 64)
    return window
end

function RenderWindow(mode::VideoMode, title::AbstractString, style::WindowStyle...=window_defaultstyle)
    return RenderWindow(mode, title, ContextSettings(), style[1])
end

function RenderWindow(title::AbstractString, width::Integer, height::Integer)
    mode = VideoMode(width, height)
    return RenderWindow(mode, title, window_defaultstyle)
end

function set_framerate_limit(window::RenderWindow, limit::Integer)
    ccall((:sfRenderWindow_setFramerateLimit, libcsfml_graphics), Void, (Ptr{Void}, UInt,), window.ptr, limit)
end

function set_vsync_enabled(window::RenderWindow, enabled::Bool)
    ccall((:sfRenderWindow_setVerticalSyncEnabled, libcsfml_graphics), Void, (Ptr{Void}, Int32,), window.ptr, enabled)
end

function isopen(window::RenderWindow)
    return Bool(ccall((:sfRenderWindow_isOpen, libcsfml_graphics), Int32, (Ptr{Void},), window.ptr))
end

function pollevent(window::RenderWindow, event::Event)
    return Bool(ccall((:sfRenderWindow_pollEvent, libcsfml_graphics), Int32, (Ptr{Void}, Ptr{Void},), window.ptr, event.ptr))
end

function waitevent(window::RenderWindow, event::Event)
    return Bool(ccall((:sfRenderWindow_waitEvent, libcsfml_graphics), Int32, (Ptr{Void}, Ptr{Void},), window.ptr, event.ptr))
end

function set_position(window::RenderWindow, pos::Vector2)
    pos = to_vec2i(pos)
    ccall((:sfRenderWindow_setPosition, libcsfml_graphics), Void, (Ptr{Void}, Vector2i,), window.ptr, pos)
end

function get_position(window::RenderWindow)
    return ccall((:sfRenderWindow_getPosition, libcsfml_graphics), Vector2i, (Ptr{Void},), window.ptr)
end

function set_size(window::RenderWindow, window_size::Vector2)
    window_size = to_vec2u(window_size)
    ccall((:sfRenderWindow_setSize, libcsfml_graphics), Void, (Ptr{Void}, Vector2u), window.ptr, window_size)
end

function get_size(window::RenderWindow)
    return to_vec2i(ccall((:sfRenderWindow_getSize, libcsfml_graphics), Vector2u, (Ptr{Void},), window.ptr))
end

function set_title(window::RenderWindow, title::AbstractString)
    ccall((:sfRenderWindow_setTitle, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar},), window.ptr, title)
end

function set_icon(window::RenderWindow, pixels::Array{UInt8}, width=64, height=64)
    ccall((:sfRenderWindow_setIcon, libcsfml_graphics), Void, (Ptr{Void}, UInt32, UInt32, Ptr{UInt8},), window.ptr, width, height, pointer(pixels))
end

function set_icon(window::RenderWindow, image::Image, width=64, height=64)
    set_icon(window, get_pixels(image), width, height)
end

function set_visible(window::RenderWindow, visible::Bool)
    ccall((:sfRenderWindow_setVisible, libcsfml_graphics), Void, (Ptr{Void}, Int32,), window.ptr, visible)
end

function set_mousecursor_visible(window::RenderWindow, visible::Bool)
    ccall((:sfRenderWindow_setMouseCursorVisible, libcsfml_graphics), Void, (Ptr{Void}, Int32,), window.ptr, visible)
end

function set_keyrepeat_enabled(window::RenderWindow, enabled::Bool)
    ccall((:sfRenderWindow_setKeyRepeatEnabled, libcsfml_graphics), Void, (Ptr{Void}, Int32,), window.ptr, enabled)
end

function set_active(window::RenderWindow, active::Bool)
    ccall((:sfRenderWindow_setActive, libcsfml_graphics), Void, (Ptr{Void}, Int32,), window.ptr, active)
end

function requestfocus(window::RenderWindow)
    ccall((:sfRenderWindow_requestFocus, libcsfml_graphics), Void, (Ptr{Void},), window.ptr)
end

function hasfocus(window::RenderWindow)
    return Bool(ccall((:sfRenderWindow_hasFocus, libcsfml_graphics), Int32, (Ptr{Void},), window.ptr))
end

function set_view(window::RenderWindow, view::View)
    ccall((:sfRenderWindow_setView, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void},), window.ptr, view.ptr)
    window._view = view
end

function get_view(window::RenderWindow)
    return View(ccall((:sfRenderWindow_getView, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), window.ptr))
end

function get_default_view(window::RenderWindow)
    return View(ccall((:sfRenderWindow_getDefaultView, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), window.ptr))
end

function get_viewport(window::RenderWindow)
    return ccall((:sfRenderWindow_getViewport, libcsfml_graphics), IntRect, (Ptr{Void},), window.ptr)
end

function draw(window::RenderWindow, object::Drawable, shader::Shader)
    state = RenderStates(shader)
    draw(window, object, state)
end

function draw(window::RenderWindow, object::Drawable, state::RenderStates)
    draw(window, object, state.ptr)
end

function draw(window::RenderWindow, object::Drawable)
    draw(window, object, C_NULL)
end

function draw(window::RenderWindow, object::CircleShape, renderStates::Ptr{Void})
    ccall((:sfRenderWindow_drawCircleShape, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), window.ptr, object.ptr, renderStates)
end

function draw(window::RenderWindow, object::RectangleShape, renderStates::Ptr{Void})
    ccall((:sfRenderWindow_drawRectangleShape, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), window.ptr, object.ptr, renderStates)
end

function draw(window::RenderWindow, object::Sprite, renderStates::Ptr{Void})
    ccall((:sfRenderWindow_drawSprite, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), window.ptr, object.ptr, renderStates)
end

function draw(window::RenderWindow, object::RenderText, renderStates::Ptr{Void})
    ccall((:sfRenderWindow_drawText, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), window.ptr, object.ptr, renderStates)
end

function draw(window::RenderWindow, object::ConvexShape, renderStates::Ptr{Void})
    ccall((:sfRenderWindow_drawConvexShape, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), window.ptr, object.ptr, renderStates)
end

function draw(window::RenderWindow, object::VertexArray, renderStates::Ptr{Void})
    ccall((:sfRenderWindow_drawVertexArray, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), window.ptr, object.ptr, renderStates)
end

function clear(window::RenderWindow, color::Color)
    ccall((:sfRenderWindow_clear, libcsfml_graphics), Void, (Ptr{Void}, Color,), window.ptr, color)
end

function display(window::RenderWindow)
    ccall((:sfRenderWindow_display, libcsfml_graphics), Void, (Ptr{Void},), window.ptr)
end

function close(window::RenderWindow)
    ccall((:sfRenderWindow_close, libcsfml_graphics), Void, (Ptr{Void},), window.ptr)
end

function destroy(window::RenderWindow)
    ccall((:sfRenderWindow_destroy, libcsfml_graphics), Void, (Ptr{Void},), window.ptr)
end

function capture(window::RenderWindow)
    return Image(ccall((:sfRenderWindow_capture, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), window.ptr))
end

function pixel2coords(window::RenderWindow, point::Vector2, targetview::View=get_view(window))
    point = to_vec2i(point)
    return ccall((:sfRenderWindow_mapPixelToCoords, libcsfml_graphics), Vector2f, (Ptr{Void}, Vector2i, Ptr{Void},), window.ptr, point, targetview.ptr)
end

function coords2pixel(window::RenderWindow, point::Vector2, targetview::View=get_view(window))
    point = to_vec2f(point)
    return ccall((:sfRenderWindow_mapCoordsToPixel, libcsfml_graphics), Vector2i, (Ptr{Void}, Vector2f, Ptr{Void},), window.ptr, point, targetview.ptr)
end


text:
using CLIMAParameters
using CLIMAParameters.Planet: planet_radius
using ..Diagnostics
using ..Mesh.Interpolation

"""
    DiagnosticsConfiguration

Container for all the `DiagnosticsGroup`s to be used for a simulation.
"""
mutable struct DiagnosticsConfiguration
    groups::Array{DiagnosticsGroup, 1}

    DiagnosticsConfiguration(
        groups::Array{DG, 1},
    ) where {DG <: DiagnosticsGroup} = new(groups)
end

"""
    InterpolationConfiguration(
        driver_config::DriverConfiguration,
        boundaries::Array,
        resolution::Tuple,
    )

Creates an `InterpolationTopology` (either an `InterpolationBrick` or an
`InterpolationCubedSphere`) to be used with a `DiagnosticsGroup`. The axes
are set up based on `boundaries` and `resolution`.
"""
function InterpolationConfiguration(
    driver_config::DriverConfiguration,
    boundaries::Array,
    resolution::Tuple,
)
    param_set = driver_config.bl.param_set
    grid = driver_config.grid
    if isa(grid.topology, StackedBrickTopology)

        axes = (
            collect(range(
                boundaries[1, 1],
                boundaries[2, 1],
                step = resolution[1],
            )),
            collect(range(
                boundaries[1, 2],
                boundaries[2, 2],
                step = resolution[2],
            )),
            collect(range(
                boundaries[1, 3],
                boundaries[2, 3],
                step = resolution[3],
            )),
        )
        return InterpolationBrick(grid, boundaries, axes[1], axes[2], axes[3])

    elseif isa(grid.topology, StackedCubedSphereTopology)

        FT = eltype(grid)
        _planet_radius::FT = planet_radius(param_set)
        info = driver_config.config_info
        vert_range = grid1d(
            _planet_radius,
            FT(_planet_radius + info.domain_height),
            nelem = info.nelem_vert,
        )
        axes = (
            collect(range(
                boundaries[1, 1],
                boundaries[2, 1],
                step = resolution[1],
            )),
            collect(range(
                boundaries[1, 2],
                boundaries[2, 2],
                step = resolution[2],
            )),
            collect(range(
                boundaries[1, 3],
                boundaries[2, 3],
                step = resolution[3],
            )),
        )
        return InterpolationCubedSphere(
            grid,
            collect(vert_range),
            info.nelem_horz,
            axes[1],
            axes[2],
            axes[3],
        )
    else
        @error "Cannot set up interpolation for this topology."
    end
end

"""
    InterpolationConfiguration(
        driver_config::DriverConfiguration,
        boundaries::Array,
        axes,
    )

Creates an `InterpolationTopology` (either an `InterpolationBrick` or an
`InterpolationCubedSphere`) to be used with a `DiagnosticsGroup`. The axes
are directly specified in lat/lon/lvl order.
"""
function InterpolationConfiguration(
    driver_config::DriverConfiguration,
    boundaries::Array,
    axes,
)
    param_set = driver_config.bl.param_set
    grid = driver_config.grid
    if isa(grid.topology, StackedBrickTopology)
        return InterpolationBrick(grid, boundaries, axes[1], axes[2], axes[3])

    elseif isa(grid.topology, StackedCubedSphereTopology)

        FT = eltype(grid)
        _planet_radius::FT = planet_radius(param_set)
        info = driver_config.config_info
        vert_range = grid1d(
            _planet_radius,
            FT(_planet_radius + info.domain_height),
            nelem = info.nelem_vert,
        )

        return InterpolationCubedSphere(
            grid,
            collect(vert_range),
            info.nelem_horz,
            axes[1],
            axes[2],
            axes[3],
        )
    else
        @error "Cannot set up interpolation for this topology."
    end
end

max_forks_repo_name: jm-c/CLIMA


text:
export replace_sjon
import Base: _replace, StringVector

# replace Char with function that compares with this Char
predicatepair(p::Pair{Char}) = equalto(p.first) => p.second

# replace collection of Char with function that checks occurrence in this collection
function predicatepair(p::Pair{<:Union{Tuple{Vararg{Char}},AbstractVector{Char},Set{Char}}})
    occursin(p.first) => p.second
end
predicatepair(p::Pair) = p

# find the first of a list of patterns.
# if several match at the same start position prefer the longer one
# if also the length is equal, prefer the first in list
# return index-range of the successful search and the pair
function findnextall!(pairs::NTuple{N,Pair}, s::AbstractString, st::Integer, rangecache) where N
    minstart=minend=endof(s)+1
    nmin=0
    for n = 1:N
        p = pairs[n]
        fr = first(rangecache[n])
        if fr ≥ st
            r = rangecache[n] # re-use previous search result
            last(r) == 0 && continue # not found
        else
            r = rangecache[n] = findnext(p.first, s, st)
            fr = first(r)
        end
        if fr > 0
            e = last(r)
            if fr < minstart || (fr == minstart && e > minend)
                minstart, minend, nmin = fr, e, n
            end
        else # not found
            rangecache[n] = endof(s)+1:0
        end
    end
    return minstart:minend, nmin
end

# cover the multiple pairs case
function replace_sjon(str::String, pat_repls_::Pair...; count::Integer=typemax(Int))
    if count == 0 || isempty(pat_repls_)
        return str
    end
    count < 0 && throw(DomainError(count, "`count` must be non-negative."))
    pat_repls = predicatepair.(pat_repls_)
    n = 1
    e = endof(str)
    i = a = start(str)
    rangecache = [0:0 for i=1:length(pat_repls)]
    r, npat = findnextall!(pat_repls, str, i, rangecache)
    j, k = first(r), last(r)
    out = IOBuffer(StringVector(floor(Int, 1.2sizeof(str))), true, true)
    out.size = 0
    out.ptr = 1
    while npat > 0 && n ≤ count
        if i == a || i <= k
            unsafe_write(out, pointer(str, i), UInt(j-i))
            pattern, repl = pat_repls[npat]
            _replace(out, repl, str, r, pattern)
        end
        if k < j
            i = j
            j > e && break
            k = nextind(str, j)
        else
            i = k = nextind(str, k)
        end
        r, npat = findnextall!(pat_repls, str, k, rangecache)
        j, k = first(r), last(r)
        n += 1
    end
    write(out, SubString(str, i))
    String(take!(out))
end


text:
@testset "Jacobi" begin

using SHIPs, Test

using SHIPs.JacobiPolys: Jacobi
using SHIPs: eval_basis, eval_basis_d

# copy-pasted from
# https://github.com/pjabardo/Jacobi.jl/blob/master/src/jac_poly.jl
# commit 9642d8f060203ddddeb17484a3aa8232022c7ba5

function jacobi(x, n, a, b)
   ox = one(x)
   zx = zero(x)
   if n==0
     return ox
   elseif n==1
     return ox/2 * (a - b + (a + b + 2)*x)
   end

   p0 = ox
   p1 = ox/2 * (a - b + (a + b + 2)*x)
   p2 = zx;

   for i = 1:(n-1)
      a1 = 2*(i+1)*(i+a+b+1)*(2*i+a+b);
      a2 = (2*i+a+b+1)*(a*a-b*b);
      a3 = (2*i+a+b)*(2*i+a+b+1)*(2*i+a+b+2);
      a4 = 2*(i+a)*(i+b)*(2*i+a+b+2);
      p2 = ox/a1*( (a2 + a3*x)*p1 - a4*p0);

      p0 = p1
      p1 = p2
   end

   return p2
end

djacobi(x, n, a, b) =  one(x)/2 * (a + b + n + 1) * jacobi(x, n-1, a+1, b+1)

@info("--------------  JACOBI TESTS ---------------")
@info("testing jacobi implementation against reference.")
for ntest = 1:30
   x = 2*rand() - 1
   α, β = rand(), rand()
   N = 30
   J = Jacobi(α, β, N, normalise=false)
   P = eval_basis(J, x)
   P1, dP = eval_basis_d(J, x)
   Ptest = [ jacobi(x, n, α, β) for n = 0:N ]
   dPtest = [ djacobi(x, n, α, β) for n = 0:N ]
   print_tf((@test P ≈ P1 ≈ Ptest))
   print_tf((@test dP ≈ dPtest))
end
println()

end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
export simple_psf
export jinc_psf

 # shift frequency to center if `shift=true`
 # and `norm` the sum to 1
function shift_and_norm(psf, shift, norm)
    if norm
        psf ./= sum(psf)
    end

    if shift
        psf = fftshift(psf)
    end
    return psf
end

"""
    simple_psf(psf_size, radius; shift=false)

Generation of an approximate 2D PSF.
`psf_size` is the output size of the PSF. The PSF will be centered
around the point [1, 1],
`radius` indicates the pupil diameter in pixel from which the PSF is generated.

# Examples
```julia-repl
julia> simple_psf([5, 5], 2)
5×5 Array{Float64,2}:
 0.36       0.104721    0.0152786    0.0152786    0.104721
 0.104721   0.0304627   0.00444444   0.00444444   0.0304627
 0.0152786  0.00444444  0.000648436  0.000648436  0.00444444
 0.0152786  0.00444444  0.000648436  0.000648436  0.00444444
 0.104721   0.0304627   0.00444444   0.00444444   0.0304627
```
"""
function simple_psf(psf_size, radius; shift=false)
    mask = rr(psf_size) .<= radius
    mask_ft = fft(mask)
    psf = abs2.(mask_ft)
    return shift_and_norm(psf, shift, true)
end

"""
    jinc_psf(psf_size, L, radius[, f]; λ=550e-9, shift=false)

Generate the normalized, incoherent 2D jinc PSF of a circular aperture.
`psf_size` is output array shape. `L` is the width of the array
expressed in meter. `radius` is the aperture radius in meter.
Keyword arguments `λ` and `f=100e-3` represent wavelength and focal
length of the lens respectively.

Reference: Mertz, J. (2019). Introduction to Optical Microscopy (2nd ed.).
"""
function jinc_psf(psf_size, L, radius, f=100e-3; λ=550e-9, shift=false)
    κ = calc_κ(λ)
    Δk⊥ = 2 * κ * radius / f
    # create real output psf

    psf = zeros(Float64, psf_size)
    # calculate each point
    for (j, x) in enumerate(fftpos(L, psf_size[2]))
        for (i, y) in enumerate(fftpos(L, psf_size[1]))
            r = sqrt(x^2 + y^2)
            psf[i, j] = jinc(π * Δk⊥ * r).^2
        end
    end
    # shift center to top left corner
    psf = ifftshift(psf)
    return shift_and_norm(psf, shift, true)
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
using Winds
using Test


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for graph conversions.

############################################################################################################################

@testset "Export" begin
   g = completegraph(10)

   setvprop!(g, :, 1:10, :p1)
   seteprop!(g, :, 1:90, :p1)

   @test export_adjacency(g) == indxs(g)
   @test export_vertex_property(g, :p1) == collect(1:10)

   sv = export_edge_property(g, :p1)
   @test nnz(sv) == 90
end

@testset "Graph Merging" begin
   g = completegraph(10)
   eit = edges(g)

   g1 = subgraph(g, eit[1:45])
   g2 = subgraph(g, eit[46:90])

   h = merge(g1, g2)
   @test h == g
end