text:
#=

Code for solving Laplacian and Diagnally Dominant Systems
This is not all of it, just the short and sweet routines.

Started by Dan Spielman
Contributors: ???

  lapWrapSolver: takes a solver for DD systems, and uses it to solve a lap system in la
  lapWrapSolver(solver, la::AbstractArray)
  lapWrapSolver(solver)
  lapWrapSolver(solver, la::AbstractArray, b) = lapWrapSolver(solver,la)(b)

  For example, to make a Cholesky-based solver for Laplacians, we created
  lapChol = lapWrapSolver(cholfact)

  augmentTree : takes a spanning tree, a graph, and adds in 2k edges of high stretch
    takes both as sparse matrices

  augTreeSolver : a solver that shouldn't suck

=#

include("sampler.jl")
include("randTrees.jl")

#==================================
  Wrapping solver for Laplacians
==================================#

#=
function lapWrapSolver(solver, la::AbstractMatrix)
    N = size(la)[1]
    lasub = la[1:(N-1),1:(N-1)]
    subSolver = solver(lasub);

    f = function(b)
        b = b - mean(b)
        bs = b[1:(N-1)]
        if isa(subSolver,Factorization)
            xs = subSolver \ bs
        else
            xs = subSolver(bs)
        end
        x = [xs;0]
        x = x - mean(x)
        return x
    end

    return f
end
=#

"""Apply the ith solver on the ith component"""
function blockSolver(comps, solvers)

    f = function(b)
        x = zeros(size(b))
        for i in 1:length(comps)
            ind = comps[i]
            bi = b[ind]
            if isa(solvers[i],Factorization)
                x[ind] = solvers[i] \ bi
            else
                x[ind] = (solvers[i])(bi)
            end
        end
        return x
    end

end

"""Takes a solver for solving nonsingular sdd systems,
and returns a solver for solving Laplacian systems.
The optional args tol and maxits are not necessarily taken by
all solvers.  But, if they are, one can pass them here"""
function lapWrapSolver(solver, la::AbstractArray; tol::Real=0.0, maxits::Integer=0, maxtime=Inf)
    N = size(la)[1]

    lasub = la[1:(N-1),1:(N-1)]

    # need to be sure that removing does not disconnect.
    # or, if it does, solve on the components!

    co = components(lasub)
    if maximum(co) == 1

        if tol > 0
            # this section is updated with maxtime, other may be not
            if maxits > 0
                if maxtime > 0
                  subSolver = solver(lasub, tol=tol, maxits=maxits, maxtime=maxtime);
                else
                  subSolver = solver(lasub, tol=tol, maxits=maxits);
                end
            else
                if maxtime > 0
                  subSolver = solver(lasub, tol=tol, maxtime=maxtime);
                else
                  subSolver = solver(lasub, tol=tol);
                end
            end
        else
            if maxits > 0
                subSolver = solver(lasub, maxits=maxits);
            else
                subSolver = solver(lasub);
            end
        end

    else

        comps = vecToComps(co)

        solvers = []
        for i in 1:length(comps)
            ind = comps[i]

            lasubsub = lasub[ind,ind]

            if (length(ind) == 1)
                ssubSolver = x->(x/lasubsub[1,1])

            elseif (length(ind) < 50)
                ssubSolver = cholfact(lasubsub)

            else

                if tol > 0
                    if maxits > 0
                        ssubSolver = solver(lasubsub, tol=tol, maxits=maxits);
                    else
                        ssubSolver = solver(lasubsub, tol=tol);
                    end
                else
                    if maxits > 0
                        ssubSolver = solver(lasubsub, maxits=maxits);
                    else
                        ssubSolver = solver(lasubsub);
                    end
                end

            end
            push!(solvers, ssubSolver)
        end

        subSolver = blockSolver(comps,solvers)

    end

    f = function(b)
        b = b - mean(b)
        bs = b[1:(N-1)]
        if isa(subSolver,Factorization)
            xs = subSolver \ bs
        else
            xs = subSolver(bs)
        end
        x = [xs;0]
        x = x - mean(x)
        return x
    end

    return f
end

#=
function lapWrapSolver(solver)
    f = function(la::AbstractMatrix)
        return lapWrapSolver(solver, la)
    end
    return f
end
=#

function lapWrapSolver(solver; tol::Real=0.0, maxits::Integer=0)
    f = function(la::AbstractArray)
        return lapWrapSolver(solver, la, tol=tol, maxits=maxits)
    end
    return f
end

"""lapChol uses a Cholesky Factorization to solver systems in Laplacians"""
function lapChol()
    return lapWrapSolver(cholfact)
end

#lapWrapSolver(solver, la::AbstractMatrix, b) = lapWrapSolver(solver,la)(b)

lapWrapSolver(solver, la::AbstractArray, b; tol::Real=0.0, maxits::Integer=0) = lapWrapSolver(solver,la, tol=tol, maxits=maxits)(b)

#=========================================
  Augmented Spanning Tree Preconditioner
=========================================#

"""Takes as input a tree and an adjacency matrix of a graph.
It then computes the stretch of every edge of the graph wrt
the tree.  It then adds back the k edges of highest stretch,
and k edges sampled according to stretch"""
function augmentTree{Tv,Ti}(tree::SparseMatrixCSC{Tv,Ti}, mat::SparseMatrixCSC{Tv,Ti}, k::Ti)

    st = compStretches(tree, mat)

    # just to be safe, remove the tree from this
    #=
    (ti,tj) = findnz(tree)
    for i in 1:length(ti)   # most of the time
        st[ti[i],tj[i]] = 0
        st[tj[i],ti[i]] = 0
    end
    =#

    (ai,aj,av) = findnz(triu(st))

    ord = sortperm(av, rev=true)

    edgeinds = zeros(Bool,length(av))
    for i in 1:k
        edgeinds[ord[i]] = true
    end

    s = sum(av[ord[(k+1):end]])
    probs = av * k / s
    probs[ord[1:k]] = 0
    edgeinds[rand(length(av)) .< probs] = true

    augi = ai[edgeinds]
    augj = aj[edgeinds]
    augm = length(augi)
    augv = zeros(Tv, augm)
    for i in 1:augm,
        augv = mat[augi[i],augj[i]]
    end

    n = size(mat)[1]
    aug = sparse(augi, augj, augv, n, n)
    aug = aug + aug'

    return tree + aug

end

"""This is an augmented spanning tree preconditioner for diagonally dominant
linear systems.  It takes as optional input a tree growing algorithm.
It adds back 2sqrt(n) edges via augmentTree: the sqrt(n) of highest stretch
and another sqrt(n) sampled according to stretch.
For most purposes, one should directly call `augTreeSolver`."""
function augTreePrecon{Tv,Ti}(ddmat::SparseMatrixCSC{Tv,Ti}; treeAlg=akpw)

  adjmat = -triu(ddmat,1)
  adjmat = adjmat + adjmat'

  tree = treeAlg(adjmat)

  n = size(ddmat)[1]

  augtree = augmentTree(tree,adjmat,convert(Int,round(sqrt(n))))

  Dx = spdiagm(ddmat*ones(n))

  augDD = Dx + spdiagm(augtree*ones(n)) - augtree

  F = cholfact(augDD)

  return x -> (F\x)

end

"""
An "augmented spanning tree" solver for positive definite diagonally dominant matrices.
It works by adding edges to a low stretch spanning tree.  It calls `augTreePrecon` to form
the preconditioner.

~~~julia
 augTreeSolver{Tv,Ti}(ddmat::SparseMatrixCSC{Tv,Ti}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, treeAlg=akpw)
~~~
"""
function augTreeSolver{Tv,Ti}(ddmat::SparseMatrixCSC{Tv,Ti}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, treeAlg=akpw)

  F = augTreePrecon(ddmat, treeAlg=treeAlg)

  f(b;maxits=maxits, maxtime=maxtime, verbose=verbose) = pcg(ddmat, b, F, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)

  return f

end

"""This is an augmented spanning tree preconditioner for Laplacians.
It takes as optional input a tree growing algorithm.
It adds back 2sqrt(n) edges via `augmentTree`: the sqrt(n) of highest stretch
and another sqrt(n) sampled according to stretch.
For most purposes, one should directly call `augTreeLapSolver`."""
function augTreeLapPrecon{Tv,Ti}(ddmat::SparseMatrixCSC{Tv,Ti}; treeAlg=akpw)

  adjmat = -triu(ddmat,1)
  adjmat = adjmat + adjmat'

  tree = treeAlg(adjmat)

  n = size(ddmat)[1]

  augtree = augmentTree(tree,adjmat,convert(Int,round(sqrt(n))))

  Dx = spdiagm(ddmat*ones(n))

  augDD = Dx + spdiagm(augtree*ones(n)) - augtree

  F = lapWrapSolver(cholfact,augDD)

  return F

end

"""
An "augmented spanning tree" solver for Laplacian matrices.
It works by adding edges to a low stretch spanning tree.  It calls `augTreeLapPrecon` to form
the preconditioner. In line with other solver, it takes as input the adjacency matrix of the system.

~~~julia
 augTreeLapSolver{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, treeAlg=akpw)
~~~
"""
function augTreeLapSolver{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, treeAlg=akpw)

  la = lap(a)

  F = augTreeLapPrecon(la, treeAlg=treeAlg)

  f(b;maxits=maxits, maxtime=maxtime, verbose=verbose) = pcg(la, b, F, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)

  return f

end

"""
A wrapper for the PyAMG solver.

~~~julia
 amgSolver{Tv,Ti}(ddmat::SparseMatrixCSC{Tv,Ti}; tol::Float64=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
~~~
"""
function AMGSolver{Tv,Ti}(ddmat::SparseMatrixCSC{Tv,Ti}; tol::Float64=1e-6, maxits=Inf, maxtime=Inf, verbose=false)

  amg = PyAMG.RugeStubenSolver(ddmat);
  M = PyAMG.aspreconditioner(amg);
  function F(b)
    return M \ b;
  end

  f(b;maxits=maxits, maxtime=maxtime, verbose=verbose) = pcg(ddmat, b, F, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)

  return f

end

"""
A wrapper for the PyAMG solver. In line with our other solvers, takes in an adjacency matrix.

~~~julia
 amgSolver{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti}; tol::Float64=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
~~~
"""
function AMGLapSolver{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti}; tol::Float64=1e-6, maxits=Inf, maxtime=Inf, verbose=false)

  la = lap(a)

  amg = PyAMG.RugeStubenSolver(la);
  M = PyAMG.aspreconditioner(amg);
  function F(b)
    return M \ b;
  end

  f(b;maxits=maxits, maxtime=maxtime, verbose=verbose) = pcg(la, b, F, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)

  return f

end


text:
using Documenter
using Optinum

makedocs(
    modules = [Optinum],
    sitename = "Optinum.jl",
    authors = "Saloua Naama, Mohamed El Waghf et Rachid ELMontassir",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
            "Accueil" => "index.md",
	    "Sujet" => "Sujet.md",
            "Algorithmes" => [
            "L'algorithme de Newton local" => "Algorithme_de_newton.md",
            "La méthode des régions de confiance" => "Regions_de_confiance.md",
            "La méthode du Lagrangien augmenté" => "Lagrangien_augmente.md"
            ],
	    "Index des fonctions" =>"fct_index.md",
            "Annexes" => "Annexes.md",
            #"Exemples d'appels" =>"Exemples.md",
            "Installation de Julia et tests unitaires" => "mise_en_place.md",
            #"Julia vs MatLab"=> "julia_vs_matlab.md",
            #"Création de Modules en Julia" => "create_package.md",
            #"Création d'un dépôt Git pour un module Julia" => "git_doc.md",
            #"Géneration de la doc" => "generation_doc.md",
            #"Intégration continue et déploiement de la Doc avec Travis" =>"Integration_continue.md",
	    #"Précompilation des modules"=>"Precompilation.md",
            #"Génération du rapport" => "generate_pdf.md",
            "Foire aux Questions" =>"FAQ.md"

            ]
    )

deploydocs(
    repo   = "github.com/mathn7/Optinum.git",
    branch = "gh-pages",
    devbranch = "master",
)

ext: jl

max_stars_repo_name: faicaltoubali/ENSEEIHT

max_issues_repo_name: faicaltoubali/ENSEEIHT

max_forks_repo_name: faicaltoubali/ENSEEIHT


text:
abstract type DOFNumbering end

"""
Creates and stores the DOF numbering associated to the approximation space
constructed using the basis `b` in the `Domain` `Ω`.

In practice, the vector `dofs` contains all unique *global* identifiers (type
`T`) of degrees of freedom present in the domain `Ω` in reference to the
global information contained in `mesh`. The length of this vector depends on
the domain `Ω` but not the values of its elements. The *local* numbering of
the DOFs are given by the indices in this vector `dofs`, so in a sense this
gives a local to global numbering mapping.

The vector `elt2dof` contains, for each element (on which the basis `b` is
defined) in the domain `Ω`, the indices in the vector `dofs` of all degrees of
freedom that are (partly) supported on the said element.

Some optimization is possible in some cases where we could avoid the
allocation of one or both attributes `dofs` and `elt2dof` altogether.

Warning: this is not robust if the underlying `mesh` is modified.
"""
struct LocalDOFNumbering{S,T} <: DOFNumbering
    mesh::AbstractMesh
    Ω::Domain
    b::PolynomialBasis
    dofs::Vector{S}
    elt2dof::Vector{T}
end

dofs(dofnb::LocalDOFNumbering) = dofnb.dofs
elt2dof(dofnb::LocalDOFNumbering) = dofnb.elt2dof

function _element_support_basis(mesh, b::LagrangeBasis{D,p}) where {D,p}
    elts = DataType[]
    for E in etypes(mesh)
        if E <: AbstractElement
            if typeof(domain(E)) == D
                push!(elts, E)
            end
        end
    end
    @assert length(elts) == 1
    return elts[1]
end

"""
One degree of freedom per mesh element in the domain `Ω`.
"""
function local_dof_numbering(mesh::GenericMesh{d,T}, Ω::Domain,
                             b::LagrangeBasis{D,0}) where {d,D,Np,T}
    E = _element_support_basis(mesh, b)
    dofs = Int64[]
    for ent in entities(Ω)
        append!(dofs, ent2tags(mesh)[ent][E])
    end
    elt2dof = collect(1:length(dofs))
    return LocalDOFNumbering(mesh, Ω, b, dofs, elt2dof)
end

"""
One degree of freedom per mesh node in the domain `Ω`.
"""
function local_dof_numbering(mesh::GenericMesh, Ω::Domain,
                             b::LagrangeBasis{D,1}) where {D}
    E = _element_support_basis(mesh, b)
    N = number_of_nodes(D)
    node_tags = Int64[]
    for ent in entities(Ω)
        ent_tags = ent2tags(mesh)[ent][E]
        node_tags_ent = unique!(sort!(elements(mesh)[E][:,ent_tags][:]))
        append!(node_tags, node_tags_ent)
    end
    dofs = unique!(sort!(node_tags))
    glob2loc = SparseVector(length(nodes(mesh)),dofs,collect(1:length(dofs)))
    elt2dof = SVector{N,Int64}[]
    for ent in entities(Ω)
        ent_tags = ent2tags(mesh)[ent][E]
        node_tags_ent = elements(mesh)[E][:,ent_tags]
        for ie in 1:size(node_tags_ent,2)
            push!(elt2dof, glob2loc[node_tags_ent[:,ie]])
        end
    end
    return LocalDOFNumbering(mesh, Ω, b, dofs, elt2dof)
end

function assembly(mesh, Ω, u, v; order=1,
                  f=(u,v)->(i,j,x̂,_,_)->u(x̂)[j]*v(x̂)[i])
    submesh = view(mesh, Ω)
    etype2qrule = Mesh._qrule_for_mesh(submesh, order)
    E = _element_support_basis(submesh, u)
    @assert E == _element_support_basis(submesh, v)
    qrule = etype2qrule[E]
    udofnb = local_dof_numbering(mesh, Ω, u)
    vdofnb = local_dof_numbering(mesh, Ω, v)
    n = length(dofs(udofnb))
    m = length(dofs(vdofnb))
    Is = Int64[]
    Js = Int64[]
    Vs = Float64[]
    M  = zeros(Float64,length(v),length(u))
    for (iel, el) in enumerate(elements(submesh, E))
        utags = elt2dof(udofnb)[iel]
        vtags = elt2dof(vdofnb)[iel]
        elementary_matrix!(M,el, u, v, qrule; f=f)
        for (iloc, iglob) in enumerate(vtags)
            for (jloc, jglob) in enumerate(utags)
                push!(Is, iglob)
                push!(Js, jglob)
                push!(Vs, M[iloc,jloc])
            end
        end
        fill!(M,0)
    end
    A = sparse(Is,Js,Vs,m,n)
    return A
end

function assembly(mesh, Ω, v; order=1,
                  f=(v)->(i,x̂,_,_)->v(x̂)[i])
    submesh = view(mesh, Ω)
    etype2qrule = Mesh._qrule_for_mesh(submesh, order)
    E = _element_support_basis(submesh, v)
    qrule = etype2qrule[E]
    vdofnb = local_dof_numbering(mesh, Ω, v)
    m = length(dofs(vdofnb))
    b = zeros(Float64, m)
    for (iel, el) in enumerate(elements(submesh, E))
        vtags = elt2dof(vdofnb)[iel]
        M = elementary_matrix(el, v, qrule; f=f)
        for (iloc, iglob) in enumerate(vtags)
            b[iglob] += M[iloc]
        end
    end
    return b
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using DiffEqPhysics, ForwardDiff, OrdinaryDiffEq
using StaticArrays, LinearAlgebra
using SafeTestsets, Test

@safetestset "Hamiltonian Test" begin include("hamiltonian_test.jl") end
#@safetestset "N-Body Test" begin include("nbody_test.jl") end


text:
"""
"""
module GPU # begin submodule PredictMD.GPU

############################################################################
# PredictMD.GPU source files ###############################################
############################################################################

end # end submodule PredictMD.GPU


text:
import Zygote

using Test: @test, @testset, @inferred
using Random: bitrand
using Statistics: mean
using LinearAlgebra: I
using RestrictedBoltzmannMachines: RBM, BinaryRBM
using RestrictedBoltzmannMachines: visible, hidden, weights
using RestrictedBoltzmannMachines: energy, interaction_energy, free_energy, ∂free_energy
using RestrictedBoltzmannMachines: inputs_v_to_h, inputs_h_to_v
using WhitenedRBMs: BinaryWhiteRBM, WhiteRBM, Affine
using WhitenedRBMs: whiten, blacken, whiten_visible, whiten_hidden
using WhitenedRBMs: whiten!, whiten_visible!, whiten_hidden!
using WhitenedRBMs: energy_shift, energy_shift_visible, energy_shift_hidden

@testset "whiten / blacken" begin
    rbm = @inferred BinaryRBM(randn(3), randn(2), randn(3,2))
    affine_v = @inferred Affine(randn(3,3), randn(3))
    affine_h = @inferred Affine(randn(2,2), randn(2))
    white_rbm = @inferred whiten(rbm, affine_v, affine_h)
    @test white_rbm.affine_v.u ≈ affine_v.u
    @test white_rbm.affine_h.u ≈ affine_h.u
    @test white_rbm.affine_v.A ≈ affine_v.A
    @test white_rbm.affine_h.A ≈ affine_h.A
    @test affine_v.A' * weights(white_rbm) * affine_h.A ≈ weights(rbm)
    @test visible(rbm).θ ≈ visible(white_rbm).θ - weights(rbm) * affine_h.u
    @test hidden(rbm).θ ≈ hidden(white_rbm).θ - weights(rbm)' * affine_v.u
    @test blacken(rbm) == rbm

    brbm = @inferred blacken(white_rbm)
    wrbm = @inferred whiten(white_rbm)
    @test visible(brbm).θ ≈ visible(wrbm).θ ≈ visible(rbm).θ
    @test hidden(brbm).θ ≈ hidden(wrbm).θ ≈ hidden(rbm).θ
    @test weights(brbm) ≈ weights(wrbm) ≈ weights(rbm)

    @test iszero(wrbm.affine_v.u)
    @test iszero(wrbm.affine_h.u)
    @test wrbm.affine_v.A ≈ I
    @test wrbm.affine_h.A ≈ I

    # whiten_visible
    affine_v_new = @inferred Affine(randn(3,3), randn(3))
    white_rbm_new = @inferred whiten_visible(white_rbm, affine_v_new)
    @test white_rbm_new.affine_v.u == affine_v_new.u
    @test white_rbm_new.affine_v.A == affine_v_new.A
    @test white_rbm_new.affine_h.u == affine_h.u
    @test white_rbm_new.affine_h.A == affine_h.A
    white_rbm_expected = whiten(white_rbm, affine_v_new, affine_h)
    @test visible(white_rbm_new).θ ≈ visible(white_rbm_expected).θ
    @test hidden(white_rbm_new).θ ≈ hidden(white_rbm_expected).θ
    @test weights(white_rbm_new) ≈ weights(white_rbm_expected)

    # whiten_hidden
    affine_h_new = @inferred Affine(randn(2,2), randn(2))
    white_rbm_new = @inferred whiten_hidden(white_rbm, affine_h_new)
    @test white_rbm_new.affine_v.u == affine_v.u
    @test white_rbm_new.affine_v.A == affine_v.A
    @test white_rbm_new.affine_h.u == affine_h_new.u
    @test white_rbm_new.affine_h.A == affine_h_new.A
    white_rbm_expected = whiten(white_rbm, affine_v, affine_h_new)
    @test visible(white_rbm_new).θ ≈ visible(white_rbm_expected).θ
    @test hidden(white_rbm_new).θ ≈ hidden(white_rbm_expected).θ
    @test weights(white_rbm_new) ≈ weights(white_rbm_expected)

    wrbm1 = whiten_visible(whiten_hidden(rbm, affine_h), affine_v)
    wrbm2 = whiten_hidden(whiten_visible(rbm, affine_v), affine_h)
    for wrbm in (wrbm1, wrbm2)
        @test wrbm.affine_v.u == affine_v.u
        @test wrbm.affine_v.A == affine_v.A
        @test wrbm.affine_h.u == affine_h.u
        @test wrbm.affine_h.A == affine_h.A
        @test visible(wrbm).θ ≈ visible(white_rbm).θ
        @test hidden(wrbm).θ ≈ hidden(white_rbm).θ
        @test weights(wrbm) ≈ weights(white_rbm)
    end
end

@testset "energy invariance" begin
    affine_v = @inferred Affine(randn(3,3), randn(3))
    affine_h = @inferred Affine(randn(2,2), randn(2))
    rbm = BinaryRBM(randn(3), randn(2), randn(3,2))
    white_rbm = @inferred whiten(rbm, affine_v, affine_h)
    v = bitrand(size(visible(rbm))..., 100)
    h = bitrand(size(hidden(rbm))..., 100)

    @test energy(visible(rbm), v) ≈ energy(visible(white_rbm), v) + v' * weights(rbm) * affine_h.u
    @test energy(hidden(rbm), h) ≈ energy(hidden(white_rbm), h) + h' * weights(rbm)' * affine_v.u

    @test inputs_v_to_h(rbm, v .- affine_v.u) ≈ @inferred inputs_v_to_h(white_rbm, v)
    @test inputs_h_to_v(rbm, h .- affine_h.u) ≈ @inferred inputs_h_to_v(white_rbm, h)

    ΔE = interaction_energy(rbm, affine_v.u, affine_h.u)::Real
    @test @inferred(energy_shift(rbm, affine_v, affine_h)) ≈ ΔE
    @test energy(rbm, v, h) ≈ energy(white_rbm, v, h) .- ΔE
    @test free_energy(rbm, v) ≈ free_energy(white_rbm, v) .- ΔE

    @test energy_shift_visible(rbm, affine_v) ≈ energy_shift(rbm, affine_v, one(affine_h))
    @test energy_shift_hidden(rbm, affine_h)  ≈ energy_shift(rbm, one(affine_v), affine_h)
end

@testset "whiten!" begin
    rbm = @inferred BinaryWhiteRBM(
        randn(3), randn(2), randn(3,2),
        Affine(randn(3,3), randn(3)),
        Affine(randn(2,2), randn(2))
    )
    @test iszero(energy_shift(rbm, rbm.affine_v, rbm.affine_h))

    v = bitrand(size(visible(rbm))..., 100)
    h = bitrand(size(hidden(rbm))..., 100)
    E = energy(rbm, v, h)
    F = free_energy(rbm, v)
    affine_v = @inferred Affine(randn(3,3), randn(3))
    affine_h = @inferred Affine(randn(2,2), randn(2))

    @test energy_shift_visible(rbm, affine_v) ≈ energy_shift(rbm, affine_v, rbm.affine_h)
    @test energy_shift_hidden(rbm, affine_h) ≈ energy_shift(rbm, rbm.affine_v, affine_h)

    ΔE = energy_shift(rbm, affine_v, affine_h)
    @test energy(whiten(rbm, affine_v, affine_h), v, h) ≈ E .+ ΔE

    whiten!(rbm, affine_v, affine_h)
    @test rbm.affine_v.u ≈ affine_v.u
    @test rbm.affine_h.u ≈ affine_h.u
    @test rbm.affine_v.A ≈ affine_v.A
    @test rbm.affine_h.A ≈ affine_h.A
    @test energy(rbm, v, h) ≈ E .+ ΔE
    @test free_energy(rbm, v) ≈ F .+ ΔE
end

@testset "free energy" begin
    affine_v = Affine(randn(3,3), randn(3))
    affine_h = Affine(randn(2,2), randn(2))
    white_rbm = BinaryWhiteRBM(randn(3), randn(2), randn(3,2), affine_v, affine_h)
    v = bitrand(size(visible(white_rbm))...)
    F = -log(sum(exp(-energy(white_rbm, v, h)) for h in [[0,0], [0,1], [1,0], [1,1]]))
    @test free_energy(white_rbm, v) ≈ F
end

@testset "∂free energy" begin
    affine_v = Affine(randn(3,3), randn(3))
    affine_h = Affine(randn(2,2), randn(2))
    white_rbm = BinaryWhiteRBM(randn(3), randn(2), randn(3,2), affine_v, affine_h)
    v = bitrand(size(visible(white_rbm))...)
    gs = Zygote.gradient(white_rbm) do white_rbm
        mean(free_energy(white_rbm, v))
    end
    ∂ = ∂free_energy(white_rbm, v)
    @test ∂.visible.θ ≈ only(gs).visible.θ
    @test ∂.hidden.θ ≈ only(gs).hidden.θ
    @test ∂.w ≈ only(gs).w
end


text:
using DataStructures: CircularBuffer
include("AlgorithmicEnv.jl")

"""
Task is to reverse content over the input tape.
http://arxiv.org/abs/1511.07275
"""
mutable struct ReverseEnv <: TapeAlgorithmicEnv
    MIN_REWARD_SHORTFALL_FOR_PROMOTION::Float32
    base::Int
    episode_total_reward::Float32
    reward_shortfalls::CircularBuffer{Float32}
    charmap::Array{Char, 1}
    min_length::Int
    action_space::TupleSpace
    observation_space::Discrete

    MOVEMENTS::NTuple
    READ_HEAD_START::Int

    target::Array{Int8, 1}  # target tape
    input_data::Array{Int8, 1}
    time::Int8
    read_head_position::Int8  # current position of the read head.
    write_head_position::Int8
    last_action
    last_reward::Float32
end

ReverseEnv(base::Int=2) =
    ReverseEnv(-1f-1, # MIN_REWARD_SHORTFALL_FOR_PROMOTION
            base, # last
            0f0, # episode_total_reward
            CircularBuffer{Float32}(50),  # reward_shortfalls
            push!(['A'+i for i=0:base-1], ' '),  # charmap
            1,  # starting min_length
            TupleSpace([Discrete(2), Discrete(2), Discrete(base)]),  # action_space
            Discrete(base+1),  # observation_space
            (:left, :right), 1,  # MOVEMENTS and READ_HEAD_START
            Int8[], Int8[], 0, 1, 1, nothing, 0.0)

target_from_input_data!(env::ReverseEnv) =
    env.target = [char for char in reverse(env.input_data)]


text:
using RegularizedRegression
using BaseTestNext
using Compat

const NumTestRepetitions = 30

using RegularizedRegression: learn, coefficients, RegressionResult

function make_dataset(pRange, nRange, coefRange)
    N = rand(pRange)
    P = rand(nRange)
    X = randn(P, N)
    sel = rand(1:P)
    coef = rand(coefRange)
    y = X[sel, :] * coef
    ds = MatrixDataSet(X, y)
    return N, P, sel, X, y, ds
end

function test_learner_with_includeIntercept(learnertype, maxMape = 10.0)
    N, P, sel, X, y, ds = make_dataset(10:50, 2:8, 1.0:0.01:10.0)
    l = learnertype()

    res = learn(l, ds)
    @test (typeof(res) <: Dict || typeof(res) <: RegressionResult)

    beta = coefficients(l, ds)
    @test typeof(beta) <: Array
    @test length(beta) == P+1

    beta2 = coefficients(l, ds; includeIntercept = false)
    @test typeof(beta2) <: Array
    @test length(beta2) == P

    yhat = X' * beta2
    @test RegularizedRegression.mape(yhat, y) < maxMape # Accuracy should be very high since no noise
end

function test_learner(learnertype, maxMape = 10.0)
    N, P, sel, X, y, ds = make_dataset(10:50, 2:8, 1.0:0.01:10.0)
    l = learnertype()

    res = learn(l, ds)
    @test (typeof(res) <: Dict || typeof(res) <: RegressionResult)

    beta = coefficients(l, ds)
    @test typeof(beta) <: Array
    @test length(beta) == P

    yhat = X' * beta
    @test RegularizedRegression.mape(yhat, y) < maxMape # Accuracy should be very high since no noise
end


text:
function runfig7()
    p = Params(Ln = 195, W = 2500, Ls = 20, Ws = 20, scale = 40, λ = 5,
    Δ = 1, d = 0, τ = 1);
    ϕlist = -2π:π/10:2π
    println("1")
    _, esa = spectrumvsphase(ϕlist,
        reconstruct(p, EZ = SA[0, 2, 0], μN = 0.6, α = 0), nev = 16)
    println("2")
    _, esb = spectrumvsphase(ϕlist,
        reconstruct(p, EZ = SA[0, 4, 0], μN = 1.192, α = 0), nev = 16)
    println("3")
    _, esc = spectrumvsphase(ϕlist,
        reconstruct(p, EZ = SA[0, 2, 0], μN = 0.6, α = 2), nev = 16)
    println("4")
    ϕs, esd = spectrumvsphase(ϕlist,
        reconstruct(p, EZ = SA[0, 4, 0], μN = 1.192, α = 2), nev = 16)
    savespectrumvsphase("fig7", p, ϕs, esa, esb, esc, esd)
end

# Extra: checking W and Ln dependence

# function Wtest(Wlist, p, ϕ0, panel; nev = 4, kw...)
#     elist = zeros(Float64, nev, length(Wlist))
#     for i in 1:length(Wlist)
#         p = reconstruct(p, W = Wlist[i])
#         ph = rectangle_weaklink(p; kw...)
#         s = spectrum(ph(phi = ϕ0), method = ArpackPackage(sigma = 1e-7, nev = nev))
#         elist[:,i] = s.energies

#         println(s.energies, " W: ",Wlist[i])
#     end
#     return Wlist, elist
# end

# function testW(Ln0, ϕ0, panel)
#     Wlist = 2000:200:2500
#     p = Params(Ln = Ln0, Ls = 20, Ws = 20, scale = 40, λ = 5, Δ = 1, d = 0, τ = 1);
#     if panel == "a"
#         p = reconstruct(p, EZ = SA[0, 2, 0], μN = 0.6, α = 0)
#     elseif panel == "b"
#         p =  reconstruct(p, EZ = SA[0, 4, 0], μN = 1.192, α = 0)
#     elseif panel == "c"
#         p = reconstruct(p, EZ = SA[0, 2, 0], μN = 0.6, α = 2)
#     else
#         p =  reconstruct(p, EZ = SA[0, 4, 0], μN = 1.192, α = 2)
#     end

#     Wtest(Wlist, p, ϕ0, panel)
# end

# function Lntest(Lnlist, p, ϕ0, panel; nev = 4, kw...)
#     elist = zeros(Float64, nev, length(Lnlist))
#     for i in 1:length(Lnlist)
#         p = reconstruct(p, Ln = Lnlist[i])
#         ph = rectangle_weaklink(p; kw...)
#         s = spectrum(ph(phi = ϕ0), method = ArpackPackage(sigma = 1e-7, nev = nev))
#         elist[:,i] = s.energies

#         println(s.energies, " Ln: ",Lnlist[i])
#     end
#     return Lnlist, elist
# end

# function testLn(W0, ϕ0, panel, Lnlist )

#     p = Params(W = W0, Ls = 20, Ws = 20, scale = 40, λ = 5, Δ = 1, d = 0, τ = 1);
#     if panel == "a"
#         p = reconstruct(p, EZ = SA[0, 2, 0], μN = 0.6, α = 0)
#     elseif panel == "b"
#         p =  reconstruct(p, EZ = SA[0, 4, 0], μN = 1.192, α = 0)
#     elseif panel == "c"
#         p = reconstruct(p, EZ = SA[0, 2, 0], μN = 0.6, α = 2)
#     else
#         p =  reconstruct(p, EZ = SA[0, 4, 0], μN = 1.192, α = 2)
#     end

#     Lntest(Lnlist, p, ϕ0, panel)
# end

max_stars_repo_name: fernandopenaranda/MBSinBLG

max_issues_repo_name: fernandopenaranda/MBSinBLG

max_forks_repo_name: fernandopenaranda/MBSinBLG


text:
function compute_thrust_direction(obj::ThrustDirectionProvider, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "computeThrustDirection", Vector3D, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end


text:
module MsgPackRpcClient

include("const.jl")
include("base.jl")
include("socks.jl")
include("session.jl")
include("future.jl")

using MsgPack

export MsgPackRpcClientSession, call, get

function call(s::MsgPackRpcClientSession.Session, method::String, params...; sync::Bool = DEFAULT_SYNC_POLICY, sock = nothing)
  if sock == nothing
    try
      sock = s.get_sock()
    catch
      MsgPackRpcClientSocks.connect_and_push!(s.socks)
      s.ptr = 1
      sock = s.get_sock()
    end
  end

  # NOTE: Get msg_id and calculate next_id never to use the same msg_id even if a request fails.
  msg_id = s.next_id
  # NOTE: For compatibility, msgid must be Int32
  s.next_id = s.next_id >= 1<<(BIT_OF_MSGID - 1) ? 0 : s.next_id + 1

  future = send_request(sock, msg_id, method, params)

  if sync == true
    receive_response(sock, future)
    return get(future)
  else
    future.task = @async receive_response(sock, future)
    return future
  end
end

function send_request(sock, msg_id::Int, method::String, args)
  params = MsgPackRpcClientBase.to_a(args)
  packed_data = MsgPack.pack({REQUEST, msg_id, method, params})
  if typeof(sock) == Base.TcpSocket
    send_data_in_tcp(sock, packed_data)
  #else if typeof(sock) == Base.UdpSocket
  #  send_data_in_udp(sock, packed_data, host, port)
  #else
  #  # error
  end
  MsgPackRpcClientFuture.Future(TIMEOUT_IN_SEC, nothing, nothing, nothing, false, nothing, nothing, nothing, msg_id, nothing)
end

function send_data_in_tcp(sock::Base.TcpSocket, data)
  write(sock, data)
  nothing
end

function send_data_in_udp(sock::Base.UdpSocket, data, host::String, port::Int)
  send(sock, host, port, data)
  nothing
end

function get(future::MsgPackRpcClientFuture.Future)
  if future.task != nothing
    wait(future.task)
  end
  if future.error == nothing
    # if future.result_handler != nothing
    #   return result_handler.call()
    # end
    return future.result
  end
  if future.result == nothing
    # TODO: raise instead of return
    return nothing
  end
  # if future.error_handler != nothing
  #   return error_handler.call()
  # end
  # TODO: raise instead of return
  nothing
end

function receive_response(sock, future::MsgPackRpcClientFuture.Future; interval = DEFAULT_INTERVAL)
  unpacked_data = {}
  timeout = future.timeout

  while 0 <= timeout
    receive_data(sock, future)
    # if future.is_set == false
    #  println("continue")
    #  sleep(1)
    #  continue
    # end
    unpacked_data = MsgPack.unpack(future.raw)
    if unpacked_data[1] != RESPONSE || unpacked_data[2] != future.msg_id # type, msgid
      timeout -= interval
      sleep(interval)
      continue
    else
      break
    end
    return nothing
  end

  if unpacked_data[3] != nothing # error
    future.error = unpacked_data[3]
  end

  if unpacked_data[4] != nothing #result
    future.result = unpacked_data[4]
  end

  future
end

function receive_data(sock, future::MsgPackRpcClientFuture.Future)
  join(sock, future)
  nothing
end

function join(sock, future::MsgPackRpcClientFuture.Future; interval = DEFAULT_INTERVAL)
  timeout = future.timeout
  while future.is_set == false && timeout > 0
    future.raw = readavailable(sock)
    if length(future.raw) > 0
      future.is_set = true
      break
    end
    sleep(interval)
    timeout -= interval
  end
  future
end

end # module MsgPackRpcClient

max_stars_repo_name: keithseahus/msgpack-rpc-julia

max_issues_repo_name: keithseahus/msgpack-rpc-julia

max_forks_repo_name: keithseahus/msgpack-rpc-julia


text:
using SafeTestsets

@safetestset "Unit tests" begin
    include("unittest.jl")
end
@safetestset "End to end test" begin
    include("full.jl")
end


text:
# Tensor product basis
# ------------------------------------------------------------------------------------------
struct TensorProductBasis{B1 <: BasisSet, B2 <: BasisSet} <: BasisSet
    b1 :: B1
    b2 :: B2
end
Base.length(b::TensorProductBasis) = length(b.b1) * length(b.b2)
leftindex(b::TensorProductBasis, idx::Integer) = rem(idx - 1, length(b.b1)) + 1
rightindex(b::TensorProductBasis, idx::Integer) = div(idx - 1, length(b.b1)) + 1

struct ProductExtensionOperator{B1 <: BasisSet, B2 <: BasisSet, O1 <: LinearOperator{B1}, O2 <: LinearOperator{B2}} <: LinearOperator{TensorProductBasis{B1, B2}}
    tpb :: TensorProductBasis{B1, B2}
    op1 :: O1
    op2 :: O2
end
basis(op::ProductExtensionOperator) = op.tpb
function apply(i::Integer, op::ProductExtensionOperator, j::Integer)
    i1, i2 = leftindex(op.tpb, i), rightindex(op.tpb, i)
    j1, j2 = leftindex(op.tpb, j), rightindex(op.tpb, j)
    apply(i1, op.op1, j1) * apply(i2, op.op2, j2)
end

# Tensor sum basis
# ------------------------------------------------------------------------------------------
struct TensorSumBasis{B1 <: BasisSet, B2 <: BasisSet} <: BasisSet
    b1 :: B1
    b2 :: B2
end
Base.length(b::TensorSumBasis) = length(b.b1) + length(b.b2)

struct SumExtensionOperator{B1 <: BasisSet,B2 <: BasisSet,O1 <: LinearOperator{B1}, O2 <: LinearOperator{B2}} <: LinearOperator{TensorSumBasis{B1, B2}}
    tsb :: TensorSumBasis{B1, B2}
    op1 :: O1
    op2 :: O2
end
basis(op::SumExtensionOperator) = op.tsb
function apply(i::Integer, op::SumExtensionOperator, j::Integer) :: ComplexF64
    N1 = length(op.tsb.b1)
    if i <= N1 && j <= N1
        apply(i, op.op1, j)
    elseif i > N1 && j > N1
        apply(i - N1, op.op2, j - N1)
    else
        0.0
    end
end


text:
# # Julia Crash Course

# Hello! This is a quick intro to programming in Julia to help you hit the ground running with the _12 Steps to Navier–Stokes_.

# There are two ways to enjoy these lessons with Julia:

# 1. You can download and install a Julia distribution on your computer. Our recommendation is via [julialang.org](julialang.org). Open them up in your favorite IDE and run them one at a time.

# 2. You can run  Julia in a Pluto notebook with `using Pluto; Pluto.run()`.

# In either case, you will probably want to download a copy of this notebook, or the whole AeroPython (?) collection. We recommend that you then follow along each lesson, experimenting with the code in the notebooks, or typing the code into a separate Julia interactive session.

#If you decided to work on your local Julia installation, you will have to navigate in the terminal to the folder that contains the `.jl` files. Then, to launch the notebook server, just type:
#`julia> using Pluto; Pluto.run()`

#You will get a new browser window or tab with a list of the notebooks available in that folder. Click on one and start working!

## Libraries

#Julia is a high-level open-source language.  But the _Julia world_ is inhabited by many packages or libraries that provide useful things like gpu kernels, plotting functions, and much more. We can import libraries of functions to expand the capabilities of Python in our programs.

#OK! We'll start by importing a few libraries to help us out. First: our favorite library is **NumPy**, ... Just kidding! We don't need NumPy. We will need a 2D plotting library which we will to plot our results.
#The following code will be at the top of most of your programs, so execute this cell first:

#```julia
## <-- comments in Julia are denoted by the pound sign, like this one

#using LinearAlgebra                # we import the array library
#using Plots    # import plotting library
#```

#We are importing one library named `LinearAlgebra` and we are importing a package called `Plots`.
#To use a function belonging to one of these libraries, we have to tell Julia where to look for it. For that, each function name is written following the library name, with a dot in between.
#So if we want to use the unexported functions in `LinearAlgebra`, we have to put a name in front of it:
#
#

myarray =  [ 1 2; 3 4 ]
myarray
LinearAlgebra.norm(myarray)
5.47722

# In this tutorial, it's not too important to know where functions are defined since `using Pkg` calls them all into the global namespace to be used, but if
# you were developing your own code to share with others some additional details about namespacing could save you some headaches.

## Variables

# Julia doesn't require explicitly declared variable types like C and other languages.

a = 5        #a is an integer 5
b = "five"   #b is a string of the word "five'
c = 5.0      #c is a floating point 5

type(a)
type(b)
type(c)

# Note that if you divide an integer by an integer that yields a float. If you want to use integer division use the `div` operator.

## Whitespace in Julia
#
# Doesn't exist!
#
## Slicing Arrays

#In Julia, you can look at portions of arrays in the same way as in `Matlab`, with a few extra tricks thrown in.  Let's take an array of values from 1 to 5.

myvals = [1, 2, 3, 4, 5]

#Julia uses a **one-based index**, so let's look at the first and last element in the array `myvals`

myvals[1], myvals[5]

#Note here, the slice is inclusive on the front end the back end.

## Assigning Array Variables

# One of the strange little quirks/features in Python that often confuses people comes up when assigning and comparing arrays of values.
# Here is a quick example.  Let's start by defining a 1-D array called $a$:

a = collect(range(start = 1.0, stop = 5.0, length = 5))

#OK, so we have an array $a$, with the values 1 through 5.  I want to make a copy of that array, called $b$, so I'll try the following:

b = a

# Great.  So $a$ has the values 1 through 5 and now so does $b$.
# Now that I have a backup of $a$, I can change its values without worrying about losing data (or so I may think!).

a[2] = 17
a

#Here, the 2nd element of $a$ has been changed to 17.  Now let's check on $b$.
b

# And that's how things go wrong!  When you use a statement like $a = b$, rather than copying all the values of $a$ into a new array called $b$,
# Julia just creates an alias (or a pointer) called $b$ and tells it to route us to $a$.
# So if we change a value in $a$ then $b$ will reflect that change (technically, this is called *assignment by reference*).
# If you want to make a true copy of the array, you have to tell Julia to copy every element of $a$ into a new array.  Let's call it $c$.

c = copy(a)

# Now, we can try again to change a value in $a$ and see if the changes are also seen in $c$.

a[2] = 3
a
c

# OK, it worked!  If the difference between `a = b` and `a = b.copy()` is unclear, you should read through this again.  This issue will come back to haunt you otherwise.

## Learn More

# There are a lot of resources online to learn more about using Julia and other libraries.
# Just for kicks, here we use Pluto's feature for embedding videos to point you to a short video on YouTube on using NumPy arrays.

## TODO
#from IPython.display import YouTubeVideo
# a short video about using NumPy arrays, from Enthought
#YouTubeVideo('vWkb7VahaXQ')

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
function surf(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.surf(gr.ptr, Data(z).ptr, stl, opt)
end

function surf(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.surf_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function mesh(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.mesh(gr.ptr, Data(z).ptr, stl, opt)
end

function mesh(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.mesh_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function fall(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.fall(gr.ptr, Data(z).ptr, stl, opt)
end

function fall(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.fall_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function belt(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.belt(gr.ptr, Data(z).ptr, stl, opt)
end

function belt(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.belt_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function boxs(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.boxs(gr.ptr, Data(z).ptr, stl, opt)
end

function boxs(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.boxs_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function tile(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.tile(gr.ptr, Data(z).ptr, stl, opt)
end

function tile(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.tile_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function dens(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.dens(gr.ptr, Data(z).ptr, stl, opt)
end

function dens(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.dens_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function cont(gr::Graph, v::Array{cmgl.Float, 1}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.cont_val(gr.ptr, Data(v).ptr, Data(z).ptr, stl, opt)
end

function cont(gr::Graph, v::Array{cmgl.Float, 1}, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.cont_xy_val(gr.ptr, Data(v).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function cont(gr::Graph, v::cmgl.Float, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.cont_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(z).ptr, stl, opt)
end

function cont(gr::Graph, v::cmgl.Float, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.cont_xy_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function cont(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.cont(gr.ptr, Data(z).ptr, stl, opt)
end

function cont(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.cont_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function contf(gr::Graph, v::Array{cmgl.Float, 1}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contf_val(gr.ptr, Data(v).ptr, Data(z).ptr, stl, opt)
end

function contf(gr::Graph, v::Array{cmgl.Float, 1}, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contf_xy_val(gr.ptr, Data(v).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function contf(gr::Graph, v::cmgl.Float, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contf_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(z).ptr, stl, opt)
end

function contf(gr::Graph, v::cmgl.Float, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contf_xy_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function contf(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contf(gr.ptr, Data(z).ptr, stl, opt)
end

function contf(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contf_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function contd(gr::Graph, v::Array{cmgl.Float, 1}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contd_val(gr.ptr, Data(v).ptr, Data(z).ptr, stl, opt)
end

function contd(gr::Graph, v::Array{cmgl.Float, 1}, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contd_xy_val(gr.ptr, Data(v).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function contd(gr::Graph, v::cmgl.Float, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contd_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(z).ptr, stl, opt)
end

function contd(gr::Graph, v::cmgl.Float, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contd_xy_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function contd(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contd(gr.ptr, Data(z).ptr, stl, opt)
end

function contd(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contd_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function contv(gr::Graph, v::Array{cmgl.Float, 1}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contv_val(gr.ptr, Data(v).ptr, Data(z).ptr, stl, opt)
end

function contv(gr::Graph, v::Array{cmgl.Float, 1}, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contv_xy_val(gr.ptr, Data(v).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function contv(gr::Graph, v::cmgl.Float, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contv_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(z).ptr, stl, opt)
end

function contv(gr::Graph, v::cmgl.Float, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contv_xy_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function contv(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contv(gr.ptr, Data(z).ptr, stl, opt)
end

function contv(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.contv_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function axial(gr::Graph, v::Array{cmgl.Float, 1}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.axial_val(gr.ptr, Data(v).ptr, Data(z).ptr, stl, opt)
end

function axial(gr::Graph, v::Array{cmgl.Float, 1}, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.axial_xy_val(gr.ptr, Data(v).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function axial(gr::Graph, v::cmgl.Float, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.axial_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(z).ptr, stl, opt)
end

function axial(gr::Graph, v::cmgl.Float, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.axial_xy_val(gr.ptr, Data(cmgl.Float[v]).ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function axial(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.axial(gr.ptr, Data(z).ptr, stl, opt)
end

function axial(gr::Graph, x::Array{cmgl.Float, 2}, y::Array{cmgl.Float, 2}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.axial_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function grid(gr::Graph, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.grid(gr.ptr, Data(z).ptr, stl, opt)
end

function grid(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float, 2}; stl::String="", opt::String="")
    cmgl.grid_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end


text:
__precompile__(true)

"""
    Pandoc

Pandoc wrapper to read JSON AST from `pandoc`

See https://hackage.haskell.org/package/pandoc-types-1.17.5.4/docs/Text-Pandoc-Definition.html
"""
module Pandoc

import JSON
import Markdown

const PANDOC_JL_EXECUTABLE = get(ENV, "PANDOC_JL_EXECUTABLE", "pandoc")

@enum Alignment AlignLeft=1 AlignRight=2 AlignCenter=3 AlignDefault=4
@enum ListNumberStyle DefaultStyle=1 Example=2 Decimal=3 LowerRoman=4 UpperRoman=5 LowerAlpha=6 UpperAlpha=7
@enum ListNumberDelim DefaultDelim=1 Period=2 OneParen=3 TwoParens=4
@enum QuoteType SingleQuote=1 DoubleQuote=2
@enum MathType DisplayMath=1 InlineMath=2
@enum CitationMode AuthorInText=1 SuppressAuthor=2 NormalCitation=3

abstract type Element end

abstract type Inline <: Element end
abstract type Block <: Element end

const Format = String
const TableCell = Vector{Block}

struct Attributes
    identifier::String
    classes::Vector{String}
    attributes::Vector{Pair{String, String}}
end
function Attributes(identifier, classes, attributes::Vector{Any})
    attributes = Pair{String, String}[attr[1] => attr[2] for attr in attributes]
    return Attributes(identifier, classes, attributes)
end
Attributes() = Attributes("", [], [],)

struct Citation
    mode::CitationMode
    prefix::Vector{Inline}
    hash::Int
    id::String
    suffix::Vector{Inline}
    note_number::Int
end

struct ListAttributes
    number::Int
    style::ListNumberStyle
    delim::ListNumberDelim
end

"""Plain text, not a paragraph"""
struct Plain <: Block
    content:: Vector{Inline}
end

"""Paragraph"""
struct Para <: Block
    content::Vector{Inline}
end

"""Multiple non-breaking lines"""
struct LineBlock <: Block
    content::Vector{Vector{Inline}}
end

"""Code block (literal) with attributes"""
struct CodeBlock <: Block
    attr::Attributes
    content::String
end
CodeBlock(content) = CodeBlock(Attributes(), content)

"""Raw block"""
struct RawBlock <: Block
    format::Format
    content::String
end

"""Block quote (list of blocks)"""
struct BlockQuote <: Block
    content::Vector{Block}
end

"""Ordered list (attributes and a list of items, each a list of blocks)"""
struct OrderedList <: Block
    attr::ListAttributes
    content::Vector{Vector{Block}}
end

"""Bullet list (list of items, each a list of blocks)"""
struct BulletList <: Block
    content::Vector{Vector{Block}}
end

"""
Definition list Each list item is a pair consisting of a term (a list of inlines) and one or more definitions (each a list of blocks)"""
struct DefinitionList <: Block
    content::Vector{Pair{Vector{Inline}, Vector{Vector{Block}}}}
end

"""Header - level (integer) and text (inlines)"""
struct Header <: Block
    level::Int
    attr::Attributes
    content::Vector{Element}
end
Header(level) = Header(level, Attributes(), [])
Header(level, content) = Header(level, Attributes(), content)

"""Horizontal rule"""
struct HorizontalRule <: Block end

"""Table, with caption, column alignments (required), relative column widths (0 = default), column headers (each a list of blocks), and rows (each a list of lists of blocks)"""
struct Table <: Block
    content::Vector{Inline}
    alignments::Vector{Alignment}
    widths::Vector{Float64}
    headers::Vector{TableCell}
    rows::Vector{Vector{TableCell}}
end

"""Generic block container with attributes"""
struct Div <: Block
    attr::Attributes
    content::Vector{Block}
end
Div(content) = Div(Attributes(), content)

struct Null <: Block end

struct Target
    url::String
    title::String
end

"""Text (string)"""
struct Str <: Inline
    content::String
end

"""Emphasized text (list of inlines)"""
struct Emph <: Inline
    content::Vector{Inline}
end

"""Strongly emphasized text (list of inlines)"""
struct Strong <: Inline
    content::Vector{Inline}
end

"""Strikeout text (list of inlines)"""
struct Strikeout <: Inline
    content::Vector{Inline}
end

"""Superscripted text (list of inlines)"""
struct Superscript <: Inline
    content::Vector{Inline}
end

"""Subscripted text (list of inlines)"""
struct Subscript <: Inline
    content::Vector{Inline}
end

"""Small caps text (list of inlines)"""
struct SmallCaps <: Inline
    content::Vector{Inline}
end

"""Quoted text (list of inlines)"""
struct Quoted <: Inline
    quote_type::QuoteType
    content::Vector{Inline}
end

"""Citation (list of inlines)"""
struct Cite <: Inline
    citations::Vector{Citation}
    content::Vector{Inline}
end

"""Inline code (literal)"""
struct Code <: Inline
    attr::Attributes
    content::String
end
Code(content) = Code(Attributes(), content)

"""Inter-word space"""
struct Space <: Inline end

"""Soft line break"""
struct SoftBreak <: Inline end

"""Hard line break"""
struct LineBreak <: Inline end

"""TeX math (literal)"""
struct Math <: Inline
    math_type::MathType
    content::String
end

"""Raw inline"""
struct RawInline <: Inline
    format::Format
    content::String
end

"""Hyperlink: alt text (list of inlines), target"""
struct Link <: Inline
    attr::Attributes
    content::Vector{Inline}
    target::Target
end
Link(content, target) = Link(Attributes(), content, target)

"""Image: alt text (list of inlines), target"""
struct Image <: Inline
    attr::Attributes
    content::Vector{Inline}
    target::Target
end
Image(content, target) = Image(Attributes(), content, target)

"""Footnote or endnote"""
struct Note <: Inline
    content::Vector{Block}
end

"""Generic inline container with attributes"""
struct Span <: Inline
    attr::Attributes
    content::Vector{Inline}
end
Span(content) = Span(Attributes(), content)

struct Unknown
    e
    t
end

pandoc_api_version() = v"1.17.5-4"
pandoc_api_version(v) = pandoc_api_version(v, Val(length(v)))
pandoc_api_version(v, length::Val{0}) = error("Version array has to be length > 0 but got `$v` instead")
pandoc_api_version(v, length::Val{1}) = VersionNumber(v[1])
pandoc_api_version(v, length::Val{2}) = VersionNumber(v[1], v[2])
pandoc_api_version(v, length::Val{3}) = VersionNumber(v[1], v[2], v[3])
pandoc_api_version(v, length) = VersionNumber(v[1], v[2], v[3], tuple(v[4:end]...))

mutable struct Document
    data::Dict{String, Any}
    pandoc_api_version::VersionNumber
    meta::Dict{String, Any}
    blocks::Vector{Element}
end

function Document(data)
    pav = pandoc_api_version(data["pandoc-api-version"])
    meta = data["meta"]
    blocks = get_elements(data["blocks"])
    return Document(data, pav, meta, blocks)
end

function get_element(t, e)
    u = Unknown(e, t)
    @warn "Unknown element parsed: $u"
    return u
end

get_element(::Type{Null}, e) = Null()
get_element(::Type{SoftBreak}, e) = SoftBreak()
get_element(::Type{LineBreak}, e) = LineBreak()
get_element(::Type{HorizontalRule}, e) = HorizontalRule()
get_element(::Type{Space}, e) = Space()

get_element(::Type{Str}, e) = Str(e["c"])

function get_element(::Type{LineBlock}, e)
    content = Vector{Inline}[]
    for inlines in e["c"]
        push!(content, Inline[get_element(i) for i in inlines])
    end
    return LineBlock(content)
end

function get_element(::Type{Superscript}, e)
    return Superscript(Inline[get_element(i) for i in e["c"]])
end

function get_element(::Type{Subscript}, e)
    return Subscript(Inline[get_element(i) for i in e["c"]])
end

function get_element(::Type{Note}, e)
    return Note(Block[get_element(i) for i in e["c"]])
end

function get_element(::Type{Math}, e)
    math_type = eval(Symbol(e["c"][1]["t"]))
    content = e["c"][2]::String
    return Math(math_type, content)
end

function get_element(::Type{RawInline}, e)
    format = e["c"][1]::String
    content = e["c"][2]::String
    return RawInline(format, content)
end

function get_element(::Type{RawBlock}, e)
    format = e["c"][1]::String
    content = e["c"][2]::String
    return RawBlock(format, content)
end

function get_element(::Type{Quoted}, e)
    quote_type = eval(Symbol(e["c"][1]["t"]))
    content = Inline[get_element(i) for i in e["c"][2]]
    return Quoted(quote_type, content)
end

function get_element(::Type{Table}, e)
    c = e["c"]

    content = Inline[get_element(i) for i in c[1]]
    alignments = Alignment[ eval(Symbol(a["t"])) for a in c[2] ]
    widths = Float64[w for w in c[3]]
    headers = TableCell[]
    rows = Vector{TableCell}[]

    for blocks in c[4]
        push!(headers, Block[ get_element(b) for b in blocks ])
    end

    for list_of_blocks in c[5]
        row = TableCell[]
        for blocks in list_of_blocks
            push!(row, Block[ get_element(b) for b in blocks ])
        end
        push!(rows, row)
    end

    return Table(content, alignments, widths, headers, rows)
end

function get_element(::Type{Span}, e)
    return Span(
                Attributes(e["c"][1]...),
                Inline[get_element(i) for i in e["c"][2]]
               )
end

function get_element(::Type{Image}, e)
    attr = Attributes(e["c"][1]...)
    content = Inline[get_element(i) for i in e["c"][2]]
    target = Target(e["c"][3]...)
    return Image(attr, content, target)
end

function get_element(::Type{Strikeout}, e)
    return Strikeout(Inline[get_element(i) for i in e["c"]])
end

function get_element(::Type{SmallCaps}, e)
    return SmallCaps(Inline[get_element(i) for i in e["c"]])
end

function get_element(::Type{Strong}, e)
    return Strong(Inline[get_element(i) for i in e["c"]])
end

function get_element(::Type{DefinitionList}, e)
    c = e["c"]
    dl = Vector{Pair{Vector{Inline}, Vector{Vector{Block}}}}( [] )
    for items in c
        vector_blocks = Vector{Block}[]
        for blocks in items[2]
            push!(Block[get_element(b) for b in blocks])
        end
        inlines = Inline[get_element(i) for i in items[1]]
        push!(dl, (inlines => vector_blocks))
    end
    return DefinitionList(dl)
end

function get_element(::Type{Div}, e)
    attr = Attributes(e["c"][1]...)
    blocks = Block[get_element(b) for b in e["c"][2]]
    return Div(attr, blocks)
end

function get_element(::Type{BlockQuote}, e)
    blocks = Block[get_element(b) for b in e["c"]]
    return BlockQuote(blocks)
end

function get_element(::Type{ListAttributes}, e)
    number = e[1]
    style = eval(Symbol(e[2]["t"]))
    delim = eval(Symbol(e[3]["t"]))
    return ListAttributes(number, style, delim)
end

function get_element(::Type{OrderedList}, e)
    attr = get_element(ListAttributes, e["c"][1])
    ol = Vector{Block}[]
    for elements in e["c"][2]
        sol = Block[]
        for element in elements
            el = get_element(element)
            push!(sol, el)
        end
        push!(ol, sol)
    end
    return OrderedList(attr, ol)
end

function get_element(::Type{Citation}, e)
    mode = eval(Symbol(e["citationMode"]["t"]))
    prefix = Inline[get_element(i) for i in e["citationPrefix"]]
    hash = e["citationHash"]::Int
    id = e["citationId"]::String
    suffix = Inline[get_element(i) for i in e["citationSuffix"]]
    note_number = e["citationNoteNum"]::Int
    return Citation(mode, prefix, hash, id, suffix, note_number)
end

function get_element(::Type{Cite}, e)
    citations = Citation[get_element(Citation, c) for c in e["c"][1]]
    inlines = Inline[get_element(i) for i in e["c"][2]]
    return Cite(citations, inlines)
end

function get_element(::Type{Code}, e)
    attr = Attributes(e["c"][1]...)
    content = e["c"][2]
    return Code(attr, content)
end

function get_element(::Type{Plain}, e)
    content = Inline[]
    for se in e["c"]
        push!(content, get_element(se))
    end
    return Plain(content)
end

function get_element(::Type{CodeBlock}, e)
    attr = Attributes(e["c"][1]...)
    content = e["c"][2]::String
    return CodeBlock(attr, content)
end

function get_element(::Type{BulletList}, e)
    bl = Vector{Block}[]
    for elements in e["c"]
        sbl = Block[]
        for element in elements
            el = get_element(element)
            push!(sbl, el)
        end
        push!(bl, sbl)
    end
    return BulletList(bl)
end

function get_element(::Type{Emph}, e)
    return Emph(Element[get_element(se) for se in e["c"]])
end

function get_element(::Type{Para}, e)
    return Para(Element[get_element(se) for se in e["c"]])
end

function get_element(::Type{Link}, e)
    c = e["c"]
    identifier = c[1][1]::String
    classes = String[s for s in c[1][2]]
    attributes = Pair[ (String(k) => String(v)) for (k,v) in c[1][3] ]

    content = Element[]
    for se in c[2]
        push!(content, get_element(se))
    end

    target = Target(c[3][1], c[3][2])

    return Link(Attributes(identifier, classes, attributes), content, target)
end

function get_element(::Type{Header}, e)
    c = e["c"]
    level = c[1]::Int
    identifier = c[2][1]::String
    classes = String[s for s in c[2][2]]
    attributes = Pair[ (String(k) => String(v)) for (k,v) in c[2][3] ]

    content = Element[]
    for se in c[3]
        push!(content, get_element(se))
    end

    return Header(level, Attributes(identifier, classes, attributes), content)
end

get_element(e) = get_element(Val{Symbol(e["t"])}(), e)
get_element(::Val{:Header}, e) = get_element(Header, e)
get_element(::Val{:Link}, e) = get_element(Link, e)
get_element(::Val{:Space}, e) = get_element(Space, e)
get_element(::Val{:Emph}, e) = get_element(Emph, e)
get_element(::Val{:Str}, e) = get_element(Str, e)
get_element(::Val{:Para}, e) = get_element(Para, e)
get_element(::Val{:HorizontalRule}, e) = get_element(HorizontalRule, e)
get_element(::Val{:BulletList}, e) = get_element(BulletList, e)
get_element(::Val{:Plain}, e) = get_element(Plain, e)
get_element(::Val{:Code}, e) = get_element(Code, e)
get_element(::Val{:CodeBlock}, e) = get_element(CodeBlock, e)
get_element(::Val{:LineBreak}, e) = get_element(LineBreak, e)
get_element(::Val{:Cite}, e) = get_element(Cite, e)
get_element(::Val{:BlockQuote}, e) = get_element(BlockQuote, e)
get_element(::Val{:OrderedList}, e) = get_element(OrderedList, e)
get_element(::Val{:DefinitionList}, e) = get_element(DefinitionList, e)
get_element(::Val{:Strong}, e) = get_element(Strong, e)
get_element(::Val{:SmallCaps}, e) = get_element(SmallCaps, e)
get_element(::Val{:Span}, e) = get_element(Span, e)
get_element(::Val{:Strikeout}, e) = get_element(Strikeout, e)
get_element(::Val{:Image}, e) = get_element(Image, e)
get_element(::Val{:Table}, e) = get_element(Table, e)
get_element(::Val{:SoftBreak}, e) = get_element(SoftBreak, e)
get_element(::Val{:Quoted}, e) = get_element(Quoted, e)
get_element(::Val{:RawInline}, e) = get_element(RawInline, e)
get_element(::Val{:RawBlock}, e) = get_element(RawBlock, e)
get_element(::Val{:Math}, e) = get_element(Math, e)
get_element(::Val{:Superscript}, e) = get_element(Superscript, e)
get_element(::Val{:Subscript}, e) = get_element(Subscript, e)
get_element(::Val{:Note}, e) = get_element(Note, e)
get_element(::Val{:LineBlock}, e) = get_element(LineBlock, e)
get_element(::Val{:Div}, e) = get_element(Div, e)
get_element(::Val{:Null}, e) = get_element(Null, e)

function get_elements(blocks)
    elements = Element[]
    for e in blocks
        push!(elements, get_element(e))
    end
    return elements
end

function parse(markdown)
    cmd = pipeline(`echo $markdown`, `$PANDOC_JL_EXECUTABLE -t json`)
    data = read(cmd, String)
    return Document(JSON.parse(data))
end

function parse_file(filename)
    cmd = `$PANDOC_JL_EXECUTABLE -t json $filename`
    data = read(cmd, String)
    return Document(JSON.parse(data))
end

exists() = run(`which $PANDOC_JL_EXECUTABLE`, wait=false) |> success

### Convert

function Base.convert(::Type{Document}, md::Markdown.MD)

    pav = v"1.17.5-4"
    meta = Dict{String, Any}()
    data = Dict{String, Any}()
    blocks = Element[]

    for e in md.content
        push!(blocks, convert(Element, e))
    end

    return Document(data, pav, meta, blocks)

end

Base.convert(::Type{Element}, e::AbstractString) = convert(Str, e)
Base.convert(::Type{Inline}, e::AbstractString) = convert(Str, e)
Base.convert(::Type{Str}, e::AbstractString) = Str(e)

Base.convert(::Type{Element}, e::Markdown.Header) = convert(Header, e)
function Base.convert(::Type{Header}, e::Markdown.Header{V}) where V
    content = Element[]

    for s in split(e.text[1])
        push!(content, convert(Str, s))
        push!(content, Space())
    end
    if length(content) > 0 && content[end] isa Space
        pop!(content) # remove last Space()
    end

    return Header(
        V #= level =#,
        Attributes(
                replace(lowercase(e.text[1]), " " => "-"),
                [],
                []
               ),
        content,
    )
end

function _convert(::Type{Vector{Inline}}, text::AbstractString)
    return content
end

Base.convert(::Type{Element}, e::Markdown.Link) = convert(Link, e)
Base.convert(::Type{Inline}, e::Markdown.Link) = convert(Link, e)
function Base.convert(::Type{Link}, e::Markdown.Link)

    content = Inline[]
    text = if e.text isa AbstractString
        e.text
    else
        length(e.text) > 0 ? e.text[1] : ""
    end
    for s in split(text)
        push!(content, convert(Str, s))
        push!(content, Space())
    end
    pop!(content) # remove last Space()
    target = Target(e.url, "")

    return Link(
                Attributes(),
                content,
                target,
               )
end

Base.convert(::Type{Element}, e::Markdown.Italic) = convert(Emph, e)
Base.convert(::Type{Inline}, e::Markdown.Italic) = convert(Emph, e)
Base.convert(::Type{Emph}, e::Markdown.Italic) = Emph(Inline[i for i in e.text])

Base.convert(::Type{Element}, e::Markdown.Bold) = convert(Strong, e)
Base.convert(::Type{Inline}, e::Markdown.Bold) = convert(Strong, e)
Base.convert(::Type{Strong}, e::Markdown.Bold) = Strong(Inline[i for i in e.text])

Base.convert(::Type{Element}, e::Markdown.Paragraph) = convert(Para, e)
function Base.convert(::Type{T}, e::Markdown.Paragraph) where T<:Union{Para,Plain}
    content = Inline[]
    for c in e.content
        if c isa AbstractString
            for s in split(c)
                push!(content, convert(Str, s))
                push!(content, Space())
            end
            if length(content) > 0 && content[end] isa Space
                pop!(content) # remove last Space()
            end
        else
            push!(content, c)
        end
        push!(content, Space())
    end
    if length(content) > 0 && content[end] isa Space
        pop!(content) # remove last Space()
    end
    return T(content)
end

Base.convert(::Type{Element}, e::Markdown.HorizontalRule) = convert(HorizontalRule, e)
Base.convert(::Type{HorizontalRule}, e::Markdown.HorizontalRule) = HorizontalRule()

Base.convert(::Type{Element}, e::Markdown.LineBreak) = convert(LineBreak, e)
Base.convert(::Type{Inline}, e::Markdown.LineBreak) = convert(LineBreak, e)
Base.convert(::Type{LineBreak}, e::Markdown.LineBreak) = LineBreak()

Base.convert(::Type{Element}, e::Markdown.BlockQuote) = convert(BlockQuote, e)
Base.convert(::Type{Block}, e::Markdown.BlockQuote) = convert(BlockQuote, e)
function Base.convert(::Type{BlockQuote}, e::Markdown.BlockQuote)

    content = Block[]

    for c in e.content
        if c isa Markdown.Paragraph
            push!(content, convert(Para, c))
        else
            push!(content, c)
        end
    end

    return BlockQuote(content)
end

Base.convert(::Type{Element}, e::Markdown.Code) = convert(CodeBlock, e)
Base.convert(::Type{Block}, e::Markdown.Code) = convert(CodeBlock, e)
Base.convert(::Type{CodeBlock}, e::Markdown.Code) = CodeBlock(Attributes("", [e.language], []), e.code)

Base.convert(::Type{Inline}, e::Markdown.Code) = convert(Code, e)
Base.convert(::Type{Code}, e::Markdown.Code) = Code(Attributes("", [e.language], []), e.code)

Base.convert(::Type{Element}, e::Markdown.List) = e.ordered == 1 ? convert(OrderedList, e) : convert(BulletList, e)
Base.convert(::Type{Block}, e::Markdown.List) = e.ordered == 1 ? convert(OrderedList, e) : convert(BulletList, e)

function Base.convert(::Type{BulletList}, e::Markdown.List)
    content = Vector{Block}[]
    for items in e.items
        block = Block[]
        for item in items
            if item isa Markdown.Paragraph
                push!(block, convert(Plain, item))
            else
                push!(block, item)
            end
        end
        push!(content, block)
    end
    return BulletList(content)
end

function Base.convert(::Type{OrderedList}, e::Markdown.List)
    content = Vector{Block}[]
    for items in e.items
        block = Block[]
        for item in items
            if item isa Markdown.Paragraph
                push!(block, convert(Plain, item))
            else
                push!(block, item)
            end
        end
        push!(content, block)
    end
    # always returns list level 1 with Decimal with Period
    return OrderedList(ListAttributes(1, Pandoc.Decimal, Pandoc.Period), content)
end

Base.convert(::Type{Element}, e::Markdown.LaTeX) = convert(Inline, e)
Base.convert(::Type{Inline}, e::Markdown.LaTeX) = convert(Math, e)
Base.convert(::Type{Math}, e::Markdown.LaTeX) = Math(InlineMath, e.formula)

Base.convert(::Type{Inline}, e::Markdown.Image) = convert(Image, e)
Base.convert(::Type{Image}, e::Markdown.Image) = Image(Attributes(), [], Target(e.url, e.alt))

Base.convert(::Type{Element}, e::Markdown.Footnote) = convert(Inline, e)
Base.convert(::Type{Inline}, e::Markdown.Footnote) = convert(Note, e)
function Base.convert(::Type{Note}, e::Markdown.Footnote)
    isnothing(e.text) && return Note([])
    content = Block[]
    for p in e.text
        push!(content, p)
    end
    return Note(content)
end

### overloads

# fallback
compare(x, y) = false

function compare(x::T, y::T) where T
    equality = true
    for fieldname in fieldnames(T)
        equality = equality && ( getfield(x, fieldname) == getfield(y, fieldname) )
    end
    return equality
end

Base.:(==)(x::T, y::T) where T<:Pandoc.Element = compare(x, y)
Base.:(==)(x::Attributes, y::Attributes) = compare(x, y)
Base.:(==)(x::Citation, y::Citation) = compare(x, y)
Base.:(==)(x::ListAttributes, y::ListAttributes) = compare(x, y)
Base.:(==)(x::Target, y::Target) = compare(x, y)

end # module


text:
using Clang.wrap_c

include("constants.jl")
headers = [joinpath(LIB_FOLDER, "iaea_phsp.h")]
wc = wrap_c.init(
    headers=headers,
    output_file="raw.jl",
    header_library = _ -> :IAEA_LIB,
)
wc.headers = headers
run(wc)


text:
module GraphUtilities

using GenericTensorNetworks, Graphs
using JSON, DelimitedFiles
using Configurations, Random
using LegibleLambdas

using Requires
function __init__()
    @require UnitDiskMapping="1b61a8d9-79ed-4491-8266-ef37f39e1727" using UnitDiskMapping
end

export GraphProblemConfig, problem_instance, foldername
export SmallGraphConfig, DiagGraphConfig, SquareGraphConfig, RegularGraphConfig, MISProjectGraphConfig, MappedRegularGraphConfig
export save_property, load_property
export unified_solve

include("config.jl")
include("fileio.jl")
include("json.jl")
include("verify.jl")

end


text:
"""
function build_variable_dictionary(data::Dict)
    Builds a dictionary which is necessary to build the H matrix.
    This dictionary allows to pass the correct variables to each measurement function.
"""
function build_variable_dictionary(data::Dict)

    bus_terms = [bus["terminals"] for (b, bus) in data["bus"]]
    ref_bus = [bus["index"] for (b, bus) in data["bus"] if bus["bus_type"]==3]
    @assert length(ref_bus) == 1 "There is more than one reference bus, double-check model"
    ref_bus_term = data["bus"]["$(ref_bus[1])"]["terminals"]

    total_vm = sum(length.(bus_terms))
    total_va = sum(length.(bus_terms)) - length(ref_bus_term)

    variable_dict = Dict{String, Any}()
    variable_dict["vm"] = Dict{String, Any}()
    variable_dict["va"] = Dict{String, Any}()

    vmcount = 1
    for (b, bus) in data["bus"]
        variable_dict["vm"][b] = Dict{String, Any}()
        for i in bus["terminals"]
            variable_dict["vm"][b]["$i"] = Dict{String, Any}()
            variable_dict["vm"][b]["$i"] = vmcount
            vmcount+=1
        end
    end

    vacount = total_vm+1
    for (b, bus) in data["bus"]
        if b != "$(ref_bus[1])"
            variable_dict["va"][b] = Dict{String, Any}()
            for i in bus["terminals"]
                variable_dict["va"][b]["$i"] = Dict{String, Any}()
                variable_dict["va"][b]["$i"] = vacount
                vacount+=1
            end
        end
    end
    return variable_dict
end

function variables_of_buses(Bi::Array, variable_dict::Dict)
    vm_indices = []
    va_indices = []
    for bus in Bi
        push!(vm_indices, variable_dict["vm"]["$bus"])
        if haskey(variable_dict["va"], "$bus") push!(va_indices, variable_dict["va"]["$bus"]) end # if it is a ref bus, the key does not exist
    end
    return vm_indices, va_indices
end
"""
Function to build multi branch input for the injection functions (as injections can occur at buses with multiple incoming/outcoming branches,
    which is not the case for 'flows', which only refer to a single branch )
"""
function build_multibranch_input(qmeas, meas, data, variable_dict)

    cmp_bus = occursin("g", string(qmeas)) ? data["gen"]["$(meas["cmp_id"])"]["gen_bus"] : data["load"]["$(meas["cmp_id"])"]["load_bus"]
    vm_cmp_bus = variable_dict["vm"]["$cmp_bus"]
    va_cmp_bus = haskey(variable_dict["va"], "$cmp_bus") ? variable_dict["va"]["$cmp_bus"] : Dict{String, Any}() # is empty if cmp_bus is the ref bus
    conn_branches = [branch for (_, branch) in data["branch"] if (branch["f_bus"] == cmp_bus || branch["t_bus"] == cmp_bus)]
    g = []
    b = []
    G_fr = []
    B_fr = []
    t_buses = []
    f_conn = []
    t_conn = []
    for branch in conn_branches
        push!(g, _PMD.calc_branch_y(branch)[1])
        push!(b, _PMD.calc_branch_y(branch)[2])
        push!(G_fr, branch["g_fr"])
        push!(B_fr, branch["b_fr"])
        push!(t_buses, branch["t_bus"])
        push!(t_buses, branch["f_bus"])
        push!(f_conn, branch["f_connections"])
        push!(t_conn, branch["t_connections"])
    end
    t_buses = filter(x-> x!=cmp_bus, t_buses) |> unique
    vm_adj_bus, va_adj_bus = variables_of_buses(t_buses, variable_dict)
    return cmp_bus, vm_cmp_bus, va_cmp_bus, conn_branches, f_conn, t_conn, g, b, G_fr, B_fr, t_buses, vm_adj_bus, va_adj_bus
end
"""
Function to build single branch input for the flow measurement functions (as 'flows' only refer to a single branch, as opposed to
    injections, which can occur at buses with multiple incoming/outcoming branches)
"""
function build_singlebranch_input(meas, data, variable_dict)

    branch = data["branch"]["$(meas["cmp_id"])"]
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    vm_fr = variable_dict["vm"]["$f_bus"]
    va_fr = haskey(variable_dict["va"], "$f_bus") ? variable_dict["va"]["$f_bus"] : Dict{String, Any}() # is empty if cmp_bus is the ref bus

    g, b = _PMD.calc_branch_y(branch)
    G_fr = branch["g_fr"]
    B_fr = branch["b_fr"]

    vm_to, va_to = variables_of_buses([t_bus], variable_dict)
    return f_bus, vm_fr, va_fr, [1], [branch["f_connections"]], [branch["t_connections"]], [g], [b], [G_fr], [B_fr], t_bus, vm_to, va_to
end
"""
    Given a state estimation or powerflow result or similar dictionary `se_sol`, and the system variables `variable_dict`,
    this function builds an array with the numerical state of the system, e.g., variables in `variable_dict`
        are assigned a scalar that correspond to their value in `se_sol`.
        The resulting array has the same index order of that of `h_functions`.
"""
function build_state_array(sol_dict::Dict, variable_dict::Dict)
    state_length = sum( length(val) for (_, val) in variable_dict["vm"]) + sum( length(val) for (_, val) in variable_dict["va"])
    state_array = zeros(Float64, (state_length,) )
    for (key, val) in variable_dict["vm"]
        bus_idx = minimum([idx for (_, idx) in val])
        for i in 1:length(val)
            state_array[bus_idx+i-1] = sol_dict["solution"]["bus"][key]["vm"][i]
        end
    end
    for (key, val) in variable_dict["va"]
        bus_idx = minimum([idx for (_, idx) in val])
        for i in 1:length(val)
            state_array[bus_idx+i-1] = sol_dict["solution"]["bus"][key]["va"][i]
        end
    end
    return state_array
end
"""
Adds "virtual" zero zib residuals to a state estimation solution dictionary `se_sol`.
The `data` dictionary must have measurement data for these zibs, which can be added
    with _PMDSE.add_zib_virtual_meas!
"""
function add_zib_virtual_residuals!(se_sol::Dict, data::Dict)
    original_meas = [m for (m, meas) in se_sol["solution"]["meas"]]
    for (m, meas) in data["meas"]
        if m ∉ original_meas
            se_sol["solution"]["meas"][m] = Dict("res"=>zeros(length(meas["dst"])))
        end
    end
end


text:
# Autogenerated wrapper script for libsingular_julia_jll for x86_64-linux-musl-cxx03-julia_version+1.6.0
export libsingular_julia

using libcxxwrap_julia_jll
using Singular_jll
JLLWrappers.@generate_wrapper_header("libsingular_julia")
JLLWrappers.@declare_library_product(libsingular_julia, "libsingular_julia.so")
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll, Singular_jll)
    JLLWrappers.@init_library_product(
        libsingular_julia,
        "lib/libsingular_julia.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
a = VecE2()
@test typeof(a) <: VecE
@test typeof(a) <: AbstractVec
@test a.x == 0.0
@test a.y == 0.0

b = VecE2(0.0,0.0)
@test b.x == 0.0
@test b.y == 0.0

@test length(a) == 2
@test a == b
@test isequal(a, b)
@test copy(a) == a
@test convert(Vector{Float64}, a) == [0.0,0.0]
@test convert(VecE2, [0.0,0.0]) == a

@test isapprox(polar(1.0,0.0), VecE2(1.0,0.0))
@test isapprox(polar(2.0,π/2), VecE2(0.0,2.0))
@test isapprox(polar(3.0,-π/2), VecE2(0.0,-3.0))
@test isapprox(polar(-0.5,1.0*π), VecE2(0.5,0.0))

a = VecE2(0.0,1.0)
b = VecE2(0.5,2.0)

@test a != b

@test a + b == VecE2(0.5, 3.0)
@test a + 1 == VecE2(1.0, 2.0)
@test a + 0.5 == VecE2(0.5, 1.5)

@test a - b == VecE2(-0.5, -1.0)
@test a - 1 == VecE2(-1.0,  0.0)
@test a - 0.5 == VecE2(-0.5, 0.5)

@test a * 2 == VecE2(0.0, 2.0)
@test a * 0.5 == VecE2(0.0, 0.5)

@test a / 2 == VecE2(0.0, 0.5)
@test a / 0.5 == VecE2(0.0, 2.0)

@test b^2 == VecE2(0.25, 4.0)
@test b^0.5 == VecE2(0.5^0.5, 2.0^0.5)

@test a % 2.0 == VecE2(0.0, 1.0)

@test isfinite(a)
@test !isfinite(VecE2(Inf,0))
@test !isfinite(VecE2(0,Inf))
@test !isfinite(VecE2(0,-Inf))
@test !isfinite(VecE2(Inf,Inf))

@test !isinf(a)
@test isinf(VecE2(Inf,0))
@test isinf(VecE2(0,Inf))
@test isinf(VecE2(0,-Inf))
@test isinf(VecE2(Inf,Inf))

@test !isnan(a)
@test isnan(VecE2(NaN,0))
@test isnan(VecE2(0,NaN))
@test isnan(VecE2(NaN,NaN))

@test round(VecE2(0.25,1.75)) == VecE2(0.0,2.0)
@test round(VecE2(-0.25,-1.75)) == VecE2(-0.0,-2.0)

@test floor(VecE2(0.25,1.75)) == VecE2(0.0,1.0)
@test floor(VecE2(-0.25,-1.75)) == VecE2(-1.0,-2.0)

@test ceil(VecE2(0.25,1.75)) == VecE2(1.0,2.0)
@test ceil(VecE2(-0.25,-1.75)) == VecE2(-0.0,-1.0)

@test trunc(VecE2(0.25,1.75)) == VecE2(0.0,1.0)
@test trunc(VecE2(-0.25,-1.75)) == VecE2(-0.0,-1.0)

@test clamp(VecE2(1.0, 10.0), 0.0, 5.0) == VecE2(1.0, 5.0)
@test clamp(VecE2(-1.0, 4.0), 0.0, 5.0) == VecE2(0.0, 4.0)

c = VecE2(3.0,4.0)

@test isapprox(abs(a), 1.0)
@test isapprox(abs(c), 5.0)
@test isapprox(hypot(a), 1.0)
@test isapprox(hypot(c), 5.0)
@test isapprox(abs2(a), 1.0)
@test isapprox(abs2(c), 25.0)
@test isapprox(norm(a), VecE2(0.0,1.0))
@test isapprox(norm(b), VecE2(0.5/sqrt(4.25),2.0/sqrt(4.25)))
@test isapprox(norm(c), VecE2(3.0/5,4.0/5))

@test isapprox(dist(a,a), 0.0)
@test isapprox(dist(b,b), 0.0)
@test isapprox(dist(c,c), 0.0)
@test isapprox(dist(a,b), hypot(0.5, 1.0))
@test isapprox(dist(a,c), hypot(3.0, 3.0))

@test isapprox(dist2(a,a), 0.0)
@test isapprox(dist2(b,b), 0.0)
@test isapprox(dist2(c,c), 0.0)
@test isapprox(dist2(a,b), 0.5^2 + 1^2)
@test isapprox(dist2(a,c), 3^2 + 3^2)

@test isapprox(dot(a, b), 2.0)
@test isapprox(dot(b, c), 1.5 + 8.0)
@test isapprox(dot(c, c), abs2(c))

@test isapprox(proj(b, a, Float64), 2.0)
@test isapprox(proj(c, a, Float64), 4.0)

@test isapprox(proj(b, a, VecE2), VecE2(0.0, 2.0))
@test isapprox(proj(c, a, VecE2), VecE2(0.0, 4.0))

@test isapprox(lerp(VecE2(0,0), VecE2(2,-2), 0.25), VecE2(0.5,-0.5))
@test isapprox(lerp(VecE2(2,-2), VecE2(3,-3), 0.5), VecE2(2.5,-2.5))

@test isapprox(rot180(b), VecE2(-0.5,-2.0))
@test isapprox(rotr90(b), VecE2( 2.0,-0.5))
@test isapprox(rotl90(b), VecE2(-2.0, 0.5))

@test isapprox(rot(b, 1.0*pi), VecE2(-0.5,-2.0))
@test isapprox(rot(b, -π/2), VecE2( 2.0,-0.5))
@test isapprox(rot(b,  π/2), VecE2(-2.0, 0.5))
@test isapprox(rot(b, 2π), b)


text:
#const XYIterGMMorLevel=Union{XYIterLevel, XYIterGMM}

iterloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterControl, RegressionType::Val;
  ) where {TM, TV, T} = sum(
    (Θ(Xv, RegressionType) .- Xv.ṽ).^2)

function currentloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterControl,
   absμₖ::TM) where {TM<:CuMatrix, TV<:CuVector, T}

  someones::CuVector{T} = CUDA.ones(T, size(Θ.A,1))
  return sum((absμₖ*someones .- Xv.ṽ).^2)
end

function testiterloss(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterControl;
  RegressionType::Val=Val(:choleskyzygote)) where {TM, TV, T}

  vhattest = iterloss(deepcopy(Θ), Xv, RegressionType)

  vhatver = sum((Θ(Xv, RegressionType) .- Xv.ṽ) .^ 2)

  @assert vhatver ≈ vhattest
end

function itermodel(panel::AbstractDataFrame,  zs::ZSpec,::Val{:ols},
    ::Type{T}=PARAM[:itertype],
    ::Type{TV}=PARAM[:itergpu] ? CUDA.CuVector{T} : Vector{T},
    ::Type{TM}=PARAM[:itergpu] ? CUDA.CuMatrix{T} : Matrix{T},
    PredictionType::Val = Val(PARAM[:predictiontype]),
    RegressionType::Val = Val(PARAM[:iterregressiontype]),
    SolveType::Val = Val(PARAM[:itersolvetype]),
    ; kwargs... #in case we want to use kwargs in other methods
    ) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}
  (T == Float32) && (PARAM[:minweightmagnitude]<10^(-7)) && @warn(
    "T=$T while PARAM[:minweightmagnitude]=$(PARAM[:minweightmagnitude]).
    Therefore the minweight is effectivelly 0- is this really what you want?")
  CUDA.allowscalar(false)

  Zms::MeasureSpec = zs.Zms
  Fw::Vector{Symbol} = Fweights(Zms, :Fw)
  FRtw::Vector{Symbol} = Fweights(Zms, :FRtw)
  FRLw::Vector{Symbol} = Fweights(Zms, :FRLw)
  FW::Vector{Symbol} = Zms.FW
  FWgroup = Fgroupedcontrols(Zms)

  #form the data into a more usable format and expand the arrays
  p = modelparts(panel, T,TV,TM,
    weightdata = (;Fw, FRtw, FRLw), Fvol=Zms.Fvol; zs, FW, FWgroup)
  dims = p.dims
  tidx::Dict = p.tidx
  cleanpanel::SubDataFrame = p.cleanpanel
  Ncleanpanel::Int = nrow(cleanpanel)
  ws::TM = p.dat[:Fw]
  Rtws::TM = p.dat[:FRtw]
  RLws::TM = p.dat[:FRLw]
  W::TM = p.dat[:FW]
  Wgroup = p.dat[:FWgroup]
  ts::Vector{Int} = p.ts

  v::TV = p.dat[:Fvol]

  #println("ts: $ts")
  Xv = XYIterControl(ws, Rtws, RLws, v, ts, W=W, Wgroup=Wgroup)

  #construct the parameter
  Θ = VolumePartsIter(dims.T, dims.K, ts, T)

  local λ⁺::T = T(-1)

  #for debugging purposes- will trhow an error at some point if true
  PARAM[:runstacktraceoniter] && solveA!(iterloss,Θ,Xv,SolveType)

  #main optimization routine
  try
    λ⁺ = solveA!(iterloss,Θ,Xv,SolveType)
  catch err
    errmessage = "$err"
    if length(errmessage) > 20_000 #sometimes we get really long error messages
      errmessage = errmessage[1:20_000]
    end
    #@warn "Solver failed with error $errmessage"
    @error "Solver failed with error $errmessage\n" exception=(err, catch_backtrace())

  finally
    results::DataFrame = formresults(panel, Zms, Θ, Xv, Val{:iter}())

    #recompute the loss if needed
    if λ⁺<0
      throw("λ⁺ < 0!! This probably means the gpu failed. If thats not what happened
      it shouldn't be too hard to fix up the code and allow for this scenario. If it is,
      Just restart Julia and rerun.")
      λ⁺ = updateA₁!(Θ.g|>TV|>deepcopy, iterloss, Θ, Xv, RegressionType, Val{:identity}())
      updateAfromGAₜ!(Θ,Xv,1) #update to the final values of A
      #currentloss(Θ, Xv, projectvolumefromA(Θ,Xv))
      λ⁺ = lossfromA(Θ, Xv, RegressionType)
    end

    λ⁺ = isfinite(λ⁺) ? Int(round(min(λ⁺, typemax(Int)÷10))) : T(-99.0)

    return λ⁺, results
  end

  @assert false #should never reach here due to the finalizer
end

function itermodel(panel::AbstractDataFrame,  zs::ZSpec,::Val{:bootstrap},
    ::Type{T}=PARAM[:itertype],
    ::Type{TV}=PARAM[:iterbootstrapgpu] ? CUDA.CuVector{T} : Vector{T},
    ::Type{TM}=PARAM[:iterbootstrapgpu] ? CUDA.CuMatrix{T} : Matrix{T},
    PredictionType::Val = Val(PARAM[:predictiontype]),
    RegressionType::Val = Val(PARAM[:iterregressiontype]);
    funds, kwargs...
    ) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}
  (T == Float32) && throw(
    "If you really want to use Float32 using the bootstrap approach,
    you can disable this.")
  CUDA.allowscalar(false)

  Zms::MeasureSpec = zs.Zms
  Fw::Vector{Symbol} = Fweights(Zms, :Fw)
  FRtw::Vector{Symbol} = Fweights(Zms, :FRtw)
  FRLw::Vector{Symbol} = Fweights(Zms, :FRLw)
  FW::Vector{Symbol} = Zms.FW
  FWgroup = Fgroupedcontrols(Zms)

  #form the data into a more usable format and expand the arrays
  p = modelparts(panel, T,TV,TM,
    weightdata = (;Fw, FRtw, FRLw), Fvol=Zms.Fvol; zs, FW, FWgroup)
  dims = p.dims
  tidx::Dict = p.tidx
  cleanpanel::SubDataFrame = p.cleanpanel
  Ncleanpanel::Int = nrow(cleanpanel)
  ws::TM = p.dat[:Fw]
  Rtws::TM = p.dat[:FRtw]
  RLws::TM = p.dat[:FRLw]
  W::TM = p.dat[:FW]
  Wgroup = p.dat[:FWgroup]
  ts::Vector{Int} = p.ts

  v::TV = p.dat[:Fvol]

  #println("ts: $ts")
  Xv = XYIterControl(ws, Rtws, RLws, v, ts, W=W, Wgroup=Wgroup)

  #construct the parameter
  Θ = VolumePartsIter(dims.T, dims.K, ts, T, TM, TV)

  G = extractGfromfunds(TM; funds, zs, tidx, ts)
  println("dims G: $(size(G))
    dims Θ.G: $(size(Θ.G))")
  Θ.G .= G

  local λ⁺::T = T(-1)

  #for debugging purposes- will trhow an error at some point if true
  PARAM[:runstacktraceoniter] && solveA!(iterloss,Θ,Xv,
    Val{:bootstrap}())

  #main optimization routine
  #try
    λ⁺ = solveA!(iterloss,Θ,Xv,Val{:bootstrap}())
  #=catch err
    errmessage = "$err"
    if length(errmessage) > 20_000 #sometimes we get really long error messages
      errmessage = errmessage[1:20_000]
    end
    #@warn "Solver failed with error $errmessage"
    @error "Solver failed with error $errmessage\n" exception=(err, catch_backtrace())

  finally
    results::DataFrame = formresults(panel, Zms, Θ, Xv, Val{:iter}())

    #recompute the loss if needed
    if λ⁺<0
      throw("λ⁺ < 0!! This probably means the gpu failed. If thats not what happened
      it shouldn't be too hard to fix up the code and allow for this scenario. If it is,
      Just restart Julia and rerun.")
      λ⁺ = updateA₁!(Θ.g|>TV|>deepcopy, iterloss, Θ, Xv, RegressionType, Val{:identity}())
      updateAfromGAₜ!(Θ,Xv,1) #update to the final values of A
      #currentloss(Θ, Xv, projectvolumefromA(Θ,Xv))
      λ⁺ = lossfromA(Θ, Xv, RegressionType)
    end=#

    results::DataFrame = formresults(panel, Zms, Θ, Xv, Val{:iter}())
    results = addfundcolstoresults(results; funds, zs, tidx, ts)

    λ⁺ = isfinite(λ⁺) ? Int(round(min(λ⁺, typemax(Int)÷10))) : T(-99.0)

    return λ⁺, results
  #end

  @assert false #should never reach here due to the finalizer
end

#extracts the G from the fund file
function extractGfromfunds(::Type{TM}; funds, zs, tidx, ts, fundsource=PARAM[:iterfundsource],
  fundassetsprefix = PARAM[:iterfundassetprefix]) where TM <: AbstractMatrix
  ms = zs.originalms
  f = funds[fundsource] |> deepcopy
  alldates = sort(keys(tidx) |> collect)

  @assert issorted(f, :date)
  funddates = f.date

  #verify we have coverage of the growth rates over the estimation period
  @assert (setdiff(alldates, funddates) |> isempty) "
    setdiff(alldates,funddates): $(setdiff(alldates,funddates))"

  cols = ms.Fξs .|> (s->Symbol(fundassetsprefix, s))

  cleanfunds = completecases(f, cols)

  A = f[f.date .|> (d->insorted(d, alldates)), cols] |> Matrix{Float64}
  G = (A[2:end, :] ./ A[1:(end-1),:])' |> TM

  #NOTE: this approach won't work with smoothing- but wouldn't smoothing
  #break the accounting relationship?
  return G
end

function addfundcolstoresults(results; funds, zs, tidx, ts, fundsource=PARAM[:iterfundsource],
  fundassetsprefix = PARAM[:iterfundassetprefix],
  additionalfundscols = PARAM[:iteradditionalfundscols]) where TM <: AbstractMatrix
  ms = zs.originalms
  f = funds[fundsource] |> deepcopy
  alldates = sort(keys(tidx) |> collect)

  @assert issorted(f, :date)
  funddates = f.date

  #verify we have coverage of the growth rates over the estimation period
  @assert (setdiff(alldates, funddates) |> isempty) "
    setdiff(alldates,funddates): $(setdiff(alldates,funddates))"

  cols = [ms.Fξs .|> (s->Symbol(fundassetsprefix, s)); additionalfundscols]

  cleanfunds = completecases(f, cols)

  colstoadd = f[f.date .|> (d->insorted(d, alldates)), [:date; cols;]] |> deepcopy

  newcolnames = cols .|> (s)->Symbol(fundsource, s)
  rename!(colstoadd, Dict(cols .=> newcolnames))

  results = leftjoin(results, colstoadd, on=:date, validate=(true,true))

  #NOTE: this approach won't work with smoothing- but wouldn't smoothing
  #break the accounting relationship?
  return results
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Random

temp = zeros(Int64, 3, 2, 4)

for i in 1:3
    for j in 1:2
        for k in 1:4
            temp[i,j,k] = i*j*k
        end
    end
end

temp

using HDF5

temp[:]

h5write("../DHC/scratch_AKS/order_test.h5", "main/data", temp[:])

for

temp = zeros(Int64, 3, 3)

for i in 1:3
    for j in 1:3
        temp[i,j] = i*(j+1)
    end
end

temp

temp[:]

temp[1,2]

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
# exports
export create, Sphere, show

mutable struct Sphere <: Shape
    radius::Float64
end

function Base.show(io::IO, x::Sphere)
    println(io, "Sphere")
    println(io, "Radius: $(x.radius)")
end

function create(c::Sphere; resolution=nothing, rand_=0.0, type::Int64=1)
    radius = c.radius
    bounds = [-radius radius; -radius radius; -radius radius]
    x, v, y, vol = create(Cuboid(bounds), resolution=resolution, rand_=rand_)
    mask = vec(sum(x.^2, dims=1) .<= radius^2)
    mesh = x[:, mask]
    # x, v, y, vol
    return mesh, zeros(size(mesh)), copy(mesh), vol, id*ones(Int64, length(vol))
end

max_stars_repo_name: ravinderbhattoo/PDMesh

max_issues_repo_name: ravinderbhattoo/PDMesh

max_forks_repo_name: ravinderbhattoo/PDMesh


text:
include("dynamic_piomelli.jl")

@par function realspace_dynamic_les_calculation!(s::A) where {A<:@par(AbstractSimulation)}
    real!(s.lesmodel.û)
    real!(s.lesmodel.M)
    if is_dynP_les(A)
        real!(s.lesmodel.ŵ)
    end
    if is_piomelliSmag_les(A)
        dt = get_dt(s)
        dtm1 = s.lesmodel.dtm1[]
        dtf = dt == 0.0 ? 1.0 : dt
        dtm1f = dtm1 == 0.0 ? dtf : dtm1

        if s.iteration[] == 0
            copyto!(s.lesmodel.cm1.rr,s.lesmodel.c.rr)
        end

        calc_cstar!(s.lesmodel.c.rr,s.lesmodel.cm1.rr,dtf,dtm1f)
        s.lesmodel.dtm1[] = dt
    end

    if is_piomelliP_les(A)
        dt = get_dt(s)
        dtm1 = s.lesmodel.dtpm1[]
        dtf = dt == 0.0 ? 1.0 : dt
        dtm1f = dtm1 == 0.0 ? dtf : dtm1

        if s.iteration[] == 0
            copyto!(s.lesmodel.cpm1.rr,s.lesmodel.cp.rr)
        end

        calc_cstar!(s.lesmodel.cp.rr,s.lesmodel.cpm1.rr,dtf,dtm1f)
        s.lesmodel.dtpm1[] = dt
    end

    @mthreads for j in TRANGE
        realspace_dynamic_les_calculation!(s,j)
    end
    return nothing
end

@par function realspace_dynamic_les_calculation!(s::A,j::Integer) where {A<:@par(AbstractSimulation)}
    u = s.u.rr
    L = s.lesmodel.L.rr
    tau = s.lesmodel.tau.rr
    Sa = s.lesmodel.S.rr
    Δ² = s.lesmodel.Δ²

    if is_piomelliSmag_les(A)
        c = s.lesmodel.c.rr
    end

    if is_dynP_les(A)
        w = s.rhs.rr
        P = s.lesmodel.P.rr
        if is_piomelliP_les(A)
            cp = s.lesmodel.cp.rr
        end
    end

    @inbounds @msimd for i in REAL_RANGES[j]
        S = tau[i]
        if is_piomelliSmag_les(A)
            Sa[i] = 2*c[i]*Δ²*norm(S)*S
        else
            Sa[i] = Δ²*norm(S)*S
        end
        v = u[i]
        L[i] = symouter(v,v)
        if is_piomelliP_les(A)
            P[i] = cp[i]*Δ²*Lie(S,AntiSymTen(0.5*w[i]))
        elseif is_dynSandP_les(A)
            P[i] = Δ²*Lie(S,AntiSymTen(0.5*w[i]))
        end
    end

    return nothing
end

@par function fourierspace_dynamic_les_calculation!(s::A) where {A<:@par(AbstractSimulation)}
    fullfourier!(s.lesmodel.S)
    fullfourier!(s.lesmodel.L)
    is_dynP_les(A) && fullfourier!(s.lesmodel.P)

    ind = CartesianIndices((Base.OneTo(NY),Base.OneTo(NZ)))
    @mthreads for jk in ind
        j,k = Tuple(jk)
        fourierspace_dynamic_les_calculation!(s,j,k)
    end

    real!(s.lesmodel.S)
    real!(s.lesmodel.L)
    is_dynP_les(A) && real!(s.lesmodel.P)
    return nothing
end

@par function fourierspace_dynamic_les_calculation!(s::A,j::Integer, k::Integer) where {A<:@par(AbstractSimulation)}
    S = s.lesmodel.S.c
    L = s.lesmodel.L.c
    Δ̂² = s.lesmodel.Δ̂²
    if is_dynP_les(A)
        P = s.lesmodel.P.c
    end

    @inbounds begin
        @msimd for i in Base.OneTo(NX)
            G = Gaussfilter(Δ̂²,i,j,k)
            S[i,j,k] *= G
            L[i,j,k] *= G
            if is_dynP_les(A)
                P[i,j,k] *= G
            end
        end
    end
end

@par function coefficient_les_calculation!(s::A) where {A<:@par(AbstractSimulation)}
    @mthreads for j in TRANGE
        coefficient_les_calculation!(s,j)
    end

    s.lesmodel.avg && average_c!(s.lesmodel.c,s.lesmodel.Δ̂²)

    if is_dynP_les(s)
        s.lesmodel.avg && average_c!(s.lesmodel.cp,s.lesmodel.Δ̂²)
    end

    return nothing
end

@par function coefficient_les_calculation!(s::A,j::Integer) where {A<:@par(AbstractSimulation)}

    rhs = s.rhs.rr
    τ = s.lesmodel.tau.rr
    Δ² = s.lesmodel.Δ²

    if haspassivescalarles(A)
        φ = s.passivescalar.φ.field.data
        fφ = s.passivescalar.flux.rr
        has_les_scalar_vorticity_model(A) && (cφ = s.passivescalar.lesmodel.c*Δ²)
    end

    if hasdensityles(A)
        ρ = s.densitystratification.ρ.field.data
        f = s.densitystratification.flux.rr
        has_les_density_vorticity_model(A) && (cρ = s.densitystratification.lesmodel.c*Δ²)
    end

    ca = s.lesmodel.c.rr
    û = s.lesmodel.û.rr
    La = s.lesmodel.L.rr
    Ma = s.lesmodel.M.rr
    Sa = s.lesmodel.S.rr
    Δ̂² = s.lesmodel.Δ̂² + Δ²
    #cmin = s.lesmodel.cmin

    if is_dynP_les(A)
        Pa = s.lesmodel.P.rr
        ŵ = s.lesmodel.ŵ.rr
        cpa = s.lesmodel.cp.rr
    end

    @inbounds @msimd for i in REAL_RANGES[j]
        w = rhs[i]
        S = τ[i]

        L = traceless(symouter(û[i],û[i]) - La[i])
        La[i] = L
        Sh = Ma[i]

        if is_dynSmag_les(A)
            M = Δ̂²*norm(Sh)*Sh - Sa[i]
            #c = max(0.0,0.5*((L:M)/(M:M)))
            c = max(-1.0,min(1.0,0.5*((L:M)/(M:M))))
        elseif is_piomelliSmag_les(A)

            τ̂ = Sa[i]

            if is_piomelliP_les(A)
                τ̂ += Pa[i]
            end

            ah = norm(Sh)

            M = 2*Δ̂²*ah*Sh
            Mn = M/(M:M)
            c = (τ̂ + L):Mn
            c = max(-1.0,min(1.0,c))
            #c = ifelse(ah<10.0,0.0,max(0.0,c))
        end
        ca[i] = c

        if is_dynP_les(A)
            wh = ŵ[i]
            Ph = Lie(Sh,AntiSymTen(0.5*wh))

            if is_dynSandP_les(A)
                Mp = Δ̂²*Ph - Pa[i]
                cp = max(-1.0,min(1.0,(L:Mp)/(Mp:Mp))) # Should I clip negative values?
                #cp = max(0.0,cp)
            elseif is_piomelliP_les(A)
                if !is_piomelliSmag_les(A)
                    τ̂ = Pa[i]
                end
                Mp = Δ̂²*Ph
                Mpn = Mp/(Mp:Mp)
                cp = max(-1.0,min(1.0,(τ̂ + L):Mpn))
            end

            cpa[i] = cp
        end

    end

    return nothing
end


text:
# Importing the physical constants:
constants = Constants()

T = 28.0 - constants.K₀ # Current temperature
Tᵣ = 25.0 - constants.K₀ # Reference temperature
A = Fvcb()

@testset "Γ_star()" begin
    @test PlantBiophysics.Γ_star(T,Tᵣ,constants.R) ==
        PlantBiophysics.arrhenius(42.75,37830.0,T,Tᵣ,constants.R)
end;

@testset "standard arrhenius()" begin
    @test PlantBiophysics.arrhenius(42.75,37830.0,T,Tᵣ,constants.R) ≈ 49.76935360399572
end;

@testset "arrhenius() with negative effect of too high T" begin
    @test PlantBiophysics.arrhenius(A.JMaxRef,A.Eₐⱼ,T,Tᵣ,A.Hdⱼ,A.Δₛⱼ) ≈ 278.5161762418
    # NB: value checked using plantecophys.
end;

@testset "arrhenius() with negative effect of too high T" begin
    @test PlantBiophysics.arrhenius(A.JMaxRef,A.Eₐⱼ,T,Tᵣ,A.Hdⱼ,A.Δₛⱼ) ≈ 278.5161762418
    # NB: value checked using plantecophys.
end;

@testset "compare arrhenius() implementations" begin
# arrhenius with negative effect of too high T should yield the same result as the standard Arrhenius
# when Δₛ = 0.0
    @test PlantBiophysics.arrhenius(A.JMaxRef,A.Eₐⱼ,28.0-constants.K₀,A.Tᵣ-constants.K₀,A.Hdⱼ,0.0) ==
        PlantBiophysics.arrhenius(A.JMaxRef,A.Eₐⱼ,28.0-constants.K₀,A.Tᵣ-constants.K₀)
end;


text:
## Utility functions

# Sum of square differences
sosdiff(a::Number,b::Number) = abs2(a-b)
sosdiff(A,B) = sum(ab -> sosdiff(ab...),zip(A,B))

# Hard thresholding
hard(x, beta) = abs(x) < beta ? zero(x) : x

# Partial objective
_obj(zlk,λ) = sum(z -> (abs(z) < sqrt(2λ) ? abs2(z)/2 : λ), zlk)

# Conversions between filter forms
_filtermatrix(hlist) =
    (hcat([vec(h) for h in hlist]...)::Matrix{eltype(first(hlist))},
     size(first(hlist)))
_filterlist(Hmatrix,R) = [reshape(h,map(n->n,R)) for h in eachcol(Hmatrix)]


text:
using SymbolicCodegen
using Test

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