text:
# aes256cbc_tests.jl

import AES256CBC
using Base.Test

test_pwd = "Secret Passphrase"
test_salt = "a3e550e89e70996c" # must be random UBytes by genRandUBytes(8)
expected_k = "e299ff9d8e4831f07e5323913c53e5f0fec3a040a211d6562fa47607244d0051"
expected_i = "7c7ed9434ddb9c2d1e1fcc38b4bf4667" # initial value

test_data = [
  ( # test short
    "Message", # auto padding
    "da8aab1b904205a7e49c1ecc7118a8f4",
    "da8aab1b904205a7e49c1ecc7118a8f4", # last iv
  ),( # test long
    "Message\t\t\t\t\t\t\t\t\t", # with self padding
    "da8aab1b904205a7e49c1ecc7118a8f4804bef7be79216196739de7845da182d",
    "804bef7be79216196739de7845da182d", # last iv
  )]

@test AES256CBC.string2bytes("abC") == AES256CBC.UBytes([0x61, 0x62, 0x43])
println(@sprintf "genPadding(0) == [%s]" AES256CBC.genPadding(0))
@test AES256CBC.genPadding(1) == AES256CBC.UBytes([1])
@test AES256CBC.genPadding(5) == Array{UInt8, 1}([5, 5, 5, 5, 5])
@test AES256CBC.genPadding(9) == Array{UInt8, 1}([9, 9, 9, 9, 9, 9, 9, 9, 9])
println(@sprintf "genPadding(16) == [%s]" AES256CBC.genPadding(16))
salt = AES256CBC.genRandUBytes(8)
println(@sprintf "genRandUBytes(8) == [%s]" bytes2hex(salt))
salt = AES256CBC.genRandUBytes(8)
println(@sprintf "genRandUBytes(8) == [%s]" bytes2hex(salt))
salt = AES256CBC.genRandUBytes(8)
println(@sprintf "genRandUBytes(8) == [%s]" bytes2hex(salt))
salt = hex2bytes(test_salt) # replace for test
println(@sprintf "salt[%s]" bytes2hex(salt))
(key32, iv16) = AES256CBC.genKey32Iv16(AES256CBC.string2bytes(test_pwd), salt)
println(@sprintf "key32[%s]\niv16[%s]" bytes2hex(key32) bytes2hex(iv16))
@test key32 == hex2bytes(expected_k)
@test iv16 == hex2bytes(expected_i)

for (i, (test_plain, expected_enc, expected_liv)) in enumerate(test_data)
  println(i)
  # i == 1 && begin println("*** skip ***"); continue end # skip short test

  encoded = AES256CBC.encryptAES256CBC(key32, iv16,
    AES256CBC.string2bytes(test_plain))
  println(@sprintf "encoded[%s]" bytes2hex(encoded))
  @test encoded == hex2bytes(expected_enc)
  decoded = AES256CBC.decryptAES256CBC(key32, iv16, encoded)
  println(@sprintf "decoded[%s]" bytes2hex(decoded))
  @test decoded == AES256CBC.string2bytes(test_plain)
  println(@sprintf "decoded[%s]" bytestring(decoded))
  @test bytestring(decoded) == test_plain
end


text:
# Simple fetch single image from MongoDB store, using Neo4j lookup

using Caesar

# include(joinpath(dirname(@__FILE__),"..","database","blandauthremote.jl"))
# addrdict["session"] = "SESSTURT45"
# cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)
cloudGraph, addrdict = standardcloudgraphsetup()

session = addrdict["session"]

@show IDs = getPoseExVertexNeoIDs(
  cloudGraph.neo4j.connection,
  sessionname=session,
  reqbackendset=false
);

vid, cvid = IDs[2]

# cv = CloudGraphs.get_vertex(cloudGraph, cvid, false )
cv = CloudGraphs.get_vertex(cloudGraph, cvid, true  )

img = Caesar.getBigDataElement(cv, "keyframe_rgb")
imf = open("test.png","w")
write(imf, img.data)
close(imf)
run(`eog test.png`)
Base.rm("test.png")

rgb = ImageMagick.readblob(img.data);

# using JSON
# using LibBSON
#
# mongk = JSON.parse(cv.properties["mongo_keys"])
# BSONOID(mongk["keyframe_rgb"])

#

ext: jl


text:
function request(api::API, request_method, endpoint;
            auth = AnonymousAuth(), handle_error = true,
            headers = Dict(), params = Dict(), allow_redirects = true,
            query = Dict())
    authenticate_headers!(headers, auth)
    headers["Accept"] = "application/json"
    api_endpoint = api_uri(api, endpoint)
    if request_method == Requests.get
        params["minorversion"] = "4"
        r = request_method(api_endpoint; headers = headers, query = params, allow_redirects = allow_redirects)
    else
        query["minorversion"] = 4
        r = request_method(api_endpoint; headers = headers, query = query,
                                            json = params, allow_redirects = allow_redirects)
    end
    handle_error && handle_response_error(r)
    return r
end

hget(api::API, endpoint; options...) = request(api, Requests.get, endpoint; options...)
post(api::API, endpoint; options...) = request(api, Requests.post, endpoint; options...)

##################
# Error Handling #
##################

function handle_response_error(r::HttpCommon.Response)
    if r.status >= 400
        dj = Requests.json(r)
        # Seem like qbo isn't consistent about this
        data = haskey(dj, "Fault") ? dj["Fault"] : dj["fault"]
        message = "Error found in QBO reponse:\n"*
                    "\tStatus Code: $(r.status)\n"
        for error in (haskey(data, "Error") ? data["Error"] : data["error"])
            message *= """
            \tMessage: $(strip(haskey(error, "Message") ? error["Message"] : error["message"]))\n
            \tDetail: $(strip(haskey(error, "Detail") ? error["Detail"] : error["detail"]))\n
            """
        end
        error(message)
    end
end


text:
include("../examples/rijke.jl")
using PyPlot
using JLD
function plot_sensitivity()
	vsh, dJds = test_lss()
	fig, ax = subplots(1,1)
	fig1, ax1 = subplots(1,1)
	vsh = vsh[:,100:end,10]'
	n = 0.005*(axes(vsh)[1] .- 1.)
	n_samples = size(dJds)[1]
	dJds_avg = sum(dJds)/n_samples
	var_dJds = sum(x -> x^2, dJds .- dJds_avg)/n_samples
	ax.plot(n, vsh[:,1], label=L"v_{\rm sh,x}")
	ax.plot(n, vsh[:,2], label=L"v_{\rm sh,y}")
	ax.plot(n, vsh[:,3], label=L"v_{\rm sh,z}")
	ax.set_xlabel("time",fontsize=18)
	ax.set_ylabel("components of the shadowing direction",
					  fontsize=18)
	ax.xaxis.set_tick_params(labelsize=18)
	ax.yaxis.set_tick_params(labelsize=18)
	fig.legend(fontsize=18)
	ax1.set_ylabel("sensitivities from NILSS",fontsize=18)
	ax1.set_xlabel("sample number",fontsize=18)
	ax1.set_title("trajectory length = 5000", fontsize=18)
	ax1.xaxis.set_tick_params(labelsize=18)
	ax1.yaxis.set_tick_params(labelsize=18)
	ax1.grid(true)
	ax1.errorbar(x=1:n_samples,y=dJds,yerr=sqrt(var_dJds),
			 linestyle="none",ms=4)
	ax1.plot(sum(dJds)/n_samples*ones(n_samples),"--")
	#mean value: 0.905328466, variance = 0.042, without
	# time dilation
end
function plot_condition_number()
	dJds, condnum = test_condition_number()
	fig, ax = subplots(1,1)
	n_samples = size(dJds)[1]
	n_arr = StepRange(500, 155, 5000)
	ax.plot(n_arr, condnum, ".", ms=10.0)
	ax.set_xlabel("trajectory length",fontsize=18)
	ax.set_ylabel("condition number",
					  fontsize=18)
	ax.xaxis.set_tick_params(labelsize=18)
	ax.yaxis.set_tick_params(labelsize=18)
end
function plot_da_err(z_prd,z_obs,fig, ax)
	n_trj, n_gd, n_exps = size(z_prd)
	z_opt_prd = view(z_prd, :, n_gd, :)
	err = abs.((z_obs .- z_opt_prd)./z_obs)
	max_err = findmax(err, dims=1)[1][:]
	max_err = log.(max_err)
	max_err = max_err .- minimum(max_err)
	norm_max_err = max_err/maximum(max_err)
	clr_map = plt.get_cmap("winter")
	clr_acc_err = clr_map(norm_max_err)
	times = dt*LinRange(0,n_trj,n_trj)
	times = times
	err = err
	for i = 1:n_exps
		ax.semilogy(times, err[:,i], ".",
					ms=0.5,color=clr_acc_err[i,:])
	end
	mean_err = (sum(err, dims=2)/n_exps)[:]
	std_err = sum((err .- mean_err).^2.0, dims=2)/n_exps
	std_err = sqrt.(std_err[:])
	ax.semilogy(times, mean_err, ".", color="k",ms=2.0)
end
function read_error_files()
	n_times_1 = 4
	n_exps_t = 15
	n_times_2 = 8
	n_exps = (n_times_1 + n_times_2)*n_exps_t
	filenames = Array{String,1}(undef, n_exps)
	for i = 1:n_times_1
		for j = 1:n_exps_t
			filenames[(i-1)*n_exps_t + j] = string(
			"../data/rijke_DA/1-to-60/rijke_exp",
			string(i+1), "_", string(j), ".jld")
		end
	end
	n_cltd = n_times_1*n_exps_t
	for i = 1:n_times_2
		for j = 1:n_exps_t
			filenames[n_cltd + (i-1)*n_exps_t + j] = string(
			"../data/rijke_DA/61-to-180/rijke_exp",
			string(i+1), "_", string(j), ".jld")
		end
	end

	data = load(filenames[1])
	prd = data["z_prd"]
	obs = data["z_obs"]
	n_trj, n_gd, x = size(prd)
	z_prd = zeros(n_trj, n_gd, n_exps)
	z_obs = zeros(n_trj, n_exps)
	z_prd[:,:,1] = prd[:,:,1]
	z_obs[:,1] = obs[:,1]
	for k = 2:n_exps
		data = load(filenames[k])
		prd = data["z_prd"]
		obs = data["z_obs"]
		z_prd[:,:,k] = prd[:,:,1]
		z_obs[:,k] = obs[:,1]
	end

	fig, ax = subplots(1,1)

	ax.xaxis.set_tick_params(labelsize=24)
	ax.yaxis.set_tick_params(labelsize=24)
	ax.set_xlabel("time", fontsize=24)
	ax.set_ylabel("Prediction error", fontsize=24)
	ax.grid(true)

	plot_da_err(z_prd,z_obs,fig, ax)

end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
using SparseDiffTools
using SparseArrays
using VertexSafeGraphs
using LightGraphs
using Random

const LG = LightGraphs

n = rand(10:100)
matrices = Array{SparseMatrixCSC, 1}(undef, 0)

for _ in 1:20
    n = rand(10:100)
    matrix = sprand(Int64, n, n, 0.5)
    push!(matrices, copy(matrix))
end

for i in 1:20
    matrix = matrices[i]
    g = matrix2graph(matrix, false)
    for e in edges(g)
        src = LG.src(e)
        dst = LG.dst(e)
        col1 = abs.(matrix[:, src])
        col2 = abs.(matrix[:, dst])
        pr = col1' * col2
        @test pr != 0
    end
end

for i in 1:20
    matrix = matrices[i]
    g = matrix2graph(matrix, true)
    for e in edges(g)
        src = LG.src(e)
        dst = LG.dst(e)
        row1 = abs.(matrix[src, :])
        row2 = abs.(matrix[dst, :])
        pr = row1' * row2
        @test pr != 0
    end
end


text:
using LinearAlgebra

# Linear Combination
a=[1,2]; b=[3,4];
alpha = -0.5; beta = 1.5;
c = alpha*a + beta*b

# Define function
function lincomb(coeff, vectors)
    n = length(vectors[1])
    a = zeros(n);
    for i = 1:length(vectors)
        a = a + coeff[i] * vectors[i];
    end
return a
end

#Alternate definition
function lincomb(coeff, vectors)
    return sum( coeff[i] * vectors[i]
        for i = 1:length(vectors) )
    end

# Use funtion
lincomb( ( -0.5, 1.5), ( [1, 2], [ 3, 4]) )

max_stars_repo_name: lynnlangit/julia-linear-algebra

max_issues_repo_name: lynnlangit/julia-linear-algebra

max_forks_repo_name: lynnlangit/julia-linear-algebra


text:
using Unitful
using SampledSignals

export aba_, ab, stimulus, audiospect_stimulus

const sr = 8000

silence(len) = SampleBuf(zeros(SampledSignals.inframes(Int,len,sr)),sr)
function tone(freq,len)
  freq_ = SampledSignals.inHz(freq)
  len_ = SampledSignals.inseconds(len)
  SampleBuf(sin.(2π.*freq_.*range(0,stop=len_,step=1/sr)), sr)
end

function ramp(xs,len)
  n = SampledSignals.inframes(Int,len,samplerate(xs))
  t = range(0,stop=1,length=n)
  xs[1:n] .*= cos.(π.*(t.-1)./2)
  xs[end-n+1:end] .*= cos.(π.*t./2)
  xs
end

normpower(x) = x ./ sqrt.(mean(x.*x,dims=1))
amplify(ratio) = x -> x.*ratio

const stim_cache = Dict{Vector{Union{Number,String}},AbstractMatrix}()
function audiospect_stimulus(params,settings;cache=false)
  settings = read_settings(settings)
  params = read_params(params)
  if cache
    get!(stim_cache,[[params[:Δt],params[:Δf],params[:f]];
                             collect(values(settings.stimulus))]) do
      audiospect_stimulus_(params,settings)
    end
  else
    audiospect_stimulus_(params,settings)
  end
end

function audiospect_stimulus_(params,settings)
  stim = stimulus(params[:Δt],params[:f],params[:Δf];
                  settings.stimulus...) |> normpower |> amplify(-10dB)
  @info "Stimulus is $(maximum(domain(stim))) seconds long."
  as = Audiospect(;settings.freqs.analyze...)
  filt(as,stim,false)
end

function stimulus(total_len,freq,delta;repeats=10,tone_len_fraction=0.5,
                  pattern="ab",ramp_len=0ms)
  if pattern == "ab"
    ab(tone_len_fraction*total_len,(1-tone_len_fraction)*total_len,1,
       repeats,freq,delta;ramp_len=asseconds(ramp_len))
  elseif pattern == "aba_"
    aba_(tone_len_fraction*total_len,(1-tone_len_fraction)*total_len,
         repeats,freq,delta;ramp_len=asseconds(ramp_len))
  else
    error("Unexpected stimulus pattern '$pattern'")
  end
end

function aba_(tone_len,spacing_len,repeats,freq,delta;ramp_len=0ms)
  a_freq=freq
  b_freq=freq*(2^(delta/12))
  a = tone(a_freq,tone_len)
  b = tone(b_freq,tone_len)
  a = ramp_len > 0ms ? ramp(a,25ms) : a
  b = ramp_len > 0ms ? ramp(b,25ms) : b

  aba_ = [a;silence(spacing_len);
          b;silence(spacing_len);
          a;silence(2spacing_len+tone_len)]

  aba_seq = silence(0s)
  for i = 1:repeats
    aba_seq = [aba_seq; aba_]
  end

  aba_seq
end

function ab(tone_len,spacing_len,offset_ratio,repeats,freq,delta;ramp_len=0ms)
  @assert 0 <= offset_ratio <= 1

  a_freq=freq
  b_freq=freq*(2^(delta/12))
  a = tone(a_freq,tone_len)
  b = tone(b_freq,tone_len)
  a = ramp_len > 0ms ? ramp(a,25ms) : a
  b = ramp_len > 0ms ? ramp(b,25ms) : b

  offset = silence(offset_ratio*(tone_len+spacing_len))
  ab = [a;offset].+[offset;b]
  ab_len = 2(tone_len+spacing_len)
  ab = [ab; silence(ab_len - duration(ab)*s)]

  ab_seq = silence(0s)
  for i = 1:repeats
    ab_seq = [ab_seq; ab]
  end

  ab_seq
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..")) # add InteratomicPotentials to environment stack

using InteratomicPotentials
using Documenter
using DocumenterCitations
using Literate

DocMeta.setdocmeta!(InteratomicPotentials, :DocTestSetup, :(using InteratomicPotentials); recursive = true)

bib = CitationBibliography(joinpath(@__DIR__, "citations.bib"))

# Generate examples

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR = joinpath(@__DIR__, "src/generated")

examples = Pair{String,String}[]

for (_, name) in examples
    example_filepath = joinpath(EXAMPLES_DIR, string(name, ".jl"))
    Literate.markdown(example_filepath, OUTPUT_DIR, documenter = true)
end

examples = [title => joinpath("generated", string(name, ".md")) for (title, name) in examples]

makedocs(bib;
    modules = [InteratomicPotentials],
    authors = "CESMIX-MIT",
    repo = "https://github.com/cesmix-mit/InteratomicPotentials.jl/blob/{commit}{path}#{line}",
    sitename = "InteratomicPotentials.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://cesmix-mit.github.io/InteratomicPotentials.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "InteratomicPotentials Interface" => "interface.md",
        "Examples" => examples,
        "API Reference" => "api.md",
    ],
    doctest = true,
    linkcheck = true,
    strict = true
)

deploydocs(;
    repo = "github.com/cesmix-mit/InteratomicPotentials.jl",
    devbranch = "main",
    push_preview = true
)


text:
module LaplaceBeltramiCubedSphereTestsSeq

  using Test
  using Gridap
  using GridapGeosciences
  using Plots

  include("../ConvergenceAnalysisTools.jl")
  include("../LaplaceBeltramiCubedSphereTests.jl")

  n_values=generate_n_values(2)
  hs=[2.0/n for n in n_values]
  model_args_series=[(n,) for n in n_values]

  @time ahs1,ak1errors,as1=convergence_study(solve_laplace_beltrami,hs,model_args_series,1,8)
  @test round(as1,digits=1) ≈ 1.0

  @time ahs2,ak2errors,as2=convergence_study(solve_laplace_beltrami,hs,model_args_series,2,12)
  @test round(as2,digits=1) ≈ 2.0

end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
import .ExpressionExplorer: SymbolsState, FunctionName

"Information container about the cells to run in a reactive call and any cells that will err."
struct TopologicalOrder
	"Cells that form a directed acyclic graph, in topological order."
	runnable::Array{Cell,1}
	"Cells that are in a directed cycle, with corresponding `ReactivityError`s."
	errable::Dict{Cell,ReactivityError}
end

"Return a `TopologicalOrder` that lists the cells to be evaluated in a single reactive run, in topological order. Includes the given roots."
function topological_order(notebook::Notebook, topology::NotebookTopology, roots::Array{Cell,1}; allow_multiple_defs=false)::TopologicalOrder
	entries = Cell[]
	exits = Cell[]
	errable = Dict{Cell,ReactivityError}()

	function dfs(cell::Cell)
		if cell in exits
			return
		elseif haskey(errable, cell)
			return
		elseif length(entries) > 0 && entries[end] == cell
			return # a cell referencing itself is legal
		elseif cell in entries
			currently_in = setdiff(entries, exits)
			cycle = currently_in[findfirst(isequal(cell), currently_in):end]
			for cell in cycle
				errable[cell] = CyclicReferenceError(topology, cycle...)
			end
			return
		end

		push!(entries, cell)
		assigners = where_assigned(notebook, topology, cell)
		if !allow_multiple_defs && length(assigners) > 1
			for c in assigners
				errable[c] = MultipleDefinitionsError(topology, c, assigners)
			end
		end
		referencers = where_referenced(notebook, topology, cell) |> Iterators.reverse
		for c in (allow_multiple_defs ? referencers : union(assigners, referencers))
			if c != cell
				dfs(c)
			end
		end
		push!(exits, cell)
	end

	# we first move cells to the front if they call `import` or `using`
    # we use MergeSort because it is a stable sort: leaves cells in order if they are in the same category
    prelim_order_1 = sort(roots, alg=MergeSort, by=c -> cell_precedence_heuristic(topology, c))
	# reversing because our search returns reversed order
	prelim_order_2 = Iterators.reverse(prelim_order_1)
	dfs.(prelim_order_2)
	ordered = reverse(exits)
	TopologicalOrder(setdiff(ordered, keys(errable)), errable)
end

function disjoint(a::Set, b::Set)
	!any(x in a for x in b)
end

"Return the cells that reference any of the symbols defined by the given cell. Non-recursive: only direct dependencies are found."
function where_referenced(notebook::Notebook, topology::NotebookTopology, myself::Cell)::Array{Cell,1}
	to_compare = union(topology[myself].definitions, topology[myself].funcdefs_without_signatures)
	where_referenced(notebook, topology, to_compare)
end
"Return the cells that reference any of the given symbols. Non-recursive: only direct dependencies are found."
function where_referenced(notebook::Notebook, topology::NotebookTopology, to_compare::Set{Symbol})::Array{Cell,1}
	return filter(notebook.cells) do cell
		!disjoint(to_compare, topology[cell].references)
	end
end

"Return the cells that also assign to any variable or method defined by the given cell. If more than one cell is returned (besides the given cell), then all of them should throw a `MultipleDefinitionsError`. Non-recursive: only direct dependencies are found."
function where_assigned(notebook::Notebook, topology::NotebookTopology, myself::Cell)::Array{Cell,1}
	self = topology[myself]
	return filter(notebook.cells) do cell
		other = topology[cell]
		!(
			disjoint(self.definitions,                 other.definitions) &&

			disjoint(self.definitions,                 other.funcdefs_without_signatures) &&
			disjoint(self.funcdefs_without_signatures, other.definitions) &&

			disjoint(self.funcdefs_with_signatures,    other.funcdefs_with_signatures)
		)
	end
end

"Return whether any cell references the given symbol. Used for the @bind mechanism."
function is_referenced_anywhere(notebook::Notebook, topology::NotebookTopology, sym::Symbol)::Bool
	any(notebook.cells) do cell
		sym ∈ topology[cell].references
	end
end

"Return whether any cell defines the given symbol. Used for the @bind mechanism."
function is_assigned_anywhere(notebook::Notebook, topology::NotebookTopology, sym::Symbol)::Bool
	any(notebook.cells) do cell
		sym ∈ topology[cell].definitions
	end
end

"""Assigns a number to a cell - cells with a lower number might run first.

This is used to treat reactive dependencies between cells that cannot be found using static code anylsis."""
function cell_precedence_heuristic(topology::NotebookTopology, cell::Cell)::Number
	if :LOAD_PATH ∈ topology[cell].references
		# https://github.com/fonsp/Pluto.jl/issues/323
		1
	elseif !isempty(cell.module_usings)
		# always do `using X` before other cells, because we don't (yet) know which cells depend on it (we only know it with `import X` and `import X: y, z`)
		2
	elseif :include ∈ topology[cell].references
		# https://github.com/fonsp/Pluto.jl/issues/193
		# because we don't (yet) know which cells depend on it
		3
	else
		4
	end
end

const md_and_friends = [Symbol("@md_str"), Symbol("@html_str")]

"""Does the cell only contain md"..." and html"..."?

This is used to run these cells first."""
function is_just_text(topology::NotebookTopology, cell::Cell)::Bool
	# https://github.com/fonsp/Pluto.jl/issues/209
	isempty(topology[cell].definitions) && isempty(topology[cell].funcdefs_with_signatures) &&
		length(topology[cell].references) <= 2 &&
		topology[cell].references ⊆ md_and_friends
end


text:
#=
Contains info specific to finite volume such as
- cell enters
- face centers
- interpolation weights
-

Also has functions useful for FV such as
- cell to node / node to cell interpolation
-
=#

struct FVInfo
    cellCenters::Array{Float64}     # Coordinates of cell centers
    faceCenters::Array{Float64}     # Coordinates of face centers

    cell2node::Array{Array{Int64,1}}            # Neighboring cell indices for each node
    cell2nodeWeight::Array{Array{Float64,1}}    # Cell to node interpolation weights
end

etypetonv = [2, 3, 4, 4, 8, 6, 5, 2, 3, 4, 4, 8, 6, 5, 1, 4, 8, 6, 5]; # number of vertices for each type
etypetodim= [1, 2, 2, 3, 3, 3, 3, 1, 2, 2, 3, 3, 3, 3, 1, 2, 2, 3, 3]; # dimension of each type
etypetonf = [2, 3, 4, 4, 6, 5, 5, 2, 3, 4, 4, 6, 5, 5, 1, 4, 6, 5, 5]; # number of faces for element types
etypetoftype=[0,1, 1, 2, 3, 3, 3, 0, 1, 1, 2, 3, 3, 3, 0, 0, 0, 0, 0]; # type of faces for this element type

# Build the FVInfo from a given grid
function build_FV_info(grid)
    nel = size(grid.loc2glb, 2);
    nface = size(grid.face2glb, 3);
    nnode = size(grid.allnodes, 2);
    dim = size(grid.allnodes, 1);

    cell_centers = zeros(dim, nel);
    face_centers = zeros(dim, nface);
    cell2node = Array{Array{Int64,1},1}(undef,nnode);
    cell2nodeWeight = Array{Array{Float64,1},1}(undef,nnode);

    # Find cell centers
    for ei=1:nel
        center = zeros(dim);
        nvtx = size(grid.loc2glb, 1);
        for ni=1:nvtx
            center .+= grid.allnodes[:, grid.loc2glb[ni, ei]];
        end
        cell_centers[:, ei] = center ./ nvtx;
    end

    # Find face centers
    for fi=1:nface
        center = zeros(dim);
        nvtx = size(grid.face2glb, 1);
        for ni=1:nvtx
            center .+= grid.allnodes[:, grid.face2glb[ni, 1, fi]];
        end
        face_centers[:, fi] = center ./ nvtx;
    end

    # Build cell2node map
    # Check all elements to find ones touching this node.
    nvtx = size(grid.loc2glb, 1);
    for ni=1:nnode
        cell2node[ni] = zeros(Int,0);
        for ei=1:nel
            for nj=1:nvtx
                if grid.loc2glb[nj, ei] == ni
                    # This cell touches this node
                    push!(cell2node[ni], ei);
                    break;
                end
            end
        end
    end

    # Compute cell2node interpolation weights
    for ni=1:nnode
        cell2nodeWeight[ni] = zeros(length(cell2node[ni]));
        sumweight = 0;
        for ei=1:length(cell2node[ni])
            dist = norm(grid.allnodes[:,ni] - cell_centers[:, cell2node[ni][ei]]);
            sumweight += 1 / dist;
            cell2nodeWeight[ni][ei] = 1 / dist;
        end
        for ei=1:length(cell2node[ni])
            cell2nodeWeight[ni][ei] /= sumweight;
        end
    end

    return FVInfo(cell_centers, face_centers, cell2node, cell2nodeWeight);
end

# Interpolate nodal values from neighboring cells.
# Returns an array of nodal values.
function FV_cell_to_node(cell_values, node_values; dofs = 1)
    Nnodes = size(grid_data.allnodes,2);
    if node_values === nothing
        node_values = zeros(Nnodes * dofs);
    end

    for ni=1:Nnodes
        if dofs > 1
            st = (ni-1)*dofs + 1;
            en = ni*dofs;
            node_values[st:en] .= 0;
            ncell = length(fv_info.cell2nodeWeight[ni]);
            for ci=1:ncell
                cind = fv_info.cell2node[ni][ci];
                cst = (cind-1)*dofs + 1;
                cen = cind*dofs;
                node_values[st:en] .+= fv_info.cell2nodeWeight[ni][ci] .* cell_values[cst:cen];
            end

        else
            node_values[ni] = 0;
            ncell = length(fv_info.cell2nodeWeight[ni]);
            for ci=1:ncell
                node_values[ni] += fv_info.cell2nodeWeight[ni][ci] * cell_values[fv_info.cell2node[ni][ci]];
            end
        end
    end

    return node_values;
end

# Find cell averages from nodal values
# Returns an array of cell values.
function FV_node_to_cell(node_values, cell_values = nothing)
    Ncells = size(grid_data.loc2glb,2);
    if cell_values === nothing
        cell_values = zeros(Ncells);
    end

    for ci=1:Ncells
        cell_values[ci] = 0;
        nnode = length(grid_data.loc2glb[:,ci]);
        for ni=1:nnode
            cell_values[ci] += node_values[grid_data.loc2glb[ni, ci]];
        end
        cell_values[ci] /= nnode;
    end

    return cell_values;
end


text:
using Polyorder
using YAML
using CxxWrap
using Test

@testset "Config" begin
    config_file = "./DIS.yml"
    config = Polyorder.parse_config(config_file)

    @test config["Version"]["version"] == 11.1
    @test config["Model"]["n_block"] == [4, 1]
    @test config["IO"]["is_display"] == true
end

@testset "Model" begin
    config_file = "./DIS.yml"
    model, config = Polyorder.init_model(config_file)

    model_alias = model
    @test model_alias === model

    @test Polyorder.Hi(model) ≈ -0.7641464471205148 atol=1e-6
    @test Polyorder.run!(model, config) ≈ 3.6134535528794847 atol=1e-6

    F, mu = Polyorder.free_energy!(model, config)
    @test F ≈ 3.6134535528794847 atol=1e-6
    @test mu ≈ 0.6744577118259163 atol=1e-6
end


text:
### A Pluto.jl notebook ###
# v0.14.3

using Markdown
using InteractiveUtils

# ╔═╡ 7dfe963f-3601-4515-a636-4930816a4e6c
using Pkg, DrWatson

# ╔═╡ 73d28c7f-8426-4c5a-b166-23dcac5c16ce
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers
	using Test
end

# ╔═╡ ac233fe5-a00f-412f-ad4a-a1d39ee20517
md" ### f-01-01.jl"

# ╔═╡ d82838e5-1591-41d9-b5e0-0c49f48d7814
begin
	a = [10. 5. 6.; 5. 20. 4.; 6. 4. 30.]
	x0 = [1., 1., 1.]

	tol = 1.0e-5
	limit = 100
end;

# ╔═╡ 570d9517-14d8-4d40-8e29-039269a5e37d
(iters, eigval, eigvec) = nmex(a, x0; tol=tol, limit=limit, show=8)

# ╔═╡ 1614c806-1bdc-40bb-ae8c-e02d29d128fe
jeigval, jeigvec = eigen(a)

# ╔═╡ 90c91c6a-1b9b-4d34-858a-7535b1090fee
@test round(maximum(jeigval), digits=3) == round(eigval, digits=3)

# ╔═╡ 47256683-e6b6-4fae-916c-f80ba57a232f
md" ### End of f-01-01.jl"

# ╔═╡ Cell order:
# ╟─ac233fe5-a00f-412f-ad4a-a1d39ee20517
# ╠═7dfe963f-3601-4515-a636-4930816a4e6c
# ╠═73d28c7f-8426-4c5a-b166-23dcac5c16ce
# ╠═d82838e5-1591-41d9-b5e0-0c49f48d7814
# ╠═570d9517-14d8-4d40-8e29-039269a5e37d
# ╠═1614c806-1bdc-40bb-ae8c-e02d29d128fe
# ╠═90c91c6a-1b9b-4d34-858a-7535b1090fee
# ╟─47256683-e6b6-4fae-916c-f80ba57a232f


text:
# The @compat macro is used to implement compatibility rules that require
# syntax rewriting rather than simply new function/constant/module definitions.

export @compat

function _compat(ex::Expr)
    if ex.head === :quote && isa(ex.args[1], Symbol)
        # Passthrough
        return ex
    end
    return Expr(ex.head, map(_compat, ex.args)...)
end

_compat(ex) = ex

macro compat(ex)
    esc(_compat(ex))
end


text:
function check_hsg()
    assert_clean_working_directory()
    run_hsg()
    assert_clean_working_directory()
    return nothing
end

function assert_clean_working_directory()
    if !isempty(strip(read(`git status --short`, String)))
        msg = "The working directory is dirty"
        @error msg

        println("Output of `git status`:")
        println(strip(read(`git status`, String)))

        run(`git add -A`)

        println("Output of `git diff HEAD`:")
        println(strip(read(`git diff HEAD`, String)))

        run(`git reset`)

        throw(ErrorException(msg))
    else
        @info "The working directory is clean"
        return nothing
    end
end

function run_hsg()
    env2 = copy(ENV)
    delete!(env2, "JULIA_DEPOT_PATH")
    delete!(env2, "JULIA_LOAD_PATH")
    delete!(env2, "JULIA_PROJECT")
    env2["JULIA_DEPOT_PATH"] = mktempdir(; cleanup = true)
    julia_binary = Base.julia_cmd().exec[1]
    hsg_directory = joinpath("ext", "HistoricaStdlibGenerator")
    hsg_generate_file = joinpath(hsg_directory, "generate_historical_stdlibs.jl")

    cmd_1 = `$(julia_binary) --project=$(hsg_directory) -e 'import Pkg; Pkg.instantiate()'`
    cmd_2 = `$(julia_binary) --project=$(hsg_directory) $(hsg_generate_file)`

    run(setenv(cmd_1, env2))
    run(setenv(cmd_2, env2))

    return nothing
end

check_hsg()


text:
using Base: Indices
using OffsetArrays: OffsetArray

abstract type AbstractDataSeries{T <: AbstractData, N, AT <: AbstractArray{T,N}} <: AbstractArray{T,N} end

function initialize!(ds::AbstractDataSeries{T,1}, q₀::T) where {T}
    for j in eachindex(ds)
        ds[j] = zero(q₀)
    end
    return ds
end

function initialize!(ds::AbstractDataSeries{T,1}, q₀::Vector{T}) where {T}
    @assert nsamples(ds) == length(q₀) == 1
    initialize!(ds, q₀[begin])
end

function initialize!(ds::AbstractDataSeries{T,2}, q₀::Vector{T}) where {T}
    @assert nsamples(ds) == length(q₀)
    for k in axes(ds,2)
        for j in axes(ds,1)
            ds[j,k] = zero(q₀[k])
        end
    end
    return ds
end

function DataSeriesConstructor(::Type{dsType}, ::Type{T}, nt::Int, ni::Int=1) where {dsType <: AbstractDataSeries, T <: Number}
    ds = dsType{T, ni == 1 ? 1 : 2}(nt, ni)
    initialize!(ds, zero(T))
end

function DataSeriesConstructor(::Type{dsType}, ::Type{T}, nd::Int, nt::Int, ni::Int) where {dsType <: AbstractDataSeries, T <: Number}
    ds = dsType{Vector{T}, ni == 1 ? 1 : 2}(nt, ni)
    initialize!(ds, [zeros(T, nd) for _ in 1:ni])
end

function DataSeriesConstructor(::Type{dsType}, q₀::T, nt::Int, ni::Int=1) where {dsType <: AbstractDataSeries, T <: AbstractData}
    @assert ni == 1
    ds = dsType{T,1}(nt, ni)
    initialize!(ds, q₀)
end

function DataSeriesConstructor(::Type{dsType}, q₀::Vector{T}, nt::Int, ni::Int=1) where {dsType <: AbstractDataSeries, T <: Number}
    @assert ni == 1 || ni == length(q₀)
    ds = ( ni == 1 ? dsType{Vector{T},1}(nt, ni) : dsType{T,2}(nt, ni) )
    initialize!(ds, q₀)
end

function DataSeriesConstructor(::Type{dsType}, q₀::Vector{T}, nt::Int, ni::Int=length(q₀)) where {dsType <: AbstractDataSeries, T <: AbstractArray{<:Number}}
    @assert ni == length(q₀)
    ds = dsType{T, ni == 1 ? 1 : 2}(nt, ni)
    initialize!(ds, q₀)
end

function DataSeriesConstructor(::Type{dsType}, d::AbstractArray{T,1}) where {dsType <: AbstractDataSeries, T}
    nt = size(d,1)-1
    ni = 1
    ds = dsType{T,1}(nt, ni)
    copy!(parent(ds), d)
    return ds
end

function DataSeriesConstructor(::Type{dsType}, d::AbstractArray{T,2}) where {dsType <: AbstractDataSeries, T}
    nt = size(d,1)-1
    ni = size(d,2)
    ds = dsType{T,2}(nt, ni)
    copy!(parent(ds), d)
    return ds
end

Base.parent(ds::DS) where {DS <: AbstractDataSeries} = error("parent() not implemented for $(DS)")

Base.eltype(ds::AbstractDataSeries{T,N}) where {T,N} = T
Base.ndims(ds::AbstractDataSeries{T,N}) where {T,N} = N

Base.:(==)(ds1::AbstractDataSeries, ds2::AbstractDataSeries) = (parent(ds1) == parent(ds2))

function Base.show(io::IO, ds::DS) where {T, N, DS <: AbstractDataSeries{T,N}}
    print(io, "$(DS) with data type ", T, " and ", N, " dimensions:\n")
    print(io, parent(ds))
end

ntime(ds::AbstractDataSeries) = lastindex(ds.d, 1) - 1
nsamples(ds::AbstractDataSeries) = size(ds.d, 2)

Base.size(ds::AbstractDataSeries) = size(ds.d)
Base.size(ds::AbstractDataSeries, d) = size(ds.d, d)

function Base.similar(ds::DS) where {T, DS <: AbstractDataSeries{T,1}}
    newds = DS(ntime(ds), nsamples(ds))
    for j in eachindex(ds)
        newds[j] = zero(ds[j])
    end
    return newds
end

function Base.similar(ds::DS) where {T, DS <: AbstractDataSeries{T,2}}
    newds = DS(ntime(ds), nsamples(ds))
    for k in axes(ds,2)
        for j in axes(ds,1)
            newds[j,k] = zero(ds[j,k])
        end
    end
    return newds
end

Base.eachindex(::IndexCartesian, ds::AbstractDataSeries) = CartesianIndices(axes(ds))
# Base.eachindex(::IndexLinear, ds::AbstractDataSeries{T,1}) where {T} = LinearIndices(axes(ds,1))

Base.firstindex(ds::AbstractDataSeries{T,1}) where {T} = 0
Base.firstindex(ds::AbstractDataSeries{T,N}) where {T,N} = firstindex(ds.d)

Base.firstindex(ds::AbstractDataSeries{T,1}, d) where {T} = d ≥ 1 && d ≤ 1 ? 0 : 1
Base.firstindex(ds::AbstractDataSeries{T,2}, d) where {T} = d ≥ 1 && d ≤ 2 ? (0, 1)[d] : 1

Base.lastindex(ds::AbstractDataSeries{T,1}) where {T} = ntime(ds)
Base.lastindex(ds::AbstractDataSeries{T,N}) where {T,N} = lastindex(ds.d)

Base.lastindex(ds::AbstractDataSeries{T,1}, d) where {T} = d ≥ 1 && d ≤ 1 ? ntime(ds) : 1
Base.lastindex(ds::AbstractDataSeries{T,2}, d) where {T} = d ≥ 1 && d ≤ 2 ? (ntime(ds), nsamples(ds))[d] : 1

@inline Base.axes(ds::AbstractDataSeries{T,1}) where {T} = (0:ntime(ds),)
@inline Base.axes(ds::AbstractDataSeries{T,2}) where {T} = (0:ntime(ds), 1:nsamples(ds))
@inline Base.axes(ds::AbstractDataSeries{T,N}, d) where {T,N} = d ≥ 1 && d ≤ N ? axes(ds)[d] : (1:1)

Base.strides(ds::AbstractDataSeries) = strides(ds.d)

function get_data!(ds::AbstractDataSeries{T,1}, x::T, n::Int, k::Int=1) where {T}
    @assert n ≥ 0 && n ≤ ntime(ds)
    @assert k == 1
    x .= ds.d[n+1]
end

function get_data!(ds::AbstractDataSeries{T,2}, x::T, n::Int, k::Int) where {T}
    @assert n ≥ 0 && n ≤ ntime(ds)
    @assert k ≥ 1 && k ≤ nsamples(ds)
    x .= ds.d[n+1,k]
end

function get_data!(ds::AbstractDataSeries{T,2}, x::Vector{T}, n::Int) where {T}
    @assert length(x) == nsamples(ds)
    @assert n ≥ 0 && n ≤ ntime(ds)
    @inbounds for k in eachindex(x)
        x[k] = ds.d[n+1,k]
    end
end

function set_data!(ds::AbstractDataSeries{T,1}, x::T, n::Int, k::Int=1) where {T}
    @assert n ≥ 0 && n ≤ ntime(ds)
    @assert k == 1
    if T <: Number
        ds.d[n+1] = x
    elseif T <: AbstractArray
        ds.d[n+1] .= x
    end
end

function set_data!(ds::AbstractDataSeries{T,1}, x::Vector{T}, n::Int, k::Int=1) where {T}
    @assert nsamples(ds) == length(x) == 1
    set_data!(ds, x[begin], n, k)
end

function set_data!(ds::AbstractDataSeries{T,2}, x::T, n::Int, k::Int) where {T}
    @assert n ≥ 0 && n ≤ ntime(ds)
    if T <: Number
        ds.d[n+1,k] = x
    elseif T <: AbstractArray
        ds.d[n+1,k] .= x
    end
end

function set_data!(ds::AbstractDataSeries{T,2}, x::Vector{T}, n::Int) where {T}
    @assert length(x) == nsamples(ds)
    @assert n ≥ 0 && n ≤ ntime(ds)
    @inbounds for k in axes(ds.d, 2)
        if T isa Number
            ds.d[n+1,k] = x[k]
        elseif T isa AbstractArray
            ds.d[n+1,k] .= x[k]
        end
    end
end

function reset!(ds::AbstractDataSeries{T,1}) where {T}
    @inbounds ds[begin] = ds[end]
end

function reset!(ds::AbstractDataSeries{T,2}) where {T}
    @inbounds for k in axes(ds,2)
        ds[begin,k] = ds[end,k]
    end
end

@inline function Base.getindex(ds::AbstractDataSeries{T,1}, I::Colon) where {T}
    OffsetArray([ds[i] for i in axes(ds,1)], axes(ds,1))
end

@inline function Base.getindex(ds::AbstractDataSeries{T,1}, j::Union{UnitRange,Int}) where {T}
    getindex(ds.d, j.+1)
end

@inline function Base.getindex(ds::AbstractDataSeries{T,1}, i::Union{Int,CartesianIndex}, j::Union{Int,CartesianIndex}) where {T}
    ds[j][i]
end

@inline function Base.getindex(ds::AbstractDataSeries{T,1}, i::Union{Int,CartesianIndex}, J::AbstractRange{Int}) where {T}
    [ds[j][i] for j in J]
end

@inline function Base.getindex(ds::AbstractDataSeries{T,1}, i::Union{Int,CartesianIndex}, J::Colon) where {T}
    OffsetArray([ds[j][i] for j in axes(ds,1)], axes(ds,1))
end

@inline Base.getindex(ds::AbstractDataSeries{T,1}, I, J::Colon) where {T} = getindex(ds, I, axes(ds,1))

@inline function Base.getindex(ds::AbstractDataSeries{T,2}, j::Union{UnitRange,Int}, k::Union{UnitRange,Int}) where {T}
    getindex(ds.d, j.+1, k)
end

@inline function Base.getindex(ds::AbstractDataSeries{T,2}, i::Union{Int,CartesianIndex}, j::Union{Int,CartesianIndex}, k::Union{Int,CartesianIndex}) where {T}
    ds[j,k][i]
end

@inline function Base.getindex(ds::AbstractDataSeries{T,2}, i::Union{Int,CartesianIndex}, J::AbstractRange{Int}, K::AbstractRange{Int}) where {T}
    [ds[j,k][i] for j in J, k in K]
end

@inline function Base.getindex(ds::AbstractDataSeries{T,2}, i::Union{Int,CartesianIndex}, J::Colon, K::AbstractRange{Int}) where {T}
    OffsetArray([ds[j,k][i] for j in axes(ds,1), k in K], axes(ds,1), K)
end

@inline function Base.getindex(ds::AbstractDataSeries{T,2}, i::Union{Int,CartesianIndex}, J::AbstractRange{Int}, k::Union{Int,CartesianIndex}) where {T}
    [ds[j,k][i] for j in J]
end

@inline function Base.getindex(ds::AbstractDataSeries{T,2}, i::Union{Int,CartesianIndex}, J::Colon, k::Union{Int,CartesianIndex}) where {T}
    OffsetArray([ds[j,k][i] for j in axes(ds,1)], axes(ds,1))
end

@inline function Base.getindex(ds::AbstractDataSeries{T,2}, i::Union{Int,CartesianIndex}, j::Union{Int,CartesianIndex}, K::AbstractRange{Int}) where {T}
    [ds[j,k][i] for k in K]
end

@inline Base.getindex(ds::AbstractDataSeries{T,2}, J::Colon, K::Colon) where {T} = getindex(ds, axes(ds,1), axes(ds,2))
@inline Base.getindex(ds::AbstractDataSeries{T,2}, J, K::Colon) where {T} = getindex(ds, J, axes(ds,2))
@inline Base.getindex(ds::AbstractDataSeries{T,2}, J::Colon, K) where {T} = getindex(ds, axes(ds,1), K)

@inline Base.getindex(ds::AbstractDataSeries{T,2}, I, J::Colon, K::Colon) where {T} = getindex(ds, I, axes(ds,1), axes(ds,2))
@inline Base.getindex(ds::AbstractDataSeries{T,2}, I, J::Colon, K) where {T} = getindex(ds, I, axes(ds,1), K)
@inline Base.getindex(ds::AbstractDataSeries{T,2}, I, J, K::Colon) where {T} = getindex(ds, I, J, axes(ds,2))

@inline function Base.setindex!(ds::AbstractDataSeries{T,1}, x, j::Union{UnitRange,Int}) where {T}
    setindex!(ds.d, x, j.+1)
end

@inline function Base.setindex!(ds::AbstractDataSeries{T,2}, x, j::Union{UnitRange,Int}, k::Union{UnitRange,Int}) where {T}
    setindex!(ds.d, x, j.+1, k)
end

#
# @inline function Base.getindex(ds::AbstractDataSeries{T,2}, I::Union{Int,CartesianIndex}, J::Union{Int,CartesianIndex,AbstractRange{Int}}, K::Union{Int,CartesianIndex,AbstractRange{Int}}) where {T}
#     x = zeros(eltype(ds.d[begin]), (1 for i in 1:length(I))..., length(J))
#     for j in J
#         x[Tuple(I)..., j] = ds[j][I]
#     end
# end
#
# similar(storagetype, shape)
#

function fromarray(::Type{dsType}, d::AbstractArray{T,1}, ni::Int=1) where {dsType <: AbstractDataSeries, T <: Number}
    @assert ni == 1
    return dsType(d)
end

function fromarray(::Type{dsType}, d::AbstractArray{T,2}, ni::Int=1) where {dsType <: AbstractDataSeries, T <: Number}
    @assert ni == 1 || ni == size(d,2)

    if ni == 1
        ds = [Vector{T}(d[:,n]) for n in axes(d,2)]
    else
        ds = Array(d)
    end

    return dsType(ds)
end

function fromarray(::Type{dsType}, d::AbstractArray{T,N}, ni::Int=1) where {dsType <: AbstractDataSeries, T <: Number, N}
    @assert ni == 1 || ni == size(d)[end]
    @assert N ≥ 3

    if ni == 1
        AT = Array{T,ndims(d)-1}
        ds = [AT(d[axes(d)[1:end-1]..., n]) for n in axes(d)[end]]
    else
        AT = Array{T,ndims(d)-2}
        ds = [AT(d[axes(d)[1:end-2]..., n, k]) for n in axes(d)[end-1], k in axes(d)[end]]
    end

    return dsType(ds)
end

struct DataSeries{T <: AbstractData, N} <: AbstractDataSeries{T, N, Array{T,N}}
    d::Array{T,N}

    function DataSeries{T,N}(nt, ni) where {T <: AbstractData, N}
        @assert nt ≥ 0
        @assert ni > 0

        @assert N ∈ (1,2)

        if N == 1
            d = Array{T,N}(undef, nt+1)
        elseif N == 2
            d = Array{T,N}(undef, nt+1, ni)
        end

        new(d)
    end
end

Base.parent(ds::DataSeries) = ds.d

DataSeries(::Type{T}, nt::Int, ni::Int=1) where {T} = DataSeriesConstructor(DataSeries, T, nt, ni)
DataSeries(::Type{T}, nd::Int, nt::Int, ni::Int) where {T} = DataSeriesConstructor(DataSeries, T, nd, nt, ni)
DataSeries(q₀::T, nt::Int, ni::Int=1) where {T <: AbstractData} = DataSeriesConstructor(DataSeries, q₀, nt, ni)
DataSeries(q₀::Vector{T}, nt::Int, ni::Int=1) where {T <: Number} = DataSeriesConstructor(DataSeries, q₀, nt, ni)
DataSeries(q₀::Vector{T}, nt::Int, ni::Int=length(q₀)) where {T <: AbstractArray{<:Number}} = DataSeriesConstructor(DataSeries, q₀, nt, ni)
DataSeries(d::AbstractArray) = DataSeriesConstructor(DataSeries, d)


text:
# Addresses https://github.com/JuliaAI/MLFlowClient.jl/issues/17
include("test_base.jl")

@testset "Issue17" begin
    @ensuremlf

    mlf_experiment = getorcreateexperiment(mlf, "issue17")
    mlf_run = createrun(
        mlf,
        mlf_experiment,
        tags=[
            Dict(
                "key" => "mlflow.runName",
                "value" => "run_name"
            )
        ]
    )
    @test isa(mlf_run, MLFlowRun)
end

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
using HTTP
using JSON

for (_id, convo) in JSON.parse(String(HTTP.get("https://raw.githubusercontent.com/BYU-PCCL/chitchat-dataset/master/chitchat_dataset/dataset.json").body))
	for message in convo["messages"]
		for utterance in message
			println(utterance["text"])
		end
	end
end

max_stars_repo_name: NancyFulda/chitchat-dataset

max_issues_repo_name: NancyFulda/chitchat-dataset

max_forks_repo_name: NancyFulda/chitchat-dataset


text:
# Autogenerated wrapper script for libcleri_jll for x86_64-w64-mingw32
export libcleri

using PCRE2_jll
JLLWrappers.@generate_wrapper_header("libcleri")
JLLWrappers.@declare_library_product(libcleri, "libcleri.dll")
function __init__()
    JLLWrappers.@generate_init_header(PCRE2_jll)
    JLLWrappers.@init_library_product(
        libcleri,
        "bin\\libcleri.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
import Distributions
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")

@halfwidth_plot_defaults()

plotlyjs()

ic = "point_mass"
function p()
    ks_data = CSV.read((@__DIR__)*"/../../data/errors/ks_"*ic*".csv",DataFrame)
    tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/ks_"*ic*".csv",DataFrame)
    ks_data.DG_limit2 = tmp.DG_limit

    ks_data_lwr = CSV.read((@__DIR__)*"/../../data/errors/ks_lwr_"*ic*".csv",DataFrame)
    tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/ks_lwr_"*ic*".csv",DataFrame)
    ks_data_lwr.DG_limit2 = tmp.DG_limit

    ks_data_upr = CSV.read((@__DIR__)*"/../../data/errors/ks_upr_"*ic*".csv",DataFrame)
    tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/ks_upr_"*ic*".csv",DataFrame)
    ks_data_upr.DG_limit2 = tmp.DG_limit

    plot()
    linestyles_vec = [:solid,:dash,:dashdot, :dot, :dash]
    markerstyles_vec = [:cross,:diamond,:circle,:dot,:diamond]
    colours = [1;3;2;4;5]
    _names = ["DG", "DG lim", "Unif","QBD-RAP", "DG lin lim"]

    for (c,col) in enumerate(names(ks_data))
        if c<5
            plot!(log10.(1:2:21),log10.(ks_data[:,col]),
                ribbon=(log10.(ks_data[:,col])-log10.(ks_data_lwr[:,col]),log10.(ks_data_upr[:,col])-log10.(ks_data[:,col])),
                label=_names[c],linestyle=linestyles_vec[c],
                # marker=markerstyles_vec[c],
                linewidth=2,color=colours[c])
        else
            plot!(log10.((1:2:21).+1),log10.(ks_data[:,col]),
                ribbon=(log10.(ks_data[:,col])-log10.(ks_data_lwr[:,col]),log10.(ks_data_upr[:,col])-log10.(ks_data[:,col])),
                label=_names[c],linestyle=linestyles_vec[c],
                # marker=markerstyles_vec[c],
                linewidth=2,color=colours[c])
        end
    end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error")
# plot!(legend=(0.8,0.8))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/ks_error_formatted.pdf")
end
p()

function p()
    l1_cdf_data = CSV.read((@__DIR__)*"/../../data/errors/l1_"*ic*".csv",DataFrame)
    tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/l1_"*ic*".csv",DataFrame)
    l1_cdf_data.DG_limit2 = tmp.DG_limit

    l1_cdf_data_lwr = CSV.read((@__DIR__)*"/../../data/errors/l1_lwr_"*ic*".csv",DataFrame)
    tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/l1_lwr_"*ic*".csv",DataFrame)
    l1_cdf_data_lwr.DG_limit2 = tmp.DG_limit

    l1_cdf_data_upr = CSV.read((@__DIR__)*"/../../data/errors/l1_upr_"*ic*".csv",DataFrame)
    tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/l1_upr_"*ic*".csv",DataFrame)
    l1_cdf_data_upr.DG_limit2 = tmp.DG_limit
    _names = ["DG", "DG lim", "Unif","QBDRAP", "DG lin lim"]

    plot()
    linestyles_vec = [:solid,:dash,:dashdot, :dot, :dash]
    markerstyles_vec = [:cross,:diamond,:circle,:cross, :diamond]
    colours = [1;3;2;4;5]

    for (c,col) in enumerate(names(l1_cdf_data))
        if c < 5
            plot!(log10.(1:2:21),log10.(l1_cdf_data[:,col]),
                ribbon=(log10.(l1_cdf_data[:,col])-log10.(l1_cdf_data_lwr[:,col]),log10.(l1_cdf_data_upr[:,col])-log10.(l1_cdf_data[:,col])),
                label=_names[c],linestyle=linestyles_vec[c],
                # marker=markerstyles_vec[c],
                linewidth=2,color=colours[c])
        else
            plot!(log10.((1:2:21).+1),log10.(l1_cdf_data[:,col]),
                ribbon=(log10.(l1_cdf_data[:,col])-log10.(l1_cdf_data_lwr[:,col]),log10.(l1_cdf_data_upr[:,col])-log10.(l1_cdf_data[:,col])),
                label=_names[c],linestyle=linestyles_vec[c],
                # marker=markerstyles_vec[c],
                linewidth=2,color=colours[c])
        end
    end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - CDF")
# plot!(legend=(0.95,0.95))
# plot!(legend=(0.215,0.535))
# plot!(ylims=ylims(plot!()).+(-0.13,0.0),xlims=xlims(plot!()).+(-0.03,0.0))
@add_lines!(l1_cdf_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/l1_cdf_error_formatted.pdf")
end
p()

function p()
colours = [1;3;2;4;5]

@fullwidth_plot_defaults()
sim_data = CSV.read((@__DIR__)*"/../../data/point_mass/sims.csv",DataFrame)
data = CSV.read((@__DIR__)*"/../../data/point_mass/order_21_model_dg1.csv",DataFrame)
T_idx = only(findall(data.t.==1.2))
sim_cdf = CSV.read((@__DIR__)*"/../../data/point_mass/sims_cdf_evaluated.csv",DataFrame).phase_1[1:T_idx]
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
plot(data.t[1:T_idx],data.phase_1[1:T_idx],label="DG - 21",color=colours[1],linestyle=linestyles_vec[1],linewidth=2)
data = CSV.read((@__DIR__)*"/../../data/point_mass/order_21_model_qbdrap4.csv",DataFrame)
plot!(data.t[1:T_idx],data.phase_1[1:T_idx],label="QBDRAP - 21",color=colours[4],linestyle=linestyles_vec[4],linewidth=2)
data = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/point_mass/order_21_model_dg1.csv",DataFrame)
plot!(data.t[1:T_idx],data.phase_1[1:T_idx],label="DG lin lim - 22",color=colours[5],linestyle=linestyles_vec[5],linewidth=2)
plot!(data.t[1:T_idx],sim_cdf,label="Sim",color=6,linestyle=:dashdot,linewidth=2)
plot!(xlabel="t")
plot!(ylabel="Probability")
plot!(title="First hitting time CDF")
plot!(legend=(0.2,0.9))
savefig((@__DIR__)*"/cdf_order21DG_and_sims.pdf")

end
p()


text:
include("am.jl")

using BenchmarkTools

# Create an instance of Household
am = Household(r = 0.03, β = 0.96, a_size = 500)
a_size = am.a_size
V, g_a, g_c, g_h = initialize(am)

println("===============")
println("TIME:")
println("===============")
println("VFI - Discrete with no refinements. Grid size = $a_size")
@btime compute_fixed_point(v -> BruteForce_Std(v, am),
                        V,
                        max_iter=1000,
                        err_tol=1e-4,
                        verbose=false)
println("===============")
println("VFI - Discrete using Monotonicity. Grid size = $a_size")
@btime compute_fixed_point(v -> BruteForce_Mtcy(v, am),
                        V,
                        max_iter=1000,
                        err_tol=1e-4,
                        verbose=false)
println("===============")
println("VFI - Continous using linear interpolation. Grid size = $a_size")
@btime compute_fixed_point(v -> BellmanOperator(v, am),
                        V,
                        max_iter=1000,
                        err_tol=1e-4,
                        verbose=false)
println("===============")
println("PFI - Continous using linear interpolation. Grid size = $a_size")
@btime compute_fixed_point(c -> ColemanOperator(c, am),
                        g_c,
                        max_iter=1000,
                        err_tol=1e-4,
                        verbose=false)
println("===============")
println("Endogeous Grid Method. Grid size = $a_size")
@btime compute_fixed_point(c -> EGMOperator_Labor(c, am),
                        g_c,
                        max_iter=1000,
                        err_tol=1e-4,
                        verbose=false)


text:
using ODBC

conn = ODBC.connect("Chinook",usr="malcolm",pwd="mypasswd")
res = ODBC.query("select count(*) from Customers", conn)
println("Number of customers: $(res[:K][1])")

sql = "select a.LastName, a.FirstName,";
sql *= " count(b.InvoiceId) as Invs, sum(b.Total) as Amt";
sql *= " from Customer a";
sql *= " join Invoice b on a.CustomerId = b.CustomerId";
sql *= " group by a.LastName having Amt >= 45.00";
sql *= " order by Amt desc;";

res = ODBC.query(sql);

for i in 1:size(res)[1]
  LastName  = res[:LastName][i]
  FirstName = res[:FirstName][i]
  Invs      = res[:Invs][i]
  Amt       = res[:Amt][i]
  @printf "%10s %10s %4d %10.2f\n" LastName FirstName Invs Amt
end

sql = "select a.LastName, a.FirstName, d.Name as TrackName";
sql *= " from Customer a";
sql *= " join Invoice b on a.CustomerId = b.CustomerId";
sql *= " join InvoiceLine c on b.InvoiceId = c.InvoiceId";
sql *= " join Track d on c.TrackId = d.TrackId";
sql *= " where a.LastName = 'Cunningham' limit 5;";

res = ODBC.query(sql);

for i in 1:size(res)[1]
  LastName  = res[:LastName][i]
  FirstName = res[:FirstName][i]
  TrackName = res[:TrackName][i]
  @printf "%15s %15s %15s\n" LastName FirstName TrackName
end


text:
using ValidatedNumerics, ValidatedNumerics.RootFinding
using Test

@testset "Bisection tests" begin
    @testset "Single variable" begin
        f(x) = sin(x) + cos(x/2) + x/5

        roots = bisection(f, -∞..∞, tolerance=1e-3)
        roots2 = [root.interval for root in roots]

        @test roots2 ==  [  Interval(-0.8408203124999999, -0.8398437499999999),
                            Interval(3.6425781249999996, 3.6435546874999996),
                            Interval(6.062499999999999, 6.063476562499999)
                         ]

    end

    @testset "Two variables" begin
        @intervalbox g(x, y) = (x^2 + y^2 - 1, y - x)

        X = (-∞..∞) × (-∞..∞)
        roots = bisection(g, X, tolerance=1e-3)
        roots2 = [root.interval for root in roots]

        @test roots2 == [IntervalBox(Interval(-0.7080078124999999, -0.7070312499999999), Interval(-0.7080078124999999, -0.7070312499999999)),
 IntervalBox(Interval(-0.7080078124999999, -0.7070312499999999), Interval(-0.7070312499999999, -0.7060546874999999)),
 IntervalBox(Interval(-0.7070312499999999, -0.7060546874999999), Interval(-0.7080078124999999, -0.7070312499999999)),
 IntervalBox(Interval(0.7060546874999999, 0.7070312499999999), Interval(0.7070312499999999, 0.7080078124999999)),
 IntervalBox(Interval(0.7070312499999999, 0.7080078124999999), Interval(0.7060546874999999, 0.7070312499999999)),
 IntervalBox(Interval(0.7070312499999999, 0.7080078124999999), Interval(0.7070312499999999, 0.7080078124999999))]

    end

end


text:
# Genetic Algorithms
# ==================
#         objfun: Objective fitness function
#              N: Search space dimensionality
# initPopulation: Generation function which produce individual population entities.
# populationSize: Size of the population
#  crossoverRate: The fraction of the population at the next generation, not including elite children,
#                 that is created by the crossover function.
#   mutationRate: Probability of chromosome to be mutated
#              ɛ: Positive integer specifies how many individuals in the current generation
#                 are guaranteed to survive to the next generation.
#                 Floating number specifies fraction of population.
#
#
function sGA(objfun::Function,
             N::Int,
             initPopulation::Function = randperm,
             populationSize::Int = 50,
             crossoverRate::Float64 = 0.8,
             mutationRate::Float64 = 0.1,
             ɛ::Real = 0,
             selection::Function = ((x,n)->1:n),
             crossover::Function = ((x,y)->(y,x)),
             mutation::Function = (x->x),
             iterations::Int = 10*N,
             tol = 0.0,
             tolIter = 10,
             vstep::Int = iterations,
             rebal = false,
             inform = true,
             robust = false,
             signature::String = "TSP$N")

    store = Dict{Symbol,Any}()

    # Setup parameters
    elite = isa(ɛ, Int) ? ɛ : round(Int, ɛ * populationSize)
    fitFunc = inverseFunc(objfun)

    # Initialize population
    fitness = zeros(Float64,populationSize)
    trueFit = zeros(Float64,populationSize)
    population = Array(Vector{Int}, populationSize)
    offspring = similar(population)

    # Generate population
    for i in 1:populationSize
      population[i] = initPopulation(N)
      fitness[i] = fitFunc(population[i])
    end
    fitidx = sortperm(fitness, rev = true)

    #store descriptive information about run
    keep(inform, :selection, string(selection), store)
    keep(inform, :crossover, string(crossover), store)
    keep(inform, :mutation, string(mutation), store)
    keep(inform, :cxRate, crossoverRate, store)
    keep(inform, :mxRate, mutationRate, store)
    keep(inform, :elite, elite, store)
    keep(inform, :popsize, populationSize, store)
    keep(inform, :id, signature, store)
    keep(inform, :ipop, string(initPopulation), store)

    # Generate and evaluate offspring
    itr = 1
    bestFitness = 0.0
    bestIndividual = 0
    fittol = 0.0
    fittolitr = 1
    while true
        # Select offspring
        selected = selection(fitness, populationSize)

        # Perform mating
        offidx = randperm(populationSize)
        for i in 1:2:populationSize
            j = (i == populationSize) ? i-1 : i+1
            if rand() < crossoverRate
                offspring[i], offspring[j] = crossover(population[selected[offidx[i]]], population[selected[offidx[j]]])
            else
                offspring[i], offspring[j] = copy(population[selected[i]]), copy(population[selected[j]])
            end
        end

        # Perform mutation
        for i in 1:populationSize
            if rand() < mutationRate
                mutation(offspring[i])
            end
        end

        # Elitism
        if elite > 0
          subs = sample(1:populationSize,elite,replace=false)
            for i in 1:elite
                offspring[subs[i]] = population[fitidx[i]]
            end
        end

        # New generation
        for i in 1:populationSize
            population[i] = copy(offspring[i])
            fitness[i] = fitFunc(offspring[i])
            trueFit[i] = objfun(offspring[i])
            rebal && (population[i] = rebalTo1(population[i]))
        end
        fitidx = sortperm(fitness, rev = true)
        ss = summarystats(trueFit)

        # Store generational Data
        keep(robust, :min, ss.min, store)
        keep(robust, :max, ss.max, store)
        keep(robust, :mean, ss.mean, store)
        keep(robust, :median, ss.median, store)

        bestIndividual = fitidx[1]

        curGenFitness = ss.min
        fittol = abs(bestFitness - curGenFitness)
        bestFitness = curGenFitness

        # Verbose step
        (mod(itr,vstep) == 0) && println("BEST:", bestFitness, " : ", "G: ", itr)

        # Terminate:
        #  if fitness tolerance is met for specified number of steps
        if fittol < tol
            if fittolitr > tolIter
                break
            else
                fittolitr += 1
            end
        else
            fittolitr = 1
        end
        # if number of iterations more then specified
        if itr >= iterations
            break
        end
        itr += 1
    end
    keep(inform, :best, bestFitness, store)
    keep(inform, :pop, population, store)

    # return population[bestIndividual], bestFitness, itr, fittol, store
    return store
end


text:
@testset "Cora" begin
    data  = Cora()
    @test length(data) == 1
    g = data[1]
    @test g isa MLDatasets.Graph

    @test g.num_nodes == 2708
    @test g.num_edges == 10556
    @test size(g.node_data.features) == (1433, g.num_nodes)
    @test size(g.node_data.targets) == (g.num_nodes,)
    @test size(g.node_data.train_indices) == (140,)
    @test size(g.node_data.val_indices) == (500,)
    @test size(g.node_data.test_indices) == (1000,)
    @test g.edge_index isa Tuple{Vector{Int}, Vector{Int}}
    s, t = g.edge_index
    for a in (s, t)
        @test a isa Vector{Int}
        @test length(a) == g.num_edges
        @test minimum(a) == 1
        @test maximum(a) == g.num_nodes
    end
end

@testset "CiteSeer" begin
    data  = CiteSeer()
    @test length(data) == 1
    g = data[1]
    @test g isa MLDatasets.Graph

    @test g.num_nodes == 3327
    @test g.num_edges == 9104
    @test size(g.node_data.features) == (3703, g.num_nodes)
    @test size(g.node_data.targets) == (g.num_nodes,)
    @test size(g.node_data.train_indices) == (120,)
    @test size(g.node_data.val_indices) == (500,)
    @test size(g.node_data.test_indices) == (1015,)
    @test g.edge_index isa Tuple{Vector{Int}, Vector{Int}}
    s, t = g.edge_index
    for a in (s, t)
        @test a isa Vector{Int}
        @test length(a) == g.num_edges
        @test minimum(a) == 1
        @test maximum(a) == g.num_nodes
    end
end

@testset "PubMed" begin
    data  = PubMed()
    @test length(data) == 1
    g = data[1]
    @test g isa MLDatasets.Graph

    @test g.num_nodes == 19717
    @test g.num_edges == 88648
    @test size(g.node_data.features) == (500, g.num_nodes)
    @test size(g.node_data.targets) == (g.num_nodes,)
    @test size(g.node_data.train_indices) == (60,)
    @test size(g.node_data.val_indices) == (500,)
    @test size(g.node_data.test_indices) == (1000,)
    @test g.edge_index isa Tuple{Vector{Int}, Vector{Int}}
    s, t = g.edge_index
    for a in (s, t)
        @test a isa Vector{Int}
        @test length(a) == g.num_edges
        @test minimum(a) == 1
        @test maximum(a) == g.num_nodes
    end
end

@testset "KarateClub" begin
    data  = KarateClub()
    @test length(data) == 1
    g = data[1]
    @test g isa MLDatasets.Graph

    @test g.num_nodes == 34
    @test g.num_edges == 156
    @test size(g.node_data.labels_clubs) == (g.num_nodes,)
    @test sort(unique(g.node_data.labels_clubs)) == [0, 1]
    @test size(g.node_data.labels_comm) == (g.num_nodes,)
    @test sort(unique(g.node_data.labels_comm)) == [0, 1, 2, 3]

    @test g.edge_index isa Tuple{Vector{Int}, Vector{Int}}
    s, t = g.edge_index
    for a in (s, t)
        @test a isa Vector{Int}
        @test length(a) == g.num_edges
        @test minimum(a) == 1
        @test maximum(a) == g.num_nodes
    end
end

@testset "PolBlogs" begin
    data  = PolBlogs()
    @test length(data) == 1
    g = data[1]
    @test g isa MLDatasets.Graph

    @test g.num_nodes == 1490
    @test g.num_edges == 19025
    @test size(g.node_data.labels) == (g.num_nodes,)
    @test sort(unique(g.node_data.labels)) == [0, 1]

    @test g.edge_index isa Tuple{Vector{Int}, Vector{Int}}
    s, t = g.edge_index
    for a in (s, t)
        @test a isa Vector{Int}
        @test length(a) == g.num_edges
        @test minimum(a) >= 1
        @test maximum(a) <= g.num_nodes
    end
end


text:
"""
    apply_filter_extend(u, filter, domain)

Extend function beyond beyond boundaries before filtering.
"""
function apply_filter_extend(u, filter, domain)
    error("Not implemented")
end

function apply_filter_extend(u_int, filter::TopHatFilter, domain::ClosedIntervalDomain)
    h = filter.width
    a, b = domain.left, domain.right

    function ū_ext(x, t)
        xₗ, xᵣ = x - h(x), x + h(x)
        1 / 2h(x) * (
            u_int(min(b, xᵣ), t) - u_int(max(a, xₗ), t) +
            g_a(t) * max(0.0, a - xₗ) +
            g_b(t) * max(0.0, xᵣ - b)
        )
    end
end

"""
    apply_filter_extend(u, filter, domain)

Extend function beyond beyond boundaries before filtering.
"""
function apply_filter_extend(u, filter::ConvolutionalFilter, domain::ClosedIntervalDomain)
    h = filter.width
    G = filter.kernel
    a, b = domain.left, domain.right

    # Exact filtered solution
    function ū_ext(x)
        x₋, x₊ = x - h(x), x + h(x)
        xₗ = max(a, x₋)
        xᵣ = min(b, x₊)

        Gfun = Fun(ξ -> G(ξ - x), x₋..x₊)
        Gfun /= sum(Gfun)

        S = sum(Fun(Gfun, xₗ..xᵣ) * Fun(u, xₗ..xᵣ))
        if x₋ < a
            S += u(a) * sum(Fun(Gfun, x₋..a))
        end
        if b < x₊
            S += u(b) * sum(Fun(Gfun, b..x₊))
        end

        S
    end
end


text:
include("../../MaximinOPF/src/MaximinOPF.jl")
using JuMP
using PowerModels
using SCS
using ProxSDP
PowerModels.silence()

testcase = Dict(
	"file" => "data/case30.m",
 	"name" => "case30K4",
 	"attack_budget" => 4,
 	"inactive_indices" => [],
 	"protected_indices" => []
	)

pm_data = PowerModels.parse_file(testcase["file"])
pm_data["attacker_budget"] = testcase["attack_budget"] ###Adding another key and entry
pm_data["inactive_branches"] = testcase["inactive_indices"] ###Adding another key and entry
pm_data["protected_branches"] = testcase["protected_indices"] ###Adding another key and entry

### "sdp_pm=[ SDPWRMPowerModel, SparseSDPWRMPowerModel ]"

conic_solver=with_optimizer(SCS.Optimizer,verbose=0,max_iters=30000)
sdp_solver=with_optimizer(SCS.Optimizer,verbose=0,max_iters=30000)
sparse_sdp_solver=with_optimizer(SCS.Optimizer,verbose=0,max_iters=30000)
#sdp_solver=with_optimizer(ProxSDP.Optimizer, max_iter=100000, log_verbose=false, initial_beta=0.005, full_eig_decomp=false )
#sparse_sdp_solver=with_optimizer(ProxSDP.Optimizer, max_iter=100000,log_verbose=false,initial_beta=0.005,full_eig_decomp=true)

for pm_form in MaximinOPF.conic_supported_pm
    println("\nFormulating and solving the dual minmax problem for the form ",pm_form)
    dual_minmax_model=MaximinOPF.MaximinOPFModel(pm_data, pm_form; enforce_int=false)
    if pm_form == SDPWRMPowerModel
        JuMP.set_optimizer(dual_minmax_model,sdp_solver)
    elseif pm_form == SparseSDPWRMPowerModel
        JuMP.set_optimizer(dual_minmax_model,sparse_sdp_solver)
    else
        JuMP.set_optimizer(dual_minmax_model,conic_solver)
    end
    JuMP.optimize!(dual_minmax_model)
    println("\tOptimal value using powerform ", pm_form, " is: ",JuMP.objective_value(dual_minmax_model), " with status ",JuMP.termination_status(dual_minmax_model))
end


text:
function check_precondition(::Planetoid, dataset::Symbol)
    @assert dataset in subdatasets(Planetoid) "`dataset` should be one of citeseer, cora, pubmed."
end

function check_precondition(::Entities, dataset::Symbol)
    @assert dataset in subdatasets(Entities) "`dataset` should be one of :aifb, :am, :mutag, :bgs."
end

traindata(::Dataset) = throw(ArgumentError("Undefined dataset."))

function traindata(pla::Planetoid, dataset::Symbol)
    check_precondition(pla, dataset)
    file = @datadep_str "Planetoid/$(dataset).train.jld2"
    @load file train_X train_y
    train_X, train_y
end

# function traindata(cora::Cora)
#     file = datadep"Cora/cora.train.jld2"
#     @load file graph train_X train_y
#     graph, train_X, train_y
# end

function traindata(::PPI)
    file = datadep"PPI/ppi.train.jld2"
    @load file graph X y ids
    graph, X, y, ids
end

# function traindata(ent::Entities, dataset::Symbol)
#     check_precondition(ent, dataset)
#     file = @datadep_str "Entities/$(dataset).train.jld2"
#     @load file graph train_X train_y
#     graph, train_X, train_y
# end

validdata(::Dataset) = throw(ArgumentError("Undefined dataset."))

function validdata(::PPI)
    file = datadep"PPI/ppi.valid.jld2"
    @load file graph X y ids
    graph, X, y, ids
end

testdata(::Dataset) = throw(ArgumentError("Undefined dataset."))

function testdata(pla::Planetoid, dataset::Symbol)
    check_precondition(pla, dataset)
    file = @datadep_str "Planetoid/$(dataset).test.jld2"
    @load file test_X test_y
    test_X, test_y
end

# function testdata(dataset::Symbol)
#     file = datadep"Cora/cora.test.jld2"
#     @load file graph test_X test_y
#     graph, test_X, test_y
# end

function testdata(::PPI)
    file = datadep"PPI/ppi.test.jld2"
    @load file graph X y ids
    graph, X, y, ids
end

# function testdata(ent::Entities, dataset::Symbol)
#     check_precondition(ent, dataset)
#     file = @datadep_str "Entities/$(dataset).test.jld2"
#     @load file graph test_X test_y
#     graph, test_X, test_y
# end

function graphdata(pla::Planetoid, dataset::Symbol)
    check_precondition(pla, dataset)
    file = @datadep_str "Planetoid/$(dataset).graph.jld2"
    @load file sg
    sg
end

function graphdata(cora::Cora)
    file = datadep"Cora/cora.graph.jld2"
    @load file sg
    sg
end

function graphdata(::Reddit)
    file = datadep"Reddit/reddit.graph.jld2"
    @load file sg
    sg
end

function alldata(pla::Planetoid, dataset::Symbol)
    check_precondition(pla, dataset)
    file = @datadep_str "Planetoid/$(dataset).all.jld2"
    @load file all_X all_y
    all_X, all_y
end

function alldata(cora::Cora)
    file = datadep"Cora/cora.all.jld2"
    @load file all_X all_y
    all_X, all_y
end

function alldata(::Reddit)
    file = datadep"Reddit/reddit.all.jld2"
    @load file all_X all_y
    all_X, all_y
end

function rawdata(pla::Planetoid, dataset::Symbol)
    check_precondition(pla, dataset)
    file = @datadep_str "Planetoid/$(dataset).raw.jld2"
    @load file raw
    raw
end

function rawdata(cora::Cora)
    file = datadep"Cora/cora.raw.jld2"
    @load file raw
    raw
end

function rawdata(::Reddit)
    file = datadep"Reddit/reddit.raw.jld2"
    @load file raw
    raw
end

function rawdata(::QM7b)
    file = datadep"QM7b/qm7b.raw.jld2"
    @load file raw
    raw
end

function metadata(pla::Planetoid, dataset::Symbol)
    check_precondition(pla, dataset)
    file = @datadep_str "Planetoid/$(dataset).metadata.jld2"
    @load file meta
    meta
end

function metadata(cora::Cora)
    file = datadep"Cora/cora.metadata.jld2"
    @load file meta
    meta
end

function metadata(::Reddit)
    file = datadep"Reddit/reddit.metadata.jld2"
    @load file meta
    meta
end


text:
export Mv, manynth, dimsnth

# """
# Multivariate distribution: Random array where each variable is ciid given
# values for parameters.

# `Mv(dist, shape)`

# # Arguments
# - `dist` a variable class, i.e. `dist(id, ω)` must b defined, e.g. `Normal(0, 1)`
# - `shape` Dimensions of Multivariate

# # Returns

# # Example
# ```julia
# x = 1 ~ Normal(0, 1)
# function f(id, ω)
#   x(ω) + Uniform(0, 1)(id, ω)
# end
# xs = 2 ~ Mv(f, (3, 3))
# randsample((x, xs))
# ```
# """
# struct Mv{T, SHAPE}
#   dist::T
#   shape::SHAPE
# end

# Mv(dist, N::Integer) = Mv(dist, (N,))

# traitlift(::Type{<:Mv}) = Lift()

# # Base.eltype(Mv{T}) where {T} =
# prim(d::Normal) = StdNormal()
# func(d::Normal, x) = x * d.σ + d.μ
# @inline Var.recurse(mv::Mv{<:Distribution}, id, ω) =
#   map(x -> func(mv.dist, x), resolve(Mv(prim(mv.dist), mv.shape), id, ω))
# f(x::Dims) = map(i->1:i, x)
# g(x::Dims) = Iterators.product(f(x)...)
# @inline Var.recurse(mv::Mv{<:T}, id, ω) where T =
#   map(id_ -> mv.dist(append(id, id_), ω), g(mv.shape))
# Base.rand(rng::AbstractRNG, mv::Mv{<:PrimDist}) =
#   rand(rng, mv.dist, mv.shape)

struct Mv{IDXS, FS} <:  AbstractVariable
  idxs::IDXS
  f::FS
end

Var.recurse(mv::Mv, ω) =  map(i -> mv.f(i, ω), mv.idxs)

# FIXME: Maybe the ids should come first, not sure
@inline manynth(f, ids) = Mv(ids, f)
@inline dimsnth(f, shape::Dims) = Mv(CartesianIndices(shape), f)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# true composite models:
abstract type DeterministicComposite <: Deterministic end
abstract type ProbabilisticComposite <: Probabilistic end
abstract type JointProbabilisticComposite <: JointProbabilistic end
abstract type IntervalComposite <: Interval end
abstract type  UnsupervisedComposite <: Unsupervised end
abstract type  StaticComposite <: Static end

# surrogate composite models:
struct DeterministicSurrogate <: Deterministic end
struct ProbabilisticSurrogate <: Probabilistic end
struct JointProbabilisticSurrogate <: JointProbabilistic end
struct IntervalSurrogate <: Interval end
struct UnsupervisedSurrogate <: Unsupervised end
struct StaticSurrogate <: Static end

Deterministic() = DeterministicSurrogate()
Probabilistic() = ProbabilisticSurrogate()
JointProbabilistic() = JointProbabilisticSurrogate()
Interval() = IntervalSurrogate()
Unsupervised() = UnsupervisedSurrogate()
Static() = StaticSurrogate()

const SupervisedComposite =
    Union{DeterministicComposite,ProbabilisticComposite,JointProbabilisticComposite,IntervalComposite}

const SupervisedSurrogate =
    Union{DeterministicSurrogate,ProbabilisticSurrogate,JointProbabilisticSurrogate,IntervalSurrogate}

const Surrogate = Union{DeterministicSurrogate,
                        ProbabilisticSurrogate,
                        JointProbabilisticSurrogate,
                        IntervalSurrogate,
                        UnsupervisedSurrogate,
                        StaticSurrogate}

const Composite = Union{SupervisedComposite,UnsupervisedComposite,
                        StaticComposite}

for T in [:DeterministicComposite,
          :ProbabilisticComposite,
          :JointProbabilisticComposite,
          :IntervalComposite,
          :UnsupervisedComposite,
          :StaticComposite,
          :DeterministicSurrogate,
          :JointProbabilisticSurrogate,
          :ProbabilisticSurrogate,
          :IntervalSurrogate,
          :UnsupervisedSurrogate,
          :StaticSurrogate]
    eval(:(MMI.is_wrapper(::Type{$T}) = true))
    eval(:(MMI.package_name(::Type{$T}) = "MLJBase"))
    str = string(T)
    eval(:(MMI.load_path(::Type{$T}) = "MLJBase."*$str))
end

# aliases for legacy code:
const DeterministicNetwork = DeterministicComposite
const ProbabilisticNetwork = ProbabilisticComposite
const UnsupervisedNetwork = UnsupervisedComposite


text:
module QXTools
using Logging
using Reexport
using QXZoo
using QXTns
import DataStructures: OrderedDict

export TensorNetworkCircuit, add_input!, add_output!

# circuit manipulation
include("circuits/circuits.jl")
# circuit to tn conversion
include("tn_conversion.jl")
# contraction planning
include("contraction_planning.jl")
# dsl and parameter files
include("compute_graph/compute_graph.jl")
# functions for output arguments and parsing
include("outputs.jl")
# simulation utilities
include("simulation.jl")

@reexport using QXTools.Circuits

end


text:
println("\n\n\nRunning tests ...")

push!(LOAD_PATH, joinpath(@__DIR__, "../.."))
using InformedDifferentialEvolution
using Test, Statistics

optimum = [-2 4]'
costf(x) = (r = sqrt(sum((x-optimum).^2)); sin(r)^2+r/2)
costf2(x, ignore) = (r = sqrt(sum((x-optimum).^2)); sin(r)^2+r/2)

best, info = DE(costf, [-10,-10], [10,10])
@test best == optimum

best, info = DE(costf2, [-10,-10], [10,10], data = 1)
@test best == optimum

f(x) = sum(abs.(x))
predictor(r,pop,costs) = r[:] = vcat(mean(pop,dims=1), mean(pop,dims=1))
mi = [-1,-1]; ma = [1,1]

best, info = DE(f, mi, ma, predictors = [predictor, :default])
@test maximum(abs.(best)) < 0.000001

function g(x, previouscost = Inf)
    r = 0.
    for i = -5:5
        r += abs.(x[1])+abs.(x[2]+i)
        if r > previouscost
            break
        end
    end
    r
end
best, info = DE(g, [-100,-100], [100,100])
@test best == [0 0]'

println("   done running tests!")


text:
#=

# Dispatch

Some API functions cannot be called directly from the Vulkan library. In particular, extension functions *must* be called through a pointer obtained via API commands. In addition, most functions can be made faster by calling directly into their function pointer, instead of going through the loader trampoline which resolves the function pointer every time.

To circumvent that, we provide a handy way for retrieving and calling into function pointers, as well as a thread-safe global dispatch table that can automate this work for you.

## Retrieving function pointers

API Function pointers can be obtained with the [`function_pointer`](@ref) function, using the API function name.

=#

using SwiftShader_jll # hide
using Vulkan
@set_driver :SwiftShader # hide

const instance = Instance([], [])
const device = Device(
    first(unwrap(enumerate_physical_devices(instance))), [DeviceQueueCreateInfo(0, [1.])],
    [],
    [],
)

#=

```@repl dispatch
function_pointer("vkCreateInstance")
function_pointer(instance, "vkDestroyInstance")
function_pointer(device, "vkCreateFence")
```

It is essentially a wrapper around `get_instance_proc_addr` and `get_device_proc_addr`, leveraging multiple dispatch to make it more intuitive to use.

## Providing function pointers

Every wrapper function or handle constructor has a signature which accepts the standard function arguments plus a function pointer to call into. For example, instead of doing

=#

foreach(println, unwrap(enumerate_instance_layer_properties()))

# you can do

fptr = function_pointer("vkEnumerateInstanceLayerProperties")
foreach(println, unwrap(enumerate_instance_layer_properties(fptr)))

# In the case of a handle constructor which calls both a creation and a destruction function, there is an argument for each corresponding function pointer:

fptr_create = function_pointer(device, "vkCreateFence")
fptr_destroy = function_pointer(device, "vkDestroyFence")
Fence(device, fptr_create, fptr_destroy)

#=

## Automatic dispatch

Querying and retrieving function pointers every time results in a lot of boilerplate code. To remedy this, we provide a thread-safe global dispatch table. Upon every API call, it will call the function pointer directly, and will automatically retrieve it if necessary. It has one dispatch table per device, so that multiple devices can be used at the same time.

This feature can be disabled by setting the [preference](@ref Preferences) `USE_DISPATCH_TABLE` to `"false"`.

=#


text:
using Test

@testset "intercpetion" begin
    for i in 1:1000
        Potential_Evaporation = rand(0.1:0.1:5)[1]
        Precipitation = rand(0.1:1:20)[1]
        Interceptioncapacity = rand(0.1:0.1:3)[1]
        Storage = rand(0:0.1:Interceptioncapacity)[1]
        Temp_Thresh =rand(-2:0.1:2)[1]
        Temp = rand(Temp_Thresh + 0.1:1:25)[1]

        #if precipitation is zero, no effective precipitation
        Precipitation = 0.
        Effective_Precipitation, Interception_Evaporation, Interceptionstorage = interception(Potential_Evaporation, Precipitation, Temp, Storage, Interceptioncapacity, Temp_Thresh)
        @test Effective_Precipitation == 0
        @test Interception_Evaporation == min(Potential_Evaporation * 0.5, Storage)
        @test Interceptionstorage == Storage - min(Potential_Evaporation * 0.5, Storage)

        # if precipitation and storage is zero, all outputs are zero
        Precipitation = 0.
        Storage = 0.
        Effective_Precipitation, Interception_Evaporation, Interceptionstorage = interception(Potential_Evaporation, Precipitation, Temp, Storage, Interceptioncapacity, Temp_Thresh)
        @test Effective_Precipitation == 0
        @test Interception_Evaporation == 0
        @test Interceptionstorage == 0

        prec = 14

        # if storage is zero and rain is present, evaporation, and outflow
        Precipitation = rand(0.1:1:20)[1]
        Storage = 0.
        Effective_Precipitation, Interception_Evaporation, Interceptionstorage = interception(Potential_Evaporation, Precipitation, Temp, Storage, Interceptioncapacity, Temp_Thresh)
        @test round(Effective_Precipitation, digits = prec) == round(Precipitation - (Interceptionstorage + Interception_Evaporation), digits = prec)
        # evaporation also occurs on rainy days
        New_Storage = min(Precipitation, Interceptioncapacity)
        @test round(Interception_Evaporation, digits = prec) == round(min(New_Storage, Potential_Evaporation * 0.5), digits = prec)
        @test round(Interceptionstorage, digits= prec) == round(min(Interceptioncapacity, Precipitation) - Interception_Evaporation, digits= prec)
        if Effective_Precipitation > 0
            @test round(Interceptionstorage, digits = prec) == round(Interceptioncapacity - Interception_Evaporation, digits = prec)
        end

        Precipitation = rand(0:1.:20)[1]
        Storage = rand(0:0.1:Interceptioncapacity)[1]
        Effective_Precipitation, Interception_Evaporation, Interceptionstorage = interception(Potential_Evaporation, Precipitation, Temp, Storage, Interceptioncapacity, Temp_Thresh)
        @test round(Effective_Precipitation, digits= prec) == round(Precipitation - (Interceptionstorage + Interception_Evaporation - Storage), digits=prec)
        New_Storage = min(Precipitation + Storage, Interceptioncapacity)
        @test round(Interception_Evaporation, digits = prec) == round(min(New_Storage, Potential_Evaporation * 0.5), digits = prec)
        @test round(Interceptionstorage, digits=prec) == round(min(Interceptioncapacity, Precipitation + Storage) - Interception_Evaporation, digits=prec)
        if Effective_Precipitation > 0
            @test round(Interceptionstorage, digits = prec) == round(Interceptioncapacity - Interception_Evaporation, digits =prec)
        end
        # if potential evaporation is 0, evaporation will always be zero
        Potential_Evaporation = 0.
        Precipitation = 0.
        Effective_Precipitation, Interception_Evaporation, Interceptionstorage = interception(Potential_Evaporation, Precipitation, Temp, Storage, Interceptioncapacity, Temp_Thresh)
        @test Interception_Evaporation == 0

        # Interceptioncapacity = 0
        # Storage = rand(0:0.1:Interceptioncapacity)[1]
        # Effective_Precipitation, Interception_Evaporation, Interceptionstorage = interception(Potential_Evaporation, Precipitation, Temp, Storage, Interceptioncapacity, Temp_Thresh)
        # @test Effective_Precipitation == Precipitation
        # @test Interceptionstorage == 0
        # @test Interception_Evaporation == 0
        Potential_Evaporation = rand(0.1:0.1:5)[1]
        Interceptioncapacity = 0.
        Storage = rand(0:0.1:Interceptioncapacity)[1]
        Precipitation = rand(0.1:1:20)[1]
        Effective_Precipitation, Interception_Evaporation, Interceptionstorage = interception(Potential_Evaporation, Precipitation, Temp, Storage, Interceptioncapacity, Temp_Thresh)
        @test Effective_Precipitation == Precipitation
        @test Interceptionstorage == 0
        @test Interception_Evaporation == 0

        # if temperature below threshold there is no precipitation and evaporation
        Potential_Evaporation = rand(0.1:0.1:10)[1]
        Precipitation = rand(0.1:1:20)[1]
        Interceptioncapacity = rand(0.1:0.1:3)[1]
        Storage = rand(0:0.1:Interceptioncapacity)[1]
        Temp_Thresh =rand(-2:0.1:2)[1]
        Temp = rand(-10:0.1:Temp_Thresh)[1]
        Effective_Precipitation, Interception_Evaporation, Interceptionstorage = interception(Potential_Evaporation, Precipitation, Temp, Storage, Interceptioncapacity, Temp_Thresh)
        @test Effective_Precipitation == 0
        @test Interception_Evaporation == 0
        @test Interceptionstorage == Storage

    end

end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
#=
ported from:
COTD Entry submitted by John W. Ratcliff [jratcliff@verant.com]

THIS IS A CODE SNIPPET WHICH WILL EFFICIEINTLY TRIANGULATE ANY
POLYGON/CONTOUR (without holes) AS A STATIC CLASS.  THIS SNIPPET
IS COMPRISED OF 3 FILES, TRIANGULATE.H, THE HEADER FILE FOR THE
TRIANGULATE BASE CLASS, TRIANGULATE.CPP, THE IMPLEMENTATION OF
THE TRIANGULATE BASE CLASS, AND TEST.CPP, A SMALL TEST PROGRAM
DEMONSTRATING THE USAGE OF THE TRIANGULATOR.  THE TRIANGULATE
BASE CLASS ALSO PROVIDES TWO USEFUL HELPER METHODS, ONE WHICH
COMPUTES THE AREA OF A POLYGON, AND ANOTHER WHICH DOES AN EFFICENT
POINT IN A TRIANGLE TEST.
SUBMITTED BY JOHN W. RATCLIFF (jratcliff@verant.com) July 22, 2000

 **********************************************************************/
 ************ HEADER FILE FOR TRIANGULATE.H ***************************/
 **********************************************************************/
 *****************************************************************/
 ** Static class to triangulate any contour/polygon efficiently **/
 ** You should replace Vector2d with whatever your own Vector   **/
 ** class might be.  Does not support polygons with holes.      **/
 ** Uses STL vectors to represent a dynamic array of vertices.  **/
 ** This code snippet was submitted to FlipCode.com by          **/
 ** John W. Ratcliff (jratcliff@verant.com) on July 22, 2000    **/
 ** I did not write the original code/algorithm for this        **/
 ** this triangulator, in fact, I can't even remember where I   **/
 ** found it in the first place.  However, I did rework it into **/
 ** the following black-box static class so you can make easy   **/
 ** use of it in your own code.  Simply replace Vector2d with   **/
 ** whatever your own Vector implementation might be.           **/
 *****************************************************************/
=#

function area{N, T}(contour::AbstractVector{Point{N, T}})
    n = length(contour)
    A = zero(T)
    p=n; q=1
    while q <= n
        A += cross(contour[p], contour[q])
        p = q; q +=1
    end
    return A*T(0.5)
end

"""
 InsideTriangle decides if a point P is Inside of the triangle
 defined by A, B, C.
"""
function InsideTriangle{T<:Point}(A::T, B::T, C::T, P::T)
    a = C-B; b = A-C; c = B-A
    ap = P-A; bp = P-B; cp = P-C
    a_bp = a[1]*bp[2] - a[2]*bp[1];
    c_ap = c[1]*ap[2] - c[2]*ap[1];
    b_cp = b[1]*cp[2] - b[2]*cp[1];
    t0 = zero(eltype(T))
    return ((a_bp >= t0) && (b_cp >= t0) && (c_ap >= t0))
end

function snip{N, T}(
        contour::AbstractVector{Point{N, T}}, u, v, w, n, V
    )
    A = contour[V[u]]
    B = contour[V[v]]
    C = contour[V[w]]
    x = (
        ((B[1]-A[1])*(C[2]-A[2])) -
        ((B[2]-A[2])*(C[1]-A[1]))
    )
    if 0.0000000001f0 > x
        return false
    end

    for p=1:n
        ((p == u) || (p == v) || (p == w)) && continue;
        P = contour[V[p]]
        if InsideTriangle(A, B, C, P)
            return false;
        end
    end
    return true;
end

"""
Triangulates a Polygon given as a `contour`::AbstractArray{Point} without holes.
It will return a Vector{`facetype`}, defining indexes into `contour`
"""
function polygon2faces{P<:Point}(
        contour::AbstractArray{P}, facetype=GLTriangle
    )
    #= allocate and initialize list of Vertices in polygon =#
    result = facetype[]

    n = length(contour)
    if n < 3
        error("Not enough points in the contour. Found: $contour")
    end
    #= we want a counter-clockwise polygon in V =#
    if 0 < area(contour)
        V = Int[i for i=1:n]
    else
        V = Int[i for i=n:-1:1]
    end

    nv = n

    #=  remove nv-2 Vertices, creating 1 triangle every time =#
    count = 2*nv   #= error detection =#
    v=nv
    while nv>2
        #= if we loop, it is probably a non-simple polygon =#
        if 0 >= count
            return result
        end
        count -= 1

        #= three consecutive vertices in current polygon, <u,v,w> =#
        u = v; (u > nv) && (u = 1) #= previous =#
        v = u+1; (v > nv) && (v = 1) #= new v =#
        w = v+1; (w > nv) && (w = 1) #= next =#
        if snip(contour, u, v, w, nv, V)
            #= true names of the vertices =#
            a = V[u]; b = V[v]; c = V[w];
            #= output Triangle =#
            push!(result, facetype(Triangle{Int}(a, b, c)))
            #= remove v from remaining polygon =#
            s = v; t = v+1
            while t<=nv
                V[s] = V[t]
                s += 1; t += 1
            end
            nv -= 1
            #= resest error detection counter =#
            count = 2*nv
        end
    end

    return result
end

function topoint{T}(::Type{Point{3, T}}, p::Point{2, T})
    Point{3, T}(p[1], p[2], T(0))
end
function topoint{T}(::Type{Point{3, T}}, p::Point{3, T})
    p
end
function topoint{T}(::Type{Point{2, T}}, p::Point{3, T})
    Point{2, T}(p[1], p[2])
end

@compat function (::Type{M}){M<:AbstractMesh, P<:Point}(
        points::AbstractArray{P}
    )
    faces = polygon2faces(points, facetype(M))
    VT = vertextype(M)
    GLNormalMesh(faces=faces, vertices=VT[topoint(VT, p) for p in points])
end


text:
#=
MIT License

Copyright (c) 2017 Joseph R. Peterson

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

# External Libraries
using DataTables
using ArgParse

# Create command line parser and parse arguments
function parse_cmd_line()
	argSettings = ArgParseSettings()
	@add_arg_table argSettings begin
		"--expressionData", "-i"
			help = "A tab/comma separated file containing gene expression data."
			arg_type = String
			required = true
		"--outputDir", "-o"
			help = "Name of directory into which results are stored."
			arg_type = String
			default = "out"
		"--threshold", "-t"
			help = "A threshold on the minimum coefficient of variation for genes to be considered by the algorithm (pre-algorithm filtering)."
			arg_type = Float64
			default = 0.0
		"--significant", "-s"
			help = "Significant interactions are those 's' standard deviations away from mean (post-algorithm filtering)."
			arg_type = Float64
			default = 2.0
		"--plot", "-p"
			help = "Plot intermediate results"
			default = false
			action = :store_true
	end

	# Actually pargse aguments
	return parse_args(argSettings)
end

# Main computes the network
function main()
	# Get command line arguments
	args = parse_cmd_line()

	# Read data files
	println("Reading inputs...")
	expressionData = readtable(args["expressionData"])

	# Run algorithm
	println("Computing couplings...")
	genes, couplingStrength = plogp(expressionData, args["threshold"])

	# Cull network
	println("Thresholding significant interactions...")

	# Plot Data (if requested)
	if args["plot"]
		println("Plotting...")
	end

	# Save resulting weight and graph files
	println("Saving results...")

	# Successfully completed
	println("Done!")
end
main()


text:
module Gaugefields

using Requires
include("./output/verboseprint_mpi.jl")
#include("./output/verboseprint.jl")
include("./SUN_generator.jl")
include("./autostaples/wilsonloops.jl")
include("./AbstractGaugefields.jl")
include("./output/io.jl")
include("./output/ildg_format.jl")
include("./output/bridge_format.jl")
include("./autostaples/Loops.jl")
include("./smearing/Abstractsmearing.jl")

include("./action/GaugeActions.jl")
include("./heatbath/heatbathmodule.jl")
include("./smearing/gradientflow.jl")

function __init__()
    @require MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" begin
        import .AbstractGaugefields_module:identityGaugefields_4D_wing_mpi,
                Gaugefields_4D_wing_mpi,calc_rank_and_indices,barrier,comm,
                setvalue!
    end
end

# Write your package code here.
import .Verboseprint_mpi:Verbose_print,println_verbose_level1,println_verbose_level2,println_verbose_level3
import .AbstractGaugefields_module:AbstractGaugefields,IdentityGauges,RandomGauges,Oneinstanton,calculate_Plaquette,
                                    calculate_Polyakov_loop,map_U!,evaluate_gaugelinks_evenodd!,normalize!,normalize3!,normalizeN!,
                                    shift_U,evaluate_gaugelinks!,
                                    TA_Gaugefields,Abstractfields,Staggered_Gaugefields,staggered_U,
                                    clear_U!,set_wing_U!,Traceless_antihermitian!,initialize_TA_Gaugefields,
                                    substitute_U!,Traceless_antihermitian!,Traceless_antihermitian

#=
                                    import .AbstractGaugefields_module:AbstractGaugefields,identitymatrix,Abstractfields,
                                    shift_U,construct_staple!,set_wing_U!,
                                    calculate_Plaquette,substitute_U!,calculate_Polyakov_loop,construct_gauges,
                                    Gaugefields_4D_wing,
                                    identityGaugefields_4D_wing,
                                    add_force!,exp_aF_U!,clear_U!,add_U!,exptU!,
                                    Traceless_antihermitian!,Traceless_antihermitian,Generator,
                                    Staggered_Gaugefields,staggered_U,
                                    Traceless_antihermitian_add!,
                                    IdentityGauges,RandomGauges,Oneinstanton,
                                    construct_Λmatrix_forSTOUT!,
                                    evaluate_gaugelinks_evenodd!,
                                    map_U!
                                    =#
import .Loops_module:Loops,evaluate_loops,calc_large_wilson_loop!,evaluate_loops!
import .Wilsonloops_module:Wilson_loop_set,make_staples,Wilson_loop_set,
                                    make_cloverloops,Tensor_derivative_set, make_loops,
                                    make_plaq_staple,make_links,make_plaq,
                                    make_loopforactions,make_plaqloops,make_rectloops,make_polyakovloops,
                                    make_plaq_staple_prime,
                                    calc_coordinate,make_plaq_staple_prime,calc_shift,
                                    Tensor_wilson_lines_set,Tensor_wilson_lines,Tensor_derivative_set,
                                    get_leftstartposition,get_rightstartposition,Wilson_loop,calc_loopset_μν_name,
                                    make_originalactions_fromloops
import .Abstractsmearing_module:Abstractsmearing,Nosmearing,Stoutsmearing,calc_smearedU,
                                construct_smearing,back_prop,CovNeuralnet#gradientflow!
import .ILDG_format:ILDG,load_gaugefield!,save_binarydata
import .heatbath_module:SU2update_KP!,SUNupdate_matrix!,SU3update_matrix!,heatbath!,
            Heatbath,Heatbath_update,overrelaxation!
import .Bridge_format:save_textdata,load_BridgeText!
import Wilsonloop:loops_staple
import .Abstractsmearing_module:STOUT_Layer,CovNeuralnet,calc_smearedU, construct_smearing,
        set_parameters,get_parameter_derivatives,apply_smearing_U
import .SUN_generator:Generator
import .Gradientflow_module:Gradientflow,Gradientflow_general,flow!,get_tempG,get_eps
#import .Verbose_print:Verbose_level,Verbose_3,Verbose_2,Verbose_1,println_verbose3,println_verbose2,println_verbose1,
#    print_verbose1,print_verbose2,print_verbose3

import .ILDG_format:ILDG,load_gaugefield,load_gaugefield!,save_binarydata,load_binarydata!
import .IOmodule:saveU,loadU,loadU!

import .AbstractGaugefields_module:AbstractGaugefields,identitymatrix,Abstractfields,
            shift_U,construct_staple!,set_wing_U!,
            calculate_Plaquette,substitute_U!,calculate_Polyakov_loop,construct_gauges,
            Gaugefields_4D_wing,
            identityGaugefields_4D_wing,
            add_force!,exp_aF_U!,clear_U!,add_U!,exptU!,
            Traceless_antihermitian!,Traceless_antihermitian,Generator,
            Staggered_Gaugefields,staggered_U,
            Traceless_antihermitian_add!,
            IdentityGauges,RandomGauges,Oneinstanton,Initialize_4DGaugefields,
            construct_Λmatrix_forSTOUT!,
            evaluate_gaugelinks_evenodd!,
            map_U!,initialize_TA_Gaugefields,gauss_distribution!,Initialize_Gaugefields,
            construct_Adjoint_rep_Gaugefields,get_myrank,getvalue,get_nprocs,
            write_to_numpyarray,map_U_sequential!
import Wilsonloop:make_loops_fromname
import .GaugeAction_module:GaugeAction,evaluate_GaugeAction_untraced!,evaluate_GaugeAction_untraced,calc_dSdUμ,calc_dSdUμ!,get_temporary_gaugefields,evaluate_GaugeAction

export IdentityGauges,RandomGauges,Oneinstanton,calculate_Plaquette,calculate_Polyakov_loop
export ILDG,load_gaugefield!,save_binarydata
export SU2update_KP!,SUNupdate_matrix!,SU3update_matrix!
export map_U!
export evaluate_gaugelinks_evenodd!,normalize!,normalize3!,normalizeN!
export loops_staple
export save_textdata,load_BridgeText!
export shift_U,evaluate_gaugelinks!,Gradientflow,flow!
export heatbath!,Heatbath
export STOUT_Layer,CovNeuralnet,calc_smearedU,make_loops_fromname
export GaugeAction,evaluate_GaugeAction_untraced!,evaluate_GaugeAction_untraced,calc_dSdUμ,calc_dSdUμ!
export initialize_TA_Gaugefields,gauss_distribution!
export exptU!,get_temporary_gaugefields,Traceless_antihermitian_add!,evaluate_GaugeAction,substitute_U!,set_wing_U!,Traceless_antihermitian!
export Initialize_Gaugefields,back_prop
export Initialize_4DGaugefields
export set_parameters,get_parameter_derivatives,apply_smearing_U
export construct_Adjoint_rep_Gaugefields
export get_myrank,getvalue,get_nprocs,Gradientflow_general
export Heatbath_update
export println_verbose_level1,println_verbose_level2,println_verbose_level3
export overrelaxation!
export AbstractGaugefields,Traceless_antihermitian
export write_to_numpyarray,map_U_sequential!
export load_binarydata!

end


text:
using Test
using ProbNumDiffEq
using DiffEqDevTools

# Simple linear problem
linear(u, p, t) = p .* u
linear_jac(u, p, t) = p
linear_analytic(u0, p, t) = @. u0 * exp(p * t)
prob = ODEProblem(
    ODEFunction(linear, jac=linear_jac, analytic=linear_analytic),
    [big(1 / 2)],
    big.((0.0, 1.0)),
    big(1.01),
)

# Step-sizes
dts1 = 1 .// 2 .^ (7:-1:4)
dts2 = 1 .// 2 .^ (10:-1:8)

@testset "EK0(order=$q) convergence" for q in 1:5
    sim = test_convergence(dts1, prob, EK0(order=q))
    @test sim.𝒪est[:final] ≈ q + 1 atol = 0.3
    @test sim.𝒪est[:l2] ≈ q + 1 atol = 0.3
    @test sim.𝒪est[:l∞] ≈ q + 1 atol = 0.3
end
@testset "EK0(order=$q) convergence" for q in [7, 8]
    sim = test_convergence(dts2, prob, EK0(order=q))
    @test sim.𝒪est[:final] ≈ q + 1 atol = 1.5
    @test sim.𝒪est[:l2] ≈ q + 1 atol = 1
    @test sim.𝒪est[:l∞] ≈ q + 1 atol = 1
end

@testset "EK1(order=$q) convergence" for q in 1:5
    sim = test_convergence(dts1, prob, EK1(order=q))
    @test sim.𝒪est[:final] ≈ q + 1 atol = 0.55
    @test sim.𝒪est[:l2] ≈ q + 1 atol = 0.5
    @test sim.𝒪est[:l∞] ≈ q + 1 atol = 0.5
end
@testset "EK1(order=$q) convergence" for q in [7, 8, 10]
    sim = test_convergence(dts2, prob, EK1(order=q))
    @test sim.𝒪est[:final] ≈ q + 1 atol = 1
    @test sim.𝒪est[:l2] ≈ q + 1 atol = 1
    @test sim.𝒪est[:l∞] ≈ q + 1 atol = 1
end