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

using Base.Test
import Base.LinAlg: BlasReal, BlasFloat

n = 10 #Size of test matrix
srand(1)

@testset for relty in (Int, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
    if relty <: AbstractFloat
        dv = convert(Vector{elty}, randn(n))
        ev = convert(Vector{elty}, randn(n-1))
        if (elty <: Complex)
            dv += im*convert(Vector{elty}, randn(n))
            ev += im*convert(Vector{elty}, randn(n-1))
        end
    elseif relty <: Integer
        dv = convert(Vector{elty}, rand(1:10, n))
        ev = convert(Vector{elty}, rand(1:10, n-1))
        if (elty <: Complex)
            dv += im*convert(Vector{elty}, rand(1:10, n))
            ev += im*convert(Vector{elty}, rand(1:10, n-1))
        end
    end

    @testset "Constructors" begin
        @test Bidiagonal(dv,ev,'U') == Bidiagonal(dv,ev,true)
        @test_throws ArgumentError Bidiagonal(dv,ev,'R')
        @test_throws DimensionMismatch Bidiagonal(dv,ones(elty,n),true)
        @test_throws ArgumentError Bidiagonal(dv,ev)
    end

    @testset "getindex, setindex!, size, and similar" begin
        ubd = Bidiagonal(dv, ev, true)
        lbd = Bidiagonal(dv, ev, false)
        # bidiagonal getindex / upper & lower
        @test_throws BoundsError ubd[n + 1, 1]
        @test_throws BoundsError ubd[1, n + 1]
        @test ubd[2, 2] == dv[2]
        # bidiagonal getindex / upper
        @test ubd[2, 3] == ev[2]
        @test iszero(ubd[3, 2])
        # bidiagonal getindex / lower
        @test lbd[3, 2] == ev[2]
        @test iszero(lbd[2, 3])
        # bidiagonal setindex! / upper
        cubd = copy(ubd)
        @test_throws ArgumentError ubd[2, 1] = 1
        @test_throws ArgumentError ubd[3, 1] = 1
        @test (cubd[2, 1] = 0; cubd == ubd)
        @test ((cubd[1, 2] = 10) == 10; cubd[1, 2] == 10)
        # bidiagonal setindex! / lower
        clbd = copy(lbd)
        @test_throws ArgumentError lbd[1, 2] = 1
        @test_throws ArgumentError lbd[1, 3] = 1
        @test (clbd[1, 2] = 0; clbd == lbd)
        @test ((clbd[2, 1] = 10) == 10; clbd[2, 1] == 10)
        # bidiagonal setindex! / upper & lower
        @test_throws BoundsError ubd[n + 1, 1] = 1
        @test_throws BoundsError ubd[1, n + 1] = 1
        @test ((cubd[2, 2] = 10) == 10; cubd[2, 2] == 10)
        # bidiagonal size
        @test_throws ArgumentError size(ubd, 0)
        @test size(ubd, 1) == size(ubd, 2) == n
        @test size(ubd, 3) == 1
        # bidiagonal similar
        @test isa(similar(ubd), Bidiagonal{elty})
        @test isa(similar(ubd, Int), Bidiagonal{Int})
        @test isa(similar(ubd, Int, (3, 2)), Matrix{Int})
    end

    BD = Bidiagonal(dv, ev, true)
    @testset "show" begin
        dstring = sprint(Base.print_matrix,BD.dv')
        estring = sprint(Base.print_matrix,BD.ev')
        @test sprint(show,BD) == "$(summary(BD)):\n diag:$dstring\n super:$estring"
        BD = Bidiagonal(dv,ev,false)
        @test sprint(show,BD) == "$(summary(BD)):\n diag:$dstring\n sub:$estring"
    end

    @testset for isupper in (true, false)
        T = Bidiagonal(dv, ev, isupper)

        @testset "Constructor and basic properties" begin
            @test size(T, 1) == size(T, 2) == n
            @test size(T) == (n, n)
            @test Array(T) == diagm(dv) + diagm(ev, isupper?1:-1)
            @test Bidiagonal(Array(T), isupper) == T
            @test big.(T) == T
            @test Array(abs.(T)) == abs.(diagm(dv)) + abs.(diagm(ev, isupper?1:-1))
            @test Array(real(T)) == real(diagm(dv)) + real(diagm(ev, isupper?1:-1))
            @test Array(imag(T)) == imag(diagm(dv)) + imag(diagm(ev, isupper?1:-1))
        end
        z = zeros(elty, n)

        @testset for func in (conj, transpose, ctranspose)
            @test func(func(T)) == T
        end

        @testset "triu and tril" begin
            @test istril(Bidiagonal(dv,ev,'L'))
            @test !istril(Bidiagonal(dv,ev,'U'))
            @test tril!(Bidiagonal(dv,ev,'U'),-1) == Bidiagonal(zeros(dv),zeros(ev),'U')
            @test tril!(Bidiagonal(dv,ev,'L'),-1) == Bidiagonal(zeros(dv),ev,'L')
            @test tril!(Bidiagonal(dv,ev,'U'),-2) == Bidiagonal(zeros(dv),zeros(ev),'U')
            @test tril!(Bidiagonal(dv,ev,'L'),-2) == Bidiagonal(zeros(dv),zeros(ev),'L')
            @test tril!(Bidiagonal(dv,ev,'U'),1)  == Bidiagonal(dv,ev,'U')
            @test tril!(Bidiagonal(dv,ev,'L'),1)  == Bidiagonal(dv,ev,'L')
            @test tril!(Bidiagonal(dv,ev,'U'))    == Bidiagonal(dv,zeros(ev),'U')
            @test tril!(Bidiagonal(dv,ev,'L'))    == Bidiagonal(dv,ev,'L')
            @test_throws ArgumentError tril!(Bidiagonal(dv,ev,'U'),n+1)

            @test istriu(Bidiagonal(dv,ev,'U'))
            @test !istriu(Bidiagonal(dv,ev,'L'))
            @test triu!(Bidiagonal(dv,ev,'L'),1)  == Bidiagonal(zeros(dv),zeros(ev),'L')
            @test triu!(Bidiagonal(dv,ev,'U'),1)  == Bidiagonal(zeros(dv),ev,'U')
            @test triu!(Bidiagonal(dv,ev,'U'),2)  == Bidiagonal(zeros(dv),zeros(ev),'U')
            @test triu!(Bidiagonal(dv,ev,'L'),2)  == Bidiagonal(zeros(dv),zeros(ev),'L')
            @test triu!(Bidiagonal(dv,ev,'U'),-1) == Bidiagonal(dv,ev,'U')
            @test triu!(Bidiagonal(dv,ev,'L'),-1) == Bidiagonal(dv,ev,'L')
            @test triu!(Bidiagonal(dv,ev,'L'))    == Bidiagonal(dv,zeros(ev),'L')
            @test triu!(Bidiagonal(dv,ev,'U'))    == Bidiagonal(dv,ev,'U')
            @test_throws ArgumentError triu!(Bidiagonal(dv,ev,'U'),n+1)
        end

        Tfull = Array(T)
        @testset "Linear solves" begin
            if relty <: AbstractFloat
                c = convert(Matrix{elty}, randn(n,n))
                b = convert(Matrix{elty}, randn(n, 2))
                if (elty <: Complex)
                    b += im*convert(Matrix{elty}, randn(n, 2))
                end
            elseif relty <: Integer
                c = convert(Matrix{elty}, rand(1:10, n, n))
                b = convert(Matrix{elty}, rand(1:10, n, 2))
                if (elty <: Complex)
                    b += im*convert(Matrix{elty}, rand(1:10, n, 2))
                end
            end
            condT = cond(map(Complex128,Tfull))
            promty = typeof((zero(relty)*zero(relty) + zero(relty)*zero(relty))/one(relty))
            if relty != BigFloat
                x = T.'\c.'
                tx = Tfull.' \ c.'
                elty <: AbstractFloat && @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf))
                @test_throws DimensionMismatch T.'\b.'
                x = T'\c.'
                tx = Tfull' \ c.'
                @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf))
                @test_throws DimensionMismatch T'\b.'
                x = T\c.'
                tx = Tfull\c.'
                @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf))
                @test_throws DimensionMismatch T\b.'
            end
            @test_throws DimensionMismatch T \ ones(elty,n+1,2)
            @test_throws DimensionMismatch T.' \ ones(elty,n+1,2)
            @test_throws DimensionMismatch T' \ ones(elty,n+1,2)

            @test_throws DimensionMismatch T \ RowVector(ones(elty,n+1))
            @test_throws DimensionMismatch T.' \ RowVector(ones(elty,n+1))
            @test_throws DimensionMismatch T' \ RowVector(ones(elty,n+1))
            @test_throws DimensionMismatch Base.LinAlg.At_ldiv_B(T, RowVector(ones(elty,n+1)))
            @test_throws DimensionMismatch Base.LinAlg.Ac_ldiv_B(T, RowVector(ones(elty,n+1)))
            let bb = b, cc = c
                for atype in ("Array", "SubArray")
                    if atype == "Array"
                        b = bb
                        c = cc
                    else
                        b = view(bb, 1:n)
                        c = view(cc, 1:n, 1:2)
                    end
                end
                x = T \ b
                tx = Tfull \ b
                @test_throws DimensionMismatch Base.LinAlg.naivesub!(T,ones(elty,n+1))
                @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf))
                @testset "Generic Mat-vec ops" begin
                    @test T*b ≈ Tfull*b
                    @test T'*b ≈ Tfull'*b
                    if relty != BigFloat # not supported by pivoted QR
                        @test T/b' ≈ Tfull/b'
                    end
                end
            end
        end

        @testset "Round,float,trunc,ceil" begin
            if elty <: BlasReal
                @test floor.(Int, T) == Bidiagonal(floor.(Int, T.dv), floor.(Int, T.ev), T.isupper)
                @test isa(floor.(Int, T), Bidiagonal)
                @test trunc.(Int,T) == Bidiagonal(trunc.(Int, T.dv), trunc.(Int, T.ev), T.isupper)
                @test isa(trunc.(Int,T), Bidiagonal)
                @test round.(Int, T) == Bidiagonal(round.(Int, T.dv), round.(Int, T.ev), T.isupper)
                @test isa(round.(Int, T), Bidiagonal)
                @test ceil.(Int,T) == Bidiagonal(ceil.(Int,T.dv), ceil.(Int,T.ev), T.isupper)
                @test isa(ceil.(Int,T), Bidiagonal)
                @test floor.(T) == Bidiagonal(floor.(T.dv), floor.(T.ev), T.isupper)
                @test isa(floor.(T), Bidiagonal)
                @test trunc.(T) == Bidiagonal(trunc.(T.dv), trunc.(T.ev), T.isupper)
                @test isa(trunc.(T), Bidiagonal)
                @test round.(T) == Bidiagonal(round.(T.dv), round.(T.ev), T.isupper)
                @test isa(round.(T), Bidiagonal)
                @test ceil.(T) == Bidiagonal(ceil.(T.dv), ceil.(T.ev), T.isupper)
                @test isa(ceil.(T), Bidiagonal)
            end
        end

        @testset "Diagonals" begin
            @test diag(T,2) == zeros(elty, n-2)
            @test_throws ArgumentError diag(T,n+1)
        end

        @testset "Eigensystems" begin
            if relty <: AbstractFloat
                d1, v1 = eig(T)
                d2, v2 = eig(map(elty<:Complex ? Complex128 : Float64,Tfull))
                @test (isupper ? d1 : reverse(d1)) ≈ d2
                if elty <: Real
                    Test.test_approx_eq_modphase(v1, isupper ? v2 : v2[:,n:-1:1])
                end
            end
        end

        @testset "Singular systems" begin
            if (elty <: BlasReal)
                @test AbstractArray(svdfact(T)) ≈ AbstractArray(svdfact!(copy(Tfull)))
                @test svdvals(Tfull) ≈ svdvals(T)
                u1, d1, v1 = svd(Tfull)
                u2, d2, v2 = svd(T)
                @test d1 ≈ d2
                if elty <: Real
                    Test.test_approx_eq_modphase(u1, u2)
                    Test.test_approx_eq_modphase(v1, v2)
                end
                @test 0 ≈ vecnorm(u2*diagm(d2)*v2'-Tfull) atol=n*max(n^2*eps(relty),vecnorm(u1*diagm(d1)*v1'-Tfull))
                @inferred svdvals(T)
                @inferred svd(T)
            end
        end

        @testset "Binary operations" begin
            @test -T == Bidiagonal(-T.dv,-T.ev,T.isupper)
            @test convert(elty,-1.0) * T == Bidiagonal(-T.dv,-T.ev,T.isupper)
            @test T * convert(elty,-1.0) == Bidiagonal(-T.dv,-T.ev,T.isupper)
            @testset for isupper2 in (true, false)
                dv = convert(Vector{elty}, relty <: AbstractFloat ? randn(n) : rand(1:10, n))
                ev = convert(Vector{elty}, relty <: AbstractFloat ? randn(n-1) : rand(1:10, n-1))
                T2 = Bidiagonal(dv, ev, isupper2)
                Tfull2 = Array(T2)
                for op in (+, -, *)
                    @test Array(op(T, T2)) ≈ op(Tfull, Tfull2)
                end
            end
        end

        @test inv(T)*Tfull ≈ eye(elty,n)
    end

    @test Matrix{Complex{Float64}}(BD) == BD

end

# Issue 10742 and similar
let A = Bidiagonal([1,2,3], [0,0], true)
    @test istril(A)
    @test isdiag(A)
end

# test construct from range
@test Bidiagonal(1:3, 1:2, true) == [1 1 0; 0 2 2; 0 0 3]

@testset "promote_rule" begin
    A = Bidiagonal(ones(Float32,10),ones(Float32,9),true)
    B = rand(Float64,10,10)
    C = Tridiagonal(rand(Float64,9),rand(Float64,10),rand(Float64,9))
    @test promote_rule(Matrix{Float64}, Bidiagonal{Float64}) == Matrix{Float64}
    @test promote(B,A) == (B, convert(Matrix{Float64}, A))
    @test promote(C,A) == (C,Tridiagonal(zeros(Float64,9),convert(Vector{Float64},A.dv),convert(Vector{Float64},A.ev)))
end

import Base.LinAlg: fillslots!, UnitLowerTriangular
@testset "fill! and fillslots!" begin
    let #fill!
        let # fillslots!
            A = Tridiagonal(randn(2), randn(3), randn(2))
            @test fillslots!(A, 3) == Tridiagonal([3, 3.], [3, 3, 3.], [3, 3.])
            B = Bidiagonal(randn(3), randn(2), true)
            @test fillslots!(B, 2) == Bidiagonal([2.,2,2], [2,2.], true)
            S = SymTridiagonal(randn(3), randn(2))
            @test fillslots!(S, 1) == SymTridiagonal([1,1,1.], [1,1.])
            Ult = UnitLowerTriangular(randn(3,3))
            @test fillslots!(Ult, 3) == UnitLowerTriangular([1 0 0; 3 1 0; 3 3 1])
        end
        let # fill!(exotic, 0)
            exotic_arrays = Any[Tridiagonal(randn(3), randn(4), randn(3)),
            Bidiagonal(randn(3), randn(2), rand(Bool)),
            SymTridiagonal(randn(3), randn(2)),
            sparse(randn(3,4)),
            Diagonal(randn(5)),
            sparse(rand(3)),
            LowerTriangular(randn(3,3)),
            UpperTriangular(randn(3,3))
            ]
            for A in exotic_arrays
                fill!(A, 0)
                for a in A
                    @test a == 0
                end
            end
        end
        let # fill!(small, x)
            val = randn()
            b = Bidiagonal(randn(1,1), true)
            st = SymTridiagonal(randn(1,1))
            for x in (b, st)
                @test Array(fill!(x, val)) == fill!(Array(x), val)
            end
            b = Bidiagonal(randn(2,2), true)
            st = SymTridiagonal(randn(3), randn(2))
            t = Tridiagonal(randn(3,3))
            for x in (b, t, st)
                @test_throws ArgumentError fill!(x, val)
                @test Array(fill!(x, 0)) == fill!(Array(x), 0)
            end
        end
    end
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
"""
    bernoulli(n)

 Calculates the first 35 Bernoulli numbers ``B_n``  (of the first-kind or NIST type)
 e.g., see

 + [http://mathworld.wolfram.com/BernoulliNumber.html](http://mathworld.wolfram.com/BernoulliNumber.html)
 + [https://en.wikipedia.org/wiki/Bernoulli_number](https://en.wikipedia.org/wiki/Bernoulli_number)
 + [http://dlmf.nist.gov/24](http://dlmf.nist.gov/24)

 N.B. Bernoulli numbers of second kind only seem to differ in that ``B_1 = + 1/2`` (instead of -1/2)

## Arguments
* ``n`` `::Integer`: the index into the series, ``n=0,1,2,3,...,35`` (for larger ``n`` use `bernoulli(n,0.0)` )

 We only provide the 1st 36 values as beyond this, we can't
 return Int64 rationals, so best to compute the real approximation
 using `bernoulli(n,0.0)`.

 Odd values for ``n>1`` are all zero.

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> bernoulli(6)
1//42
```
"""
function bernoulli(n::Integer)
    # this just does a lookup -- seemed like it would be easier to code and faster
    # for the size of numbers I am working with
    if n<0
        throw(DomainError(n))
    elseif n > 1 && isodd(n)
        return 0 // 1
    elseif n > 35
        throw(DomainError(n, "If n > 35, then the numerator needs Int128 at least, and worse, so this code is not the code you want. Try using bernoulli(n, 0.0) to get a floating point approximation to the result."))
    end

    # Denominator of Bernoulli number B_n
    #   http://oeis.org/A027642
    D = [2, 6, 1, 30, 1, 42, 1, 30, 1, 66, 1, 2730, 1, 6, 1, 510, 1, 798, 1, 330, 1, 138, 1, 2730, 1, 6, 1, 870, 1, 14322, 1, 510, 1, 6, 1, 1919190, 1, 6, 1, 13530, 1, 1806, 1, 690, 1, 282, 1, 46410, 1, 66, 1, 1590, 1, 798, 1, 870, 1, 354, 1, 56786730]

    # Numerator of Bernoulli number B_n (storing 62 of these because they are easy)
    #   http://oeis.org/A027641
    N = [-1, 1, 0, -1, 0, 1, 0, -1, 0, 5, 0, -691, 0, 7, 0, -3617, 0, 43867, 0, -174611, 0, 854513, 0, -236364091, 0, 8553103, 0, -23749461029, 0, 8615841276005, 0, -7709321041217, 0, 2577687858367, 1]

    if n==0
        return 1
    else
        return N[n] // D[n]
    end
end

ext: jl


text:
using Test
using FdeSolver

@testset "FdeSolver.jl" begin
    b = FdeSolver.greet()
    @test @isdefined(b)
end


text:
export compute_error

function compute_error(u, v, data::DataSet)

    str3 = "donnee_"
    str5 = ".txt"

    epsilon = data.epsilon

    if (epsilon == 10       )  str4 = "10"        end
    if (epsilon == 5        )  str4 = "5"         end
    if (epsilon == 2.5      )  str4 = "2_5"       end
    if (epsilon == 1        )  str4 = "1"         end
    if (epsilon == 0.5      )  str4 = "0_5"       end
    if (epsilon == 0.2      )  str4 = "0_2"       end
    if (epsilon == 0.25     )  str4 = "0_25"      end
    if (epsilon == 0.1      )  str4 = "0_1"       end
    if (epsilon == 0.05     )  str4 = "0_05"      end
    if (epsilon == 0.025    )  str4 = "0_025"     end
    if (epsilon == 0.01     )  str4 = "0_01"      end
    if (epsilon == 0.005    )  str4 = "0_005"     end
    if (epsilon == 0.0025   )  str4 = "0_0025"    end
    if (epsilon == 0.001    )  str4 = "0_001"     end
    if (epsilon == 0.0005   )  str4 = "0_0005"    end
    if (epsilon == 0.00025  )  str4 = "0_00025"   end
    if (epsilon == 0.0001   )  str4 = "0_0001"    end
    if (epsilon == 0.00005  )  str4 = "0_00005"   end
    if (epsilon == 0.000025 )  str4 = "0_000025"  end
    if (epsilon == 0.00001  )  str4 = "0_00001"   end
    if (epsilon == 0.000005 )  str4 = "0_000005"  end
    if (epsilon == 0.0000025)  str4 = "0_0000025" end
    if (epsilon == 0.000001 )  str4 = "0_000001"  end

    ref_file = joinpath("test", "donnees_data3_128_micmac/", str3 * str4 * str5)

    ndata = 128
    uv    = zeros(Float64, (ndata, 4))

    open(ref_file) do f

        for (j,line) in enumerate(eachline(f))
            for (i, val) in enumerate( [ parse(Float64, val) for val in split(line)])
                uv[j, i] = val
            end
        end

    end

    nx   = data.nx
    xmin = data.xmin
    xmax = data.xmax
    T    = data.T
    x    = data.x
    dx   = (xmax - xmin) / nx
    L    = xmax - xmin
    kx   = zeros(Float64, nx)
    kx   = 2π / (xmax - xmin) * vcat(0:nx÷2-1,-nx÷2:-1)

    ua = zeros(ComplexF64, (nx, 4))
    va = fft(uv, 1) / ndata
    k  = zeros(Float64, ndata)
    k .= 2π / L * vcat(0:ndata÷2-1,-ndata÷2:-1)

    for j in 1:ndata
        vv  = transpose(va[j, :])
	    ua .= ua .+ vv .* exp.(1im * k[j] * (x .- xmin))
    end

    uref = ua[:, 1] .+ 1im * ua[:, 2]
    vref = ua[:, 3] .+ 1im * ua[:, 4]

    refH1 = sqrt(dx * norm(ifft(1im * sqrt.(1 .+ kx.^2) .* fft(uref,1),1))^2
               + dx * norm(ifft(1im * sqrt.(1 .+ kx.^2) .* fft(vref,1),1))^2)

    err  = (sqrt(dx * norm(ifft(1im * sqrt.(1 .+ kx.^2) .* fft(u .- uref,1),1))^2
               + dx * norm(ifft(1im * sqrt.(1 .+ kx.^2) .* fft(v .- vref,1),1))^2)) / refH1

    err

end


text:
"""
    DataDistribution(X; kwargs...)

A probability distribution representing the frequency of codewords in a given
sample. The sample must be a `Matrix{Bool}` or a `BitMatrix` with columns
representing codewords. The number of bits = `size(X,1)` must be less than 64.

First sets `X = X[I,:]` then proceeds

"""
mutable struct DataDistribution <: AbstractBinaryVectorDistribution
    P::SparseVector{Float64}
    N::Int
    X::BitMatrix
    # I::Vector{Int} # moving this to metadata
    metadata::Dict{Any,Any}
    cache::Dict{Any,Any}

    function DataDistribution(X::AbstractMatrix{Bool}; kwargs...)
        (N,m) = size(X)
        if N >= 64
            error("Attempting to compute data distribution of data with >= 64 bit words")
        end
        P = sparsevec(_compute_Xtilde(X)[:], ones(m), 2^N)
        new(P/m, N, BitMatrix(X), Dict(kwargs), Dict(kwargs))
    end
end
_compute_Xtilde(X::AbstractMatrix{Bool}) = 1 .+ sum(X .* [2^i for i = 0:(size(X,1)-1)], 1)
_compute_Xtilde(X::BitMatrix) = [1 + Int(X[:,k].chunks[1]) for k = 1:size(X,2)]

n_bits(DD::DataDistribution) = DD.N

==(P1::DataDistribution, P2::DataDistribution) = P1.P == P2.P

# pdf(Pr::DataDistribution, x::Vector{Bool}) = length(x) == n_bits(Pr) ? Pr.P[1 + dot([2^i for i = 0:(Pr.N - 1)], x)] : error("DataDistribution pdf: out of domain error")
# pdf(Pr::DataDistribution, x::BitVector) = length(x) == n_bits(Pr) ? Pr.P[1 + Int(x.chunks[1])] : error("DataDistribution pdf: out of domain error")
pdf(Pr::DataDistribution, x) = length(x) == n_bits(Pr) ? Pr.P[1 + _binary_to_int(x)] : error("DataDistribution pdf: out of domain error")

get_pdf(DD::DataDistribution) = full(DD.P)

function show(io::IO, P::DataDistribution)
    println(io, "Data Distribution")
    println(io, "N_neurons: $(n_bits(P))")
    # println(io, "Indices:   $(P.I)")
    println(io, "N_samples: $(size(P.X,2))")
    show_metadata(io, P)
end

expectation_matrix(DD::DataDistribution) = DD.X * DD.X' / size(DD.X,2)

entropy(DD::DataDistribution) = -sum_kbn([p * log(p) for p in nonzeros(DD.P)])
entropy2(DD::DataDistribution) = -sum_kbn([p * log2(p) for p in nonzeros(DD.P)])

# random(DD::DataDistribution, n_samples::Int=1) = _random_exact(DD, n_samples)

max_stars_repo_name: sekunder/BinaryVectorProbability

max_issues_repo_name: sekunder/BinaryVectorProbability

max_forks_repo_name: sekunder/BinaryVectorProbability


text:
#=
Sections of the code below were proposed by Simon Danisch:
    https://gist.github.com/SimonDanisch/ae046f3a0c78b26242e78fa9b139aa11#file-benchmark-jl
=#

using BenchmarkTools

"""
    regular_time_step(u::Matrix)

Take a time step in the Jacobi numerical approximation u using un-optimized loops.
"""
function regular_time_step(u::Matrix)
    n, m = size(u)
    error = 0.0
    for i in 2:n-1
        for j in 2:m-1
            temp = u[i, j]
            u[i, j] = ((u[i-1, j] + u[i+1, j] +
                       u[i, j-1] + u[i, j+1])*4.0 +
                       u[i-1, j-1] + u[i-1, j+1] +
                       u[i+1, j-1] + u[i+1, j+1])/20.0

            difference = u[i, j] - temp
            error = error + difference*difference
        end
    end
    return sqrt(error)
end

"""
    optimized_time_step{T}(u::Matrix{T})

Take a time step in the Jacobi numerical approximation u using optimized loops.
"""
function optimized_time_step{T}(u::Matrix{T})
    n, m = size(u)
    error = T(0.0) # not an optimization, but makes it work for all element types
    @inbounds for i = 2:n-1
        for j = 2:m-1
            temp = u[j, i];
            u[j, i] = ( T(4) *(u[j-1, i] + u[j+1, i] +
                        u[j, i-1] + u[j, i+1]) +
                        u[j-1, i-1] + u[j+1, i+1] +
                        u[j+1, i-1] + u[j-1, i+1]) / T(20)
            difference = u[j, i] - temp
            error += difference*difference;
        end
    end
    return sqrt(error)
end

"""
    optimized_time_step{T}(u::Matrix{T})

Take a time step in the Jacobi numerical approximation u using optimized loops and SIMD.
"""
function optimized_time_step_simd{T}(u::Matrix{T})
    n, m = size(u)
    error = T(0.0) # not an optimization, but makes it work for all element types
    @inbounds for i = 2:n-1
        @simd for j = 2:m-1
            temp = u[j, i];
            u[j, i] = (T(4)*(u[j-1, i] + u[j+1, i] +
                        u[j, i-1] + u[j, i+1]) +
                        u[j-1, i-1] + u[j+1, i+1] +
                        u[j+1, i-1] + u[j-1, i+1])/T(20)
            difference = u[j, i] - temp
            error += difference*difference;
        end
    end
    return sqrt(error)
end

"""
    optimized_vectorized_time_step(u::Matrix)

Take a time step in the Jacobi numerical approximation u using optimized vectorization.
"""
function optimized_vectorized_time_step(u::Matrix)
    (n, m) = size(u)
    u_old = copy(u)
    #@views u[2:n-1, 2:m-1] = ((u[1:n-2, 2:m-1] + u[3:n, 2:m-1] +
    u[2:n-1, 2:m-1] = ((u[1:n-2, 2:m-1] + u[3:n, 2:m-1] +
                        u[2:n-1, 1:m-2] + u[2:n-1, 3:m])*4.0 +
                        u[1:n-2, 1:m-2] + u[1:n-2, 3:m] +
                        u[3:n, 1:m-2] + u[3:n, 3:m])/20.0
    v = vec(u - u_old)
    return sqrt(dot(v, v))
end

# C version
const clib_name5 = tempname()
source = """
#include <math.h>
double c_time_step(int N, double *u[])
{
    double temp, difference;
    int i, j;
    double error;
    error = 0.0;
    for (i = 1; i < N-1; i++) {
        for (j = 1; j < N-1; j++) {
            temp = u[j][i];
            u[j][i] = ( 4.0 *(u[j-1][i] + u[j+1][i] +
                        u[j][i-1] + u[j][i+1]) +
                        u[j-1][i-1] + u[j+1][i+1] +
                        u[j+1][i-1] + u[j-1][i+1])/20.;
            difference = u[j][i] - temp;
            error += difference*difference;
        }
    }
    return sqrt(error);
}
"""
run(`gcc --version`)
# gcc (Ubuntu 5.4.1-2ubuntu1~16.04) 5.4.1 20160904
open(`gcc -fPIC -Ofast -O3 -xc -shared -o $(clib_name5).so -`, "w") do pipe
    print(pipe, source)
end

"""
    c_time_step(X::Vector{Vector{Float64}})

Take a time step in the Jacobi numerical approximation u for the C version.
"""
function c_time_step(X::Vector{Vector{Float64}})
    ccall(("c_time_step", clib_name5), Float64, (Cint, Ptr{Ptr{Float64}}), length(X), X)
end

"""
    solver(u, time_step)

Find the numerical approximation u by iterating the time steps.
"""
function solver(u, time_step)
    iteration = 0
    error = 2.0
    while iteration < 100000 && error > 1e-6
        error = time_step(u)
        iteration += 1
    end
    u, error, iteration
end

# Get the number of points from the command line
n, = size(ARGS)
if (n < 1)
    println("Usage: jacobiDriver4-benchmark.jl N")
    println("       --->  Please specify the number of points.")
    exit()
end

num_points = parse(Int, ARGS[1])

println("-------------------------------")
println(@sprintf "Jacobi Iterations: %d" num_points)
println("-------------------------------")

"""
    initialize_array(numPoints = num_points)

Create and populate a square array with side length num_points.
"""
#function initialize_array(numPoints = 100)
function initialize_array(numPoints = num_points)
    pi_c = 4.0 * atan(1.0)
    v = zeros(Float64, numPoints, numPoints)
    x = [(i-1)*pi_c/(numPoints-1) for i in 1:numPoints]
    v[1, :] = [sin(a) for a in x]
    v[numPoints, :] = [sin(a)*exp(-pi_c) for a in x]
    v
end

initialize_c() = (v = initialize_array(); v_2 = [v[i, :] for i=1:size(v, 1)])

functions = (regular_time_step, optimized_time_step, optimized_vectorized_time_step, optimized_time_step_simd) #
benchmarks = map(functions) do time_step
    v = initialize_array()
    solver(v, time_step) # compile
    v = initialize_array()
    # should be enough iterations already to be okay with just a simple elapsed
    result = @elapsed begin
        u, error, iteration = solver(v, time_step)
    end
    Symbol(time_step) => (result, iteration)
end

result = @elapsed begin
    v = initialize_c()
    u, error, iteration = solver(v, c_time_step)
end
benchmarks = (benchmarks..., :c_time_step => (result, iteration))
for (name, (t, sweeps)) in benchmarks
    println(name, ":")
    println("          Number of sweeps: ", sweeps)
    println("              Elapsed Time: $(t) Seconds")
end


text:
module AssociatedLegendreFunction

import ..OffsetVec, ..OffsetMat

### The code in this module is based on a paper by Xing et al.:
### https://doi.org/10.1007/s00190-019-01331-0.  All references to equation numbers are for that
### paper.

# These functions implement Eqs. (13), absorbing a factor of 1/2 into b, and converting to the
# relevant float type where needed.
@inline a(n, ::Type{T}) where {T<:Real} = √((2n+1)/T(2n-1))
@inline b(n, ::Type{T}) where {T<:Real} = √((n-1)*(2n+1)/T(2n*(2n-1)))
@inline c(n, m, ::Type{T}) where {T<:Real} = √((n+m)*(n-m)*(2n+1) / T(2n-1)) / n
@inline d(n, m, ::Type{T}) where {T<:Real} = √((n-m)*(n-m-1)*(2n+1) / T(2n-1)) / (2n)
@inline e(n, m, ::Type{T}) where {T<:Real} = √((n+m)*(n+m-1)*(m!=1 ? 1 : 2)*(2n+1) / T(2n-1)) / (2n)

# These are more efficient to compute than the above, because they involve no division and fewer T-multiplications
@inline b̂(n, ::Type{T}) where {T<:Real} = √T(2n*(n-1))  # b = b̂ * a / 2n
@inline ĉ(n, m, ::Type{T}) where {T<:Real} = √T(4*(n+m)*(n-m))  # c = ĉ * a / 2n
@inline d̂(n, m, ::Type{T}) where {T<:Real} = √T((n-m)*(n-m-1))  # d = d̂ * a / 2n
@inline ê(n, m, ::Type{T}) where {T<:Real} = √T((n+m)*(n+m-1)*(m!=1 ? 1 : 2))  # e = ê * a / 2n

struct ALFRecursionCoefficients{T<:Real}
    nmax::Int
    a::Vector{T}
    b::Vector{T}
    cde::Matrix{T}
end
function ALFRecursionCoefficients(nmax, ::Type{T}=Float64) where {T<:Real}
    avalues = Vector{T}(undef, nmax+1)
    bvalues = Vector{T}(undef, nmax+1)
    cdevalues = Matrix{T}(undef, 3, (nmax*(nmax+3))÷2+1)
    i = 0
    @inbounds for n in 1:(nmax+1)
        avalues[n] = a(n, T)
        bvalues[n] = b(n, T)
        for m in 1:n
            i+=1
            cdevalues[1, i] = c(n, m, T)
            cdevalues[2, i] = d(n, m, T)
            cdevalues[3, i] = e(n, m, T)
        end
    end
    ALFRecursionCoefficients{T}(nmax, avalues, bvalues, cdevalues)
end
Base.show(io::IO, r::ALFRecursionCoefficients) = print(io, "$(typeof(r))($(r.nmax))")

function ALFrecurse!(P̄ₙ::AbstractVector{T}, P̄ₙ₋₁::AbstractVector{T}, n::Int, t::T, u::T, aₙ::T, bₙ::T, cdeₙ::AbstractMatrix{T}) where {T<:Real}
    @fastmath @inbounds begin
        let m=0
            # Eq. (12a); note that we absorbed a factor of 1/2 into the definition of b
            P̄ₙ[m] = t * aₙ * P̄ₙ₋₁[m] - u * bₙ * P̄ₙ₋₁[m+1]
        end
        for m in 1:n-2
            # Eq. (12b)
            cₙₘ = cdeₙ[1, m]
            dₙₘ = cdeₙ[2, m]
            eₙₘ = cdeₙ[3, m]
            P̄ₙ[m] = t * cₙₘ * P̄ₙ₋₁[m] - u * (dₙₘ * P̄ₙ₋₁[m+1] - eₙₘ * P̄ₙ₋₁[m-1])
        end
        let m=n-1
            # Eq. (12b)
            cₙₘ = cdeₙ[1, m]
            eₙₘ = cdeₙ[3, m]
            P̄ₙ[m] = t * cₙₘ * P̄ₙ₋₁[m] + u * eₙₘ * P̄ₙ₋₁[m-1]
        end
        let m=n
            # Eq. (12b)
            eₙₘ = cdeₙ[3, m]
            P̄ₙ[m] = u * eₙₘ * P̄ₙ₋₁[m-1]
        end
    end
end

function ALFrecurse!(P̄ₙ::AbstractVector{T}, P̄ₙ₋₁::AbstractVector{T}, n::Int, t::T, u::T) where {T<:Real}
    @fastmath @inbounds begin
        aₙ = a(n, T)
        aₙ_over_2n = aₙ / T(2n)
        bₙ = b̂(n, T) * aₙ_over_2n
        let m=0
            # Eq. (12a); note that we absorbed a factor of 1/2 into the definition of b
            P̄ₙ[m] = t * aₙ * P̄ₙ₋₁[m] - u * bₙ * P̄ₙ₋₁[m+1]
        end
        for m in 1:n-2
            # Eq. (12b)
            ĉₙₘ = ĉ(n, m, T)
            d̂ₙₘ = d̂(n, m, T)
            êₙₘ = ê(n, m, T)
            P̄ₙ[m] = (t * ĉₙₘ * P̄ₙ₋₁[m] - u * (d̂ₙₘ * P̄ₙ₋₁[m+1] - êₙₘ * P̄ₙ₋₁[m-1])) * aₙ_over_2n
        end
        let m=n-1
            # Eq. (12b)
            ĉₙₘ = ĉ(n, m, T)
            êₙₘ = ê(n, m, T)
            P̄ₙ[m] = (t * ĉₙₘ * P̄ₙ₋₁[m] + u * êₙₘ * P̄ₙ₋₁[m-1]) * aₙ_over_2n
        end
        let m=n
            # Eq. (12b)
            êₙₘ = ê(n, m, T)
            P̄ₙ[m] = u * êₙₘ * P̄ₙ₋₁[m-1] * aₙ_over_2n
        end
    end
end

"""
    ALFcompute(expiβ, nmax)
    ALFcompute!(P̄, expiβ, nmax)
    ALFcompute(expiβ, nmax, recursion_coefficients)
    ALFcompute!(P̄, expiβ, nmax, recursion_coefficients)

Compute the "fully normalized" Associated Legendre Functions up to some maximum `n` value `nmax`.

These functions can take a vector P̄, to store the data, stored in order of increasing `m` most
rapidly varying and then increasing `n`.  If not supplied, P̄ will be constructed for you and
returned.

The optional `recursion_coefficients` argument must be an `ALFRecursionCoefficients`, which stores
various constant coefficients used in the recursion.  This object requires more than 3x the memory
and more than 20x the time to compute a single P̄ vector without this argument, but passing it will
typically speed up the calculation of each P̄ by a factor of 8x or so.  Thus, if you expect to
compute P̄ more than a few times, it will take less time to pre-compute those factors, and pass them
to this function.

Note that the base real type will be inferred from the (complex) type of `expiβ`.  If present, the
base types of `P̄` and `recursion_coefficients` must agree.

"""
function ALFcompute!(P̄::Vector{T}, expiβ::Complex{T}, nmax::Int, recursion_coefficients::ALFRecursionCoefficients{T}) where {T<:Real}
    min_length = (nmax * (nmax + 3)) ÷ 2 + 1
    if length(P̄) < min_length
        throw(ArgumentError("Length of P̄ ($(length(P̄))) must be at least $min_length for nmax=$nmax"))
    end

    P̄ₙ₋₁ = OffsetVec(P̄, -1)  # -(n*(n-1))÷2-1
    P̄ₙ = OffsetVec(P̄, -1)  # -(n*(n+1))÷2-1
    cdeₙ = OffsetMat(recursion_coefficients.cde, 0, 0)  # -(n*(n-1))÷2

    @inbounds begin
        t = expiβ.re  # = cosβ
        u = expiβ.im  # = sinβ

        # Initialize with Eq. (14)
        if nmax ≥ 0
            P̄ₙ[0] = one(T)  # n=0
        end
        if nmax ≥ 1
            P̄ₙ.offset -= 1
            sqrt3 = √T(3)
            P̄ₙ[0] = sqrt3 * t
            P̄ₙ[1] = sqrt3 * u
        end
        for n in 2:min(nmax, recursion_coefficients.nmax)
            aₙ = recursion_coefficients.a[n]
            bₙ = recursion_coefficients.b[n]
            cdeₙ.offset2 -= (n-1)
            P̄ₙ₋₁.offset = P̄ₙ.offset
            P̄ₙ.offset -= n
            ALFrecurse!(P̄ₙ, P̄ₙ₋₁, n, t, u, aₙ, bₙ, cdeₙ)
        end
        if nmax > recursion_coefficients.nmax
            for n in max(2,recursion_coefficients.nmax+1):nmax
                P̄ₙ₋₁.offset = P̄ₙ.offset
                P̄ₙ.offset -= n
                ALFrecurse!(P̄ₙ, P̄ₙ₋₁, n, t, u)
            end
        end
    end
end
function ALFcompute!(P̄::Vector{T}, expiβ::Complex{T}, nmax::Int) where {T<:Real}
    recursion_coefficients = ALFRecursionCoefficients(0, T)
    ALFcompute!(P̄, expiβ, nmax, recursion_coefficients)
end

function ALFcompute(expiβ::Complex{T}, nmax::Int, recursion_coefficients::ALFRecursionCoefficients{T}) where {T<:Real}
    min_length = (nmax * (nmax + 3)) ÷ 2 + 1
    P̄ = Vector{T}(undef, min_length)
    ALFcompute!(P̄, expiβ, nmax, recursion_coefficients)
    P̄
end
function ALFcompute(expiβ::Complex{T}, nmax::Int) where {T<:Real}
    recursion_coefficients = ALFRecursionCoefficients(0, T)
    ALFcompute(expiβ, nmax, recursion_coefficients)
end

end  # module ALF


text:
using DrWatson
@quickactivate :ConditionalNFs4Imaging

using InvertibleNetworks
using ProgressMeter
using PyPlot
using Seaborn
using LinearAlgebra
using Random
using Statistics
using Flux
include(scriptsdir("create_OOD_shots_and_RTM.jl"))

# Random seed
Random.seed!(19)

args = read_config("amortized_MAP.json")
args = parse_input_args(args)

max_epoch = args["max_epoch"]
lr = args["lr"]
idx = args["idx"]
sigma = args["sigma"]
nsrc = args["nsrc"]
nrec = args["nrec"]
pretrained_model_config = args["pretrained_model_config"]
sim_name = args["sim_name"]
resume = args["resume"]

# Loading the existing weights, if any.
loaded_keys =
    resume_from_checkpoint(args, ["J", "Z0", "src_enc", "opt", "epoch", "fval", "mval"])
J = loaded_keys["J"]
Z0 = loaded_keys["Z0"]
src_enc = loaded_keys["src_enc"]
opt = loaded_keys["opt"]
init_epoch = loaded_keys["epoch"]
fval = loaded_keys["fval"]
mval = loaded_keys["mval"]

# Create shot data and RTM image given a index for out-of-distribution seismic image.
X_OOD, lin_data, Y_OOD = create_OOD_shots_and_RTM(args)

# Create simultaneous source operators and return the source encoding.
if J == nothing
    J, _, src_enc =
        create_operator(; nrec = nrec, nsrc = nsrc, sigma = sigma, sim_src = true)
end

# Mix the data according to the source encoding to create simultaneous data.
mix_data!(lin_data, src_enc)

# Loading the pretrained model.
pretrained_model_args = read_config(pretrained_model_config)
loaded_keys = load_experiment(
    pretrained_model_args,
    ["Params", "AN_params_x", "AN_params_y", "n_hidden", "depth"],
)
Params = loaded_keys["Params"]
AN_params_x = loaded_keys["AN_params_x"]
AN_params_y = loaded_keys["AN_params_y"]
n_hidden = loaded_keys["n_hidden"]
depth = loaded_keys["depth"]

AN_x = ActNorm(1)
AN_y = ActNorm(1)
put_params!(AN_x, convert(Array{Any,1}, AN_params_x))
put_params!(AN_y, convert(Array{Any,1}, AN_params_y))

X_OOD = AN_x.forward(X_OOD)
Y_OOD = AN_y.forward(Y_OOD)
X_OOD = wavelet_squeeze(X_OOD)
Y_OOD = wavelet_squeeze(Y_OOD)

nx, ny, n_in = size(X_OOD)[1:3]

# Create network
CH = NetworkConditionalHINT(n_in, n_hidden, depth, logdet = false)
put_params!(CH, convert(Array{Any,1}, Params))
CH = CH |> gpu

# Corresponding data latent variable
Y_OOD = Y_OOD |> gpu
X_OOD = X_OOD |> gpu
Zy = CH.forward_Y(Y_OOD)

CHrev = reverse(CH)
AN_x_rev = reverse(deepcopy(AN_x))

# Draw new Zx, while keeping Zy fixed
Z0 == nothing && (Z0 = zeros(Float32, nx, ny, n_in, 1))
Z0 = Z0 |> gpu

X0 = CHrev.forward(Z0, Zy)[1]

data_loader = Flux.DataLoader(range(1, nsrc, step = 1), batchsize = 1, shuffle = true)
num_batches = length(data_loader)

# Optimizer
opt == nothing && (opt = Flux.ADAM(lr))

p = Progress(num_batches * (max_epoch - init_epoch + 1))
fval == nothing && (fval = zeros(Float32, num_batches * max_epoch))
mval == nothing && (mval = zeros(Float32, num_batches * max_epoch))
ΔZx = zeros(Float32, nx, ny, n_in, 1) |> gpu

for epoch = init_epoch:max_epoch
    for (itr, idx) in enumerate(data_loader)
        Base.flush(Base.stdout)

        fval[(epoch-1)*num_batches+itr], X_hat = amortized_MAP_loss!(
            ΔZx,
            CHrev,
            Z0,
            Zy,
            lin_data[idx..., :],
            J[idx...],
            sigma,
            AN_x_rev,
        )

        Flux.update!(opt, vec(Z0), vec(ΔZx))

        mval[(epoch-1)*num_batches+itr] = norm(X_hat - X_OOD)^2 / length(X_hat)

        ProgressMeter.next!(
            p;
            showvalues = [
                (:Epoch, epoch),
                (:Itreration, itr),
                (:Objective, fval[(epoch-1)*num_batches+itr]),
                (:Error, mval[(epoch-1)*num_batches+itr]),
            ],
        )

    end
    global Z0 = Z0 |> cpu
    save_dict =
        @strdict sim_name idx sigma lr max_epoch nsrc nrec src_enc fval mval Z0 J epoch opt pretrained_model_config
    @tagsave(
        datadir(sim_name, savename(save_dict, "jld2"; digits = 6)),
        save_dict;
        safe = true
    )
    global Z0 = Z0 |> gpu
end

X_hat = CHrev.forward(Z0, Zy)[1]

X_OOD = X_OOD |> cpu
X_hat = X_hat |> cpu
X0 = X0 |> cpu
Y_OOD = Y_OOD |> cpu
Zy = Zy |> cpu
Z0 = Z0 |> cpu

X_OOD = wavelet_unsqueeze(X_OOD)
X_hat = wavelet_unsqueeze(X_hat)
X0 = wavelet_unsqueeze(X0)
Y_OOD = wavelet_unsqueeze(Y_OOD)
Zy = wavelet_unsqueeze(Zy)
Z0 = wavelet_unsqueeze(Z0)

X_OOD = AN_x.inverse(X_OOD)
X_hat = AN_x.inverse(X_hat)
X0 = AN_x.inverse(X0)
Y_OOD = AN_y.inverse(Y_OOD)

Y_OOD[1:10, :, :, :] .= 0.0f0
X_hat[1:10, :, :, :] .= 0.0f0
X0[1:10, :, :, :] .= 0.0f0

# Saving parameters and logs
epoch = max_epoch
save_dict =
    @strdict sim_name idx sigma lr max_epoch nsrc nrec src_enc fval mval X_hat Z0 pretrained_model_args pretrained_model_config X_OOD Y_OOD X0 epoch
@tagsave(datadir(sim_name, savename(save_dict, "jld2"; digits = 6)), save_dict; safe = true)


text:
#These types must be the same as in the QDLDL header
c_int   = Int64;
c_float = Float64;
c_bool  = Bool;

function LDLsolve(A::SparseMatrixCSC{c_float,c_int},b::Array{c_float})

    LDLlib = "../../build/out/libqdldl.dylib"

    work   = zeros(c_int,A.n)
    Lnz    = zeros(c_int,A.n)
    etree  = zeros(c_int,A.n)

    #Julia matrices are 1-index.  C code is 0 indexed.
    Ap = A.colptr-1;
    Ai = A.rowval-1;
    Ax = A.nzval;

    #Construct the elimination tree and column counts for A
    sumLnz = ccall((:QDLDL_etree,LDLlib),c_int,(c_int,Ptr{c_int},Ptr{c_int},Ptr{c_int},Ptr{c_int},Ptr{c_int}),A.n,Ap,Ai,work,Lnz,etree)

    if(sumLnz < 0)
        return "Abort on etree";
    end

    #allocate target matrix data for L
    Ln = A.n;
    Lp = zeros(c_int,A.n+1);
    Li = zeros(c_int,sumLnz);
    Lx = zeros(c_float,sumLnz);

    #allocate D and its inverse
    D      = zeros(c_float,A.n)
    Dinv   = zeros(c_float,A.n)

    #scratch workspace values
    fwork            = zeros(c_float,Ln)
    iwork            = zeros(c_int,Ln*3)
    bwork            = zeros(c_bool,Ln)

    #call the C factorization function
    posCountD = ccall((:QDLDL_factor,LDLlib),c_int,
         (     c_int,              #n
               Ptr{c_int},         #Ap
               Ptr{c_int},         #Ai
               Ptr{c_float},       #Ax
               Ptr{c_int},         #Lp
               Ptr{c_int},         #Li
               Ptr{c_float},       #Lx
               Ptr{c_float},       #D
               Ptr{c_float},       #Dinv
               Ptr{c_int},         #Lnz
               Ptr{c_int},         #etree
               Ptr{c_bool},        #bwork
               Ptr{c_int},         #iwork
               Ptr{c_float}),      #fwork
          A.n,Ap,Ai,Ax,Lp,Li,Lx,D,Dinv,Lnz,etree,bwork,iwork,fwork)

    if(posCountD < 0)
        return "Abort on factor";
    end

    #Julia indexing is 1 based
    L = SparseMatrixCSC(Ln,Ln,Lp+1,Li+1,Lx)

    #call the C solve function
    x = copy(b)
    ccall((:QDLDL_solve,LDLlib),Void,
         (     c_int,              #n
               Ptr{c_int},         #Lp
               Ptr{c_int},         #Li
               Ptr{c_float},       #Lx
               Ptr{c_float},       #Dinv
               Ptr{c_float}),      #x (starts as b)
            Ln,Lp,Li,Lx,Dinv,x)

    return x
end

function qdMatExample()

    i = [1, 7, 2, 3, 6, 10, 2, 3, 10, 4, 8, 5, 2, 6, 1, 7, 9, 4, 8, 7, 9, 2, 3, 10]
    j = [1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 5, 6, 6, 7, 7, 7, 8, 8, 9, 9, 10, 10, 10]
    s = [1.000000e+00, 3.503210e-01, 4.606408e-01, -1.211894e-01, -2.900576e-02, 1.824522e-01, -1.211894e-01, 4.179283e-01, -1.565056e+00, 1.778279e-01, -8.453948e-02, 1.000000e-01, -2.900576e-02, -1.000000e+00, 3.503210e-01, -4.410924e-01, 1.786632e-01, -8.453948e-02, -3.162278e-01, 1.786632e-01, -2.990768e-01, 1.824522e-01, -1.565056e+00, -1.000000e-01]

    A = sparse(i,j,s)
end

function qdRandom(m,n)

    density = 0.1
    A =  pdRandom(m,density)
    C = -pdRandom(n,density)
    B = sprandn(m,n,density)
    H = [A B;B' C]

end

function pdRandom(n,density)

    X = sprandn(n,n,density)
    X = 0.5*(X+X') + n*speye(n)

end

println("Things that should fail")
println("----------------------------")
A = sparse([0.0 1.0; 1.0 0.0])
b = randn(c_float,A.n);
println("non LDL decomposable   : ",  LDLsolve(triu(A),b))

A = sparse([1.0 1.0; 1.0 1.0])
b = randn(c_float,A.n);
println("Tril structure         : ",  LDLsolve(tril(A),b))

A = sparse([5.0 1.0; 1.0 5.0])
b = randn(c_float,A.n);
println("TRIL and TRIL parts    : ",  LDLsolve(A,b))

A = sparse([0.0 5.0; 5.0 1.0])
b = randn(c_float,A.n);
println("Zero at A[0,0]         : ",  LDLsolve(triu(A),b))

A = sparse([1.0 1.0; 1.0 1.0])
b = randn(c_float,A.n);
println("Rank deficient matrix  : ",  LDLsolve(triu(A),b))

A = qdMatExample();
A[5,5] = 0.0;
b = randn(c_float,A.n);
println("Zero with allocation   : ",  LDLsolve(triu(A),b))

println("\n\n")
println("Things that should work   ")
println("(Errors are relative to A\\b)")
println("-----------------------------")

A = qdMatExample();
b = randn(c_float,A.n);
println("Basic example                : ",  norm(A\b-LDLsolve(triu(A),b)))

A = sparse([4.0 1.0 2.0; 1.0 0.0 1.0;2.0 1.0 -3.0])
b = [6.0,9.0,12.0];
println("Zero mid matrix              : ",  norm(A\b-LDLsolve(triu(A),b)))

A = sparse([1.0 5.0; 5.0 0.0])
b = randn(c_float,A.n);
println("Zero at A[end,end]           : ",  norm(A\b-LDLsolve(triu(A),b)))

A = sparse([1.0 1.0;1.0 -1.0]);
b = cumsum(ones(c_float(A.n)));
println("2x2                          : ", norm(A\b- LDLsolve(triu(A),b)))

A = speye(100);
b = ones(c_float(A.n));
println("Big identity matrix          : ", norm(A\b- LDLsolve(triu(A),b)))

A = qdMatExample();
b = randn(c_float,A.n);
println("example QD                   : ", norm(A\b- LDLsolve(triu(A),b)))

A = speye(1)/5.0;
b = randn(c_float,A.n);
println("singleton                    : ", norm(A\b- LDLsolve(triu(A),b)))

A = sparse([1.0 1.0;1.0 (1.0-1e-8)]);
b = randn(c_float,A.n);
println("ill conditioned              : ", norm(A\b- LDLsolve(triu(A),b)))

A = sparse([1.0 1.0;1.0 -(1.0 -1e-8)]);
b = randn(c_float,A.n);
println("ill conditioned QD           : ", norm(A\b- LDLsolve(triu(A),b)))

A = sparse([1e-8 1.0;1.0 -(1.0 -1e-8)]);
b = randn(c_float,A.n);
println("+Epsilon at A[0,0]           : ", norm(A\b- LDLsolve(triu(A),b)))

A = qdRandom(100,90);
b = randn(c_float,A.n);
println("Big QD                       : ", norm(A\b- LDLsolve(triu(A),b)))

A = qdRandom(0,90);
b = randn(c_float,A.n);
println("Big ND                       : ", norm(A\b- LDLsolve(triu(A),b)))

A = qdRandom(100,0);
b = randn(c_float,A.n);
println("Big PD                       : ", norm(A\b- LDLsolve(triu(A),b)))

A = sparse([1.0 5.0; 5.0 1.0])
b = randn(c_float,A.n);
println("indefinite matrix            : ", norm(A\b- LDLsolve(triu(A),b)))

Ai = [0; 1; 2; 1; 0; 3; 4; 5; 5; 6; 4; 3] + 1
Aj = [0; 1; 2; 2; 2; 3; 4; 5; 6; 6; 6; 6] + 1
Ax = [-0.25; -0.25; 1; 0.513578; 0.529142; -0.25; -0.25; 1.10274; 0.15538; 1.25883; 0.13458; 0.621134]
A = sparse(Ai, Aj, Ax)
b = randn(c_float,A.n);
println("tricky permuted osqp matrix  : ", norm((A + A' - diagm(diag(A)))\b- LDLsolve(A, b)))

max_stars_repo_name: ANYbotics/qdldl

max_issues_repo_name: ANYbotics/qdldl

max_forks_repo_name: ANYbotics/qdldl


text:
export move_agent!, add_agent!, add_agent_single!,
move_agent_single!, kill_agent!, coord2vertex, vertex2coord

"""
    kill_agent!(agent::AbstractAgent, model::ABM)

Remove an agent from the list of agents and from the space.
"""
function kill_agent!(agent::AbstractAgent, model::ABM)
  if typeof(agent.pos) <: Tuple
    agentnode = coord2vertex(agent.pos, model)
  else
    agentnode = agent.pos
  end
   # remove from the space
  splice!(agent_positions(model)[agentnode],
          findfirst(a->a==agent.id, agent_positions(model)[agentnode]))
  delete!(model.agents, agent.id)
  return
end

"""
    move_agent!(agent::AbstractAgent [, pos], model::ABM) → agent

Add `agentID` to the new position `pos` (either tuple or integer) in the model
and remove it from the old position
(also update the agent to have the new position).

If `pos` is a tuple, it represents the coordinates of the grid node.
If `pos` is an integer, it represents the node number in the graph.
If `pos` is not given, the agent is moved to a random position on the grid.
"""
function move_agent!(agent::AbstractAgent, pos::Tuple, model::ABM)
  nodenumber = coord2vertex(pos, model)
  move_agent!(agent, nodenumber, model)
end

function move_agent!(agent::AbstractAgent, pos::Integer, model::ABM)
  push!(model.space.agent_positions[pos], agent.id)
  # remove agent from old position
  if typeof(agent.pos) <: Tuple
    oldnode = coord2vertex(agent.pos, model)
    splice!(model.space.agent_positions[oldnode], findfirst(a->a==agent.id, model.space.agent_positions[oldnode]))
    agent.pos = vertex2coord(pos, model)  # update agent position
  else
    splice!(model.space.agent_positions[agent.pos], findfirst(a->a==agent.id, model.space.agent_positions[agent.pos]))
    agent.pos = pos
  end
  return agent
end

function move_agent!(agent::AbstractAgent, model::ABM)
  nodenumber = rand(1:nv(model.space))
  move_agent!(agent, nodenumber, model)
  return agent
end

"""
    move_agent_single!(agent::AbstractAgent, model::ABM)

Move agent to a random nodes on the grid while respecting a maximum of one agent
per node. If there are no empty nodes, the agent wont move.

Return the agent's new position.
"""
function move_agent_single!(agent::AbstractAgent, model::ABM)
  # TODO: this inefficient
  empty_cells = [i for i in 1:length(model.space.agent_positions) if length(model.space.agent_positions[i]) == 0]
  if length(empty_cells) > 0
    random_node = rand(empty_cells)
    move_agent!(agent, random_node, model)
  end
  return agent
end

"""
    add_agent!(agent::AbstractAgent [, pos], model::ABM) → agent

Adds the agent to the `pos` in the space and to the list of agents.
If `pos` is not given, the agent is added to a random position.
The agent's position is always updated to match `pos`.
"""
function add_agent!(agent::AbstractAgent, pos::Tuple, model::ABM)
  # node number from x, y, z coordinates
  nodenumber = coord2vertex(pos, model)
  add_agent!(agent, nodenumber, model)
end

function add_agent!(agent::AbstractAgent, pos::Integer, model::ABM)
  push!(model.space.agent_positions[pos], agent.id)
  model.agents[agent.id] = agent
  # update agent position
  if typeof(agent.pos) <: Integer
    agent.pos = pos
  elseif typeof(agent.pos) <: Tuple
    agent.pos = vertex2coord(pos, model)
  else
    throw("Unknown type of agent.pos.")
  end
  return agent
end

function add_agent!(agent::AbstractAgent, model::ABM)
  nodenumber = rand(1:nv(model.space))
  add_agent!(agent, nodenumber, model)
  return agent
end

"""
    add_agent_single!(agent::AbstractAgent, model::ABM) → agent

Add agent to a random node in the space while respecting a maximum one agent per node.
This function does not do anything if there are no empty nodes.

Return the new agent.
"""
function add_agent_single!(agent::AbstractAgent, model::ABM)
  msa = model.space.agent_positions
  empty_cells = [i for i in 1:length(msa) if length(msa[i]) == 0]
  if length(empty_cells) > 0
    random_node = rand(empty_cells)
    add_agent!(agent, random_node, model)
  end
  return agent
end

function biggest_id(model)
    if isempty(model.agents)
        return 0
    else
        return maximum(keys(model.agents))
    end
end

"""
    add_agent!(node, model::ABM, properties...)
Add a new agent in the given `node`, by constructing the agent type of
the `model` and propagating all extra `properties` to the constructor.
"""
function add_agent!(node, model::ABM, properties...)
    id = biggest_id(model) + 1
    A = agenttype(model)
    cnode = correct_pos_type(node, model)
    agent = A(id, cnode, properties...)
    add_agent!(agent, cnode, model)
end

"""
    add_agent!(model::ABM, properties...)
Similar with `add_agent!(node, model, properties...)`, but adds the
created agent to a random node.
This function also works for models without a spatial structure.
"""
function add_agent!(model::ABM{A, Nothing}, properties...) where {A}
  @assert model.space == nothing
  id = biggest_id(model) + 1
  model.agents[id] = A(id, properties...)
  return model.agents[id]
end

function add_agent!(model::ABM{A, S}, properties...) where {A, S<:AbstractSpace}
  id = biggest_id(model) + 1
  n = rand(1:nv(model))
  cnode = correct_pos_type(n, model)
  model.agents[id] = A(id, cnode, properties...)
  return model.agents[id]
end

"""
    add_agent_single!(model::ABM, properties...)
Same as `add_agent!(model, properties...)` but ensures that it adds an agent
into a node with no other agents (does nothing if no such node exists).
"""
function add_agent_single!(model::ABM, properties...)
  msa = model.space.agent_positions
  id = biggest_id(model) + 1
  A = agenttype(model)
  empty_cells = [i for i in 1:length(msa) if length(msa[i]) == 0]
  if length(empty_cells) > 0
    random_node = rand(empty_cells)
    cnode = correct_pos_type(node, model)
    agent = A(id, cnode, properties...)
    add_agent!(agent, random_node, model)
    return agent
  end
end


text:
"""
    dct2!(a, fft_type::FFTType=Radix2Type())

Discrete cosine transform of type II

Y = DCT2(X) returns the discrete cosine transform of X,
based on the staggered-grid definition
           N
   Y(k) =  ∑  X(j) cos (π*(k-1)*(j-1/2)/N)
          j=1
The vector Y is the same size as X and contains the
discrete cosine transform coefficients.

Y = DCT2(X,N) pads or truncates the vector X to length N
before transforming.

If X is a matrix, the DCT2 operation is applied to each
column. This transform can be inverted using IDCT2.

Verified correct within machine precision for 2D array input 2/28/2021

Traditional dct2:
          N-1
   Y(k) =  ∑  X(j) cos (π*(j+1/2)*k/N)
          j=0

"""
function dct2!(x, fft_type=Radix2Type())
    x .= dct2(x, fft_type)
end
function dct2(a, fft_type::FFTType=Radix2Type())
    FT = eltype(a)
    if min(size(a)...)==1
        if size(a,2)>1
            do_trans = true;
        else
            do_trans = false;
        end
        a = a[:];
    else
        do_trans = false;
    end
    n = size(a,1);
    m = size(a,2);
    # Pad or truncate a if necessary
    if size(a,1)<n
      aa = zeros(n,m);
      aa[1:size(a,1),:] = a;
    else
      aa = a[1:n,:];
    end
    y = zeros(2*n, m);
    y[1:n,:] .= aa;
    i = Complex(0,1);
    e = 0.5 * exp.(-i*FT(0.5) * π * (0:n-1) ./ n);
    y[n+1:2*n,:] .= reverse(aa[1:n,:]; dims=1);
    yy = fft(y, 1);
    b = zeros(n,m);
    for l = 1:m
        b[1:n,l] .= real.(yy[1:n,l] .* e);
    end
    if isreal(a)
        b = real.(b);
    end
    if do_trans
        b = b'
    end
    return b
end


text:
abstract type TrustRegionUpdater end
struct TrustRegion{M, SP, D}
	scheme::M
	spsolve::SP
	Δupdate::D
end
summary(tr::TrustRegion) = "$(summary(modelscheme(tr))) with $(summary(algorithm(tr)))"
function initial_preconditioner(approach::TrustRegion, x)
  nothing
end
"""
  BTR() <: TrustRegionUpdater

Basic trust region updater following, and named after [CGT].
"""
struct BTR{T}
  Δmin::T
end
TrustRegion(; deltamin=nothing) = TrustRegion(Newton(), NTR(), BTR(deltamin))
TrustRegion(m, sp=NTR(); deltamin=nothing) = TrustRegion(m, sp, BTR(deltamin))
modelscheme(tr::TrustRegion) = tr.scheme
algorithm(tr::TrustRegion) = tr.spsolve

solve(problem::OptimizationProblem, x0, scheme, options::OptimizationOptions) =
  solve(problem, (x0, nothing), TrustRegion(scheme, NWI()), options)

solve(problem::OptimizationProblem{<:Any, <:Nothing, <:Nothing, <:Nothing}, x0, approach::TrustRegion, options::OptimizationOptions) =
  solve(problem, (x0, nothing), approach, options)

solve(problem::OptimizationProblem{<:Any, <:Nothing, <:Nothing, <:Nothing}, s0::Tuple, approach::TrustRegion, options::OptimizationOptions) =
  solve(problem, s0, approach, options)

function solve(problem::OptimizationProblem, x0, scheme::Newton, options::OptimizationOptions)
    solve(problem, (x0, nothing), TrustRegion(scheme, NTR()), options)
end
function solve(problem::OptimizationProblem, x0, approach::TrustRegion, options::OptimizationOptions)
    solve(problem, (x0, nothing), approach, options)
end
include("optimize/inplace_loop.jl")
include("optimize/outofplace_loop.jl")


text:
__precompile__()

module XYZ

# export
#     # cloud
#     Cloud, BoundingBox, positions, attributes, boundingbox, calc_bbox, update_bbox!,
#     # pointfilters
#     ground, water, notground, notoutlier, lastreturn,
#     # classifiers
#     classify_ground!, classify_ground_approx!, classify_outliers!, classify_water!,
#     classify_min!, classify_max!, classify!, reset_class!,
#     classify_below_surface!, classify_above_surface!, classify_water_approx!,
#     # raster
#     Raster, rasterize, define_raster,
#     # profile
#     define_profile, filter_profile,
#     # reducers
#     reduce_min, reduce_max, reduce_pointfilter,
#     # writers
#     to_tif, to_las, to_xyz, grid2tif,
#     # readers
#     read_pointcloud,
#     # utils
#     describe, paint!, create_kdtree, point_density!, grid2cloud_attribute!, copy_classification!

using Compat
using Nullables
using ColorTypes
using FixedPointNumbers
using Colors
using GDAL
using FileIO
using LasIO
using NearestNeighbors
using PeatUtils
using ProgressMeter
using SortingAlgorithms
using StaticArrays

include("cloud.jl")
include("pointfilters.jl")
include("rasterize.jl")
include("cloudutils.jl")
include("profiler.jl")
include("cloudfilter.jl")
include("cloudclassifier.jl")
include("lasio.jl")
include("writers.jl")

end


text:
using BlossomMatching
using Test

@testset "BlossomMatching.jl" begin
    # Write your own tests here.
end


text:
mdp = GridWorld()

sim(mdp, max_steps=100) do state
    @assert isa(state, GridWorldState)
    acts = actions(mdp)
    return rand(acts)
end

pomdp = BabyPOMDP()

sim(pomdp, max_steps=100) do obs
    @assert isa(obs, Bool)
    acts = actions(pomdp)
    return rand(acts)
end

sim(pomdp, false, max_steps=100) do obs
    @assert isa(obs, Bool)
    acts = actions(pomdp)
    return rand(acts)
end

sim(pomdp, initial_state=true, max_steps=100) do obs
    @assert isa(obs, Bool)
    acts = actions(pomdp)
    return rand(acts)
end


text:
#
# The filter combinator.
#

function FilterQuery(base::Query, pred::Query)
    @assert fits(output(base), input(pred))
    @assert fits(output(pred), Domain(Bool))
    @assert !isplural(output(pred))
    dom = domain(output(base))
    return (
        base >>
        RecordQuery(ItQuery(dom), pred) >>
        Query(
            FilterSig(),
            Input([dom, output(pred)]),
            Output(dom) |> setoptional()))
end

immutable FilterSig <: AbstractPrimitive
end

describe(io::IO, ::FilterSig) = print(io, "filter′")

ev(sig::FilterSig, dv::DataVector) =
    isplain(column(dv, 2)) ?
        plain_filter_impl(values(dv, 1), values(dv, 2)) :
        filter_impl(values(dv, 1), offsets(dv, 2), values(dv, 2))

function plain_filter_impl(vals::AbstractVector, predvals::AbstractVector{Bool})
    len = length(vals)
    size = 0
    for pred in predvals
        if pred
            size += 1
        end
    end
    if size == len
        return OptionalColumn(OneTo(len+1), vals)
    end
    offs = Vector{Int}(len+1)
    offs[1] = 1
    idxs = Vector{Int}(size)
    n = 1
    for k in eachindex(predvals)
        if predvals[k]
            idxs[n] = k
            n += 1
        end
        offs[k+1] = n
    end
    return OptionalColumn(offs, vals[idxs])
end

function filter_impl(vals::AbstractVector, predoffs::AbstractVector{Int}, predvals::AbstractVector{Bool})
    len = length(vals)
    size = 0
    for pred in predvals
        if pred
            size += 1
        end
    end
    if size == len
        return OptionalColumn(OneTo(len+1), vals)
    end
    offs = Vector{Int}(len+1)
    offs[1] = 1
    idxs = Vector{Int}(size)
    n = 1
    for k in 1:len
        l = predoffs[k]
        r = predoffs[k+1]
        if l < r && predvals[l]
            idxs[n] = k
            n += 1
        end
        offs[k+1] = n
    end
    return OptionalColumn(offs, vals[idxs])
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

#=
- read meshes from different formats
- reorder connectivity, create element sets, node sets, ...
- create partitions for parallel runs
- renumber elements / nodes
- maybe precheck for bad elements
- check surface normal direction in boundary elements
- orientation of 2d elements
- etc only topology related stuff
=#

importall Base

using JuliaFEM

type Mesh
    nodes :: Dict{Int64, Vector{Float64}}
    node_sets :: Dict{Symbol, Set{Int64}}
    elements :: Dict{Int64, Vector{Int64}}
    element_types :: Dict{Int64, Symbol}
    element_codes :: Dict{Int64, Symbol}
    element_sets :: Dict{Symbol, Set{Int64}}
    surfaces :: Dict{Symbol, Vector{Tuple{Int64, Symbol}}}
    surface_types :: Dict{Symbol, Symbol}
end

function Mesh()
    return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict())
end

function add_node!(mesh::Mesh, nid::Int, ncoords::Vector{Float64})
    mesh.nodes[nid] = ncoords
end

function add_nodes!(mesh::Mesh, nodes::Dict{Int64, Vector{Float64}})
    for (nid, ncoords) in nodes
        add_node!(mesh, nid, ncoords)
    end
end

function add_node_to_node_set!(mesh::Mesh, set_name, nids...)
    if !haskey(mesh.node_sets, set_name)
        mesh.node_sets[set_name] = Set{Int64}()
    end
    push!(mesh.node_sets[set_name], nids...)
end

function add_element!(mesh::Mesh, elid::Int, eltype::Symbol, connectivity::Vector{Int64})
    mesh.elements[elid] = connectivity
    mesh.element_types[elid] = eltype
end

function add_elements!(mesh::Mesh, elements::Dict{Int64, Tuple{Symbol, Vector{Int64}}})
    for (elid, (eltype, elcon)) in elements
        add_element!(mesh, elid, eltype, elcon)
    end
end

function add_element_to_element_set!(mesh::Mesh, set_name, elids...)
    if !haskey(mesh.element_sets, set_name)
        mesh.element_sets[set_name] = Set{Int64}()
    end
    push!(mesh.element_sets[set_name], elids...)
end

function copy(mesh::Mesh)
    mesh2 = Mesh()
    mesh2.nodes = copy(mesh.nodes)
    mesh2.node_sets = copy(mesh.node_sets)
    mesh2.elements = copy(mesh.elements)
    mesh2.element_types = copy(mesh.element_types)
    mesh2.element_sets = copy(mesh.element_sets)
    return mesh2
end

function filter_by_element_id(mesh::Mesh, element_ids::Vector{Int64})
    mesh2 = copy(mesh)
    mesh2.elements = Dict()
    for elid in element_ids
        if haskey(mesh.elements, elid)
            mesh2.elements[elid] = mesh.elements[elid]
        end
    end
    return mesh2
end

function filter_by_element_set(mesh::Mesh, set_name)
    filter_by_element_id(mesh::Mesh, collect(mesh.element_sets[set_name]))
end

function create_element(mesh::Mesh, id::Int)
    connectivity = mesh.elements[id]
    element_type = JuliaFEM.(mesh.element_types[id])
    element = Element(element_type, connectivity)
    update!(element, "geometry", mesh.nodes)
    element.id = id
    return element
end

function create_elements(mesh::Mesh; element_type=nothing)
    element_ids = collect(keys(mesh.elements))
    if element_type != nothing
        filter!(id -> mesh.element_types[id] == element_type, element_ids)
    end
    elements = [create_element(mesh, id) for id in element_ids]
    return elements
end

function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothing)
    if isempty(element_sets)
        element_ids = collect(keys(mesh.elements))
    else
        element_ids = Set{Int64}()
        for set_name in element_sets
            element_ids = union(element_ids, mesh.element_sets[set_name])
        end
    end

    if element_type != nothing
        filter!(id -> mesh.element_types[id] == element_type, element_ids)
    end

    elements = [create_element(mesh, id) for id in element_ids]
    return elements
end

function create_elements(mesh::Mesh, element_sets::AbstractString...; element_type=nothing)
    element_sets = map(parse, element_sets)
    return create_elements(mesh, element_sets...; element_type=element_type)
end

""" find npts nearest nodes from mesh and return id numbers as list. """
function find_nearest_nodes(mesh::Mesh, coords::Vector, npts=1)
    dist = Dict{Int64, Float64}()
    for (nid, c) in mesh.nodes
        dist[nid] = norm(coords-c)
    end
    s = sort(collect(dist), by=x->x[2])
    nd = s[1:npts] # [(id1, dist1), (id2, dist2), ..., (id_npts, dist_npts)]
    node_ids = [n[1] for n in nd]
    return node_ids
end

"""
Apply new node ordering to elements. In JuliaFEM same node ordering is used
than in ABAQUS and if mesh is parsed from FEM format with other node ordering
this can be used to reorder nodes.

Parameters
----------
mapping :: Dict{Symbol, Vector{Int}}
    e.g. :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

"""
function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol, Vector{Int}})
    for (elid, eltype) in mesh.element_types
        haskey(mapping, eltype) || continue
        new_order = mapping[eltype]
        element_connectivity = mesh.elements[elid]
        new_element_connectivity = element_connectivity[new_order]
        mesh.elements[elid] = new_element_connectivity
    end
end

function JuliaFEM.Problem{P<:FieldProblem}(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int64)
    problem = Problem{P}(name, dimension, "none", [], Dict(), Assembly(), P())
    problem.elements = create_elements(mesh, name)
    return problem
end

function JuliaFEM.Problem{P<:BoundaryProblem}(mesh::Mesh, ::Type{P}, name, dimension, parent_field_name)
    problem = Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
    problem.elements = create_elements(mesh, name)
    return problem
end

"""
Swap surface element connectivity s.t. normals point outward
"""
function check_orientation!
    # TODO
end

"""
Partition model using METIS
"""
function partition_model!
    # TODO
end


text:
tutorial_path(path="") = joinpath("src", path)

function scan_path!(c, path; ext=".md") # dispatches based on path type
    tut_path = tutorial_path(path)
    if isdir(tut_path)
        add_dir!(c, path; ext=ext)
    elseif isfile(tut_path)
        add_page!(c, path; ext=ext)
    else
        error("$(path) is neither a file nor a directory")
    end
    if length(c) == 1
        return c[1]
    else
        return c
    end
end

function add_page!(c, path; ext)
    tut_path = tutorial_path(path)
    if endswith(tut_path, ext)
        title = open(tut_path) do f
            firstline = readline(f)
            strip(firstline, ['#', ' '])
        end
        push!(c, string(title) => string(path))
    end
    return c
end

function add_dir!(c, path; ext)
    new_c = []
    for p in readdir(tutorial_path(path))
        scan_path!(new_c, joinpath(path, p); ext=ext)
    end
    if length(new_c) == 1
        push!(c, new_c[1])
    elseif length(new_c) > 1
        title = basename(path)
        push!(c, string(title) => new_c)
    end
    return c
end

tutorials = scan_path!([], "Tutorials")


text:
# MIT License: Copyright (c) 2016: Andy Ferris.
# See LICENSE.md for further licensing test

@inline function LinearAlgebra.eigen(S::SymmetricTensor{2, 1, T}) where {T}
    @inbounds Eigen(Vec{1, T}((S[1, 1],)), one(Tensor{2, 1, T}))
end

function LinearAlgebra.eigen(S::SymmetricTensor{2, 2, T}) where {T}
    @inbounds begin
        if S[2,1] == 0 # diagonal tensor
            S11 = S[1,1]
            S22 = S[2,2]
            if S11 < S22
                λ = Vec{2}((S11, S22))
                Φ = Tensor{2,2,T}((T(1), T(0), T(0), T(1)))
            else # S22 <= S11
                λ = Vec{2}((S22, S11))
                Φ = Tensor{2,2,T}((T(0), T(1), T(1), T(0)))
            end
        else
            # eigenvalues from quadratic formula
            trS_half = tr(S) / 2
            tmp2 = trS_half * trS_half - det(S)
            tmp2 < 0 ? tmp = zero(tmp2) : tmp = sqrt(tmp2) # Numerically stable for identity matrices, etc.
            λ = Vec{2}((trS_half - tmp, trS_half + tmp))

            Φ11 = λ[1] - S[2,2]
            n1 = sqrt(Φ11 * Φ11 + S[2,1] * S[2,1])
            Φ11 = Φ11 / n1
            Φ12 = S[2,1] / n1

            Φ21 = λ[2] - S[2,2]
            n2 = sqrt(Φ21 * Φ21 + S[2,1] * S[2,1])
            Φ21 = Φ21 / n2
            Φ22 = S[2,1] / n2

            Φ = Tensor{2, 2}((Φ11, Φ12,
                              Φ21, Φ22))
        end
        return Eigen(λ, Φ)
    end
end

# A small part of the code in the following method was inspired by works of David
# Eberly, Geometric Tools LLC, in code released under the Boost Software
# License. See LICENSE.md
function LinearAlgebra.eigen(S::SymmetricTensor{2, 3, T}) where {T}
    @inbounds begin
        R = typeof((one(T)*zero(T) + zero(T))/one(T))
        SR = convert(SymmetricTensor{2, 3, R}, S)

        S11 = SR[1, 1]; S22 = SR[2, 2]; S33 = SR[3, 3]
        S12 = SR[1, 2]; S13 = SR[1, 3]; S23 = SR[2, 3]

        p1 = abs2(S12) + abs2(S13) + abs2(S23)
        if (p1 == 0) # diagonal tensor
            v1, v2, v3 = basevec(Vec{3, R})
            if S11 < S22
                if S22 < S33
                    return Eigen(Vec{3, R}((S11, S22, S33)), Tensor{2, 3, R}((v1[1], v1[2], v1[3], v2[1], v2[2], v2[3], v3[1], v3[2], v3[3])))
                elseif S33 < S11
                    return Eigen(Vec{3, R}((S33, S11, S22)), Tensor{2, 3, R}((v3[1], v3[2], v3[3], v1[1], v1[2], v1[3], v2[1], v2[2], v2[3])))
                else
                    return Eigen(Vec{3, R}((S11, S33, S22)), Tensor{2, 3, R}((v1[1], v1[2], v1[3], v3[1], v3[2], v3[3], v2[1], v2[2], v2[3])))
                end
            else #S22 < S11
                if S11 < S33
                    return Eigen(Vec{3, R}((S22, S11, S33)), Tensor{2, 3, R}((v2[1], v2[2], v2[3], v1[1], v1[2], v1[3], v3[1], v3[2], v3[3])))
                elseif S33 < S22
                    return Eigen(Vec{3, R}((S33, S22, S11)), Tensor{2, 3, R}((v3[1], v3[2], v3[3], v2[1], v2[2], v2[3], v1[1], v1[2], v1[3])))
                else
                    return Eigen(Vec{3, R}((S22, S33, S11)), Tensor{2, 3, R}((v2[1], v2[2], v2[3], v3[1], v3[2], v3[3], v1[1], v1[2], v1[3])))
                end
            end
        end

        q = (S11 + S22 + S33) / 3
        p2 = abs2(S11 - q) + abs2(S22 - q) + abs2(S33 - q) + 2 * p1
        p = sqrt(p2 / 6)
        invp = inv(p)
        b11 = (S11 - q) * invp
        b22 = (S22 - q) * invp
        b33 = (S33 - q) * invp
        b12 = S12 * invp
        b13 = S13 * invp
        b23 = S23 * invp
        B = SymmetricTensor{2, 3, R}((b11, b12, b13, b22, b23, b33))
        r = det(B) / 2

        # In exact arithmetic for a symmetric matrix -1 <= r <= 1
        # but computation error can leave it slightly outside this range.
        if (r <= -1)
            phi = R(pi) / 3
        elseif (r >= 1)
            phi = zero(R)
        else
            phi = acos(r) / 3
        end

        λ3 = q + 2 * p * cos(phi)
        λ1 = q + 2 * p * cos(phi + (2*R(pi)/3))
        λ2 = 3 * q - λ1 - λ3 # since tr(S) = λ1 + λ2 + λ3

        if r > 0 # Helps with conditioning the eigenvector calculation
            (λ1, λ3) = (λ3, λ1)
        end

        # Calculate the first eigenvector
        # This should be orthogonal to these three rows of A - λ1*I
        # Use all combinations of cross products and choose the "best" one
        r₁ = Vec{3, R}((S11 - λ1, S12, S13))
        r₂ = Vec{3, R}((S12, S22 - λ1, S23))
        r₃ = Vec{3, R}((S13, S23, S33 - λ1))
        n₁ = r₁ ⋅ r₁
        n₂ = r₂ ⋅ r₂
        n₃ = r₃ ⋅ r₃

        r₁₂ = r₁ × r₂
        r₂₃ = r₂ × r₃
        r₃₁ = r₃ × r₁
        n₁₂ = r₁₂ ⋅ r₁₂
        n₂₃ = r₂₃ ⋅ r₂₃
        n₃₁ = r₃₁ ⋅ r₃₁

        # we want best angle so we put all norms on same footing
        # (cheaper to multiply by third nᵢ rather than divide by the two involved)
        if n₁₂ * n₃ > n₂₃ * n₁
            if n₁₂ * n₃ > n₃₁ * n₂
                Φ1 = r₁₂ / sqrt(n₁₂)
            else
                Φ1 = r₃₁ / sqrt(n₃₁)
            end
        else
            if n₂₃ * n₁ > n₃₁ * n₂
                Φ1 = r₂₃ / sqrt(n₂₃)
            else
                Φ1 = r₃₁ / sqrt(n₃₁)
            end
        end

        # Calculate the second eigenvector
        # This should be orthogonal to the previous eigenvector and the three
        # rows of A - λ2*I. However, we need to "solve" the remaining 2x2 subspace
        # problem in case the cross products are identically or nearly zero

        # The remaing 2x2 subspace is:
        if abs(Φ1[1]) < abs(Φ1[2]) # safe to set one component to zero, depending on this
            orthogonal1 = Vec{3, R}((-Φ1[3], zero(R), Φ1[1])) / sqrt(abs2(Φ1[1]) + abs2(Φ1[3]))
        else
            orthogonal1 = Vec{3, R}((zero(R), Φ1[3], -Φ1[2])) / sqrt(abs2(Φ1[2]) + abs2(Φ1[3]))
        end
        orthogonal2 = Φ1 × orthogonal1

        # The projected 2x2 eigenvalue problem is C x = 0 where C is the projection
        # of (A - λ2*I) onto the subspace {orthogonal1, orthogonal2}
        a_orth1_1 = S11 * orthogonal1[1] + S12 * orthogonal1[2] + S13 * orthogonal1[3]
        a_orth1_2 = S12 * orthogonal1[1] + S22 * orthogonal1[2] + S23 * orthogonal1[3]
        a_orth1_3 = S13 * orthogonal1[1] + S23 * orthogonal1[2] + S33 * orthogonal1[3]

        a_orth2_1 = S11 * orthogonal2[1] + S12 * orthogonal2[2] + S13 * orthogonal2[3]
        a_orth2_2 = S12 * orthogonal2[1] + S22 * orthogonal2[2] + S23 * orthogonal2[3]
        a_orth2_3 = S13 * orthogonal2[1] + S23 * orthogonal2[2] + S33 * orthogonal2[3]

        c11 = orthogonal1[1]*a_orth1_1 + orthogonal1[2]*a_orth1_2 + orthogonal1[3]*a_orth1_3 - λ2
        c12 = orthogonal1[1]*a_orth2_1 + orthogonal1[2]*a_orth2_2 + orthogonal1[3]*a_orth2_3
        c22 = orthogonal2[1]*a_orth2_1 + orthogonal2[2]*a_orth2_2 + orthogonal2[3]*a_orth2_3 - λ2

        # Solve this robustly (some values might be small or zero)
        c11² = abs2(c11)
        c12² = abs2(c12)
        c22² = abs2(c22)
        if c11² >= c22²
            if c11² > 0 || c12² > 0
                if c11² >= c12²
                    tmp = c12 / c11
                    p2 = inv(sqrt(1 + abs2(tmp)))
                    p1 = tmp * p2
                else
                    tmp = c11 / c12 # TODO check for compex input
                    p1 = inv(sqrt(1 + abs2(tmp)))
                    p2 = tmp * p1
                end
                Φ2 = p1*orthogonal1 - p2*orthogonal2
            else # c11 == 0 && c12 == 0 && c22 == 0 (smaller than c11)
                Φ2 = orthogonal1
            end
        else
            if c22² >= c12²
                tmp = c12 / c22
                p1 = inv(sqrt(1 + abs2(tmp)))
                p2 = tmp * p1
            else
                tmp = c22 / c12
                p2 = inv(sqrt(1 + abs2(tmp)))
                p1 = tmp * p2
            end
            Φ2 = p1*orthogonal1 - p2*orthogonal2
        end

        # The third eigenvector is a simple cross product of the other two
        Φ3 = Φ1 × Φ2 # should be normalized already

        # Sort them back to the original ordering, if necessary
        if r > 0
            (λ1, λ3) = (λ3, λ1)
            (Φ1, Φ3) = (Φ3, Φ1)
        end

        λ = Vec{3}((λ1, λ2, λ3))
        Φ = Tensor{2, 3}((Φ1[1], Φ1[2], Φ1[3],
                          Φ2[1], Φ2[2], Φ2[3],
                          Φ3[1], Φ3[2], Φ3[3]))
        return Eigen(λ, Φ)
    end
end


text:
using Random
Random.seed!()

passLength, numMatchesForLog = 8, 1
possibleChars = ['a':'z' ; 'A':'Z' ; '0':'9']

correctPassword = "3xyZu4vN"

numMatch(loginPassword) =
    sum([loginPassword[i] == correctPassword[i] for i in 1:passLength])

N = 10^7

passwords = [String(rand(possibleChars,passLength)) for _ in 1:N]
numLogs   = sum([numMatch(p) >= numMatchesForLog for p in passwords])
println("Number of login attempts logged: ", numLogs)
println("Proportion of login attempts logged: ", numLogs/N)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using DiffEqBase, DiffEqJump, DiffEqProblemLibrary, Plots, Statistics
using DiffEqProblemLibrary.JumpProblemLibrary: importjumpproblems; importjumpproblems()
import DiffEqProblemLibrary.JumpProblemLibrary: prob_jump_dnarepressor
gr()
fmt = :png

methods = (Direct(),DirectFW(),FRM(),FRMFW(),SortingDirect(),NRM(),DirectCR(),RSSA())
shortlabels = [string(leg)[12:end-2] for leg in methods]
prob    = prob_jump_dnarepressor.discrete_prob
tf      = prob_jump_dnarepressor.tstop
rn      = prob_jump_dnarepressor.network
ploth   = plot(reuse=false)
for (i,method) in enumerate(methods)
    jump_prob = JumpProblem(prob, method, rn, save_positions=(false,false))
    sol = solve(jump_prob, SSAStepper(), saveat=tf/1000.)
    plot!(ploth,sol.t,sol[3,:],label=shortlabels[i], format=fmt)
end
plot(ploth, title="Protein level", xlabel="time",format=fmt)

p = []
for (i,method) in enumerate(methods)
    jump_prob = JumpProblem(prob, method, rn, save_positions=(false,false))
    sol = solve(jump_prob, SSAStepper(), saveat=tf/1000.)
    push!(p, plot(sol,title=shortlabels[i],leg=false,format=fmt))
end
plot(p...,format=fmt)

function run_benchmark!(t, jump_prob, stepper)
    sol = solve(jump_prob, stepper)
    @inbounds for i in 1:length(t)
        t[i] = @elapsed (sol = solve(jump_prob, stepper))
    end
end

nsims = 500
benchmarks = Vector{Vector{Float64}}()
for method in methods
    jump_prob = JumpProblem(prob, method, rn, save_positions=(false,false))
    stepper = SSAStepper()
    t = Vector{Float64}(undef,nsims)
    run_benchmark!(t, jump_prob, stepper)
    push!(benchmarks, t)
end

medtimes = Vector{Float64}(undef,length(methods))
stdtimes = Vector{Float64}(undef,length(methods))
avgtimes = Vector{Float64}(undef,length(methods))
for i in 1:length(methods)
    medtimes[i] = median(benchmarks[i])
    avgtimes[i] = mean(benchmarks[i])
    stdtimes[i] = std(benchmarks[i])
end
println(medtimes/medtimes[1])

using DataFrames
df = DataFrame(names=shortlabels,medtimes=medtimes,relmedtimes=(medtimes/medtimes[1]),avgtimes=avgtimes, std=stdtimes, cv=stdtimes./avgtimes)
sa = [text(string(round(mt,sigdigits=2),"s"),:center,10) for mt in df.medtimes]
bar(df.names,df.relmedtimes,legend=:false, fmt=fmt)
scatter!(df.names, .05 .+ df.relmedtimes, markeralpha=0, series_annotations=sa, fmt=fmt)
ylabel!("median relative to Direct")
title!("Negative Feedback Gene Expression Model")

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
module mmass0
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using LinearAlgebra
using Test
function test()
    a = SysmatAssemblerSparse(0.0)
    start!(a, 7, 7)
    m = [0.24406   0.599773    0.833404  0.0420141
    0.786024  0.00206713  0.995379  0.780298
    0.845816  0.198459    0.355149  0.224996]
    gi = [1 7 5]
    gj = [5 2 1 4]
    for j in 1:size(m, 2), i in 1:size(m, 1)
        assemble!(a, gi[i], gj[j], m[i, j])
    end
    m = [0.146618  0.53471   0.614342    0.737833
    0.479719  0.41354   0.00760941  0.836455
    0.254868  0.476189  0.460794    0.00919633
    0.159064  0.261821  0.317078    0.77646
    0.643538  0.429817  0.59788     0.958909]
    gi =  [2 3 1 7 5]
    gj = [6 7 3 4]
    for j in 1:size(m, 2), i in 1:size(m, 1)
        assemble!(a, gi[i], gj[j], m[i, j])
    end
    A = finish!(a)
    B = [ 0.833404  0.599773    0.460794    0.0512104  0.24406   0.254868  0.476189
    0         0          0.614342    0.737833    0        0.146618  0.53471
    0         0          0.00760941  0.836455    0        0.479719  0.41354
    0         0           0           0          0         0         0
    0.355149  0.198459    0.59788     1.1839     0.845816  0.643538  0.429817
    0         0           0           0          0         0         0
    0.995379  0.00206713  0.317078    1.55676    0.786024  0.159064  0.261821  ]
    @test norm(A - B) / norm(B) < 1.0e-5
    true
end
end
using .mmass0
mmass0.test()

module mmass0a
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using LinearAlgebra
using Test
function test()
    a = SysmatAssemblerSparse(0.0)
    start!(a, 7, 7)
    m = [0.24406   0.599773    0.833404  0.0420141
    0.786024  0.00206713  0.995379  0.780298
    0.845816  0.198459    0.355149  0.224996]
    gi = [1 7 5]
    gj = [5 2 1 4]
    la = LocalMatrixAssembler(3, 4, 0.0)
    init!(la, gi, gj)
    for j in 1:size(m, 2), i in 1:size(m, 1)
        la[i, j] += m[i, j]
    end
    assemble!(a, la)
    m = [0.146618  0.53471   0.614342    0.737833
    0.479719  0.41354   0.00760941  0.836455
    0.254868  0.476189  0.460794    0.00919633
    0.159064  0.261821  0.317078    0.77646
    0.643538  0.429817  0.59788     0.958909]
    gi =  [2 3 1 7 5]
    gj = [6 7 3 4]
    la = LocalMatrixAssembler(5, 4, 0.0)
    init!(la, gi, gj)
    for j in 1:size(m, 2), i in 1:size(m, 1)
        la[i, j] += m[i, j]
    end
    assemble!(a, la)
    A = finish!(a)
    B = [ 0.833404  0.599773    0.460794    0.0512104  0.24406   0.254868  0.476189
    0         0          0.614342    0.737833    0        0.146618  0.53471
    0         0          0.00760941  0.836455    0        0.479719  0.41354
    0         0           0           0          0         0         0
    0.355149  0.198459    0.59788     1.1839     0.845816  0.643538  0.429817
    0         0           0           0          0         0         0
    0.995379  0.00206713  0.317078    1.55676    0.786024  0.159064  0.261821  ]
    @test norm(A - B) / norm(B) < 1.0e-5
    true
end
end
using .mmass0a
mmass0a.test()

module mmass0t
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using LinearAlgebra
using Test
function test()
    a = SysmatAssemblerSparse(0.0)
    start!(a, 7, 7)
    m = transpose([0.24406   0.599773    0.833404  0.0420141
    0.786024  0.00206713  0.995379  0.780298
    0.845816  0.198459    0.355149  0.224996])
    gj = [1 7 5]
    gi = [5 2 1 4]
    la = LocalMatrixAssembler(4, 3, 0.0)
    @test size(la) == (4, 3)
    init!(la, gi, gj)
    for j in 1:size(m, 2), i in 1:size(m, 1)
        la[i, j] += m[i, j]
    end
    assemble!(a, transpose(la))
    m = transpose([0.146618  0.53471   0.614342    0.737833
    0.479719  0.41354   0.00760941  0.836455
    0.254868  0.476189  0.460794    0.00919633
    0.159064  0.261821  0.317078    0.77646
    0.643538  0.429817  0.59788     0.958909]   )
    gj =  [2 3 1 7 5]
    gi = [6 7 3 4]
    la = LocalMatrixAssembler(4, 5, 0.0)
    init!(la, gi, gj)
    for j in 1:size(m, 2), i in 1:size(m, 1)
        la[i, j] += m[i, j]
    end
    assemble!(a, transpose(la))
    A = finish!(a)
    B = [ 0.833404  0.599773    0.460794    0.0512104  0.24406   0.254868  0.476189
    0         0          0.614342    0.737833    0        0.146618  0.53471
    0         0          0.00760941  0.836455    0        0.479719  0.41354
    0         0           0           0          0         0         0
    0.355149  0.198459    0.59788     1.1839     0.845816  0.643538  0.429817
    0         0           0           0          0         0         0
    0.995379  0.00206713  0.317078    1.55676    0.786024  0.159064  0.261821  ]
    @test norm(A - B) / norm(B) < 1.0e-5
    la = LocalVectorAssembler(4, 0.0)
    @test size(la)  == (4, )
    true
end
end
using .mmass0t
mmass0t.test()

module mvass1
using Elfel.Assemblers: SysvecAssembler, start!, finish!, assemble!
using LinearAlgebra
using Test
function test()
    N = 30
    v = fill(0.0, N)
    a = SysvecAssembler()
    start!(a,  N)
    vec, dofnums = rand(3), [1, 7, 2]
    for i in 1:length(vec)
        assemble!(a, dofnums[i], vec[i])
    end
    for (i, p) in zip(dofnums, vec)
        v[i] += p
    end
    vec, dofnums = rand(7), [29, 15, 1, 7, 3, 6, 2]
    for i in 1:length(vec)
        assemble!(a, dofnums[i], vec[i])
    end
    for (i, p) in zip(dofnums, vec)
        v[i] += p
    end
    w = finish!(a)
    @test norm(v-w) / norm(w) <= 1.0e-9
end
end
using .mvass1
mvass1.test()

module mass1
using Elfel
using Elfel.RefShapes: RefShapeTriangle, manifdim, RefShapeInterval
using Elfel.FElements: FE, refshape
using Elfel.FElements: bfun, bfungradpar
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Test
function test()
    a = SysmatAssemblerSparse(0.0)
    start!(a, 7, 7)
    m = [0.24406   0.599773    0.833404  0.0420141
    0.786024  0.00206713  0.995379  0.780298
    0.845816  0.198459    0.355149  0.224996]
    gi = [1 7 5]
    gj = [5 2 1 4]
    for j in 1:size(m, 2), i in 1:size(m, 1)
          assemble!(a, gi[i], gj[j], m[i, j])
    end

    m = [0.146618  0.53471   0.614342    0.737833
    0.479719  0.41354   0.00760941  0.836455
    0.254868  0.476189  0.460794    0.00919633
    0.159064  0.261821  0.317078    0.77646
    0.643538  0.429817  0.59788     0.958909]
    gi =  [2 3 1 7 5]
    gj = [6 7 3 4]
    for j in 1:size(m, 2), i in 1:size(m, 1)
          assemble!(a, gi[i], gj[j], m[i, j])
    end
    A = finish!(a)
    @test abs.(maximum([0.833404 0.0 0.460794 0.05121043 0.24406 0.254868 0.476189;
        0.0 0.0 0.614342 0.737833 0.0 0.146618 0.53471; 0.0 0.0 0.00760941 0.836455 0.0 0.479719 0.41354; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.355149 0.0 0.59788 1.183905 0.845816 0.643538 0.429817; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.995379 0.0 0.0 0.780298 0.786024 0.0 0.0]  - A)) < 1.0e-5
    #     a = SysmatAssemblerSparse(0.0)
    #     startassembly!(a, 5, 5, 3, 7, 7)
    #     m = [0.24406   0.599773    0.833404  0.0420141
    #     0.786024  0.00206713  0.995379  0.780298
    #     0.845816  0.198459    0.355149  0.224996]
    #     assemble!(a, m, [1 7 5], [5 0 1 4])
    #     m = [0.146618  0.53471   0.614342    0.737833
    #  0.479719  0.41354   0.00760941  0.836455
    #  0.254868  0.476189  0.460794    0.00919633
    #  0.159064  0.261821  0.317078    0.77646
    #  0.643538  0.429817  0.59788     0.958909]
    #     assemble!(a, m, [2 3 1 0 5], [6 7 3 4])
    #     A = makematrix!(a)
    #     # @show Matrix(A)
    #     @test abs.(maximum([0.833404 0.0 0.460794 0.05121043 0.24406 0.254868 0.476189;
    # 0.0 0.0 0.614342 0.737833 0.0 0.146618 0.53471; 0.0 0.0 0.00760941 0.836455 0.0 0.479719 0.41354; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.355149 0.0 0.59788 1.183905 0.845816 0.643538 0.429817; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.995379 0.0 0.0 0.780298 0.786024 0.0 0.0]  - A)) < 1.0e-5

end
end
using .mass1
mass1.test()


text:
module NaiveNASlib

import Statistics
import Logging
using Logging: LogLevel, @logmsg

import JuMP
using JuMP: @variable, @constraint, @objective, @expression, MOI, MOI.INFEASIBLE, MOI.FEASIBLE_POINT
import Cbc
import Functors
using Functors: @functor, functor

# Computation graph
export CompGraph, nvertices, vertices, findvertices, inputs, outputs, name

# Vertex size operations
export nin, nout, Δnin!, Δnout!, Δsize!, relaxed

# Connectivity mutation
export remove!, insert!, create_edge!, remove_edge!

# Create vertices
export inputvertex, immutablevertex, absorbvertex, invariantvertex, conc

include("vertex.jl")
include("compgraph.jl")
include("prettyprint.jl")

include("mutation/vertex.jl")
include("mutation/graph.jl")

include("mutation/jumpnorm.jl")
include("mutation/sizestrategies.jl")

include("mutation/size.jl")
include("mutation/select.jl")
include("mutation/structure.jl")

include("api/vertex.jl")
include("api/size.jl")
include("api/Advanced.jl")
include("api/Extend.jl")

end # module


text:
# find julia files to be converted into README.md

function gernerate_README(examplestoconvert, folder)

    ps = map(examplestoconvert) do str
        fs = filter(s -> split(s, ".")[end] == "jl", readdir(folder*str))
        if length(fs) > 0
            [str, fs[1]]
        else
            String[]
        end
    end
    filter!(p -> p != String[], ps)

    # This code below writes the README.md
    for p in ps
        s1 = string(folder,p[1],"/",p[2])
        str_arr = split(read(open(s1,"r"), String), "\n")

        # convert to md format
        str_arr = map(str_arr) do s
            if length(s) > 2 && s[1] == '#' && s[2] == ' '
              string(s[3:end],"\n")
            else
              string(s,"\n")
            end
        end

        s2 = string(folder,p[1],"/README.md")
        f = open(s2,"w")
        write(f, string(str_arr...))
        close(f)
    end
end


text:
using Flux

struct MLP
    layers
end

@Flux.functor MLP

"""
    MLP(inopts::Int, layeropts::Union{Int,Pair{Int}}...; activation = relu, outactivation = nothing)

Build a MultiLayer Perceptron.

# Examples
```jldoctest
julia> model = MLP(3, 5, 10 => [:act => relu, :dropout => 0.2], 5 => [:act => relu, :batchnorm],  2; activation = σ, outactivation = identity)
MLP(Chain(Dense(3, 5, σ), Dense(5, 10, relu), Dropout{Float64}(0.2, Colon(), true), Dense(10, 5, relu), BatchNorm(5), Dense(5, 2)))
```
"""
function MLP(inopts::Int, layeropts::Union{Int,Pair{Int}}...; activation = relu, outactivation = nothing)
    outactivation = outactivation === nothing ? activation : outactivation
    function extract_layeropts(l, i)
        act = i < length(layeropts) ? activation : outactivation
        drop = 0.0
        bn = false
        if l isa Pair{Int}
            opts = last(l)
            opts = opts isa Union{AbstractArray,Tuple} ? opts : [opts]
            for o in opts
                o = o isa Symbol ? o => true : o
                if first(o) == :act
                    act = last(o)
                elseif first(o) == :dropout
                    drop = last(o)
                elseif first(o) == :batchnorm
                    bn = true
                else
                    error("unknown layer option: $(first(o))")
                end
            end
        end
        (size = first(l), activation = act, batchnorm = bn, dropout = drop)
    end

    layers = []
    nprev = inopts
    for (i, l) in enumerate(layeropts)
        opts = extract_layeropts(l, i)
        push!(layers, Dense(nprev, opts.size, opts.activation))
        opts.batchnorm && push!(layers, BatchNorm(opts.size))
        opts.dropout > 0.0 && push!(layers, Dropout(opts.dropout))
        nprev = first(l)
    end
    MLP(Chain(layers...))
end

"""
    (m::MLP)(x::AbstractArray)

Evaluate MLP.

# Examples
```jldoctest
julia> using Random; Random.seed!(1);

julia> X = randn(Float32, 3, 4)
3×4 Array{Float32,2}:
  0.297288  -0.0104452   2.29509   0.431422
  0.382396  -0.839027   -2.26709   0.583708
 -0.597634   0.311111    0.529966  0.963272

julia> model = MLP(3, 5, 2)
MLP(Chain(Dense(3, 5, relu), Dense(5, 2, relu)))

julia> model(X)
Tracked 2×4 Array{Float32,2}:
 0.0  0.0  0.0  0.377914
 0.0  0.0  0.0  0.687224
```
"""
(m::MLP)(X::AbstractArray) = m.layers(X)

max_stars_repo_name: dhonza/FluxGoodies

max_issues_repo_name: dhonza/FluxGoodies

max_forks_repo_name: dhonza/FluxGoodies


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Represents a Fibre Channel volume. Fibre Channel volumes can only be mounted as read/write once. Fibre Channel volumes support ownership management and SELinux relabeling.

    IoK8sApiCoreV1FCVolumeSource(;
        fsType=nothing,
        lun=nothing,
        readOnly=nothing,
        targetWWNs=nothing,
        wwids=nothing,
    )

    - fsType::String : Filesystem type to mount. Must be a filesystem type supported by the host operating system. Ex. \"ext4\", \"xfs\", \"ntfs\". Implicitly inferred to be \"ext4\" if unspecified.
    - lun::Int32 : Optional: FC target lun number
    - readOnly::Bool : Optional: Defaults to false (read/write). ReadOnly here will force the ReadOnly setting in VolumeMounts.
    - targetWWNs::Vector{String} : Optional: FC target worldwide names (WWNs)
    - wwids::Vector{String} : Optional: FC volume world wide identifiers (wwids) Either wwids or combination of targetWWNs and lun must be set, but not both simultaneously.
"""
mutable struct IoK8sApiCoreV1FCVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    lun::Any # spec type: Union{ Nothing, Int32 } # spec name: lun
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    targetWWNs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: targetWWNs
    wwids::Any # spec type: Union{ Nothing, Vector{String} } # spec name: wwids

    function IoK8sApiCoreV1FCVolumeSource(;fsType=nothing, lun=nothing, readOnly=nothing, targetWWNs=nothing, wwids=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("lun"), lun)
        setfield!(o, Symbol("lun"), lun)
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("targetWWNs"), targetWWNs)
        setfield!(o, Symbol("targetWWNs"), targetWWNs)
        validate_property(IoK8sApiCoreV1FCVolumeSource, Symbol("wwids"), wwids)
        setfield!(o, Symbol("wwids"), wwids)
        o
    end
end # type IoK8sApiCoreV1FCVolumeSource

const _property_map_IoK8sApiCoreV1FCVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("lun")=>Symbol("lun"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("targetWWNs")=>Symbol("targetWWNs"), Symbol("wwids")=>Symbol("wwids"))
const _property_types_IoK8sApiCoreV1FCVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("lun")=>"Int32", Symbol("readOnly")=>"Bool", Symbol("targetWWNs")=>"Vector{String}", Symbol("wwids")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1FCVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1FCVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1FCVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1FCVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1FCVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1FCVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1FCVolumeSource)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1FCVolumeSource }, name::Symbol, val)
end


text:
using JultraDark
using Test, Documenter
using SafeTestsets
using MPI

MPI.Init()

doctest(JultraDark)

@safetestset "grids.jl" begin
    include("grids.jl")
end

@testset "Phase grad Grids" begin
    include("phase_gradient.jl")
end

@safetestset "Evolution" begin
    include("evolution.jl")
end

# Run examples
@safetestset "full sim" begin
    include("sims/full.jl")
end

# Put notebook in module to emulate SafeTestsets
# https://github.com/YingboMa/SafeTestsets.jl/issues/3
module GrowthNotebook
    using NBInclude
    @nbinclude("../examples/growth.ipynb")
end


text:
@testset "conv_bias_act" begin
    x = rand(4,4,3,3)
    w = rand(2,2,3,3)
    b = rand(1,1,1,3)
    cdims = DenseConvDims(x, w; stride=2)
    @test NNlib.conv_bias_act(x, w, cdims, b, relu) ≈ relu.(conv(x, w, cdims) .+ b) atol=1e-5
end