text:
export GHType,
       GHInteger,
       GHNumber,
       GHPoint,
       GHBoolean

GHType(type, name, Nickname=name[1:1], Description=name, Access=:item, Default=0;
    nickname=Nickname, description=Description, access=Access, default=Default)=
    "$(type)|$(name)|$(nickname)|$(description)|$(access)|$(default)"
GHInteger(name, Nickname=name[1:1], Description=name, Access=:item, Default=0;
    nickname=Nickname, description=Description, access=Access, default=Default) =
    GHType("Integer", name, nickname, description, access, default)
GHNumber(name, Nickname=name[1:1], Description=name, Access=:item, Default=0;
    nickname=Nickname, description=Description, access=Access, default=Default) =
    GHType("Number", name, nickname, description, access, default)
GHPoint(name, Nickname=name[1:1], Description=name, Access=:item, Default=u0();
    nickname=Nickname, description=Description, access=Access, default=Default) =
    GHType("Point", name, nickname, description, access, default)
GHBoolean(name, Nickname=name[1:1], Description=name, Access=:item, Default=true;
    nickname=Nickname, description=Description, access=Access, default=Default) =
    GHType("Boolean", name, nickname, description, access, default)


text:
using DSP, JetPackDSP, Jets, Test

n1,n2,n3 = 80,70,10

@testset "JopFilter linearity test, Lowpass, dimension=$(length(N)), T=$(T)" for T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    A = JopFilter(JetSpace(T,N), Lowpass(T(0.25)), Butterworth(T,5))
    lhs, rhs = linearity_test(A)
    @test lhs ≈ rhs
    lhs, rhs = linearity_test(A')
    @test lhs ≈ rhs
end

@testset "JotOpFilter linearity test, Highpass, dimension=$(length(N)), T=$(T)" for T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    A = JopFilter(JetSpace(T,N), Highpass(T(0.25)), Butterworth(T,5))
    lhs, rhs = linearity_test(A)
    @test lhs ≈ rhs
    lhs, rhs = linearity_test(A')
    @test lhs ≈ rhs
end

@testset "JotOpFilter linearity test, Bandpass, dimension=$(length(N)), T=$(T)" for T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    A = JopFilter(JetSpace(T,N), Bandpass(T(0.2),T(0.3)), Butterworth(T,5))
    lhs, rhs = linearity_test(A)
    @test lhs ≈ rhs
    lhs, rhs = linearity_test(A')
    @test lhs ≈ rhs
end

@testset "JotOpFilter dot product test, Lowpass, dimension=$(length(N)), T=$(T)" for T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    A = JopFilter(JetSpace(T,N), Lowpass(T(0.1)), Butterworth(5))
    lhs,rhs = dot_product_test(A, -1 .+ 2 .* rand(domain(A)), -1 .+ 2 .* rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JotOpFilter dot product test, Highpass, dimension=$(length(N)), T=$(T)" for T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    A = JopFilter(JetSpace(T,N), Highpass(T(0.4)), Butterworth(5))
    lhs,rhs = dot_product_test(A, -1 .+ 2 .* rand(domain(A)), -1 .+ 2 .* rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JotOpFilter dot product test, Bandpass, dimension=$(length(N)), T=$(T)" for T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    A = JopFilter(JetSpace(T,N), Bandpass(T(0.1),T(0.4)), Butterworth(5))
    lhs,rhs = dot_product_test(A, -1 .+ 2 .* rand(domain(A)), -1 .+ 2 .* rand(range(A)))
    @test lhs ≈ rhs
end

ext: jl


text:
function encode(n)
  sort(split(string(n), ""))
end

n = BigInt(1)
while true
  global n
  if encode(n) == encode(2n) == encode(3n) == encode(4n) == encode(5n) == encode(6n)
    println(n)
    break
  end
  n += 1
end

max_stars_repo_name: pdbartlett/misc-stuff

max_issues_repo_name: pdbartlett/misc-stuff

max_forks_repo_name: pdbartlett/misc-stuff


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

function addext(filename::String, extension::String)::String
    filename = rstrip(filename, '.')
    extension = extension[1] == '.' ? extension : '.' * extension
    endswith(filename, extension) ? filename : filename * extension
end

function prepare_path(filename::String, extension::String, directory = pwd())::String
    # make sure the filename has the right extension
    filename = addext(filename, extension)
    # normalize path
    normpath(joinpath(directory, filename))
end

function skipStarting(stream::IOStream, chars::AbstractVector{Char})
    c = read(stream, Char)
    while c in chars
        readline(stream)
        c = read(stream, Char)
    end
    skip(stream, -ncodeunits(c))
end

module RCompat

function r2j_weights(rweights::AbstractMatrix{Float64}, rows::Int, columns::Int)
    result = Array{Float64, 3}(undef, size(rweights)[2], rows, columns)
    for r in 1:rows, c in 1:columns
        ind = j2r_ind(r,c,columns)
        result[:,r,c] = rweights[ind,:]
    end
    result
end

function j2r_weights(jweights::AbstractArray{Float64, 3})
    (n, rows, columns) = size(jweights)
    result = Array{Float64, 2}(undef, rows*columns, n)
    for r in 1:rows, c in 1:columns
        ind = j2r_ind(r,c,columns)
        result[ind,:] = jweights[:,r,c]
    end
    result
end

@inline j2r_ind(_, r, c, columns) = j2r_ind(r, c, columns)

function j2r_ind(i::CartesianIndex, columns)
    j2r_ind(i.I...,columns)
end

function j2r_ind(r, c, columns)
    (r-1) * columns  + c
end

function r2j_ind(i, col)
    (div(i-1, col) + 1,
     mod(i-1, col) + 1)
end

end # module

max_stars_repo_name: ckafi/DataIO

max_issues_repo_name: ckafi/DataIO

max_forks_repo_name: ckafi/DataIO


text:
include("./functions.jl")

using .Functions: softmax, cross_entropy_error

abstract type AbstractLayer end

mutable struct MatMul <: AbstractLayer
    params
    grads
    x
    function MatMul(W)
        layer = new()
        layer.params = [W]
        layer.grads = [zero(W)]
        return layer
    end
end

function forward!(layer::MatMul, x)
    W, = layer.params
    layer.x = x
    return out = x * W
end

# Keyword Argのデフォルトが指定されていない場合，必ず呼出側で引数指定する
function backward!(layer::MatMul, dout)
    W, = layer.params
    dx = dout * W'
    dW = layer.x' * dout
    layer.grads[1] = deepcopy(dW)
    return dx
end

mutable struct Affine <: AbstractLayer
    params
    grads
    x
    function Affine(W, b)
        layer = new()
        layer.params = [W, b]
        layer.grads = [zero(W), zero(b)]
        return layer
    end
end

function forward!(layer::Affine, x)
    W, b = layer.params
    layer.x = x
    return out = x * W .+ b
end

function backward!(layer::Affine; dout)
    W, b = layer.params
    dx = dout * W'
    dW = layer.x' * dout
    db = sum(dout, dims = 1)

    layer.grads[1] = dW
    layer.grads[2] = db
    return dx
end

mutable struct Softmax <: AbstractLayer
    params
    grads
    out
    function Softmax()
        layer = new()
        layer.params = []
        layer.grads = []
        return layer
    end
end

function forward!(layer::Softmax, x)
    return layer.out = softmax(x)
end

function backward!(layer::Softmax; dout)
    dx = layer.out .* dout
    sumdx = sum(dx, dims = 2)
    return dx -= layer.out .* sumdx
end

mutable struct SoftmaxWithLoss <: AbstractLayer
    params
    grads
    y # SoftmaxレイヤからCrossEntropyErrorレイヤへのデータ
    t
    function SoftmaxWithLoss()
        layer = new()
        layer.params = []
        layer.grads = []
        return layer
    end
end

function forward!(layer::SoftmaxWithLoss, x, t)
    layer.t = t
    layer.y = softmax(x)
    return loss = cross_entropy_error(layer.y, layer.t)
end

function backward!(layer::SoftmaxWithLoss; dout = 1.0)
    batch_size = size(layer.t, 1)
    return dx = (layer.y .- layer.t) .* dout ./ batch_size
end

mutable struct Sigmoid <: AbstractLayer
    params
    grads
    out
    function Sigmoid()
        layer = new()
        layer.params = []
        layer.grads = []
        return layer
    end
end

function forward!(layer::Sigmoid, x)
    return layer.out = 1 ./ (1 .+ exp.(-x))
end

function backward!(layer::Sigmoid; dout)
    return dx = dout .* (1.0 .- layer.out) .* layer.out
end

mutable struct Embedding <: AbstractLayer
    params
    grads
    idx
    function Embedding(W)
        layer = new()
        layer.params = [W]
        layer.grads = [zero(W)]
        layer.idx = nothing
        return layer
    end
end

function forward!(layer::Embedding, idx)
    W, = layer.params
    layer.idx = idx
    out = selectdim(W, 1, idx)
    return out
end

function backward!(layer::Embedding, dout)
    dW, = self.grads
    dW = deepcopy(zero(dW))
    selectdim(dW, 1, layer.idx) .+= dout
end

mutable struct SigmoidWithLoss <: AbstractLayer
    params
    grads
    loss
    y
    t
    function SigmoidWithLoss()
        self = new()
        self.loss = nothing
        self.y = nothing
        self.t = nothing
        return self
    end
end

function forward!(layer::SigmoidWithLoss, x, t)
    layer.t = t
    layer.y = 1 ./ (1 .+ exp.(-x))

    layer.loss = cross_entropy_error(hcat(1 - layer.y, layer.y), layer.t)
    return layer.loss
end

function backward!(layer::SigmoidWithLoss; dout = 1)
    batch_size = size(layer.t, 1)
    dx = (layer.y - layer.t) .* dout ./ batch_size
    return dx
end


text:
"""
    objective_function_template(nobs::Function,
                                obj_contribution::Function)

A constructor of objects of [composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) for defining an `objective_function_template`.

Arguments
===
+ `nobs`: a function of `data` that computes the number of observations of the particular data type,
+ `obj_contribution`: a function of the parameters `theta`, the `data` and the observation index `i` that returns a `Float64`.

Result
===
An `objective_function_template` object with fields `nobs` and `obj_contributions`.

"""
struct objective_function_template
    nobs::Function
    obj_contribution::Function
end

"""
    objective_function(theta::Vector,
                       data::Any,
                       template::objective_function_template,
                       br::Bool = false)

Evaluates the objective function at `theta` by adding up all contributions in
`data`, according to [`objective_function_template`](@ref).

Arguments
===
+ `theta`: a `Vector` of parameter values at which to evaluate the objective function
+ `data`: typically an object of [composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) with all the data required to compute the `objective_function`.
+ `template`: an [`objective_function_template`](@ref) object
+ `br`: a `Bool`. If `false` (default), the objective function is constructed by adding up all contributions in
`data`, according to [`objective_function_template`](@ref), before it is evaluated at `theta`. If `true` then the bias-reducing penalty in [Kosmidis & Lunardon, 2020](http://arxiv.org/abs/2001.03786) is computed and added to the objective function.

Result
===
A `Float64`.

Details
===
`data` can be used to pass additional constants other than the actual data to the objective.
"""
function objective_function(theta::Vector,
                            data::Any,
                            template::objective_function_template,
                            br::Bool = false)
    p = length(theta)
    n_obs = template.nobs(data)
    objective = 0
    for i in 1:n_obs
        objective += template.obj_contribution(theta, data, i)
    end
    if (br)
        objective + obj_quantities(theta, data, template, true)[1]
    else
        objective
    end
end

function obj_quantities(theta::Vector,
                        data::Any,
                        template::objective_function_template,
                        penalty::Bool = false)
    objective_i = (pars, i) -> template.obj_contribution(pars, data, i)
    results_i = (x, i) -> begin
        out = DiffResults.HessianResult(x)
        ForwardDiff.hessian!(out, pars -> objective_i(pars, i), x)
    end
    p = length(theta)
    n_obs = template.nobs(data)
    psi = zeros(p)
    emat = zeros(p, p)
    jmat = zeros(p, p)
    for i in 1:n_obs
        cdiffres = results_i(theta, i)
        cpsi = DiffResults.gradient(cdiffres)
        psi += cpsi
        emat += cpsi * cpsi'
        jmat += - DiffResults.hessian(cdiffres)
    end
    jmat_inv = try
        inv(jmat)
    catch
        fill(NaN, p, p)
    end
    vcov = jmat_inv * (emat * jmat_inv)
    if (penalty)
        br_penalty = - tr(jmat_inv * emat) / 2
        # br_penalty = n_obs * log(det(Matrix{Float64}(I * n_obs, p, p) -
        #                              jmat_inv * emat)) / 2
        # br_penalty = + log(det(sum(njmats))) / 2 - log(det(emat)) / 2
        [br_penalty, jmat_inv, emat, psi]
    else
        [vcov, jmat_inv, emat, psi]
    end
end

function obj_quantities_old(theta::Vector,
                            data::Any,
                            template::objective_function_template,
                            penalty::Bool = false)
    function gr_i(eta::Vector, i::Int)
        out = similar(eta)
        ForwardDiff.gradient!(out, beta -> template.obj_contribution(beta, data, i), eta)
    end
    function he_i(eta::Vector, i::Int)
        out = similar(eta, p, p)
        ForwardDiff.hessian!(out, beta -> template.obj_contribution(beta, data, i), eta)
    end
    p = length(theta)
    n_obs = template.nobs(data)
    psi = zeros(p)
    emat = zeros(p, p)
    jmat = zeros(p, p)
    for i in 1:n_obs
        cpsi = gr_i(theta, i)
        psi += cpsi
        emat += cpsi * cpsi'
        jmat += - he_i(theta, i)
    end
    jmat_inv = try
        inv(jmat)
    catch
        fill(NaN, p, p)
    end
    vcov = jmat_inv * (emat * jmat_inv)
    if (penalty)
        br_penalty = - tr(jmat_inv * emat) / 2
        # br_penalty = n_obs * log(det(Matrix{Float64}(I * n_obs, p, p) -
        #                              jmat_inv * emat)) / 2
        # br_penalty = + log(det(sum(njmats))) / 2 - log(det(emat)) / 2
        [br_penalty, jmat_inv, emat, psi]
    else
        [vcov, jmat_inv, emat, psi]
    end
end

"""
    estimating_function_template(object::objective_function_template)

Construct an [`estimating_function_template`](@ref) from an [`objective_function_template`](@ref).

Arguments
===
+ `object`: an [`objective_function_template`](@ref)

Details
===
The field `ef_contribution` of the result is computed by differentiating (using the [**ForwardDiff**](https://github.com/JuliaDiff/ForwardDiff.jl) package) `object.obj_contribution` with respect to `theta`.

Result
===
A [`estimating_function_template`](@ref) object.

"""
function estimating_function_template(object::objective_function_template)
    function ef_contribution(theta::Vector, data::Any, i::Int64)
        out = similar(theta)
        ForwardDiff.gradient!(out, b -> object.obj_contribution(b, data, i), theta)
    end
    estimating_function_template(object.nobs, ef_contribution)
end


text:
export graph_newton_schulz, graph_newton_schulz_degopt

"""
     (graph,crefs)=graph_newton_schulz(
        k,
        T = ComplexF64;
        input = :A,
        Z = :Z,
        Q = :Q,
        V = :V,
    )

Does `k` iterations of the Newton–Schulz iteration for approximating the inverse
of a matrix A (name given by `input`), i.e., the recursion

    V_k+1 = V_k*(2*I - A*V_k),

with `V_0=A`.
The recursion is implemented using the graph operations

    Z_i=A*V_i
    Q_i=2*I-Z_i
    V_{i+1}=V_i*Q_i,

and the kwargs `Z`, `Q`, and `V` determines the naming of the corresponding
intermediate steps.

**References**

1. G. Schulz. "Iterative Berechnung der reziproken Matrix". Zeitschrift für
   Angewandte Mathematik und Mechanik, 13:57–59, 1933.
   DOI: [10.1002/zamm.19330130111](https://doi.org/10.1002/zamm.19330130111)

2. N. J. Higham. "Functions of Matrices". SIAM, Philadelphia, PA, 2008.
   DOI: [10.1137/1.9780898717778](https://doi.org/10.1137/1.9780898717778)
"""
function graph_newton_schulz(
    k,
    T = ComplexF64;
    input = :A,
    Z = :Z,
    Q = :Q,
    V = :V,
)
    graph = Compgraph(T)
    cref = Vector{Tuple{Symbol,Int}}(undef, 0)

    for i = 1:k
        if (i == 1)
            Vi = input
        else
            Vi = Symbol("$(V)$i")
        end
        Zi = Symbol("$(Z)$i")
        add_mult!(graph, Zi, input, Vi)
        Qi = Symbol("$(Q)$i")
        add_lincomb!(graph, Qi, 2, :I, -1.0, Zi)
        Vii = Symbol("$(V)$(i+1)")
        add_mult!(graph, Vii, Vi, Qi)
        add_output!(graph, Vii)
    end
    Vii = Symbol("$(V)$(k+1)")

    return (graph, cref)
end

"""
     (graph,crefs)=graph_newton_schulz_degopt(k, T=ComplexF64; input=:A)

Does `k` iterations of the Newton–Schulz iteration for approximating the inverse,
using the recursion

    Z_i=A*V_i
    V_{i+1}=V_i*(2*I-Z_i).

The function makes a call to [`graph_degopt`](@ref), resulting in more
degrees of freedom in `crefs`. See also [`graph_newton_schulz`](@ref).
"""
function graph_newton_schulz_degopt(k, T = ComplexF64; input = :A)
    x = Vector{Tuple{Vector{T},Vector{T}}}(undef, 2 * k)
    # Z_i=A*V_i --- Odd numbers
    for i = 1:2:(2*k-1)
        x[i] =
            (vcat(zero(T), one(T), zeros(T, i - 1)), vcat(zeros(T, i), one(T)))
    end
    # V_{i+1}=V_i*(2*I-Z_i) --- Even numbers
    for i = 2:2:(2*k)
        x[i] = (
            vcat(zeros(T, i - 1), one(T), zero(T)),
            vcat(2, zeros(T, i - 1), -one(T)),
        )
    end

    return graph_degopt(x, vcat(zeros(T, 2 * k + 1), one(T)), input = input)
end


text:
using FujitsuBLAS
using Test
using LinearAlgebra

@testset "FujitsuBLAS.jl" begin
    @test !iszero(FujitsuBLAS.__init__())
    blases = BLAS.get_config().loaded_libs
    fjblas = findfirst(x -> contains(x.libname, r"libfjlapackexsve_ilp64.so$"), blases)
    @test !isnothing(fjblas)
    @test blases[fjblas].interface === :ilp64
end


text:
using Test

@test 1 ≈ 0.999999999999

@test 1 ≈ 0.99999

#           Test.@inferred [AllowedType] f(x)

f(a) = a > 1 ? 1 : 1.0

typeof(f(2))
typeof(f(1))
@code_warntype f(1)
@code_warntype f(2)

@inferred f(2)
@inferred f(1)

@inferred max(1,2)

g(a) = a < 10 ? missing : 1.0

g(2)
g(20)

@inferred g(2)
@inferred g(20)

@inferred Missing g(2)
@inferred Missing g(20)


text:
using OED
using Base.Test
using jInv.Mesh
using jInv.Utils

domain = [-1 1 -1 1.]; n = [8 9]
Mr     = getRegularMesh(domain,n)
xc     = .8*getCellCenteredGrid(Mr)

dataT     = exp.(-xc[:,1].^2- xc[:,2].^2)
dataT     = reshape(dataT,tuple(n...))

Img = randn(prod(n))
At,dAt = getLinearInterMatrixTranspose(domain,n,xc,doDerivative=true,Img=Img)
@test norm(sum(At,1)-1,Inf)<1e-12

function linearInterMatTranspTestFun(x,v=[])
	At,dAT = getLinearInterMatrixTranspose(domain,n,x,doDerivative=true,Img=Img)

	if isempty(v)
		return At*Img
	else
		return At*Img,dAT*v
	end
end
x0 = vec(xc+1e-2*randn(size(xc)))
chkDer, = checkDerivative(linearInterMatTranspTestFun,x0)
@test chkDer

x,y = getCellCenteredAxes(Mr)
A0 = getLinearInterMatrix(domain,n,copy(x0))
A0t= getLinearInterMatrixTranspose(domain,n,x0)
@test norm(A0-A0t',Inf) <= 1e-10


text:
using DSGE
using Test, BenchmarkTools
using JLD2

import DSGE: klein_transition_matrices, n_model_states, n_backward_looking_states

# What do you want to do?
check_steady_state = true
check_jacobian = true
check_solution = true
check_measurement = true

path = dirname(@__FILE__)

m = RealBond()

# Steady-state computation
if check_steady_state
    steadystate!(m)
    # @btime steadystate!(m)

    file = jldopen("$path/reference/steady_state.jld2", "r")
    saved_ell  = read(file, "ell")
    saved_c    = read(file, "c")
    saved_η    = read(file, "eta")
    saved_μ    = read(file, "mu")
    saved_β    = read(file, "beta")
    saved_χss  = read(file, "chi_ss")
    close(file)

    @testset "Check steady state outputs" begin
        @test saved_ell ≈ m[:lstar].value
        @test saved_c   ≈ m[:cstar].value
        @test saved_η   ≈ m[:ηstar].value
        @test saved_μ   ≈ m[:μstar].value
        @test saved_β   ≈ m[:βstar].value
        @test saved_χss ≈ m[:χstar].value
    end
end

# Jacobian computation
if check_jacobian
    m.testing = true        # So that it will test against the unnormalized Jacobian
    JJ = DSGE.jacobian(m)
    # @btime JJ = DSGE.jacobian(m)

    file = jldopen("$path/reference/jacobian.jld2", "r")
    saved_JJ  = read(file, "JJ")

    MUP = read(file, "MUP")
    ZP  = read(file, "ZP")
    MONP = read(file, "MONP")
    ELLP = read(file, "ELLP")
    RRP = read(file, "RRP")
    IIP = read(file, "IIP")
    WWP = read(file, "WWP")
    PIP = read(file, "PIP")
    TTP = read(file, "TTP")
    MU  = read(file, "MU")
    Z   = read(file, "Z")
    MON = read(file, "MON")
    ELL = read(file, "ELL")
    RR  = read(file, "RR")
    II  = read(file, "II")
    WW  = read(file, "WW")
    PI  = read(file, "PI")
    TT  = read(file, "TT")

    F1 = read(file, "F1")
    F2 = read(file, "F2")
    F3 = read(file, "F3")
    F4 = read(file, "F4")
    F5 = read(file, "F5")
    F6 = read(file, "F6")
    F7 = read(file, "F7")
    F8 = read(file, "F8")
    F9 = read(file, "F9")
    close(file)

    @testset "Check jacobian outputs" begin
        @testset "Euler Equation" begin
            @test saved_JJ[F1, ELLP] ≈ JJ[F1, ELLP]
            @test saved_JJ[F1, RRP]  ≈ JJ[F1, RRP]
            @test saved_JJ[F1, WWP]  ≈ JJ[F1, WWP]
            @test saved_JJ[F1, TTP]  ≈ JJ[F1, TTP]
            @test saved_JJ[F1, ELL]  ≈ JJ[F1, ELL]
            @test saved_JJ[F1, RR]   ≈ JJ[F1, RR]
            @test saved_JJ[F1, WW]   ≈ JJ[F1, WW]
            @test saved_JJ[F1, TT]   ≈ JJ[F1, TT]
        end

        @testset "Kolmogorov Forward Equation" begin
            @test saved_JJ[F2, MUP] ≈ JJ[F2, MUP]
            @test saved_JJ[F2, MU]  ≈ JJ[F2, MU]
            @test saved_JJ[F2, ELL] ≈ JJ[F2, ELL]
            @test saved_JJ[F2, WW]  ≈ JJ[F2, WW]
            @test saved_JJ[F2, RR]  ≈ JJ[F2, RR]
            @test saved_JJ[F2, TT]  ≈ JJ[F2, TT]
        end

        @testset "Market clearing condition" begin
            @test saved_JJ[F3, MU]  ≈ JJ[F3, MU]
            @test saved_JJ[F3, Z]   ≈ JJ[F3, Z]
            @test saved_JJ[F3, ELL] ≈ JJ[F3, ELL]
            @test saved_JJ[F3, WW]  ≈ JJ[F3, WW]
            @test saved_JJ[F3, RR]  ≈ JJ[F3, RR]
            @test saved_JJ[F3, TT]  ≈ JJ[F3, TT]
        end

        @testset "Technology process" begin
            @test saved_JJ[F4, ZP] ≈ JJ[F4, ZP]
            @test saved_JJ[F4, Z] ≈ JJ[F4, Z]
        end

        @testset "Phillips equation" begin
            @test saved_JJ[F5, PIP] ≈ JJ[F5, PIP]
            @test saved_JJ[F5, Z]   ≈ JJ[F5, Z]
            @test saved_JJ[F5, WW]  ≈ JJ[F5, WW]
            @test saved_JJ[F5, PI]  ≈ JJ[F5, PI]
        end

        @testset "Taylor equation" begin
            @test saved_JJ[F6, II]  ≈ JJ[F6, II]
            @test saved_JJ[F6, PI]  ≈ JJ[F6, PI]
            @test saved_JJ[F6, MON] ≈ JJ[F6, MON]
        end

        @testset "Fisher equation" begin
            @test saved_JJ[F7, II]  ≈ JJ[F7, II]
            @test saved_JJ[F7, PIP] ≈ JJ[F7, PIP]
            @test saved_JJ[F7, RR]  ≈ JJ[F7, RR]
        end

        @testset "Transfers equation" begin
            @test saved_JJ[F8, TT] ≈ JJ[F8, TT]
            @test saved_JJ[F8, Z]  ≈ JJ[F8, Z]
            @test saved_JJ[F8, WW] ≈ JJ[F8, WW]
        end
    end
end

# Solve
if check_solution
    m.testing = false      # So the Jacobian will be normalized within the klein solution
    gx, hx = klein(m)
    # @btime klein(m)

    file = jldopen("$path/reference/solve.jld2", "r")
    saved_gx   = read(file, "gx")
    saved_hx   = read(file, "hx")
    close(file)

    @testset "Check solve outputs" begin
        @test saved_gx  ≈ gx
        @test saved_hx  ≈ hx
    end
end

# Measurement
if check_measurement
    TTT_jump, TTT_state = klein(m)

    # Transition
    TTT, RRR = klein_transition_matrices(m, TTT_state, TTT_jump)
    CCC = zeros(n_model_states(m))

    transition_equation = Transition(TTT, RRR, CCC)

    # Measurement (needs the additional TTT_jump argument)
    measurement_equation = measurement(m, TTT, TTT_jump, RRR, CCC)

    @load "$path/reference/measurement_equation.jld2" obsmat

    @testset "Check measurement equation" begin
        @test obsmat ≈ measurement_equation[:ZZ][:, 1:n_backward_looking_states(m)]
    end
end


text:
using LightGraphs
using Base.Meta
import TikzGraphs
using DataStructures

type FlowGraph
    ast::CodeInfo
    statements::Vector{Any}
    in_edges::Vector{Vector{Int}}
end

function FlowGraph(ast::CodeInfo)

    # Build a list of statements and goto labels:
    statements = []
    targets = Dict{Int,Int}()
    for (i, stmt) in enumerate(ast.code)
        if isa(stmt, LineNumberNode) || isexpr(stmt, :inbounds)
            # skip these, they complicate flow analysis
            # (they may prevent detecting some jumps as break, return or continue)
            continue
        end
        if isexpr(stmt, :meta) && stmt.args[1] in (:pop_loc, :push_loc)
            # ditto
            continue
        end
        if isa(stmt, LabelNode)
            targets[stmt.label] = length(statements)+1
            continue
        end
        push!(statements, stmt)
    end

    in_edges = [Int[] for _ in eachindex(statements)]
    for (i, stmt) in enumerate(statements)
        if isa(stmt, GotoNode)
            target = targets[stmt.label]
            statements[i] = GotoNode(target)
            push!(in_edges[target], i)
            continue
        end
        if isexpr(stmt, :return)
            continue
        end
        if isexpr(stmt, :gotoifnot)
            test, label = stmt.args
            target = targets[label]
            statements[i] = Expr(:gotoifnot, test, target)
            push!(in_edges[target], i)
        end
        if i < length(statements)
            push!(in_edges[i+1], i)
        end
    end

    for (i, stmt) in enumerate(statements)
        println("$i $stmt")
    end

    FlowGraph(ast, statements, in_edges)
end

function showgraph(flow::FlowGraph)
    g = DiGraph(length(flow.statements))
    for i in eachindex(flow.statements)
        for j in flow.in_edges[i]
            add_edge!(g, j, i)
        end
    end
    display(TikzGraphs.plotHelper(g, "layered", "layered layout", ""))
end

"""
Note: this function does not handle the general case of converting
arbitrary programs with gotos to structured programs. It only knows
how to handle gotos generated by Julia from what was once a structured
program!
"""
function raise_flow(istart, flow::FlowGraph, seen=Set(Int[]), while_start=-1, while_end=-1)
    dest = []

    i = istart
    while 0 < i <= length(flow.statements)
        stmt = flow.statements[i]

        # Don't overrun the end of while loops:
        # ======================================
        if i == while_end+1
            break
        end

        # Check for flow interruptions:
        # ===================================
        if isexpr(stmt, :return)
            push!(seen, i)
            push!(dest, stmt)
            return dest, -1
        end
        if isa(stmt, GotoNode)
            target = stmt.label
            if target-1 == while_end
                push!(seen, i)
                push!(dest, Expr(:break))
                return dest, -1
            elseif target == while_start
                push!(seen, i)
                push!(dest, Expr(:continue))
                return dest, -1
            end
        end

        # Stop at reconvergence points until
        # all branches have been followed:
        # =======================================
        for j in flow.in_edges[i]
            if j < i
#                print("$istart-$i: checking $j...")
                if !(j in seen)
#                    println("WAIT")
                    return dest, i
#                else
#                    println("OK")
                end
            end
        end

        # OK, we're processing that node
#        println("$istart-$i: proceed")
        push!(seen, i)

        # Follow forward goto nodes:
        # ===============================
        if isa(stmt, GotoNode)
            target = stmt.label
            if target > i
                i = target
                continue
            end
            error("Unhandled backwards goto at $i")
        end

        # Detect while loops
        # ======================
        backwards = [j for j in flow.in_edges[i] if j > i]
        if length(backwards) > 1
            error("Too many backwards jump to $i")
        end

        # Handle gotoifnot nodes:
        # ============================
        if isexpr(stmt, :gotoifnot)
            test, else_target = stmt.args
#            println("$i: 'if' $test $else_target")

            if !isempty(backwards) && else_target == backwards[1] + 1
                # 'while $test' loop.
#                println("While loop from $i to $backwards "
#                        "$while_start $while_end")
                while_body, tail = raise_flow(i+1, flow, seen, i, backwards[1])
                if tail != -1 && tail != i
                    error("Bad 'while' reconvergence $i:$tail")
                end
                while_stmt = Expr(:while, test, Expr(:block, while_body...))
                push!(dest, while_stmt)
                i = backwards[1] + 1
                continue
            else
                # 'if then else' branch.
                # Visit both branches without letting each branch know
                # the nodes seen by the other branch:
                seen_if = Set(seen)
                if_branch, if_tail = raise_flow(i+1, flow, seen_if, while_start, while_end)

                seen_else = Set(seen)
                else_branch, else_tail = raise_flow(else_target, flow, seen_else, while_start, while_end)

#                println("$i: 'if' branches reached $if_tail/$else_tail")

                # Check for proper reconvergence:
                if else_tail == if_tail || else_tail == -1 || if_tail == -1
                    # Good reconvergence, mark the branches as seen:
                    push!(seen, seen_if...)
                    push!(seen, seen_else...)

                    # Build the if statement:
                    if_stmt = if length(else_branch) > 0
                        Expr(:if, test, Expr(:block, if_branch...),
                                        Expr(:block, else_branch...))
                    else
                        Expr(:if, test, Expr(:block, if_branch...))
                    end
                    push!(dest, if_stmt)

                    # Decide where to go next!
                    i = max(if_tail, else_tail)
                    continue
                else
                    error("Bad 'if' reconvergence $i:$if_tail/$else_tail")
                end
            end
        elseif !isempty(backwards)
            # 'while true' statement
            if i != while_start # otherwise we're already there
#                println("While loop from $i to $backwards "
#                        "$while_start $while_end")
                while_body, tail = raise_flow(i, flow, seen, i, backwards[1])
                if tail != -1 && tail != i
                    error("Bad 'while' reconvergence $i:$tail")
                end
                while_stmt = Expr(:while, :true, Expr(:block, while_body...))
                push!(dest, while_stmt)
                i = backwards[1] + 1
                continue
            end
        end

        # Just a regular statement:
#        push!(dest, LineNumberNode(i))
        push!(dest, stmt)
        i += 1
    end

    dest, -1
end


text:
function vocab(file)
    cut = 10
    counter = Dict{String, Int}()
    open(file, "r") do io
        for line in eachline(io)
            for word in split(line)
                counter[word] = 1 + get(counter, word, 0)
            end
        end
    end

    # Sort and display
    guys = collect(counter)
    sort!(guys, by=(x) -> x[2], rev=true)

    for x in guys
        if x[2] > cut
            println(x[1], '\t',  x[2])
        end
    end
end

vocab(ARGS[1])

max_stars_repo_name: alexalemi/textbench

max_issues_repo_name: alexalemi/textbench

max_forks_repo_name: alexalemi/textbench


text:
using CombUtils

function main()
    f1 = big(1)
    f2 = big(1)
    index = 2
    num_digits = 1
    while num_digits < 1000
        index += 1
        tmp = f2
        f2 += f1
        f1 = tmp
        num_digits = ndigits(f2)
    end
    @printf("The first index where the nth fib has 1000 digits is %d\n", index)
end

main()


text:
module ChooseOptimizer

using Cbc, JuMP

export set_solver, get_solver, get_solver_name
export set_solver_options, clear_solver_options, get_solver_options
export set_solver_verbose

_SOLVER = Cbc
_SOLVER_OPTS = Dict{String,Any}()

"""
`set_solver(OPT_NAME::Module=Cbc)` sets the optimization solver to be used.
Note: This automatically invokes `clear_solver_options()`.
"""
function set_solver(OPT_NAME::Module = Cbc)
    clear_solver_options()
    global _SOLVER = OPT_NAME
    # @info "Solver set to $_SOLVER"
    nothing
end

"""
`set_solver_options(kwd::Symbol, val)` set an option to be used
by the solver.

`set_solver_options(d::Dict)` adds all the options in `d`.
"""
function set_solver_options(kwd::String, val)
    global _SOLVER_OPTS[kwd] = val
end

function set_solver_options(d::Dict)
    for k in keys(d)
        set_solver_options(k, d[k])
    end
    nothing
end

"""
`clear_solver_options()` clears all solver options.
"""
function clear_solver_options()
    global _SOLVER_OPTS = Dict{String,Any}()
    nothing
end

"""
`get_solver_options()` returns the dictionary of current solver options.
"""
function get_solver_options()::Dict
    global _SOLVER_OPTS
    return _SOLVER_OPTS
end

"""
`get_solver()` generates the necessary arguments for creating a `JuMP`
`Model`. That is, one just uses `Model(get_solver())`.
"""
function get_solver()
    global _SOLVER
    global _SOLVER_OPTS

    d = get_solver_options()
    return JuMP.optimizer_with_attributes(_SOLVER.Optimizer, d...)

    # return JuMP.with_optimizer(_SOLVER.Optimizer; get_solver_options()...)
end

function get_solver_name()
    global _SOLVER
    return _SOLVER
end

"""
`set_solver_verbose(verb = true)` sets how verbose the solver is
during its work.
"""
function set_solver_verbose(verb::Bool = true)
    global _SOLVER
    val = verb ? 1 : 0
    key = "Unknown"
    # @info "Setting verbose option for $_SOLVER to $verb"
    try
        if _SOLVER == Cbc
            key = "LogLevel"
        end
    catch
    end

    try
        if _SOLVER == Main.Cbc
            key = "LogLevel"
        end
    catch
    end

    try
        if _SOLVER == Gurobi
            key = "OutputFlag"
        end
    catch
    end

    try
        if _SOLVER == Main.Gurobi
            key = "OutputFlag"
        end
    catch
    end

    try
        if _SOLVER == GLPK
            key = "msg_lev"
            val = verb ? 2 : 0 # GLPK uses 2 for "normal" messaging
        end
    catch
    end

    try
        if _SOLVER == Main.GLPK
            key = "msg_lev"
            val = verb ? 2 : 0 # GLPK uses 2 for "normal" messaging
        end
    catch
    end

    if key == "Unknown"
        @warn("Unable to set verbose option for $_SOLVER")
    else
        set_solver_options(key, val)
    end
    nothing
end

end


text:
export Segmenter

using Merlin

type Segmenter
    xdict::Dict{String,Int}
    ydict::Dict{String,Int}
    nn
end

function Segmenter()
    Segmenter(Dict{String,Int}(), Dict{String,Int}(), Model(4))
end

function (seg::Segmenter)(str::String)
    ids = map(collect(str)) do c
        id = seg.dict_x[string(c)]
        #id < 0 && throw("Invalid character: $(c).")
        id
    end
    x = reshape(ids, 1, length(ids))
    z = seg(Var(x))
    map(id -> getkey(seg.dict_y,id), z)
end

type Model
    g::Graph
end

function Model(outsize::Int)
    T = Float32
    x = Var([1,2])
    h = Lookup(T,300,10)(x)
    h = window(h, (70,), pads=(30,), strides=(10,))
    h = Linear(T,70,70)(h)
    h = tanh(h)
    h = Linear(T,70,outsize)(h)
    g = Graph([x], [h])
    Model(g)
end

function (m::Model)(x::Var, y=nothing)
    z = m.g(x)
    y == nothing ? argmax(z,1) : crossentropy(y,z)
end

function train(seg::Segmenter, trainpath::String, nepochs::Int)
    train_x, train_y = read(trainpath, seg.xdict, seg.ydict)
    n = floor(Int, length(train_x)*0.9)
    train_x, test_x = train_x[1:n], train_x[n+1:end]
    train_y, test_y = train_y[1:n], train_y[n+1:end]
    info("# training samples: $(length(train_x))")
    info("# test samples: $(length(test_x))")
    info("# x: $(length(seg.xdict))")
    info("# y: $(length(seg.ydict))")

    opt = SGD(0.0005)
    for epoch = 1:nepochs
        println("epoch: $epoch")
        loss = fit(train_x, train_y, seg.nn, opt)
        println("loss: $loss")

        ys = cat(1, map(y -> vec(y.data), test_y)...)
        zs = cat(1, map(x -> vec(seg.nn(x).data), test_x)...)
        acc = mean(i -> ys[i] == zs[i] ? 1.0 : 0.0, 1:length(ys))
        println("test accuracy: $acc")
        println()
    end
end

function read(path::String, xdict, ydict)
    data_x, data_y = Var[], Var[]
    buf_x, buf_y = Int[], Int[]
    lines = open(readlines, path)
    for i = 1:length(lines)
        line = chomp(lines[i])
        if isempty(line) || i == length(lines)
            push!(data_x, Var(reshape(buf_x,1,length(buf_x))))
            push!(data_y, Var(reshape(buf_y,1,length(buf_y))))
            buf_x = Int[]
            buf_y = Int[]
        else
            items = split(line, "\t")
            id = get!(xdict, items[1], length(xdict)+1)
            push!(buf_x, id)
            id = get!(ydict, items[2], length(ydict)+1)
            push!(buf_y, id)
        end
    end
    data_x, data_y
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
# SArray.jl

struct Values{N,T} <: TupleVector{N,T}
    v::NTuple{N,T}
    Values{N,T}(x::NTuple{N,T}) where {N,T} = new{N,T}(x)
    Values{N,T}(x::NTuple{N,Any}) where {N,T} = new{N,T}(convert_ntuple(T, x))
end

@pure @generated function (::Type{Values{N,T}})(x::Tuple) where {T, N}
    return quote
        @_inline_meta
        Values{N,T}(x)
    end
end

@inline Values(a::TupleVector{N}) where N = Values{N}(Tuple(a))
@propagate_inbounds Base.getindex(v::Values, i::Int) = v.v[i]
@inline Tuple(v::Values) = v.v
Base.dataids(::Values) = ()

# See #53
Base.cconvert(::Type{Ptr{T}}, a::Values) where {T} = Base.RefValue(a)
Base.unsafe_convert(::Type{Ptr{T}}, a::Base.RefValue{SA}) where {N,T,SA<:Values{N,T}} = Ptr{T}(Base.unsafe_convert(Ptr{Values{N,T}}, a))

# SVector.jl

@inline Values(x::NTuple{N,Any}) where N = Values{N}(x)
@inline Values{N}(x::NTuple{N,T}) where {N,T} = Values{N,T}(x)
@inline Values{N}(x::T) where {N,T<:Tuple} = Values{N,promote_tuple_eltype(T)}(x)

# Some more advanced constructor-like functions
@pure @inline Base.zeros(::Type{Values{N}}) where N = zeros(Values{N,Float64})
@pure @inline Base.ones(::Type{Values{N}}) where N = ones(Values{N,Float64})

# Converting a CartesianIndex to an SVector
Base.convert(::Type{Values}, I::CartesianIndex) = Values(I.I)
Base.convert(::Type{Values{N}}, I::CartesianIndex{N}) where {N} = Values{N}(I.I)
Base.convert(::Type{Values{N,T}}, I::CartesianIndex{N}) where {N,T} = Values{N,T}(I.I)

@pure Base.promote_rule(::Type{Values{N,T}}, ::Type{CartesianIndex{N}}) where {N,T} = Values{N,promote_type(T,Int)}


text:
using Base.Test
using QuantumOptics

shape1 = [5]
shape2 = [2, 3]

b1 = GenericBasis(shape1)
b2 = GenericBasis(shape2)

@test b1.shape == shape1
@test b2.shape == shape2
@test b1 != b2
@test b1 == b1

comp_b1 = tensor(b1, b2)
comp_b2 = tensor(b1, b1, b2)
@test comp_b1.shape == [prod(shape1), prod(shape2)]
@test comp_b2.shape == [prod(shape1), prod(shape1), prod(shape2)]

comp_b1_b2 = tensor(comp_b1, comp_b2)
@test comp_b1_b2.shape == [prod(shape1), prod(shape2), prod(shape1), prod(shape1), prod(shape2)]
@test comp_b1_b2 == CompositeBasis(b1, b2, b1, b1, b2)

@test comp_b2.shape == tensor(b1, comp_b1).shape
@test comp_b2 == tensor(b1, comp_b1)

@test ptrace(comp_b2, [1]) == ptrace(comp_b2, [2]) == comp_b1
@test ptrace(comp_b2, [1, 2]) == ptrace(comp_b1, [1])
@test ptrace(comp_b2, [2, 3]) == ptrace(comp_b1, [2])


text:
using MPI
using ClimateMachine
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes

using Dates
using ClimateMachine.Diagnostics
using CLIMAParameters

macro oldnew()
 return eval(newstyle)
end

import  ClimateMachine.SystemSolvers:
    BatchedGeneralizedMinimalResidual, linearsolve!

# Needed to control diagnostics dispatch
using ClimateMachine.ConfigTypes
struct IVDCConfigType <: ClimateMachineConfigType end

include("MYDIAGSModel.jl")
include("IVDCModel.jl")
include("ivdc_diagnostics.jl")

 #####
 # Basic initialization
 #####
 ClimateMachine.init()
 ArrayType = ClimateMachine.array_type()
 mpicomm = MPI.COMM_WORLD
 FT = Float64

 #####
 # Create mesh
 #####
 const N = 4
 # const Nˣ = 20
 # const Nʸ = 20
 const Nˣ = 2
 const Nʸ = 2
 const Nᶻ = 20
 const Lˣ = 4e6  # m
 const Lʸ = 4e6  # m
 const H = 1000  # m

 xrange = range(FT(0); length = Nˣ + 1, stop = Lˣ)
 yrange = range(FT(0); length = Nʸ + 1, stop = Lʸ)
 zrange = range(FT(-H); length = Nᶻ + 1, stop = 0)

 gconf  = (N=N,Nˣ=Nˣ,Nʸ=Nʸ,Lˣ=Lˣ,Lʸ=Lʸ,H=H,xrange=xrange,yrange=yrange,zrange=zrange)

 brickrange_2D = (xrange, yrange)
 topl_2D       =
   BrickTopology( mpicomm, brickrange_2D, periodicity = (false, false) );
 grid_2D = DiscontinuousSpectralElementGrid(
   topl_2D,
   FloatType = FT,
   DeviceArray = ArrayType,
   polynomialorder = N,
 )
 brickrange_3D = (xrange, yrange, zrange)
 topl_3D = StackedBrickTopology(
   mpicomm,
   brickrange_3D;
   periodicity = (false, false, false),
   boundary = ((1, 1), (1, 1), (2, 3)),
 )
 grid_3D = DiscontinuousSpectralElementGrid(
   topl_3D,
   FloatType = FT,
   DeviceArray = ArrayType,
   polynomialorder = N,
 )

 # Wave speeds for numerical flux
 const cʰ = 1  # typical of ocean internal-wave speed
 const cᶻ = 0

 # Set a timestep for implicit solve
 dt = 5400/5

 # Create a container balance law for working with multiple DG models in Diagnostics

 # Create balance law and RHS arrays for diffusion equation
 ivdc_dg = IVDCDGModel(
  IVDCModel{FT}(;dt=dt,cʰ=cʰ,cᶻ=cᶻ,gconf=gconf),
  grid_3D,
  RusanovNumericalFlux(),
  CentralNumericalFluxSecondOrder(),
  CentralNumericalFluxGradient();
  direction=VerticalDirection(),
 );

 ivdc_Q   = init_ode_state(ivdc_dg, FT(0); init_on_cpu=true)
 ivdc_RHS = init_ode_state(ivdc_dg, FT(0); init_on_cpu=true)

 # Instantiate a batched GM-RES solver uaing balance law as its operator
 ivdc_bgm_solver=BatchedGeneralizedMinimalResidual(
   ivdc_dg,
   ivdc_Q;
   max_subspace_size=10);

 # Save initial state
 θ₀=deepcopy( ivdc_Q.θ )

 # Setup some callbacks and diagnostics
 diagsconf=(myivdcdg=ivdc_dg,myivdcq=ivdc_Q)
 mydiags_dg=MYDIAGSModel(
  MYDIAGSModel{FT}(;diagsconf=diagsconf),
  grid_3D,
  RusanovNumericalFlux(),
  CentralNumericalFluxSecondOrder(),
  CentralNumericalFluxGradient();
  direction=VerticalDirection(),
 );

 callbacks = ()
 struct EarthParameterSet <: AbstractEarthParameterSet end
 pset = EarthParameterSet()
 dgn_starttime = replace(string(now()), ":" => ".")
 od=ClimateMachine.ClimateMachine_Settings().output_dir
 Diagnostics.init(
   mpicomm,
   pset,
   mydiags_dg,
   ivdc_Q,
   dgn_starttime,
   od
  )

 # dgn_parms=(interval="1steps",type=IVDCConfigType(),name="default")
 dgngrp=setup_single_column_default_diagnostics( IVDCConfigType(), "1steps", "default_output_group")
 # dgngrp=setup_single_column_default_diagnostics(dgn_parms.type,dgn_parms.interval,dgn_parms.name)
 # dgn_config=ClimateMachine.DiagnosticsConfiguration([dgngrp])

 # Set up right hand side
 # for i=1:200000
 for i=1:50
  ivdc_RHS.θ   .= ivdc_Q.θ/dt
  #### ivdc_RHS.θ   .= ivdc_Q.θ
  ivdc_dg.state_auxiliary.θ_init .= ivdc_Q.θ

  println("Before maximum ", maximum(ivdc_Q.θ) )
  println("Before minimum ", minimum(ivdc_Q.θ) )

  # Evaluate operator
  #### ivdc_dg(ivdc_Q,ivdc_RHS,nothing,0;increment=false);
  #### println( maximum(ivdc_Q) )

  # Now try applying batched GM res solver
  lm!(y,x)=ivdc_dg(y,x,nothing,0;increment=false)
  solve_time = @elapsed iters = linearsolve!(lm!, ivdc_bgm_solver, ivdc_Q, ivdc_RHS);
  println("solver iters, time: ",iters, ", ", solve_time)

  println("After maximum ", maximum(ivdc_Q.θ) )
  println("After minimum ", minimum(ivdc_Q.θ) )
 end

 ## using Plots
 ## ni=2;ei=1
 ## nj=3;ej=1

 ## θᶠ=reshape( ivdc_Q.θ,(N+1,N+1,N+1,1,Nᶻ,Nˣ,Nʸ) )[ni,nj,:,1,:,ei,ej]
 ## θ⁰=reshape( θ₀,(N+1,N+1,N+1,1,Nᶻ,Nˣ,Nʸ) )[ni,nj,:,1,:,ei,ej]
 ## zc=reshape( grid_3D.vgeo,(N+1,N+1,N+1,16,Nᶻ,Nˣ,Nʸ) )[ni,nj,:,15,:,ei,ej]
 ## plot(θᶠ[:],zc[:],label="")
 ## plot!(θ⁰[:],zc[:],label="")

max_stars_repo_name: christophernhill/clima-impldiff-testing

max_issues_repo_name: christophernhill/clima-impldiff-testing

max_forks_repo_name: christophernhill/clima-impldiff-testing


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

mutable struct ResourceNavigationLinkFormat2 <: SwaggerModel
    linkedResourceType::Any # spec type: Union{ Nothing, String } # spec name: linkedResourceType
    link::Any # spec type: Union{ Nothing, String } # spec name: link
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ResourceNavigationLinkFormat2(;linkedResourceType=nothing, link=nothing, provisioningState=nothing)
        o = new()
        validate_property(ResourceNavigationLinkFormat2, Symbol("linkedResourceType"), linkedResourceType)
        setfield!(o, Symbol("linkedResourceType"), linkedResourceType)
        validate_property(ResourceNavigationLinkFormat2, Symbol("link"), link)
        setfield!(o, Symbol("link"), link)
        validate_property(ResourceNavigationLinkFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ResourceNavigationLinkFormat2

const _property_map_ResourceNavigationLinkFormat2 = Dict{Symbol,Symbol}(Symbol("linkedResourceType")=>Symbol("linkedResourceType"), Symbol("link")=>Symbol("link"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ResourceNavigationLinkFormat2 = Dict{Symbol,String}(Symbol("linkedResourceType")=>"String", Symbol("link")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ResourceNavigationLinkFormat2 }) = collect(keys(_property_map_ResourceNavigationLinkFormat2))
Swagger.property_type(::Type{ ResourceNavigationLinkFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ResourceNavigationLinkFormat2[name]))}
Swagger.field_name(::Type{ ResourceNavigationLinkFormat2 }, property_name::Symbol) =  _property_map_ResourceNavigationLinkFormat2[property_name]

function check_required(o::ResourceNavigationLinkFormat2)
    true
end

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


text:
function randArray()::Array{Int}
	return rand(1:1000, rand(1:10, rand(1:4))...);
end

function offsetArray()::OffsetArray{Int}
	return OffsetArray(1:10, -4);
end

function randStrArray()::Array{Int}
	return rand(1:100, rand(1:4, rand(1:3))...);
end

function lipsum()::String
	return "Lorem ipsum dolor sit amet, consectetur adipiscing elit. Vivamus dignissim, lacus et aliquet pellentesque, nibh turpis suscipit justo, non gravida massa lorem eu urna. Curabitur suscipit sapien vitae odio accumsan, quis mollis lectus tincidunt. Vestibulum tincidunt hendrerit nisl ac mattis. Praesent tristique magna quis tortor porta egestas. Morbi orci purus, posuere sit amet magna ac, sagittis pretium ipsum. Nam dictum turpis eu vulputate accumsan. Aliquam scelerisque ante sit amet justo facilisis, varius ornare orci posuere. In eu quam et ex sagittis fermentum et eu mauris. Sed porttitor massa in aliquam bibendum. Nam eleifend sed arcu a lobortis. Pellentesque semper, sapien non auctor consectetur, est nunc sollicitudin dolor, quis imperdiet dui turpis ac dui. Praesent nec vestibulum sem, ut placerat arcu. Quisque rhoncus placerat turpis, nec dictum mauris pulvinar vitae. Sed consectetur sem et ante tempor rutrum id et arcu. Maecenas sit amet finibus enim. Sed maximus mollis mauris in pretium. Suspendisse vehicula dui convallis aliquet mollis. Nunc justo felis, dictum ut efficitur at, hendrerit in est. Vestibulum faucibus quam mi, eget vulputate nisi iaculis ac. Nunc sem purus, convallis quis interdum vel, laoreet non orci. Praesent sed scelerisque mi. Aliquam sit amet lectus eu odio tempus tristique in et augue. Aliquam ornare augue ac tellus bibendum luctus. Morbi et ultricies magna. Pellentesque id libero tortor. Aenean bibendum dui nec bibendum suscipit. Maecenas in odio sodales, condimentum elit sit amet, euismod nisl. Nunc varius nulla sed elit posuere, vel fringilla mauris dictum. Nulla id condimentum lorem, eget efficitur odio. Nunc fermentum eros a sollicitudin laoreet. Maecenas suscipit, ipsum tempor lacinia aliquam, urna dui sagittis nisi, vel semper sapien orci quis nunc. Proin a rhoncus nunc, dapibus viverra dolor. Aliquam erat volutpat. Sed aliquam leo in sem mattis lobortis. Morbi vestibulum lectus nibh, sit amet gravida velit efficitur nec. Nam eget ligula sit amet neque mattis malesuada cras amet.";
end

function rand2Arrays()::Tuple{Array{Int}, Array{Int}}
	dims = rand(1:10, rand(1:4));
	return rand(1:1000, dims...), rand(1:1000, dims...);
end

function rand2StrArrays()::Tuple{Array{Int}, Array{Int}}
	dims = rand(1:4, rand(1:3));
	return rand(1:100, dims...), rand(1:100, dims...);
end

function randInit()::Int
	return rand(1:1000);
end


text:
# Setup functions for MathGL.jl

# resetting all parameters, clear graph
function reset(gr::Graph)
    cmgl.set_def_param(gr.ptr)
end

# Transparency options
function alpha(gr::Graph, enable::Bool)
    cmgl.set_alpha(gr.ptr, enable)
end

function alpha(gr::Graph, val::Real)
    cmgl.set_alpha_default(gr.ptr, val)
end

function transptype(gr::Graph, typ::Int)
    cmgl.set_transp_type(gr.ptr, typ)
end

# Lightning options
function light(gr::Graph, enable::Bool)
    cmgl.set_light(gr.ptr, enable)
end

function light(gr::Graph, n::Int, enable::Bool)
    cmgl.set_light_n(gr.ptr, n, enable)
end

function light(gr::Graph, n::Int, dx::Real, dy::Real, dz::Real;
               c::String="w", bright::Real=0.5, ap::Real=0.)
    cmgl.add_light_ext(gr.ptr, n, dx, dy, dz, c, bright, ap)
end

function light(gr::Graph, n::Int, rx::Real, ry::Real, rz::Real, dx::Real, dy::Real, dz::Real;
               c::String="w", bright::Real=0.5, ap::Real=0.)
    cmgl.add_light_loc(gr.ptr, n, rx, ry, rz, dx, dy, dz, c, bright, ap)
end

function diffuse(gr::Graph, bright::Real)
    cmgl.set_difbr(gr.ptr, bright)
end

function ambient(gr::Graph, bright::Real)
    cmgl.set_ambbr(gr.ptr, bright)
end

# Fog
function fog(gr::Graph, d::Real; dz::Real=0.25)
    cmgl.set_fog(gr.ptr, d, dz)
end

# Default widths, sizes, numbers, names...
function barwidth(gr::Graph, val::Real)
    cmgl.set_bar_width(gr.ptr, val)
end

function marksize(gr::Graph, val::Real)
    cmgl.set_mark_size(gr.ptr, val)
end

function arrowsize(gr::Graph, val::Real)
    cmgl.set_arrow_size(gr.ptr, val)
end

function meshnum(gr::Graph, n::Int)
    cmgl.set_meshnum(gr.ptr, n)
end

function facenum(gr::Graph, n::Int)
    cmgl.set_facenum(gr.ptr, n)
end

function plotid(gr::Graph, id::String)
    cmgl.set_plotid(gr.ptr, id)
end

function plotid(gr::Graph)
    bytestring(cmgl.get_plotid(gr.ptr))
end

# Cutting (excluding) points from drawing:
function cut(gr::Graph, enable::Bool)
    cmgl.set_cut(gr.ptr, enable)
end

function cut(gr::Graph, x1::Real, y1::Real, z1::Real, x2::Real, y2::Real, z2::Real)
    cmgl.set_cut_box(gr.ptr, x1, y1, z1, x2, y2, z2)
end

function cut(gr::Graph, cond::String)
    cmgl.set_cutoff(gr.ptr, enable)
end

# Text and font settings
function rotatetext(gr::Graph, enable::Bool)
    cmgl.set_rotated_text(gr.ptr, enable)
end

function loadfont(gr::Graph, name::String, path::String="")
    cmgl.load_font(gr.ptr, name, path)
end

function font(gr::Graph, fnt::String)
    cmgl.set_font_def(gr.ptr, fnt)
end

function fontsize(gr::Graph, s::Real, dpi::Int; unit=:pt)
    # Found the conversions in base.h
    if str(unit) == "pt"
        cmgl.set_font_size(s*27/dpi)
    elseif str(unit) == "cm"
        cmgl.set_font_size(s*27*28.45/dpi)
    elseif str(unit) == "in"
        cmgl.set_font_size(s*27*72.27/dpi)
    elseif str(unit) == "mm"
        cmgl.set_font_size(s*27*2.845/dpi)
    else
        cmgl.set_font_size(s)
    end
end

function font(gr::Graph, fnt::String, val::Real)
    cmgl.set_font_def(gr.ptr, fnt)
    cmgl.set_font_size(gr.ptr, val)
end

function font(gr::Graph, from::Graph)
    cmgl.copy_font(gr.ptr, from.ptr)
end

function restorefont(gr::Graph)
    cmgl.restore_font(gr.ptr)
end

function defaultfont(name::String, path::String="")
    cmgl.def_font(name, path)
end

# Color palette, colors, ...
function palette(gr::Graph, colors::String)
    cmgl.set_palette(gr.ptr, colors)
end

function scheme(gr::Graph, sch::String)
    cmgl.set_def_sch(gr.ptr, sch)
end

function colorid(id::Char, r::Real, g::Real, b::Real)
    cmgl.set_color(id, r, g, b)
end

# Masks
function mask(id::Char, hex::String)
    cmgl.set_mask(id, hex)
end

function mask(id::Char, hex::Integer)
    cmgl.set_mask(id, hex)
end

function mask(gr::Graph, angle::Integer)
    cmgl.set_mask_angle(gr.ptr, angle)
end

# TODO: Error handling
# TODO: Stop drawing
function stop(gr::Graph, enable::Bool=true)
    cmgl.ask_stop(gr.ptr, enable)
end

# Ranges (partially done) TODO: allow ranges to be set according to data
function xrange(gr::Graph, x1::Real, x2::Real; add::Bool=false)
    if add
        cmgl.add_range_val(gr.ptr, 'x', x1, x2)
    else
        cmgl.set_range_val(gr.ptr, 'x', x1, x2)
    end
end

function xrange(gr::Graph, x::Array{cmgl.Float})
    cmgl.set_range_val(gr.ptr, 'x', minimum(x), maximum(x))
end

function yrange(gr::Graph, y1::Real, y2::Real; add::Bool=false)
    if add
        cmgl.add_range_val(gr.ptr, 'y', y1, y2)
    else
        cmgl.set_range_val(gr.ptr, 'y', y1, y2)
    end
end

function yrange(gr::Graph, y::Array{cmgl.Float})
    cmgl.set_range_val(gr.ptr, 'y', minimum(y), maximum(y))
end

function zrange(gr::Graph, z1::Real, z2::Real; add::Bool=false)
    if add
        cmgl.add_range_val(gr.ptr, 'z', z1, z2)
    else
        cmgl.set_range_val(gr.ptr, 'z', z1, z2)
    end
end

function zrange(gr::Graph, z::Array{cmgl.Float})
    cmgl.set_range_val(gr.ptr, 'z', minimum(z), maximum(z))
end

function crange(gr::Graph, c1::Real, c2::Real; add::Bool=false)
    if add
        cmgl.add_range_val(gr.ptr, 'c', c1, c2)
    else
        cmgl.set_range_val(gr.ptr, 'c', c1, c2)
    end
end

function crange(gr::Graph, c::Array{cmgl.Float})
    cmgl.set_range_val(gr.ptr, 'c', minimum(c), maximum(c))
end

function ranges(gr::Graph, x1::Real, x2::Real, y1::Real, y2::Real, z1::Real=0., z2::Real=0.)
    cmgl.set_ranges(gr.ptr, x1, x2, y1, y2, z1, z2)
end

function ranges(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}=[0.])
    cmgl.set_ranges(gr.ptr, minimum(x), maximum(x), minimum(y), maximum(y), minimum(z), maximum(z))
end

function origin(gr::Graph, x0::Real, y0::Real, z0::Real=NaN)
    cmgl.set_origin(gr.ptr, x0, y0, z0)
end

function zoomaxis(gr::Graph, x1::Real, y1::Real, z1::Real, x2::Real, y2::Real, z2::Real)
    cmgl.zoom_axis(gr.ptr, x1, y1, z1, x2, y2, z2)
end

# Curvilinear coordinates
function coords(gr::Graph, how::Integer=0; X::String="", Y::String="", Z::String="", C::String="")
    if X == Y == Z == C == ""
        cmgl.set_coor(gr.ptr, how)
    else
        cmgl.set_func(gr.ptr, X, Y, Z, C)
    end
end

function ternary(gr::Graph, tern::Integer)
    cmgl.set_ternary(gr.ptr, tern)
end

# Ticks TODO: many options for ticks
function adjust(gr::Graph, dir::String="xyzc")
    cmgl.adjust_ticks(gr.ptr, dir)
end

function xtick(gr::Graph, d::Real=0; ns::Int=0, org::Real=NaN, fact::String="")
    cmgl.set_ticks_fact(gr.ptr, 'x', d, ns, org, fact)
end

function ytick(gr::Graph, d::Real=0; ns::Int=0, org::Real=NaN, fact::String="")
    cmgl.set_ticks_fact(gr.ptr, 'y', d, ns, org, fact)
end

function ztick(gr::Graph, d::Real=0; ns::Int=0, org::Real=NaN, fact::String="")
    cmgl.set_ticks_fact(gr.ptr, 'z', d, ns, org, fact)
end

function ctick(gr::Graph, d::Real=0; ns::Int=0, org::Real=NaN, fact::String="")
    cmgl.set_ticks_fact(gr.ptr, 'c', d, ns, org, fact)
end

function xtick(gr::Graph, lbl::Array{String}; add=false)
    cmgl.set_ticks_str(gr.ptr, 'x', reduce(*, Base.map(x->x*"\n", lbl)), add)
end
function ytick(gr::Graph, lbl::Array{String}; add=false)
    cmgl.set_ticks_str(gr.ptr, 'y', reduce(*, Base.map(x->x*"\n", lbl)), add)
end
function ztick(gr::Graph, lbl::Array{String}; add=false)
    cmgl.set_ticks_str(gr.ptr, 'z', reduce(*, Base.map(x->x*"\n", lbl)), add)
end
function ctick(gr::Graph, lbl::Array{String}; add=false)
    cmgl.set_ticks_str(gr.ptr, 'c', reduce(*, Base.map(x->x*"\n", lbl)), add)
end

function xtick(gr::Graph, val::Array{cmgl.Float}, lbl::Array{String}; add=false)
    cmgl.set_ticks_val(gr.ptr, 'x', Data(val).ptr, reduce(*, Base.map(x->x*"\n", lbl)), add)
end
function ytick(gr::Graph, val::Array{cmgl.Float}, lbl::Array{String}; add=false)
    cmgl.set_ticks_val(gr.ptr, 'y', Data(val).ptr, reduce(*, Base.map(x->x*"\n", lbl)), add)
end
function ztick(gr::Graph, val::Array{cmgl.Float}, lbl::Array{String}; add=false)
    cmgl.set_ticks_val(gr.ptr, 'y', Data(val).ptr, reduce(*, Base.map(x->x*"\n", lbl)), add)
end
function ctick(gr::Graph, val::Array{cmgl.Float}, lbl::Array{String}; add=false)
    cmgl.set_ticks_val(gr.ptr, 'c', Data(val).ptr, reduce(*, Base.map(x->x*"\n", lbl)), add)
end

# For Resetting <- would trigger overwrite-warning
#=function xtick(gr::Graph)=#
    #=cmgl.set_ticks(gr.ptr, 'x', 0, 0, NaN)=#
#=end=#
#=function ytick(gr::Graph)=#
    #=cmgl.set_ticks(gr.ptr, 'y', 0, 0, NaN)=#
#=end=#
#=function ztick(gr::Graph)=#
    #=cmgl.set_ticks(gr.ptr, 'z', 0, 0, NaN)=#
#=end=#

function xtick(gr::Graph, val::Real, lbl::String)
    cmgl.add_tick(gr.ptr, 'x', val, lbl)
end
function ytick(gr::Graph, val::Real, lbl::String)
    cmgl.add_tick(gr.ptr, 'y', val, lbl)
end
function ztick(gr::Graph, val::Real, lbl::String)
    cmgl.add_tick(gr.ptr, 'z', val, lbl)
end
function ctick(gr::Graph, val::Real, lbl::String)
    cmgl.add_tick(gr.ptr, 'c', val, lbl)
end

function xtick(gr::Graph, templ::String)
    cmgl.set_tick_templ(gr.ptr, 'x', templ)
end
function ytick(gr::Graph, templ::String)
    cmgl.set_tick_templ(gr.ptr, 'y', templ)
end
function ztick(gr::Graph, templ::String)
    cmgl.set_tick_templ(gr.ptr, 'z', templ)
end
function ctick(gr::Graph, templ::String)
    cmgl.set_tick_templ(gr.ptr, 'c', templ)
end

function tuneticks(gr::Graph, tune::Int, pos::Real)
    cmgl.tune_ticks(gr.ptr, tune, pos)
end

function tickshift(gr::Graph, dx::Real, dy::Real, dz::Real, dc::Real)
    cmgl.tickshift(gr.ptr, dx, dy, dz, dc)
end
function tickshift(gr::Graph; dx::Real=0., dy::Real=0., dz::Real=0., dc::Real=0.)
    cmgl.tickshift(gr.ptr, dx, dy, dz, dc)
end

function tickrotate(gr::Graph, val::Bool=true)
    cmgl.set_tick_rotate(gr.ptr, val)
end

function tickskip(gr::Graph, val::Bool=true)
    cmgl.set_tick_skip(gr.ptr, val)
end

#TODO: time formats for tics etc.

function origintick(gr::Graph, val::Bool=true)
    cmgl.set_flag(gr.ptr, val, cmgl.MGL_NO_ORIGIN)
end

function ticklen(gr::Graph, val::Real, stt::Real=1.)
    cmgl.set_tick_len(gr.ptr, val, stt)
end

function axisstl(gr::Graph, stl::String="", tck::String="", sub::String="")
    cmgl.set_axis_stl(gr.ptr, stl, tck, sub)
end

# Subplots and Rotations
# Preferred Order: subplot,multiplot,inplot -> title -> rotate,aspect -> plotting
# or:              columnplot, gridplot, stickplot, relative inplot -> ...

function subplot(gr::Graph, nx::Int, ny::Int, m::Int; stl::String="<>_^", dx::Real=0, dy::Real=0)
    cmgl.subplot_d(gr.ptr, nx, ny, m-1, stl, dx, dy)
end

function multiplot(gr::Graph, nx::Int, ny::Int, m::Int, dx::Int, dy::Int; stl::String="<>_^")
    cmgl.multiplot(gr.ptr, nx, ny, m-1, dx, dy, stl)
end

function inplot(gr::Graph, x1::Real, x2::Real, y1::Real, y2::Real, rel::Bool=true)
    if rel
        cmgl.inplot(gr.ptr, x1, x2, y1, y2)
    else
        cmgl.relplot(gr.ptr, x1, x2, y1, y2)
    end
end

function columnplot(gr::Graph, num::Int, ind::Int, d::Real=0.)
    cmgl.columnplot_d(gr.ptr, num, ind-1, d)
end

function gridplot(gr::Graph, nx::Int, ny::Int, ind::Int, d::Real=0.)
    cmgl.gridplot_d(gr.ptr, nx, ny, ind-1, d)
end

function stickplot(gr::Graph, num::Int, ind::Int, tet::Real, phi::Real)
    cmgl.stickplot(gr.ptr, num, ind-1, tet, phi)
end

function title(gr::Graph, txt::String; stl::String="", size::Real=-2)
    cmgl.title(gr.ptr, txt, stl, size)
end
#=function title(gr::Graph, txt::UTF8String; stl::String="", size::Int=-2)=#
    #=cmgl.titlew(gr.ptr, txt, stl, size)=#
#=end=#

function rotate(gr::Graph, tetx::Real, tetz::Real, tety::Real=0.)
    cmgl.rotate(gr.ptr, tetx, tetz, tety)
end

function rotate(gr::Graph, tet::Real, x::Real, y::Real, z::Real)
    cmgl.rotate_vector(gr.ptr, tet, x, y, z)
end

function xrotate(gr::Graph, tet::Real)
    cmgl.rotate_vector(gr.ptr, tet, 1., 0., 0.)
end
function yrotate(gr::Graph, tet::Real)
    cmgl.rotate_vector(gr.ptr, tet, 0., 1., 0.)
end
function zrotate(gr::Graph, tet::Real)
    cmgl.rotate_vector(gr.ptr, tet, 0., 0., 1.)
end

function aspect(gr::Graph, ax::Real, ay::Real, az::Real=1.)
    cmgl.aspect(gr.ptr, ax, ay, az)
end

#TODO Push and Pop for matrices on stack
function plotfactor(gr::Graph, val::Real)
    cmgl.plotfactor(gr.ptr, val)
end

function perspective(gr::Graph, a::Real)
    cmgl.perspective(gr.ptr, a)
end

function view(gr::Graph, tetx::Real, tetz::Real, tety::Real=0.)
    cmgl.view(gr.ptr, tetx, tetz, tety)
end

function zoom(gr::Graph, x1::Real, y1::Real, x2::Real, y2::Real)
    cmgl.set_zoom(gr.ptr, x1, y1, x2, y2)
end

# Additional helper functions
function setsize(gr::Graph, width::Int, height::Int)
    cmgl.set_size(gr.ptr, width, height)
end

function quality(gr::Graph, val::Int)
    cmgl.set_quality(gr.ptr, val)
end

function quality(gr::Graph)
    cmgl.get_quality(gr.ptr)
end

function startgroup(gr::Graph, name::String)
    cmgl.start_group(gr.ptr, name)
end

function endgroup(gr::Graph)
    cmgl.end_group(gr.ptr)
end

# File Export
function write(gr::Graph, fname::String="", descr::String="")
    cmgl.write_frame(gr.ptr, fname, descr)
end

function writepng(gr::Graph, fname::String, descr::String="", compr::Int=""; alpha::Bool=true)
    if alpha
        cmgl.write_png(gr.ptr, fname, descr)
    else
        cmgl.write_png_solid(gr.ptr, fname, descr)
    end
end

function writejpg(gr::Graph, fname::String, descr::String="")
    cmgl.write_jpg(gr.ptr, fname, descr)
end

function writegif(gr::Graph, fname::String, descr::String="")
    cmgl.write_gif(gr.ptr, fname, descr)
end

function writebmp(gr::Graph, fname::String, descr::String="")
    cmgl.write_bmp(gr.ptr, fname, descr)
end

function writetga(gr::Graph, fname::String, descr::String="")
    cmgl.write_tga(gr.ptr, fname, descr)
end

function writeeps(gr::Graph, fname::String, descr::String="")
    cmgl.write_eps(gr.ptr, fname, descr)
end

function writebps(gr::Graph, fname::String, descr::String="")
    cmgl.write_bps(gr.ptr, fname, descr)
end

function writesvg(gr::Graph, fname::String, descr::String="")
    cmgl.write_svg(gr.ptr, fname, descr)
end

function writetex(gr::Graph, fname::String, descr::String="")
    cmgl.write_tex(gr.ptr, fname, descr)
end

function writeprc(gr::Graph, fname::String, descr::String=""; pdf=true)
    cmgl.write_prc(gr.ptr, fname, descr, pdf)
end

function writeobj(gr::Graph, fname::String, descr::String="")
    cmgl.write_obj(gr.ptr, fname, descr)
end

function writexyz(gr::Graph, fname::String, descr::String="")
    cmgl.write_xyz(gr.ptr, fname, descr)
end

function writestl(gr::Graph, fname::String, descr::String="")
    cmgl.write_stl(gr.ptr, fname, descr)
end

function writeoff(gr::Graph, fname::String, descr::String=""; colored::Bool=true)
    cmgl.write_off(gr.ptr, fname, descr, colored)
end

function writejson(gr::Graph, fname::String, descr::String="")
    cmgl.write_json(gr.ptr, fname, descr)
end

function exportmgld(gr::Graph, fname::String, descr::String="")
    cmgl.export_mgld(gr.ptr, fname, descr)
end

function importmgld(gr::Graph, fname::String; add::Bool=true)
    cmgl.import_mgld(gr.ptr, fname, add)
end

#=#TODO:show image function=#
function showwith(gr::Graph, viewer::String; wait::Bool=true)
    cmgl.show_image(gr.ptr, viewer, wait)
end

function newframe(gr::Graph)
    cmgl.new_frame(gr.ptr)
end

function endframe(gr::Graph)
    cmgl.end_frame(gr.ptr)
end

function getnumframe(gr::Graph)
    cmgl.get_num_frame(gr.ptr)
end

function setnumframe(gr::Graph, n::Int)
    cmgl.set_num_frame(gr.ptr, n)
end

function getframe(gr::Graph, n::Int)
    cmgl.get_frame(gr.ptr, n)
end

function showframe(gr::Graph, n::Int)
    cmgl.show_frame(gr.ptr, n)
end

function delframe(gr::Graph, n::Int)
    cmgl.del_frame(gr.ptr, n)
end

function resetframes(gr::Graph)
    cmgl.reset_frames(gr.ptr)
end

function clearframe(gr::Graph, n::Int)
    cmgl.clear_frame(gr.ptr, n)
end

function startgif(gr::Graph, fname::String, ms::Int=100)
    cmgl.start_gif(gr.ptr, fname, ms)
end

function closegif(gr::Graph)
    cmgl.close_gif(gr.ptr)
end

function width(gr::Graph)
    cmgl.get_width(gr.ptr)
end
function height(gr::Graph)
    cmgl.get_height(gr.ptr)
end

function getrgb(gr::Graph)
    pointer_to_array(cmgl.get_rgb(gr.ptr), (3, width(gr), height(gr)), true)
end

function getrgba(gr::Graph)
    pointer_to_array(cmgl.get_rgba(gr.ptr), (4, width(gr), height(gr)), true)
end

function getwidth(gr::Graph)
    cmgl.get_width(gr.ptr)
end

function getheight(gr::Graph)
    cmgl.get_height(gr.ptr)
end

function calcxyz(gr::Graph, xs::Int, ys::Int)
    x = Int[0]
    y = Int[0]
    z = Int[0]
    cmgl.calc_xyz(gr.ptr, xs, ys, x, y, z)
    [x,y,z]
end

function calcscr(gr::Graph, x::Int, y::Int, z::Int)
    xs = Int[0]
    ys = Int[0]
    cmgl.calc_scr(gr.ptr, x, y, z, xs, ys)
    [xs, ys]
end

function objectid(gr::Graph, id::Int)
    cmgl.set_obj_id(gr.ptr, id)
end

function objectid(gr::Graph, xs::Int, ys::Int)
    cmgl.get_obj_id(gr.ptr, xs, ys)
end

function subplotid(gr::Graph, xs::Int, ys::Int)
    cmgl.get_spl_id(gr.ptr, xs, ys)
end

function highlight(gr::Graph, id::Int)
    cmgl.highlight(gr.ptr, id)
end

function isactive(gr::Graph, xs::Int, ys::Int, d::Int=1)
    cmgl.is_active(gr.ptr, xs, ys, d)
end

function drawregion(gr::Graph, nx::Int, ny::Int, m::Int)
    cmgl.set_draw_reg(gr.ptr, nx, ny, m)
end

## NOT TO BE EXPORTED ###
function numthreads(n::Int)
    cmgl.set_num_thr(n)
end

function combine(gr::Graph, gr2::Graph)
    cmgl.combine_gr(gr.ptr, gr2.ptr)
end

function mpisend(gr::Graph, id::Int)
    cmgl.mpi_send(gr.ptr, id)
end

function mpirecv(gr::Graph, id::Int)
    cmgl.mpi_recv(gr.ptr, id)
end

function clf(gr::Graph)
    cmgl.clf(gr.ptr)
end

function clf(gr::Graph, col::String)
    cmgl.clf_str(gr.ptr, col)
end

function clf(gr::Graph, col::Char)
    cmgl.clf_chr(gr.ptr, col)
end

function clf(gr::Graph, r::Real, g::Real, b::Real)
    cmgl.clf_rgb(gr.ptr, r, g, b)
end

function rasterize(gr::Graph)
    cmgl.rasterize(gr.ptr)
end

function background(gr::Graph, fname::String, alpha::Float64=1)
    cmgl.load_background(gr.ptr, fname, alpha)
end


text:
module ReactiveMPModelsTemplateTest

using Test, InteractiveUtils
using Rocket, ReactiveMP, GraphPPL, Distributions
using BenchmarkTools, Random, Plots, Dates, LinearAlgebra, StableRNGs

# Please use StableRNGs for random number generators

## Model definition
## -------------------------------------------- ##
@model function template_model(n, x0, c_, P_)
    error("define model")
end
## -------------------------------------------- ##
## Inference definition
## -------------------------------------------- ##
function template_inference(data, x0, c, P)
    error("define inference")
end

@testset "Model template" begin

    @testset "Use case template" begin
        ## -------------------------------------------- ##
        ## Data creation
        ## -------------------------------------------- ##
        _
        ## -------------------------------------------- ##
        ## Inference execution
        -
        ## -------------------------------------------- ##
        ## Test inference results
        _
        ## -------------------------------------------- ##
        ## Form debug output
        base_output = joinpath(pwd(), "_output", "models")
        mkpath(base_output)
        timestamp        = Dates.format(now(), "dd-mm-yyyy-HH-MM")
        plot_output      = joinpath(base_output, "template_model_plot_$(timestamp)_v$(VERSION).png")
        benchmark_output = joinpath(base_output, "template_model_benchmark_$(timestamp)_v$(VERSION).txt")
        ## -------------------------------------------- ##
        ## Create output plots

        savefig(p, plot_output)
        ## -------------------------------------------- ##
        ## Create output benchmarks
        benchmark = @benchmark 1 + 1#
        open(benchmark_output, "w") do io
            show(io, MIME("text/plain"), benchmark)
            versioninfo(io)
        end
        ## -------------------------------------------- ##
    end

end

end


text:
# Objective functions
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export log_likelihood, ∇log_likelihood, Hlog_likelihood, mse, ∇mse, Hmse

###################################################################################################
# Mean squared error
"""
    f = mse(X, Y)

Mean squared error between arrays/tensor X and Y

See also: [`∇mse`](@ref)
"""
mse(X::AbstractArray{Float32, N}, Y::AbstractArray{Float32, N}) where N = .5f0/size(X, N)*norm(X - Y, 2)^2

"""
    ∇f = ∇mse(X, Y)

Gradient of the MSE loss with respect to input tensors X and Y.

See also: [`mse`](@ref)
"""
∇mse(X::Array{Float32, N}, Y::Array{Float32, N}) where N = 1f0/size(X, N)*(X - Y)

"""
    Hf = Hmse(X, Y)

Hessian of the MSE loss with respect to input tensors X.

See also: [`mse`](@ref)
"""
function Hmse(X::Array{Float32, N}, Y::Array{Float32, N}) where N
    return InvertibleNetworkLinearOperator{Array{Float32, N},Array{Float32, N}}(
        ΔX -> 1f0/size(X, N)*ΔX,
        ΔX -> 1f0/size(X, N)*ΔX)
end

###################################################################################################
# Log-likelihood

"""
    f = log_likelihood(X; μ=0f0, σ=1f0)

Log-likelihood of X for a Gaussian distribution with given mean μ and variance σ. All elements of X are assumed to be iid.

See also: [`∇log_likelihood`](@ref)
"""
log_likelihood(X::AbstractArray{Float32, N}; μ=0f0, σ=1f0) where N = 1f0/size(X, N)*sum(-.5f0*((X .- μ)/σ).^2)

"""
    Hf = Hlog_likelihood(X; μ=0f0, σ=1f0)

Hessian of the log-likelihood function with respect to the input tensor X.

See also: [`log_likelihood`](@ref)
"""
∇log_likelihood(X::AbstractArray{Float32, N}; μ=0f0, σ=1f0) where N = -1f0/size(X, N)*(X .- μ)/σ^2

"""
    Hf = Hlog_likelihood(X; μ=0f0, σ=1f0)

Hessian of the log-likelihood function with respect to the input tensor X.

See also: [`log_likelihood`](@ref)
"""
function Hlog_likelihood(X::AbstractArray{Float32, N}; μ=0f0, σ=1f0) where N
    return InvertibleNetworkLinearOperator{AbstractArray{Float32, N},AbstractArray{Float32, N}}(
        ΔX -> -1f0/size(X, N)*ΔX/σ^2,
        ΔX -> -1f0/size(X, N)*ΔX/σ^2)
end


text:
using Parser
for n in names(CSTParser, true, true)
    eval(:(import CSTParser.$n))
end

applychange(str, s1, s2, t) = string(str[1:s1], t, str[s2 + 1:end])

str, s1, s2, txt = "a + b", 0, 0, "x"
str, s1, s2, txt = "a + b", 1, 1, "x"
str, s1, s2, txt = "a + b", 2, 2, "x"
str, s1, s2, txt = "a + b", 4, 4, "b"
str, s1, s2, txt = "a + b", 5, 5, "b"
str, s1, s2, txt = "a + b", 4, 5, "c"
str, s1, s2, txt = "function fname(arg) end", 9, 10, "g"
str, s1, s2, txt = "function fname(arg) end", 15, 18, "g"
str, s1, s2, txt = """
function fname(arg)
    if cond == true
        for i = 1:20
            out += arg
        end
    end
end""", 73, 76, "d"
str, s1, s2, txt = "name", 2, 2, "+"

str1 = applychange(str, s1, s2, txt);
f = File([], [], "", Parser.parse(str, true));
reparse(f, str1, s1, s2, txt);

function reparse(f::File, newstr, startbyte, endbyte, text)
    size_new = sizeof(newstr)
    size_old = size_new - sizeof(text) - startbyte + endbyte

    if isempty(f.ast.args)
        f.ast = parse(newstr, true)
        return
    # elseif startbyte == size_old
    #     startbyte1 = 0
    #     for i in 1:length(f.ast)-1
    #         startbyte1 += f.ast[i].span
    #     end
    #     ps = ParseState(newstr)
    #     while ps.nt.startbyte < startbyte1 && !ps.done
    #         next(ps)
    #     end
    #     newblocks, ps = parse(ps, true)
    #     pop!(f.ast.args)
    #     append!(f.ast.args, newblocks.args)
    #     f.ast.span = sum(x.span for x in f.ast)
    #     return
    end

    # find old token
    # if on the last byte of a node we need to check the boundary
    first_tok, first_tree, first_ind = find(f.ast, startbyte)
    if !(first_tok isa IDENTIFIER)
        first_tok, first_tree, first_ind = find(f.ast, startbyte + 1)
    end
    last_tok, last_tree, last_ind = find(f.ast, endbyte)
    if !(last_tok isa IDENTIFIER)
        last_tok, last_tree, last_ind = find(f.ast, endbyte + 1)
    end

    if first_tok == last_tok && first_tok isa LITERAL || first_tok isa IDENTIFIER
        ps = ParseState(newstr)
        while ps.nnt.startbyte < startbyte && !ps.done
            next(ps)
        end
        if ps.nt.kind == Tokens.COMMENT
            new_tok, ps = parse(ps)
        elseif ps.nt.kind == Tokens.IDENTIFIER
            new_tok = INSTANCE(next(ps))
        else
            return
        end

        insert_range = endbyte - startbyte
        if new_tok.span == first_tok.span + sizeof(text) - insert_range
            first_tree[end][first_ind[end]] = new_tok
            Dspan = new_tok.span - first_tok.span
            for x in first_tree
                x.span += 1
            end
        end
    end

    if Expr(parse(newstr, true)) != Expr(f.ast)
        info("reparse failed")
        info(Expr(parse(newstr, true)))
        info(Expr(f.ast))
        # f.ast = parse(newstr, true)
    end
end

function get_byteposition(doc, line, character)
    line_offsets = get_line_offsets(doc)

    current_offset = line == 0 ? 0 : line_offsets[line]
    for i = 1:character
        current_offset = nextind(doc._content, current_offset)
    end
    return current_offset
end


text:
@testset "Issue 156" begin
    @test length(children(node(:div, rand(3, 3)))) == 1

    n = rand(10, 10) |> dom"div"
    @test length(children(n)) == 1  # should have one child, and that child should be the given array
    @test size(first(children(n))) == (10, 10)
end

# Nodes don't support equality tests currently, so
# we have to do our own comparison of their contents
function compare_nodes(n1::Node, n2::Node)
    (instanceof(n1) == instanceof(n2) &&
     all(compare_nodes(c1, c2) for (c1, c2) in zip(children(n1), children(n2))) &&
     props(n1) == props(n2))
end
compare_nodes(x1, x2) = x1 == x2

@testset "Constructing Nodes from Strings" begin
    @test compare_nodes(node("div", "hello", "world"),
                        node(:div, "hello", "world"))
end

@testset "Setting node fields" begin
    n1 = node(:div, node(:div, "hello"), node(:div, "world"), foo="bar")
    n2 = node(:div, node(:div, "hello"), node(:div, "world"), foo="bar")

    n3 = setinstanceof(n1, :ul)
    # set... methods shouldn't mutate the original object, so
    # n1 and n2 are still the same
    @test compare_nodes(n1, n2)
    # but n1 and n3 are now different
    @test !compare_nodes(n1, n3)
    @test instanceof(n3) == WebIO.promote_instanceof(:ul)

    n3 = setchildren(n1, [node(:ul, "list"), node(:li, "item")])
    # set... methods shouldn't mutate the original object, so
    # n1 and n2 are still the same
    @test compare_nodes(n1, n2)
    # but n1 and n3 are now different
    @test !compare_nodes(n1, n3)
    @test length(children(n3)) == 2
    @test compare_nodes(children(n3)[1], node(:ul, "list"))
    @test compare_nodes(children(n3)[2], node(:li, "item"))

    n3 = setprops(n1, Dict(:one => 1, :two => 2))
    # set... methods shouldn't mutate the original object, so
    # n1 and n2 are still the same
    @test compare_nodes(n1, n2)
    # but n1 and n3 are now different
    @test !compare_nodes(n1, n3)
    @test props(n3) == Dict(:one => 1, :two => 2)
end


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

"""
Convenience functions for metaprogramming.
"""
module Meta

using ..CoreLogging

export quot,
       isexpr,
       show_sexpr,
       @dump

quot(ex) = Expr(:quote, ex)

isexpr(ex::Expr, head)          = ex.head === head
isexpr(ex::Expr, heads::Union{Set,Vector,Tuple}) = in(ex.head, heads)
isexpr(ex,       head)          = false

isexpr(ex,       head, n::Int)  = isexpr(ex, head) && length(ex.args) == n

# ---- show_sexpr: print an AST as an S-expression ----

show_sexpr(ex) = show_sexpr(stdout, ex)
show_sexpr(io::IO, ex) = show_sexpr(io, ex, 0)
show_sexpr(io::IO, ex, indent::Int) = show(io, ex)

const sexpr_indent_width = 2

function show_sexpr(io::IO, ex::QuoteNode, indent::Int)
    inner = indent + sexpr_indent_width
    print(io, "(:quote, #QuoteNode\n", " "^inner)
    show_sexpr(io, ex.value, inner)
    print(io, '\n', " "^indent, ')')
end
function show_sexpr(io::IO, ex::Expr, indent::Int)
    inner = indent + sexpr_indent_width
    print(io, '(')
    show_sexpr(io, ex.head, inner)
    for arg in ex.args
        print(io, ex.head === :block ? ",\n"*" "^inner : ", ")
        show_sexpr(io, arg, inner)
    end
    if isempty(ex.args); print(io, ",)")
    else print(io, (ex.head === :block ? "\n"*" "^indent : ""), ')')
    end
end

"""
    @dump expr

Show every part of the representation of the given expression. Equivalent to
[`dump(:(expr))`](@ref dump).
"""
macro dump(expr)
    dump(expr)
end

"""
    lower(m, x)

Takes the expression `x` and returns an equivalent expression in lowered form
for executing in module `m`.
See also [`code_lowered`](@ref).
"""
lower(m::Module, @nospecialize(x)) = ccall(:jl_expand, Any, (Any, Any), x, m)

"""
    @lower [m] x

Return lowered form of the expression `x` in module `m`.
By default `m` is the module in which the macro is called.
See also [`lower`](@ref).
"""
macro lower(code)
    return :(lower($__module__, $(QuoteNode(code))))
end
macro lower(mod, code)
    return :(lower($(esc(mod)), $(QuoteNode(code))))
end

## interface to parser ##

"""
    ParseError(msg)

The expression passed to the [`parse`](@ref) function could not be interpreted as a valid Julia
expression.
"""
struct ParseError <: Exception
    msg::AbstractString
end

"""
    parse(str, start; greedy=true, raise=true, depwarn=true)

Parse the expression string and return an expression (which could later be passed to eval
for execution). `start` is the index of the first character to start parsing. If `greedy` is
`true` (default), `parse` will try to consume as much input as it can; otherwise, it will
stop as soon as it has parsed a valid expression. Incomplete but otherwise syntactically
valid expressions will return `Expr(:incomplete, "(error message)")`. If `raise` is `true`
(default), syntax errors other than incomplete expressions will raise an error. If `raise`
is `false`, `parse` will return an expression that will raise an error upon evaluation. If
`depwarn` is `false`, deprecation warnings will be suppressed.

```jldoctest
julia> Meta.parse("x = 3, y = 5", 7)
(:(y = 5), 13)

julia> Meta.parse("x = 3, y = 5", 5)
(:((3, y) = 5), 13)
```
"""
function parse(str::AbstractString, pos::Int; greedy::Bool=true, raise::Bool=true,
               depwarn::Bool=true)
    # pos is one based byte offset.
    # returns (expr, end_pos). expr is () in case of parse error.
    bstr = String(str)
    # For now, assume all parser warnings are depwarns
    ex, pos = with_logger(depwarn ? current_logger() : NullLogger()) do
        ccall(:jl_parse_string, Any,
              (Ptr{UInt8}, Csize_t, Int32, Int32),
              bstr, sizeof(bstr), pos-1, greedy ? 1 : 0)
    end
    if raise && isa(ex,Expr) && ex.head === :error
        throw(ParseError(ex.args[1]))
    end
    if ex === ()
        raise && throw(ParseError("end of input"))
        ex = Expr(:error, "end of input")
    end
    return ex, pos+1 # C is zero-based, Julia is 1-based
end

"""
    parse(str; raise=true, depwarn=true)

Parse the expression string greedily, returning a single expression. An error is thrown if
there are additional characters after the first expression. If `raise` is `true` (default),
syntax errors will raise an error; otherwise, `parse` will return an expression that will
raise an error upon evaluation.  If `depwarn` is `false`, deprecation warnings will be
suppressed.

```jldoctest
julia> Meta.parse("x = 3")
:(x = 3)

julia> Meta.parse("x = ")
:($(Expr(:incomplete, "incomplete: premature end of input")))

julia> Meta.parse("1.0.2")
ERROR: Base.Meta.ParseError("invalid numeric constant \\\"1.0.\\\"")
Stacktrace:
[...]

julia> Meta.parse("1.0.2"; raise = false)
:($(Expr(:error, "invalid numeric constant \"1.0.\"")))
```
"""
function parse(str::AbstractString; raise::Bool=true, depwarn::Bool=true)
    ex, pos = parse(str, 1, greedy=true, raise=raise, depwarn=depwarn)
    if isa(ex,Expr) && ex.head === :error
        return ex
    end
    if !done(str, pos)
        raise && throw(ParseError("extra token after end of expression"))
        return Expr(:error, "extra token after end of expression")
    end
    return ex
end

end # module

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
push!(LOAD_PATH, "../src")
using Base.Test
using PumiInterface
using ODLCommonTools
using SummationByParts
using PdePumiInterface

function not_isapprox(args...; kwargs...)
  return !isapprox(args...; kwargs...)
end

include("defs.jl")
include("common_functions.jl")
include("test_funcs.jl")
include("test_adapt.jl")

  dmg_name = "parallel2.dmg"
  smb_name = "./parallel2.smb"
  order = 1
@testset "Testing PUMIInterface.jl" begin

  m_ptr, dim = apf.loadMesh(dmg_name, smb_name, order)
  mshape_ptr, num_Entities, n_arr = apf.initMesh(m_ptr)
#  num_Entities, m_ptr, mshape_ptr = apf.init2(dmg_name, smb_name, order)

  myrank = MPI.Comm_rank(MPI.COMM_WORLD)
  @test ( num_Entities[1] )== 6
  @test ( num_Entities[2] )== 9
  @test ( num_Entities[3] )== 4

  @test ( apf.countPeers(m_ptr, apf.VERTEX) )== 1
  peers = Array{Int32}(1)
  apf.getPeers(m_ptr, peers)

  @test ( peers[1] )== 1-myrank

  verts = Array{Ptr{Void}}(num_Entities[1])
  edges = Array{Ptr{Void}}(num_Entities[2])
  faces = Array{Ptr{Void}}(num_Entities[3])

  for (i, e) in enumerate(apf.MeshIterator(m_ptr, 0))
    verts[i] = e
  end

  for (i, e) in enumerate(apf.MeshIterator(m_ptr, 1))
    edges[i] = e
  end

  for (i, e) in enumerate(apf.MeshIterator(m_ptr, 2))
    faces[i] = e
  end

  # check the apf.countRemotes function
  # the only thing that can really be tested is the number of total remotes
  remote_cnt = 0 # count the total number of remotes
  remote_sum = 0 # sum of the remote values
  part_nums = Array{Cint}(1)
  entities = Array{Ptr{Void}}(1)

  sleep(5*myrank)
  for i=1:length(verts)
    nremotes = apf.countRemotes(m_ptr, verts[i])
    remote_cnt += Int(nremotes)
    if nremotes != 0
      apf.getRemotes(part_nums, entities)
      remote_sum += Int(part_nums[1])
    end
  end

  @test ( remote_cnt )== 3
  @test ( remote_sum )== 3*(1 - myrank)

  remote_cnt = 0
  remote_sum = 0
  for i=1:length(edges)
    nremotes =  apf.countRemotes(m_ptr, edges[i])
    remote_cnt += Int(nremotes)
    if nremotes != 0
      apf.getRemotes(part_nums, entities)
      remote_sum += Int(part_nums[1])
    end
  end

  @test ( remote_cnt )== 2
  @test ( remote_sum )== 2*(1-myrank)
end

@testset "----- Testing PdePumiInterfaceDG -----" begin
  opts = Dict{Any, Any}(
    "dmg_name" => dmg_name,
    "smb_name" => smb_name,
    "coloring_distance" => 2,
    "order" => order,
    "numBC" => 1,
    "BC1" =>  [0,1,2,3],
    "run_type" => 1,
    "verify_coloring" => true,
    "use_edge_res" => false,
    "write_edge_vertnums" => true,
    "write_face_vertnums" => true,
    "write_boundarynums" => true,
    "write_dxidx" => false,
    "write_coords" => false,
    "write_sparsity" => false,
    "write_offsets" => false,
    "write_counts" => false,
    "write_sparsity_nodebnds" => false,
    "write_dofs" => false,
    "write_interfaces" => false,
    "write_boundaries" => false,
    "write_sharedboundaries" => false,
    "use_linear_metrics" => true,
    )

  sbp = getTriSBPOmega(degree=order)
#  sbp = TriSBP{Float64}(degree=order, internal=true)

  vtx = [-1. 1 -1; -1 -1 1]
#  interp_op = SummationByParts.buildinterpolation(sbp, vtx)
  sbpface = TriFace{Float64}(order, sbp.cub, vtx.')
  mesh_c = PumiMeshDG2(Complex128, sbp, opts, sbpface)
#  mesh_c = PumiMeshDG2{Complex128, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface)
  mesh = PumiMeshDG2(Float64, sbp, opts, sbpface)
#  mesh = PumiMeshDG2{Float64, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface)

  @test ( mesh.numVert )== 6
  @test ( mesh.numEdge )== 9
  @test ( mesh.numEl )== 4
  @test ( mesh.numInterfaces )== 3
  @test ( mesh.myrank )== MPI.Comm_rank(mesh.comm)
  @test ( mesh.commsize )== 2  # these tests only work for 2 ranks
  @test ( mesh.peer_parts[1] )== 1-(mesh.myrank)
  @test ( mesh.peer_face_counts[1] )== 2
  @test ( length(mesh.bndries_local[1]) )== 2
  @test ( length(mesh.bndries_remote[1]) )== 2
  @test ( length(mesh.shared_interfaces[1]) )== 2

  # check the interfaces
  for i=1:length(mesh.shared_interfaces[1])
    interface_i = mesh.shared_interfaces[1][i]
    @test  interface_i.elementL  > 0
    @test  interface_i.elementL  < mesh.numEl + 1
    @test  interface_i.elementR  > mesh.numEl
    @test  interface_i.elementR  < mesh.numEl + 3
    @test ( interface_i.orient )== 1
  end

  @test ( mesh.shared_element_offsets[1] )== (mesh.numEl + 1)
  @test ( mesh.shared_element_offsets[2] )== (mesh.numEl + 3)

  # check the gathering of the elements on the boundaries
  for i=1:mesh.npeers
    peernum = mesh.peer_parts[i]
    local_els = mesh.local_element_lists[i]
    for j = 1:length(local_els)
      el = local_els[j]
      el_ptr = mesh.elements[el]
      down, tmp = PumiInterface.apf.getDownward(mesh.m_ptr, el_ptr, apf.EDGE)
      cnt = 0 # count number of shared edges
      for k=1:3
        if apf.isShared(mesh.m_ptr, down[k])
          nremotes = PumiInterface.apf.countRemotes(mesh.m_ptr, down[k])
          partnums = Array{Cint}(nremotes)
          entities = Array{Ptr{Void}}(nremotes)
          PumiInterface.apf.getRemotes(partnums, entities)
          for p=1:nremotes
            if partnums[p] == peernum
              cnt += 1
            end  # end if
          end  # end loop p
        end  # end if apf.isShared
      end  # end loop k

      @test  cnt  > 0
    end  # end loop j

    # check element numbers are unique
    local_elnums = copy(local_els)
    sort!(local_elnums)
    @test ( local_elnums )== unique(local_elnums)
  end  # end loop i

  # check the interpolation
  # use the fact that dxidx, jac are constant across the element
  for p=1:mesh.npeers
    interfaces_p = mesh.shared_interfaces[p]
    dxidx_p = mesh.dxidx_sharedface[p]
    jac_p = mesh.jac_sharedface[p]
    for i=1:mesh.peer_face_counts[p]
      interface_i = interfaces_p[i]
      dxidx_ref = mesh.dxidx[:, :, 1, interface_i.elementL]
      jac_ref= mesh.jac[1, interface_i.elementL]

      for j=1:mesh.sbpface.numnodes
        for k=1:2
          for m=1:2
            @test isapprox( dxidx_p[k, m, j, i], dxidx_p[k,m]) atol=1e-12
          end
        end
        @test isapprox( jac_p[j, i], jac_ref) atol=1e-12
      end
    end
  end

  # check the adjacency dictonary
  adj_dict, revadj = PdePumiInterface.getLocalAdjacency(mesh)
  @test ( length(adj_dict) )== mesh.numSharedEl
  # now check that no element appears more than twice
  count_dict = Dict{Int, Int}()
  for (key, val) in adj_dict
    el1 = val[1]
    el2 = val[2]
    for elnums = 1:2
      val_i = val[elnums]
      if !haskey(count_dict, val_i)
        count_dict[val_i] = 1
      else
        count_dict[val_i] = count_dict[val_i] + 1
      end
    end
  end

  for (key, val) in count_dict
    @test  val   < 3
    @test ( val[1] )!=0
  end

  # check the revadj
  count_dict2 = Dict{Int, Int}()
  for i=1:mesh.numSharedEl
    val1 = revadj[i, 1]
    val2 = revadj[i, 2]
    @test ( val1 )!=0
    if !haskey(count_dict2, val1)
      count_dict2[val1] = 1
    else
      count_dict2[val1] += 1
    end

    if val2 != 0
      if !haskey(count_dict2, val2)
        count_dict2[val1] = 1
      else
        count_dict2[val1] += 1
      end
    end

  end

  # check no value appears more than twice
  for (key, val) in count_dict2
    @test  val  < 3
    @test ( val )!=0
  end

  # check neighbor_colors for uniqueness
  myrank = mesh.myrank
  for i=1:mesh.numEl
    colors_i = mesh.neighbor_colors[:, i]
    sort!(colors_i)
    nzs = countnz(colors_i)
    nz_start = length(colors_i) - nzs + 1
    colors_i2 = colors_i[nz_start:end]
    @test ( colors_i2 )== unique(colors_i2)
  end

  for i=1:mesh.npeers
    numel_i = mesh.shared_element_offsets[i+1] - mesh.shared_element_offsets[i]
    masks = mesh.shared_element_colormasks[i]
    for j = 1:numel_i
      val = 0
      for k = 1:mesh.numColors
        val += masks[k][j]
      end
      @test ( val )== 1  # each element should have exactly 1 color
    end
  end

  # check dof numbers for non-local elements
  npeer_els = mesh.shared_element_offsets[end] = mesh.shared_element_offsets[1]
  dof_min = mesh.numDofPerNode*mesh.numNodesPerElement*mesh.numEl
  dof_max = mesh.numDofPerNode*mesh.numNodesPerElement*(npeer_els + mesh.numEl)  + 1
  for i=(mesh.numEl + 1):size(mesh.dofs, 3)
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.numDofPerNode
        if mesh.myrank == 0
          @test  (mesh.dofs[k, j, i] + mesh.dof_offset)  > dof_min
          @test  (mesh.dofs[k, j, i] + mesh.dof_offset)  < dof_max
        else
          @test  mesh.dofs[k, j, i]  < 1
        end
      end
    end
  end

  vshare = mesh.vert_sharing
  @test ( vshare.npeers )== 1
  @test ( vshare.peer_nums )== [1 - mesh.myrank]
  @test ( vshare.counts[1] )== 3
  @test ( length(vshare.vert_nums[1]) )== 3
  @test ( length(keys(vshare.rev_mapping)) )== 3

  # check vert_nums and rev_mapping in more detail
  for i=1:3
    @test  vshare.vert_nums[1][i]  > 0
    @test  vshare.vert_nums[1][i]  < mesh.numVert + 1
  end

  for (key, val) in vshare.rev_mapping
    @test ( length(val.first) )== 1
    @test ( length(val.second) )== 1
    @test ( val.first[1] )== 1-mesh.myrank  # part number
    @test  val.second[1]  > 0
    @test  val.second[1]  < 4# number of shared vertices + 1
  end

  compare_meshes(mesh, mesh_c)
  test_metrics_rev(mesh, mesh_c, sbp, opts)

  MPI.Barrier(mesh.comm)

end

@testset "----- Testing PdePumiInterface3DG -----" begin

  degree = 1
  Tsbp = Float64
  sbp = getTetSBPOmega(degree=degree)
#  sbp = TetSBP{Tsbp}(degree=degree, internal=true)
  ref_verts = sbp.vtx
  interp_op = SummationByParts.buildinterpolation(sbp, ref_verts.')
  face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
  topo = ElementTopology{3}(face_verts)
  sbpface = TetFace{Tsbp}(degree, sbp.cub, ref_verts)

  dmg_name = ".null"
  smb_name = "pcube2.smb"

  opts = PdePumiInterface.get_defaults()
  opts["dmg_name"] = dmg_name
  opts["smb_name"] = smb_name
  opts["use_linear_metrics"] = true
  opts["coloring_distance"] = 2
  opts["order"] = degree
  opts["numBC"] = 1
  opts["BC1"] = [0,1,2,3,4,5]

#  interp_op = eye(4)
  mesh_c = PumiMeshDG3(Complex128, sbp, opts, sbpface, topo)
#  mesh_c = PumiMeshDG3{Complex128, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)
  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  mesh = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  @test ( mesh.numEl )== 12
  @test ( mesh.numGlobalEl )== 16
#  @test ( mesh.numBoundaryFaces )==#  @fact mesh.numBoundaryFaces --> 4
  @test ( mesh.numInterfaces )== 14
  @test ( mesh.numBoundaryFaces )== 16
  @test ( mesh.myrank )== MPI.Comm_rank(mesh.comm)
  @test ( mesh.commsize )== MPI.Comm_size(mesh.comm)

  if mesh.myrank == 0
    @test ( mesh.peer_parts )== [1]
  else
    @test ( mesh.peer_parts )== [0]
  end

  @test ( mesh.peer_face_counts[1] )== 4

  # TODO: test dxidx, jac sharedface
  # check interface array
  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    @test  iface_i.elementL  > 0
    @test  iface_i.elementL  < mesh.numEl+1
    @test  iface_i.elementR  > 0
    @test  iface_i.elementR  < mesh.numEl+1
    @test  iface_i.faceL  > 0
    @test  iface_i.faceL  < 5
    @test  iface_i.faceR  > 0
    @test  iface_i.faceR  < 5
    @test  iface_i.orient  > 0
    @test  iface_i.orient  < 4
  end

  # check boundary array
  for i=1:mesh.numBoundaryFaces
    bndry_i = mesh.bndryfaces[i]
    @test  bndry_i.element  > 0
    @test  bndry_i.element  < mesh.numEl + 1
    @test  bndry_i.face  > 0
    @test  bndry_i.face  < 5
  end

  # check mapping interpolation
  # should be constant within an element for straight-sided elements
  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    el_i = iface_i.elementL
    dxidx_el = mesh.dxidx[:, :, 1, el_i]
    jac_el = mesh.jac[:, el_i]
    jac_face = mesh.jac_face[:, i]

    for j=1:mesh.numNodesPerFace
      dxidx_face = mesh.dxidx_face[:, :, j, i]
      for k=1:3
        for p=1:3
          @test isapprox( dxidx_face[p, k], dxidx_el[p, k]) atol=1e-13
        end
      end
      @test isapprox( jac_face[j], jac_el[j]) atol=1e-13
    end
  end  # end loop over interfaces

  # check mapping interpolation
  for i=1:mesh.npeers

    dxidx_peer = mesh.dxidx_sharedface[i]
    jac_peer = mesh.jac_sharedface[i]
    for j=1:mesh.peer_face_counts[i]
      iface_i = mesh.shared_interfaces[i][j]
      el_j = iface_i.elementL
      dxidx_el = mesh.dxidx[:, :, 1, el_j]
      jac_el = mesh.jac[:, el_j]
      jac_face = jac_peer[:, j]

      for k=1:mesh.numNodesPerFace
        dxidx_face = dxidx_peer[:, :, k, j]
        for p=1:3
          for n=1:3
            @test isapprox( dxidx_face[n, p], dxidx_el[n, p]) atol=1e-13
          end
        end
        @test isapprox( jac_face[k], jac_el[k]) atol=1e-13
      end
    end   # end loop over peer faces
  end  # end loop over peers

  # test curvilinear metrics reverse mode
  compare_meshes(mesh, mesh_c)
  test_metrics_rev(mesh, mesh_c, sbp, opts)

  # test periodic things
  smb_name = "tet2_pxz_p2_.smb"

  opts = PdePumiInterface.get_defaults()
  opts["smb_name"] = smb_name
  opts["dmg_name"] = dmg_name
  opts["use_linear_metrics"] = true
  opts["coloring_distance"] = 2
  opts["order"] = degree
  opts["numBC"] = 1
  opts["BC1"] = [0,2,4,5]

  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  mesh = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  @test ( mesh.numPeriodicInterfaces )== 4
  @test ( mesh.numInterfaces )== (mesh.numFace - 8 - 3*4 - 8)
  @test ( mesh.peer_face_counts[1] )== 8

  smb_name = "tet2_pxy_p2_.smb"
  opts["smb_name"] = smb_name
  opts["BC1"] = [1,2,3,4]
  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  mesh = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  @test ( mesh.numPeriodicInterfaces )== 0
  @test ( mesh.numInterfaces )== (mesh.numFace - 2*8 - 4*4)
  @test ( mesh.peer_face_counts[1] )== 16

  vshare = mesh.vert_sharing
  @test ( vshare.npeers )== 1
  @test ( vshare.peer_nums )== [1 - mesh.myrank]
  @test ( vshare.counts[1] )== 18
  @test ( length(vshare.vert_nums[1]) )== 18
  @test ( length(keys(vshare.rev_mapping)) )== 18

  # check vert_nums and rev_mapping in more detail
  for i=1:18
    @test  vshare.vert_nums[1][i]  > 0
    @test  vshare.vert_nums[1][i]  < mesh.numVert + 1
  end

  for (key, val) in vshare.rev_mapping
    @test ( length(val.first) )== 1
    @test ( length(val.second) )== 1
    @test ( val.first[1] )== 1-mesh.myrank  # part number
    @test  val.second[1]  > 0
    @test  val.second[1]  < 19# number of shared vertices + 1
  end

  MPI.Barrier(mesh.comm)

  # check assertion for ambiguous parallel things
  dmg_name = ".null"
  smb_name = "tet111_.smb"

  opts = PdePumiInterface.get_defaults()
  opts["dmg_name"] = dmg_name
  opts["smb_name"] = smb_name
  opts["order"] = degree
  opts["numBC"] = 0
  opts["coloring_distance"] = 2

  @test_throws Exception  PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  @test_throws Exception #  @fact_throws PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  # check 2 x 2 x 2 mesh works
  dmg_name = ".null"
  smb_name = "tet222_.smb"

  opts["dmg_name"] = dmg_name
  opts["smb_name"] = smb_name

  PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  # this doesn't work because each partition is a 1 x 1 x 1 cube with periodic faces
  dmg_name = ".null"
  smb_name = "tet211_.smb"
  opts["dmg_name"] = dmg_name
  opts["smb_name"] = smb_name

  @test_throws Exception  PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  @test_throws Exception #  @fact_throws PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

end

MPI.Barrier(MPI.COMM_WORLD)
if MPI.Initialized()
  MPI.Finalize()
end