text:
# 1-D Calculations
"""
    mass_matrix(b::Basis)

    Calculate the mass-matrix for the basis `b` with Gauss quadrature.
"""
function mass_matrix(basis::BSplineBasis{dirichlet})
    x, w = gausslegendre(basis.k + 10)
    a = basis.interval.left
    b = basis.interval.right
    N = length(basis.B)
    Δx = (b-a)/basis.N

    M = zeros(N, N)
    # Do Quadrature over every element -> exact
    for el in 1:basis.N
        b = a + Δx
        xx = a .+ (b - a) .*(x .+ 1) ./ 2

        # Otherwise sometimes the last element is wrong
        if xx[end] > basis.interval.right
            xx[end] = basis.interval.right
        end

        for i in 1:N
            for j in 1:i
                M[i, j] = M[j, i] += 2/Δx * dot(basis[xx, i] .* basis[xx, j], w)
            end
        end
        a = b
    end

    return  M
end

# Todo more precisely, defined on the same domain
function mass_matrix(basis1::BSplineBasis{T}, basis2::BSplineBasis{T}) where {T}
    x, w = gausslegendre(basis1.k + 2)
    a1 = basis1.interval.left
    b1 = basis1.interval.right
    N1 = length(basis1.B)
    Δx1 = (b1-a1)/basis1.N

    a2 = basis2.interval.left
    b2 = basis2.interval.right
    N2 = length(basis2.B)
    Δx2 = (b2-a2)/basis2.N

    M = zeros(N1, N2)
    # Do Quadrature over every element -> exact
    for el1 in 1:basis1.N
        b1 = a1 + Δx1
        xx = a1 .+ (b1 - a1) .*(x .+ 1) ./ 2

        if xx[end] > basis1.interval.right
            xx[end] = basis1.interval.right
        end

        for i in 1:N1
            for j in 1:N2
                M[i, j] += 2/(b1 - a1) .* dot(basis1[xx, i] .* basis2[xx, j], w)
            end
        end
        a1 = b1
    end

    return  M
end

# For the stiffness_matrix of basis b, call stiffness_matrix(b')
# Or should we take b as an input?
# Care for the factor in the integral, the derived basis functions, i.e. b', have already 2/Δx as an factor!
"""
    stiffness_matrix(b'::Basis)

    Calculate the stiffnes-matrix for the basis `b` with Gauss quadrature. Caution: We directly input b' here as a basis.
"""
function stiffness_matrix(basis::BSplineBasis{dirichlet})
    x, w = gausslegendre(basis.k + 1)
    a = basis.interval.left
    b = basis.interval.right
    N = length(basis.B)
    Δx = (b-a)/basis.N

    M = zeros(N+1, N+1)
    # The derivative of the basis functions is the known formula for BSplines
    # B'_i(x) = D_i(x) - D_{i+1}(x), where D_i are the basisfunctions of b'
    for el in 1:basis.N
        b = a + Δx
        xx = a .+ (b - a) .*(x .+ 1) ./ 2

        if xx[end] > basis.interval.right
            xx[end] = basis.interval.right
        end

        M[1, 1] += 2/(b - a) .* dot((-basis[xx, 1]) .* (-basis[xx, 1] ), w)

        for j in 1:N-1
            M[1, j+1] = M[j+1, 1] += 2/(b - a) .* dot((-basis[xx, 1]) .* (basis[xx, j] - basis[xx, j+1] ), w)
        end

        for i in 1:N-1
            for j in 1:i

                M[i+1, j+1] = M[j+1, i+1] += 2/(b - a) .* dot((basis[xx, i] .- basis[xx, i+1]) .* (basis[xx, j] .- basis[xx, j+1]), w)
            end
        end

        for j in 1:N-1
            M[N+1, j+1] = M[j+1, N+1] += 2/(b - a) .* dot((basis[xx, N]) .* (basis[xx, j] - basis[xx, j+1] ), w)
        end

        M[N+1, N+1] += 2/(b - a) .* dot((basis[xx, N]) .* (basis[xx, N] ), w)

        a = b
    end

    return  M
end

# For periodic boundary conditions on the basis use the circulant methods
function mass_matrix(basis::BSplineBasis{periodic})
    x, w = gausslegendre(basis.k+1)
    a = basis.interval.left
    b = basis.interval.right
    N = length(basis.B)
    Δx = (b-a)/basis.N
    c = zeros(N)

    for el in 1:basis.N
        b = a + Δx
        xx = a .+ (b - a) .*(x .+ 1) ./ 2

        for i in 1:N
            c[i] += 2/(b - a) .* dot(basis[xx, i] .* basis[xx, 1], w)
        end
        a = b
    end

    return Circulant(c)
end

# Input derivative basis
function stiffness_matrix(basis::BSplineBasis{periodic})
    x, w = gausslegendre(basis.k+1)
    a = basis.interval.left
    b = basis.interval.right
    N = length(basis.B)
    Δx = (b-a)/basis.N
    c = zeros(N)

    for el in 1:basis.N
        b = a + Δx
        xx = a .+ (b - a) .*(x .+ 1) ./ 2

        for i in 1:N-1
            c[i] += 2/(b - a) .* dot((basis[xx, i] .- basis[xx, i+1]) .* (basis[xx, 1] .- basis[xx, 2]), w)
        end
        c[N] += 2/(b - a) .* dot((basis[xx, N] .- basis[xx, 1]) .* (basis[xx, 1] .- basis[xx, 2]), w)
        a = b
    end

    return Circulant(c)
end

# Since for now we only have a periodic basis, f has to be as well
function L2_prod(f::Function, basis::Basis, order::Int64=basis.k+1)
    # what degree do we want here?
    x, w = gausslegendre(order+1)
    a = basis.interval.left
    b = basis.interval.right
    N = length(basis.B)
    Δx = (b-a)/basis.N
    fh = zeros(N)

    for el in 1:basis.N
        b = a + Δx
        xx = a .+ (b - a) .*(x .+ 1) ./ 2

        if xx[end] > basis.interval.right
            xx[end] = basis.interval.right
        end

        for i in 1:N
            fh[i] += 2/(b - a) .* dot(basis[xx, i] .* f.(xx), w)
        end
        a = b
    end

    return fh
end

#maybe find a better name
function col_mat(basis::Basis, m::Int)
    n = length(basis.B)
    x = range(basis.interval.left+1e-14, basis.interval.right-1e-14, length = m)
    M = zeros(m, n)

    for i in 1:n
        M[:,i] = basis[x, i]
    end

    return x, M
end

#n-D Calculations
function mass_matrix(basis::TensorProductBasis{N}) where N
    matrices = []
    for i in 1:N
        push!(matrices, mass_matrix(basis[i]))
    end

    return ⊗(matrices...)
end

function stiffness_matrix(basis::TensorProductBasis{N}) where N
    matrices = []
    for i in 1:N
        push!(matrices,  stiffness_matrix(basis[i]))
    end

    return ⊗(matrices...)
end

#Todo: Generalize the following functions and make them more efficient
function L2_prod(f::Function, basis::BSplineTensorProductBasis{2, T}) where T
    # what degree do we want here?
    X,Y = basis.B

    x1, w1 = gausslegendre(X.k + 1)
    x2, w2 = gausslegendre(Y.k + 1)

    a1 = X.interval.left
    b1 = X.interval.right
    N1 = length(X.B)
    Δx1 = (b1-a1)/X.N

    a2 = Y.interval.left
    b2 = Y.interval.right
    N2 = length(Y.B)
    Δx2 = (b2-a2)/Y.N

    vectors = zeros(N1, N2)

    #loop over basis functions
    for i in 1:N1
        for j in 1:N2
            #loop over elements in x-direction
            for el_x in 1:X.N
                b1 = a1 + Δx1
                xx1 = a1 .+ (b1 - a1) .*(x1 .+ 1) ./ 2

                if xx1[end] > X.interval.right
                    xx1[end] = X.interval.right
                end

                val_y = zeros(length(xx1))
                # at every quadrature point in x-direction, do a full quadrature in y-direction
                for xn in 1:length(xx1)
                    val = 0

                    for el_y in 1:Y.N
                        b2 = a2 + Δx2
                        xx2 = a2 .+ (b2 - a2) .*(x2 .+ 1) ./ 2
                        if xx2[end] > Y.interval.right
                            xx2[end] = Y.interval.right
                        end

                        val += 2/(b2 - a2) .* dot(Y[xx2, j] .* f.(xx1[xn], xx2), w2)

                        a2 = b2
                    end

                    val_y[xn] = val
                    a2 = Y.interval.left
                end

                vectors[i, j] += 2/(b1 - a1) .* dot(X[xx1, i] .* val_y, w1)

                a1 = b1
            end

            a1 = X.interval.left
        end

    end

    return reshape(permutedims(vectors, [2 ,1]), N1*N2)
end

function L2_prod(f::Function, basis::BSplineTensorProductBasis{3, T}) where T
    X, Y, Z = basis.B
    N = (length(X.B), length(Y.B), length(Z.B))

    x1, w1 = gausslegendre(X.k + 1)
    x2, w2 = gausslegendre(Y.k + 1)
    x3, w3 = gausslegendre(Z.k + 1)

    a1 = X.interval.left
    b1 = X.interval.right
    N1 = length(X.B)
    Δx1 = (b1-a1)/X.N

    a2 = Y.interval.left
    b2 = Y.interval.right
    N2 = length(Y.B)
    Δx2 = (b2-a2)/Y.N

    a3 = Z.interval.left
    b3 = Z.interval.right
    N3 = length(Z.B)
    Δx3 = (b3-a3)/Z.N

    vectors = zeros(Float64, N1, N2, N3)

    xx1 = zeros(Float64, length(x1))
    xx2 = zeros(Float64, length(x2))
    xx3 = zeros(Float64, length(x3))
    #loop over basis functions
    for i in 1:N1
        for j in 1:N2
            for l in 1:N3
                #loop over elements in x-direction
                for el_x in 1:X.N
                    b1 = a1 + Δx1
                    xx1 = a1 .+ (b1 - a1)/2 * (x1 .+ 1)

                    if xx1[end] > X.interval.right
                        xx1[end] = X.interval.right
                    end

                    val_y = zeros(Float64, length(xx1))
                    # at every quadrature point in x-direction, do a full quadrature in y-direction
                    for xn in 1:length(xx1)
                        val2 = 0.0

                        for el_y in 1:Y.N
                            b2 = a2 + Δx2
                            xx2 = a2 .+ (b2 - a2)/2 * (x2 .+ 1)
                            if xx2[end] > Y.interval.right
                                xx2[end] = Y.interval.right
                            end
                                val_z = zeros(Float64, length(xx2))
                                # at every quadrature point in x-direction, do a full quadrature in y-direction
                                for yn in 1:length(xx2)
                                    val3 = 0.0

                                    for el_z in 1:Z.N
                                        b3 = a3 + Δx3
                                        xx3 = a3 .+ (b3 - a3)/2 * (x3 .+ 1)
                                        if xx3[end] > Z.interval.right
                                            xx3[end] = Z.interval.right
                                        end

                                        val3 += 2/(b3 - a3) * dot(Z[xx3, l] .* f.(xx1[xn], xx2[yn], xx3), w3)

                                        a3 = b3
                                    end

                                    val_z[yn] = val3
                                    a3 = Z.interval.left
                                end

                            val2 += 2/(b2 - a2) * dot(Y[xx2, j] .* val_z, w2)

                            a2 = b2
                        end

                        val_y[xn] = val2
                        a2 = Y.interval.left
                    end
                    vectors[i, j, l] += 2/(b1 - a1) * dot(X[xx1, i] .* val_y, w1)

                    a1 = b1
                end

                a1 = X.interval.left
            end
        end

    end

    return reshape(permutedims(vectors, [3, 2 ,1]), N1*N2*N3)
end

#x = (x, y)  1-d interval discretizations
function col_mat(basis::TensorProductBasis{2}, mm::Int)
    X,Y = basis.B
    N = (length(X.B), length(Y.B))

    x = range(X.interval.left+1e-14, X.interval.right-1e-14, length = mm)
    y = range(Y.interval.left+1e-14, Y.interval.right-1e-14, length = mm)
    xy = [(a,b) for a in x for b in y]

    M1 = zeros(mm, N[1])
    M2 = zeros(mm, N[2])

    for i in 1:N[1]
        M1[:, i] = X[x, i]
    end

    for i in 1:N[2]
        M2[:, i] = Y[y, i]
    end

    return xy, M1 ⊗ M2
end

function col_mat(basis::BSplineTensorProductBasis{3, T}, mm::Int) where T
    X, Y, Z = basis.B
    N = (length(X.B), length(Y.B), length(Z.B))

    x = range(X.interval.left+1e-14, X.interval.right-1e-14, length = mm)
    y = range(Y.interval.left+1e-14, Y.interval.right-1e-14, length = mm)
    z = range(Z.interval.left+1e-14, Z.interval.right-1e-14, length = mm)
    xyz = [(a, b, c) for a in x for b in y for c in z]

    M1 = zeros(mm, N[1])
    M2 = zeros(mm, N[2])
    M3 = zeros(mm, N[3])

    for i in 1:N[1]
        M1[:, i] = X[x, i]
    end

    for i in 1:N[2]
        M2[:, i] = Y[y, i]
    end

    for i in 1:N[3]
        M3[:, i] = Z[z, i]
    end

    return xyz, M1 ⊗ M2 ⊗ M3
end

function L2_prod(f::Function, basis::SArray{Tuple{3}, BSplineTensorProductBasis{3, T}}) where {T}
    fh = []

    for i in 1:3
        push!(fh, L2_prod((x,y,z) -> f(x,y,z)[i], basis[i]))
    end

    return SVector(fh...)
end

function L2_prod(f::Function, basis::SArray{Tuple{2}, BSplineTensorProductBasis{2, T}}) where {T}
    fh = []

    for i in 1:2
        push!(fh, L2_prod((x,y) -> f(x,y)[i], basis[i]))
    end

    return SVector(fh...)
end

function L2_proj(f::Function, basis::SArray{Tuple{3}, BSplineTensorProductBasis{3, T}}) where {T}
    fh = []
    M = mass_matrix.(basis)

    for i in 1:3
        push!(fh, M[i]\L2_prod((x,y,z) -> f(x,y,z)[i], basis[i]))
    end

    return SVector(fh...)
end

function L2_proj(f::Function, basis::BSplineTensorProductBasis{3, T}) where {T}
    M = mass_matrix(basis)

    fh = M\L2_prod(f, basis)
    return SVector{length(fh)}(fh)
end

function L2_proj(f::Function, basis::SArray{Tuple{2}, BSplineTensorProductBasis{2, T}}) where {T}
    fh = []
    M = mass_matrix.(basis)

    for i in 1:2
        push!(fh, M[i]\L2_prod((x,y) -> f(x,y)[i], basis[i]))
    end

    return SVector(fh...)
end

function L2_proj(f::Function, basis::BSplineTensorProductBasis{2, T}) where {T}
    M = mass_matrix(basis)
    fh = M\L2_prod(f, basis)
    return SVector{length(fh)}(fh)
end

function L2_proj(f::Function, basis::BSplineBasis{T}) where {T}
    M = mass_matrix(basis)
    fh = M\L2_prod(f, basis)
    return SVector{length(fh)}(fh)
end

##weighted mass matrices
function weighted_mass_matrix(basis::Basis, ω::Function)
    x, w = gausslegendre(basis.k + 1)
    a = basis.interval.left
    b = basis.interval.right
    N = length(basis.B)
    Δx = (b-a)/basis.N

    M = zeros(N, N)

    for el in 1:basis.N
        b = a + Δx
        xx = a .+ (b - a) .*(x .+ 1) ./ 2
        if xx[end] > basis.interval.right
            xx[end] = basis.interval.right
        end

        for i in 1:N
            for j in 1:N
                M[i, j] += 2/(b - a) .* dot(basis[xx, i] .* basis[xx, j] .* ω.(xx), w)
            end
        end
        a = b
    end

    return  M
end

#ω(x,y)
function weighted_mass_matrix(basis::BSplineTensorProductBasis{2, T}, ω::Function) where {T}
    X,Y = basis.B

    #maybe we need higher order
    x1, w1 = gausslegendre(X.k + 1)
    x2, w2 = gausslegendre(Y.k + 1)

    a1 = X.interval.left
    b1 = X.interval.right
    N1 = length(X.B)
    Δx1 = (b1-a1)/X.N

    a2 = Y.interval.left
    b2 = Y.interval.right
    N2 = length(Y.B)
    Δx2 = (b2-a2)/Y.N

    M = zeros(N1*N2, N1*N2)

    vectors = zeros(N1, N2)
    for k in 1:N1
        for l in 1:N2

            vectors = zeros(N1, N2)
            for i in 1:N1
                for j in 1:N2
                    #loop over elements in x-direction
                    for el_x in 1:X.N
                        b1 = a1 + Δx1
                        xx1 = a1 .+ (b1 - a1) .*(x1 .+ 1) ./ 2
                        if xx1[end] > X.interval.right
                            xx1[end] = X.interval.right
                        end

                        val_y = zeros(length(xx1))
                        # at every quadrature point in x-direction, do a full quadrature in y-direction
                        for xn in 1:length(xx1)
                            val = 0

                            for el_y in 1:Y.N
                                b2 = a2 + Δx2
                                xx2 = a2 .+ (b2 - a2) .*(x2 .+ 1) ./ 2
                                if xx2[end] > Y.interval.right
                                    xx2[end] = Y.interval.right
                                end
                                val += 2/(b2 - a2) .* dot(Y[xx2, j] .* Y[xx2, l] .* ω.(xx1[xn], xx2), w2)

                                a2 = b2
                            end

                            val_y[xn] = val
                            a2 = Y.interval.left
                        end
                        vectors[i, j] += 2/(b1 - a1) .* dot(X[xx1, i] .* X[xx1, k] .* val_y, w1)

                        a1 = b1
                    end
                    a1 = X.interval.left
                end

            end
            M[:, N2*(k-1) + l ] = reshape(permutedims(vectors, [2 ,1]), N1*N2)

        end
    end

    return  M
end

#special structure because of V1
#w(x, y)
function interior_product_matrix(basis::SArray{Tuple{2}, BSplineTensorProductBasis{2, T}}, ω::Function) where T
    M = weighted_mass_matrix(basis[1], basis[2], ω)

    return SArray(M, M')
end

#ω(x,y)
#V1 = SVector(X' ⊗ Y, X ⊗ Y')
function weighted_mass_matrix(basis1::BSplineTensorProductBasis{2, periodic}, basis2::BSplineTensorProductBasis{2, periodic}, ω::Function)
    X, Y = basis1.B
    X2, Y2 = basis2.B

    x1, w1 = gausslegendre(X.k + 1)
    x2, w2 = gausslegendre(Y.k + 1)

    # Assume both component bases are defined on the same domain
    a1 = X.interval.left
    b1 = X.interval.right
    N1 = length(X.B)
    Δx1 = (b1-a1)/X.N

    a2 = Y.interval.left
    b2 = Y.interval.right
    N2 = length(Y.B)
    Δx2 = (b2-a2)/Y.N

    M = zeros(N1*N2, N1*N2)

    vectors = zeros(N1, N2)
    for k in 1:length(X2.B) # build the full matrix
        for l in 1:length(Y2.B)

            vectors = zeros(N1, N2)
            for i in 1:N1
                for j in 1:N2
                    #loop over elements in x-direction
                    for el_x in 1:X.N
                        b1 = a1 + Δx1
                        xx1 = a1 .+ (b1 - a1) .*(x1 .+ 1) ./ 2
                        if xx1[end] > X.interval.right
                            xx1[end] = X.interval.right
                        end
                        val_y = zeros(length(xx1))
                        # at every quadrature point in x-direction, do a full quadrature in y-direction
                        for xn in 1:length(xx1)
                            val = 0

                            for el_y in 1:Y.N
                                b2 = a2 + Δx2
                                xx2 = a2 .+ (b2 - a2) .*(x2 .+ 1) ./ 2
                                if xx2[end] > Y.interval.right
                                    xx2[end] = Y.interval.right
                                end
                                val += 2/(b2 - a2) .* dot(Y[xx2, j] .* Y2[xx2, l] .* ω.(xx1[xn], xx2), w2)

                                a2 = b2
                            end

                            val_y[xn] = val
                            a2 = Y.interval.left
                        end
                        vectors[i, j] += 2/(b1 - a1) .* dot(X[xx1, i] .* X2[xx1, k] .* val_y, w1)

                        a1 = b1
                    end
                    a1 = X.interval.left
                end

            end
            M[:, N2*(k-1) + l ] = reshape(permutedims(vectors, [2 ,1]), N1*N2)

        end
    end

    return  M
end

function interior_product_matrix(basis::SArray{Tuple{3}, BSplineTensorProductBasis{3, T}}, ω::SArray{Tuple{3}, Function}) where T

    M1 = weighted_mass_matrix(V1[1], V1[2], ω[3])
    M2 = weighted_mass_matrix(V1[1], V1[3], ω[2])
    M3 = weighted_mass_matrix(V1[2], V1[3], ω[1])

    #M = [0 -M1 M2;
    #     M1 0 -M3;
    #     -M2 M3 0]

    #we could take it as a curl matrix
    return SVector(M1, M2, M3)
end

function L2_penalty(uh::T, basis::BSplineBasis, ϵ) where T

    M = mass_matrix(basis)
    M0 = copy(M)

    M0[1, 1] *= (1 + 1/ϵ)
    M0[end, end] *= (1 + 1/ϵ)

    return M0\(M*uh)

end

#returns invertible matrix, to get coefficients from point values
function col_mat_proj(basis::BSplineTensorProductBasis{2, periodic})
    X,Y = basis.B
    N = (length(X.B), length(Y.B))
    k = (X.k, Y.k)

    #periodic
    x = X.t[k[1]+2:N[1]+k[1]+1]
    y = Y.t[k[2]+2:N[2]+k[2]+1]

    xy = [(a,b) for a in x for b in y]

    M1 = zeros(N[1], N[1])
    M2 = zeros(N[2], N[2])

    for i in 1:N[1]
        M1[:, i] = X[x, i]
    end

    for i in 1:N[2]
        M2[:, i] = Y[y, i]
    end

    return xy, M1 ⊗ M2
end

#returns invertible matrix, to get coefficients from point values
function col_mat_proj(basis::BSplineTensorProductBasis{2, dirichlet})
    X,Y = basis.B
    N = (length(X.B), length(Y.B))
    k = (X.k, Y.k)

    x = range(X.interval.left+1e-14, X.interval.right-1e-14, length= N[1])
    y = range(Y.interval.left+1e-14, Y.interval.right-1e-14, length= N[2])
    xy = [(a,b) for a in x for b in y]

    M1 = zeros(N[1], N[1])
    M2 = zeros(N[2], N[2])

    for i in 1:N[1]
        M1[:, i] = X[x, i]
    end

    for i in 1:N[2]
        M2[:, i] = Y[y, i]
    end

    return xy, M1 ⊗ M2
end

#returns invertible matrix, to get coefficients from point values
function col_mat_proj(basis::BSplineBasis{dirichlet})
    N = length(basis.B)
    #k = basis.k

    x = range(basis.interval.left+1e-14, basis.interval.right-1e-14, length= N)

    M1 = zeros(N, N)

    for i in 1:N
        M1[:, i] = basis[x, i]
    end

    return x, M1
end


text:
module F16Model

using ForwardDiff, HDF5, Interpolations
include("NonlinearF16Model.jl");
end # module

ext: jl


text:
using SpatialFields
using Base.Test
# using PyPlot

function hrbf_2d_5th_power()
	points = SVector{2, Float64}[[1; 0], [0; 1], [-1; 0], [0;-1]]
	normals = SVector{2, Float64}[[1; 1], [0; 1], [-1; 1], [0; -1]]
	num_points = size(points, 2)

	field = HermiteRadialField(points, normals, x -> x^5)

	X = linspace(-2, 2)
	Y = linspace(-2, 2)
	Z = zeros(length(X), length(Y))

    for i = 1:length(X)
	    for j = 1:length(Y)
	        Z[j,i] = field([X[i], Y[j]])
	    end
	end

	for i in 1:size(points, 2)
		@test isapprox(field(points[i]), 0, atol=1e-6)
		eps = 1e-4
		nudged_point = points[i] + eps * SVector(normals[i])
		@test isapprox(dot(SpatialFields.gradient(field)(points[i]), SVector(normals[i])) * eps, field(nudged_point), atol=1e-6)
	end

	# clf()
	# hold(true)
	# PyPlot.contour(X, Y, Z, [-0.1, 0.0, 0.1])
	# for i = 1:num_points
	#     PyPlot.quiver(points[:,i]..., normals[:,i]...)
	# end
	# axis("equal")
	# show()
end

hrbf_2d_5th_power()


text:
function cd3_dissociation!(du,u,p,t,dens,cons)
    # unpack rates and constants
    nᵣ,nₓ,n₃ = u
    k₋₃ = p[1]
    k₋₁,k₋₂ = cons
    # mᵣ,mₗ,mₓ,A = dens

    # model
    du[1] = dnᵣ = - k₋₁*nᵣ
    du[2] = dnₓ = - k₋₂*nₓ
    du[3] = dn₃ = - k₋₃*n₃

end

function bimolecular_diss_bell_activated_state_w_cons!(du,u,p,t,dens,cons;f=0.0)
    # unpack rates and constants
    Pᵣ,Pₐ = u
    x₋ᵒf,
        k₋ᵒs,x₋ᵒs,
        kᵒa,xᵒa,
        k₋ᵒa,x₋ᵒa = p
    k₋ᵒf = cons[1] # zero force diss, generally obtained from thermal fluct.

    k₋f = bell_diss_kT(k₋ᵒf, x₋ᵒf, f) # fast diss form L-R
    k₋s = bell_diss_kT(k₋ᵒs, x₋ᵒs, f) # slow diss from L-R*
    ka  = bell_diss_kT(kᵒa, xᵒa, f)   # activation from L-R to L-R*
    k₋a = bell_ass_kT(k₋ᵒa, x₋ᵒa, f)  # deactivation from L-R* to L-R

    # model
    du[1] = dPᵣ = - k₋f.*Pᵣ - ka.*Pᵣ + k₋a.*Pₐ
    du[2] = dPₐ = + ka.*Pᵣ - k₋a.*Pₐ - k₋s.*Pₐ
end

function cd3_tri_dissociation_bell!(du,u,p,t,dens,cons;f=0.0)
    # unpack rates and constants
    Pᵣ₁,Pₐ₁,
        Pᵣ₂,Pₐ₂,
        Pᵣ₃,Pₐ₃ = u

    x₋₃ᵒf,
        k₋₃ᵒs, x₋₃ᵒs,
        k₃ᵒa, x₃ᵒa,
        k₋₃ᵒa, x₋₃ᵒa = p

    #  for now declare the kᵢ from bimolecular as global variables
    k₋₁ᵒf, x₋₁ᵒf,
        k₋₁ᵒs, x₋₁ᵒs,
        k₁ᵒa, x₁ᵒa,
        k₋₁ᵒa, x₋₁ᵒa,
        k₋₂ᵒf, x₋₂ᵒf,
        k₋₂ᵒs, x₋₂ᵒs,
        k₂ᵒa, x₂ᵒa,
        k₋₂ᵒa, x₋₂ᵒa,
        k₋₃ᵒf         = cons

    k₋₁f = bell_diss_kT(k₋₁ᵒf, x₋₁ᵒf, f) # fast diss form L-R
    k₋₁s = bell_diss_kT(k₋₁ᵒs, x₋₁ᵒs, f) # slow diss from L-R*
    k₁a  = bell_diss_kT(k₁ᵒa, x₁ᵒa, f)   # activation from L-R to L-R*
    k₋₁a = bell_ass_kT(k₋₁ᵒa, x₋₁ᵒa, f)  # deactivation from L-R* to L-R

    k₋₂f = bell_diss_kT(k₋₂ᵒf, x₋₂ᵒf, f) # fast diss form L-R
    k₋₂s = bell_diss_kT(k₋₂ᵒs, x₋₂ᵒs, f) # slow diss from L-R*
    k₂a  = bell_diss_kT(k₂ᵒa, x₂ᵒa, f)   # activation from L-R to L-R*
    k₋₂a = bell_ass_kT(k₋₂ᵒa, x₋₂ᵒa, f)  # deactivation from L-R* to L-R

    k₋₃f = bell_diss_kT(k₋₃ᵒf, x₋₃ᵒf, f) # fast diss form L-R
    k₋₃s = bell_diss_kT(k₋₃ᵒs, x₋₃ᵒs, f) # slow diss from L-R*
    k₃a  = bell_diss_kT(k₃ᵒa, x₃ᵒa, f)   # activation from L-R to L-R*
    k₋₃a = bell_ass_kT(k₋₃ᵒa, x₋₃ᵒa, f)  # deactivation from L-R* to L-R

    # model
    # γϵ
    du[1] = dPᵣ₁ = - k₋₁f*Pᵣ₁ - k₁a*Pᵣ₁ + k₋₁a*Pₐ₁
    du[2] = dPₐ₁ = + k₁a*Pᵣ₁ - k₋₁a*Pₐ₁ - k₋₁s*Pₐ₁
    # δϵ
    du[3] = dPᵣ₂ = - k₋₂f*Pᵣ₂ - k₂a*Pᵣ₂ + k₋₂a*Pₐ₂
    du[4] = dPₐ₂ = + k₂a*Pᵣ₂ - k₋₂a*Pₐ₂ - k₋₂s*Pₐ₂
    # tri
    du[5] = dPᵣ₃ = - k₋₃f*Pᵣ₃ - k₃a*Pᵣ₃ + k₋₃a*Pₐ₃
    du[6] = dPₐ₃ = + k₃a*Pᵣ₃ - k₋₃a*Pₐ₃ - k₋₃s*Pₐ₃

end

## functions to calculate initial bond distribution vector (at t=0)

#=
if the delay between binding and beginning of the measurement is sufficiently
long, the probability of occupying either state will reach an equilibrium. This
equilibrium is independent of force
=#
function initial_state(kᵒa,k₋ᵒa)
        u₀1 = k₋ᵒa / (kᵒa + k₋ᵒa)
        u₀ = [u₀1, (1.0 - u₀1)]
end

#=
assumes that the time between binding and measurement is short relative to the
transition rates between the two bound states, then the bond will stay in the
initial state it occupied upon binding. The probabilities of occupying each of
these states can be derived from the principle of detailed balance.
=#
function initial_state(k₋ᵒf,k₋ᵒs,kᵒa,k₋ᵒa)
        u₀1 = (k₋ᵒa * k₋ᵒf) / (k₋ᵒa * k₋ᵒf + kᵒa * k₋ᵒs)
        u₀ = [u₀1, (1.0 - u₀1)]
end

#=
assumes enough time to reach equilibrium, and that the bond experiences some
fraction a of the measured force.
IMPORTANT: must add a*f to function call: initial_state(p, a * f)
=#

function initial_state(p::Vector,f)
        k₋ᵒf,x₋ᵒf,
        k₋ᵒs,x₋ᵒs,
        kᵒa,xᵒa,
        k₋ᵒa,x₋ᵒa = p

        k₋f = bell_diss_kT(k₋ᵒf, x₋ᵒf, f) # fast diss form L-R
        k₋s = bell_diss_kT(k₋ᵒs, x₋ᵒs, f) # slow diss from L-R*
        ka  = bell_diss_kT(kᵒa, xᵒa, f)   # activation from L-R to L-R*
        k₋a = bell_ass_kT(k₋ᵒa, x₋ᵒa, f)

        u₀1 = k₋a / (ka + k₋a)
        u₀ = [u₀1, (1.0 - u₀1)]
end

# function
# u = []
# for f in Fs
#         u₀ = initial_state(rates,f)
#         push!(u,u₀)
# end
#
# plot(reduce(hcat,u))

function initial_u(k,M)
    if ~isempty(M.cons)
        k = vcat(M.cons,k)
    end

    if M.u₀_type == "equilibrium"
        k₋ᵒf,x₋ᵒf,k₋ᵒs,x₋ᵒs,kᵒa,xᵒa,k₋ᵒa,x₋ᵒa = k
        new_u₀ = initial_state(kᵒa,k₋ᵒa)

    elseif M.u₀_type == "detailed_balance"
        k₋ᵒf,x₋ᵒf,k₋ᵒs,x₋ᵒs,kᵒa,xᵒa,k₋ᵒa,x₋ᵒa = k
        new_u₀ = initial_state(k₋ᵒf,k₋ᵒs,kᵒa,k₋ᵒa)

    elseif M.u₀_type == "force"
        f = M.cons
        a = k[end]
        new_u₀ = initial_state(k,a*f)

    else
        new_u₀ = []
    end

    if isempty(M.u₀)
        u₀ = new_u₀
    else
        u₀ = M.u₀
        u₀ = convert.(eltype(new_u₀),u₀)
        u₀ = append!(u₀, new_u₀) # append new initial states to provided initial states
    end
    return u₀
end


text:
"""This is a function to calcuate (x+y)*(x+y)"""
f(x,y)=x^2+2*x*y+y^2
"""This calculates the derivative with respect to x of the
function f(x,y)"""
derivative_of_f1(x,y)=ForwardDiff.derivative(x->f(x,y),x)


text:
include("triangulate_cell.jl")
include("triangulate_unstruct.jl")
include("triangulate_struct.jl")

triangulate(dataset::AbstractVTKMultiblockData) = triangulate(VTKUnstructuredData(dataset))

function triangulate(dataset::AbstractTimeSeriesVTKData)
    return VTKTimeSeriesData(dataset.timemarkers, [triangulate(dataset[i]) for
        i in 1:length(dataset)])
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule PMIx_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("PMIx")
JLLWrappers.@generate_main_file("PMIx", UUID("32165bc3-0280-59bc-8c0b-c33b6203efab"))
end  # module PMIx_jll


text:
using Graphs

function wordindices{T <: String}(word::T)
  wordInd = map((x) -> int(x), collect(word))
  wordLen = length(wordInd)
  return wordInd, wordLen
end

function flooredmod(a::Int, n::Int)
  return int(a - n * floor(a/n))
end

type MinimalPerfectHash
  t1::Array{Int64,1}
  t2::Array{Int64,1}
  G::GenericGraph{
    ExVertex,
    ExEdge{ExVertex},
    Array{ExVertex,1},
    Array{ExEdge{ExVertex},1},
    Array{Array{ExEdge{ExVertex},1},1}}
  #Array{Array{ExVertex,1},1},
  g::Array{Int64,1}
  m::Int64
  n::Int64

  function MinimalPerfectHash{T <: String}(rawwords::Vector{T})
    words = unique(rawwords)
    numwords = int64(length(words))

    G, t1, t2, n = mapping(words)

    # Set all edges to unvisited state
    for v in vertices(G)
      v.attributes["visited"] = false
    end

    # Initialize g functions
    g = Array(Int64, length(vertices(G)))

    mphbase = new(t1, t2, G, g, numwords, n)

    for v in vertices(G)
      if v.attributes["visited"] == false
        g[v.index] = 0
        traverse(v, mphbase)
      end
    end

    return mphbase
  end
end

function hash{T <: String}(word::T, mph::MinimalPerfectHash)
  wordInd, wordLen = wordindices(word)
  u = flooredmod((mph.t1[1:wordLen]' * wordInd)[1], mph.n)
  v = flooredmod((mph.t2[1:wordLen]' * wordInd)[1], mph.n)
  hashval = flooredmod((mph.g[u+1] + mph.g[v+1]), mph.m) + 1
  return hashval
end

function buildgraph(n::Int)
  G = graph(Array(ExVertex,0), Array(ExEdge{ExVertex},0), is_directed = false)
  for i = 1:n
    add_vertex!(G, "edge:" * string(i))
  end
  return G
end

function mapping{T <: String}(words::Array{T})
  numwords = length(words)
  uniqChars = unique(collect(join(words)))
  maxwl = maximum(collect(map((x) -> length(x), words)))

  n = int(25/12 * numwords)
  t1 = Array(Int64, maxwl)
  t2 = Array(Int64, maxwl)

  G = buildgraph(n)

  isCyclic = true
  loopCounter = 0

  while isCyclic
    G = buildgraph(n)
    V = vertices(G)

    rand!(0:(n-1), t1)
    rand!(0:(n-1), t2)

    for i = 1:numwords
      wordInd, wordLen = wordindices(words[i])
      u = flooredmod((t1[1:wordLen]' * wordInd)[1], n)
      v = flooredmod((t2[1:wordLen]' * wordInd)[1], n)
      edg = add_edge!(G, V[u+1], V[v+1])
      edg.attributes["h"] = i - 1
    end

    isCyclic = test_cyclic_by_dfs(G)

    if loopCounter > 20
      throw(ArgumentError("[MPHF:mapping] Could not find acyclic graph for input.\n"))
    end

    loopCounter += 1
  end

  return G, t1, t2, n
end

function traverse(u::ExVertex, mph::MinimalPerfectHash)
  u.attributes["visited"] = true
  for edg in out_edges(u, mph.G)
    w = edg.target
    if !w.attributes["visited"]

      mph.g[w.index] = flooredmod(edg.attributes["h"] - mph.g[u.index], mph.m)
      traverse(w, mph)
    end
  end
end


text:
midpoint(lo::T, hi::T) where T<:Integer = lo + ((hi - lo) >>> 0x01)
midpoint(lo::Integer, hi::Integer) = midpoint(promote(lo, hi)...)
struct HeapSortAlg  <: Base.Sort.Algorithm end
const DEFAULT_UNSTABLE = QuickSort
# const DEFAULT_STABLE   = MergeSort
const SMALL_ALGORITHM  = InsertionSort
const SMALL_THRESHOLD  = 20
const HeapSort = HeapSortAlg()

# Base on Julia sort

# Ordering in the following algorithms always is with by = identity
function ds_sort!(v, idx::Vector{<:Integer}, lo::Integer, hi::Integer, ::InsertionSortAlg, o::Ordering)
    @inbounds for i = lo+1:hi
        j = i
        x = v[i]
        y = idx[i]
        while j > lo
            if lt(o, x, v[j-1])
                v[j] = v[j-1]
                idx[j] = idx[j-1]
                j -= 1
                continue
            end
            break
        end
        v[j] = x
        idx[j] = y
    end
end

@inline function _selectpivot!(v, idx::Vector{<:Integer}, lo::Integer, hi::Integer, o::Ordering)
    @inbounds begin
        mi = midpoint(lo, hi)

        # sort v[mi] <= v[lo] <= v[hi] such that the pivot is immediately in place
        if lt(o, v[lo], v[mi])
            v[mi], v[lo] = v[lo], v[mi]
            idx[mi], idx[lo] = idx[lo], idx[mi]
        end

        if lt(o, v[hi], v[lo])
            if lt(o, v[hi], v[mi])
                v[hi], v[lo], v[mi] = v[lo], v[mi], v[hi]
                idx[hi], idx[lo], idx[mi] = idx[lo], idx[mi], idx[hi]
            else
                v[hi], v[lo] = v[lo], v[hi]
                idx[hi], idx[lo] = idx[lo], idx[hi]
            end
        end

        # return the location of pivot
        return lo
    end
end

function _partition!(v, idx::Vector{<:Integer}, lo::Integer, hi::Integer, o::Ordering)
    pv = _selectpivot!(v, idx, lo, hi, o)
    pivot = v[pv]
    pidx = idx[pv]
    # pivot == v[lo], v[hi] > pivot
    i, j = lo, hi
    @inbounds while true
        i += 1; j -= 1
        while lt(o, v[i], pivot); i += 1; end;
        while lt(o, pivot, v[j]); j -= 1; end;
        i >= j && break
        v[i], v[j] = v[j], v[i]
        idx[i], idx[j] = idx[j], idx[i]
    end
    v[j], v[lo] = pivot, v[j]
    idx[j], idx[lo] = pidx, idx[j]

    # v[j] == pivot
    # v[k] >= pivot for k > j
    # v[i] <= pivot for i < j
    return j
end

function ds_sort!(v, idx::Vector{<:Integer}, lo::Integer, hi::Integer, a::QuickSortAlg, o::Ordering)
    @inbounds while lo < hi
        hi-lo <= SMALL_THRESHOLD && return ds_sort!(v, idx, lo, hi, SMALL_ALGORITHM, o)
        j = _partition!(v, idx, lo, hi, o)
        if j-lo < hi-j
            # recurse on the smaller chunk
            # this is necessary to preserve O(log(n))
            # stack space in the worst case (rather than O(n))
            lo < (j-1) && ds_sort!(v, idx, lo, j-1, a, o)
            lo = j+1
        else
            j+1 < hi && ds_sort!(v, idx, j+1, hi, a, o)
            hi = j-1
        end
    end
end

# simple parallel implementation of QuickSort
# the assumption is that x[lo:mid] is sorted and x[mid+1:hi] is also sorted,
# the function uses this information to sort x[lo:hi]
# x_cpy is a copy of the x, idx_cpy is a copy of idx
function _sort_two_sorted_half!(x, x_cpy, idx::Vector{<:Integer}, idx_cpy, lo, mid, hi, o; cpy_offset = 0)
    st1 = lo
    en1 = mid
    st2 = mid+1
    en2 = hi
    cnt = lo
    @inbounds while true
        if lt(o, x_cpy[st1-cpy_offset], x_cpy[st2-cpy_offset])
            x[cnt] = x_cpy[st1-cpy_offset]
            idx[cnt] = idx_cpy[st1-cpy_offset]
            st1 += 1
            cnt += 1
            st1 > en1 && break
        else
            x[cnt] = x_cpy[st2-cpy_offset]
            idx[cnt] = idx_cpy[st2-cpy_offset]
            st2 += 1
            cnt += 1
            st2 > en2 && break
        end
    end
    @inbounds if st1 > en1
        while cnt <= hi
            x[cnt] = x_cpy[st2-cpy_offset]
            idx[cnt] = idx_cpy[st2-cpy_offset]
            st2 += 1
            cnt += 1
        end
    elseif st2 > en2
        while cnt <= hi
            x[cnt] = x_cpy[st1-cpy_offset]
            idx[cnt] = idx_cpy[st1-cpy_offset]
            st1 += 1
            cnt += 1
        end
    end
end

# to simplify the problem we assume number_of_chunks is 2^n for some n
function _sort_chunks!(x, idx::Vector{<:Integer}, lo, hi, number_of_chunks, a::Base.Sort.Algorithm, o::Ordering)
    rangelen = hi - lo + 1
    st_offset = lo - 1
    cz = div(rangelen, number_of_chunks)
    en = hi
    Threads.@threads for i in 1:number_of_chunks
        ds_sort!(x, idx, (i-1)*cz+1+st_offset, i*cz+st_offset, a, o)
    end
    # take care of the last few observations
    if number_of_chunks*div(rangelen, number_of_chunks)+st_offset < en
        ds_sort!(x, idx, number_of_chunks*div(rangelen, number_of_chunks)+1+st_offset, en, a, o)
    end
end

function _sort_multi_sorted_chunk!(x, idx::Vector{<:Integer}, lo, hi, number_of_chunks, a::Base.Sort.Algorithm, o::Ordering)
    rangelen = hi - lo + 1
    st_offset = lo - 1
    cz = div(rangelen, number_of_chunks)
    en = hi
    current_numberof_chunks = number_of_chunks
    x_cpy = x[lo:hi]
    idx_cpy = idx[lo:hi]
    while true
        Threads.@threads for i in 1:2:current_numberof_chunks
            _sort_two_sorted_half!(x, x_cpy, idx, idx_cpy, (i-1)*cz+1+st_offset, i*cz+st_offset, (i+1)*cz+st_offset, o; cpy_offset = lo-1)
        end
        cz *= 2
        current_numberof_chunks = current_numberof_chunks  >> 1
        current_numberof_chunks < 2 && break
        copyto!(x_cpy, 1, x, lo, rangelen)
        copyto!(idx_cpy, 1, idx, lo, rangelen)
    end
    # take care of the last few (less than number_of_chunks) observations
    if number_of_chunks*div(rangelen, number_of_chunks)+st_offset < en
        copyto!(x_cpy, 1, x, lo, rangelen)
        copyto!(idx_cpy, 1, idx, lo, rangelen)
        _sort_two_sorted_half!(x, x_cpy, idx, idx_cpy, lo, number_of_chunks*div(rangelen, number_of_chunks)+st_offset, en, o; cpy_offset = lo-1)
    end
end

# sorting a vector using parallel quick sort
# it uses a simple algorithm for doing this, and to make it even simpler the number of threads must be in the form of 2^n
function hp_ds_sort!(x, idx, a::Base.Sort.Algorithm, o::Ordering; lo = 1, hi = length(x))
    cpucnt = Threads.nthreads()
    @assert cpucnt >= 2 "we need at least 2 cpus for parallel sorting"
    cpucnt = 2 ^ floor(Int, log2(cpucnt))
    _sort_chunks!(x , idx, lo, hi, cpucnt, a, o)
    _sort_multi_sorted_chunk!(x, idx, lo, hi, cpucnt, a, o)
end

# Heapsort
# modified from DataStructures.jl, SortingAlgorithms.jl
# Binary heap indexing
heapleft(i::Integer) = 2i
heapright(i::Integer) = 2i + 1
heapparent(i::Integer) = div(i, 2)

# Binary min-heap percolate down.
function percolate_down!(xs::AbstractArray, idx, i::Integer, x=xs[i], idval = idx[i], o::Ordering=Forward, len::Integer=length(xs))
    @inbounds while (l = heapleft(i)) <= len
        r = heapright(i)
        j = r > len || lt(o, xs[l], xs[r]) ? l : r
        if lt(o, xs[j], x)
            xs[i] = xs[j]
            idx[i] = idx[j]
            i = j
        else
            break
        end
    end
    xs[i] = x
    idx[i] = idval
end

percolate_down!(xs::AbstractArray, idx, i::Integer, o::Ordering, len::Integer=length(xs)) = percolate_down!(xs, idx, i, xs[i], idx[i], o, len)

function heapify!(xs::AbstractArray, idx, o::Ordering=Forward)
    for i in heapparent(length(xs)):-1:1
        percolate_down!(xs, idx, i, o)
    end
    return xs
end

function ds_sort!(v::AbstractVector, idx::AbstractVector{<:Integer}, lo::Integer, hi::Integer, a::HeapSortAlg, o::Ordering)
    hi-lo <= SMALL_THRESHOLD && return ds_sort!(v, idx, lo, hi, SMALL_ALGORITHM, o)
    if lo > 1 || hi < length(v)
        return ds_sort!(view(v, lo:hi), view(idx, lo:hi), 1, length(v), a, o)
    end
    r = ReverseOrdering(o)
    heapify!(v, idx, r)
    @inbounds for i = length(v):-1:2
        # Swap the root with i, the last unsorted position
        x = v[i]
        idxval = idx[i]
        v[i] = v[1]
        idx[i] = idx[1]
        # The heap portion now ends at position i-1, but needs fixing up
        # starting with the root
        percolate_down!(v, idx, 1, x, idxval, r, i-1)
    end
    v
end

function percolate_down2!(xs::AbstractArray, i::Integer, x=xs[i], o::Ordering=Forward, len::Integer=length(xs))
    @inbounds while (l = heapleft(i)) <= len
        r = heapright(i)
        j = r > len || lt(o, xs[l], xs[r]) ? l : r
        if lt(o, xs[j], x)
            xs[i] = xs[j]
            i = j
        else
            break
        end
    end
    xs[i] = x
end

percolate_down2!(xs::AbstractArray, i::Integer, o::Ordering, len::Integer=length(xs)) = percolate_down2!(xs, i, xs[i], o, len)

function heapify2!(xs::AbstractArray, o::Ordering=Forward)
    for i in heapparent(length(xs)):-1:1
        percolate_down2!(xs, i, o)
    end
    return xs
end

function Base.sort!(v::AbstractVector, lo::Integer, hi::Integer, a::HeapSortAlg, o::Ordering = Forward)
    hi-lo <= SMALL_THRESHOLD && return sort!(v, lo, hi, SMALL_ALGORITHM, o)
    if lo > 1 || hi < length(v)
        return sort!(view(v, lo:hi), 1, length(v), a, o)
    end
    r = ReverseOrdering(o)
    heapify2!(v, r)
    @inbounds for i = length(v):-1:2
        # Swap the root with i, the last unsorted position
        x = v[i]
        v[i] = v[1]
        # The heap portion now ends at position i-1, but needs fixing up
        # starting with the root
        percolate_down2!(v, 1, x, r, i-1)
    end
    v
end


text:
module CSVanalyzer

import DataFrames: DataFrame
import HypothesisTests: pvalue, SignedRankTest
import Printf.@printf
import Statistics: mean, median, std
import DelimitedFiles.readdlm

readcsv(fname) = readdlm(fname, ',', '\n')

export statitistic,printSummary,statsToLatex,getComparison,comparisonToLatex

function statitistic(M::Matrix{Float64}; mapping::Function=identity)
    # rows = funs
    # cols = runs

    result = zeros(size(M,1), 5)
    for i = 1:size(M,1)
        x = mapping(M[i,:])

        result[i, 1] = minimum(x)
        result[i, 2] = median(x)
        result[i, 3] = mean(x)
        result[i, 4] = maximum(x)
        result[i, 5] = std(x)
    end

    return result
end

function printSummary(stats::Matrix{Float64}; mapping::Function=identity)
    for i = 1:size(stats,1)
        @printf("Best = %.4e  Median = %.4e  Mean = %.4e  Worst = %.4e  std = %.4e\n",
                    stats[i, 1],
                    stats[i, 2],
                    stats[i, 3],
                    stats[i, 4],
                    stats[i, 5])
    end
end

function statsToLatex(stats::Matrix{Float64}; mapping::Function=identity)

    println("fn & Best &  Median &  Mean &  Worst &  Std. \\\\ \n")

    for i = 1:size(stats,1)
        @printf("%d & %.4e & %.4e & %.4e & %.4e & %.4e \\\\ \\hline \n", i,
                    stats[i, 1],
                    stats[i, 2],
                    stats[i, 3],
                    stats[i, 4],
                    stats[i, 5])
    end
end

function getComparison(algorithms::Dict{String, String}, algName::String)

    outputTable = DataFrame()

    compTable = readcsv(algorithms[algName])

    i = 1
    for name = keys(algorithms)
        thetable = readcsv(algorithms[name])
        stats =  statitistic(thetable)

        outputTable[Symbol("$(name)_mean")] =  stats[ :, 3 ]
        # outputTable[Symbol("$(name)_std")]  =  stats[ :, 5 ]

        if name == algName
            continue
        end

        res = String[]

        for i = 1:size(compTable,1)
            x, y = compTable[i,:], thetable[i,:]

            p = pvalue(SignedRankTest(x, y))

            if p < 0.05
                mx, my = mean(x), mean(y)
                if mx < my
                    push!(res, "+")
                elseif mx == my
                    push!(res, "≈")
                else
                    push!(res, "-")
                end

            else
                push!(res, "≈")
            end
        end

        outputTable[Symbol("$(algName)_vs_$(name)")]  =  res

        i += 1

    end

    return outputTable
end

function comparisonToLatex(algorithms::Dict{String, String}, algName::String)

    thetable = getComparison(algorithms, algName)

    print("fn & ")
    for n = names(thetable)
        print(string(n), " & ")
    end

    println("")

    for i = 1:size(thetable, 1)
        print(i, " & ")
        for j = 1:size(thetable, 2)
            d = thetable[i, j]

            if typeof(d) == String

                if d == "≈"
                    d = "\\approx"
                end

                print("\$$d\$")

                if j == size(thetable, 2)
                    println("\\\\ \\hline")
                else
                    print(" & ")
                end
                continue
            end

            @printf(" %.4e & ", d)
        end

    end
end

statitistic(fname::String; mapping::Function=identity) = statitistic( readcsv(fname); mapping=mapping )
printSummary(fname::String; mapping::Function=identity) = printSummary( statitistic(readcsv(fname)); mapping=mapping )
statsToLatex(fname::String; mapping::Function=identity) = statsToLatex(statitistic( readcsv(fname); mapping=mapping ))
printComparison(algorithms::Dict{String, String}, algName::String) = println(getComparison(D))

end # module


text:
"""
    spatialsubset(C::ClimGrid, poly::Array{N, 2})

Returns the spatial subset of ClimGrid C. The spatial subset is defined by the polygon poly, defined on a -180, +180 longitude reference.

"""
function spatialsubset(C::ClimGrid, poly)

    # Some checks for polygon poly

    if size(poly, 1) != 2 && size(poly, 2) == 2
        # transpose
        poly = poly'
    end

    lonsymbol = Symbol(C.dimension_dict["lon"])
    latsymbol = Symbol(C.dimension_dict["lat"])

    longrid = C.longrid
    latgrid = C.latgrid

    lon = C[1][Axis{lonsymbol}][:]
    lat = C[1][Axis{latsymbol}][:]

    msk = inpolygrid(longrid, latgrid, poly)

    begin
        I = findall(!isnan, msk)
        idlon, idlat = (getindex.(I, 1), getindex.(I, 2))
    end
    minXgrid = minimum(idlon)
    maxXgrid = maximum(idlon)
    minYgrid = minimum(idlat)
    maxYgrid = maximum(idlat)

    # Get DATA
    data = C[1].data

    if ndims(data) == 3
        data = data[minXgrid:maxXgrid, minYgrid:maxYgrid, :]
    elseif ndims(data) == 4
        data = data[minXgrid:maxXgrid, minYgrid:maxYgrid, :, :]
    end

    #new mask (e.g. representing the region of the polygon)
    msk = msk[minXgrid:maxXgrid, minYgrid:maxYgrid]
    data = applymask(data, msk)

    # Get lon-lat for such region
    longrid = longrid[minXgrid:maxXgrid, minYgrid:maxYgrid]
    latgrid = latgrid[minXgrid:maxXgrid, minYgrid:maxYgrid]
    lon = lon[minXgrid:maxXgrid]
    lat = lat[minYgrid:maxYgrid]

    if ndims(data) == 3
      dataOut = AxisArray(data, Axis{lonsymbol}(lon), Axis{latsymbol}(lat),  Axis{:time}(C[1][Axis{:time}][:]))
    elseif ndims(data) == 4
      dataOut = AxisArray(data, Axis{lonsymbol}(lon), Axis{latsymbol}(lat), Axis{:level}(C[1][Axis{:plev}][:]), Axis{:time}(C[1][Axis{:time}][:]))
    end

    return ClimGrid(dataOut, longrid=longrid, latgrid=latgrid, msk=msk, grid_mapping=C.grid_mapping, timeattrib=C.timeattrib, dimension_dict=C.dimension_dict, model=C.model, frequency=C.frequency, experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable=C.variable, typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)

end

"""
    shapefile_coords_poly(poly::Shapefile.Polygon)

Return the polygons contained in shp.shapes[i]. It returns an array containing the polygons.

See also [`shapefile_coords`](@ref), which returns vectors as opposed to array. Returned polygon is consistent with the data extraction of the [`load`](@ref) function.

"""
function shapefile_coords_poly(poly::Shapefile.Polygon)
    x, y = shapefile_coords(poly)
    return [x y]'
end

"""
    extractpoly(file::String; n::Int=1)

Returns the n-th polygon contained in *file*.
"""
function extractpoly(file::String; n::Int=1)

    shp = open(file) do fd
        read(fd, Shapefile.Handle)
    end

    poly = shapefile_coords_poly(shp.shapes[n])
    return poly
end

"""
    shapefile_coords(poly::Shapefile.Polygon)

This function return the polygons contained in shp.shapes[i]. It returns the x and y coordinates vectors.

See also [`shapefile_coords_poly`](@ref), which returns a polygon that ca be used for data extraction of the [`load`](@ref).

"""
function shapefile_coords(poly::Shapefile.Polygon)
    start_indices = poly.parts .+ 1
    end_indices = vcat(poly.parts[2:end], length(poly.points))
    x, y = zeros(0), zeros(0)
    for (si,ei) in zip(start_indices, end_indices)
        push!(x, NaN)
        push!(y, NaN)
        for pt in poly.points[si:ei]
            push!(x, pt.x)
            push!(y, pt.y)
        end
    end
    x, y
end

function shiftgrid_180_west_east(longrid)

    for ilon in eachindex(longrid)
        if longrid[ilon] >= 180
            longrid[ilon] = longrid[ilon] - 360
        end
    end

    ieast = longrid .>= 0
    iwest = longrid .< 0

    grideast = reshape(longrid[ieast], :, size(longrid, 2))
    gridwest = reshape(longrid[iwest], :, size(longrid, 2))

    longrid_flip = vcat(gridwest, grideast)

    return longrid_flip

end

function shiftgrid_180_east_west(longrid)

    for ilon in eachindex(longrid)
        if longrid[ilon] >= 180
            longrid[ilon] = longrid[ilon] - 360
        end
    end

    ieast = longrid .>= 0
    iwest = longrid .< 0

    grideast = reshape(longrid[ieast], :, size(longrid, 2))
    gridwest = reshape(longrid[iwest], :, size(longrid, 2))

    longrid_flip = vcat(grideast, gridwest)

    return longrid_flip

end

function shiftarray_west_east(data, longrid_flip)

    ieast = longrid_flip .>= 0
    iwest = longrid_flip .< 0

    msk = permute_west_east2D(data, iwest, ieast)
    return msk

end

function shiftarray_east_west(data, longrid_flip)

    ieast = longrid_flip .>= 0
    iwest = longrid_flip .< 0

    msk = permute_east_west2D(data, iwest, ieast)
    return msk

end

function shiftvector_180_east_west(lon_raw)

    for ilon in eachindex(lon_raw)
        if lon_raw[ilon] >= 180
            lon_raw[ilon] = lon_raw[ilon] - 360
        end
    end
    ieast_vec = lon_raw .>= 0
    iwest_vec = lon_raw .< 0
    lon_raw_flip = [lon_raw[ieast_vec];lon_raw[iwest_vec]]
    return lon_raw_flip

end

function shiftvector_180_west_east(lon_raw)

    for ilon in eachindex(lon_raw)
        if lon_raw[ilon] >= 180
            lon_raw[ilon] = lon_raw[ilon] - 360
        end
    end

    ieast_vec = lon_raw .>= 0
    iwest_vec = lon_raw .< 0

    lon_raw_flip = [lon_raw[iwest_vec];lon_raw[ieast_vec]]

    return lon_raw_flip

end

function meridian_check(poly)

    minpoly = minimum(poly[1, .!isnan.(poly[1, :])])
    maxpoly = maximum(poly[1, .!isnan.(poly[1, :])])

    meridian = false
    if sign(minpoly)*sign(maxpoly) == sign(-1.0) # polygon crosses the meridian
       meridian = true
    end

    return meridian

end

function findmindist(p1::Tuple{Float64, Float64}, p2::Tuple{Array{Float64,1},Array{Float64,1}}; R=6372.8)

    dist = Array{Float64}(undef, length(p2[1]), 1)

    for idx = 1:length(p2[1])
        p3 = ( p2[1][idx], p2[2][idx] )
        dist[idx] = haversine(p1, p3, 6372.8)
    end

    return findmin(dist)

end


text:
export LocalHamiltonian, LHIterator

mutable struct LocalHamiltonian{B, R, M<:AbstractSparseMatrix}
    data::M
    nzbase::Vector{B}
end

function LocalHamiltonian(::Type{B}, ::Type{R}, mat::M) where {B, R, M}
    nzbase = [convert(B, i-1) for i in mat.rowval]
    LocalHamiltonian{B, R, M}(mat, nzbase)
end

LocalHamiltonian(::Type{T}, ::Type{R}, mat::AbstractSparseMatrix) where {T <: SiteLabel, R} =
    LocalHamiltonian(SubSites{T, eltype(T), Int(log2(size(mat, 1)))}, R, mat)
LocalHamiltonian(::Type{R}, mat::AbstractSparseMatrix) where {R <: AbstractRegion} =
    LocalHamiltonian(Bit, R, mat)
LocalHamiltonian(::Type{T}, mat::AbstractSparseMatrix) where {T <: SiteLabel} =
    LocalHamiltonian(Bit, Region{1}, mat)
LocalHamiltonian(mat::AbstractSparseMatrix) =
    LocalHamiltonian(Bit, Region{1}, mat)

# LocalHamiltonian shares some similar interface
# with SparseMatrix, for convinience.
# But it is not a subtype of AbstractArray
import Base: eltype, length, ndims, size, eachindex, stride, getindex
import Compat: axes
eltype(h::LocalHamiltonian) = eltype(h.data)
length(h::LocalHamiltonian) = length(h.data)
ndims(h::LocalHamiltonian) = ndims(h.data)
size(h::LocalHamiltonian) = size(h.data)
size(h::LocalHamiltonian, n::Integer) = size(h.data, n)
axes(h::LocalHamiltonian) = axes(h.data)
axes(h::LocalHamiltonian, d::Integer) = axes(h.data)
eachindex(h::LocalHamiltonian) = eachindex(h.data)
stride(h::LocalHamiltonian, k::Integer) = stride(h.data, k)
strides(x::LocalHamiltonian) = strides(x.data)

# TODO: support indexing with Sites
getindex(x::LocalHamiltonian, index::Integer...) = getindex(x.data, index...)
setindex!(x::LocalHamiltonian, val, index::Integer...) = setindex!(x.data, val, index...)

import Base: show

function show(io::IO, h::LocalHamiltonian{B, R, M}) where {B, R, M}
    print(io, "LocalHamiltonian{$B, $R, $M}")
    print(io, h.data)
end

mutable struct LHIterator{B, R, M} <: AbstractHamiltonianIterator
    data::M
    nzbase::Vector{B}
    index::Int
end

LHIterator(::Type{R}, data::M, nzbase::Vector{B}, index::Int) where {R, M, B} =
    LHIterator{B, R, M}(data, nzbase, index)

LHIterator(h::LocalHamiltonian{B, R, M}, rhs::SubSites) where {B, R, M} =
    LHIterator(R, h.data, h.nzbase, convert(Int, rhs) + 1)
LHIterator(h::LocalHamiltonian{B}, rhs::Number...) where B =
    LHIterator(h, B(rhs))

import Base: start, next, done, eltype, length

start(itr::LHIterator) = itr.data.colptr[itr.index]
next(itr::LHIterator, state) = (itr.data.nzval[state], itr.nzbase[state]), state+1
done(itr::LHIterator, state) = state >= itr.data.colptr[itr.index+1]

eltype(itr::LHIterator{Ti, B}) where {Ti, B} = Tuple{eltype(itr.data), B}
length(itr::LHIterator{Ti, B}) where {Ti, B} =
    itr.data.colptr[itr.index] - itr.data.colptr[itr.index]


text:
using DSGE
using HDF5, Base.Test

path = dirname(@__FILE__)

h5 = h5open("$path/../reference/posterior.h5")
data = read(h5, "data")
lh_expected = read(h5, "lnpy")
post_expected = read(h5, "obj")
close(h5)

m = Model990()

lh, _ = likelihood(m, data)
@test_approx_eq lh_expected lh

post = posterior(m, data)[:post]
@test_approx_eq post_expected post

x = map(α->α.value, m.parameters)
post_at_start = posterior!(m, x, data)[:post]
@test_approx_eq post_expected post_at_start

# Ensure if we are not evaluating at start vector, then we do not get the reference
# posterior
y = x .+ 0.01
post_not_at_start = posterior!(m, y, data)[:post]
ϵ = 1.0
@test abs(post_at_start - post_not_at_start) > ϵ


text:
using Cxx

cxx"""class Bill {
         public:
            int64_t a;
            int64_t b;
       };"""

 type Bill
          a::Int
          b::Int
       end

s = Bill(1, 2)

cxx"""#include <stdio.h>"""

cxx"""void myfn(Bill * m) { printf("%ld, %ld\n", m->a, m->b);}"""

cxxptr = @cxxnew Bill()
@assert isa(cxxptr, CppPtr)
@cxx myfn(cxxptr)

@cxx myfn(jpcpp"Bill"(s))


text:
using StructJuMP

######
# DATA
######
NUMSCEN = 3
TESTTIME = GEN = WIND = BUS = LIN = GEN = LOAD = 1:2
NODES = HORIZON = [1]
FIXEDGENCOST = THETASCALE = COSTSCALE = 1
Pmax = zeros(2,2)
np_cap = wind_share = X = V = max_ur = Pgen_init = gen_cost = ones(2)
rec_bus = snd_bus = bus_gen = bus_load = bus_wind = 1:2
wind_total = [1 1]
node = NODES[1]
Pload = ones(2,2)
ref_bus = 1

m = StructuredModel()

@variable(m, dummyV >= 0)

s0 = StructuredModel(parent=m, id=1)
@variable(s0, 0 <= y[TESTTIME,GEN] <= 1)
@objective(s0, :Min, FIXEDGENCOST*sum(y[t,j] for t=TESTTIME,j=GEN))

for it=1:NUMSCEN
    st = StructuredModel(parent=s0, id=1+it)
    @variable(st, 0 <= Pgen[TESTTIME,j=GEN] <= np_cap[j])
    @variable(st, 0 <= Pwind[t=TESTTIME,j=WIND] <= wind_total[node,t]*wind_share[j])
    @variable(st, -THETASCALE*π/2 <= theta[TESTTIME,j=BUS] <= THETASCALE*π/2)
    @variable(st, -Pmax[i] <= P[TESTTIME,i=LIN] <= Pmax[i])
    for t in TESTTIME
        for j in BUS
            @constraint(st, ( sum(P[t,i] for i=LIN if j==rec_bus[i])
                                -sum(P[t,i] for i=LIN if j==snd_bus[i])
                                +sum(Pgen[t,i] for i=GEN if j==bus_gen[i])
                                -sum(Pload[t,i] for i=LOAD if j==bus_load[i])
                                +sum(wind_total[node,t]*wind_share[i]-Pwind[t,i] for i=WIND if j==bus_wind[i])
                                == 0))
        end
        @constraint(st, theta[t,ref_bus] == 0)
        for i in LIN
            @constraint(st, X[i]*P[t,i] - V[i]^2*(theta[t,snd_bus[i]]-theta[t,rec_bus[i]])/THETASCALE == 0)
        end
    end
    # ramp constraints
    for t in HORIZON
        for i in GEN
            @constraint(st, -max_ur[i] <= Pgen[t,i] - Pgen[t+1,i] <= max_ur[i]) # range constraint, beware...
        end
    end

    for i in GEN
        @constraint(st, -max_ur[i] <= Pgen[1,i] - Pgen_init[i] <= max_ur[i])
    end

    # linking constraints
    # gen[s,t,j] <= np_cap[j]*y[t,j]
    for t in TESTTIME
        for j in GEN
            @constraint(st, -np_cap[j] <= Pgen[t,j] - np_cap[j]*y[t,j] <= 0) # should wrap the call to y in an ancestor-type thing
        end
    end

    # need Exp here for earlier sml versions
    @objective(st, :Min, sum(COSTSCALE*gen_cost[i]*Pgen[t,i] for t=TESTTIME, i=GEN))

end


text:
using SHIPs, JuLIP, BenchmarkTools
using JuLIP: alloc_temp, alloc_temp_d

trans = PolyTransform(2, 1.0)
fcut = PolyCutoff2s(Val(2), 0.5, 3.0)
ships = [SHIPBasis(SparseSHIP(n, 15, 1.5), trans, fcut)
         for n = 2:4]

NR = 50
Rs = 1.0 .+ 2*(rand(JVecF, NR) .- 0.5)
Zs = zeros(Int16, NR)
z0 = 0

@info("profile `precompute_A!`")
tmp = alloc_temp(ships[1], NR)
@btime SHIPs.precompute_A!($tmp, $(ships[1]), $Rs, $Zs)
@btime SHIPs.precompute_A!($tmp, $(ships[1]), $Rs, $Zs)
@info("profile `precompute_grads!`")
tmpd = alloc_temp_d(ships[1], NR)
@btime SHIPs.precompute_grads!($tmpd, $(ships[1]), $Rs, $Zs)
@btime SHIPs.precompute_grads!($tmpd, $(ships[1]), $Rs, $Zs)

@info("profile basis computation")
for n = 2:4
   @info("  body-order $(n+1):")
   🚢 = ships[n-1]
   B = SHIPs.alloc_B(🚢)
   @info("     eval_basis:")
   @btime SHIPs.eval_basis!($B, $tmp, $🚢, $Rs, $Zs, $z0)
   @btime SHIPs.eval_basis!($B, $tmp, $🚢, $Rs, $Zs, $z0)
   @info("     eval_basis_d:")
   dB = SHIPs.alloc_dB(🚢, Rs)
   @btime SHIPs.eval_basis_d!($B, $dB, $tmpd, $🚢, $Rs, $Zs, $z0)
   @btime SHIPs.eval_basis_d!($B, $dB, $tmpd, $🚢, $Rs, $Zs, $z0)
end

# ##
# using Profile
# 🚢 = ships[2]
# B = SHIPs.alloc_B(🚢)
# dB = SHIPs.alloc_dB(🚢, Rs)
# @btime SHIPs.eval_basis!($B, $tmp, $🚢, $Rs, $Zs, $z0)
# @btime SHIPs.eval_basis_d!($B, $dB, $tmpd, $🚢, $Rs, $Zs, $z0)
#
# @code_warntype SHIPs._eval_basis!(B, tmp, 🚢, Val{3}(), 1, 🚢.NuZ[3,1])
#
# SHIPs._eval_basis!(B, tmp, 🚢, Val{3}(), 1, 🚢.NuZ[3,1])
#
# ##
#
# function runn(N, f, args...)
#    for n = 1:N
#       f(args...)
#    end
# end
#
# runn(10, SHIPs.eval_basis!, B, tmp, 🚢, Rs, Zs, z0)
# runn(10, SHIPs.eval_basis_d!, B, dB, tmpd, 🚢, Rs, Zs, z0)
#
# ##
#
# Profile.clear()
# @profile runn(100, SHIPs.eval_basis!, B, tmp, 🚢, Rs, Zs, z0)
# Profile.print()

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
# -*- coding: utf-8 -*-
# ![logo](https://julialang.org/v2/img/logo.svg)

# # Installation
#
# - Download for your plaform at http://julialang.org/downloads/
#
# - You can also try online via https://www.juliabox.com/ which provides hosted Julia notebooks on Google Cloud.
#
# - If you need corporate support, https://juliacomputing.com/products/juliapro.html provides a _batteries included_ installation of Julia. The company is run by some of the creators of the language and provides a single install with all of the most important add-ons.
#
# - Most Linux distributions come with Julia packages in their repositories. However, these may lag somewhat behind the current rather fast development cycle.
#
# - Notebooks are available on https://github.com/pnavaro/math-julia

# ## Installing Julia on Linux
#
# Download the tar archive from https://julialang.org/downloads and rename the directory like so:
# ```
# cd ~
# tar -xzf Downloads/julia-1.0.1-linux-x86_64.tar.gz
# ```
#
# You may also want to add `julia` to your path
#
# ```
# export PATH=$PATH:~/julia-1.0.1/bin
# echo 'export PATH=$PATH:~/julia-1.0.1/bin' >>~/.bash_profile
# ```

# # Julia shell
# - Type julia in a terminal
# ![julia shell](images/julia_prompt.png)
#
#
#

# - "Tab" for autocompletion
# - pwd() (Current directory)
# - cd() (Change directory),
# - homedir() (Home directory)
# - Use ; to switch to unix shell.
# - help("sin") ou ?sin toggles the help mode.

# ## Run julia code.
#
# Julia programs use ".jl" extension by convention and can be executed in the julia prompt with:
# ```julia
# include("my_program.jl")
# ```
# - Julia is first translated into an intermediate representation.
# - Then LLVM compiles it for your machine.
#
# This means that
# - re-running the same code is faster the second time around
# - it runs at speeds comparable to compiled C or Fortran code

# # IDE
#
# - Julia shell
# - Command line : `julia my_program.jl`
# - Julia notebooks (Jupyter)
# - No dedicated IDE but plugins for existing tools.
# - IntellJ, Eclipse, PyCharm
# - [VS Code](https://github.com/JuliaEditorSupport/julia-vscode)
# - Atom and [Juno](http://junolab.org/)
# - [Sublime Text 3](https://github.com/quinnj/Sublime-IJulia)
#
# I personally use Jupyter for development but Atom seems to be the most used environment.

# # Packages
#
# - Use ] to switch to package manager.
#
# ```
# pkg> add IJulia
# ```
#
# will install the package. Type `help` to display all available commands.
#
# Use "backspace" to go back to the `julia` shell.
#
# https://pkg.julialang.org
#
# The package is installed in directory `~/.julia/`
#
# To import the package, type:
# ```julia
# using IJulia
# ```

# # Jupyter notebook
#
# The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text. Uses include: data cleaning and transformation, numerical simulation, statistical modeling, data visualization, machine learning, and much more.
#
# type:
# ```julia
# using IJulia
# notebook()  # use notebook(detached=true)
# ```
#
# To convert a notebook file `.ipynb` in julia program `.jl`:
# * In the top menu `File->Download as`
# * Command line : ```ipython nbconvert --to script my_notebook.ipynb```.
#

# # Why Julia?
#
# - Increase the programmer productivity.
# - Python is a very nice language to code a prototype but i had to use Cython to make it faster. And had to
# code in Fortran if the code had to be deployed on HPC cluster
# - High-level languages like python and R let one explore and experiment rapidly, but can run slow.
# - Low-level languages like Fortran/C++ tend to take longer to develop, but run fast.
# - This is sometimes called the "two language problem" and is something the Julia developers set out to eliminate.
# - My code runs much faster than Python, and typically similar to Fortran.
#
# Julia provides a "best of both worlds" experience for programmers who need to develop novel algorithms and bring them into production environments with minimal effort.
#

# # Julia features
#
# - High-level language for numerical computing.
# - Born in 2009 and version 1.0 was released in August 2018.
# - Dynamically-typed with optional types, feels like a scripting language, and has good support for interactive use.
# - Easy to learn for people that comes from R and Python.
# - Designed from the beginning to be fast as Fortran and C
# - Julia programs compile to efficient native code via [LLVM](https://llvm.org).
# - Designed for parallelism, and provides built-in primitives for parallel computing
# - can call C and Fortran libraries natively
# - can call Python libraries via PyCall package
# - can call R functions via RCall package

# # Julia is a language made for Science.
#
# http://www.stochasticlifestyle.com/some-state-of-the-art-packages-in-julia-v1-0
#
# Mathematics
# * JuliaDiff – Differentiation tools
# * JuliaDiffEq – Differential equation solving and analysis
# * JuliaGeometry – Computational Geometry
# * JuliaGraphs – Graph Theory and Implementation
# * JuliaIntervals - Rigorous numerics with interval arithmetic & applications
# * JuliaMath – Mathematics made easy in Julia
# * JuliaOpt – Optimization (Gitter)
# * JuliaPolyhedra – Polyhedral computation
# * JuliaSparse – Sparse matrix solvers
#
# Data Science
# * JuliaML – Machine Learning
# * JuliaStats – Statistics
# * JuliaImages – Image Processing
# * JuliaText – Natural Language Processing (NLP), Computational Linguistics and (textual) Information Retrieval
# * JuliaDatabases – Various database drivers for Julia
# * JuliaData – Data manipulation, storage, and I/O in Julia
#
#

max_stars_repo_name: pnavaro/math-julia

max_issues_repo_name: pnavaro/math-julia

max_forks_repo_name: pnavaro/math-julia


text:
# ====================
# Testing the partial regex implementation in XUnit
module RegexTests
using Test  # For testing the behavior
import XUnit  # The behavior to test

@testset "partial regexes" begin
    date_re = XUnit.partial(raw"^\d\d-\d\d-\d\d\d\d")  # DD-MM-YYYY Date specifier

    # Full match
    @test XUnit.pmatch(date_re, "22-02-20") !== nothing
    # Partial matches
    @test XUnit.pmatch(date_re, "22-02") !== nothing
    @test XUnit.pmatch(date_re, "22-") !== nothing
    # Non-matches
    @test XUnit.pmatch(date_re, "-2") === nothing
    @test XUnit.pmatch(date_re, "-22-") === nothing
    @test XUnit.pmatch(date_re, "") === nothing
end
end
# ====================

using XUnit

function test()
    @testset "t" begin
        for ch in ["a","b","c"]
            @testset "$ch" begin
                @test 1==1
                for i in 1:3
                    @testset "$i" begin
                        @test 1==1
                        @test_broken 0==1
                    end
                end
            end
        end
    end
end

function closure(expected::Vector{String})
    cl = Set{String}()
    for path in expected
        segs = split(path, "/")
        s = ""
        first = true
        for seg in segs
            if first
                s = seg
                first = false
            else
                s = s*"/"*seg
            end
            push!(cl, s)
        end
    end
    sort(collect(cl))
end

function case(str::String, test::Function, depth::Int64,
        expected::Vector{String}, args...)
    println(str)
    println("======")
    (_, state) = XUnit.runtests_return_state(test, depth, args...)
    println()
    seen = sort(collect(keys(filter(kv -> kv.second, state.seen))))
    expected = sort(closure(expected))
    @test seen == expected
    state
end

case(str::String, test::Function, expected::Vector{String}, args...) =
    case(str, test, typemax(Int64), expected, args...)

@testset "XUnit Filter" begin
    case("Run all tests", test,
        [ "t/a/1", "t/a/2", "t/a/3",
          "t/b/1", "t/b/2", "t/b/3",
          "t/c/1", "t/c/2", "t/c/3" ]
    )

    case("Run everything under 't/a'", test,
        [ "t/a/1", "t/a/2", "t/a/3" ],
         "t/a/.*"
    )

    # "Run 'a/2'"
    case("Run 't/a/2'", test,
        [ "t/a/2" ],
         "t/a/2")

    # "Run 'a/2' and 'b/3'"
    case("Run 't/a/2' and 't/b/3'", test,
        [ "t/a/2", "t/b/3" ],
        "t/a/2", "t/b/3")

    # "Run all except 'b/2' and 'b/3'"
    case("Run all except 't/b/2' and 't/b/3'", test,
        [ "t/a/1", "t/a/2", "t/a/3",
          "t/b/1",
          "t/c/1", "t/c/2", "t/c/3" ],
        "¬t/b/2", "-t/b/3")

    # "Print names of top-level test sets"
    case("Show top-level test sets", test, 2,
        [ "t/a", "t/b", "t/c" ]
    )

    # "Print names of top-level and second-level test sets"
    case("Print names of top-level and second-level test sets", test, 3,
        [ "t/a/1", "t/a/2", "t/a/3",
          "t/b/1", "t/b/2", "t/b/3",
          "t/c/1", "t/c/2", "t/c/3" ]
    )
end


text:
module Junctions

import Base: convert, getindex, reduce
import Base.Broadcast

export @all, @any,
    eigenstates,
    Junction

struct Shift{F<:Function}
    continuation::F
end

getindex(s::Shift) = s.continuation()

struct Junction{T}
    reducer
end

struct And end
struct Or end

function construct_junction(typ, values)
    expr = :(nothing)
    for v in Iterators.reverse(values)
        expr = quote
            let v = $v
                acc = op(acc, v)
                if cont(v)
                    $expr
                end
            end
        end
    end

    expr = quote
        acc = init
        $expr
        return acc
    end

    :(Junction{$typ}((init, op, cont) -> $expr))
end

macro all(values...)
    construct_junction(:And, values)
end

macro any(values...)
    construct_junction(:Or, values)
end

reduce(op, j::Junction, init) = j.reducer(init, op, _ -> true)

eigenstates(j::Junction, ::Type{T} = Any) where T = reduce(push!, j, T[])

convert(::Type{Bool}, j::Junction{And}) = j.reducer(true, &, identity)
convert(::Type{Bool}, j::Junction{Or}) = j.reducer(false, |, !)

# broadcasting
struct JunctionStyle <: Broadcast.BroadcastStyle end
Base.BroadcastStyle(::Type{<:Junction}) = JunctionStyle()
Broadcast.BroadcastStyle(s::JunctionStyle, ::Broadcast.BroadcastStyle) = s

Broadcast.broadcastable(x::Junction) = x

function Broadcast.broadcasted(::JunctionStyle, f, arg1::Junction{T}) where T
    function reducer(init, op, cont)
        arg1.reducer(init, (acc, v) -> op(acc, f(v)), v -> cont(f(v)))
    end
    Junction{T}(reducer)
end

end # module


text:
using Base.Test

@testset "No config default tests" begin
    nothing
end


text:
export StrongStabilityPreservingRungeKutta
export SSPRK22Heuns, SSPRK22Ralstons, SSPRK33ShuOsher, SSPRK34SpiteriRuuth

"""
    StrongStabilityPreservingRungeKutta(f, RKA, RKB, RKC, Q; dt, t0 = 0)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

The constructor builds a strong-stability-preserving Runge--Kutta scheme
based on the provided `RKA`, `RKB` and `RKC` coefficient arrays.

The available concrete implementations are:

  - [`SSPRK33ShuOsher`](@ref)
  - [`SSPRK34SpiteriRuuth`](@ref)
"""
mutable struct StrongStabilityPreservingRungeKutta{T, RT, AT, Nstages} <:
               AbstractODESolver
    "time step"
    dt::RT
    "time"
    t::RT
    "elapsed time steps"
    steps::Int
    "rhs function"
    rhs!
    "Storage for RHS during the `StrongStabilityPreservingRungeKutta` update"
    Rstage::AT
    "Storage for the stage state during the `StrongStabilityPreservingRungeKutta` update"
    Qstage::AT
    "RK coefficient vector A (rhs scaling)"
    RKA::Array{RT, 2}
    "RK coefficient vector B (rhs add in scaling)"
    RKB::Array{RT, 1}
    "RK coefficient vector C (time scaling)"
    RKC::Array{RT, 1}

    function StrongStabilityPreservingRungeKutta(
        rhs!,
        RKA,
        RKB,
        RKC,
        Q::AT;
        dt = 0,
        t0 = 0,
    ) where {AT <: AbstractArray}
        T = eltype(Q)
        RT = real(T)
        new{T, RT, AT, length(RKB)}(
            RT(dt),
            RT(t0),
            0,
            rhs!,
            similar(Q),
            similar(Q),
            RKA,
            RKB,
            RKC,
        )
    end
end

"""
    ODESolvers.dostep!(Q, ssp::StrongStabilityPreservingRungeKutta, p,
                       time::Real, [slow_δ, slow_rv_dQ, slow_scaling])

Use the strong stability preserving Runge--Kutta method `ssp` to step `Q`
forward in time from the current time `time` to final time `time + getdt(ssp)`.

If the optional parameter `slow_δ !== nothing` then `slow_rv_dQ * slow_δ` is
added as an additional ODE right-hand side source. If the optional parameter
`slow_scaling !== nothing` then after the final stage update the scaling
`slow_rv_dQ *= slow_scaling` is performed.
"""
function dostep!(
    Q,
    ssp::StrongStabilityPreservingRungeKutta,
    p,
    time,
    slow_δ = nothing,
    slow_rv_dQ = nothing,
    in_slow_scaling = nothing,
)
    dt = ssp.dt

    RKA, RKB, RKC = ssp.RKA, ssp.RKB, ssp.RKC
    rhs! = ssp.rhs!
    Rstage, Qstage = ssp.Rstage, ssp.Qstage

    rv_Q = realview(Q)
    rv_Rstage = realview(Rstage)
    rv_Qstage = realview(Qstage)
    groupsize = 256

    rv_Qstage .= rv_Q
    for s in 1:length(RKB)
        rhs!(Rstage, Qstage, p, time + RKC[s] * dt, increment = false)

        slow_scaling = nothing
        if s == length(RKB)
            slow_scaling = in_slow_scaling
        end
        event = Event(array_device(Q))
        event = update!(array_device(Q), groupsize)(
            rv_Rstage,
            rv_Q,
            rv_Qstage,
            RKA[s, 1],
            RKA[s, 2],
            RKB[s],
            dt,
            slow_δ,
            slow_rv_dQ,
            slow_scaling;
            ndrange = length(rv_Q),
            dependencies = (event,),
        )
        wait(array_device(Q), event)
    end
    rv_Q .= rv_Qstage
end

@kernel function update!(
    dQ,
    Q,
    Qstage,
    rka1,
    rka2,
    rkb,
    dt,
    slow_δ,
    slow_dQ,
    slow_scaling,
)
    i = @index(Global, Linear)
    @inbounds begin
        if slow_δ !== nothing
            dQ[i] += slow_δ * slow_dQ[i]
        end
        Qstage[i] = rka1 * Q[i] + rka2 * Qstage[i] + dt * rkb * dQ[i]
        if slow_scaling !== nothing
            slow_dQ[i] *= slow_scaling
        end
    end
end

"""
    SSPRK22Heuns(f, Q; dt, t0 = 0)

This function returns a [`StrongStabilityPreservingRungeKutta`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the second-order, 2-stage, strong-stability-preserving, Runge--Kutta scheme
of Shu and Osher (1988) (also known as Heun's method.)
Exact choice of coefficients from wikipedia page for Heun's method :)

### References
    @article{shu1988efficient,
      title={Efficient implementation of essentially non-oscillatory shock-capturing schemes},
      author={Shu, Chi-Wang and Osher, Stanley},
      journal={Journal of computational physics},
      volume={77},
      number={2},
      pages={439--471},
      year={1988},
      publisher={Elsevier}
    }
    @article {Heun1900,
       title = {Neue Methoden zur approximativen Integration der
       Differentialgleichungen einer unabh\"{a}ngigen Ver\"{a}nderlichen}
       author = {Heun, Karl},
       journal = {Z. Math. Phys},
       volume = {45},
       pages = {23--38},
       year = {1900}
    }
"""
function SSPRK22Heuns(F, Q::AT; dt = 0, t0 = 0) where {AT <: AbstractArray}
    T = eltype(Q)
    RT = real(T)
    RKA = [RT(1) RT(0); RT(1 // 2) RT(1 // 2)]
    RKB = [RT(1), RT(1 // 2)]
    RKC = [RT(0), RT(1)]
    StrongStabilityPreservingRungeKutta(F, RKA, RKB, RKC, Q; dt = dt, t0 = t0)
end

"""
    SSPRK22Ralstons(f, Q; dt, t0 = 0)

This function returns a [`StrongStabilityPreservingRungeKutta`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the second-order, 2-stage, strong-stability-preserving, Runge--Kutta scheme
of Shu and Osher (1988) (also known as Ralstons's method.)
Exact choice of coefficients from wikipedia page for Heun's method :)

### References
    @article{shu1988efficient,
      title={Efficient implementation of essentially non-oscillatory shock-capturing schemes},
      author={Shu, Chi-Wang and Osher, Stanley},
      journal={Journal of computational physics},
      volume={77},
      number={2},
      pages={439--471},
      year={1988},
      publisher={Elsevier}
    }
    @article{ralston1962runge,
      title={Runge-Kutta methods with minimum error bounds},
      author={Ralston, Anthony},
      journal={Mathematics of computation},
      volume={16},
      number={80},
      pages={431--437},
      year={1962},
      doi={10.1090/S0025-5718-1962-0150954-0}
    }
"""
function SSPRK22Ralstons(F, Q::AT; dt = 0, t0 = 0) where {AT <: AbstractArray}
    T = eltype(Q)
    RT = real(T)
    RKA = [RT(1) RT(0); RT(5 // 8) RT(3 // 8)]
    RKB = [RT(2 // 3), RT(3 // 4)]
    RKC = [RT(0), RT(2 // 3)]
    StrongStabilityPreservingRungeKutta(F, RKA, RKB, RKC, Q; dt = dt, t0 = t0)
end

"""
    SSPRK33ShuOsher(f, Q; dt, t0 = 0)

This function returns a [`StrongStabilityPreservingRungeKutta`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the third-order, 3-stage, strong-stability-preserving, Runge--Kutta scheme
of Shu and Osher (1988)

### References
    @article{shu1988efficient,
      title={Efficient implementation of essentially non-oscillatory shock-capturing schemes},
      author={Shu, Chi-Wang and Osher, Stanley},
      journal={Journal of computational physics},
      volume={77},
      number={2},
      pages={439--471},
      year={1988},
      publisher={Elsevier}
    }
"""
function SSPRK33ShuOsher(F, Q::AT; dt = 0, t0 = 0) where {AT <: AbstractArray}
    T = eltype(Q)
    RT = real(T)
    RKA = [RT(1) RT(0); RT(3 // 4) RT(1 // 4); RT(1 // 3) RT(2 // 3)]
    RKB = [RT(1), RT(1 // 4), RT(2 // 3)]
    RKC = [RT(0), RT(1), RT(1 // 2)]
    StrongStabilityPreservingRungeKutta(F, RKA, RKB, RKC, Q; dt = dt, t0 = t0)
end

"""
    SSPRK34SpiteriRuuth(f, Q; dt, t0 = 0)

This function returns a [`StrongStabilityPreservingRungeKutta`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the third-order, 4-stage, strong-stability-preserving, Runge--Kutta scheme
of Spiteri and Ruuth (1988)

### References
    @article{spiteri2002new,
      title={A new class of optimal high-order strong-stability-preserving time discretization methods},
      author={Spiteri, Raymond J and Ruuth, Steven J},
      journal={SIAM Journal on Numerical Analysis},
      volume={40},
      number={2},
      pages={469--491},
      year={2002},
      publisher={SIAM}
    }
"""
function SSPRK34SpiteriRuuth(
    F,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    T = eltype(Q)
    RT = real(T)
    RKA = [RT(1) RT(0); RT(0) RT(1); RT(2 // 3) RT(1 // 3); RT(0) RT(1)]
    RKB = [RT(1 // 2); RT(1 // 2); RT(1 // 6); RT(1 // 2)]
    RKC = [RT(0); RT(1 // 2); RT(1); RT(1 // 2)]
    StrongStabilityPreservingRungeKutta(F, RKA, RKB, RKC, Q; dt = dt, t0 = t0)
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
include("vdw.jl")
include("rk.jl")
include("rks.jl")
include("pr.jl")
const VdW = VanDerWaals
const RK = RedlichKwong
const PR = PengRobinson
const SRK = SoaveRedlichKwong

export VdW
export RK
export SRK
export PR