text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).
struct MCMCDiagnostics
    samples::DensitySampleVector{<:AbstractVector{<:Real}}
    chainresults::Array{}
end

MCMCDiagnostics(samples::DensitySampleVector, chainresults = []) =
    MCMCDiagnostics(unshaped.(samples), chainresults)

@recipe function f(mcmc::MCMCDiagnostics;
                    vsel = collect(1 : size(mcmc.samples.v.data, 1)),
                    chains = collect(1 : size(mcmc.chainresults, 1)),
                    diagnostics = [:histogram, :kde, :trace, :acf],
                    trace = Dict(),
                    kde = Dict(),
                    acf = Dict(),
                    histogram = Dict(),
                    description = true
                    )

    seriestype = get(plotattributes, :seriestype, :histogram)

    nparams = length(vsel)
    nchains = length(chains)
    ndiagnostics = length(diagnostics)
    description ? ndescription = 1 : ndescription = 0

    size --> (ndiagnostics*350+100, nparams*nchains*200)
    layout --> Main.Plots.grid(nparams*nchains, ndiagnostics+ndescription)

    ctr = 1

	unique_chain_ids = sort(unique(mcmc.samples.info.chainid))

    for chain in chains

        #indices of samples from current chain
        r = mcmc.samples.info.chainid .== unique_chain_ids[chain]

        for p in vsel
            if description
                @series begin
                    subplot := ctr
                    label --> ""
                    annotations --> (1.5, 2.5, "Chain $chain  \n Parameter $p \n $(nsamples(mcmc.chainresults[chain])) Samples")
                    grid --> false
                    xaxis --> false
                    yaxis --> false
                    [1,1], [2, 2]
                end
                ctr +=1
            end

            # samples from current chain
            s = flatview(mcmc.samples[r].v)[p, :]

            for d in diagnostics

                # samples histogram
                if d == :histogram

                    chain_samples = mcmc.samples[r]

                    @series begin
                        subplot := ctr
                        seriestype --> get(histogram, "seriestype", :stephist)
                        intervals --> get(histogram, "intervals", default_credibilities)
                        bins --> get(histogram, "bins", 200)
                        normalize --> get(histogram, "normalize", true)
                        colors --> get(histogram, "colors", default_colors)
                        mean --> get(histogram, "mean", false)
                        std --> get(histogram, "std", false)
                        globalmode --> get(histogram, "globalmode", false)
                        marginalmode --> get(histogram, "marginalmode", false)
                        legend --> get(histogram, "legend", true)

                        label --> get(histogram, "label", "")
                        title --> get(histogram, "title", "Histogram")
                        xguide --> get(histogram, "xlabel", "\$\\theta_$(p)\$")
                        yguide --> get(histogram, "ylabel", "\$p(\\theta_$(p))\$")

                        (chain_samples, p)
                    end

                # trace plot
                elseif d == :trace
                    x = collect(1:1:length(s))

                    @series begin
                        subplot := ctr
                        seriestype --> get(trace, "seriestype", :line)
                        label --> get(trace, "label", "")
                        title --> get(trace, "title", "Trace")
                        xguide --> get(trace, "xlabel", "Iterations")
                        yguide --> get(trace, "ylabel", "")
                        seriescolor --> get(trace, "seriescolor", :dodgerblue)
                        linecolor --> get(trace, "linecolor", :dodgerblue)
                        linestyle --> get(trace, "linestyle", :solid)
                        linewidth --> get(trace, "linewidth", 1)
                        linealpha --> get(trace, "linealpha", 1)
                        seriesalpha --> get(trace, "seriesalpha", 1)
                        markershape --> get(trace, "markershape", :none)
                        markersize --> get(trace, "markersize", 1)
                        markeralpha --> get(trace, "markeralpha", 1)
                        markercolor --> get(trace, "markercolor", :dodgerblue)
                        markerstrokealpha --> get(trace, "markerstrokealpha", 1)
                        markerstrokecolor --> get(trace, "markerstrokecolor", :dodgerblue)
                        markerstrokestyle --> get(trace, "markerstrokestyle", :solid)
                        markerstrokewidth --> get(trace, "markerstrokewidth", 1)
                        legend --> get(trace, "legend", true)

                        x, s
                    end

                # kernel density estimate
                elseif d == :kde

                    npoints = get(kde, "npoints", 2048)
                    bandwidth = get(kde, "bandwidth", KernelDensity.default_bandwidth(s))
                    boundary = get(kde, "boundary", KernelDensity.kde_boundary(s, bandwidth))
                    kernel = get(kde, "kernel", Distributions.Normal)

                    k = KernelDensity.kde(
                        s,
                        bandwidth = bandwidth,
                        boundary = boundary,
                        npoints = npoints,
                        kernel = kernel
                    )

                    @series begin
                        subplot := ctr
                        seriestype --> get(kde, "seriestype", :line)
                        label --> get(kde, "label", "")
                        title --> get(kde, "title", "KDE")
                        xguide --> get(kde, "xlabel", "\$\\theta_$(p)\$")
                        yguide --> get(kde, "ylabel", "\$p(\\theta_$(p))\$")
                        seriescolor --> get(kde, "seriescolor", :dodgerblue)
                        linecolor --> get(kde, "linecolor", :dodgerblue)
                        linestyle --> get(kde, "linestyle", :solid)
                        linewidth --> get(kde, "linewidth", 1)
                        linealpha --> get(kde, "linealpha", 1)
                        seriesalpha --> get(kde, "seriesalpha", 1)
                        markershape --> get(kde, "markershape", :none)
                        markersize --> get(kde, "markersize", 1)
                        markeralpha --> get(kde, "markeralpha", 1)
                        markercolor --> get(kde, "markercolor", :dodgerblue)
                        markerstrokealpha --> get(kde, "markerstrokealpha", 1)
                        markerstrokecolor --> get(kde, "markerstrokecolor", :dodgerblue)
                        markerstrokestyle --> get(kde, "markerstrokestyle", :solid)
                        markerstrokewidth --> get(kde, "markerstrokewidth", 1)
                        legend --> get(kde, "legend", true)

                        k.x, k.density
                    end

                # autocorrelation function
                elseif d == :acf

                    lags = get(acf, "lags", :none)
                    demean = get(acf, "demean", true)
                    lags != :none ? autocorr = StatsBase.autocor(s, lags, demean=demean) : autocorr = StatsBase.autocor(s, demean=demean)

                    @series begin
                        subplot := ctr

                        seriestype --> get(acf, "seriestype", :bar)
                        label --> get(acf, "label", "")
                        title --> get(acf, "title", "ACF")
                        xguide --> get(acf, "xlabel", "Lags")
                        yguide --> get(acf, "ylabel", "")
                        seriescolor --> get(acf, "seriescolor", :dodgerblue)
                        linecolor --> get(acf, "linecolor", :dodgerblue)
                        linestyle --> get(acf, "linestyle", :solid)
                        linewidth --> get(acf, "linewidth", 0)
                        linealpha --> get(acf, "linealpha", 0)
                        seriesalpha --> get(acf, "seriesalpha", 1)
                        markershape --> get(acf, "markershape", :none)
                        markersize --> get(acf, "markersize", 1)
                        markeralpha --> get(acf, "markeralpha", 1)
                        markercolor --> get(acf, "markercolor", :dodgerblue)
                        markerstrokealpha --> get(acf, "markerstrokealpha", 1)
                        markerstrokecolor --> get(acf, "markerstrokecolor", :dodgerblue)
                        markerstrokestyle --> get(acf, "markerstrokestyle", :solid)
                        markerstrokewidth --> get(acf, "markerstrokewidth", 1)
                        legend --> get(acf, "legend", true)
                        lags != :none ? (lags, autocorr) : (0-0.5:1:length(autocorr)-0.5, autocorr)
                    end
                end
                ctr +=1
            end
        end
    end
end


text:
"""
    table_raw_data(lines, eachlines)

From a raw data txt file creates a table and group data counting pokes
"""
function table_raw_data(lines, eachlines)
    prov = convert_in_table(lines, eachlines)
    # count_pokes!(prov) ##complex calculation look up functions on separate file
    prov[!,:PokeIn] = ismatch.(r"poke_\d$", prov.Event)
    prov[!,:PokeIn_count] = WaltonForaging.jumpcount(prov.PokeIn)
    WaltonForaging.findpokesout!(prov)
    maximum(prov.PokeIn_count) == maximum(prov.PokeOut_count)
    prov[!,:Poke] = Vector{Union{Missing,Int}}(missing,nrow(prov))
    clean_iterator = filter(x -> !ismissing(x) && x!=0 ,prov.PokeIn_count)
    for x in clean_iterator#1:maximum(prov.PokeIn_count)
        prov[findfirst(prov.PokeIn_count .== x) : findfirst(prov.PokeOut_count .== x),
            :Poke] .= x
    end
    prov[!,:Rsync_count] = accumulate(+,[ismatch(r"rsync$",ev) for ev in prov.Event], init = 0)
    return prov
end
"""
    convert_in_table(session_path::String)

From a raw data txt file creates a table that is filled in 3 ways
1 - Read first 5 lines containing session's global information
2 - Identify Plines to globally update subsequent T, AFT, IFT values in separate columns
3 - reads line by line to identify: Time, EventName and EventType in separate columns
"""

function convert_in_table(lines, eachlines)
    prov = DataFrame(Time=Float64[],Event = String[],
        Type = String[], P = Float64[],  T = Float64[],
        AFT = Float64[], IFT = Float64[])
    P, T, AFT, IFT = WaltonForaging.initialvalues(lines)
    for l in eachlines
        isempty(l) && continue
        if WaltonForaging.ismatch.(r"AFT:",l) #if the line contain AFT read it using special parsing
            P, T, AFT, IFT = WaltonForaging.parsePline(l)
        elseif l[1] == 'D'
            r = split(l, " ")
            time = parse(Float64,r[2])
            code = r[3]
            push!(prov,(Time = time,
                Event = get(WaltonForaging.InfoDict,code,"nope: $code"),
                Type = parse(Float64,code) <=16 ? "State" : "Action",
                P = P, T = T, AFT = AFT, IFT = IFT))
        end
    end
    return prov
end

"""
    initialvalues(session_path::String)

Opens a raw data file as a vector of string and identifies the first P-line
containing T, AFT and IFT info to initiate `table_raw_data(session_path)`
"""
function initialvalues(lines)
    first_Pline = lines[findfirst(ismatch.(r"AFT:",lines))]
    P, T, AFT, IFT = parsePline(first_Pline)
    return P, T, AFT, IFT
end

"""
    parsePline(line::String)

When feeded a P line in the format
P 36861 R#:5 P#:1 RP#:5 P:-1 T:3694 AFT:3445 IFT:227.5831
    (numbers are just an example)
split the string using space character in a vector
takes the last 3 elements [T:3694,AFT:3445,IFT:227.5831]
Iterates the elements split them using column character
Parse the second element after column in a Float values
returns the 3 values in a tuple (3694,3445,227.5831)
Therefore results have to be assign in the returning variables in the
correct order: T, AFT, IF
"""
function parsePline(line::String)
    pre = split(line," ")
    starting_vals = vcat(pre[4],pre[end-2:end])
    tuple([parse(Float64,split(v,":")[2]) for v in starting_vals]...)
end

ext: jl


text:
for N in [Float64, Rational{Int}, Float32]
    # random hyperplane
    rand(Hyperplane)

    # normal constructor
    a = ones(N, 3)
    b = N(5)
    hp = Hyperplane(a, b)

    # corner case: zero normal vector
    @test_throws AssertionError Hyperplane(N[0, 0], N(1))

    # dimension
    @test dim(hp) == 3

    # support function
    hp2 = Hyperplane(N[1, 0], N(1))
    @test ρ(N[2, 0], hp2) == N(2)
    @test ρ(N[-2, 0], hp2) == N(-2)
    @test ρ(N[1, 1], hp2) == N(Inf)

    # support vector and membership
    function test_svec(hp)
        d1 = copy(hp.a)
        @test σ(d1, hp) ∈ hp
        @test σ(N(2) * d1, hp) ∈ hp
        d2 = N[1, 0, 0]
        @test_throws ErrorException σ(d2, hp)
        d3 = N[1, 1, 2]
        @test_throws ErrorException σ(d3, hp)
        d4 = zeros(N, 3)
        @test σ(d4, hp) ∈ hp
    end
    # tests 1
    test_svec(hp)
    # tests 2
    test_svec(Hyperplane(N[0, 0, 1], b))

    # support vector in opposite direction
    hp2 = Hyperplane(N[1], N(1))
    @test σ(N[-1], hp2) ∈ hp2
    # support vector in other directions throws an error (but see #750)
    @test_throws ErrorException σ(N[1, 1], Hyperplane(N[1, 0], N(1)))

    # boundedness
    @test !isbounded(hp)

    # universality
    @test !isuniversal(hp)
    res, w = isuniversal(hp, true)
    @test !res && w ∉ hp

    # isempty
    @test !isempty(hp)

    # an_element and membership
    @test an_element(hp) ∈ hp

    # constrained dimensions
    @test constrained_dimensions(Hyperplane(N[1, 0, 1], N(1))) == [1, 3]
    @test constrained_dimensions(Hyperplane(N[0, 1, 0], N(1))) == [2]

    # constraints_list
    @test ispermutation(constraints_list(hp),
                        [HalfSpace(a, b), HalfSpace(-a, -b)])

    # test concrete linear map of a hyerplane
    H = Hyperplane(N[1, -1], N(0)) # x = y
    M = N[1 0; 0 0] # non-invertible matrix
    @test_throws ArgumentError linear_map(M, H)
    M = N[2 2; 0 1] # invertible matrix
    @test linear_map(M, H) == Hyperplane(N[0.5, -2.0], N(0.0))

    # translation
    @test translate(hp, N[1, 2, 3]) == Hyperplane(ones(N, 3), N(11))

    # intersection emptiness
    u = Universe{N}(3)
    res, v = is_intersection_empty(u, hp, true)
    @test !is_intersection_empty(u, hp) && !res && v ∈ hp && v ∈ u
    res, v = is_intersection_empty(hp, u, true)
    @test !is_intersection_empty(hp, u) && !res && v ∈ hp && v ∈ u
end

# Polyhedra tests that only work with Float64
for N in [Float64]
    hp = Hyperplane(ones(N, 3), N(5))

    # intersection emptiness
    b = BallInf(zeros(N, 3), N(1))
    empty_intersection, v = is_intersection_empty(b, hp, true)
    @test empty_intersection && is_intersection_empty(b, hp) && v == N[]
    b = BallInf(N[3, 3, 3], N(1))
    empty_intersection, v = is_intersection_empty(b, hp, true)
    @test empty_intersection && is_intersection_empty(b, hp) && v == N[]
    b = BallInf(N[2, 2, 2], N(1))
    empty_intersection, v = is_intersection_empty(b, hp, true)
    @test !empty_intersection && !is_intersection_empty(b, hp) && v ∈ hp && v ∈ b
end


text:
################################################################################
#
#    NfOrd/NfOrdFracIdl.jl : Fractional ideals of generic
#                               orders in number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export frac_ideal

################################################################################
#
#  Construction
#
################################################################################

doc"""
***
    frac_ideal(O::NfOrd, A::FakeFmpqMat) -> NfOrdFracIdl

> Creates the fractional ideal of $\mathcal O$ with basis matrix $A$.
"""
function frac_ideal(O::NfOrd, x::FakeFmpqMat)
  y = hnf(x)
  z = NfOrdFracIdl(O, y)
  return z
end

doc"""
***
    frac_ideal(O::NfOrd, A::fmpz_mat, b::fmpz) -> NfOrdFracIdl

> Creates the fractional ideal of $\mathcal O$ with basis matrix $A/b$.
"""
function frac_ideal(O::NfOrd, x::fmpz_mat, y::fmpz)
  y = FakeFmpqMat(x, y)
  z = NfOrdFracIdl(O, y)
  return z
end

frac_ideal(O::NfOrd, x::fmpz_mat, y::Integer) = frac_ideal(O, x, fmpz(y))

doc"""
***
    frac_ideal(O::NfOrd, I::NfOrdIdl) -> NfOrdFracIdl

> Turns the ideal $I$ into a fractional ideal of $\mathcal O$.
"""
function frac_ideal(O::NfOrd, x::NfOrdIdl)
  z = NfOrdFracIdl(O, x, fmpz(1))
  return z
end

doc"""
***
    frac_ideal(O::NfOrd, I::NfOrdIdl, b::fmpz) -> NfOrdFracIdl

> Creates the fractional ideal $I/b$ of $\mathcal O$.
"""
function frac_ideal(O::NfOrd, x::NfOrdIdl, y::fmpz)
  z = NfOrdFracIdl(O, x, deepcopy(y)) # deepcopy x?
  return z
end

frac_ideal(O::NfOrd, x::NfOrdIdl, y::Integer) = frac_ideal(O, x, fmpz(y))

doc"""
***
    frac_ideal(O::NfOrd, a::nf_elem) -> NfOrdFracIdl

> Creates the principal fractional ideal $(a)$ of $\mathcal O$.
"""
function frac_ideal(O::NfOrd, x::nf_elem)
  z = NfOrdFracIdl(O, deepcopy(x))
  return z
end

doc"""
***
    frac_ideal(O::NfOrd, a::NfOrdElem) -> NfOrdFracIdl

> Creates the principal fractional ideal $(a)$ of $\mathcal O$.
"""
function frac_ideal(O::NfOrd, x::NfOrdElem)
  z = NfOrdFracIdl(O, elem_in_nf(x))
  return z
end

################################################################################
#
#  String I/O
#
################################################################################

export basis_mat, norm, inv, ==, *, integral_split

export parent, order, basis_mat, basis_mat_inv, basis, norm,
       ring_of_multipliers, ==

################################################################################
#
#  Parent
#
################################################################################

parent(a::NfOrdFracIdl) = NfOrdFracIdlSet(order(a), false)

function Base.hash(a::NfOrdFracIdl, h::UInt)
  b = simplify(a)
  return hash(b.num, hash(b.den, h))
end

################################################################################
#
#  Order
#
################################################################################

order(a::NfOrdFracIdlSet) = a.order

order(a::NfOrdFracIdl) = a.order

################################################################################
#
#  (Inverse) basis matrix
#
################################################################################

doc"""
***
    basis_mat_inv(I::NfOrdFracIdl) -> FakeFmpqMat

> Returns the inverse of the basis matrix of $I$.
"""
function basis_mat_inv(a::NfOrdFracIdl)
  if isdefined(a, :basis_mat_inv)
    return deepcopy(a.basis_mat_inv)
  else
    a.basis_mat_inv = inv(basis_mat(a))
    return deepcopy(a.basis_mat_inv)
  end
end

################################################################################
#
#  Basis
#
################################################################################

doc"""
***
    basis(I::NfOrdFracIdl) -> Array{nf_elem, 1}

> Returns the $\mathbf Z$-basis of $I$.
"""
function basis(a::NfOrdFracIdl)
  B = basis_mat(a)
  d = degree(order(a))
  O = order(a)
  K = nf(O)
  Oba = O.basis_nf
  res = Array{nf_elem}(d)

  for i in 1:d
    z = K()
    for j in 1:d
      z = z + B.num[i, j]*Oba[j]
    end
    z = divexact(z, B.den)
    res[i] = z
  end

  return res
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(x::NfOrdFracIdl, dict::ObjectIdDict)
  z = NfOrdFracIdl(deepcopy(x.num), deepcopy(x.den))
  if isdefined(x, :basis_mat)
    z.basis_mat = deepcopy(x.basis_mat)
  end
  return z
end

################################################################################
#
#  Basis matrix
#
################################################################################

doc"""
***
    basis_mat(I::NfOrdFracIdl) -> FakeFmpqMat

> Returns the basis matrix of $I$ with respect to the basis of the order.
"""
function basis_mat(x::NfOrdFracIdl)
  if isdefined(x, :basis_mat)
    return deepcopy(x.basis_mat)
  else
    x.basis_mat = FakeFmpqMat(basis_mat(numerator(x)), denominator(x))
    return deepcopy(x.basis_mat)
  end
end

################################################################################
#
#  Numerator and denominator
#
################################################################################

numerator(x::NfOrdFracIdl) = x.num

denominator(x::NfOrdFracIdl) = deepcopy(x.den)

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, s::NfOrdFracIdlSet)
   print(io, "Set of fractional ideals of ")
   print(io, s.order)
end

function show(io::IO, id::NfOrdFracIdl)
  if isdefined(id, :num) && isdefined(id, :den)
    print(io, "1//", id.den, " * ")
    print(io, id.num)
  else
    print(io, "Fractional ideal with basis matrix\n")
    print(io, id.basis_mat)
  end
end

################################################################################
#
#  Norm
#
################################################################################

doc"""
***
    norm(I::NfOrdFracIdl) -> fmpq

> Returns the norm of $I$
"""
function norm(A::NfOrdFracIdl)
  if isdefined(A, :norm)
    return deepcopy(A.norm)
  else
    A.norm = norm(A.num)//A.den^degree(order(A))
    return deepcopy(A.norm)
  end
end

################################################################################
#
#  Minimum
#
################################################################################

function minimum(A::NfOrdFracIdl)
  return minimum(A.num)//A.den
end

################################################################################
#
#  Inverse
#
################################################################################

doc"""
***
    inv(A::NfOrdFracIdl) -> NfOrdFracIdl

> Returns the fractional ideal $B$ such that $AB = \mathcal O$.
"""
function inv(A::NfOrdFracIdl)
  B = inv(A.num)
  g = gcd(B.den, A.den)
  B.den = divexact(B.den, g)
  a = divexact(A.den, g)
  return prod_by_int(B, a)
end

################################################################################
#
#  Simplification
#
################################################################################

function simplify_exact!(A::NfOrdFracIdl)
  g = A.den

  A.den = fmpz(1)
  A.num = divexact(A.num, g)
end

function simplify(A::NfOrdFracIdl)
  simplify(A.num)

  if has_2_elem(A.num)
    ZK = order(A)
    g = Base.reduce(gcd, elem_in_basis(ZK(A.num.gen_two)))
    g = gcd(g, A.den)
    g = gcd(g, A.num.gen_one)
  else
    b = basis_mat(A.num, Val{false})
    g = gcd(denominator(A), content(b))
  end

  if g != 1
    A.num = divexact(A.num, g)
    A.den = divexact(A.den, g)
  end

  return A
end

################################################################################
#
#  Primeness
#
################################################################################

isprime_known(A::NfOrdFracIdl) = isprime_known(A.num)

################################################################################
#
#  Equality
#
################################################################################

doc"""
***
    ==(x::NfOrdFracIdl, y::NfOrdFracIdl) -> Bool

> Returns whether $x$ and $y$ are equal.
"""
function ==(A::NfOrdFracIdl, B::NfOrdFracIdl)
  #return B.den * basis_mat(A.num) == A.den * basis_mat(B.num)
  if (!ismaximal_known(order(A)) || (ismaximal_known(order(A)) && !ismaximal(order(A))))
    return basis_mat(A) == basis_mat(B)
  end

  if !isdefined(A, :num) || !isdefined(B, :num)
    return basis_mat(A) == basis_mat(B)
  end
  D = inv(B)
  E = prod(A, D)
  C = simplify(E)
  return isone(C.den) && isone(norm(C))
end

################################################################################
#
#  Binary operations
#
################################################################################

function prod(a::NfOrdFracIdl, b::NfOrdFracIdl)
  A = a.num* b.num
  return NfOrdFracIdl(A, a.den*b.den)
end

doc"""
***
    *(I::NfOrdFracIdl, J::NfOrdFracIdl) -> NfOrdFracIdl

> Returns $IJ$.
"""
*(A::NfOrdFracIdl, B::NfOrdFracIdl) = prod(A, B)

################################################################################
#
#  Powering
#
################################################################################

function ^(A::NfOrdFracIdl, a::Int)
  if a == 0
    B = NfOrdFracIdl(ideal(order(A), 1), fmpz(1))
    return B
  end

  if a == 1
    return A # copy?
  end

  if a < 0
    return inv(A^(-a))
  end

  if a == 2
    return A*A
  end

  if mod(a, 2) == 0
    return (A^div(a, 2))^2
  else
    return A * A^(a - 1)
  end
end

function //(A::NfOrdFracIdl, B::NfOrdFracIdl)
  C = prod(A, inv(B))
  return C
end

################################################################################
#
#  Ad hoc binary operations
#
################################################################################

function prod_by_int(A::NfOrdFracIdl, a::fmpz)
  return NfOrdFracIdl(A.num * a, A.den)
end

*(A::NfOrdFracIdl, a::fmpz) = prod_by_int(A, a)
*(a::fmpz, A::NfOrdFracIdl) = prod_by_int(A, a)

function *(A::NfOrdIdl, B::NfOrdFracIdl)
  z = NfOrdFracIdl(A*B.num, B.den)
  return z
end

*(A::NfOrdFracIdl, B::NfOrdIdl) = NfOrdFracIdl(A.num*B, A.den)

function *(A::NfOrdFracIdl, a::nf_elem)
  C = *(A, a*order(A))
  return C
end

*(a::nf_elem, A::NfOrdFracIdl) = A*a

function *(A::NfOrdIdl, a::nf_elem)
  C = *(A, a*order(A))
  return C
end

*(a::nf_elem, A::NfOrdIdl) = A*a

function //(A::NfOrdFracIdl, B::NfOrdIdl)
  C = prod(A, inv(B))
  return C
end

function //(A::NfOrdIdl, B::NfOrdIdl)
  return A*inv(B)
end

function //(A::NfOrdIdl, B::NfOrdFracIdl)
  return A*inv(B)
end

function //(A::NfOrdFracIdl, a::nf_elem)
  C = prod(A, Idl((order(A), inv(a))))
  return C
end

function //(A::NfOrdIdl, d::fmpz)
  return Hecke.NfOrdFracIdl(A, d)
end

function //(A::NfOrdIdl, d::Integer)
  return A//fmpz(d)
end

function +(A::NfOrdIdl, B::NfOrdFracIdl)
  return (A*denominator(B)+numerator(B))//denominator(B)
end

+(A::NfOrdFracIdl, B::NfOrdIdl) = B+A

function +(A::NfOrdFracIdl, B::Hecke.NfOrdFracIdl)
  d = lcm(denominator(A), denominator(B))
  ma = div(d, denominator(A))
  mb = div(d, denominator(B))
  return (numerator(A)*ma + numerator(B)*mb)//d
end

function *(x::nf_elem, y::NfOrd)
  b, z = _check_elem_in_order(denominator(x, y)*x, y)
  return NfOrdFracIdl(ideal(y, y(z)), denominator(x, y))
end
################################################################################
#
#  Conversion
#
################################################################################

function (ord::NfOrdIdlSet)(b::NfOrdFracIdl)
   b.den > 1 && error("not integral")
   return b.num
end

function ideal(O::NfOrd, a::nf_elem)
  return a*O
end

doc"""
***
    integral_split(A::NfOrdFracIdl) -> NfOrdIdl, NfOrdIdl
> Computes the unique coprime integral ideals $N$ and $D$ s.th. $A = ND^{-1}$
"""
function integral_split(A::NfOrdFracIdl)
  d = simplify(inv(A + ideal(order(A), fmpz(1))))
  @assert denominator(d) == 1
  n = simplify(A*d)
  @assert denominator(n) == 1
  return numerator(n), numerator(d)
end

doc"""
    factor(I::NfOrdFracIdl) -> Dict{NfOrdIdl, Int}
> The factorisation of $I$.
"""
function factor(I::NfOrdFracIdl)
  if iszero(norm(I))
    error("Cannot factor zero ideal")
  end
  n, d = integral_split(I)
  fn = factor(n)
  fd = factor(d)
  for (k, v) = fd
    if haskey(fn, k)
      fn[k] -= v
    else
      fn[k] = -v
    end
  end
  return fn
end

function one(A::NfOrdFracIdl)
  return NfOrdFracIdl(ideal(order(A), 1), fmpz(1))
end

doc"""
    valuation(A::NfOrdFracIdl, p::NfOrdIdl)
> The valuation of $A$ at $p$.
"""
function valuation(A::NfOrdFracIdl, p::NfOrdIdl)
  return valuation(A.num, p) - valuation(A.den, p)
end


text:
@testset "De Bruijn Neighbors" begin
    @test collect(neighbors(DNAKmer("ACG")))  ==
        map(DNAKmer, ["CGA",  "CGC",  "CGG",  "CGT"])
    @test collect(neighbors(DNAKmer("GGGG"))) ==
        map(DNAKmer, ["GGGA", "GGGC", "GGGG", "GGGT"])
    @test collect(neighbors(RNAKmer("ACG")))  ==
        map(RNAKmer, ["CGA",  "CGC",  "CGG",  "CGU"])
    @test collect(neighbors(RNAKmer("GGGG"))) ==
        map(RNAKmer, ["GGGA", "GGGC", "GGGG", "GGGU"])
end


text:
"""A node which just wraps a block."""
struct BlockNode{T} <: NodeOp{T}
    block::Block{T}
end

function Base.show(io::IO, ::BlockNode{T}) where {T}
    return print(io, "BlockNode{$T}")
end

Base.hash(x::BlockNode, h::UInt64) = hash(x.block, hash(:BlockNode, h))
Base.:(==)(x::BlockNode{T}, y::BlockNode{T}) where {T} = x.block == y.block
Base.isequal(x::BlockNode{T}, y::BlockNode{T}) where {T} = isequal(x.block, y.block)

stateless(::BlockNode) = true

function run_node!(
    op::BlockNode{T}, ::EmptyNodeEvaluationState, time_start::DateTime, time_end::DateTime
) where {T}
    return _slice(op.block, time_start, time_end)
end

"""
    block_node(block::Block)

Construct a node whose values are read directly from the given `block`.
"""
block_node(block::Block) = obtain_node((), BlockNode(block))

# TODO Identity mapping... probably just want a cache of empty blocks by T somewhere?
"""
    empty_node(T)

Construct a node with value type `T` which, if evaluated, will never tick.
"""
empty_node(T) = block_node(Block{T}())

"""
A node op which ticks once a day at the specified local time of day.
"""
struct Iterdates <: NodeOp{DateTime}
    time_of_day::Time
    tz::TimeZone
end

stateless(::Iterdates) = true

"""
    _to_utc(dt::DateTime, tz::TimeZone) -> DateTime

Interpret naive `dt` to be in timezone `tz`. Return a naive `DateTime`, but in UTC.
"""
function _to_utc(dt::DateTime, tz::TimeZone)
    # If we are already in UTC, there's nothing to do.
    tz == tz"UTC" && return dt
    return ZonedDateTime(dt, tz).utc_datetime
end

function run_node!(
    op::Iterdates, ::EmptyNodeEvaluationState, time_start::DateTime, time_end::DateTime
)
    times = DateTime[]

    # Add some padding for safety — we will only emit those times that are in the correct
    # time range.
    d1 = Date(time_start) - Day(2)
    d2 = Date(time_end) + Day(2)
    for d in d1:Day(1):d2
        # Get the appropriate time of day in this timezone.
        t = _to_utc(d + op.time_of_day, op.tz)
        t < time_start && continue
        t >= time_end && break
        push!(times, t)
    end

    return Block(unchecked, times, times)
end

# TODO add the equivalent of an RDate step & offset, for valid steps. Not sure what the
#   equivalent of these in Julia would be...
"""
    iterdates(time_of_day::Time=Time(0), tz::TimeZone=tz"UTC", occurrence=1)

Create a node which ticks exactly once a day at `time_of_day` in timezone `tz`.

This defaults to midnight in UTC. If `tz` is set otherwise, then each knot will appear at
`time_of_day` in that timezone.

Note that:
    * All knot times in `TimeDag` are considered to be in UTC.
    * It is possible to select a `time_of_day` that does not exist for every day. This will
        lead to an exception being raised during evaluation.

In a given knot, each value will be of type `DateTime`, and equal the time of the knot.
"""
function iterdates(time_of_day::Time=Time(0), tz::TimeZone=tz"UTC")
    return obtain_node((), Iterdates(time_of_day, tz))
end

"""
A node op which ticks every `delta`, such that a knot would appear on `epoch`.
"""
struct Pulse <: NodeOp{DateTime}
    delta::Millisecond
    epoch::DateTime
end

stateless(::Pulse) = true

function run_node!(
    op::Pulse, ::EmptyNodeEvaluationState, time_start::DateTime, time_end::DateTime
)
    rem_ = (time_start - op.epoch) % op.delta
    pulse_start = rem_ == Millisecond(0) ? time_start : time_start + (op.delta - rem_)

    rem_ = (time_end - op.epoch) % op.delta
    pulse_end = rem_ == Millisecond(0) ? time_end - op.delta : time_end - rem_

    times = collect(pulse_start:(op.delta):pulse_end)
    return Block(unchecked, times, times)
end

# The Julia epoch. This is the DateTime whose internal representation is "0".
const _JULIA_EPOCH = DateTime(0, 1, 1) - Millisecond(Dates.DATETIMEEPOCH)

"""
    pulse(delta::TimePeriod[; epoch::DateTime])

Obtain a node which ticks every `delta`. Each value will equal the time of the knot.

Knots will be placed such that the difference between its time and `epoch` will always be an
integer multiple of `delta`. By default `epoch` is set to the Julia `DateTime` epoch, which
is `DateTime(0, 12, 31)`.
"""
function pulse(delta::TimePeriod; epoch::DateTime=_JULIA_EPOCH)
    delta > Millisecond(0) || throw(ArgumentError("delta must be positive, got $delta"))
    return obtain_node((), Pulse(Millisecond(delta), epoch))
end

# TODO We may want to generalise or otherwise refactor to allow reading multiple value
# fields.
# TODO we need some nicer APIs here, which will necessitate a new TeaFiles version.
"""A node which will read from an appropriately formatted TeaFile."""
struct TeaFileNode{T} <: NodeOp{T}
    path::String
    i_time_field::Int64
    i_value_field::Int64

    function TeaFileNode(path::AbstractString, value_field_name::AbstractString)
        # Verify that the file has a valid time index.
        metadata = open(path) do f
            read(f, TeaFiles.Header.TeaFileMetadata)
        end
        time_section = TeaFiles.Header.get_section(metadata, TeaFiles.Header.TimeSection)
        time_field = TeaFiles.Header.get_primary_time_field(metadata)
        compatible_time = (
            TeaFiles.Header.field_type(time_field) == Int64 &&
            TeaFiles.Header.is_julia_time_compatible(time_section)
        )
        if !compatible_time
            throw(ArgumentError("Got incompatible time field."))
        end

        item_section = TeaFiles.Header.get_section(metadata, TeaFiles.Header.ItemSection)
        i_value_field = only(findall(x -> x.name == value_field_name, item_section.fields))
        T = TeaFiles.Header.field_type(item_section.fields[i_value_field])

        # This is guaranteed to be unique.
        i_time_field = findfirst(x -> x.offset == time_field.offset, item_section.fields)

        return new{T}(path, i_time_field, i_value_field)
    end
end

function Base.hash(x::TeaFileNode, h::UInt64)
    return hash(x.path, hash(x.i_time_field, hash(x.i_value_field, hash(:TeaFileNode, h))))
end
function Base.:(==)(x::TeaFileNode{T}, y::TeaFileNode{T}) where {T}
    return (
        x.path == y.path &&
        x.i_time_field == y.i_time_field &&
        x.i_value_field == y.i_value_field
    )
end

# TODO One could imagine a more optimal reading process where the state contains the
#   location in the file we're up to, to avoid a binary search.
#   The current reading process in `run_node!` does wasteful allocation, which could also
#   be improved by dealing with a lower level API in TeaFiles (when exposed).
stateless(::TeaFileNode) = true

function run_node!(
    op::TeaFileNode{T}, ::EmptyNodeEvaluationState, time_start::DateTime, time_end::DateTime
) where {T}
    rows = TeaFiles.read(op.path; lower=time_start, upper=time_end)
    times = Vector{DateTime}(undef, length(rows))
    values = Vector{T}(undef, length(rows))

    @inbounds for (i, row) in enumerate(rows)
        times[i] = getfield(row, op.i_time_field)
        values[i] = getfield(row, op.i_value_field)
    end

    # Note that this is *not* unchecked, since we should not trust data that we retrieve
    # from disk - e.g. it may have been corrupted.
    return Block(times, values)
end

"""
    tea_file(path::AbstractString, value_field_name)

Get a node that will read data from the tea file at `path`.

Such a tea file must observe the following properties, which will be verified at runtime:
* Have a primary time field which is compatible with a Julia `DateTime`.
* Have exactly one column with name `value_field_name`.
* Have *strictly* increasing times.

Upon node creation, the metadata section of the file will be parsed to infer the value type
of the resulting node. However, the bulk of the data will only be read at evaluation time.

# See also
* The [tea file spec](http://discretelogics.com/resources/teafilespec/)
* [TeaFiles.jl](https://github.com/tpgillam/TeaFiles.jl)
"""
function tea_file(path::AbstractString, value_field_name::AbstractString)
    return obtain_node((), TeaFileNode(path, value_field_name))
end
tea_file(path::AbstractString, field_name::Symbol) = tea_file(path, string(field_name))


text:
# C Language Interface.

"""
    @gc_safe @ccall f(arg::Type, ...)::Type

Wrap a `@ccall` with `jl_gc_safe_[enter|leave]`.

This allows the GC to run during the `@ccall`.
Without this all threads will be blocked if the GC tries to run
during the `@ccall`.

The `@ccall` must not call back into Julia.
"""
macro gc_safe(ex)
    quote
        old_state = @ccall jl_gc_safe_enter()::Int8
        n = $(esc(ex))
        @ccall jl_gc_safe_leave(old_state::Int8)::Cvoid
        n
    end
end

# End of file: ccall.jl


text:
using DataFrames
using NetCDF
using IcanProj
using ProgressMeter

function aggregate_forcings(path, var, res)

    # Symbol for resolution

    ind_res = Symbol("ind_$(res)")

    # Metadata from table

    df_meta = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "df_links.csv")) |> DataFrame

    # Metadata from input netcdfs

    file_in = joinpath(path, "forcings_ts/$(var)_1km.nc")

    nc_in = Dict("ind_senorge" => ncread(file_in, "id"),
                 "time_str" => ncread(file_in, "time_str"),
                 "var_atts" => ncgetatt(file_in, var, "units"))

    # Create netcdfs

    file_out = joinpath(path, "forcings_ts/$(var)_$(res).nc")

    var_atts = Dict("units" => nc_in["var_atts"])

    dim_time = length(nc_in["time_str"])

    dim_space = length(unique(df_meta[ind_res]))

    time_str = nc_in["time_str"]

    id = convert(Array{Int64}, unique(df_meta[ind_res]))

    id_desc = String(ind_res)

    create_netcdf(file_out, var, var_atts, dim_time, dim_space, time_str, id, id_desc)

    # Aggregate forcings

    icol_out = 1

    @showprogress 1 "Running..." for ind_unique in unique(df_meta[ind_res])

        itarget = find(df_meta[ind_res] .== ind_unique)

        data_agg = fill(0.0, dim_time)

        tmp = fill(0.0, dim_time)

        for ind_senorge in df_meta[:ind_senorge][itarget]

            icol_in = find(nc_in["ind_senorge"] .== ind_senorge)[1]

            ncread!(file_in, var, tmp, start=[1,icol_in], count=[-1,1])

            data_agg = data_agg + tmp / length(itarget)

        end

        ncwrite(data_agg, file_out, var, start=[1,icol_out], count=[-1,1])

        icol_out += 1

    end

    ncclose()

    return nothing

end

# Settings

path = "/data04/jmg/fsm_simulations/netcdf"

var_all = ["ilwr", "iswr", "pres", "rainf", "snowf", "wind", "rhum", "tair"]

res_all = ["50km", "25km", "10km", "5km"]

# Loop over all variables and resolutions

for var in var_all, res in res_all

    aggregate_forcings(path, var, res)

    print("Finished $(var) for $(res)\n")

end


text:
include_dependency("guess_task_layout.jl")
const TASK_TID_OFFSET = let script = joinpath(@__DIR__, "guess_task_layout.jl")
    out = read(`$(Base.julia_cmd()) --startup-file=no $script`, String)
    parse(Int, out)
end

include("taskhack_impl.jl")
check_set_raw_tid()


text:
#!/usr/bin/env julia
# extractvars.jl is a script to print out the variables contained in a julia program
# Example Usage:
#
#      julia ../bin/extractvars.jl epicookbook/src/SIRModel.jl epicookbook/src/ScalingModel.jl
#
# Example Output:
# loaded
# epicookbook/src/SIRModel.jl, "sir_sol"
# epicookbook/src/SIRModel.jl, "tspan"
# epicookbook/src/SIRModel.jl, "du[2]"
# epicookbook/src/SIRModel.jl, "(b, g)"
# epicookbook/src/SIRModel.jl, "S"
# epicookbook/src/SIRModel.jl, "pram"
# epicookbook/src/SIRModel.jl, "init"
# epicookbook/src/SIRModel.jl, "sir_prob2"
# epicookbook/src/SIRModel.jl, "(S, I, R)"
# epicookbook/src/SIRModel.jl, "du[3]"
# epicookbook/src/SIRModel.jl, "parms"
# epicookbook/src/SIRModel.jl, "β"
# epicookbook/src/SIRModel.jl, "γ"
# epicookbook/src/SIRModel.jl, "I"
# epicookbook/src/SIRModel.jl, "du[1]"
# epicookbook/src/SIRModel.jl, "sir_prob"
# epicookbook/src/ScalingModel.jl, "α"
# epicookbook/src/ScalingModel.jl, "sir_sol"
# epicookbook/src/ScalingModel.jl, "tspan"
# epicookbook/src/ScalingModel.jl, "βs[:, i]"
# epicookbook/src/ScalingModel.jl, "w"
# epicookbook/src/ScalingModel.jl, "μ"
# epicookbook/src/ScalingModel.jl, "K"
# epicookbook/src/ScalingModel.jl, "init"
# epicookbook/src/ScalingModel.jl, "(β, r, μ, K, α)"
# epicookbook/src/ScalingModel.jl, "dS"
# epicookbook/src/ScalingModel.jl, "ws"
# epicookbook/src/ScalingModel.jl, "(S, I)"
# epicookbook/src/ScalingModel.jl, "parms"
# epicookbook/src/ScalingModel.jl, "du"
# epicookbook/src/ScalingModel.jl, "r"
# epicookbook/src/ScalingModel.jl, "dI"
# epicookbook/src/ScalingModel.jl, "β"
# epicookbook/src/ScalingModel.jl, "m"
# epicookbook/src/ScalingModel.jl, "βs"
# epicookbook/src/ScalingModel.jl, "sir_prob"
# epicookbook/src/ScalingModel.jl, "i"
try
    using SemanticModels
catch

    import Pkg;
    Pkg.add("SemanticModels")
end

using SemanticModels.Parsers
import SemanticModels.Parsers.walk

"""    findassigns(expr::Expr)

findassign walks the AST of `expr` to find the assignments to variables.

This function returns a reference to the original expression so that you can modify it inplace
and is intended to help users rewrite expressions for generating new models.

See also: [`findassign`](@ref).
"""
function findassigns(expr::Expr)
    # g(y) = filter(x->x!=nothing, y)
    matches = Dict{Any,Set{Expr}}()
    g(y::Any) = :()
    f(x::Any) = :()
    f(x::Expr) = begin
        if x.head == :(=)
            s = get(matches, x.args[1], Set{Expr}())
            matches[x.args[1]] = push!(s, x)
            return x
        end
        walk(x, f, g)
    end
    walk(expr, f, g)
    return matches
end

"""    findvars(file)

finds all the variable names from a julia script.

Example:
    file = "epicookbook/src/DiscreteTimeSIR.jl"
    matches = findvars(file)
"""
function findvars(file)
    ex = parsefile(file)
    matches = findassigns(ex)
    matches = matches
    return matches
end

if length(ARGS) > 0
    for file in ARGS
        try
            matches = findvars(file)
            for k in keys(matches)
                println("$file,\"$k\"")
            end
        catch
            @warn "Processing failed" file=file
        end
    end
end


text:
module QuantileRegressions

    import StatsModels: TableRegressionModel, FormulaTerm, coef, @formula
    export qreg, coef, vcov, stderr, quantiles, IP, IRLS, @formula

    using DataFrames, Distributions, LinearAlgebra, LinearAlgebra.BLAS, StatsModels, StatsBase

    mutable struct QRegModel
        beta::Vector{Float64}
        vcov::Matrix{Float64}
        stderr::Vector{Float64}
        q
    end

    abstract type Solver end

    struct IP <: Solver
    end

    struct IRLS{T} <: Solver
        tol::T
        threshold::T
        maxiter::Integer
    end
    IRLS(;tol::Real = 1e-12, maxiter::Integer = 1_000, threshold::Real = 1e-5) = IRLS(tol, threshold, maxiter)

    include("InteriorPoint.jl")
    include("IRLS.jl")
    include("Covariance.jl")
    include("npqreg.jl")

    function qreg(f::FormulaTerm, df::AbstractDataFrame, q, weights::AbstractVector, s::Solver = IP())
        mf = ModelFrame(f, df)
        mm = ModelMatrix(mf)
        mm = ModelMatrix(mm.m.*weights, mm.assign)
        mr = response(mf).*weights

        coef = qreg_coef(mr, mm.m, q, s)
        vcov = qreg_vcov(mr, mm.m, coef, q)
        stderror = sqrt.(diag(vcov))
        return TableRegressionModel(QRegModel(coef, vcov, stderror, q), mf, mm)
    end
    function qreg(f::FormulaTerm, df::AbstractDataFrame, q, s::Solver = IP())
        mf = ModelFrame(f, df)
        mm = ModelMatrix(mf)
        mr = response(mf)
        coef = qreg_coef(mr, mm.m, q, s)
        vcov = qreg_vcov(mr, mm.m, coef, q)
        stderror = sqrt.(diag(vcov))
        return TableRegressionModel(QRegModel(coef, vcov, stderror, q), mf, mm)
    end
    qreg(f::FormulaTerm, df::AbstractDataFrame, s::Solver) = qreg(f, df, 0.5, s)

    StatsBase.coef(x::QRegModel) = x.beta
    StatsBase.coef(rm::TableRegressionModel) = coef(rm.model)

    StatsBase.vcov(x::QRegModel) = x.vcov

    StatsBase.stderror(x::QRegModel) = x.stderr

    quantiles(x::QRegModel) = x.q
    quantiles(rm::TableRegressionModel) = quantiles(rm.model)

    function StatsBase.coeftable(mm::QRegModel)
        cc = coef(mm)
        se = stderror(mm)
        tt = cc./se
        CoefTable(hcat(kron(mm.q,ones(length(se))), cc,se,tt),
                 ["Quantile",  "Estimate","Std.Error","t value"],
                 ["x$i" for i = 1:length(cc)])
    end
end


text:
module DotNET

import Pkg.Artifacts:@artifact_str

export CLRObject,null,isnull,CLRException,@T_str,
    clrtypeof,isclrtype,isassignable,
    clreltype,
    delegate

include("typedef.jl")

include("CoreCLR.jl")
using .CoreCLR

include("CLRBridge.jl")
using .CLRBridge

include("typeinfo.jl")

include("marshalling.jl")

include("callback.jl")

include("reflection.jl")

include("operators.jl")

struct DummyCLRHost <: CLRHost end

const CURRENT_CLR_HOST = Ref{CLRHost}(DummyCLRHost())

function __init__()
    if CURRENT_CLR_HOST[] != DummyCLRHost() return end
    coreclr = detect_runtime(CoreCLRHost)
    if !isempty(coreclr)
        init_coreclr(first(coreclr))
        return
    end
    @error """
    No .NET Core runtime found on this system.
    Try specifying DOTNET_ROOT environment variable or adding 'dotnet' executable to PATH, then run 'DotNET.__init__()' again.
    Note that .NET Framework is current not supported.
    """
end

function init_coreclr(runtime)
    CoreCLR.init(runtime)
    dir = artifact"clrbridge"
    clrbridge = joinpath(dir, "CLRBridge.dll")
    if !isfile(clrbridge)
        error("""
        Artifact is present but CLRBridge.dll not found, possibly due to a broken package installation.
        You may need to delete the artifact directory '$dir' and try 'DotNET.__init__()' again.
        """)
    end
    tpalist = build_tpalist(dirname(runtime.path))
    push!(tpalist, clrbridge)
    CURRENT_CLR_HOST[] = create_host(CoreCLRHost;tpalist = tpalist)
    CLRBridge.init(CURRENT_CLR_HOST[])
end

function build_tpalist(dir)
    joinpath.(dir, filter(x->splitext(x)[2] == ".dll", readdir(dir)))
end

end # module


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Extension of ColorTypes.RGB for color manipulation

#############
# Operations

"""
    ==(c1::RGB, c2::RGB)

Elementwise comparison of two colors.
"""
(==)(c1::RGB, c2::RGB) = all(isequal.(c1, c2))

"""
    +(c1::RGB, c2::RGB)

Return the elementwise sum of two colors.

# Examples

```jldoctest
julia> RGB(1f0, 2f0, 3f0) + RGB(4f0, 5f0, 6f0)
RGB color with eltype Float32
R: 5.0, G: 7.0, B: 9.0
```
"""
(+)(c1::RGB, c2::RGB) = RGB((+).(c1, c2))

"""
    -(c1::RGB, c2::RGB)

Return the elementwise difference of two colors.

# Examples

```jldoctest
julia> RGB(1f0, 2f0, 3f0) - RGB(4f0, 5f0, 6f0)
RGB color with eltype Float32
R: -3.0, G: -3.0, B: -3.0
```
"""
(-)(c1::RGB, c2::RGB) = RGB((-).(c1, c2))

"""
    *(scalar::Number, c::RGB)
    *(c::RGB, scalar::Number)

Return a `RGB` color with each component multiplied by `scalar`.

# Examples

```jldoctest
julia> c = 2.0 * RGB(4f0, 5f0, 6f0)
RGB color with eltype Float32
R: 8.0, G: 10.0, B: 12.0

julia> RGB(4f0, 5f0, 6f0) * 2.0 == c
true
```

```jldoctest
julia> c = 2.0 * RGB(4.0, 5.0, 6.0)
RGB color with eltype Float64
R: 8.0, G: 10.0, B: 12.0

julia> RGB(4.0, 5.0, 6.0) * 2.0 == c
true
```

Note that the eltype of RGB is mantained.
"""
(*)(scalar::Number, c::RGB{T}) where {T} = RGB{T}((*).(scalar, c))
(*)(c::RGB, scalar::Number) = scalar * c

"""
    *(c1::RGB, c2::RGB)

Return the elementwise product of two colors.

# Examples

```jldoctest
julia> RGB(1f0, 2f0, 3f0) * RGB(4f0, 5f0, 6f0)
RGB color with eltype Float32
R: 4.0, G: 10.0, B: 18.0
```
"""
(*)(c1::RGB, c2::RGB) = RGB((*).(c1, c2))

"""
    ≈(c1::RGB, c2::RGB)

Check if two colors are close.

# Examples

```jldoctest
julia> c = RGB(1f0, 2f0, 3f0) * RGB(4f0, 5f0, 6f0)
RGB color with eltype Float32
R: 4.0, G: 10.0, B: 18.0

julia> c ≈ RGB(4f0, 10f0, 18f0)
true

julia> c ≈ RGB(0f0, 0f0, 0f0)
false
```
"""
(≈)(c1::RGB, c2::RGB) = all((≈).(c1, c2))

#############
# Iterations

length(::RGB) = 3
firstindex(::RGB) = 1
lastindex(c::RGB) = 3

# Since there is no standard that specifies the order in which the colors
# should be reported, here we use the convention whereby the colors should
# be reported in the RGB order (first R, then G, and finally B).
function getindex(c::RGB, i::Integer)
    if i == 1
        return c.r
    elseif i == 2
        return c.g
    elseif i == 3
        return c.b
    else
        throw(BoundsError(c, i))
    end
end

getindex(c::RGB, i::CartesianIndex{1}) = getindex(c, Tuple(i)[1])
# setindex! not implemented since RGB is immutable

iterate(c::RGB, state = 1) = state > 3 ? nothing : (c[state], state +1)

###############
# Broadcasting

# Broadcasting a function `f` on a series of arguments `xs` (ie. `f.(xs...)`) is equivalent to writing
# the following: `materialize(broadcasted(combine_styles(map(broadcastable, xs)), f, xs...))`
# Let's analyze what each of these functions does and why each of the following methods is needed

# Let's start with broadcastable: broadcasting maps its arguments so that each argument is transformed
# into a type supporting indexing and the method `axes`, and thus being able to be converted into an
# array-like type. We'll trick the broadcasting process into thinking this is the case for `RGB` so that
# `broadcastable` will return its argument when it is applied to an `RGB`

axes(::RGB) = (OneTo(3),)

broadcastable(c::RGB) = c
broadcastable(::Type{RGB}) = RGB # broadcastable is also applied to types

# Then the result of the previous mapping is passed onto the function `combine_styles`, which
# relies on calls to the constructor of `BroadcastStyle`. We need to create a `BroadcastStyle`
# exclusive to `RGB` so that we can then specialize other methods that are fed `BroadcastStyle`s

struct RGBBroadcastStyle <: BroadcastStyle end

# Then we specialize the constructor from instances of `RGB`

BroadcastStyle(::Type{<:RGB}) = RGBBroadcastStyle()

# And the constructor that combines the `RGBBroadcastStyle` with any other `BroadcastStyle`
# we want our style to have precedence over any other so that if an `RGB` type is present among
# the arguments of broadcasting the result will be of type `RGB`

BroadcastStyle(::RGBBroadcastStyle, ::BroadcastStyle) = RGBBroadcastStyle()

# The call to `materialize` returns a call to `copy`, which by default is specialized for array-like types
# we then have to implement a method that treats any `Broadcasted` type with our custom `RGBBroadcastStyle`
# in an appropriate way.

@inline function copy(bc::Broadcasted{RGBBroadcastStyle})
    ElType = combine_eltypes(bc.f, bc.args)
    if ElType <: Fractional # IF ElType instances can be stored in a RGB instance
        # the call to `convert` to a `Broadcasted` type of style `Nothing` computes
        # the result of the broadcasting and stores it into an array. splatting this
        # array into an `RGB` constructor gives us the desired result
        return RGB{ElType}(convert(Broadcasted{Nothing}, bc)...)
    else #IF ElType instances cannot be stored in a RGB instance
        return copy(convert(Broadcasted{Broadcast.DefaultArrayStyle{ElType}}, bc))
    end
end

################
# Miscellaneous

show(io::IO, c::RGB) = print(io, "($(c.r) $(c.g) $(c.b))")

function show(io::IO, ::MIME"text/plain", c::RGB{T}) where {T}
    print(io, "RGB color with eltype $T\n", "R: $(c.r), G: $(c.g), B: $(c.b)")
end

@doc raw"""
    luminosity(c::RGB)

Return the mean value between the maximum component and the minumum component of a color:

```math
\frac{\mathrm{max}(c) + \mathrm{min}(c)}{2}
```

# Examples

```jldoctest
julia> luminosity(RGB(1f0, 2f0, 3f0))
2.0f0
```
"""
luminosity(c::RGB) = (max(c...) + min(c...)) / 2

@doc raw"""
    clamp(c::RGB)

Return a clamped `RGB` color, with each component `x` obtained with the formula:

```math
\frac{x}{1 + x}
```

# Examples

```julia
julia> clamp(RGB(1f0, 2f0, 3f0))
RGB color with eltype Float32
R: 0.5, G: 0.6666667, B: 0.75
```
"""
clamp(c::RGB{T}) where {T} = RGB{T}(map(x -> x / (1f0 + x), c)...)

@doc raw"""
    γ_correction(c::RGB, γ::Float32)

Return a `RGB` color, with each component `x` corrected with the formula:

```math
x^{1 / \gamma}
```

# Examples

```julia
julia> c = RGB(1f0, 2f0, 3f0);

julia> γ_correction(c, 1f0)
RGB color with eltype Float32
R: 1.0, G: 2.0, B: 3.0

julia> γ_correction(c, 0.8f0)
RGB color with eltype Float32
R: 1.0, G: 2.3784142, B: 3.948222

julia> γ_correction(c, 2.4f0)
RGB color with eltype Float32
R: 1.0, G: 1.3348398, B: 1.580522
```
"""
γ_correction(c::RGB{T}, γ::Float32) where {T} = RGB{T}(map(x -> x^(1f0 / γ), c)...)

eltype(::RGB{T}) where {T} = T

function zero(T::Type{<:RGB})
    z = zero(eltype(T))
    RGB(z, z, z)
end

zero(c::RGB) = zero(typeof(c))

function one(T::Type{<:RGB})
    z = one(eltype(T))
    RGB(z, z, z)
end

one(c::RGB) = one(typeof(c))

"""
    BLACK

A `RGB{Float32}` representing a black color.
"""
const BLACK = zero(RGB{Float32})

"""
    WHITE

A `RGB{Float32}` representing a white color.
"""
const WHITE = one(RGB{Float32})

"""
    RED

A `RGB{Float32}` representing a red color.
"""
const RED = RGB(1f0, 0f0, 0f0)

"""
    GREEN

A `RGB{Float32}` representing a green color.
"""
const GREEN = RGB(0f0, 1f0, 0f0)

"""
    BLUE

A `RGB{Float32}` representing a blue color.
"""
const BLUE = RGB(0f0, 0f0, 1f0)

"""
    CYAN

A `RGB{Float32}` representing a cyan color.
"""
const CYAN = RGB(0f0, 1f0, 1f0)

"""
    MAGENTA

A `RGB{Float32}` representing a magenta color.
"""
const MAGENTA = RGB(1f0, 0f0, 1f0)

"""
    YELLOW

A `RGB{Float32}` representing a yellow color.
"""
const YELLOW = RGB(1f0, 1f0, 0f0)


text:
import ElectromagneticFields

import GeometricIntegrators.Equations: ODE, SODE

export guiding_center_4d_ode, guiding_center_4d_sode

function guiding_center_4d_periodicity(q, periodic=true)
    period = zeros(eltype(q), size(q,1))

    if periodic
        period[1:3] .= periodicity(q)
    end

    return period
end

transform_v(V) = (t,q,q̇,parameters) -> V(t, transform_q̃_to_q(t, q, parameters), q̇, parameters)

function guiding_center_4d_ode(qᵢ, params; periodic=true)
    ODE(transform_v(v), qᵢ; parameters=params, h=hamiltonian,
            periodicity=guiding_center_4d_periodicity(qᵢ, periodic))
end

function guiding_center_4d_sode(qᵢ, params; periodic=true)
    SODE(Tuple(transform_v(V) for V in (v₁, v₂, v₃, v₄, v₅, v₆)), qᵢ; parameters=params,
            periodicity=guiding_center_4d_periodicity(qᵢ, periodic))
end


text:
function Base.:!(bc::CS.BasicConstraint)
    if bc.fct != equal
        throw(ErrorException("!BasicConstraint is only implemented for !equal"))
    end
    if length(bc.indices) != 2
        throw(ErrorException("!BasicConstraint is only implemented for !equal with exactly 2 variables"))
    end
    bc.fct = not_equal
    return bc
end

function not_equal(com::CS.CoM, constraint::BasicConstraint; logs = true)
    indices = constraint.indices

    search_space = com.search_space
    # we always have only two variables
    v1 = search_space[indices[1]]
    v2 = search_space[indices[2]]
    fixed_v1 = isfixed(v1)
    fixed_v2 = isfixed(v2)
    if !fixed_v1 && !fixed_v2
        return true
    elseif fixed_v1 && fixed_v2
        if value(v1) == value(v2)
            return false
        end
        return true
    end
    # one is fixed and one isn't
    if fixed_v1
        prune_v = v2
        prune_v_idx = 2
        other_val = value(v1)
    else
        prune_v = v1
        prune_v_idx = 1
        other_val = value(v2)
    end
    if has(prune_v, other_val)
        if !rm!(com, prune_v, other_val)
            return false
        end
    end
    return true
end

"""
    not_equal(com::CoM, constraint::Constraint, value::Int, index::Int)

Returns whether the `not_equal` constraint can be still fulfilled.
"""
function not_equal(com::CoM, constraint::Constraint, value::Int, index::Int)
    if index == constraint.indices[1]
        other_var = com.search_space[constraint.indices[2]]
    else
        other_var = com.search_space[constraint.indices[1]]
    end
    return !issetto(other_var, value)
end


text:
"""
    ArrayOutput <: Output

    ArrayOutput(init; tspan::AbstractRange, [aux, mask, padval])

A simple output that stores each step of the simulation in a vector of arrays.

# Arguments

- `init`: initialisation `AbstractArrayArray` or `NamedTuple` of `AbstractArrayArray`.

# Keywords (passed to [`Extent`](@ref))

$EXTENT_KEYWORDS

An `Extent` object can be also passed to the `extent` keyword, and other keywords will be ignored.
"""
mutable struct ArrayOutput{T,F<:AbstractVector{T},E} <: Output{T,F}
    frames::F
    running::Bool
    extent::E
end
function ArrayOutput(; frames, running, extent, kw...)
    append!(frames, _zerogrids(init(extent), length(tspan(extent))-1))
    ArrayOutput(frames, running, extent)
end

"""
    ResultOutput <: Output

    ResultOutput(init; tspan::AbstractRange, kw...)

A simple output that only stores the final result, not intermediate frames.

# Arguments

- `init`: initialisation `Array` or `NamedTuple` of `Array`

# Keywords (passed to [`Extent`](@ref))

$EXTENT_KEYWORDS

An `Extent` object can be also passed to the `extent` keyword, and other keywords will be ignored.
"""
mutable struct ResultOutput{T,F<:AbstractVector{T},E} <: Output{T,F}
    frames::F
    running::Bool
    extent::E
end
ResultOutput(; frames, running, extent, kw...) = ResultOutput(frames, running, extent)

isstored(o::ResultOutput) = false
storeframe!(o::ResultOutput, data::AbstractSimData) = nothing

function finalise!(o::ResultOutput, data::AbstractSimData)
    # Only store after the last frame
    _storeframe!(o, eltype(o), data)
end


text:
#=
KTruss:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra

#A is the original matrix, k is the requested truss, n the nodes in the original matrix, Output: K-Truss graph
function ktruss(A, k, n)

    C = A
    print("STARTING ALGORITHM ...\n\n")
    for i = 1:n-1
       C = mul(C, C, Semirings.PLUS_LAND, mask=C)
       C = select(SelectOps.GE, C, k-2)
       print(C)
       #***
       print("\n\n\n")
       print("ITERATION FINISHED")
       print("\n\n\n")
    end

    return C

end

#***
#optional: we could also tell the algorithm to stop if the non zero values didn't change from the previous iteration


text:
function _prepare(nodes::Matrix{T},
                  kernel::RK
                 ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    n = size(nodes, 1)
    n_1 = size(nodes, 2)
    min_bound = Vector{T}(undef, n)
    compression::T = 0
    @inbounds for i = 1:n
        min_bound[i] = nodes[i,1]
        maxx::T = nodes[i,1]
        for j = 2:n_1
            min_bound[i] = min(min_bound[i], nodes[i,j])
            maxx = max(maxx, nodes[i,j])
        end
        compression = max(compression, maxx - min_bound[i])
    end
    if compression <= eps(T(1.0))
        error("Cannot prepare the spline: `nodes` data are not correct.")
    end

    t_nodes = similar(nodes)
    @inbounds for j = 1:n_1
        for i = 1:n
            t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression
        end
    end

    if T(kernel.ε) == T(0.0)
        ε = _estimate_ε(t_nodes)
        if isa(kernel, RK_H0)
            kernel = RK_H0(ε)
        elseif isa(kernel, RK_H1)
            ε *= T(1.5)
            kernel = RK_H1(ε)
        elseif isa(kernel, RK_H2)
            ε *= T(2.0)
            kernel = RK_H2(ε)
        else
            error("incorrect `kernel` type.")
        end
    end

    gram = _gram(t_nodes, kernel)
    chol = nothing
    try
        chol = cholesky(gram)
    catch
        error("Cannot prepare the spline: Gram matrix is degenerate.")
    end

    cond = _estimate_cond(gram, chol)

    spline = NormalSpline(kernel,
                          compression,
                          t_nodes,
                          nothing,
                          nothing,
                          nothing,
                          nothing,
                          min_bound,
                          gram,
                          chol,
                          nothing,
                          cond
                         )
    return spline
end

function _construct(spline::NormalSpline{T, RK},
                    values::Vector{T},
                    cleanup::Bool = false
                   ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    if(length(values) != size(spline._nodes, 2))
        error("Number of data values does not correspond to the number of nodes.")
    end
    if isnothing(spline._chol)
        error("Gram matrix was not factorized.")
    end

    mu = Vector{T}(undef, size(spline._gram, 1))
    ldiv!(mu, spline._chol, values)

    spline = NormalSpline(spline._kernel,
                          spline._compression,
                          spline._nodes,
                          values,
                          nothing,
                          nothing,
                          nothing,
                          spline._min_bound,
                          cleanup ? nothing : spline._gram,
                          cleanup ? nothing : spline._chol,
                          mu,
                          spline._cond
                         )
    return spline
end

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

function _prepare(nodes::Matrix{T},
                  d_nodes::Matrix{T},
                  es::Matrix{T},
                  kernel::RK
                 ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
    n = size(nodes, 1)
    n_1 = size(nodes, 2)
    n_2 = size(d_nodes, 2)

    if(size(es, 2) != n_2)
        error("Number of derivative directions does not correspond to the number of derivative nodes.")
    end

    t_es = similar(es)
    try
        @inbounds for i = 1:n_2
            t_es[:,i] = es[:,i] ./ norm(es[:,i])
        end
    catch
        error("Cannot normalize derivative direction: zero direction vector.")
    end

    min_bound = Vector{T}(undef, n)
    compression::T = 0
    @inbounds for i = 1:n
        min_bound[i] = nodes[i,1]
        maxx::T = nodes[i,1]
        for j = 2:n_1
            min_bound[i] = min(min_bound[i], nodes[i,j])
            maxx = max(maxx, nodes[i,j])
        end
        for j = 1:n_2
            min_bound[i] = min(min_bound[i], d_nodes[i,j])
            maxx = max(maxx, d_nodes[i,j])
        end
        compression = max(compression, maxx - min_bound[i])
    end

    if compression <= eps(T(1.0))
        error("Cannot prepare the spline: `nodes` data are not correct.")
    end

    t_nodes = similar(nodes)
    t_d_nodes = similar(d_nodes)

    @inbounds for j = 1:n_1
        for i = 1:n
            t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression
        end
    end

    @inbounds for j = 1:n_2
        for i = 1:n
            t_d_nodes[i,j] = (d_nodes[i,j] - min_bound[i]) / compression
        end
    end

    if T(kernel.ε) == T(0.0)
        ε = _estimate_ε(t_nodes, t_d_nodes)
        if isa(kernel, RK_H1)
            ε *= T(2.0)
            kernel = RK_H1(ε)
        elseif isa(kernel, RK_H2)
            ε *= T(2.5)
            kernel = RK_H2(ε)
        else
            error("incorrect `kernel` type.")
        end
    end

    gram = _gram(t_nodes, t_d_nodes, t_es, kernel)
    chol = nothing
    try
        chol = cholesky(gram)
    catch
        error("Cannot prepare the spline: Gram matrix is degenerate.")
    end

    cond = _estimate_cond(gram, chol)

    spline = NormalSpline(kernel,
                          compression,
                          t_nodes,
                          nothing,
                          t_d_nodes,
                          t_es,
                          nothing,
                          min_bound,
                          gram,
                          chol,
                          nothing,
                          cond
                         )
    return spline
end

function _construct(spline::NormalSpline{T, RK},
                    values::Vector{T},
                    d_values::Vector{T},
                    cleanup::Bool = false
                   ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    if(length(values) != size(spline._nodes, 2))
        error("Number of data values does not correspond to the number of nodes.")
    end
    if(length(d_values) != size(spline._d_nodes, 2))
        error("Number of derivative values does not correspond to the number of derivative nodes.")
    end
    if isnothing(spline._chol)
        error("Gram matrix was not factorized.")
    end

    mu = Vector{T}(undef, size(spline._gram, 1))
    ldiv!(mu, spline._chol, [values; spline._compression .* d_values])

    spline = NormalSpline(spline._kernel,
                          spline._compression,
                          spline._nodes,
                          values,
                          spline._d_nodes,
                          spline._es,
                          d_values,
                          spline._min_bound,
                          cleanup ? nothing : spline._gram,
                          cleanup ? nothing : spline._chol,
                          mu,
                          spline._cond
                         )
    return spline
end

function _evaluate(spline::NormalSpline{T, RK},
                   points::Matrix{T},
                   do_parallel::Bool = false
                  ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    if isnothing(spline)
        error("Spline was not prepared.")
    end
    if isnothing(spline._mu)
        error("Spline coefficients were not calculated.")
    end

    if size(points, 1) != size(spline._nodes, 1)
        if size(points, 1) == 1 && size(points, 2) > 1
            error("Incorrect first dimension of the `points` parameter (use 'evaluate_one' function for evaluating the spline at one point).")
        else
            error("Incorrect first dimension of the `points` parameter (the spline was built in the space of different dimension).")
        end
    end

    n = size(spline._nodes, 1)
    n_1 = size(spline._nodes, 2)
    m = size(points, 2)

    pts = similar(points)
    @inbounds for j = 1:m
        for i = 1:n
            pts[i,j] = (points[i,j] - spline._min_bound[i]) / spline._compression
        end
    end

    spline_values = Vector{T}(undef, m)
    iend1 = iend2 = iend3 = 0
    istart2 = istart3 = istart4 = 0
    mu = spline._mu[1:n_1]
    d_mu = spline._mu[(n_1 + 1):end]
    if do_parallel && m >= 1000 && Threads.nthreads() >= 4
        step = m ÷ 4
        iend1 = 1 + step
        istart2 = iend1 + 1
        iend2 = istart2 + step
        istart3 = iend2 + 1
        iend3 = (istart3 + step) < m ? (istart3 + step) : m
        istart4 = iend3 + 1
        @inbounds Threads.@threads for t = 1:4
            if t == 1
                _do_work(1, iend1, pts, spline._nodes, mu, spline._kernel, spline_values)
                if !isnothing(spline._d_nodes)
                    _do_work_d(1, iend1, pts, spline._d_nodes, spline._es, d_mu, spline._kernel, spline_values)
                end
            elseif t == 2
                _do_work(istart2, iend2, pts, spline._nodes, mu, spline._kernel, spline_values)
                if !isnothing(spline._d_nodes)
                    _do_work_d(istart2, iend2, pts, spline._d_nodes, spline._es, d_mu, spline._kernel, spline_values)
                end
            elseif t == 3
                _do_work(istart3, iend3, pts, spline._nodes, mu, spline._kernel, spline_values)
                if !isnothing(spline._d_nodes)
                    _do_work_d(istart3, iend3, pts, spline._d_nodes, spline._es, d_mu, spline._kernel, spline_values)
                end
            elseif t == 4
                if istart4 <= m
                    _do_work(istart4, m, pts, spline._nodes, mu, spline._kernel, spline_values)
                    if !isnothing(spline._d_nodes)
                        if istart4 <= m
                            _do_work_d(istart4, m, pts, spline._d_nodes, spline._es, d_mu, spline._kernel, spline_values)
                        end
                    end
                end
            end
        end
    else
        _do_work(1, m, pts, spline._nodes, mu, spline._kernel, spline_values)
        if !isnothing(spline._d_nodes)
            _do_work_d(1, m, pts, spline._d_nodes, spline._es, d_mu, spline._kernel, spline_values)
        end
    end

    return spline_values
end

function _do_work(istart::Int,
                  iend::Int,
                  points::Matrix{T},
                  nodes::Matrix{T},
                  mu::Vector{T},
                  kernel::RK,
                  spline_values::Vector{T}
                ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    n_1 = size(nodes, 2)
    h_values = Vector{T}(undef, n_1)
    @inbounds for p = istart:iend
        for i = 1:n_1
            h_values[i] = _rk(kernel, points[:,p], nodes[:,i])
        end
        spline_values[p] = sum(mu .* h_values)
    end
end

function _do_work_d(istart::Int,
                    iend::Int,
                    pts::Matrix{T},
                    d_nodes::Matrix{T},
                    es::Matrix{T},
                    d_mu::Vector{T},
                    kernel::RK,
                    spline_values::Vector{T}
                   ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
   n_2 = size(d_nodes, 2)
   d_h_values = Vector{T}(undef, n_2)
   @inbounds for p = istart:iend
       for i = 1:n_2
           d_h_values[i] = _∂rk_∂e(kernel, pts[:,p], d_nodes[:,i], es[:,i])
       end
       spline_values[p] += sum(d_mu .* d_h_values)
    end
end

function _evaluate_gradient(spline::NormalSpline{T, RK},
                            point::Vector{T}
                           ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    if isnothing(spline._mu)
        error("Spline coefficients were not calculated.")
    end

    n = size(spline._nodes, 1)
    n_1 = size(spline._nodes, 2)

    pt = Vector{T}(point)
    @inbounds for i = 1:n
        pt[i] = (point[i] - spline._min_bound[i]) / spline._compression
    end

    d_h_values = Vector{T}(undef, n_1)
    grad = Vector{T}(undef, n)
    mu = spline._mu[1:n_1]
    @inbounds for k = 1:n
        for i = 1:n_1
            d_h_values[i] = _∂rk_∂η_k(spline._kernel, pt, spline._nodes[:,i], k)
        end
        grad[k] = sum(mu .* d_h_values)
    end

    if !isnothing(spline._d_nodes)
        n_2 = size(spline._d_nodes, 2)
        d_h_values = Vector{T}(undef, n_2)
        d_mu = spline._mu[n_1+1:end]
        @inbounds for k = 1:n
            for i = 1:n_2
                d_h_values[i] = T(0.0)
                for l = 1:n
                    d_h_values[i] += (_∂²rk_∂η_r_∂ξ_k(spline._kernel, pt, spline._d_nodes[:,i], k, l) * spline._es[l,i])
                end
            end
            grad[k] += sum(d_mu .* d_h_values)
        end
    end
    @inbounds for k = 1:n
        grad[k] /= spline._compression
    end
    return grad
end


text:
################################################################################
#
#  Small groups
#
################################################################################

export number_of_small_groups, small_group, small_groups_limit

const small_groups_limit = 63

include("small_groups_1_63")

function number_of_small_groups(i)
  return length(small_groups_1_63[i])
end

isfrom_db(G::GrpGen) = G.isfromdb

# Data in the DB:
# <gens, rels, nontrivrels, orderdis, ordersubdis, IsAbelian(G), IsCyclic(G),
# IsSolvable(G), IsNilpotent(G), #AutomorphismGroup(G);

function small_group(i, j)
  i < 1 || i > 63 && error("Group order ($i) must be between 1 and $small_groups_limit")
  j < 1 || j > number_of_small_groups(i) && error("Index ($j) must be between 1 and $(number_of_small_groups(i))")
  P = PermGroup(i)
  Gens = [P()]
  for p in small_groups_1_63[i][j][1]
    push!(Gens, P(p))
  end
  G, PtoG, GtoP = generic_group(closure(Gens, *, P()), *)

  # The small groups have generators, for which we know
  # generators for the relations.
  # This is used in the homomorphism computation.
  # To make this work, the generators must much exactly
  # with the generators in the database.
  G.gens = Vector{Int}(undef, length(small_groups_1_63[i][j][1]))
  for k in 1:length(small_groups_1_63[i][j][1])
    p = small_groups_1_63[i][j][1][k]
    G.gens[k] = PtoG[P(p)][]
  end

  G.issolvable = small_groups_1_63[i][j][8] ? 1 : 2
  G.isnilpotent = small_groups_1_63[i][j][9] ? 1 : 2
  G.isfromdb = true
  G.small_group_id = (i, j)
  return G
end


text:
using Reactive

import Base: IndexStyle, size, getindex, setindex!, _unsafe_getindex!, push!, similar

using Base.Broadcast: map_newindexer, _broadcast_eltype, broadcast_indices

struct ArrayNode{T, N, F, AT <: AbstractArray} <: AbstractArray{T, N}
    signal::Signal{AT}
    convert::F
end
const AbstractNode{T} = Union{ArrayNode{T}, Node{T}}

function Base.show(io::IO, node::AbstractNode{T}) where T
    print(io, "<Node: ")
    show(io, to_value(node))
    print(io, ">")
end

function children(x::Signal)
    filter(x-> x != nothing, (x-> x.value).(Reactive.nodes[Reactive.edges[x.id]]))
end
children(x::AbstractNode) = children(to_signal(x))

function disconnect!(x::Signal, toremove::Signal)
    x.parents = (Iterators.filter(x-> x != toremove, x.parents)...,)
end

disconnect_all_children!(x::AbstractNode) = disconnect_all_children!(to_signal(x))

function disconnect_all_children!(x::Signal)
    for child in children(x)
        disconnect!(child, x)
    end
end

function closenode!(x::AbstractNode)
    s = to_signal(x)
    disconnect_all_children!(s)
    close(s, false)
end

"""
Lift nodes of a scene by supplying fetching nodes via keys into the scene.
A tuple of symboles will be interpreted as `scene[:tupl_elem_1, :tupl_elem_2, ...]`
"""
function lift_node(f, scene::Scene, keys::Union{Symbol, Tuple}...)
    args = map(x-> getindex(scene, x), keys)
    lift_node(f, args...)
end

"""
Registers a callback to `nodes`, which calls function `f` whenever any node in `nodes` updates.
`f` will get the values of each `node` as an argument, so basically `f(to_value.(nodes))`.
Returns a new node which is the result of `f` applied to the updated `nodes`.
"""
function lift_node(f, nodes::AbstractNode...)
    args = to_signal.(nodes)
    s = foreach(f, args...)
    to_node(s)
end

function push!(x::AbstractNode, value)
    val = x.convert(to_value(value))
    s = to_signal(x)
    Reactive.set_value!(s, val) # make the value available already!
    push!(s, val)
    val
end

IndexStyle(x::ArrayNode{T, N, F, AT}) where {T, N, F, AT} = IndexStyle(AT)
Array(x::ArrayNode) = value(to_signal(x))
similar(A::ArrayNode) = to_node(similar(Array(A)))
similar(A::ArrayNode, ::Type{T}) where T = to_node(similar(Array(A), T))
similar(A::ArrayNode, dims::NTuple{N, Int}) where N = to_node(similar(Array(A), dims))
similar(A::ArrayNode, ::Type{T}, dims::NTuple{N, Int}) where {T, N} = to_node(similar(Array(A), T, dims))

size(x::ArrayNode) = size(value(to_signal(x)))

# Index Linear
function getindex(x::ArrayNode, i::Integer)
    Array(x)[i]
end
function setindex!(x::ArrayNode{T, N}, val, i::Integer) where {T, N}
    Array(x)[i] = val
    push!(to_signal(x), Array(x)) # update array
    val
end
function getindex(x::ArrayNode, i::Vararg{Int, N}) where N
    Array(x)[i...]
end

function setindex!(x::ArrayNode{T, N}, val, i::Vararg{Int, N}) where {T, N}
    Array(x)[i...] = val
    push!(to_signal(x), Array(x))
    val
end

function _unsafe_getindex!(dest::ArrayNode, src::AbstractArray, Is::Union{Real, AbstractArray}...)
    Array(dest)[Is...] = Array(src)
    push!(to_signal(dest), Array(dest))
    return dest
end

to_node(obj::AbstractNode) = obj
function to_node(obj::AbstractNode, f)
    to_node(map(f, to_signal(obj)), f)
end
to_node(obj, f = identity) = to_node(Signal(f(obj)), f)
to_node(obj::Signal, f = identity) = Node(map(f, obj), f)
function to_node(obj::Signal{AT}, f::F = identity) where {AT <: AbstractArray, F}
    A = value(obj)
    ArrayNode{eltype(A), ndims(A), F, AT}(obj, f)
end

to_value(obj::AbstractNode) = value(to_signal(obj))
to_value(obj) = obj

to_signal(obj::AbstractNode) = obj.signal
to_signal(obj) = obj

# Called by Base broadcasting mechanisms (in place and out of place)
Base.Broadcast._containertype(::Type{<:ArrayNode}) = ArrayNode
Base.Broadcast.promote_containertype(::Type{Any}, ::Type{ArrayNode}) = ArrayNode
Base.Broadcast.promote_containertype(::Type{ArrayNode}, ::Type{Any}) = ArrayNode

@inline function Base.Broadcast.broadcast_c!(f, ::Type{ArrayNode}, ::Type, C, A, Bs::Vararg{Any, N}) where N
    lift_node(to_node(C), to_node(A), to_node.(Bs)...) do c, a, bs...
        broadcast!(f, c, a, bs)
    end
    C
end

@inline function Base.Broadcast.broadcast_c(f, ::Type{ArrayNode}, A, Bs...)
    lift_node(to_node(A), to_node.(Bs)...) do a, bs...
        broadcast(f, a, bs...)
    end
end
function Base.append!(node::ArrayNode, values)
    v = to_value(node)
    append!(v, values)
    push!(node, v)
    node
end


text:
"""
    abstract type ProjectiveTransform <: Transform

Abstract supertype for projective transformations. See
[Projective transformations](../docs/projective/interface.md).
"""
abstract type ProjectiveTransform <: Transform end

"""
    getprojection(tfm, bounds; randstate)

Create a projection for an item with spatial bounds `bounds`.
The projection should be a `CoordinateTransformations.Transformation`.
See [CoordinateTransformations.jl](https://github.com/JuliaGeometry/CoordinateTransformations.jl)
"""
function getprojection end

struct Bounds{N}
    rs::NTuple{N, UnitRange{Int}}
end

function Base.show(io::IO, bounds::Bounds{N}) where N
    print(io, "Bounds(")
    for (i, r) in enumerate(bounds.rs)
        print(io, first(r), ':', last(r))
        i == N || print(io, '×')
    end
    print(io, ")")
end

Bounds(sz::NTuple{N, <:Integer}) where N = Bounds(Tuple(1:n for n in sz))
Bounds(axes::NTuple{N, <:Base.OneTo}) where N = Bounds(convert.(UnitRange{Int}, axes))

Base.:(==)(bs1::Bounds, bs2::Bounds) = false
function Base.:(==)(bs1::Bounds{N}, bs2::Bounds{N}) where N
    c1 = all((first(bs1.rs[i]) == first(bs2.rs[i])) for i in 1:N)
    c2 = all((last(bs1.rs[i]) == last(bs2.rs[i])) for i in 1:N)
    return c1 && c2
end

"""
    transformbounds(bounds, P)

Apply `CoordinateTransformations.Transformation` to `bounds`.
"""
function transformbounds(bounds::Bounds, P::CoordinateTransformations.Transformation)
    return Bounds(ImageTransformations.autorange(CartesianIndices(bounds.rs), P))
end

"""
    getbounds(item)

Return the spatial bounds of `item`. For a 2D-image (`Image{2}`)
the bounds are the 4 corners of the bounding rectangle. In general,
for an N-dimensional item, the bounds are a vector of the N^2 corners
of the N-dimensional hypercube bounding the data.
"""
function getbounds end

getbounds(wrapper::ItemWrapper) = getbounds(getwrapped(wrapper))

"""
    project(P, item, indices)

Project `item` using projection `P` and crop to `indices` if given.
"""
function project end

"""
    project!(bufitem, P, item, indices)

Project `item` using projection `P` and crop to `indices` if given.
Store result in `bufitem`. Inplace version of [`project`](#).

Default implementation falls back to `project`.
"""
function project!(bufitem, P, item, indices)
    titem = project(P, item, indices)
    copyitemdata!(bufitem, titem)

    return bufitem
end

"""
    projectionbounds(tfm, P, bounds, randstate)
"""
function projectionbounds(tfm, P, bounds; randstate = getrandstate(tfm))
    return transformbounds(bounds, P)
end

# With the interface defined, any `ProjectiveTransform` can be `apply`ed
# like this:

function apply(tfm::ProjectiveTransform, item::Item; randstate = getrandstate(tfm))
    bounds = getbounds(item)
    P = getprojection(tfm, bounds; randstate = randstate)
    bounds_ = projectionbounds(tfm, P, bounds; randstate = randstate)
    return P isa IdentityTransformation ? item : project(P, item, bounds_)
end

# For the buffered version, `project!` is used. Of course the size
# of the data produced by `tfm` needs to be the same every time it
# is applied for this to work.

function apply!(
        bufitem::AbstractItem,
        tfm::ProjectiveTransform,
        item::AbstractItem;
        randstate = getrandstate(tfm))
    bounds = getbounds(item)
    P = getprojection(tfm, bounds; randstate = randstate)
    bounds_ = projectionbounds(tfm, P, bounds; randstate = randstate)
    res = project!(bufitem, P, item, bounds_)
    return res
end

# The simplest `ProjectiveTransform` is the static `Project`.

struct Project{T<:CoordinateTransformations.Transformation} <: ProjectiveTransform
    P::T
end

getprojection(tfm::Project, bounds; randstate = nothing) = tfm.P


text:
using NDTensors,
      Test
using LinearAlgebra

@testset "random_orthog" begin
  n,m = 10,4
  O1 = random_orthog(n,m)
  @test norm(transpose(O1)*O1 - Diagonal(fill(1.,m))) < 1E-14

  n,m = 4,10
  O2 = random_orthog(n,m)
  @test norm(O2*transpose(O2) - Diagonal(fill(1.,n))) < 1E-14
end

@testset "random_unitary" begin
  n,m = 10,4
  U1 = random_unitary(n,m)
  @test eltype(U1) <: ComplexF64
  @test norm(U1'*U1 - Diagonal(fill(1.,m))) < 1E-14

  n,m = 4,10
  U2 = random_unitary(n,m)
  @test norm(U2*U2' - Diagonal(fill(1.,n))) < 1E-14
end

nothing


text:
# Sum-Product Networks as described by Poon&Domingos, Gens&Domingos etc.

abstract SPNode{T} # SPN node with names of type T

immutable Sum{T} <: SPNode{T}
  children::Vector{SPNode{T}}
  weights::Vector{Float64}
end

immutable Prod{T} <: SPNode{T}
  children::Vector{SPNode{T}}
end

# Variables have names of a given type.
immutable Var{T} <: SPNode{T}
  name::T
end

immutable NVar{T} <: SPNode{T}
  name::T
  negname::T
end
negname(name::Symbol) = symbol("__N__" * string(name))
negname{I<:Integer}(name::I) = -name
NVar(name) = NVar{typeof(name)}(name, negname(name))

# A SPN with variable names of type T.
type SumProductNetwork{T}
  root::SPNode{T}
end
SPN{T}(r::SPNode{T}) = SumProductNetwork{T}(r)

# Direct calculation of probabilities
calc{T, N <: Number}(spn::SumProductNetwork{T}, vals::Dict{T,N}) = calc(spn.root, vals)
function calc{T, N <: Number}(s::Sum{T}, vals::Dict{T,N})
  sum(calcchildren(s.children, vals) .* s.weights)
end
calcchildren{T, N <: Number}(children::Vector{SPNode{T}}, vals::Dict{T,N}) = map(c->calc(c, vals), children)
calc{T, N <: Number}(s::Prod{T}, vals::Dict{T,N}) = prod(calcchildren(s.children, vals))
calc{T, N <: Number}(v::Var{T}, vals::Dict{T,N}) = haskey(vals, v.name) ? vals[v.name] : 1
function calc{T, N <: Number}(v::NVar{T}, vals::Dict{T,N})
  haskey(vals, v.negname) ? vals[v.negname] : (haskey(vals, v.name) ? (1-vals[v.name]) : 1)
end

# But we also want to be able to calculate probabilities
# when values of some variables is not known. This is done by setting both
# the vars and negvars of the unknown variables to 1. Eas

# Simple examples
s1 = SPN(Var(1))
calc(s1, Dict{Int64,Int64}(1 => 1))

s2 = SPN(NVar(1))
calc(s2, Dict{Int64,Int64}(1 => 1))
calc(s2, Dict{Int64,Int64}(1 => 0))
calc(s2, Dict{Int64,Int64}(-1 => 1))
calc(s2, Dict{Int64,Int64}(-1 => 0))
calc(s2, Dict{Int64,Int64}())

s3 = SPN(Var(:a))
calc(s3, Dict{Symbol,Int64}(:a => 1))

# Fig 1 from cs.ubc.ca paper:
at = Sum(SPNode{Symbol}[Var(:a), NVar(:a)], Float64[0.5, 0.5])
bt = Sum(SPNode{Symbol}[Var(:b), NVar(:b)], Float64[0.5, 0.5])
ct = Sum(SPNode{Symbol}[Var(:c), NVar(:c)], Float64[0.5, 0.5])
dt = Sum(SPNode{Symbol}[Var(:d), NVar(:d)], Float64[0.5, 0.5])
left = Prod(SPNode{Symbol}[at, bt])
right = Prod(SPNode{Symbol}[ct, dt])
r = Sum(SPNode{Symbol}[left, right], Float64[0.5, 0.5])
@time calc(r, Dict{Symbol,Int64}(:a => 1, :b => 0, :c => 1, :d => 1))
calc(r, Dict{Symbol,Int64}(:a => 1, :b => 0, :c => 1))


text:
"""
    Bitmap2D(grid)

Create an empty `Bitmap2D` from an grid.
"""
Bitmap2D(grid::Grid2D) = Bitmap2D(falses(size(grid)...), grid)

Base.copy(M::Bitmap2D) = Bitmap2D(copy(M.data), M.grid)
Base.length(M::Bitmap2D) = length(M.grid)
Base.size(M::Bitmap2D) = size(M.grid)
Base.size(M::Bitmap2D, i::Int64) = Base.size(M.grid, i)
Base.ndims(M::Bitmap2D) = 2
Base.getindex(M::Bitmap2D, i) = getindex(M.data, i)
Base.getindex(M::Bitmap2D, i, j) = getindex(M.data, i, j)
Base.setindex!(M::Bitmap2D, x::Bool, i) = setindex!(M.data, x, i)
Base.setindex!(M::Bitmap2D, x::Bool, i, j) = setindex!(M.data, x, i, j)

Base.intersect(M₁::Bitmap2D, M₂::Bitmap2D) = Bitmap2D(M₁.data .& M₂.data, M₂.grid)

function safe_getindex(M::Bitmap2D, i, j)
    if 1 ≤ i ≤ size(M.data, 1) && 1 ≤ j ≤ size(M.data, 2)
        return getindex(M.data, i, j)
    end
    return false
end

function safe_setindex!(M::Bitmap2D, x::Bool, i, j)
    if 1 ≤ i ≤ size(M.data, 1) && 1 ≤ j ≤ size(M.data, 2)
        setindex!(M, x, i, j)
    end
end

@recipe function plot(M::Bitmap2D; transparent=true, scatter=false)
    seriescolor --> :dodgerblue

    if scatter
        xs, ys = _scatterdata(M)

        @series begin
            st = :scatter
            xs, ys
        end

    else
        clr = Colors.parse(Colors.Colorant, plotattributes[:seriescolor])
        cmap = PlotUtils.ColorGradient([clr, clr])
        fillcolor := cmap
        colorbar := false
        aspect_ratio --> :equal
        data = map(x -> x ? 1.0 : NaN, M.data)
        if !transparent
            cmap = PlotUtils.ColorGradient([Colors.parse(Colors.Colorant, :white), clr])
            fillcolor := cmap
            data = map(x -> x ? 1.0 : 0.0, M.data)
        end

        (xmin, xmax, ymin, ymax) = limits(M.grid)

        @series begin
            st := :heatmap
            grid  --> false
            xlims --> (xmin, xmax)
            ylims --> (ymin, ymax)
            xcoordinates(M.grid), ycoordinates(M.grid), data'
        end
    end
end

function _scatterdata(D::Bitmap2D)
    xs = Float64[]
    ys = Float64[]

    for l in eachindex(D.data)
        if D[l]
            i, j = ind2sub(D.data, l)
            push!(xs, D.grid.xrange[i])
            push!(ys, D.grid.yrange[j])
        end
    end
    (xs, ys)
end

fillrect!(M::Bitmap2D, px, py) = fillrect!(M.data, px, py)

#
# function Base.show(io::IO, M::Bitmap2D)
#     #println(typeof(A), ":")
#     println(io, "$(size(M, 1))×$(size(M, 2)) Bitmap2D with window [$(M.xmin), $(M.xmax), $(M.ymin), $(M.ymax)]")
# end


text:
@testset "basic message handling" begin
    cb = 0  # callback count
    function handle(type, args...)
        Atom.handlemsg(Dict("type"     => type,
                            "callback" => (cb += 1)),
                       args...)
    end
    function expect(expected)
        @test readmsg() == ["cb", cb, expected]
    end

    # pingpong
    handle("ping")
    expect("pong")

    # echo
    handle("echo", "echome!")
    expect("echome!")

    # cd
    oldpath = pwd()
    handle("cd", joinpath(oldpath, ".."))
    expect(nothing)
    @test pwd() == realpath(joinpath(oldpath, ".."))
    cd(oldpath)

    # evalsimple
    handle("evalsimple", "1 + 1")
    expect(2)
    handle("evalsimple", "sin(pi)")
    expect(sin(pi))
end


text:
export Mesh, get_ngauss, get_area, Mesh_v2, Mesh_v3

macro exported_enum(name, args...)
    esc(quote
        @enum($name, $(args...))
        export $name
        $([:(export $arg) for arg in args]...)
        end)
end

@exported_enum FiniteElementType P1 P2 BDM1

@doc raw"""
`Mesh` holds data structures for an unstructured mesh.

- `nodes`: a $n_v \times 2$ coordinates array
- `edges`: a $n_{\text{edge}} \times 2$ integer array for edges
- `elems`: a $n_e \times 3$ connectivity matrix, 1-based.
- `nnode`, `nedge`, `nelem`: number of nodes, edges, and elements
- `ndof`: total number of degrees of freedoms
- `conn`: connectivity matrix, `nelems × 3` or `nelems × 6`, depending on whether a linear element or a quadratic element is used.
- `lorder`: order of quadrature rule for line integrals (default = 6, 4 gauss points per line segment)
- `elem_type`: type of the element (P1, P2 or BDM1)

# Constructors

```
Mesh(m::Int64, n::Int64, h::Float64; order::Int64 = -1,
            degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = -1,
            version::Int64 = 1)
```

Constructs a mesh of a rectangular domain. The rectangle is split into $m\times n$ cells, and each cell is further split into two triangles.
`order` specifies the quadrature rule order. `degree` determines the degree for finite element basis functions.

!!! info
    AdFem provides three types of triangulations for a rectangular domain. The different types of meshes can be used to validate numerical schemes.
    For example, we can change to different meshes to verify that bugs of our program do not originate from mesh types.

    ![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/AdFem/mesh_types.png?raw=true)

# Integration Order : `order` and `lorder`

TODO
"""
mutable struct Mesh
    nodes::Array{Float64, 2}
    edges::Array{Int64, 2}
    elems::Array{Int64, 2}
    nnode::Int64
    nedge::Int64
    nelem::Int64
    ndof::Int64
    conn::Array{Int64, 2}
    lorder::Int64
    elem_ndof::Int64
    elem_type::FiniteElementType
end

function Mesh(coords::Array{Float64, 2}, elems::Array{Int64, 2}, order::Int64 = -1,
        degree::Union{FiniteElementType,Int64} = 1, lorder::Int64 = -1)
    if degree==P1
        degree = 1
    elseif degree==P2
        degree = 2
    end
    if !(degree in [1, 2]) && degree!=BDM1
        error("Only degree = 1 or 2 is supported.")
    end

    if order==-1
        if degree == 1 || degree==BDM1
            order = 2
        elseif degree == 2
            order = 4
        end
    end

    if lorder==-1
        if degree == 1 || degree==BDM1
            lorder = 6
        elseif degree == 2
            lorder = 6
        end
    end

    if degree==BDM1
        degree = -1
    end
    nnode = size(coords, 1)
    nelem = size(elems,1)
    nedges = zeros(Int64, 1)
    c = [coords zeros(size(coords, 1))]'[:]
    e = Int32.(elems'[:].- 1)
    edges_ptr = @eval ccall((:init_nnfem_mesh, $LIBMFEM), Ptr{Clonglong}, (Ptr{Cdouble}, Cint,
            Ptr{Cint}, Cint, Cint, Cint, Cint, Ptr{Clonglong}), $c, Int32(size($coords, 1)), $e, Int32(size($elems,1)),
            Int32($order), Int32($lorder), Int32($degree), $nedges)
    nedges = nedges[1]
    edges = unsafe_wrap(Array{Int64,1}, edges_ptr, (2nedges,), own=true)
    edges = reshape(edges, nedges, 2)
    elem_dof = Int64(@eval ccall((:mfem_get_elem_ndof, $LIBMFEM), Cint, ()))
    conn = zeros(Int64, elem_dof * size(elems, 1))
    @eval ccall((:mfem_get_connectivity, $LIBMFEM), Cvoid, (Ptr{Clonglong}, ), $conn)
    ndof = Int64(@eval ccall((:mfem_get_ndof, $LIBMFEM), Cint, ()))
    conn = reshape(conn, elem_dof, size(elems, 1))'|>Array
    @eval ccall((:mfem_get_element_to_vertices, $LIBMFEM), Cvoid, (Ptr{Clonglong}, ), $elems)

    elem_type = missing
    if degree==1
        elem_type = P1
    elseif degree==2
        elem_type = P2
    elseif degree==-1
        elem_type = BDM1
    end
    Mesh(coords, edges,  elems, nnode, nedges, nelem, ndof, conn, lorder, elem_dof, elem_type)
end

Base.:copy(mesh::Mesh) = Mesh(copy(mesh.nodes),
                            copy(mesh.edges),
                            copy(mesh.elems),
                            copy(mesh.nnode),
                            copy(mesh.nedge),
                            copy(mesh.nelem),
                            copy(mesh.ndof),
                            copy(mesh.conn),
                            copy(mesh.lorder),
                            copy(mesh.elem_ndof),
                            mesh.elem_type)

function Mesh(m::Int64, n::Int64, h::Float64; order::Int64 = -1,
            degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = -1,
            version::Int64 = 1)
    coords = zeros((m+1)*(n+1), 2)
    elems = zeros(Int64, 2*m*n, 3)
    for i = 1:n
        for j = 1:m
            e = 2*((i-1)*m + j - 1)+1
            if version == 1
                elems[e, :] = [(i-1)*(m+1)+j; (i-1)*(m+1)+j+1; i*(m+1)+j ]
                elems[e+1, :] = [(i-1)*(m+1)+j+1; i*(m+1)+j; i*(m+1)+j+1]
            elseif version == 2
                elems[e, :] = [(i-1)*(m+1)+j; i*(m+1)+j+1; i*(m+1)+j ]
                elems[e+1, :] = [(i-1)*(m+1)+j; (i-1)*(m+1)+j+1; i*(m+1)+j+1]
            elseif version == 3
                if rand()>0.5
                    elems[e, :] = [(i-1)*(m+1)+j; (i-1)*(m+1)+j+1; i*(m+1)+j ]
                    elems[e+1, :] = [(i-1)*(m+1)+j+1; i*(m+1)+j; i*(m+1)+j+1]
                else
                    elems[e, :] = [(i-1)*(m+1)+j; i*(m+1)+j+1; i*(m+1)+j ]
                    elems[e+1, :] = [(i-1)*(m+1)+j; (i-1)*(m+1)+j+1; i*(m+1)+j+1]
                end
            end

        end
    end
    k = 1
    for i = 1:n+1
        for j = 1:m+1
            x = (j-1)*h
            y = (i-1)*h
            coords[k, :] = [x;y]
            k += 1
        end
    end
    Mesh(coords, elems, order, degree, lorder)
end

Mesh_v2(args...;kwargs...) = Mesh(args...; version = 2, kwargs...)
Mesh_v3(args...;kwargs...) = Mesh(args...; version = 3, kwargs...)

"""
    Mesh(; order::Int64 = -1,
    degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = -1)

Creates a mesh with a reference triangle.

![](https://raw.githubusercontent.com/ADCMEMarket/ADCMEImages/master/AdFem/mapping.png)

"""
function Mesh(; order::Int64 = -1,
    degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = -1)
    coords = [
        1.0 0.0
        0.0 1.0
        0.0 0.0
    ]
    elems = [
        1 2 3
    ]
    Mesh(coords, elems, order, degree, lorder)
end

function Base.:length(mesh::Mesh)
    return maximum(mesh.conn)
end

"""
    get_ngauss(mesh::Mesh)

Return the total number of Gauss points.
"""
function get_ngauss(mesh::Mesh)
    return Int64(@eval ccall((:mfem_get_ngauss, $LIBMFEM), Cint, ()))
end

"""
    get_ngauss(mesh::Mesh)

Return the areas of triangles as an array.
"""
function get_area(mesh::Mesh)
    a = zeros(size(mesh.elems,1))
    @eval ccall((:mfem_get_area, $LIBMFEM), Cvoid, (Ptr{Cdouble}, ), $a)
    a
end

"""
    gauss_nodes(mesh::Mesh)
"""
function gauss_nodes(mesh::Mesh)
    ngauss = get_ngauss(mesh)
    x = zeros(ngauss)
    y = zeros(ngauss)
    @eval ccall((:mfem_get_gauss, $LIBMFEM), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}), $x, $y)
    [x y]
end

"""
    fem_nodes(mesh::Mesh)
"""
function fem_nodes(mesh::Mesh)
    mesh.nodes
end

"""
    fvm_nodes(mesh::Mesh)
"""
function fvm_nodes(mesh::Mesh)
    nnode = size(mesh.nodes, 1)
    nelem = size(mesh.elems, 1)
    out = zeros(nelem, 2)
    for i = 1:nelem
        idx = mesh.elems[i, :]
        out[i, 1] = mean(mesh.nodes[idx, 1])
        out[i, 2] = mean(mesh.nodes[idx, 2])
    end
    return out
end

function _edge_dict(mesh::Mesh)
    D = Dict{Tuple{Int64, Int64}, Int64}()
    for i = 1:mesh.nedge
        D[(mesh.edges[i,1], mesh.edges[i,2])] = i
        D[(mesh.edges[i,2], mesh.edges[i,1])] = i
    end
    D
end