text:
using JSON
using Dates

include("struct.jl")
include("generateLatexTables.jl")

function autoExpe(expeJsonPath::String)

    parameters = readExpeFormat(expeJsonPath)

    # Number of resolution performed in this experiment
    resolutionCount = length(parameters.resolutionMethods) * length(parameters.instancesPath)
    for (key, value) in parameters.parametersToCombine
        resolutionCount *= length(value)
    end

    resolutionId = 0
    isFirstResolution = false
    savedResultsFound = false

    if !isdir(parameters.outputPath)
        mkpath(parameters.outputPath)
    end

    for instancePath in parameters.instancesPath

        while occursin("//", instancePath)
            instancePath = replace(instancePath, "//" => "/")
        end

        instanceName = splitext(basename(instancePath))[1]

        outputFile = parameters.outputPath * "/" * instanceName * ".json"
        while occursin("//", outputFile)
            outputFile = replace(outputFile, "//" => "/")
        end

        # Variable which contains all the results of the instance
        # Each result is a Dict which contains:
        # - one entry with the key "resolutionMethodName"
        # - one entry for each parameter
        # - one entry for each result
        instanceResults = Vector{Dict{String, Any}}()

        # Get the results previously computed for this instance if any
        if isfile(outputFile)
            println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS"), " Reading previous results from file \"", outputFile, "\"")
            stringdata=join(readlines(outputFile))

            try
                instanceResults = JSON.parse(stringdata)
            catch e
                id = 1
                newName = outputFile * "_" * string(id)

                while isfile(newName)
                    id += 1
                    newName = outputFile * "_" * string(id)
                end
                println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS") ,
                    " Warning: Unable to read the JSON result file. A new one will be created and the older one is moved to \"", newName, "\"")
            end

        end

        if parameters.verbosity != :None
            println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS"), " Solving instance ", instancePath, " (", floor(Int, 100 * resolutionId / resolutionCount),"%)")
        end

        # Associate one index to each parameter in the dictionary of parameters to combine
        # combinationId[i] is the index of the value of parameter i that will be considered in the next combination
        combinationId = Vector{Int}()

        # Variable used for experiments which do not have combination of parameters
        isFirstCombination = true

        # Set the time of the next compilation (another may happen before if saved results are found)
        nextCompilationTime = Dates.now() + Dates.Minute(parameters.minLatexCompilationInterval)
        isFirstLatexCompiled = false

        # While there are combinations of parameters to consider for this instance
        while getNextCombinationId!(combinationId, parameters.parametersToCombine) || isFirstCombination

            # Get the value of the parameters in the combination from their ids
            dictCombination = getCombination(combinationId, parameters.parametersToCombine)
            dictCombination["instancePath"] = instancePath

            # Get these parameters in an array
            combinationKeys = collect(keys(dictCombination))

            if parameters.verbosity != :None
                println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS"),"\t Parameters: ", getString(dictCombination), " (", floor(Int, 100 * resolutionId / resolutionCount),"%)")
            end

            # For each method
            for methodId in 1:length(parameters.resolutionMethods)

                # Reinitialise the combination
                # (otherwise if a parameter p is specified for a method A and not for a method B, then method B would still have this parameter specified)
                dictCombination = getCombination(combinationId, parameters.parametersToCombine)
                dictCombination["instancePath"] = instancePath

                currentMethodName = string(nameof(parameters.resolutionMethods[methodId]))

                methodLog = string(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS")) * "\t\t Resolution method \"" * currentMethodName * "\" (" * string(floor(Int, 100 * resolutionId / resolutionCount)) * "%)... "
                resolutionId += 1

                ## Test if the resolution must be done or not
                mustSolve = parameters.computeUnsolvedResults

                # If unsolved results must be solved and if results already saved must not be recomputed, then test if this result is already solved
                if parameters.computeUnsolvedResults

                    resultId = 1
                    combinationFoundInResults = false

                    # While we have not tested that all the solved combinations are different from the one we want to compute
                    while resultId <= length(instanceResults) && !combinationFoundInResults

                        resultDict = instanceResults[resultId]

                        # If the result corresponds to the current resolution method
                        if resultDict["resolutionMethodName"] == currentMethodName

                            # Test all the parameters of the combination to see if they are equal to the parameter of the result
                            combinationFoundInResults = true
                            paramId = 1
                            while paramId <= length(combinationKeys) && combinationFoundInResults

                                # If any parameter does not correspond, the result is not a match
                                if !haskey(resultDict, combinationKeys[paramId]) || resultDict[combinationKeys[paramId]] != dictCombination[combinationKeys[paramId]]
                                    combinationFoundInResults = false
                                end
                                paramId += 1
                            end
                        end
                        if !combinationFoundInResults
                            resultId += 1
                        end
                    end

                    if combinationFoundInResults
                        savedResultsFound = true

                        # Remove the entry if the result must be recomputed
                        if parameters.recomputeSavedResults
                            deleteat!(instanceResults, resultId)
                        else
                            mustSolve = false
                            methodLog *= "Skipped (already solved)"
                        end
                    end
                else
                    methodLog *= "Skipped (not solving unsolved results)"
                end

                # If the combination must be solved for this method
                if mustSolve

                    # Create the latex result table if:
                    # 1    - a latex format is specified; and either
                    # 2.1  - the next compilation time is reached; or
                    # 2.2 (- no latex table has yet been created; and
                    #      - results have been obtained in previous run(s))
                    # Note: Condition 2.2 enables the user to see all the results of previous runs immediatly after starting the experiment
                    if parameters.latexFormatPath != [] && (Dates.now()  > nextCompilationTime  || (!isFirstLatexCompiled &&  savedResultsFound))
                        println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS") , "\t\t Creating the latex table(s) in file ", parameters.latexOutputFile)
                        createTexTables(parameters)
                        nextCompilationTime = Dates.now() + Dates.Minute(parameters.minLatexCompilationInterval)
                        isFirstLatexCompiled = true
                    end

                    println(methodLog)

                    # Solve it
                    startingTime = time()
                    dictResults = parameters.resolutionMethods[methodId](dictCombination)
                    resolutionTime = time() - startingTime

                    dictResults["auto_expe_resolution_time"] = resolutionTime

                    # The results saved both contain the experiment results and the combination of parameters considered
                    outputResult = merge(dictCombination, dictResults)
                    outputResult["resolutionMethodName"] = currentMethodName
                    push!(instanceResults, outputResult)

                    # Display the results
                    if parameters.verbosity != :None
                        println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS"), "\t\t Results: ")
                        for (key, value) in dictResults
                            print(" "^22, "\t\t\t", key, " = ")

                            try
                                # Try to round numerical values
                                roundedValue = string(round.(value, digits = 2))

                                displayedValue = value

                                # If the rounding remove decimals
                                if length(roundedValue) < length(string(value))
                                    displayedValue = roundedValue
                                end

                                if typeof(displayedValue) <: Array && length(displayedValue) > 10
                                    println("[", displayedValue[1], ", ", displayedValue[2], ", ..., ", displayedValue[end-1], ", ", displayedValue[end], "]")
                                else
                                    println(displayedValue)
                                end
                            catch e
                                # If the value cannot be rounded (i.e., if it contains non numerical values)
                                println(value)
                            end
                        end
                    end

                    # and save them
                    open(outputFile, "w") do fout
                        JSON.print(fout, instanceResults, 4)
                    end
                    println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS"), " Wrote in ", outputFile)
                else
                    if parameters.verbosity != :None
                        println(methodLog)
                    end
                end
            end

            isFirstCombination = false
        end
    end

    if length(parameters.instancesPath) == 0
        println("Warning: no instance found")
    elseif parameters.latexFormatPath != nothing
        println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS") , " Creating the latex table(s) in file ", parameters.latexOutputFile)
        createTexTables(parameters)
    end
end

function getNextCombinationId!(combinationId::Vector{Int}, parameters::Dict{String, Vector{Any}})

    isFirstCombination = isempty(combinationId)

    nextCombinationFound = false

    paramKeys = collect(keys(parameters))
    keyId = 1

    while !nextCombinationFound && keyId <= length(paramKeys)

        keyName = paramKeys[keyId]
        valueArray = parameters[keyName]

        if isFirstCombination
            if length(valueArray) > 0
                push!(combinationId, 1)
            else
                println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS"), " Parameter ", key, "has an empty size (value: ", valueArray, ").")
            end
        else
            # If the id keyId can be incremented
            if combinationId[keyId] < length(valueArray)

                # Increment it
                combinationId[keyId] += 1

                # And set the previous ids to 1
                for secondKeyId in 1:keyId - 1
                    combinationId[secondKeyId] = 1
                end

                nextCombinationFound = true
            end
        end

        keyId += 1
    end

    return nextCombinationFound
end

function getCombination(combinationId::Vector{Int}, parameters::Dict{String, Vector{Any}})

    combination = Dict{String, Any}()

    paramKeys = collect(keys(parameters))

    for keyId in 1:length(paramKeys)

        keyName = paramKeys[keyId]
        valueArray = parameters[keyName]
        value = valueArray[combinationId[keyId]]
        combination[keyName] = value
    end

    return combination

end

function getString(dict::Dict{String, Any})
    result = ""

    for (key, val) in dict
        result *= string(key) * " = \"" * string(val) * "\", "
    end

    return result
end

"""
Get the path of the instances in the json experiment

Input
- paths: each path can correspond to an instance file or to a folder. A folder is specified either by a String (its path) or a Dictionary with key "path" (and optionally key "isRecursive"). This argument either contain one folder or an array of folders;
- extensions: list of valid instance name extensions.
"""
function readInstances(paths, extensions::Vector{String})

    individualInstances = Vector{String}()

    # If there is only one instance entry
    if typeof(paths) <: Dict || typeof(paths) == String
        paths = [paths]
    end

    for pathEntry in paths

        isRecursive = false

        cPath = ""
        if typeof(pathEntry) == String
            cPath = pathEntry
        elseif typeof(pathEntry) <: Dict && haskey(pathEntry, "path")
            cPath = pathEntry["path"]
            if haskey(pathEntry, "isRecursive")
                isRecursive = pathEntry["isRecursive"]
            end
        else
            println("Warning: instance path ignored since it is not represented by a String or a Dictionary with a key \"path\":\n\tpath: ", pathEntry)
        end

        if isfile(cPath)
            if isExtensionValid(cPath, extensions)
                push!(individualInstances, cPath)
            else
                println("Warning: File ignored since its extension is not valid.\n\tpath: ", cPath, "\n\tValid extensions: ", extensions)
            end
        elseif isdir(cPath)
            append!(individualInstances, readInstanceFolder(cPath, extensions, isRecursive))
        else
            println("Warning: Instance path ignored since it does not correspond to an existing file or folder:\n\tinstance path: ", cPath)
        end
    end

    return individualInstances
end

"""
Find all the instances with a valid extension in a folder

Input:
- folderPath: path of the folder
- extensions: list of the valid extensions
- isRecursive: true if instances are sought recursively
"""
function readInstanceFolder(folderPath::String, extensions::Vector{String}, isRecursive::Bool)

    instances = Vector{String}()
    for file in readdir(folderPath)
        if isdir(file) && isRecursive
            append!(instances, readInstanceFolder(file, extensions, isRecursive))
        elseif isExtensionValid(file, extensions)
            push!(instances, folderPath * "/" * file)
        else
            println("Warning: data file ignored since it is either a folder or a file with an invalid extension.\n\tfile path: ", file)
        end
    end

    return instances
end

"""
Test if a file extension is included in a list of extensions

Input
- path: the path to test;
- validExtensions: list of the valid extensions.
"""
function isExtensionValid(path::String, validExtensions::Vector{String})

    extension = ""
    if occursin('.', path)
        extension = path[findlast(isequal('.'), path)+1:end]
    end

    return extension in validExtensions || "." * extension in validExtensions
end

"""
Get several julia functions.

Input:
- iMethods: Dictionary with an entry "name" (the name of the julia
function) and an optional entry "module" (the module in which the
function can be found). To specify several functions, iMethods can be
an array of such dictionaries. If the module is not specified, the
Main module will be considered in priority. If it is not in the Main
module, all the other loaded modules will be considered.

Output:
- vector of Function which correspond to the found functions.
"""
function getMethods(iMethods)

    # If there is only one method
    if !(typeof(iMethods) <: Vector)
        iMethods = [iMethods]
    end

    methods = Vector{Function}()

    # For each method
    for method in iMethods
        cName = nothing
        cModule = nothing

        # If it is defined by a string
        if typeof(method) == String
            cName = method

        # If it is defined by a Dict
        elseif typeof(method) <: Dict && haskey(method, "name")
            cName = method["name"]
            if haskey(method, "module")
                cModule = method["module"]
            else
                println("A resolution method must be a String or a Dictionary which contains key \"name\". Ignoring the following method:\n\tmethod: ", method, "\n\ttype: ", typeof(method))
            end
        end

        # If the function is defined properly
        if cName != nothing

            # If its module is specified
            if cModule != nothing
                try
                    push!(methods, getfield(cModule, Symbol(cName)))
                catch e
                    println("Warning: unable to find method ", cName, " in module ", cModule)
                end
            else
                println("Debug: no module specified for method ", cName)

                moduleFound = false

                # First look for the method in the Main module
                try
                    push!(methods, getfield(Main, Symbol(cName)))
                    moduleFound = true
		    println("Debug: method found in package Main.")
                catch e
                    println("Debug: no function of name ", cName, " found in module Main")
                end

                # If the method is not in the main module test all loaded modules
                if !moduleFound

                    # Get all the loaded modules
                    modules =  Base.loaded_modules_array()
                    validModules = Vector{Tuple{Module, Function}}()
                    for moduleId in 1:length(modules)
                        try
                            f = getfield(modules[moduleId], Symbol(cName))
                            push!(modules, (modules[moduleId], f))
                        catch e
                        end
                    end

                    if length(validModules) == 0
                        println("Warning: no function ", cName, " found for any loaded module. This method is ignored")
                    else
                        if length(validModules) > 1
                            println("Warning: several loaded modules contain a function ", cName, ". We consider the one from module ", validModules[1][1])
                        else
                            println("Debug: loading function ", cName, " from module ", validModules[1][1])
                        end
                        push!(methods, validModules[1][2])
                    end
                end
            end
        end
    end

    return methods

end

"""
Read the format of the experiment from a json file
"""
function readExpeFormat(jsonFilePath::String)

    stringdata=join(readlines(jsonFilePath))
    expeDict = JSON.parse(stringdata)

    if !isa(expeDict, Dict)
       println("Error: The json file describing the experiment (", jsonFilePath, ") does not contain a dictionary.")
       return nothing
    end

    # Test if the experiment format is valid
    if !haskey(expeDict, "instancesPaths") || !haskey(expeDict, "resolutionMethods")
        println("Error: The json file describing the experiment must contain a dictionary with at least an entry \"resolutionMethods\" and an entry \"instancesPaths\"")
        return nothing
    end

    # Get the valid extensions (txt and json by default)
    extensions = Vector{String}(["txt", "json"])

    if haskey(expeDict, "instanceExtensions")
        extensions = expeDict["instanceExtensions"]

        # If there is only one extension
        if !(typeof(extensions) <: Vector)
            extensions = [extensions]
        end
    end

    # Get all the instances
    instancesPath = Vector{String}()

    if haskey(expeDict, "instancesPaths")
        append!(instancesPath, readInstances(expeDict["instancesPaths"], extensions))
    end

    # Get all the methods
    methods = getMethods(expeDict["resolutionMethods"])

    parameters = ExpeParameters(instancesPath, methods)

    # Get all the optional parameters
    if haskey(expeDict, "outputPath")
        if typeof(expeDict["outputPath"]) != String
            println("Warning: the outputPath of the experiment must be a String (current value: ", expeDict["outputPath"], " current type: ", typeof(expeDict["outputPath"]), ")")
        else
            parameters.outputPath = expeDict["outputPath"]
        end
    end

    if haskey(expeDict, "recomputeSavedResults")
        if typeof(expeDict["recomputeSavedResults"]) != Bool
            println("Warning: the recomputeSavedResults of the experiment must be a boolean (current value: ", expeDict["recomputeSavedResults"], " current type: ", typeof(expeDict["recomputeSavedResults"]), ")")
        else
            parameters.recomputeSavedResults = expeDict["recomputeSavedResults"]
        end
    end

    if haskey(expeDict, "computeUnsolvedResults")
        if typeof(expeDict["computeUnsolvedResults"]) != Bool
            println("Warning: the computeUnsolvedResults of the experiment must be a boolean (current value: ", expeDict["computeUnsolvedResults"], " current type: ", typeof(expeDict["computeUnsolvedResults"]), ")")
        else
            parameters.computeUnsolvedResults = expeDict["computeUnsolvedResults"]
        end
    end

    if haskey(expeDict, "latexOutputFile")
        if typeof(expeDict["latexOutputFile"]) != String
            println("Warning: the latexOutputFile of the experiment must be a String (current value: ", expeDict["latexOutputFile"], " current type: ", typeof(expeDict["latexOutputFile"]), ")")
        else
            parameters.latexOutputFile = expeDict["latexOutputFile"]
        end
    end

    if haskey(expeDict, "latexFormatPath")
        if typeof(expeDict["latexFormatPath"]) != String && !(typeof(expeDict["latexFormatPath"]) <: Vector)
            println("Warning: the latexFormatPath of the experiment must be a String or a vector of Strings (current value: ", expeDict["latexFormatPath"], " current type: ", typeof(expeDict["latexFormatPath"]), ")")
        else
            if !(typeof(expeDict["latexFormatPath"]) <: Vector)
                parameters.latexFormatPath = [expeDict["latexFormatPath"]]
            else
                parameters.latexFormatPath = expeDict["latexFormatPath"]
            end
        end
    end

    if haskey(expeDict, "minLatexCompilationInterval")
        if typeof(expeDict["minLatexCompilationInterval"]) != Int
            println("Warning: the minLatexCompilationInterval of the experiment must be an integer (current value: ", expeDict["minLatexCompilationInterval"], " current type: ", typeof(expeDict["minLatexCompilationInterval"]), ")")
        else
            parameters.minLatexCompilationInterval = expeDict["minLatexCompilationInterval"]
        end
    end

    if haskey(expeDict, "parametersToCombine")
        if !(typeof(expeDict["parametersToCombine"]) <: Dict)
            println("Warning: the parametersToCombine of the experiment must be a Dictionary (current value: ", expeDict["parametersToCombine"], " current type: ", typeof(expeDict["parametersToCombine"]), ")")
        else
            for (key, value) in expeDict["parametersToCombine"]
                if !(typeof(value) <:Vector)
                    expeDict["parametersToCombine"][key] = [value]
                end
            end
            parameters.parametersToCombine = expeDict["parametersToCombine"]
        end
    end

    return parameters
end

function addValuesToResults(expeJsonPath::String, addingFunctionJsonPath::String)

    # Read the experiment parameters in order to find the instances path
    parameters = readExpeFormat(expeJsonPath)

    # Get the function which get the additional result entries
    stringdata=join(readlines(addingFunctionJsonPath))
    functionDict = JSON.parse(stringdata)
    addingFunction = getMethods(functionDict)

    if length(addingFunction) != 1
        println("Warning: when adding entries in existing saved results, the json file must only contain one dictionary with an entry \"name\" and optionally an entry \"module\". However, ", length(addingFunction), " valid functions have been found from file ", addingFunctionJsonPath, ".")
    end

    addingFunction = addingFunction[1]

    for instancePath in parameters.instancesPath

        instanceName = splitext(basename(instancePath))[1]
        outputFile = parameters.outputPath * "/" * instanceName * ".json"
        outputFile = replace(outputFile, "//" => "/")

        # Variable which contains all the results of the current instance
        # Each result is a Dict which contains:
        # - one entry with the key "resolutionMethodName"
        # - one entry for each parameter
        # - one entry for each result
        instanceResults = Vector{Dict{String, Any}}()

        # Get the results previously computed for this instance if any
        if isfile(outputFile)
            println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS"), " Reading previous results from file \"", outputFile, "\"")
            stringdata=join(readlines(outputFile))

            try
                instanceResults = JSON.parse(stringdata)
                additionalResults = addingFunction(instancePath)

                # For each result of the current instance, add it the additional results
                for singleResult in instanceResults
                    merge!(singleResult, additionalResults)
                end

                # Save the updated results
                open(outputFile, "w") do fout
                    JSON.print(fout, instanceResults, 4)
                end

            catch e
                println(Dates.format(now(), "yyyy/mm/dd - HHhMM:SS") ,
                " Warning: Unable to read the JSON result file.")
            end

        end
    end
end

max_stars_repo_name: ZacharieALES/AutoExpe

max_issues_repo_name: ZacharieALES/AutoExpe

max_forks_repo_name: ZacharieALES/AutoExpe


text:
using Revise

using POMDPs
using ARDESPOT
using ProfileView
using POMDPModels

pomdp = BabyPOMDP()
bounds = IndependentBounds(reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)
solver = DESPOTSolver(epsilon_0=0.0,
                      K=100,
                      lambda=0.01,
                      bounds=bounds,
                      max_trials=3,
                      rng=MersenneTwister(4),
                      random_source=FastMersenneSource(1000, 1)
                     )
@show solver.lambda
p = solve(solver, pomdp)
b0 = initial_state_distribution(pomdp)
@show b0
@time D = ARDESPOT.build_despot(p, b0)
@show action(p, b0)
println(TreeView(D, 1, 90))

ext: jl


text:
#include("")
#positive definite
"""Matrix.jl
this is a stub
```input
Positive definite matrix #todo: how to define
-positive eigenvalues ( > 0)

```

```output:
`
2×2 Matrix{Int64}:
 5  2
 2  3

 """

A = [5 2; 2 3]
imp

# smoothFunction(A) #UncommentMe

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
# # Helmholtz equation
#
# In this example, we want to solve a (variant of) of the [Helmholtz equation](https://en.wikipedia.org/wiki/Helmholtz_equation).
# The example is inspired by an [dealii step_7](https://www.dealii.org/8.4.1/doxygen/deal.II/step_7.html) on the standard square.
#
# ```math
#  - \Delta u + u = f
# ```
#
# With boundary conditions given by
# ```math
# u = g_1 \quad x \in \Gamma_1
# ```
# and
# ```math
# n \cdot \nabla u = g_2 \quad x \in \Gamma_2
# ```
#
# Here Γ₁ is the union of the top and the right boundary of the square,
# while Γ₂ is the union of the bottom and the left boundary.
#
# ![](helmholtz.png)
#
# We will use the following weak formulation:
# ```math
# \int \nabla δu \cdot \nabla u d\Omega
# + \int δu \cdot u d\Omega
# - \int δu \cdot f d\Omega
# + \int δu \cdot (n \cdot \nabla u - g_2) d\Gamma_2
# ```
#
# where $δu$ is a suitable test function that satisfies:
# ```math
# δu = 0 \quad x \in \Gamma_1
# ```
# and $u$ is a suitable function that satisfies:
# ```math
# u = g_1 \quad x \in \Gamma_1
# ```
# The example highlights the following interesting features:
#
# * There are two kinds of boundary conditions, "Dirichlet" and "Von Neumann"
# * The example contains boundary integrals
# * The Dirichlet condition is imposed strongly and the Von Neumann condition is imposed weakly.
#
using Ferrite
using Tensors
using SparseArrays
using LinearAlgebra

const ∇ = Tensors.gradient
const Δ = Tensors.hessian;

grid = generate_grid(Quadrilateral, (150, 150))

dim = 2
ip = Lagrange{dim, RefCube, 1}()
qr = QuadratureRule{dim, RefCube}(2)
qr_face = QuadratureRule{dim-1, RefCube}(2)
cellvalues = CellScalarValues(qr, ip);
facevalues = FaceScalarValues(qr_face, ip);

dh = DofHandler(grid)
push!(dh, :u, 1)
close!(dh)

# We will set things up, so that a known analytic solution is approximately reproduced.
# This is a good testing strategy for PDE codes and known as the method of manufactured solutions.

function u_ana(x::Vec{2, T}) where {T}
    xs = (Vec{2}((-0.5,  0.5)),
          Vec{2}((-0.5, -0.5)),
          Vec{2}(( 0.5,  -0.5)))
    σ = 1/8
    s = zero(eltype(x))
    for i in 1:3
        s += exp(- norm(x - xs[i])^2 / σ^2)
    end
    return max(1e-15 * one(T), s) # Denormals, be gone
end;

dbcs = ConstraintHandler(dh)
# The (strong) Dirichlet boundary condition can be handled automatically by the Ferrite library.
dbc = Dirichlet(:u, union(getfaceset(grid, "top"), getfaceset(grid, "right")), (x,t) -> u_ana(x))
add!(dbcs, dbc)
close!(dbcs)
update!(dbcs, 0.0)

K = create_sparsity_pattern(dh);

function doassemble(cellvalues::CellScalarValues{dim}, facevalues::FaceScalarValues{dim},
                         K::SparseMatrixCSC, dh::DofHandler) where {dim}
    b = 1.0
    f = zeros(ndofs(dh))
    assembler = start_assemble(K, f)

    n_basefuncs = getnbasefunctions(cellvalues)
    global_dofs = zeros(Int, ndofs_per_cell(dh))

    fe = zeros(n_basefuncs) # Local force vector
    Ke = zeros(n_basefuncs, n_basefuncs) # Local stiffness mastrix

    @inbounds for (cellcount, cell) in enumerate(CellIterator(dh))
        fill!(Ke, 0)
        fill!(fe, 0)
        coords = getcoordinates(cell)

        reinit!(cellvalues, cell)
        # First we derive the non boundary part of the variation problem from the destined solution `u_ana`
        # ```math
        # \int \nabla δu \cdot \nabla u d\Omega
        # + \int δu \cdot u d\Omega
        # - \int δu \cdot f d\Omega
        # ```
        #+

        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            coords_qp = spatial_coordinate(cellvalues, q_point, coords)
            f_true = -LinearAlgebra.tr(hessian(u_ana, coords_qp)) + u_ana(coords_qp)
            for i in 1:n_basefuncs
                δu = shape_value(cellvalues, q_point, i)
                ∇δu = shape_gradient(cellvalues, q_point, i)
                fe[i] += (δu * f_true) * dΩ
                for j in 1:n_basefuncs
                    u = shape_value(cellvalues, q_point, j)
                    ∇u = shape_gradient(cellvalues, q_point, j)
                    Ke[i, j] += (∇δu ⋅ ∇u + δu * u) * dΩ
                end
            end
        end

        # Now we manually add the von Neumann boundary terms
        # ```math
        # \int δu \cdot (n \cdot \nabla u - g_2) d\Gamma_2
        # ```
        #+
        for face in 1:nfaces(cell)
            if onboundary(cell, face) &&
                   ((cellcount, face) ∈ getfaceset(grid, "left") ||
                    (cellcount, face) ∈ getfaceset(grid, "bottom"))
                reinit!(facevalues, cell, face)
                for q_point in 1:getnquadpoints(facevalues)
                    coords_qp = spatial_coordinate(facevalues, q_point, coords)
                    n = getnormal(facevalues, q_point)
                    g = gradient(u_ana, coords_qp) ⋅ n
                    dΓ = getdetJdV(facevalues, q_point)
                    for i in 1:n_basefuncs
                        δu = shape_value(facevalues, q_point, i)
                        fe[i] += -(δu * g) * dΓ
                        for j in 1:n_basefuncs
                            ∇u = shape_gradient(cellvalues, q_point, j)
                            Ke[i, j] += (δu * ∇u ⋅ n) * dΓ
                        end
                    end
                end
            end
        end

        celldofs!(global_dofs, cell)
        assemble!(assembler, global_dofs, fe, Ke)
    end
    return K, f
end;

K, f = doassemble(cellvalues, facevalues, K, dh);
apply!(K, f, dbcs)
u = Symmetric(K) \ f;

vtkfile = vtk_grid("helmholtz", dh)
vtk_point_data(vtkfile, dh, u)
vtk_save(vtkfile)
using Test #src
#src this test catches unexpected changes in the result over time.
#src the true maximum is slightly bigger then 1.0
@test maximum(u) ≈ 0.995280389173959 #src
@test u_ana(Vec{2}((-0.5, -0.5))) ≈ 1 #src
@test u_ana(Vec{2}((0.5, -0.5)))  ≈ 1 #src
@test u_ana(Vec{2}((-0.5, 0.5)))  ≈ 1 #src
println("Helmholtz successful")


text:
# -----------------------------------------------------------
#
# Construction of multidimensional DG derivative matrices
# based on the 1-D derivative matrix
#
# -----------------------------------------------------------

# Efficiency criticality: HIGH
# Matrix computations are only performed once,
# but this can be the main bottleneck if not done right

# Accuracy criticality: HIGH
# Critical for accurate PDE evolution

function D_matrix(d::Int, k::Int, n::Int, srefVD::Array{NTuple{3, CartesianIndex{D}}, 1},
        srefDV::Dict{NTuple{3, CartesianIndex{D}}, Int}; scheme = "sparse") where D
    D_matrix(Val(d), k, n, srefVD, srefDV, Val(Symbol(scheme)))
end

function D_matrix(::Val{d}, k::Int, n::Int, srefVD::Array{NTuple{3, CartesianIndex{D}}, 1},
        srefDV::Dict{NTuple{3, CartesianIndex{D}}, Int}, scheme::Val{Scheme}) where {d, D, Scheme}

    # len = length(srefVD[:,1])
    len = length(srefVD)
    @assert length(srefVD[:,1]) == length(srefVD)
    V2D_1D = V2Dref(Val(1),k,n,Val(:sparse)) # Why not 'scheme'?
    D2V_1D = D2Vref(Val(1),k,n,Val(:sparse))
    I = Int[]; J = Int[]; V = Float64[]

    # 1-dimensional derivative matrix - this will give all coefficient info
    Dmat_1D = periodic_DLF_matrix(k, n) #in hierarchischer Darstellung

    for j in 1:len
        lcm = srefVD[j]
        l = lcm[1][d]
        c = lcm[2][d]
        m = lcm[3][d]
        j_1D = D2V_1D[(CartesianIndex(l),CartesianIndex(c),CartesianIndex(m))]
        derivs = Dmat_1D[:, j_1D]::SparseVector{Float64,Int}

        for i_1D in derivs.nzind #1D Matrix auf Zeile der großen Matrix aufteilen (Linearisieren)
            lcm_1D = V2D_1D[i_1D]
            # using tuple constructor independent of Julia's metaprogramming:
            # make_cartesian_index(d, arr1, arr2) takes arr::CartesianIndex{1}
            # and makes a new CartesianIndex{D} using arr2::CartesianIndex{D}
            # with the dth value replaced by arr1[1]
            level2 = make_cartesian_index(d, lcm_1D[1], lcm[1])
            cutoff(scheme, level2, n) && continue
            cell2 = make_cartesian_index(d, lcm_1D[2], lcm[2])
            mode2 = make_cartesian_index(d, lcm_1D[3], lcm[3])
            i = srefDV[(level2, cell2, mode2)]
            push!(I, i)
            push!(J, j)
            push!(V, derivs[i_1D])
        end
    end
    # dropzeros! does not seem helpful for this matrix:
    return sparse(I, J, V, len, len, +)
end

function D_matrix(D::Int, d::Int, k::Int, n::Int; scheme="sparse")
    VD = V2Dref(Val(D), k, n, Val(Symbol(scheme)))
    DV = D2Vref(Val(D), k, n, Val(Symbol(scheme)))
    return D_matrix(d, k, n, VD, DV; scheme=scheme)
end

# Same functions as above but for anisotropic grids, changed to default scheme 'full'
function D_matrix(d::Int, k::Int, n::Vector{Int}, srefVD::Array{NTuple{3, CartesianIndex{D}}, 1},
                srefDV::Dict{NTuple{3, CartesianIndex{D}}, Int}; scheme = "full") where D
    D_matrix(Val(d), k, n, srefVD, srefDV, Val(Symbol(scheme)))
end

function D_matrix(::Val{d}, k::Int, n::Vector{Int}, srefVD::Array{NTuple{3, CartesianIndex{D}}, 1},
                srefDV::Dict{NTuple{3, CartesianIndex{D}}, Int}, scheme::Val{Scheme}) where {d, D, Scheme}

    # len = length(srefVD[:,1])
    # max_level = maximum(n)
    n_akt = n[d] # +1

    len = length(srefVD)
    @assert length(srefVD[:,1]) == length(srefVD)

    V2D_1D = V2Dref(Val(1),k,n_akt,Val(:full)) # Why not 'scheme'? changed from Val(:sparse) to Val(:full) # changed max_level to n
    D2V_1D = D2Vref(Val(1),k,n_akt,Val(:full)) # same here # changed max_level to n

    if n_akt == 0
        num_cells = 1
    else
        num_cells = 2^(n_akt)
    end
    # num_cells = 2^(n_akt) # -1
    # println("length(D2V_1D) = ", length(D2V_1D))
    # println("num_cells*k = ", num_cells*k)
    @assert length(D2V_1D) == num_cells*k #*(n_akt+1)

    I = Int[]; J = Int[]; V = Float64[]

    # 1-dimensional derivative matrix - this will give all coefficient info
    Dmat_1D = periodic_DLF_matrix(k, n_akt) #in hierarchischer Darstellung mit maximalem level # changed max_level to n

    len_akt = length(V2D_1D)^length(n)
    @assert length(V2D_1D[:,1]) == length(V2D_1D)

    for j in 1:len
        lcm = srefVD[j]
        # println(lcm)
        l = lcm[1][d]
        c = lcm[2][d]
        m = lcm[3][d]
        j_1D = D2V_1D[(CartesianIndex(l),CartesianIndex(c),CartesianIndex(m))]
        derivs = Dmat_1D[:, j_1D]::SparseVector{Float64,Int}

        for i_1D in derivs.nzind #1D Matrix auf Zeile der großen Matrix aufteilen (Linearisieren)
            lcm_1D = V2D_1D[i_1D]
            # using tuple constructor independent of Julia's metaprogramming:
            # make_cartesian_index(d, arr1, arr2) takes arr::CartesianIndex{1}
            # and makes a new CartesianIndex{D} using arr2::CartesianIndex{D}
            # with the dth value replaced by arr1[1]
            level2 = make_cartesian_index(d, lcm_1D[1], lcm[1])
            cutoff(scheme, level2, n) && continue ##cutoff mehrdimensional darstellen
            # cutoff der D dimensionales n übernimmt also ganzer vektor
            # aber trotzdem entweder true oder false zurückgeben
            cell2 = make_cartesian_index(d, lcm_1D[2], lcm[2])
            mode2 = make_cartesian_index(d, lcm_1D[3], lcm[3])
            i = srefDV[(level2, cell2, mode2)]
            push!(I, i)
            push!(J, j)
            push!(V, derivs[i_1D])
        end
    end
    # dropzeros! does not seem helpful for this matrix:
    return sparse(I, J, V, len, len, +)
end

# changed default scheme from sparse to full
function D_matrix(D::Int, d::Int, k::Int, n::Vector{Int}; scheme="full")
    # max_level = maximum(n) # alle D_matrizen sollen die gleich Größe haben
    # n_akt = n[d]
    VD = V2Dref(Val(D), k, n, Val(Symbol(scheme)))
    DV = D2Vref(Val(D), k, n, Val(Symbol(scheme)))
    return D_matrix(d, k, n, VD, DV; scheme=scheme)
end

function grad_matrix(D::Int, k::Int, n::Int; scheme="sparse")
    return [D_matrix(D, d, k, n; scheme=scheme) for d in 1:D]
end

function laplacian_matrix(D::Int, k::Int, n::Int; scheme="sparse")
    len = get_size(Val(D), k, n, Val(Symbol(scheme)))
    lap = spzeros(len, len)
    for i in 1:D
        D_op = D_matrix(D, i, k, n; scheme=scheme) ##
        lap += D_op * D_op
    end
    return lap
end

# Same function but for anisotropic grids
function laplacian_matrix(D::Int, k::Int, n::Vector{Int}; scheme="full")
    len = get_size(Val(D), k, n, Val(Symbol(scheme)))
    lap = spzeros(len, len)
    for i in 1:D
        D_op = D_matrix(D, i, k, n; scheme=scheme)
        lap += D_op * D_op
    end
    return lap
    # max_n = argmax(n)
    # len = get_size(Val(D), k, n[max_n], Val(Symbol(scheme)))
    # lap = spzeros(len, len)
    # for i in 1:D
    #     lap += laplacian_matrix(D,k,n[i],scheme=scheme)
    # end
    # return lap
    # lap_all = Float64[]
    # for i in 1:D
    #     len = get_size(Val(D), k, n[i], Val(Symbol(scheme)))
    #     lap = spzeros(len, len)
    #     D_op = D_matrix(D, i, k, n[i]; scheme=scheme) ##
    #     lap += D_op * D_op
    #     lap_all
    # end
    # return lap_all ##auch wieder liste
end


text: include("transferentropy.jl")


text:
if VERSION < v"0.4-"
    startswith = beginswith
end

const url_reg = r"^(([^:/?#]+):)?(//([^/?#]*))?([^?#]*)(\?([^#]*))?(#(.*))?"
const gh_path_reg_git=r"^/(.*)?/(.*)?.git$"

const releasejuliaver = v"0.4" # Current release version of Julia
const minjuliaver = v"0.4.0" # Oldest Julia version allowed to be registered
const minpkgver = v"0.0.1"   # Oldest package version allowed to be registered

print_list_3582 = false # set this to true to generate the list of grandfathered
                        # packages permitted under Issue #3582
list_3582 = Any[]
include("list_3582.jl")

# Issue 2064 - check that all listed packages have at least one tagged version
#2064## Uncomment the #2064# code blocks to generate the list of grandfathered
#2064## packages permitted
for pkg in readdir("METADATA")
    startswith(pkg, ".") && continue
    isfile(joinpath("METADATA", pkg)) && continue
    pkg in [
        "AuditoryFilters",
        "CorpusTools",
        "ErrorFreeTransforms",
        "Evapotranspiration",
        "GtkSourceWidget",
        "HiRedis",
        "KyotoCabinet",
        "LatexPrint",
        "MathLink",
        "ObjectiveC",
        "Processing",
        "RaggedArrays",
        "RationalExtensions",
        "SolveBio",
        "SortPerf",
    ] && continue
    if !("versions" in readdir(joinpath("METADATA", pkg)))
        #2064#println("        \"", pkg, "\","); continue
        error("Package $pkg has no tagged versions")
    end
end

# return julia version requirement for pkg, version
# if check is true, throw errors for metadata policy violations
function juliaver_in_require(pkg, version; check=true)
    requires_file = joinpath("METADATA", pkg, "versions", string(version), "requires")
    if !isfile(requires_file)
        if check
            error("File not found: $requires_file")
        else
            return v"0.0.0"
        end
    end
    open(requires_file) do f
        juliaver = v"0.0.0"
        hasjuliaver = false
        for line in eachline(f)
            if startswith(line, "julia")
                tokens = split(line)
                if length(tokens) <= 1
                    if check
                        error("$requires_file: oldest allowed julia version not specified (>= $minjuliaver needed)")
                    end
                else
                    juliaver = max(juliaver, convert(VersionNumber, tokens[2]))
                    hasjuliaver = true
                end
                if check && juliaver < minjuliaver
                    error("$requires_file: oldest allowed julia version $juliaver too old (>= $minjuliaver needed)")
                end
                if check && (juliaver < releasejuliaver && juliaver.patch==0 &&
                        (juliaver.prerelease != () || juliaver.build != ()))
                    # No prereleases older than current release allowed
                    error("$requires_file: prerelease $juliaver not allowed (>= $releasejuliaver needed)")
                end
            end
        end
        if check && !hasjuliaver
            error("$requires_file: no julia entry (>= $minjuliaver needed)")
        end
        return juliaver
    end
end

majmin(x::VersionNumber) = VersionNumber(x.major, x.minor, 0)

for (pkg, versions) in Pkg.Read.available()
    # Issue #2057 - naming convention check
    if endswith(pkg, ".jl")
        error("Package name $pkg should not end in .jl")
    end
    url = Pkg.Read.url(pkg)
    if length(versions) <= 0
        error("Package $pkg has no tagged versions.")
    end
    sortedversions = sort(collect(keys(versions)))
    maxver = sortedversions[end]

    m = match(url_reg, url)
    if m === nothing || length(m.captures) < 5
        error("Invalid url $url for package $pkg. Should satisfy $url_reg")
    end
    host = m.captures[4]
    if host === nothing
        error("Invalid url $url for package $pkg. Cannot extract host")
    end
    path = m.captures[5]
    if path === nothing
        error("Invalid url $url for package $pkg. Cannot extract path")
    end
    scheme = m.captures[2]
    if !(ismatch(r"git", scheme) || ismatch(r"https", scheme))
        error("Invalid url scheme $scheme for package $pkg. Should be 'git' or 'https'")
    end
    if ismatch(r"github\.com", host)
        m2 = match(gh_path_reg_git, path)
        if m2 == nothing
            error("Invalid GitHub url pattern $url for package $pkg. Should satisfy $gh_path_reg_git")
        end
        user = m2.captures[1]
        repo = m2.captures[2]

        # Issue #2057 - naming convention check
        if !endswith(repo, ".jl")
            error("Repository name $repo does not end in .jl")
        end
    end

    for (ver, avail) in versions
        # Check that all sha1 files have the correct version hashes
        sha1_file = joinpath("METADATA", pkg, "versions", string(ver), "sha1")
        if !isfile(sha1_file)
            error("File not found: $sha1_file")
        end
        sha1fromfile = open(readchomp, sha1_file)
        @assert sha1fromfile == avail.sha1

        # Issue #3582 - check that all versions of a package newer than the grandfathered
        # list in list_3582.jl are at least minpkgver and furthermore have a
        # requires file listing a minimum Julia version that is at least minjuliaver
        if print_list_3582 || !haskey(maxver_list_3582, pkg) || (ver > maxver_list_3582[pkg])
            try
                if ver < minpkgver
                    error("$pkg: version $ver no longer allowed (>= $minpkgver needed)")
                end
                # run with check=true for all versions more recent than grandfathered list
                juliaver = juliaver_in_require(pkg, ver; check=true)
                if ver == maxver
                    # check if minimum minor julia version has changed within the same
                    # minor package version (only check this for latest tag, it's okay
                    # for tags in old minor series to have too-high minimum requirements)
                    same_minor(x::VersionNumber) = (majmin(x) == majmin(ver) &&
                        juliaver_in_require(pkg, x; check=false) < juliaver)
                    ind_same_minor = findfirst(same_minor, sortedversions)
                    ind_same_minor == 0 && continue
                    first_same_minor = sortedversions[ind_same_minor]
                    juliaver_prev = juliaver_in_require(pkg, first_same_minor; check=false)
                    if majmin(juliaver) > majmin(juliaver_prev)
                        nextminor = VersionNumber(ver.major, ver.minor+1, 0)
                        error("New tag $ver of package $pkg requires julia $juliaver, ",
                            "but version $first_same_minor of $pkg requires julia ",
                            "$juliaver_prev. Use a new minor package version when support ",
                            "for an old version of Julia is dropped. Re-tag the package ",
                            "as $nextminor using `PkgDev.tag(\"$pkg\", :minor)`.")
                    end
                end
            catch err
                if print_list_3582
                    if isempty(list_3582) || list_3582[end][1] != pkg
                        push!(list_3582, (pkg, ver))
                    else
                        maxv = max(list_3582[end][2], ver)
                        list_3582[end] = (pkg, maxv)
                    end
                else
                    rethrow(err)
                end
            end
        end
    end
end
if print_list_3582
    sort!(list_3582, by=first)
    open(joinpath(dirname(@__FILE__), "list_3582.jl"), "w") do f
        println(f, """
# Issue #3582 - check that all versions of a package newer than the
# grandfathered list below are at least minpkgver and furthermore have a
# requires file listing a minimum Julia version that is at least minjuliaver
maxver_list_3582 = Dict([ # List of grandfathered packages""")
        for npkg in 1:length(list_3582)
            pkg, maxv = list_3582[npkg]
            println(f, """    ("$pkg", v"$maxv"),""")
        end
        println(f, "    ])")
    end
end

info("Checking that all entries in METADATA are recognized packages...")

# Scan all entries in METADATA for possibly unrecognized packages
const pkgs = [pkg for (pkg, versions) in Pkg.Read.available()]

for pkg in readdir("METADATA")
    # Traverse the 'versions' directory and make sure that we understand its contents
    # The only allowed subdirectories must be semvers and the only allowed
    # files within are 'sha1' and 'requires'
    #
    # Ref: #2040
    verinfodir = joinpath("METADATA", pkg, "versions")
    isdir(verinfodir) || continue # Some packages are registered but have no tagged versions. See #2064

    for verdir in readdir(verinfodir)
        version = try
            convert(VersionNumber, verdir)
        catch ArgumentError
            error("Invalid version number $verdir found in $verinfodir")
        end

        versions = Pkg.Read.available(pkg)
        if version in keys(versions)
           for filename in readdir(joinpath(verinfodir, verdir))
               if !(filename=="sha1" || filename=="requires")
                   relpath = joinpath(verinfodir, verdir, filename)
                   error("Unknown file $relpath encountered. Valid filenames are 'sha1' and 'requires'.")
               end
           end
        else
            relpath = joinpath("METADATA", pkg, "versions", verdir, "sha1")
            error("Version v$verdir of $pkg is not configured correctly. Check that $relpath exists.")
        end
    end
end

if haskey(ENV, "TRAVIS_PULL_REQUEST") && ENV["TRAVIS_PULL_REQUEST"] != "false"
    info("Checking repository tags...")
    include("check-pr.jl")
end

info("Verifying METADATA...")
if isdefined(Pkg.Entry, :check_metadata)
    Pkg.Entry.check_metadata()
else
    Pkg.add("PkgDev")
    import PkgDev
    PkgDev.Entry.check_metadata()
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct HelloRequest <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function HelloRequest(; kwargs...)
        obj = new(meta(HelloRequest), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        return obj
    end
end # mutable struct HelloRequest
const __meta_HelloRequest = Ref{ProtoMeta}()
function meta(::Type{HelloRequest})
    ProtoBuf.metalock() do
        if !isassigned(__meta_HelloRequest)
            __meta_HelloRequest[] = target = ProtoMeta(HelloRequest)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString]
            meta(
                target,
                HelloRequest,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_HelloRequest[]
    end
end
function Base.getproperty(obj::HelloRequest, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct HelloReply <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function HelloReply(; kwargs...)
        obj = new(meta(HelloReply), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        return obj
    end
end # mutable struct HelloReply
const __meta_HelloReply = Ref{ProtoMeta}()
function meta(::Type{HelloReply})
    ProtoBuf.metalock() do
        if !isassigned(__meta_HelloReply)
            __meta_HelloReply[] = target = ProtoMeta(HelloReply)
            allflds = Pair{Symbol,Union{Type,String}}[:message => AbstractString]
            meta(
                target,
                HelloReply,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_HelloReply[]
    end
end
function Base.getproperty(obj::HelloReply, name::Symbol)
    if name === :message
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

# service methods for Greeter
const _Greeter_methods = MethodDescriptor[
MethodDescriptor("SayHello", 1, HelloRequest, HelloReply)
] # const _Greeter_methods
const _Greeter_desc = ServiceDescriptor("helloworld.Greeter", 1, _Greeter_methods)

Greeter(impl::Module) = ProtoService(_Greeter_desc, impl)

mutable struct GreeterStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    GreeterStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_Greeter_desc, channel))
end # mutable struct GreeterStub

mutable struct GreeterBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    GreeterBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_Greeter_desc, channel))
end # mutable struct GreeterBlockingStub

SayHello(stub::GreeterStub, controller::ProtoRpcController, inp, done::Function) =
    call_method(stub.impl, _Greeter_methods[1], controller, inp, done)
SayHello(stub::GreeterBlockingStub, controller::ProtoRpcController, inp) =
    call_method(stub.impl, _Greeter_methods[1], controller, inp)

export HelloRequest, HelloReply, Greeter, GreeterStub, GreeterBlockingStub, SayHello


text:
function multiplyexistingvar(v::Vector{PolyVar{C}}, x::PolyVar{C}, i::Int) where {C}
    updatez = z -> begin
        newz = copy(z)
        newz[i] += 1
        newz
    end
    copy(v), updatez
end
function insertvar(v::Vector{PolyVar{C}}, x::PolyVar{C}, i::Int) where {C}
    n = length(v)
    I = 1:i-1
    J = i:n
    K = J.+1
    w = Vector{PolyVar{C}}(undef, n+1)
    w[I] = v[I]
    w[i] = x
    w[K] = v[J]
    updatez = z -> begin
        newz = Vector{Int}(undef, n+1)
        newz[I] = z[I]
        newz[i] = 1
        newz[K] = z[J]
        newz
    end
    w, updatez
end

include("cmult.jl")
include("ncmult.jl")

MP.multconstant(α, x::Monomial)   = Term(α, MA.mutable_copy(x))
MP.mapcoefficientsnz(f::Function, p::Polynomial) = Polynomial(map(f, p.a), MA.mutable_copy(p.x))
function MP.mapcoefficientsnz_to!(output::Polynomial, f::Function, t::MP.AbstractTermLike)
    MP.mapcoefficientsnz_to!(output, f, polynomial(t))
end
function MP.mapcoefficientsnz_to!(output::Polynomial, f::Function, p::Polynomial)
    resize!(output.a, length(p.a))
    @. output.a = f(p.a)
    Future.copy!(output.x.vars, p.x.vars)
    # TODO reuse the part of `Z` that is already in `output`.
    resize!(output.x.Z, length(p.x.Z))
    for i in eachindex(p.x.Z)
        output.x.Z[i] = copy(p.x.Z[i])
    end
    return output
end

# I do not want to cast x to TermContainer because that would force the promotion of eltype(q) with Int
function Base.:(*)(x::DMonomialLike, p::Polynomial)
    Polynomial(MA.mutable_copy(p.a), x*p.x)
end
function Base.:(*)(x::DMonomialLike{false}, p::Polynomial)
    # Order may change, e.g. y * (x + y) = y^2 + yx
    Polynomial(monovec(MA.mutable_copy(p.a), [x*m for m in p.x])...)
end
function Base.:(*)(p::Polynomial, x::DMonomialLike)
    Polynomial(MA.mutable_copy(p.a), p.x*x)
end

function _term_poly_mult(t::Term{C, S}, p::Polynomial{C, T}, op::Function) where {C, S, T}
    U = MA.promote_operation(op, S, T)
    if iszero(t)
        zero(Polynomial{C, U})
    else
        n = nterms(p)
        allvars, maps = mergevars([t.x.vars, p.x.vars])
        nv = length(allvars)
        # Necessary to annotate the type in case it is empty
        Z = Vector{Int}[zeros(Int, nv) for i in 1:n]
        for i in 1:n
            Z[i][maps[1]] = t.x.z
            Z[i][maps[2]] += p.x.Z[i]
        end
        Polynomial(op.(t.α, p.a), MonomialVector(allvars, Z))
    end
end
Base.:(*)(p::Polynomial, t::Term) = _term_poly_mult(t, p, (α, β) -> β * α)
Base.:(*)(t::Term, p::Polynomial) = _term_poly_mult(t, p, *)
_sumprod(a, b) = a * b + a * b
function _mul(::Type{T}, p::AbstractPolynomialLike, q::AbstractPolynomialLike) where T
    return _mul(T, polynomial(p), polynomial(q))
end
function _mul(::Type{T}, p::Polynomial{true}, q::Polynomial{true}) where T
    samevars = _vars(p) == _vars(q)
    if samevars
        allvars = copy(_vars(p))
    else
        allvars, maps = mergevars([_vars(p), _vars(q)])
    end
    N = length(p) * length(q)
    Z = Vector{Vector{Int}}(undef, N)
    a = Vector{T}(undef, N)
    i = 0
    for u in p
        for v in q
            if samevars
                z = u.x.z + v.x.z
            else
                z = zeros(Int, length(allvars))
                z[maps[1]] += u.x.z
                z[maps[2]] += v.x.z
            end
            i += 1
            Z[i] = z
            a[i] = u.α * v.α
        end
    end
    return allvars, a, Z
end
function Base.:(*)(p::Polynomial{true, S}, q::Polynomial{true, T}) where {S, T}
    if iszero(p) || iszero(q)
        zero(MA.promote_operation(*, typeof(p), typeof(q)))
    else
        polynomialclean(_mul(MA.promote_operation(MA.add_mul, S, T), p, q)...)
    end
end
function MA.mutable_operate_to!(p::Polynomial{false, T}, ::typeof(*), q1::MP.AbstractPolynomialLike, q2::MP.AbstractPolynomialLike) where T
    if iszero(q1) || iszero(q2)
        MA.mutable_operate!(zero, p)
    else
        ts = Term{false, T}[]
        MP.mul_to_terms!(ts, q1, q2)
        # TODO do better than create tmp
        tmp = polynomial!(ts)
        Future.copy!(p.a, tmp.a)
        Future.copy!(p.x.vars, tmp.x.vars)
        Future.copy!(p.x.Z, tmp.x.Z)
        return p
    end
end
function MA.mutable_operate_to!(p::Polynomial{true, T}, ::typeof(*), q1::MP.AbstractPolynomialLike, q2::MP.AbstractPolynomialLike) where T
    if iszero(q1) || iszero(q2)
        MA.mutable_operate!(zero, p)
    else
        polynomialclean_to!(p, _mul(T, q1, q2)...)
    end
end
function MA.mutable_operate!(::typeof(*), p::Polynomial{C}, q::Polynomial{C}) where C
    return MA.mutable_operate_to!(p, *, p, q)
end


text:
t1 = ([1], [0])
t2 = ([0], [1])
t3 = ([0], [-1, 1])
t4 = ([-1, 1], [0])

t5 = ([1:2], [0:1])

## pav(y) returns the isotonic regression of the predictor y.
## after W.J. Wilbur, L. Yeganova and W. Kim (2005)
function pav{T<:Real}(y::Vector{T})
    yy = convert(Vector{Float64}, y)
    N = length(y)
    F = [1:N]+1                 # forward pointer to intervals
    B = F-2                     # backward pointer to intervals
    Q1 = [1:N]                # queues
    Q2 = Int[]
    flag = false
    while length(Q1)>0
        m = 1
        while length(Q1)>0
            k = shift!(Q1)
            if m <= k
                flag=false
                println(k, " ", F[k])
                while F[k] != N && yy[k]> yy[F[k]]
                    ## pool intervals beginning at k and F[k]
                    first = k>0 ? F[k-1] : 1
                    pooled = vcat(first:F[k], F[F[k]-1]:F[F[k]])
                    println("pool ", pooled')
                    yy[pooled] = mean(y[pooled])
                    u = F[F[k]]
                    F[k] = u
                    if u<=N B[u] = k end
                    m = u
                    flag = true
                end
            end
            if flag && B[k]>=0
                push!(Q2, B[k])
            end
        end
        h = Q2; Q2 = Q1; Q1 = h
    end
    yy
end


text:
using Pkg; Pkg.activate(".")
using Luxor, Colors

function chemellia_old(loc::Point, size::Real; petalcolor="plum3", period=16pi, step=0.5)
    for pt in reverse(spiral(size, 1.3, period=period, stepby=step))
        rad = distance(loc, loc+pt)
        blendstart = between(loc, loc+pt, 0.7)
        blendend = between(loc, loc+pt, 1.4)
        setblend(blend(blendstart, blendend, petalcolor, "white"))
        circle(loc+pt, size+0.37*rad, :fill)
    end
end

# thanks to cormullion!
function chemellia(loc::Point=O;size::Real=110, npetals=105, petalcolor="maroon")
    origin(loc)
    #sethue(petalcolor)
    setline(0.75)
    #squircle(O, 2*size, 2*size, :fill, rt = 0.3)
    #for n in 110:-1:5
    local totalrot = 0
    for n in range(size, 5, length=npetals)
        rotate(deg2rad(137.5))
        totalrot = totalrot+137.5
        petal = makebezierpath([ngon(O, n, 3, 0, vertices=true)...], smoothing=1.2)
        @layer begin
            translate(n, 0)
            setblend(blend(O, O + (n, 0),  petalcolor, "white"))
            drawbezierpath(petal, :fill)
            sethue("grey40")
            drawbezierpath(petal, :stroke)
        end
    end
    rotate(deg2rad(-totalrot%360))
end

function bee(loc::Point, size::Real; aspect=1.6, angle=0)
    rotate(angle)
    ellipse(loc, aspect*size, size, :clip)
    sethue("gold")
    for f=-0.5aspect:0.4*aspect:0.5*aspect
        rect(loc+f.*(size,0)-(0,0.5*size), aspect/5*size, size, :fill)
    end
    sethue("black")
    for f=-0.3*aspect:0.4*aspect:0.3*aspect
        rect(loc+f.*(size,0)-(0,0.5*size), aspect/5*size, size, :fill)
    end
    clipreset()
    f1 = loc + 0.4*size .* (0,-2.5)
    f2 = loc + 0.4*size .* (-1.5,-1.5)
    f3 = loc + 0.4*size .* (1.5,-1.5)
    sethue("lightgrey")
    ellipse(f1, f2, 0.3*size, :fill)
    ellipse(f1, f3, 0.3*size, :fill)
end

# size in inches
function make_squarelogo(fname="squarelogo.png", size=1.5, dpi=300)
    np = size * dpi
    Drawing(np, np, fname)
    origin()
    sethue("olivedrab")
    setline(0.75) # scale this?
    squircle(O, 0.48*np, 0.48*np, :fill, rt = 0.3)
    chemellia(O + Point(0.5*np, 0.5*np), size=0.22*np)
    finish()
end

function make_biglogo(fname="./biglogo.eps", width=4.875, dpi=300)
    np_w = round(width * dpi)
    np_h = round(np_w / 3.25)
    Drawing(np_w, np_h, fname)
    origin()
    fontface("Josefin Sans Semibold")
    fontsize((11.0/80) * np_w)
    setcolor("olivedrab")
    text("Chemellia", Point(-.23*np_w,0), valign=:middle)
    bee(Point(0.266*np_w, -0.458*np_h), 0.024*np_w, aspect=1.5, angle=π/8)
    chemellia(Point(0.1375*np_w, 0.5*np_h), size=0.0563*np_w)
    finish()
end

make_squarelogo()
make_squarelogo("./squarelogo.eps")
make_biglogo()
make_biglogo("./biglogo.png")


text:
"Robust approximate comparison operator for two reals."
function ≅(a::Real, b::Real; ϵscale=1000, ϵpow=0.5)
    abs(a-b) < ϵscale*(max(eps(a),eps(b))^ϵpow)*(1+abs(a)+abs(b))
end

"""
Draw `n` random `θ` coefficients. Second parameter adds additional
dispersion on the log scale.
"""
rand_θ(n, log_σ=4) = normalize!(randn(n) .* [exp(log_σ*randn()) for _ in n], 2)

"""
Accuracy of various f(θ, x), comparing the original and BigFloat
values at various `xs`. Report relative accuracy in digits, as quantiles q
"""
function accuracy(f, θ, xs, q=[0.0, 0.1, 0.5])
    y_bigf = f.([BigFloat.(θ)], BigFloat.(xs))
    y = f.([θ], xs)         # FIXME remove [] in v0.6
    acc(a,b) = -Float64(log10(abs(a-b)/(1+abs(Float64(a)))))
    quantile(acc.(y_bigf, y), q)
end

function test_univariate(fam, expected_degf, expected_domain;
                         xs = nothing, test_basis = true, f = nothing,
                         test_partial = true, test_valuepartial = test_partial)
    @test degf(fam) == expected_degf
    @test domain(fam) == expected_domain
    @test all(p ∈ expected_domain for p in points(fam))

    # when xs given, test evaluation using basis
    if xs ≠ nothing
        θ = rand(expected_degf)
        for x in xs
            if test_basis
                @test evaluate(fam, θ, x) == dot(basis(fam, x), θ)
            else
                @test isa(evaluate(fam, θ, x), typeof(x))
            end
        end
    end

    # when f given, test approximation
    if f ≠ nothing
        zs = points(fam)
        ys = f.(zs)
        θ0 = zeros(degf(fam))
        fit!(fam, f, θ0)
        θ1 = fit(fam, f)

        @test θ0 == θ1

        for z in zs
            @test evaluate(fam, θ1, z) ≅ f(z)
        end

        if test_basis
            θ2 = basis_matrix(fam, zs) \ ys
            @test θ1 ≈ θ2
        end
    end

    # when f given, test function fitting
    if f ≠ nothing
        zs = points(fam)
        pf = fitfun(fam, f.(zs))

        @test points(pf) == points(fam)
        @test degf(pf) == degf(fam)
        @test domain(pf) == domain(fam)
        @test family(pf) == fam

        for z in zs
            fz = f(z)
            f′z = derivative(f, z)
            @test pf(z) ≅ fz
            if test_partial
                @test pf(Partial(z)) ≅ f′z
            end
            if test_valuepartial
                v, p = pf(ValuePartial(z))
                @test v ≅ fz
                @test p ≅ f′z
            end
        end
    end
end


text:
using DataStructures
using JLSponge
using Test
const SegmentArrives_Result_OK = true
const SegmentArrives_Result_NOT_SYN = false
const DEFAULT_TEST_WINDOW = UInt(137)

function build_seg(; data="", ack=false, rst=false, syn=false, fin=false,
                   seqno=WrappingInt32(0), ackno=WrappingInt32(0), win=UInt16(0))
    return seg = TCPSegment(TCPHeader(; ack, rst, syn, fin, seqno, ackno, win), data)
    return seg
end

function segment_arrives(receiver::TCPReceiver; data="", with_ack=nothing, with_syn=false,
                         with_fin=false, seqno=0, result=nothing)
    seg = build_seg(; data, syn=with_syn, fin=with_fin, seqno=WrappingInt32(seqno),
                    ack=with_ack !== nothing,
                    ackno=WrappingInt32(with_ack !== nothing ? with_ack : 0))
    segment_received!(receiver, seg)
    res = ackno(receiver) === nothing ? SegmentArrives_Result_NOT_SYN :
          SegmentArrives_Result_OK

    result !== nothing && @test result == res
end

function expect_seg(segments::Queue{TCPSegment};
    no_flag = nothing,
    ack=nothing, rst=nothing, syn=nothing, fin=nothing,
    payload_size=nothing, seqno=nothing, ackno=nothing,
    data=nothing, win = nothing)
    if no_flag !== nothing
        ack = isnothing(ack) ? false : ack
        rst = isnothing(rst) ? false : rst
        syn = isnothing(syn) ? false : syn
        fin = isnothing(fin) ? false : fin
    end

    isempty(segments) && error("No segs")
    seg = dequeue!(segments)
    hd = seg.header
    !isnothing(ack) && @test hd.ack == ack
    !isnothing(rst) && @test hd.rst == rst
    !isnothing(syn) && @test hd.syn == syn
    !isnothing(fin) && @test hd.fin == fin
    !isnothing(seqno) && @test hd.seqno == seqno
    !isnothing(ackno) && @test hd.ackno == ackno

    !isnothing(win) && @test hd.win == win
    !isnothing(payload_size) && @test length(seg.payload) == payload_size
    !isnothing(data) && @test seg.payload.storage[] == data
    seg
end

function expect_seg(sender::TCPSender;
    no_flag = nothing,
    ack=nothing, rst=nothing, syn=nothing, fin=nothing,
    payload_size=nothing, seqno=nothing, ackno=nothing,
    data=nothing, win = nothing)

    expect_seg(sender.segments_out; no_flag, ack, rst, syn, fin, payload_size, seqno, ackno, data, win)
end

function expect_seg(conn::TCPConnection;
    no_flag = nothing,
    ack=nothing, rst=nothing, syn=nothing, fin=nothing,
    payload_size=nothing, seqno=nothing, ackno=nothing,
    data=nothing, win = nothing)

    expect_seg(conn.segments_out; no_flag, ack, rst, syn, fin, payload_size, seqno, ackno, data, win)
end

expect_no_seg(sender::TCPSender) = @test isempty(sender.segments_out)

function ack_received_test(sender::TCPSender, ackno, win_size::Integer=UInt16(137))
    ack_received!(sender, ackno, UInt16(win_size))
    fill_window!(sender)
end

function write_bytes!(sender::TCPSender, bytes="", end_input=false)
    write!(stream_in(sender), bytes)
    end_input && end_input!(stream_in(sender))
    fill_window!(sender)
end

function expect_state(sender::TCPSender, state)
    @test state_summary(sender) == state
end

sender_close(sender::TCPSender) = write_bytes!(sender, "", true)

fsm_in_listen(;cap::Int=64000, retx_timeout::UInt16=UInt16(1000), fixed_isn=nothing) =
    TCPConnection(;cap, retx_timeout, fixed_isn)

function fsm_in_syn_sent(fixed_isn; cap::Int=64000, retx_timeout::UInt16=UInt16(1000))
   conn = TCPConnection(;cap, retx_timeout, fixed_isn)
   connect!(conn)
   expect_one_seg(conn; no_flag=true, syn=true, seqno=fixed_isn, payload_size=0)
   conn
end

function fsm_in_established(tx_isn::WrappingInt32=WrappingInt32(0),
    rx_isn::WrappingInt32=WrappingInt32(0); args...)

    conn = fsm_in_syn_sent(tx_isn; args...)
    send_syn!(conn, rx_isn, tx_isn + 1)
    expect_one_seg(conn; no_flag=true, ack=true, ackno=rx_isn + 1, payload_size=0)
    conn
end

function fsm_in_close_wait(tx_isn::WrappingInt32=WrappingInt32(0),
    rx_isn::WrappingInt32=WrappingInt32(0))

    conn = fsm_in_established()
    send_fin!(conn, rx_isn + 1, tx_isn + 1)
    expect_one_seg(conn, no_flag=true, ack=true, ackno=rx_isn+2)
    conn
end

function fsm_in_last_ack(tx_isn::WrappingInt32=WrappingInt32(0),
    rx_isn::WrappingInt32=WrappingInt32(0))

    conn = fsm_in_close_wait(rx_isn, rx_isn)
    conn |> end_input_stream!
    expect_one_seg(conn, no_flag=true, fin=true, ack=true, seqno=rx_isn+1, ackno=tx_isn+2)
    conn
end

function fsm_in_fin_wait_1(tx_isn::WrappingInt32=WrappingInt32(0),
    rx_isn::WrappingInt32=WrappingInt32(0))

    conn = fsm_in_established(tx_isn, rx_isn)
    conn |> end_input_stream!
    expect_one_seg(conn, no_flag=true, fin=true, ack=true, ackno=rx_isn+1, seqno=tx_isn+1)
    conn
end

function fsm_in_fin_wait_2(tx_isn::WrappingInt32=WrappingInt32(0),
    rx_isn::WrappingInt32=WrappingInt32(0))

    conn = fsm_in_fin_wait_1(tx_isn, rx_isn)
    send_ack!(conn, rx_isn + 1, tx_isn + 2)
    conn
end

function fsm_in_closing(tx_isn::WrappingInt32=WrappingInt32(0),
    rx_isn::WrappingInt32=WrappingInt32(0))

    conn = fsm_in_fin_wait_1(tx_isn, rx_isn)
    send_fin!(conn, rx_isn + 1, tx_isn + 1)
    expect_one_seg(conn, no_flag=true, ack=true, ackno=rx_isn+2)
    conn
end

function fsm_in_time_wait(tx_isn::WrappingInt32=WrappingInt32(0),
    rx_isn::WrappingInt32=WrappingInt32(0))

    conn = fsm_in_fin_wait_1(tx_isn, rx_isn)
    send_fin!(conn, rx_isn + 1, tx_isn + 2)
    expect_one_seg(conn, no_flag=true, ack=true, ackno=rx_isn+2)
    conn
end

function expect_one_seg(conn::TCPConnection;
    no_flag=nothing,
    ack=nothing, rst=nothing, syn=nothing, fin=nothing,
    payload_size=nothing, seqno=nothing, ackno=nothing,
    data=nothing, win = nothing)

    seg = expect_seg(conn.segments_out; no_flag, ack, rst, syn, fin, payload_size, seqno, ackno, data, win)
    @test conn.segments_out |> isempty
    seg
end

function send_syn!(conn::TCPConnection, seqno::WrappingInt32, ackno=nothing)
    ackno === nothing ?
    send_seg!(conn; syn=true, seqno) :
    send_seg!(conn; syn=true, seqno, ack=true, ackno)
end

function send_seg!(conn::TCPConnection;
    ack=false, rst=false, syn=false, fin=false,
    seqno=WrappingInt32(0), ackno=WrappingInt32(0),
    data="", win = UInt16(137))

    seg = build_seg(;data, ack, rst, syn, fin, seqno, ackno, win)
    segment_received!(conn, seg)
    nothing
end

function send_seg!(conn::TCPConnection, seg::TCPSegment)
    hdr = seg.header
    seg = build_seg(;data=seg.payload.storage[],
        ack=hdr.ack, rst=hdr.rst, syn=hdr.syn,
        fin=hdr.fin, seqno=hdr.seqno,
        ackno=hdr.ackno, win=hdr.win)

    segment_received!(conn, seg)
    nothing
end

function send_ack!(conn::TCPConnection, seqno::WrappingInt32, ackno::WrappingInt32, swin=DEFAULT_TEST_WINDOW)
    win = swin
    send_seg!(conn; ack=true, win, seqno, ackno)
end

function send_rst!(conn::TCPConnection, seqno::WrappingInt32, ackno=nothing)
    ackno === nothing ?
    send_seg!(conn; seqno, rst=true) :
    send_seg!(conn; seqno, rst=true, ackno=ackno, ack=true)
end

function send_fin!(conn::TCPConnection, seqno::WrappingInt32, ackno=nothing)
    ackno === nothing ?
    send_seg!(conn; seqno, fin=true) :
    send_seg!(conn; seqno, fin=true, ack=true, ackno=ackno)
end

function send_byte!(conn::TCPConnection,
    seqno::WrappingInt32,
    ackno::Union{WrappingInt32, Nothing}, val::Char)

    ackno === nothing ?
    send_seg!(conn; seqno, data=string(val)) :
    send_seg!(conn; seqno, data=string(val), ackno=ackno, ack=true)
end

function send_data!(conn::TCPConnection,
    seqno::WrappingInt32,
    ackno::WrappingInt32, data::AbstractString)
    send_seg!(conn; ack=true, ackno, data, seqno)
end

expect_state(conn::TCPConnection, state::JLSponge.State) = @test TCPState(conn) == TCPState(state)

expect_notin_state(conn::TCPConnection, state::JLSponge.State) = @test TCPState(conn) != TCPState(state)

expect_no_seg(conn::TCPConnection) = @test isempty(conn.segments_out)

function expect_data!(conn::TCPConnection, data::Nothing=nothing)
    bytes_avail = inbound_stream(conn) |> buffer_size
    actual_data = read!(inbound_stream(conn), bytes_avail)
end

function expect_data!(conn::TCPConnection, data::String)
    actual_data = expect_data!(conn)
    @test data == actual_data
end

function expect_no_data(conn::TCPConnection)
    @test inbound_stream(conn) |> buffer_size == 0
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using FractionalDiffEq, SpecialFunctions
using Test

@testset "FractionalDiffEq.jl" begin
    include("FODETests.jl")
    include("FDDETests.jl")
    include("FPDETests.jl")
    include("FDifferenceTests.jl")
    include("FIETests.jl")
    include("DODETests.jl")
    include("auxillary.jl")
    include("models.jl")
end


text:
convert(::Type{String}, s::T) where {T<:Union{TAMQPShortStr,TAMQPLongStr,TAMQPByteArray}} = String(copy(s.data))
convert(::Type{Bool}, b::TAMQPBit) = Bool(b.val & 0x1)

simplify(val::T) where {T <: Union{TAMQPShortStr,TAMQPLongStr,TAMQPByteArray}} = String(copy(val.data))
simplify(val::TAMQPFieldArray) = [simplify(elem) for elem in val.data]
simplify(table::TAMQPFieldTable) = Dict{String,Any}(simplify(f.name)=>simplify(f.val) for f in table.data)
simplify(val::TAMQPFieldValue) = simplify(val.fld)
simplify(x) = x


text:
mutable struct ScalarExponentiator <: Exponentiator
end
function ScalarExponentiator(v::AbstractVector)
    size(v,1) == 1 || error("Scalar exponentiator only applicable for 1×1 operators")
    ScalarExponentiator()
end

function (SE::ScalarExponentiator)(Ω::LinearMap{T}, τ::Number,
                                   v, w) where T<:Number
    a = Ω([one(T)])[1]
    w[:] = exp(τ*a)*v
end

export ScalarExponentiator


text:
draw(PNG("portrait.png", 100, 140),
     Paper("lightgrey")) do
    move(0, 0)
    line(1, 0)
    move(0, 0)
    line(0, 1)
end


text:
module MdfReader
import InteractiveUtils.subtypes
export MDF
export read

struct CHAR
    x::UInt8
end
struct nCHAR{n}
    chars::NTuple{n,CHAR}
end
function Base.convert(::Type{NTuple{N,CHAR}},str::AbstractString) where N
    @assert N==length(str) "string must have length $N"
    NTuple{N,CHAR}(CHAR(c) for c in str)
end
nCHAR(str::AbstractString) = nCHAR{length(str)}(str)
Base.Char(c::CHAR) = Base.Char(c.x)

Base.show(io::IO, c::CHAR) = print(io, Char(c))
Base.string(c::nCHAR{x}) where {x} = join(Char.(c.chars))
Base.show(io::IO, c::nCHAR{x}) where {x} = print(io, string(c))
function Base.convert(::Type{nCHAR{N}},str::AbstractString) where {N}
    nCHAR{length(str)}(str)
end

struct BOOL
    x::UInt16
end
struct REAL
    x::Float64
end
struct LINK
    target::UInt32
end
isnill(b::LINK) = iszero(b.target)
Base.convert(::Type{LINK}, x::Int64) = LINK(x)
Base.show(io::IO, b::LINK) = print(io, "byteNr: $(b.target)")
Base.seek(io::IO, b::LINK) = seek(io, b.target)

struct BlockHeader
    type::nCHAR{2}
    size::UInt16
    function BlockHeader(t,s)
        @assert isBlockType(t) "$t is not a valid block type"
        new(t,s)
    end
end

Base.string(hd::BlockHeader) = string(hd.type)
typestr(hd) = uppercase(string(hd)) * "BLOCK"

abstract type BLOCK end
header(b::BLOCK) = b.header

struct IDBLOCK <: BLOCK
    fileIdentifier::nCHAR{8}
    formatIdentifier::nCHAR{8}
    programIdentifier::nCHAR{8}
    defaultByteOrder::UInt16
    defaultFloatingPointFormat::UInt16
    versionNumber::UInt16
    codePageNumber::UInt16
    reserved::nCHAR{28}
    standardFlags::UInt16
    customFlags::UInt16
    IDBLOCK() = new()
end

struct HDBLOCK <: BLOCK
    header::BlockHeader
    dataGroupBlock::LINK
    fileComment::LINK
    programBlock::LINK
    numberDataGroups::UInt16
    date::nCHAR{10}
    time::nCHAR{8}
    author::nCHAR{32}
    organization::nCHAR{32}
    project::nCHAR{32}
    subject::nCHAR{32}
    timeStamp::UInt64
    utcTimeOffset::Int16
    timeQuality::UInt16
    timerIdentification::nCHAR{32}
    HDBLOCK() = new()
end

struct DGBLOCK <: BLOCK
    header::BlockHeader
    nextDataGroupBlock::LINK
    firstChannelGroupBlock::LINK
    triggerBlock::LINK
    dataBlock::LINK
    numberChannelGroups::UInt16
    numverRecordIDs::UInt16
    Reserved::UInt32
    DGBLOCK() = new()
end

struct CGBLOCK <: BLOCK
    header::BlockHeader
    nextChannelGroupBlock::LINK
    firstChannelBlock::LINK
    comment::LINK
    recordID::UInt16
    numberChannels::UInt16
    sizeDataRecord::UInt16
    numberRecords::UInt32
    firstSampleReductionBlock::LINK
    CGBLOCK() = new()
end

struct CNBLOCK <: BLOCK
    header::BlockHeader
    nextChannelBlock::LINK
    conversionFormula::LINK
    sourceDependingExtension::LINK
    dependencyBlock::LINK
    comment::LINK
    channelType::UInt16
    shortSignalName::nCHAR{32}
    signalDescription::nCHAR{128}
    startOffset::UInt16
    numberOfBits::UInt16
    signalDataType::UInt16
    valueRangeValid::UInt16
    minimumSignalValue::REAL
    maximumSignalValue::REAL
    samplingRate::REAL
    longSignalName::LINK
    displayName::LINK
    additionalByteOffset::UInt16
    CNBLOCK() = new()
end

struct TXBLOCK <: BLOCK
    header::BlockHeader
    text::String
    TXBLOCK() = new()
end
struct PRBLOCK <: BLOCK
    header::BlockHeader
    text::String
    PRBLOCK() = new()
end
function Base.read(io::IO, b::Union{TXBLOCK,PRBLOCK})
    hd = read(io, BlockHeader)
    nb = size(hd) - sizeof(hd)
    str = String(Char.(read(io, nb)))
end

# TODO: make read block such that it only reads the blocksize number of
# fields and makes the other fields without content. (e.g. zero)

linkfields(b::DataType) =
    [fieldname(b, i) for i in 1:fieldcount(b) if fieldtype(b, i) == LINK]
linkfields(b) = linkfields(typeof(b))
type(hd) = eval(Symbol(typestr(hd)))
isBlockType(hd) = isdefined(@__MODULE__,Symbol(typestr(hd))) && (type(hd) in subtypes(BLOCK))

"""
	linkoffset(b)

returns the link field offset for the given block type
"""
function linkoffset(b::Type{T}) where {T<:BLOCK}
    [fieldoffset(b, i) for i in 1:fieldcount(b) if fieldtype(b, i) == LINK]
end

function readlink(io, pos)
    seek(io, pos)
    read(io, LINK)
end
function readlinks(io, hd)
    sz = sizeof(hd)
    loffs = linkoffset(type(hd))
    pos = loffs .+ position(io) .- sz
    links = [readlink(io, p) for p in pos]
end
Base.size(hd::BlockHeader) = hd.size

type(b::BLOCK) = type(header(b))

Base.read(io::IO, d::BlockHeader) = Base.read(io,type(d),size(d))

function Base.read(io::IO, d::Type{T}, size=packed_sizeof(d)) where {T<:Union{BLOCK,BlockHeader}}
    a = Type{d}()
    ptr = Ptr{UInt8}(pointer_from_objref(a))
    unsafe_copyto!()
    read!(io,a)
    a[]
end
zero
function zero!(p::Ptr{T}) where T
    nb = sizeof(T)
    puint = Ptr{UInt8}(p)
    unsafe_copyto!()
end
Base.show(io::IO, b::BlockHeader) = print(io, "Block
      type: $(typestr(b))
      size: $(size(b))
    ")

## next functions for iterating over linked blocks
next(b::DGBLOCK) = b.nextDataGroupBlock
next(b::CGBLOCK) = b.nextChannelGroupBlock
next(b::CNBLOCK) = b.nextChannelBlock
first(b::HDBLOCK) = b.dataGroupBlock
first(b::DGBLOCK) = b.firstChannelGroupBlock
first(b::CGBLOCK) = b.firstChannelBlock

struct MDF
    identifier
    header
    dataGroups
    channelGroups
    channels
end
function readBlock(io,T)
    hd = read(io, BlockHeader)
    skip(io, -sizeof(hd))
    @assert isBlockType(hd) "blocktype $(hd.type) is not valid"
    sz = size(hd)
    b  = read(io, type(hd), sz) # call a sized block reader
end
function readBlock(io,T, link::LINK)
    seek(io, link)
    readBlock(io,T)
end
function readGroup(io, start, T)
    b = readBlock(io,T, start)
    blocks = Array{T,1}(undef, 0)
    push!(blocks,b)
    nextBlock = next(b)
    while ~isnill(nextBlock)
        b = readBlock(io,T, nextBlock)
        push!(blocks, b)
        nextBlock = next(b)
    end
    blocks
end

function MDF(filepath)
    open(filepath) do io
        id = read(io, IDBLOCK)
        hd = read(io, HDBLOCK)
        dgblocks = readGroup(io, first(hd), DGBLOCK)
        cgblocks = [readGroup(io, first(dg), CGBLOCK) for dg in dgblocks]
        cnblocks = [readGroup(io, first(cg), CNBLOCK) for cgs in cgblocks for cg in cgs]
        MDF(id, hd, dgblocks, cgblocks, cnblocks)
    end
end

end


text:
State = Int64
Action = Int64
DPTable = Array{Dict{State, Float64}, 1}

StatesAvailable = Dict(
    1 => [1],
    2 => [2, 3, 4],
    3 => [5, 6],
    4 => [7]
)

ActionsAvailble = Dict(
    (1, 1) => [2, 3, 4],
    (2, 2) => [5],
    (2, 3) => [5, 6],
    (2, 4) => [5, 6],
    (3, 5) => [7],
    (3, 6) => [7],
    (4, 7) => [7]
)

Transitions = Dict(
    (1, 2) => (2, 7),
    (1, 3) => (3, 8),
    (1, 4) => (4, 5),
    (2, 5) => (5, 12),
    (3, 5) => (5, 8),
    (3, 6) => (6, 9),
    (4, 5) => (5, 7),
    (4, 6) => (6, 13),
    (5, 7) => (7, 9),
    (6, 7) => (7, 6)
#    (7, 7) => (7, 0)
)

EndState = Dict(7 => 0.0)

function state_func(state)
    return StatesAvailable[state]
end

function actions_func(stage, state)
    return ActionsAvailble[(stage, state)]
end

function transitions(state, action)
    return Transitions[(state, action)]
end

function is_better(minimize::Bool, x::Float64, y::Float64)
    if minimize
        return x < y
    end
    return x > y

end

function compute_solution_by_br(current_stage::Int64,
        minimize::Bool,
        table::DPTable,
        state_func,
        actions_func,
        transition)
    if(current_stage == 0)
        return reverse(table)
    end
    factor = 1

    if !minimize
        factor = -1
    end

    states = state_func(current_stage)
    table_next = table[end]
    table_curr = Dict{State, Float64}()
    for state ∈ states
        table_curr[state] = factor * Inf
        actions = actions_func(current_stage, state)
        for action ∈ actions
            transition = transitions(state, action)
            new_state = transition[1]
            val = transition[2]
            if new_state ∈ keys(table_next)
                value_of_action = table_next[new_state] + val
                if is_better(minimize, value_of_action, table_curr[state])
                    table_curr[state] = value_of_action
                end
            end
        end
    end
    push!(table, table_curr)
    return compute_solution_by_br(current_stage - 1, minimize,
        table, state_func, actions_func, transition)
end

table = [Dict(7 => 0.0)]
compute_solution_by_br(3, true, table, state_func, actions_func, transitions)


text:
using GeometryTypes
using SailingPolars
using SailingNavigation
using Test

#=
The data is selected so the best time path is slightly longer than the
shortest length path. The forbidden regions are x, representing land or reef.
The allowed sailing points are . and start and finish are S and F.

x  .  .  F  .  .  x  .  x
.  .  .  .  .  .  .  x  x
x  .  .  x  x  x  .  .  .
.  .  x  x  x  x  .  x  x
x  .  .  .  x  x  .  x  .
x  .  .  .  x  x  .  x  .
.  .  .  .  x  .  .  x  .
x  .  .  .  .  .  .  x  .
.  .  .  S  .  .  .  .  .
=#
const forbidden = Point2{Int}.([
    [1, 8], [2, 1], [2, 8], [3, 5], [3, 8], [4, 1], [4, 5], [4, 6], [4, 8], [5, 1],
    [5, 5], [5, 6], [5, 8], [6, 3], [6, 4], [6, 5], [6, 6], [6, 8], [6, 9], [7, 1],
    [7, 4], [7, 5], [7, 6], [8, 8], [8, 9], [9, 1], [9, 7], [9, 9],
])

# Create regional wind patterns on the map.
function surfacebylongitude(lon)
    return lon < -155.03 ? SurfaceParameters(-5.0, 8, 150, 0.5) :
           lon < -155.99 ? SurfaceParameters(-90.0, 20, 150, 0.4) :
                           SurfaceParameters(180.0, 25, 150, 0.3)
end

# Vary wind speeds over time.
function mutatetimeslices!(slices)
    for (i, slice) in enumerate(slices), x in slice
        x.sp = SurfaceParameters(x.sp.winddeg, x.sp.windkts * (1 + 0.002 * i),
            x.sp.currentdeg, x.sp.currentkts)
    end
end

const startpos = Point2{Int}(1, 4)
const endpos = Point2{Int}(9, 4)
const pmat  = [Position(19.78 - 1/60 + i/60, -155.0 - 5/60 + j/60) for i in 0:8, j in 0:8]
const gpoints = map(pt -> GridPoint(pt, surfacebylongitude(lon(pt))), pmat)
const slices = [deepcopy(gpoints) for _ in 1:200]
mutatetimeslices!(slices)

const routeprob = RoutingProblem(10.0, slices, forbidden, 1, startpos, endpos, false)
const filename = "polar.csv"
const sp = getpolardata(filename)
const tp = minimumtimeroute(routeprob, sp)

@test tp.duration ≈ 4.0 * 60 + 43.69787966871
@test tp.path == Point{2,Int64}[[1, 4], [1, 5], [2, 6], [3, 7], [4, 7], [5, 7], [6, 7], [7, 7], [8, 6], [8, 5], [9, 4]]


text:
using MATLAB, Test

is_ci() = lowercase(get(ENV, "CI", "false")) == "true"

if !is_ci() # only test if not CI

    @testset "jl.eval" begin

        mat"""
            $result = jl.eval('2+2');
        """
        @test result == 4

        mat"""
            [$s, $c] = jl.eval('sin(pi/3), cos(pi/3)');
        """
        @test s ≈ sin(pi/3)
        @test c ≈ cos(pi/3)

        mat"""
            jleval 1 + 1;
            $result = ans;
        """
        @test result == 2
    end

    @testset "jl.call" begin
        mat"""
            $result = jl.call('factorial', int64(10));
        """
        @test result == factorial(10)
    end

    @testset "jl.mex" begin
        mat"""
        a = rand(5,5);
        jleval import MATLAB;
        jleval double_it(args::Vector{MATLAB.MxArray}) = [2*MATLAB.jvalue(arg) for arg in args];
        result = jl.mex('double_it', a);
        $result = result{1};
        $expected_result = 2*a;
        """
        @test result ≈ expected_result
    end

end

max_stars_repo_name: taylormcd/mexjulia

max_issues_repo_name: taylormcd/mexjulia

max_forks_repo_name: taylormcd/mexjulia


text:
using AMQPClient
const VIRTUALHOST = "/"
const HOST = "127.0.0.1"

function receive()
    # 1. Create a connection to the localhost or 127.0.0.1 of virtualhost '/'
    connection(; virtualhost=VIRTUALHOST, host=HOST) do conn
        # 2. Create a channel to send messages
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan
            # 3. Declare a exchange
            exchange = "direct_logs"
            exchange_declare(chan, exchange, EXCHANGE_TYPE_DIRECT)

            result, queue_name, _, _  = queue_declare(chan, "", exclusive=true)

            # 4. Receive queues to bind
            if length(Base.ARGS) <= 0
                println(Base.stdout, "Usage: [info] [warning] [error]\n")
                Base.exit(1)
            end

            # 4.1 Bind queues
            for severity in Base.ARGS[1:end]
                queue_bind(chan, queue_name, exchange,
                           severity)
            end

            println(" [*] Waiting for messages. To exit press CTRL+C")
            on_receive = (msg) -> begin
                data = String(msg.data)
                println("Received the message: $data")
                basic_ack(chan, msg.delivery_tag)
            end

            success, consumer_tag = basic_consume(chan, queue_name, on_receive)

            while true
                sleep(1)
            end
        end
    end
end

# Don't exit on Ctrl-C
Base.exit_on_sigint(false)
try
    receive()
catch ex
    if ex isa InterruptException
        println("Interrupted")
    else
        println("Exception: $ex")
    end
end

max_forks_repo_name: tomashq/rabbitmq-tutorials