text:
#declare sparse matrix based DAWG type
struct SparseMatrixDAWG
    m::SparseMatrixCSC
    rownum::Int64
    colnum::Int64
    edgedict::Dict{Char,Int32}
end

#declare double array based DAWG type
struct DoubleArrayDAWG
    b::Array{Int32,1}
    check::Array{Int32,1}
    rownum::Int64
    colnum::Int64
    edgedict::Dict{Char,Int32}
end

#DAWG object definition
mutable struct DawgNode
    id::Int32
    final::Bool
    edges::Dict
    count::Int32
end

function DawgNode(id = 0)
    return DawgNode(id, false, Dict{Char,DawgNode}(), -1)
end

function str(self::DawgNode)
    arr = IOBuffer()
    if self.final
        write(arr, '1')
    else
        write(arr, '0')
    end
    for (label, node) in sort(collect(zip(keys(self.edges), values(self.edges))))
        write(arr, label, '_', node.id, '_')
    end
    return String(take!(arr)) # takebuf_string
end

function numReachable!(self::DawgNode)
    if self.count > 0
        return self.count
    end
    count = 0
    if self.final
        count += 1
    end
    for (label, node) in sort(collect(zip(keys(self.edges), values(self.edges))))
        count += numReachable!(node)
    end
    self.count = count
    return count
end

#declare DAWG type
mutable struct Dawg
    NextId::Int32
    previousWord::AbstractString
    root::DawgNode
    uncheckedNodes::Array
    data::Array
    minimizedNodes::Dict
end

# Dawg functions
function Dawg()
    return Dawg(1, "", DawgNode(0), [], [], Dict())
end

function insert!(self::Dawg, word)
    word = utf32(word)
    if word <= self.previousWord
        throw("Error: Words must be inserted in alphabetical order: ", word)
    end
    commonPrefix = 0
    for i = 1:min(length(word), length(self.previousWord))
        word[i] == self.previousWord[i] ? commonPrefix += 1 : break
    end
    minimize!(self, commonPrefix) # HERE
    node = length(self.uncheckedNodes) == 0 ? self.root : self.uncheckedNodes[end][3]
    for letter in word[(commonPrefix+1):end]
        nextNode = DawgNode(self.NextId)
        self.NextId += 1
        node.edges[letter] = nextNode
        push!(self.uncheckedNodes, (node, letter, nextNode))
        node = nextNode
    end

    node.final = true
    self.previousWord = word
end

function minimize!(self::Dawg, downTo)
    if length(self.uncheckedNodes) < 1
        return
    end
    for i = length(self.uncheckedNodes):-1:(downTo+1)
        (parent, letter, child) = self.uncheckedNodes[i]
        if haskey(self.minimizedNodes, str(child))
            parent.edges[letter] = self.minimizedNodes[str(child)]
        else
            self.minimizedNodes[str(child)] = child
        end
        pop!(self.uncheckedNodes)
    end
end

function finish!(self::Dawg)
    minimize!(self, 0)
    numReachable!(self.root)
    self.minimizedNodes[str(self.root)] = self.root
end

function lookup(self::Dawg, word)
    word = utf32(word)
    node = self.root
    skipped = 0
    for letter in word
        if !haskey(node.edges, letter)
            return ""
        end
        for (label, child) in sort(collect(zip(keys(node.edges), values(node.edges))))
            if label == letter
                if node.final
                    skipped += 1
                end
                node = child
                break
            end
            skipped += child.count
        end
    end

    if node.final
        println("<--- This is word. ($skipped)")
    end
end

function nodeCount(self::Dawg)
    return length(self.minimizedNodes)
end

function edgeCount(self::Dawg)
    count = 0
    for (key, node) in zip(keys(self.minimizedNodes), values(self.minimizedNodes))
        count += length(node.edges)
    end
    return count
end

#recursively collect states with paths in DAWG
#used when constructing sparse matrix DAWG
#returns nothing but fills statedict, edgedict, donepathdict, and statecollect
function travrec!(node, statedict, edgedict, donepathdict, statecollect)
    for path in keys(node.edges)
        if !haskey(statedict, node.id)
            statedict[node.id] = length(statedict) + 1
            statecollect[node.id] = node
        end
        if !haskey(edgedict, path)
            edgedict[path] = length(edgedict) + 1
        end

        if get(donepathdict, "$(node.id)$path", 0) == node.edges[path].id
            return
        end
        donepathdict["$(node.id)$path"] = node.edges[path].id

        if node.edges[path].final
            if length(node.edges[path].edges) > 0
                travrec!(node.edges[path], statedict, edgedict, donepathdict, statecollect)
            end
        else
            travrec!(node.edges[path], statedict, edgedict, donepathdict, statecollect)
        end
    end
end

#recursively collect states with paths in DAWG
#used when constructing double array DAWG
#returns nothing but fills statedict, edgedict, donepathdict, and statecollect
function travrecBC!(node, statedict, edgedict, donepathdict, statecollect, termnode)
    if !haskey(statedict, node.id)
        statedict[node.id] = length(statedict) + 1
        statecollect[node.id] = node
    end

    if node.final && length(node.edges) == 0
        termnode[node.id] = node.id
        return
    end

    for path in keys(node.edges)
        if !haskey(edgedict, path)
            edgedict[path] = length(edgedict) + 1
        end

        if get(donepathdict, (node.id, path), 0) == node.edges[path].id
            return
        end
        donepathdict[(node.id, path)] = node.edges[path].id

        travrecBC!(
            node.edges[path],
            statedict,
            edgedict,
            donepathdict,
            statecollect,
            termnode,
        )
    end
    return termnode
end

#directly fill sparse matrix with DAWG states and transition information
function fillmatrixStraightIJV(rownum, colnum, statedict, edgedict, statecollect)
    block = colnum
    I = zeros(Int32, block)
    J = zeros(Int32, block)
    V = zeros(Int32, block)

    cnt = 0
    element = 0
    statesize = length(statecollect)
    for k in sort(collect(keys(statecollect)))
        node = statecollect[k]
        cnt += 1
        if cnt % 100000 == 0
            println("Done: $cnt / $statesize")
        end

        if length(I) < (element + block)
            I = vcat(I, zeros(Int32, block))
            J = vcat(J, zeros(Int32, block))
            V = vcat(V, zeros(Int32, block))
        end

        for path in keys(node.edges)
            element += 1
            I[element] = statedict[node.id]
            J[element] = edgedict[path]

            if node.edges[path].final
                if length(node.edges[path].edges) > 0
                    V[element] = -statedict[node.edges[path].id]
                else
                    V[element] = -rownum
                end
                continue
            end
            V[element] = statedict[node.edges[path].id]
        end
    end
    len = findlast(x -> x > 0, I)
    return sparse(I[1:len], J[1:len], V[1:len], rownum, colnum)
end

#perform dictionary match using sparse matrix DAWG
#input: sentence as a valid string, dawgdict as a sparse matrix based DAWG
#set extract = true (default) will display matched word as progressing
function dictionaryMatch(sentence, dawgdict::SparseMatrixDAWG; extract::Bool = true)
    rettuple = tokenizetuple(
        sentence,
        dawgdict.m,
        dawgdict.edgedict,
        dawgdict.rownum,
        dawgdict.colnum,
    )
    if extract
        usentence = utf32(sentence)
        for (from, to) in rettuple
            println(usentence[from:to])
        end
    end
    rettuple
end

#perform dictionary matching scan with spare matrix based DAWG(inner use only)
function tokenizetuple(sentence, m::SparseMatrixCSC, edgedict, rownum, def)
    curstate = 1
    cutat = []
    curpos = 1
    lastbegin = 1
    lastgood = 1
    sentence = utf32(sentence)
    senlen = length(sentence)

    while curpos <= senlen
        @inbounds col = get(edgedict, sentence[curpos], def)
        @inbounds curstate = m[curstate, col]

        if curstate > 0
            curpos += 1
        elseif curstate == 0
            if lastgood > lastbegin
                push!(cutat, (lastbegin, lastgood))
            end
            curstate = 1
            lastbegin = curpos = lastgood += 1
        elseif curstate < 0
            if curstate == -rownum
                push!(cutat, (lastbegin, curpos))
                curstate = 1
                lastbegin = lastgood = curpos += 1
            else
                curstate = -curstate
                lastgood = curpos
                curpos += 1
            end
        end
    end

    ##flush
    if lastgood > lastbegin
        push!(cutat, (lastbegin, lastgood))
    end
    return cutat
end

#perform dictionary match using double array DAWG
#input: sentence as a valid string, dawgdict as a sparse matrix based DAWG
#set extract = true (default) will display matched word as progressing
function dictionaryMatch(sentence, dawgdict::DoubleArrayDAWG; extract::Bool = true)
    rettuple =
        tokenizedoublearrayBC(sentence, dawgdict.edgedict, dawgdict.b, dawgdict.check)
    if extract
        usentence = utf32(sentence)
        for (from, to) in rettuple
            println(usentence[from:to])
        end
    end
    rettuple
end

#perform dictionary matching scan with double array based DAWG(inner use only)
function tokenizedoublearrayBC(sentence, edgedict, b, check)
    curstate = 2
    cutat = []
    curpos = 1
    lastbegin = 1
    lastgood = 1
    sentence = utf32(sentence)
    senlen = length(sentence)

    while curpos <= senlen
        @inbounds col = get(edgedict, sentence[curpos], 0)
        if col == 0
            curstate = -1
        else
            if curstate > 1
                @inbounds p = b[curstate] + col
                @inbounds curstate = check[p] == b[curstate] ? p : -1
            end
        end

        if curstate < 1
            if lastgood > lastbegin
                push!(cutat, (lastbegin, lastgood))
            end
            curstate = 2
            lastbegin = curpos = lastgood += 1
        elseif b[b[curstate]] == -1
            push!(cutat, (lastbegin, curpos))
            curstate = 2
            lastbegin = lastgood = curpos += 1
        elseif b[b[curstate]] == 0
            lastgood = curpos
            curpos += 1
        else
            curpos += 1
        end
    end

    ##flush
    if lastgood > lastbegin
        push!(cutat, (lastbegin, lastgood))
    end
    return cutat
end

#constrct a sparse matrix based DAWG from word list(dictionary)
#return SparseMatrixDAWG type
function makeSparseMatrixDAWG(list)
    println("Making DAWG... Please wait...")
    d = Dawg()
    for w in sort(list)
        insert!(d, w)
    end
    finish!(d)
    statedict = Dict()
    statenum = nodeCount(d)

    edgedict = Dict{Char,Int32}()
    edgenum = edgeCount(d)

    donepathdict = Dict()
    statecollect = Dict{Int32,DawgNode}()

    travrec!(d.root, statedict, edgedict, donepathdict, statecollect)
    colnum = 1 + length(edgedict)
    rownum = 1 + length(statedict)

    m = fillmatrixStraightIJV(rownum, colnum, statedict, edgedict, statecollect)
    println("Done!")
    return SparseMatrixDAWG(m, rownum, colnum, edgedict)
end

#find vacant space in double array (inner use only)
function x_checkBC(check, c_set)
    minc = minimum(c_set)
    firstvac = findnext(x -> x == 0, check, minc + 2)
    bf = firstvac - minc
    runningbf = bf
    if maximum(bf .+ c_set) >= length(check)
        return 0
    end

    while sum(abs.(check[bf.+c_set])) > 0 || bf < 1
        firstvac = findnext(x -> x == 0, check, runningbf)
        if firstvac == 0
            return 0
        end
        #runningbf += 1
        runningbf = firstvac + 1
        bf = firstvac - minc
        if maximum(bf .+ c_set) >= length(check)
            return 0
        end
    end
    #@show "$bf, $c_set"
    check[bf.+c_set] .= bf
    return bf
end

#constrct double array based on data collected from tree based DAWG
function getdoublearrayBC(statedict, edgedict, statecollect, termnode)

    invn = Dict()
    fv = Dict()
    maxchmap = 5 * (1 + length(edgedict))
    check = zeros(Int32, maxchmap)
    b = zeros(Int32, maxchmap)

    check[1] = 1
    b[1] = -1
    check[2] = 0
    b[2] = 2

    for k in keys(termnode)
        fv[abs(statedict[k])] = 1
    end

    #filling check
    cnt = 0
    statesize = length(statecollect)
    for k in sort(collect(keys(statecollect))) #values(statecollect)
        node = statecollect[k]
        if haskey(termnode, node.id)
            continue
        end
        cnt += 1
        if cnt % 10000 == 0
            println("1st Done: $cnt / $statesize")
        end

        c_set = [0]
        if length(node.edges) > 0
            c_set = [c_set; map(x -> edgedict[x], collect(keys(node.edges)))]
        end

        firstvac = 0
        while firstvac <= 0
            firstvac = x_checkBC(check, c_set)
            if firstvac <= 0
                b = vcat(b, zeros(Int32, maxchmap))
                check = vcat(check, zeros(Int32, maxchmap))
            end
        end

        fv[abs(statedict[node.id])] = firstvac
        if !node.final
            b[firstvac] = 2
        end
    end

    #filling base
    cnt = 0
    for node in values(statecollect)
        cnt += 1
        if cnt % 10000 == 0
            println("2nd Done: $cnt / $statesize")
        end

        bf = fv[abs(statedict[node.id])]

        for path in keys(node.edges)
            if haskey(fv, abs(statedict[node.edges[path].id]))
                b[bf+edgedict[path]] = fv[abs(statedict[node.edges[path].id])]
            end
        end
    end

    maxb = maximum(b) + length(edgedict)
    sizediff = maxb - length(b)
    if maxb > length(b)
        b = [b; zeros(Int32, sizediff)]
        check = [check; zeros(Int32, sizediff)]
    end
    return (b, check, fv)
end

#constrct a sparse matrix based DAWG from word list(dictionary)
#return DoubleArrayDAWG type
function makeDoubleArrayDAWG(list)
    println("Making DAWG... Please wait...")
    d = Dawg()
    for w in sort(list)
        insert!(d, w)
    end
    finish!(d)
    statedict = Dict()
    statenum = nodeCount(d)

    edgedict = Dict{Char,Int32}()
    edgenum = edgeCount(d)

    donepathdict = Dict()
    statecollect = Dict{Int32,DawgNode}()

    termnode = Dict()
    travrecBC!(d.root, statedict, edgedict, donepathdict, statecollect, termnode)
    colnum = 1 + length(edgedict)
    rownum = 1 + length(statedict)

    (b, check, fv) = getdoublearrayBC(statedict, edgedict, statecollect, termnode)
    println("Done!")
    return DoubleArrayDAWG(b, check, rownum, colnum, edgedict)
end


text:
# ----------------------------------------------------
# --- Import modules
# ----------------------------------------------------
using DigitalComm
using Test
# ----------------------------------------------------
# --- Tests
# ----------------------------------------------------
# Note --> The mapping system is described in bitMapping.jl

println("Tests for Hard symbol demapper");
for alphabet = [2;4;16;64;256]
	@testset  "$alphabet-QAM" begin
		# Minimal Tx
		mcs	  = alphabet;		  # Constellation size
		n	  = Int(log2(mcs));		  # Bit per symbol
		N	  = 4096;	  # Numbver of test symbols
		bitSeq	= genBitSequence(N*n);
		qamSeq	= bitMappingQAM(mcs,bitSeq);
		# Minimal Rx
		bitDec	= bitDemappingQAM(mcs,qamSeq);
		# Checking size of decoding vector
		@test length(bitDec) == N*n;
		# Checking that bang and classic methods are the same
		buff  = zeros(UInt8,N*n);
		bitDemappingQAM!(buff,mcs,qamSeq);
		@test all(bitDec .== buff);
		# Checking that decoded method is initial method
		@test all( bitDec .≈ bitSeq)
	end
end

ext: jl


text:
function onecoldencoding(M::Matrix)
    n = size(M,2)
    onecold = zeros(n)
    for sample in 1:n
        onecold[sample] = findall(M[:,sample] .== maximum(M[:,sample]))[1]
    end
    return onecold
end

function accuracy(predictor, dataset)
    return mean(onecoldencoding(predictor(dataset[1][1])) .== onecoldencoding(dataset[1][2]))
end

function nettrain(mi::Float64, x_train::Matrix, y_train::Matrix, x_test::Matrix, y_test::Matrix, niter::Int)

    traindataset = [(x_train, y_train)] # batch
    testdataset = [(x_test, y_test)]
    opt = Descent(mi)
    predictor = Chain(Dense(size(x_train, 1), size(y_train, 1), sigmoid), softmax)
    parameters = params(predictor)
    loss(x, y) = Flux.Losses.mse(predictor(x), y)

    train!(loss, parameters, ncycle(traindataset, niter), opt, cb = throttle(() -> println("Training -- ", "Loss: ", loss(x_train, y_train), " train accuracy: ", accuracy(predictor, traindataset), "; test accuracy: ", accuracy(predictor, testdataset)), 2))

    return predictor, parameters
end


text:
using EvoDynamics
using BenchmarkTools

param_file = "../test/params.yml"

a = @benchmark runmodel(param_file)
display(a)


text:
using ChainRulesTestUtils
using Documenter

makedocs(
    modules=[ChainRulesTestUtils],
    format=Documenter.HTML(prettyurls=false, assets=["assets/chainrules.css"]),
    sitename="ChainRulesTestUtils",
    authors="JuliaDiff contributors",
    strict=true,
    checkdocs=:exports,
)

const repo = "github.com/JuliaDiff/ChainRulesTestUtils.jl.git"
deploydocs(
    repo=repo,
    push_preview=true,
)


text:
using Gadfly, DataArrays, RDatasets
using Base.Dates

economics = dataset("HistData", "Prostitutes")
# NOTE: I know these aren't unix times, but I'm not sure what they are, and this
# is just a test so it doesn't matter.
dates = DateTime[unix2datetime(d) for d in economics[:Date]]
economics[:Date] = dates

p = plot(economics, x=:Date, y=:Count, Geom.line)


text:
@testset "Cycle Basis" begin

    function evaluate(x,y)
        @test length(x) == length(y)
        for i=1:length(x)
            @test isempty(setdiff(x[i],y[i]))
            @test isempty(setdiff(x[i],y[i]))
        end
    end

    # No Edges
    ex = Graph(1)
    expected_cyclebasis = Array{Int64,1}[]
    for g in testgraphs(ex)
        ex_cyclebasis = @inferred cycle_basis(g)
        @test isempty(ex_cyclebasis)
    end

    # Only one self-edge
    elist = [(1,1)]
    ex = Graph(SimpleEdge.(elist))
    expected_cyclebasis = Array{Int64,1}[[1]]
    for g in testgraphs(ex)
        ex_cyclebasis = cycle_basis(g)
        evaluate(ex_cyclebasis,expected_cyclebasis)
    end

    # Graph with one cycle
    elist = [(1,2),(2,3),(3,4),(4,1),(1,5)]
    ex = Graph(SimpleEdge.(elist))
    expected_cyclebasis = Array{Int64,1}[
        [1,2,3,4] ]
    for g in testgraphs(ex)
        ex_cyclebasis = cycle_basis(g)
        evaluate(ex_cyclebasis,expected_cyclebasis)
    end

    # Graph with 2 of 3 cycles forming a basis
    elist = [(1,2),(1,3),(2,3),(2,4),(3,4)]
    ex = Graph(SimpleEdge.(elist))
    expected_cyclebasis = Array{Int64,1}[
        [2,3,4],
        [2,1,3] ]
    for g in testgraphs(ex)
        ex_cyclebasis = cycle_basis(g)
        evaluate(ex_cyclebasis,expected_cyclebasis)
    end

    # Testing root argument
    elist = [(1,2),(1,3),(2,3),(2,4),(3,4),(1,5),(5,6),(6,4)]
    ex = Graph(SimpleEdge.(elist))
    expected_cyclebasis = Array{Int64,1}[
        [2, 4, 3],
        [1, 5, 6, 4, 3],
        [1, 2, 3] ]
    for g in testgraphs(ex)
        ex_cyclebasis = @inferred cycle_basis(g,3)
        evaluate(ex_cyclebasis,expected_cyclebasis)
    end
end


text:
# Evolution Strategy
# ==================
#
# Implementation: (μ/ρ(+/,)λ)-ES
#
# μ is the number of parents
# ρ ≤ μ the mixing number (i.e., the number of parents involved in the procreation of an offspring)
# λ is the number of offspring.
#
# Comma-selection (μ<λ must hold): parents are deterministically selected from the set of the offspring
# Plus-selection: parents are deterministically selected from the set of both the parents and offspring
#
function es(  objfun::Function, N::Int;
    initPopulation::Individual = ones(N),
    initStrategy::Strategy = strategy(),
    recombination::Function = (rs->rs[1]),
    srecombination::Function = (ss->ss[1]),
    mutation::Function = ((r,m)->r),
    smutation::Function = (s->s),
    termination::Function = (x->false),
    μ::Integer = 1,
    ρ::Integer = μ,
    λ::Integer = 1,
    selection::Symbol = :plus,
    iterations::Integer = N*100,
    verbose = false, debug = false,
    interim = false,
    showBar = false)

    @assert ρ <= μ "Number of parents involved in the procreation of an offspring should be no more then total number of parents"
    if selection == :comma
        @assert μ < λ "Offspring population must be larger then parent population"
    end

    store = Dict{Symbol,Any}()

    if showBar
        bar(0.0, iterations)
    end

    # Initialize parent population
    individual = getIndividual(initPopulation, N)
    population = fill(individual, μ)
    fitness = zeros(μ)
    for i in 1:μ
        if isa(initPopulation, Vector)
            population[i] = initPopulation.*rand(N)
        elseif isa(initPopulation, Matrix)
            population[i] = initPopulation[:, i]
        else # Creation function
            population[i] = initPopulation(N)
        end
        fitness[i] = objfun(population[i])
        debug && println("INIT $(i): $(population[i]) : $(fitness[i])")
    end
    offspring = Array{typeof(individual)}(undef, λ)
    fitoff = fill(Inf, λ)
    stgpop = fill(initStrategy, μ)
    stgoff = fill(initStrategy, λ)

    keep(interim, :fitness, copy(fitness), store)

    # Generation cycle
    count = 0
    while true

        for i in 1:λ
            # Recombine the ρ selected parents to form a recombinant individual
            if ρ == 1
                j = rand(1:μ)
                recombinantStrategy = stgpop[j]
                recombinant = copy(population[j])
            else
                idx = randperm(μ)[1:ρ]
                recombinantStrategy = srecombination(stgpop[idx])
                recombinant = recombination(population[idx])
            end

            # Mutate the strategy parameter set of the recombinant
            stgoff[i] = smutation(recombinantStrategy)

            # Mutate the objective parameter set of the recombinant using the mutated strategy parameter set
            # to control the statistical properties of the object parameter mutation
            offspring[i] = mutation(recombinant, stgoff[i])

            # Evaluate fitness
            fitoff[i] = objfun(offspring[i])
        end

        # Select new parent population
        if selection == :plus
            idx = sortperm(vcat(fitness, fitoff))[1:μ]
            skip = idx[idx.<=μ]
            for i = 1:μ
                if idx[i] ∉ skip
                    ii = idx[i] - μ
                    population[i] = offspring[ii]
                    stgpop[i] = stgoff[ii]
                    fitness[i] = fitoff[ii]
                end
            end
        else
            idx = sortperm(fitoff)[1:μ]
            population = offspring[idx]
            stgpop = stgoff[idx]
            fitness = fitoff[idx]
        end
        keep(interim, :fitness, copy(fitness), store)
        keep(interim, :bestFitness, copy(fitness[1]), store)
        keep(interim, :fitoff, copy(fitoff), store)
        keep(interim, :bestFitoff, copy(fitoff[1]), store)
        keep(interim, :individual, copy(population[1]), store)

        if showBar
            bar(count, iterations)
        end

        # termination condition
        count += 1
        if count == iterations || termination(stgpop[1])
            if showBar
                bar(iterations, iterations)
            end
            break
        end
        verbose && println("BEST: $(fitness[1]): $(stgpop[1])")
    end
    return population[1], fitness[1], count, store
end


text:
module Transient2nOrderFEOperatorsTests

using Gridap
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using Test

# Analytical functions
u(x,t) = (1.0-x[1])*x[1]*(t^2+3.0)
u(t::Real) = x -> u(x,t)
v(t::Real) = ∂t(u)(t)
a(t::Real) = ∂tt(u)(t)
f(t) = x -> ∂tt(u)(x,t) + ∂t(u)(x,t) - Δ(u(t))(x)

domain = (0,1)
partition = (2,)
model = CartesianDiscreteModel(domain,partition)

order = 2
reffe = ReferenceFE(lagrangian,Float64,order)
V0 = FESpace(
  model,
  reffe,
  conformity=:H1,
  dirichlet_tags="boundary")
U = TransientTrialFESpace(V0,u)

Ω = Triangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)

m(utt,v) = ∫(v*utt)dΩ
c(ut,v) = ∫(v*ut)dΩ
a(u,v) = ∫(∇(v)⊙∇(u))dΩ
b(v,t) = ∫(v*f(t))dΩ

res(t,(u,ut,utt),v) = m(utt,v) + c(ut,v) + a(u,v) - b(v,t)
jac(t,(u,ut,utt),du,v) = a(du,v)
jac_t(t,(u,ut,utt),dut,v) = c(dut,v)
jac_tt(t,(u,ut,utt),dutt,v) = m(dutt,v)

op = TransientFEOperator(res,jac,jac_t,jac_tt,U,V0)

t0 = 0.0
tF = 1.0
dt = 0.1
γ = 0.5
β = 0.25

U0 = U(0.0)
uh0 = interpolate_everywhere(u(0.0),U0)
vh0 = interpolate_everywhere(v(0.0),U0)
ah0 = interpolate_everywhere(a(0.0),U0)

ls = LUSolver()
odes = Newmark(ls,dt,γ,β)
solver = TransientFESolver(odes)

sol_t = solve(solver,op,(uh0,vh0,ah0),t0,tF)

l2(w) = w*w

tol = 1.0e-6
_t_n = t0

for (uh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  @test tn≈_t_n
  e = u(tn) - uh_tn
  el2 = sqrt(sum( ∫(l2(e))dΩ ))
  @test el2 < tol
end

end


text:
# # Hypercube

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Extension/hypercube.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Extension/hypercube.ipynb)
# **Contributed by**: Benoît Legat

# Given a Sum-of-Squares constraint on an algebraic set:
# ```math
# g_1(x) = , \ldots, g_m(x) = 0 \Rightarrow p(x) \ge 0.
# ```
# We can either use the certificate:
# ```math
# p(x) = s(x) + \lambda_1(x) g_1(x) + \cdots + \lambda_m(x) g_m(x), s_0(x) \text{ is SOS},
# ```
# or
# ```math
# p(x) \equiv s(x) \pmod{\la g_1(x), \ldots, g_m(x) \ra}, s_0(x) \text{ is SOS}.
# ```
# the second one leads to a *simpler* SDP but needs to compute a *Gr\"obner* basis:
#  * SemialgebraicSets implements Buchberger's algorithm.
#  * The `@set` macro recognizes variable fix, e.g., `x = 1`
#     and provides shortcut.
#  * If you know a \alert{better} way to take modulo,
#    better create your \alert{own} type of algebraic set!

# We illustrate this in this example.

using DynamicPolynomials
@polyvar x[1:3]
p = sum(x)^2
using SumOfSquares
S = algebraicset([xi^2 - 1 for xi in x])

# We will now search for the minimum of `x` over `S` using Sum of Squares Programming.
# We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices.

import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)

function min_algebraic(S)
    model = SOSModel(solver)
    @variable(model, α)
    @objective(model, Max, α)
    @constraint(model, c, p >= α, domain = S)
    optimize!(model)
    @show termination_status(model)
    @show objective_value(model)
end

min_algebraic(S)

# Note that the minimum is in fact `1`.
# Indeed, since each variables is odd (it is either `-1` or `1`)
# and there is an odd number of variables, their sum is odd.
# Therefore it cannot be zero!

# We can see that the Gröbner basis of `S` was computed

@show S.I.gröbnerbasis
S.I.algo

# The Gröbner basis is simple to compute in this case as the vector
# of `xi^2 - 1` is already a Gröbner basis.
# However, we still need to divide polynomials by the Gröbner basis
# which can be simplified in this case.

const MP = MultivariatePolynomials
const SS = SemialgebraicSets
struct HypercubeIdeal{V} <: SS.AbstractPolynomialIdeal
    variables::Vector{V}
end
struct HypercubeSet{V} <: SS.AbstractAlgebraicSet
    ideal::HypercubeIdeal{V}
end
MP.variables(set::HypercubeSet) = MP.variables(set.ideal)
MP.variables(ideal::HypercubeIdeal) = ideal.variables
MP.changecoefficienttype(set::HypercubeSet, ::Type) = set
SS.ideal(set::HypercubeSet) = set.ideal
function Base.rem(p, set::HypercubeIdeal)
    return MP.polynomial(map(MP.terms(p)) do term
        mono = MP.monomial(term)
        new_mono = one(mono)
        for (var, exp) in powers(mono)
            if var in set.variables
                exp = rem(exp, 2)
            end
            new_mono *= var^exp
        end
        MP.coefficient(term) * new_mono
    end)
end

H = HypercubeSet(HypercubeIdeal(x))

min_algebraic(H)

# Let's now try to find the correct lower bound:

function min_algebraic_rational(S, d)
    model = SOSModel(solver)
    @variable(model, q, SOSPoly(MP.monomials(x, 0:d)))
    deno = q + 1
    @constraint(model, c, deno * p >= deno, domain = S)
    optimize!(model)
    @show termination_status(model)
end

# With `d = 0`, it's the same as previously

min_algebraic_rational(H, 0)

# But with `d = 1`, we can find the correct lower bound

min_algebraic_rational(H, 1)

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
atomicnames = Dict(
    1 => "Hydrogen",
    2 => "Helium",
    3 => "Lithium",
    4 => "Beryllium",
    5 => "Boron",
    6 => "Carbon",
    7 => "Nitrogen",
    8 => "Oxygen",
    9 => "Fluorine",
    10 => "Neon",
    11 => "Sodium",
    12 => "Magnesium",
    13 => "Aluminium",
    14 => "Silicon",
    15 => "Phosphorus",
    16 => "Sulfur",
    17 => "Chlorine",
    18 => "Argon",
    19 => "Potassium",
    20 => "Calcium",
    21 => "Scandium",
    22 => "Titanium",
    23 => "Vanadium",
    24 => "Chromium",
    25 => "Manganese",
    26 => "Iron",
    27 => "Cobalt",
    28 => "Nickel",
    29 => "Copper",
    30 => "Zinc",
    31 => "Gallium",
    32 => "Germanium",
    33 => "Arsenic",
    34 => "Selenium",
    35 => "Bromine",
    36 => "Krypton",
    37 => "Rubidium",
    38 => "Strontium",
    39 => "Yttrium",
    40 => "Zirconium",
    41 => "Niobium",
    42 => "Molybdenum",
    43 => "Technetium",
    44 => "Ruthenium",
    45 => "Rhodium",
    46 => "Palladium",
    47 => "Silver",
    48 => "Cadmium",
    49 => "Indium",
    50 => "Tin",
    51 => "Antimony",
    52 => "Tellurium",
    53 => "Iodine",
    54 => "Xenon",
    55 => "Cesium",
    56 => "Barium",
    57 => "Lanthanum",
    58 => "Cerium",
    59 => "Praseodymium",
    60 => "Neodymium",
    61 => "Promethium",
    62 => "Samarium",
    63 => "Europium",
    64 => "Gadolinium",
    65 => "Terbium",
    66 => "Dysprosium",
    67 => "Holmium",
    68 => "Erbium",
    69 => "Thulium",
    70 => "Ytterbium",
    71 => "Lutetium",
    72 => "Hafnium",
    73 => "Tantalum",
    74 => "Tungsten",
    75 => "Rhenium",
    76 => "Osmium",
    77 => "Iridium",
    78 => "Platinum",
    79 => "Gold",
    80 => "Mercury",
    81 => "Thallium",
    82 => "Lead",
    83 => "Bismuth",
    84 => "Polonium",
    85 => "Astatine",
    86 => "Radon",
    87 => "Francium",
    88 => "Radium",
    89 => "Actinium",
    90 => "Thorium",
    91 => "Protactinium",
    92 => "Uranium",
    93 => "Neptunium",
    94 => "Plutonium",
    95 => "Americium",
    96 => "Curium",
    97 => "Berkelium",
    98 => "Californium",
    99 => "Einsteinium",
    100 => "Fermium",
    101 => "Mendelevium",
    102 => "Nobelium",
    103 => "Lawrencium",
    104 => "Rutherfordium",
    105 => "Dubnium",
    106 => "Seaborgium",
    107 => "Bohrium",
    108 => "Hassium",
    109 => "Meitnerium",
    110 => "Darmstadtium",
    111 => "Roentgenium",
    112 => "Copernicium",
    113 => "Ununtrium",
    114 => "Ununquadium",
    115 => "Ununpentium",
    116 => "Ununhexium",
    117 => "Ununseptium",
    118 => "Ununoctium"
)

atomicsymbols = Dict(
    1 => "H",
    2 => "He",
    3 => "Li",
    4 => "Be",
    5 => "B",
    6 => "C",
    7 => "N",
    8 => "O",
    9 => "F",
    10 => "Ne",
    11 => "Na",
    12 => "Mg",
    13 => "Al",
    14 => "Si",
    15 => "P",
    16 => "S",
    17 => "Cl",
    18 => "Ar",
    19 => "K",
    20 => "Ca",
    21 => "Sc",
    22 => "Ti",
    23 => "V",
    24 => "Cr",
    25 => "Mn",
    26 => "Fe",
    27 => "Co",
    28 => "Ni",
    29 => "Cu",
    30 => "Zn",
    31 => "Ga",
    32 => "Ge",
    33 => "As",
    34 => "Se",
    35 => "Br",
    36 => "Kr",
    37 => "Rb",
    38 => "Sr",
    39 => "Y",
    40 => "Zr",
    41 => "Nb",
    42 => "Mo",
    43 => "Tc",
    44 => "Ru",
    45 => "Rh",
    46 => "Pd",
    47 => "Ag",
    48 => "Cd",
    49 => "In",
    50 => "Sn",
    51 => "Sb",
    52 => "Te",
    53 => "I",
    54 => "Xe",
    55 => "Cs",
    56 => "Ba",
    57 => "La",
    58 => "Ce",
    59 => "Pr",
    60 => "Nd",
    61 => "Pm",
    62 => "Sm",
    63 => "Eu",
    64 => "Gd",
    65 => "Tb",
    66 => "Dy",
    67 => "Ho",
    68 => "Er",
    69 => "Tm",
    70 => "Yb",
    71 => "Lu",
    72 => "Hf",
    73 => "Ta",
    74 => "W",
    75 => "Re",
    76 => "Os",
    77 => "Ir",
    78 => "Pt",
    79 => "Au",
    80 => "Hg",
    81 => "Tl",
    82 => "Pb",
    83 => "Bi",
    84 => "Po",
    85 => "At",
    86 => "Rn",
    87 => "Fr",
    88 => "Ra",
    89 => "Ac",
    90 => "Th",
    91 => "Pa",
    92 => "U",
    93 => "Np",
    94 => "Pu",
    95 => "Am",
    96 => "Cm",
    97 => "Bk",
    98 => "Cf",
    99 => "Es",
    100 => "Fm",
    101 => "Md",
    102 => "No",
    103 => "Lr",
    104 => "Rf",
    105 => "Db",
    106 => "Sg",
    107 => "Bh",
    108 => "Hs",
    109 => "Mt",
    110 => "Ds",
    111 => "Rg",
    112 => "Cp",
    113 => "Uut",
    114 => "Uuq",
    115 => "Uup",
    116 => "Uuh",
    117 => "Uus",
    118 => "Uuo"
)

atomicmasses = Dict(
    1 => 1.0079,
    2 => 4.0026,
    3 => 6.941,
    4 => 9.0122,
    5 => 10.811,
    6 => 12.0107,
    7 => 14.0067,
    8 => 15.9994,
    9 => 18.9984,
    10 => 20.1797,
    11 => 22.9897,
    12 => 24.305,
    13 => 26.9815,
    14 => 28.0855,
    15 => 30.9738,
    16 => 32.065,
    17 => 35.453,
    19 => 39.0983,
    18 => 39.948,
    20 => 40.078,
    21 => 44.9559,
    22 => 47.867,
    23 => 50.9415,
    24 => 51.9961,
    25 => 54.938,
    26 => 55.845,
    28 => 58.6934,
    27 => 58.9332,
    29 => 63.546,
    30 => 65.39,
    31 => 69.723,
    32 => 72.64,
    33 => 74.9216,
    34 => 78.96,
    35 => 79.904,
    36 => 83.8,
    37 => 85.4678,
    38 => 87.62,
    39 => 88.9059,
    40 => 91.224,
    41 => 92.9064,
    42 => 95.94,
    43 => 98,
    44 => 101.07,
    45 => 102.9055,
    46 => 106.42,
    47 => 107.8682,
    48 => 112.411,
    49 => 114.818,
    50 => 118.71,
    51 => 121.76,
    53 => 126.9045,
    52 => 127.6,
    54 => 131.293,
    55 => 132.9055,
    56 => 137.327,
    57 => 138.9055,
    58 => 140.116,
    59 => 140.9077,
    60 => 144.24,
    61 => 145,
    62 => 150.36,
    63 => 151.964,
    64 => 157.25,
    65 => 158.9253,
    66 => 162.5,
    67 => 164.9303,
    68 => 167.259,
    69 => 168.9342,
    70 => 173.04,
    71 => 174.967,
    72 => 178.49,
    73 => 180.9479,
    74 => 183.84,
    75 => 186.207,
    76 => 190.23,
    77 => 192.217,
    78 => 195.078,
    79 => 196.9665,
    80 => 200.59,
    81 => 204.3833,
    82 => 207.2,
    83 => 208.9804,
    84 => 209,
    85 => 210,
    86 => 222,
    87 => 223,
    88 => 226,
    89 => 227,
    91 => 231.0359,
    90 => 232.0381,
    93 => 237,
    92 => 238.0289,
    95 => 243,
    94 => 244,
    96 => 247,
    97 => 247,
    98 => 251,
    99 => 252,
    100 => 257,
    101 => 258,
    102 => 259,
    104 => 261,
    103 => 262,
    105 => 262,
    107 => 264,
    106 => 266,
    109 => 268,
    111 => 272,
    108 => 277
)

# tabulated nmax,lmax for a given atomic number Z

nlmax = Dict(
    0 => (1,0),
    1 => (1,0),
    2 => (2,0),
    3 => (2,0),
    4 => (2,1),
    5 => (2,1),
    6 => (2,1),
    7 => (2,1),
    8 => (2,1),
    9 => (2,1),
    10 => (2,1),
    11 => (3,1),
    12 => (3,1),
    13 => (3,2),
    14 => (3,2),
    15 => (3,2),
    16 => (3,2),
    17 => (3,2),
    18 => (3,2)
)

# log grid specified by tuple (xmin, xmax, num of points)
# dE = error compared to reference dftdata

atomicloggrid = Dict(
    1 => (-15.2,2.7,436),  # dE = 5e-6
    2 => (-17.2,2.3,475),  # dE = 3.8e-4
    3 => (-18.4,3.2,526),  # dE = 5.8e-4
    4 => (-19.4,2.9,543),  # dE = 7.4e-4
    5 => (-19.9,3.0,558),  # dE = 1e-3
    6 => (-20.8,2.75,717), # dE = 1.4e-3
    7 => (-21.2,2.65,726), # dE = 1.9e-3
    8 => (-21.6,2.55,735), # dE = 2.4e-3
    9 => (-21.6,2.45,732), # dE = 3e-3
    10 => (-23.25,2.35,654), # dE = 3.6e-3
    11 => (-24.0,3.1,2000),  # dE = 3.9e-3
    12 => (-23.75,2.95,1706),  # dE = 4.1e-3
    13 => (-24.0,3.1,1940),  # dE = 4.3e-3
    14 => (-24.25,2.95,2425),  # dE = 4.6e-3
    15 => (-24.25,2.85,2223),  # dE = 5e-3
    16 => (-24.75,2.75,1971),  # dE = 5.4e-3
    17 => (-24.75,2.65,1769),  # dE = 5.9e-3
    18 => (-24.75,2.55,1587)  # dE = 6.4e-3
)

# nao basis functions specified by tuple (n,l,Z)

nao_basis = Dict()
nao_basis[1] = ([(2,0,2.1),(2,1,3.5)],[(1,0,0.85),(2,1,3.7),(2,0,1.2),(3,2,7.0)])
nao_basis[2] = ([(1,0,0.85),(2,1,3.5)],)
nao_basis[3] = ([(2,0,2.),(3,2,2.6)],)
nao_basis[4] = ([(3,0,2.9),(3,2,3.5)],)
nao_basis[5] = ([(2,1,1.4),(3,2,4.8),(2,0,4.)],)
nao_basis[6] = ([(2,1,1.7),(3,2,6.),(2,0,4.9)],)
nao_basis[7] = ([(2,1,1.8),(3,2,6.8),(3,0,5.8)],)
nao_basis[8] = ([(2,1,1.8),(3,2,7.6),(3,0,6.4)],)
nao_basis[9] = ([(2,1,1.7),(3,2,7.4),(3,0,6.8)],)
nao_basis[10] = ([(3,2,6.),(2,1,1.5),(3,0,7.6)],)
nao_basis[11] = ([(2,1,1.2),(3,0,1.8),(3,2,3.8)],)
nao_basis[12] = ([(2,1,1.5),(3,0,2.4)],)
nao_basis[13] = ([],)
nao_basis[14] = ([(3,2,4.2),(2,1,1.4)],)
nao_basis[15] = ([],)
nao_basis[16] = ([(2,1,1.8)],)
nao_basis[17] = ([(2,1,1.9)],)
nao_basis[18] = ([(3,0,4.5)],)

# nao ionic basis functions specified by tuple (n,l) with configuration given below

nao_ionic = Dict()
nao_ionic[4] = ([(2,1)],)
nao_ionic[12] = ([(3,2)],)
nao_ionic[13] = ([(3,2),(3,1),(3,0)],)
nao_ionic[14] = ([(3,0)],)
nao_ionic[15] = ([(3,2),(3,1),(3,0)],)
nao_ionic[16] = ([(3,2),(3,0)],)
nao_ionic[17] = ([(3,2),(3,0)],)
nao_ionic[18] = ([(3,2),(4,1)],)

# nao ionic configuration specified by tuple of shell occupation numbers

nao_ionic_config = Dict()
nao_ionic_config[4] = [2,1]  # Be+, 1s2 2s1
nao_ionic_config[12] = [2,2,6]  # Mg(+2), 1s2 2s2 2p6
nao_ionic_config[13] = [2,2,6,1]  # Al(+2), 1s2 2s2 2p6 3s1
nao_ionic_config[14] = [2,2,6,1,1]  # Si(+2), 1s2 2s2 2p6 3s1 3p1
nao_ionic_config[15] = [2,2,6,1,2]  # P(+2), 1s2 2s2 2p6 3s1 3p2
nao_ionic_config[16] = [2,2,6,1,3]  # S(+2), 1s2 2s2 2p6 3s1 3p3
nao_ionic_config[17] = [2,2,6,1,4]  # Cl(+2), 1s2 2s2 2p6 3s1 3p4
nao_ionic_config[18] = [2,2,6,1,5]  # Ar(+2), 1s2 2s2 2p6 3s1 3p5

max_stars_repo_name: johncsnyder/Element


text:
const DEFAULT_INPUTS = 20:20:240
const DEFAULT_SECONDS = 1
const DEFAULT_SAMPLES = 200

struct RunReport
    names::Array{String,1}
    inputsizes::Array{Int, 1}
    times::Array{Real,2}
    gctimes::Array{Real,2}
    allocs::Array{Int,2}
    allocsizes::Array{Int,2}
end

"""
    RunReport(args...; kwargs...) -> Manager

Creates profiles for functions of runtime and memory consumption for different input sizes.

# Arguments

- `funcarray::Array{Base.Callable}` or `func::Base.Callable`: the function(s) to profile
- `genfunc::Base.Callable`: function that generates inputs of different sizes for
the function(s)
- `insizes`: iterable of integers that are used as inputs for `genfunc`

# Keywords

- `seconds::Number`: parameter of `@benchmarkable`
- `samples::Integer`: parameter of `@benchmarkable`
"""
function RunReport(
    funcarray::Array,
    genfunc::Base.Callable,
    insizes=DEFAULT_INPUTS;
    seconds=DEFAULT_SAMPLES,
    samples=DEFAULT_SAMPLES
    )
    funnames = String.(Symbol.(funcarray))
    times = Array{Real,2}(undef, size(insizes,1), size(funnames,1))
    gctimes = times |> copy
    allocs = Array{Int, 2}(undef, size(times)...)
    allocsizes = allocs |> copy

    for (i, n) in enumerate(insizes)
        orig_in = genfunc(n)
        for (j, func) in enumerate(funcarray)
            trial = run(@benchmarkable $func(x) setup=(x=copy($orig_in)) seconds=seconds samples=samples)
            times[i, j] = median(trial.times)
            gctimes[i, j] = median(trial.gctimes)
            allocsizes[i, j] = trial.memory
            allocs[i, j] = trial.allocs
        end
    end
    RunReport(funnames, insizes, times, gctimes, allocs, allocsizes)
end

function RunReport(
    func,
    genfunc::Base.Callable,
    insizes=DEFAULT_INPUTS;
    seconds=DEFAULT_SECONDS,
    samples=DEFAULT_SAMPLES)
    return RunReport([func], genfunc, insizes, seconds=seconds, samples=samples)
end

@recipe function f(report::RunReport)
    legend := :topleft
    label := reshape(report.names, 1, size(report.names, 1))
    layout := (2, 2)
    size --> (800, 600)

    for (i, (spy, sptitle)) in zip(
            [report.times, report.gctimes, report.allocs, report.allocsizes],
            ["Runtimes", "GC times", "Allocation counts", "Allocation sizes"]
        ) |> enumerate
        @series begin
            subplot := i
            title := sptitle
            (report.inputsizes, spy)
        end
    end
end


text:
function output_writer(u::EvolVars, chart2D::Chart, atlas, tinfo::Jecco.TimeInfo,
                       io::InOut, potential::Potential, phi0)
    Nsys = length(atlas)

    # output structures
    out_bdry  = Jecco.Output(io.out_dir, "boundary_", tinfo)
    out_gauge = Jecco.Output(io.out_dir, "gauge_", tinfo)
    out_bulk  = Jecco.Output(io.out_dir, "bulk_", tinfo)

    boundary  = getboundary(u)
    gauge     = getgauge(u)
    bulkevols = getbulkevolvedpartition(u)

    # NamedTuple with potential parameters
    params = parameters(potential)
    params = merge((phi0=phi0,), params)

    #=
    for analysis, it's useful to have the boundary behaviour of the bulk fields,
    so let's output those as well in the boundary file. note that the slicing
    operation [1,:,:] converts 3D arrays into 2D arrays; since we want the same
    shape as in the remaining boundary fields -- 3D arrays of size (1,Nx,Ny) --
    we must reshape them first.
    =#
    b14  = reshape(bulkevols[1].B1[1,:,:],  size(chart2D))
    b24  = reshape(bulkevols[1].B2[1,:,:],  size(chart2D))
    g4   = reshape(bulkevols[1].G[1,:,:],   size(chart2D))
    phi3 = reshape(bulkevols[1].phi[1,:,:], size(chart2D))

    # phi2 = phi0^3 phi3 - phi0 xi^2
    phi2 = phi0*phi0*phi0 .* phi3 .- phi0 .* gauge.xi .* gauge.xi

    # output fields
    boundary_fields = (
        Jecco.Field("a4",   boundary.a4,  chart2D),
        Jecco.Field("fx2",  boundary.fx2, chart2D),
        Jecco.Field("fy2",  boundary.fy2, chart2D),
        Jecco.Field("b14",  b14,          chart2D),
        Jecco.Field("b24",  b24,          chart2D),
        Jecco.Field("g4",   g4,           chart2D),
        Jecco.Field("phi2", phi2,         chart2D),
    )
    gauge_fields = Jecco.Field("xi", gauge.xi, chart2D)
    bulkevols_fields = ntuple(i -> (
        Jecco.Field("B1 c=$i",  bulkevols[i].B1,  atlas[i]),
        Jecco.Field("B2 c=$i",  bulkevols[i].B2,  atlas[i]),
        Jecco.Field("G c=$i",   bulkevols[i].G,   atlas[i]),
        Jecco.Field("phi c=$i", bulkevols[i].phi, atlas[i])
    ), Nsys)

    last_output_boundary_t = -io.out_boundary_every_t
    last_output_gauge_t    = -io.out_gauge_every_t
    last_output_bulk_t     = -io.out_bulk_every_t

    function (u::EvolVars)
        boundary  = getboundary(u)
        gauge     = getgauge(u)
        bulkevols = getbulkevolvedpartition(u)

        it = tinfo.it
        tt = tinfo.t

        do_output_boundary = false
        do_output_gauge    = false
        do_output_bulk     = false

        if (io.out_boundary_every > 0 && it % io.out_boundary_every == 0)
            do_output_boundary = true
        end
        if (io.out_gauge_every > 0 && it % io.out_gauge_every == 0)
            do_output_gauge = true
        end
        if (io.out_bulk_every > 0 && it % io.out_bulk_every == 0)
            do_output_bulk = true
        end

        if (io.out_boundary_every_t > 0 &&
            tt >= last_output_boundary_t + io.out_boundary_every_t - 1e-12)
            do_output_boundary     = true
            last_output_boundary_t = tt
        end
        if (io.out_gauge_every_t > 0 &&
            tt >= last_output_gauge_t + io.out_gauge_every_t - 1e-12)
            do_output_gauge     = true
            last_output_gauge_t = tt
        end
        if (io.out_bulk_every_t > 0 &&
            tt >= last_output_bulk_t + io.out_bulk_every_t - 1e-12)
            do_output_bulk     = true
            last_output_bulk_t = tt
        end

        if do_output_boundary
            boundary_fields[1].data = boundary.a4
            boundary_fields[2].data = boundary.fx2
            boundary_fields[3].data = boundary.fy2
            @views copyto!(boundary_fields[4].data, bulkevols[1].B1[1,:,:])
            @views copyto!(boundary_fields[5].data, bulkevols[1].B2[1,:,:])
            @views copyto!(boundary_fields[6].data, bulkevols[1].G[1,:,:])

            @views copyto!(phi2, bulkevols[1].phi[1,:,:])
            # phi2 = phi0^3 phi[1,:,:] - phi0 xi^2
            phi2 .= phi0*phi0*phi0 .* phi2 .- phi0 .* gauge.xi .* gauge.xi

            boundary_fields[7].data = phi2

            # write data
            out_bdry(boundary_fields, params=params)
        end

        if do_output_gauge
            gauge_fields.data = gauge.xi
            # write data
            out_gauge(gauge_fields, params=params)
        end

        if do_output_bulk
            @inbounds for i in 1:Nsys
                bulkevols_fields[i][1].data = bulkevols[i].B1
                bulkevols_fields[i][2].data = bulkevols[i].B2
                bulkevols_fields[i][3].data = bulkevols[i].G
                bulkevols_fields[i][4].data = bulkevols[i].phi
            end
            # write data
            out_bulk.(bulkevols_fields, params=params)
        end

        nothing
    end
end

function output_writer(bulkconstrains::BulkPartition{Nsys,BulkConstrained{T}}, atlas,
                       tinfo::Jecco.TimeInfo, io::InOut, potential::Potential,
                       phi0) where {Nsys,T}
    @assert Nsys == length(atlas)

    # NamedTuple with potential parameters
    params = parameters(potential)
    params = merge((phi0=phi0,), params)

    # output structure
    out  = Jecco.Output(io.out_dir, "constrained_", tinfo)

    # output fields
    fields = ntuple(i -> (
        Jecco.Field("S c=$i",    bulkconstrains[i].S,    atlas[i]),
        Jecco.Field("Fx c=$i",   bulkconstrains[i].Fx,   atlas[i]),
        Jecco.Field("Fy c=$i",   bulkconstrains[i].Fy,   atlas[i]),
        Jecco.Field("B1d c=$i",  bulkconstrains[i].B1d,  atlas[i]),
        Jecco.Field("B2d c=$i",  bulkconstrains[i].B2d,  atlas[i]),
        Jecco.Field("Gd c=$i",   bulkconstrains[i].Gd,   atlas[i]),
        Jecco.Field("phid c=$i", bulkconstrains[i].phid, atlas[i]),
        Jecco.Field("Sd c=$i",   bulkconstrains[i].Sd,   atlas[i]),
        Jecco.Field("A c=$i",    bulkconstrains[i].A,    atlas[i])
    ), Nsys)

    last_output_t = -io.out_bulkconstrained_every_t

    function (bulkconstrains)
        it = tinfo.it
        tt = tinfo.t

        do_output = false

        if (io.out_bulkconstrained_every > 0 && it % io.out_bulkconstrained_every == 0)
            do_output = true
        end

        if (io.out_bulkconstrained_every_t > 0 &&
            tt >= last_output_t + io.out_bulkconstrained_every_t - 1e-12)
            do_output = true
            last_output_t = tt
        end

        if do_output
            @inbounds for i in 1:Nsys
                fields[i][1].data = bulkconstrains[i].S
                fields[i][2].data = bulkconstrains[i].Fx
                fields[i][3].data = bulkconstrains[i].Fy
                fields[i][4].data = bulkconstrains[i].B1d
                fields[i][5].data = bulkconstrains[i].B2d
                fields[i][6].data = bulkconstrains[i].Gd
                fields[i][7].data = bulkconstrains[i].phid
                fields[i][8].data = bulkconstrains[i].Sd
                fields[i][9].data = bulkconstrains[i].A
            end
            # write data
            out.(fields, params=params)
        end

        nothing
    end
end

function checkpoint_writer(u::EvolVars, chart2D::Chart, atlas, tinfo::Jecco.TimeInfo,
                           io::InOut)
    Nsys = length(atlas)

    # output structure
    out  = Jecco.Output(io.checkpoint_dir, "checkpoint_it", tinfo)

    boundary  = getboundary(u)
    gauge     = getgauge(u)
    bulkevols = getbulkevolvedpartition(u)

    # output fields
    boundary_fields = (
        Jecco.Field("a4",   boundary.a4,  chart2D),
        Jecco.Field("fx2",  boundary.fx2, chart2D),
        Jecco.Field("fy2",  boundary.fy2, chart2D),
    )
    gauge_fields = Jecco.Field("xi", gauge.xi, chart2D)
    bulkevols_fields = ntuple(i -> (
        Jecco.Field("B1 c=$i",  bulkevols[i].B1,  atlas[i]),
        Jecco.Field("B2 c=$i",  bulkevols[i].B2,  atlas[i]),
        Jecco.Field("G c=$i",   bulkevols[i].G,   atlas[i]),
        Jecco.Field("phi c=$i", bulkevols[i].phi, atlas[i])
    ), Nsys)

    function (u::EvolVars)
        boundary  = getboundary(u)
        gauge     = getgauge(u)
        bulkevols = getbulkevolvedpartition(u)

        boundary_fields[1].data = boundary.a4
        boundary_fields[2].data = boundary.fx2
        boundary_fields[3].data = boundary.fy2

        gauge_fields.data = gauge.xi

        @inbounds for i in 1:Nsys
            bulkevols_fields[i][1].data = bulkevols[i].B1
            bulkevols_fields[i][2].data = bulkevols[i].B2
            bulkevols_fields[i][3].data = bulkevols[i].G
            bulkevols_fields[i][4].data = bulkevols[i].phi
        end

        # write data
        out(boundary_fields)
        out(gauge_fields)
        out.(bulkevols_fields)

        nothing
    end
end

# restore all bulk evolved fields
function restore!(bulkevols::BulkPartition{Nsys}, ts::OpenPMDTimeSeries,
                  it::Int) where {Nsys}
    for i in 1:Nsys
        B1, chart = get_field(ts, it=it, field="B1 c=$i")
        @assert size(B1) == size(bulkevols[i].B1)
        copyto!(bulkevols[i].B1, B1)

        B2, chart = get_field(ts, it=it, field="B2 c=$i")
        @assert size(B2) == size(bulkevols[i].B2)
        copyto!(bulkevols[i].B2, B2)

        G, chart = get_field(ts, it=it, field="G c=$i")
        @assert size(G) == size(bulkevols[i].G)
        copyto!(bulkevols[i].G, G)

        phi, chart = get_field(ts, it=it, field="phi c=$i")
        @assert size(phi) == size(bulkevols[i].phi)
        copyto!(bulkevols[i].phi, phi)
    end

    nothing
end

# restore all boundary fields
function restore!(boundary::Boundary, ts::OpenPMDTimeSeries, it::Int)

    a4, chart = get_field(ts, it=it, field="a4")
    @assert size(a4) == size(boundary.a4)
    copyto!(boundary.a4, a4)

    fx2, chart = get_field(ts, it=it, field="fx2")
    @assert size(fx2) == size(boundary.fx2)
    copyto!(boundary.fx2, fx2)

    fy2, chart = get_field(ts, it=it, field="fy2")
    @assert size(fy2) == size(boundary.fy2)
    copyto!(boundary.fy2, fy2)

    nothing
end

# restore gauge field
function restore!(gauge::Gauge, ts::OpenPMDTimeSeries, it::Int)
    xi, chart = get_field(ts, it=it, field="xi")
    @assert size(xi) == size(gauge.xi)
    copyto!(gauge.xi, xi)
    nothing
end

# recover all evolved variables from a checkpoint file
function recover(bulkevols::BulkPartition, boundary::Boundary, gauge::Gauge,
                 recovery_dir::String)
    ts = try
        OpenPMDTimeSeries(recovery_dir; prefix="checkpoint_it")
    catch e
        throw(e)
    end
    # grab last iteration of the timeseries, which should be the latest checkpoint
    it = ts.iterations[end]
    println("INFO: Recovering from checkpoint file $(ts.files[end])")

    restore!(bulkevols, ts, it)
    restore!(boundary, ts, it)
    restore!(gauge, ts, it)

    ts.current_iteration, ts.current_t
end

# for non-verbose output
vprint = x -> nothing

# for verbose output
# vprint(x) = println(x)


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "Point.jl" begin

	 # function norm(p::Point2d)
   @testset "norm" begin
      @test
      @test
      @test
      @test
   end

	 # function norm(p::Point3d)
   @testset "p::Point3d" begin
      @test
      @test
      @test
      @test
   end

	 # function norm(p::Point4d)
   @testset "norm" begin
      @test
      @test
      @test
      @test
   end

	 # function normalized(p::Point2d)
   @testset "normalized" begin
      @test
      @test
      @test
      @test
   end

	 # function normalized(p::Point3d)
   @testset "normalized" begin
      @test
      @test
      @test
      @test
   end

	 # function normalized(p::Point4d)
   @testset "normalized" begin
      @test
      @test
      @test
      @test
   end

	 # function dropW(p::Point4d)
   @testset "dropWdropW" begin
      @test
      @test
      @test
      @test
   end

end


text:
using CUDArt
import Base: isequal, convert, similar, copy, copy!, eltype, length, ndims, size, isempty, issparse, stride, strides, full
import CUDArt: to_host

# I want to make the base array explicit in the type signature of
# sparse arrays: So instead of SparseMatrixCSC{T,Int32} we use the
# equivalent KUsparse{Array,T}. Sparse arrays on the gpu can then be
# represented with KUsparse{CudaArray,T}.  The index type is always
# Int32, unlike SparseMatrixCSC, KUsparse does not support Int64
# indices.

type KUsparse{A,T}
    m                       # Number of rows
    n                       # Number of columns
    colptr                  # Column i is in colptr[i]:(colptr[i+1]-1)
    rowval                  # Row values of nonzeros
    nzval                   # Nonzero values
end

# construct KUsparse
KUsparse{A<:Array,T}(::Type{A}, ::Type{T}, m::Integer, n::Integer)=KUsparse{A,T}(m,n,ones(Int32,n+1),Array(Int32,0),Array(T,0))
KUsparse{A<:CudaArray,T}(::Type{A}, ::Type{T}, m::Integer, n::Integer)=KUsparse{A,T}(m,n,CudaArray(ones(Int32,n+1)),CudaArray(Int32,0),CudaArray(T,0))
KUsparse{A,T}(::Type{A}, ::Type{T}, d::NTuple{2,Int})=KUsparse(A,T,d...)
similar{A,T}(s::KUsparse{A}, ::Type{T}, m::Integer, n::Integer)=KUsparse(A,T,m,n)
similar{A,T}(s::KUsparse{A}, ::Type{T}, d::NTuple{2,Int})=KUsparse(A,T,d...)

isequal(a::KUsparse,b::KUsparse)=((typeof(a)==typeof(b)) && (sizeof(a)==sizeof(b)) && isequal(a.colptr,b.colptr) && isequal(a.rowval,b.rowval) && isequal(a.nzval,b.nzval))

# convert to KUsparse
convert{A<:Array,T}(::Type{KUsparse{A,T}}, s::SparseMatrixCSC{T,Int32})=KUsparse{Array,T}(s.m,s.n,s.colptr,s.rowval,s.nzval)
convert{A<:CudaArray,T}(::Type{KUsparse{A,T}}, s::SparseMatrixCSC{T,Int32})=KUsparse{CudaArray,T}(s.m,s.n,CudaArray(s.colptr),CudaArray(s.rowval),CudaArray(s.nzval))
convert{A<:BaseArray,T}(::Type{KUsparse{A,T}}, s::AbstractArray{T})=convert(KUsparse{A,T}, convert(SparseMatrixCSC{T,Int32}, s))
convert{A<:BaseArray,T}(::Type{KUsparse{A}}, s::AbstractArray{T})=convert(KUsparse{A,T}, convert(SparseMatrixCSC{T,Int32}, s))
convert{T}(::Type{KUsparse}, s::AbstractArray{T})=convert(KUsparse{Array,T}, s)
convert{T}(::Type{KUsparse}, s::AbstractCudaArray{T})=convert(KUsparse{CudaArray,T}, s)

# convert to SparseMatrixCSC
convert{A<:SparseMatrixCSC}(::Type{A}, s::KUsparse)=SparseMatrixCSC(s.m,s.n,convert(Array,s.colptr),convert(Array,s.rowval),convert(Array,s.nzval))
convert{A<:SparseMatrixCSC}(::Type{A}, a::CudaMatrix)=convert(SparseMatrixCSC, to_host(a))

# convert to Array/CudaArray
convert{A<:Array}(::Type{A}, s::SparseMatrixCSC)=full(s)
convert{A<:CudaArray}(::Type{A}, s::SparseMatrixCSC)=CudaArray(convert(Array, s))
convert{A<:Array}(::Type{A}, s::KUsparse)=convert(Array, convert(SparseMatrixCSC, s))
convert{A<:CudaArray}(::Type{A}, s::KUsparse)=convert(CudaArray, convert(SparseMatrixCSC, s))

# Now we can construct a KUsparse{CudaArray,T} using gpucopy:

copy(a::KUsparse)=deepcopy(a)
# copy{A<:BaseArray,T}(a::KUsparse{A,T})=KUsparse{A,T}(a.m,a.n,copy(a.colptr),copy(a.rowval),copy(a.nzval))

cpucopy_internal{A<:CudaArray,T}(s::KUsparse{A,T},d::ObjectIdDict)=
    (haskey(d,s)||(d[s] = KUsparse{Array,T}(s.m, s.n, cpucopy_internal(s.colptr,d), cpucopy_internal(s.rowval,d), cpucopy_internal(s.nzval,d))); d[s])

gpucopy_internal{A<:Array,T}(s::KUsparse{A,T},d::ObjectIdDict)=
    (haskey(d,s)||(d[s] = KUsparse{CudaArray,T}(s.m, s.n, gpucopy_internal(s.colptr,d), gpucopy_internal(s.rowval,d), gpucopy_internal(s.nzval,d))); d[s])

atype{A}(::KUsparse{A})=A
itype(::KUsparse)=Int32
eltype{A,T}(::KUsparse{A,T})=T
length(s::KUsparse)=(s.m*s.n)
ndims(::KUsparse)=2
size(s::KUsparse)=(s.m,s.n)
size(s::KUsparse,i)=(i==1?s.m:i==2?s.n:error("Bad dimension"))
strides(s::KUsparse)=(1,s.m)
stride(s::KUsparse,i)=(i==1?1:i==2?s.m:error("Bad dimension"))
isempty(s::KUsparse)=(length(s)==0)
to_host(s::KUsparse{CudaArray})=cpucopy(s)
issparse(::KUsparse)=true
# This won't work with CudaArray
# nnz(s::KUsparse)=s.colptr[s.n+1]-1  # TODO: shall we keep colptr in RAM?

atype(::SparseMatrixCSC)=Array
itype{T,I}(::SparseMatrixCSC{T,I})=I
full(s::KUsparse)=convert(KUdense, full(convert(SparseMatrixCSC, s)))
copy!{A,T}(a::KUsparse{A,T}, b::SparseMatrixCSC{T,Int64})=copy!(a, convert(SparseMatrixCSC{T,Int32},b))

function copy!{A,T}(a::KUsparse{A,T}, b::SparseMatrixCSC{T,Int32})
    a.m=b.m
    a.n=b.n
    for n in (:colptr, :rowval, :nzval)
        length(a.(n)) != length(b.(n)) && resize!(a.(n), length(b.(n)))
        copy!(a.(n), b.(n))
    end
    return a
end

# The final prediction output y should match the input x as closely as
# possible except for being dense.  These functions support obtaining
# the dense version of x.

dtype(x)=typeof(x)
dtype{T}(x::SparseMatrixCSC{T})=Array{T}
dtype{A,T}(x::KUsparse{A,T})=KUdense{A,T}
dtype{A,T}(x::KUdense{A,T})=KUdense{A,T}  # this is necessary, otherwise 1-dim x does not match 2-dim y.

dsimilar(x,d::Dims)=similar(x,d)
dsimilar{T}(x::SparseMatrixCSC{T},d::Dims)=Array(T,d)
dsimilar{A,T}(x::KUsparse{A,T},d::Dims)=KUdense(A,T,d)

function dsimilar!(l, n, x, dims=size(x))
    if (!isdefined(l,n) || !isa(l.(n), dtype(x)))
        l.(n) = dsimilar(x, dims)
    elseif (size(l.(n)) != dims)
        l.(n) = resize!(l.(n), dims)
    end
    return l.(n)
end

# type Sparse{A,T,I<:Integer}; m; n; colptr; rowval; nzval; end

# convert{T,I}(::Type{Sparse}, s::SparseMatrixCSC{T,I})=convert(Sparse{Array}, s)

# convert{A,T,I}(::Type{Sparse{A}}, s::SparseMatrixCSC{T,I})=
#     Sparse{A,T,I}(s.m,s.n,
#                   convert(A, s.colptr),
#                   convert(A, s.rowval),
#                   convert(A, s.nzval))

# convert{A,T,I}(::Type{SparseMatrixCSC}, s::Sparse{A,T,I})=
#     SparseMatrixCSC(s.m,s.n,
#                     convert(Array, s.colptr),
#                     convert(Array, s.rowval),
#                     convert(Array, s.nzval))

# convert{T}(::Type{Sparse}, a::Array{T,2})=convert(Sparse{Array}, a)
# convert{T}(::Type{Sparse}, a::CudaArray{T,2})=convert(Sparse{CudaArray}, a)

# convert{A<:Array,T}(::Type{Sparse{A}}, a::BaseArray{T,2})=
#     convert(Sparse, convert(SparseMatrixCSC{T,Int32}, sparse(convert(Array, a))))

# convert{A<:CudaArray,T}(::Type{Sparse{A}}, a::BaseArray{T,2})=
#     gpucopy(convert(Sparse{Array}, a))

# # Now we can construct a Sparse{CudaArray,T,I} using gpucopy:

# cpucopy_internal{A<:CudaArray,T,I}(s::Sparse{A,T,I},d::ObjectIdDict)=
#     (haskey(d,s) ? d[s] :
#      Sparse{Array,T,I}(s.m, s.n,
#                        cpucopy_internal(s.colptr,d),
#                        cpucopy_internal(s.rowval,d),
#                        cpucopy_internal(s.nzval,d)))

# gpucopy_internal{A<:Array,T,I}(s::Sparse{A,T,I},d::ObjectIdDict)=
#     (haskey(d,s) ? d[s] :
#      Sparse{CudaArray,T,I}(s.m,s.n,
#                            gpucopy_internal(s.colptr,d),
#                            gpucopy_internal(s.rowval,d),
#                            gpucopy_internal(s.nzval,d)))

# # And we can construct KUsparse which uses resizeable KUdense arrays for members:

# type KUsparse0{A,T,I<:Integer}
#     m::Int                      # Number of rows
#     n::Int                      # Number of columns
#     colptr::KUdense{A,I,1}      # Column i is in colptr[i]+1:colptr[i+1], note that this is 0 based on cusparse
#     rowval::KUdense{A,I,1}      # Row values of nonzeros
#     nzval::KUdense{A,T,1}       # Nonzero values
# end

# KUsparse0{A<:Array,T,I}(::Type{A}, ::Type{T}, ::Type{I}, m::Integer, n::Integer)=
#     KUsparse0{A,T,I}(m,n,KUdense(ones(I,n+1)),KUdense(Array(I,0)),KUdense(Array(T,0)))

# KUsparse0{A<:CudaArray,T,I}(::Type{A}, ::Type{T}, ::Type{I}, m::Integer, n::Integer)=
#     KUsparse0{A,T,I}(m,n,KUdense(CudaArray(ones(I,n+1))),KUdense(CudaArray(I,0)),KUdense(CudaArray(T,0)))

# KUsparse0{A,T,I}(::Type{A}, ::Type{T}, ::Type{I}, d::NTuple{2,Int})=KUsparse0(A,T,I,d...)

# similar{A,T,I}(s::KUsparse0{A}, ::Type{T}, ::Type{I}, m::Integer, n::Integer)=KUsparse0(A,T,I,m,n)

# similar{A,T,I,U}(s::KUsparse0{A,T,I}, ::Type{U}, d::NTuple{2,Int})=KUsparse0(A,U,I,d...)

# convert{A,T,I}(::Type{KUsparse0}, s::Sparse{A,T,I})=
#     KUsparse0{A,T,I}(s.m,s.n,convert(KUdense,s.colptr),convert(KUdense,s.rowval),convert(KUdense,s.nzval))

# convert{T,I}(::Type{KUsparse0}, s::SparseMatrixCSC{T,I})=
#     KUsparse0{Array,T,I}(s.m,s.n,convert(KUdense,s.colptr),convert(KUdense,s.rowval),convert(KUdense,s.nzval))

# convert{A,T,I}(::Type{Sparse}, s::KUsparse0{A,T,I})=
#     Sparse{A,T,I}(s.m,s.n,convert(A,s.colptr),convert(A,s.rowval),convert(A,s.nzval))

# convert{A,T,I}(::Type{SparseMatrixCSC}, s::KUsparse0{A,T,I})=
#     SparseMatrixCSC(s.m,s.n,convert(Array,s.colptr),convert(Array,s.rowval),convert(Array,s.nzval))

# convert{A<:KUsparse0{Array},T}(::Type{A}, a::BaseArray{T,2})=
#     convert(KUsparse0, convert(SparseMatrixCSC{T,Int32}, sparse(convert(Array, a))))

# convert{A<:KUsparse0{CudaArray},T}(::Type{A}, a::BaseArray{T,2})=
#     gpucopy(convert(KUsparse0{Array}, a))

# convert{A<:KUsparse0{CudaArray},T,I}(::Type{A}, s::SparseMatrixCSC{T,I})=
#     KUsparse0{CudaArray,T,I}(s.m,s.n,KUdense(CudaArray(s.colptr)),KUdense(CudaArray(s.rowval)),KUdense(CudaArray(s.nzval)))

# convert{A<:KUsparse0{Array},T,I}(::Type{A}, s::SparseMatrixCSC{T,I})=
#     KUsparse0{Array,T,I}(s.m,s.n,KUdense(s.colptr),KUdense(s.rowval),KUdense(s.nzval))

# ### BASIC COPY

# copy!{A,B,T,I}(a::KUsparse0{A,T,I}, b::KUsparse0{B,T,I})=
#     (a.m=b.m;a.n=b.n;for f in (:colptr,:rowval,:nzval); copy!(a.(f),b.(f)); end; a)

# copy!{A,T,I}(a::KUsparse0{A,T,I}, b::SparseMatrixCSC{T,I})=
#     (a.m=b.m;a.n=b.n;for f in (:colptr,:rowval,:nzval); copy!(a.(f),b.(f)); end; a)

# copy(a::KUsparse0)=KUsparse0(a.m,a.n,copy(a.colptr),copy(a.rowval),copy(a.nzval))

# # We need to fix cpu/gpu copy so the type changes appropriately:

# cpucopy_internal{A<:CudaArray,T,I}(s::KUsparse0{A,T,I},d::ObjectIdDict)=
#     (haskey(d,s) ? d[s] :
#      KUsparse0{Array,T,I}(s.m,s.n,
#                          cpucopy_internal(s.colptr, d),
#                          cpucopy_internal(s.rowval, d),
#                          cpucopy_internal(s.nzval, d)))

# gpucopy_internal{A<:Array,T,I}(s::KUsparse0{A,T,I},d::ObjectIdDict)=
#     (haskey(d,s) ? d[s] :
#      KUsparse0{CudaArray,T,I}(s.m,s.n,
#                              gpucopy_internal(s.colptr, d),
#                              gpucopy_internal(s.rowval, d),
#                              gpucopy_internal(s.nzval, d)))

# ### BASIC ARRAY OPS

# for S in (:KUsparse0, :Sparse)
#     @eval begin
#         atype{A}(::$S{A})=A
#         itype{A,T,I}(::$S{A,T,I})=I
#         eltype{A,T}(::$S{A,T})=T
#         length(s::$S)=(s.m*s.n)
#         ndims(::$S)=2
#         size(s::$S)=(s.m,s.n)
#         size(s::$S,i)=(i==1?s.m:i==2?s.n:error("Bad dimension"))
#         strides(s::$S)=(1,s.m)
#         stride(s::$S,i)=(i==1?1:i==2?s.m:error("Bad dimension"))
#         isempty(s::$S)=(length(s)==0)
#         to_host(s::$S{CudaArray})=cpucopy(s)
#         issparse(::$S)=true
#     end
# end

# atype(::SparseMatrixCSC)=Array
# itype{T,I}(::SparseMatrixCSC{T,I})=I
# full(s::KUsparse0)=convert(KUdense, full(convert(SparseMatrixCSC, s)))
# full{A}(s::Sparse{A})=convert(A, full(convert(SparseMatrixCSC, s)))

# convert{A<:CudaArray,T,I}(::Type{A}, s::SparseMatrixCSC{T,I})=CudaArray(full(s))
# convert{A<:SparseMatrixCSC,T,N}(::Type{A}, a::CudaArray{T,N})=sparse(to_host(a))

# convert{A<:Array,B,T,I}(::Type{A}, s::KUsparse0{B,T,I})=convert(Array, convert(SparseMatrixCSC, s))
# convert{A<:CudaArray,B,T,I}(::Type{A}, s::KUsparse0{B,T,I})=convert(CudaArray, convert(SparseMatrixCSC, s))

# # These two already defined in sparsematrix.jl:
# # convert(::Type{Matrix}, S::SparseMatrixCSC) = full(S)
# # sparse{Tv}(A::Matrix{Tv}) = convert(SparseMatrixCSC{Tv,Int}, A)
# # But they don't cover general array conversion:
# convert(::Type{Array}, S::SparseMatrixCSC) = full(S)


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
using Sparso
using CompilerTools
using CompilerTools.CFGs
using CompilerTools.OptFramework
using CompilerTools.LivenessAnalysis

function dump_liveness(func_ast :: Expr, func_arg_types :: Tuple, func_args)
    assert(typeof(func_ast)== Expr)
    assert(func_ast.head == :lambda)

    LivenessAnalysis.set_use_inplace_naming_convention()
    liveness = LivenessAnalysis.from_expr(func_ast)#, no_mod = Sparso.create_unmodified_args_dict())
    println("Liveness:\n", liveness)

    func_ast
end

sparse_pass = OptFramework.optPass(dump_liveness, true)
OptFramework.setOptPasses([sparse_pass])
#Sparso.set_debug_level(2)

include("./cg.jl")

A       = matrix_market_read(ARGS[1], true, true)
m       = size(A, 1)
x       = zeros(Float64, m)
b       = ones(Float64, m)
tol     = 1e-10
maxiter = 1000
@acc x, k, rel_err = cg(x, A, b, tol, maxiter)

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
module BGZF

using Zlib
import Zlib: z_stream, gz_header, Reader, libz, fillbuf, Z_OK
import Base: read, readbytes!, readuntil, close, eof, seek, position

export read, readbytes!, readuntil, close, eof, seek, position

type BGZFReader <: IO
    header::gz_header
    reader::Reader
    extra ::Array{Uint16}
end

function BGZFReader(io::IO; bufsize=65536)

    reader = Reader(io; bufsize=65536)

    header = gz_header()
    extra = Array(Uint16, 3)
    header.extra = pointer(extra)
    header.extra_max = 6

    fill_gzheader(header, reader, extra)

    BGZFReader(header, reader, extra)
end

function fill_gzheader(header::gz_header, reader::Reader, extra::Array{Uint16})

    # callable after inflateInit2() or inflateReset()
    ret = ccall((:inflateGetHeader, libz),
                Int32, (Ptr{z_stream}, Ptr{gz_header}),
                &reader.strm, &header)

    if ret != Z_OK
        error("BGZFReader: Error getting header")
    end

    # Really just want to read the header, but not yet implemented
    fillbuf(reader, 1)

    # Make sure this is a BGZF file

    if header.extra_len > 6
        error("Anticipated bug: this gzip compressed file contains more than one extra field, but we're only extracting the first.  Please file a bug report with the BGZF.jl package.")
    end

    extra8 = reinterpret(Uint8, extra)
    if header.extra_len < 6 ||
        extra8[1] != uint8(66) ||
        extra8[2] != uint8(67) ||
        extra[2] != 2
        error("BGZF extra field incorrect/missing; this does not seem to be a BGZF file")
    end
end

read{T}(r::BGZFReader, a::Array{T}) = read(r.reader, a)
read(r::BGZFReader, ::Type{Uint8}) = read(r.reader, Uint8)
readbytes!(r::BGZFReader, b::AbstractArray{Uint8}, nb=length(b)) = readbytes!(r.reader, b, nb)
readuntil(r::BGZFReader, delim::Uint8) = readuntil(r.reader, delim)
close(r::BGZFReader) = close(r.reader)
eof(r::BGZFReader) = eof(r.reader)
position(r::BGZFReader) = position(r.reader)

# seek using virtual offsets into file
function seek(r::BGZFReader, n::Integer)
    coffset = int64(n>>>16)
    uoffset = int64(n & 0xFFFF)
    seek(r, coffset, uoffset)
end

function seek(r::BGZFReader, coffset::Int64, uoffset::Int64)
    seek(r.reader.io, coffset)
    inflateReset(r.reader.strm)
    fill_gzheader(r.header, r.reader, r.extra)
    skip(r.reader, uoffset)
end

end # module


text:
using Test
using ParallelMapAccumulate
using Distributed

if nprocs() == 1
    addprocs(2)
elseif nworkers() > 2
    error("Too many processes")
end

@everywhere using Pkg
@everywhere pkg"activate ."
@everywhere pkg"precompile"
@everywhere using ParallelMapAccumulate

res = pmapaccumulate(x -> x^2, +, 1:10 |> collect; distributed = false) |> collect
@test res == accumulate(+, map(x -> x^2, 1:10))

res = pmapaccumulate(x -> x^2, +, 1:10) |> collect
@test res == accumulate(+, map(x -> x^2, 1:10))


text:
module LRUSets

import Base: push!, length, values, eltype, iterate

export LRUSet

struct LRUSet{K}
    max_size::Int
    map::Dict{K, UInt64}
    last::Threads.Atomic{UInt64}

    function LRUSet{K}(max_size::Int, map::Dict{K, UInt64}) where {K}
        last = length(map) > 0 ? maximum(values(map)) : UInt64(0)
        new(max_size, map, Threads.Atomic{UInt64}(last))
    end
end

function LRUSet{K}(max_size::Int) where {K}
    LRUSet{K}(
        max_size,
        Dict{K, UInt64}(),
    )
end

function LRUSet(max_size::Int, map::Dict{K, UInt64}) where {K}
    LRUSet{K}(
        max_size,
        map,
    )
end

LRUSet(max_size::Int) = LRUSet{Any}(max_size)

length(ts::LRUSet) = length(ts.map)
eltype(::Type{LRUSet{K}}) where {K} = K

values(ts::LRUSet) = keys(ts.map)

iterate(ts::LRUSet, state = 1) = iterate(keys(ts.map), state)

function push!(ts::LRUSet{K}, e::K) where {K}
    ts.map[e] = Threads.atomic_add!(ts.last, UInt64(1))
    if length(ts) > ts.max_size
        v, k = findmin(ts.map)
        delete!(ts.map, k)
    end
    nothing
end

end


text:
# Test Multiscale Conditional HINT network from Kruse et al. (2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test, Random
Random.seed!(2)

# Define network
nx = 16
ny = 16
n_in = 2
n_hidden = 4
batchsize = 2
L = 2
K = 2

function inv_test(nx, ny, n_in, batchsize, logdet, squeeze_type, split_scales)
    print("\nMultiscale Conditional HINT invertibility test with squeeze_type=$(squeeze_type), split_scales=$(split_scales), logdet=$(logdet)\n")
    CH = NetworkMultiScaleConditionalHINT(n_in, n_hidden, L, K; squeezer = squeeze_type(), logdet=logdet, split_scales=split_scales)

    # Input image and data
    X = randn(Float32, nx, ny, n_in, batchsize)
    Y = randn(Float32, nx, ny, n_in, batchsize)

    # Test inverse
    Zx, Zy = CH.forward(X, Y)[1:2]
    X_, Y_ = CH.inverse(Zx, Zy)

    @test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)
    @test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-3)

    # Test backward
    ΔZx = randn(Float32, size(Zx))  # random derivative
    ΔZy = randn(Float32, size(Zx))
    ΔX_, ΔY_, X_, Y_ = CH.backward(ΔZx, ΔZy, Zx, Zy)

    @test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)
    @test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-3)

    # Test inverse Y only
    Zy = CH.forward_Y(Y)
    Y_ = CH.inverse_Y(Zy)

    @test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-3)
end

# Loss
function loss(CH, X, Y)
    if CH.logdet
        Zx, Zy, logdet = CH.forward(X, Y)
        f = -log_likelihood(tensor_cat(Zx, Zy)) - logdet
    else
        Zx, Zy = CH.forward(X, Y)
        f = -log_likelihood(tensor_cat(Zx, Zy))
    end
    ΔZ       = -∇log_likelihood(tensor_cat(Zx, Zy))
    ΔZx, ΔZy = tensor_split(ΔZ)
    ΔX, ΔY   = CH.backward(ΔZx, ΔZy, Zx, Zy)[1:2]
    return f, ΔX, ΔY
end

function grad_test_X(nx, ny, n_channel, batchsize, logdet, squeeze_type, split_scales)
    print("\nMultiscale Conditional HINT gradient test with squeeze_type=$(squeeze_type), split_scales=$(split_scales), logdet=$(logdet)\n")
    CH = NetworkMultiScaleConditionalHINT(n_in, n_hidden, L, K; squeezer = squeeze_type(), logdet=logdet, split_scales=split_scales)

    # Input image
    X0 = randn(Float32, nx, ny, n_channel, batchsize)
    dX = randn(Float32, nx, ny, n_channel, batchsize)

    # Input data
    Y0 = randn(Float32, nx, ny, n_channel, batchsize)
    dY = 10*randn(Float32, nx, ny, n_channel, batchsize)

    f0, gX, gY = loss(CH, X0, Y0)[1:3]

    maxiter = 5
    h = 0.1f0
    err1 = zeros(Float32, maxiter)
    err2 = zeros(Float32, maxiter)

    for j=1:maxiter
        f = loss(CH, X0 + h*dX, Y0 + h*dY)[1]
        err1[j] = abs(f - f0)
        err2[j] = abs(f - f0 - h*dot(dX, gX) - h*dot(dY, gY))
        print(err1[j], "; ", err2[j], "\n")
        h = h/2f0
    end

    @test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)
end

for squeeze_i in [ShuffleLayer, WaveletLayer, HaarLayer]
    for split_scales in [true, false]
        for logdet in [false, true]
            inv_test(nx, ny, n_in, batchsize, logdet, squeeze_i, split_scales)
            grad_test_X(nx, ny, n_in, batchsize, logdet, squeeze_i, split_scales)
        end
    end
end

###################################################################################################
# Jacobian-related tests: NetworkMultiScaleConditionalHINT

# Gradient test

# Initialization
CH = NetworkMultiScaleConditionalHINT(n_in, n_hidden, L, K; split_scales=false, k1=3, k2=1, p1=1, p2=0, squeezer = ShuffleLayer());
CH.forward(randn(Float32, nx, ny, n_in, batchsize), randn(Float32, nx, ny, n_in, batchsize))
θ = deepcopy(get_params(CH))

CH0 = NetworkMultiScaleConditionalHINT(n_in, n_hidden, L, K; split_scales=false, k1=3, k2=1, p1=1, p2=0, squeezer = ShuffleLayer());

CH0.forward(randn(Float32, nx, ny, n_in, batchsize), randn(Float32, nx, ny, n_in, batchsize))
θ0 = deepcopy(get_params(CH0))

X = randn(Float32, nx, ny, n_in, batchsize)
Y = randn(Float32, nx, ny, n_in, batchsize)

# Perturbation (normalized)
dθ = θ-θ0; dθ .*= (norm.(θ).+1f-6)./(norm.(dθ).+1f-6)
dX = randn(Float32, nx, ny, n_in, batchsize); dX *= norm(X)/norm(dX)
dY = randn(Float32, nx, ny, n_in, batchsize); dY *= norm(Y)/norm(dY)

# Jacobian eval
dZx, dZy, Zx, Zy, _, _ = CH.jacobian(dX, dY, dθ, X, Y)

# Test
print("\nJacobian test\n")
h = 0.1f0
maxiter = 5
err5 = zeros(Float32, maxiter)
err6 = zeros(Float32, maxiter)
for j=1:maxiter
    set_params!(CH, θ+h*dθ)
    Zx_, Zy_, _ = CH.forward(X+h*dX, Y+h*dY)
    err5[j] = sqrt(norm(Zx_ - Zx)^2+norm(Zy_ - Zy)^2)
    err6[j] = sqrt(norm(Zx_ - Zx - h*dZx)^2+norm(Zy_ - Zy - h*dZy)^2)
    print(err5[j], "; ", err6[j], "\n")
    global h = h/2f0
end

@test isapprox(err5[end] / (err5[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err6[end] / (err6[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Adjoint test

set_params!(CH, θ)
dZx, dZy, Zx, Zy, _, _ = CH.jacobian(dX, dY, dθ, X, Y)
dZx_ = randn(Float32, size(dZx)); dZy_ = randn(Float32, size(dZy))
dX_, dY_, dθ_, _, _, _ = CH.adjointJacobian(dZx_, dZy_, Zx, Zy)
a = dot(dZx, dZx_)+dot(dZy, dZy_)
b = dot(dX, dX_)+dot(dY, dY_)+dot(dθ, dθ_)
@test isapprox(a, b; rtol=1f-1)  ####### need to check low accuracy here


text:
module FiberNlse
export nm, ns, ps, pm, km, mW, mm, GHz, THz, m, W, c, smf28,edfa, simulate, configure, inputSignal, Fiber, Bundle
using FFTW
using ProgressBars

include("datatypes.jl")
# Physical units & constants

nm = ns = 1e-9
ps = pm = 1e-12
km = kHz = 1e3
mW = mm = 1e-3
GHz = 1e9
Thz = 1e12
m = 1
W = 1
c = 299792458

function configure(Nₜ::Int,
    Nₗ::Int,
    fib::Fiber,
    T,
    λ, conf::SimulationConfig=default_config)

    return configure(Nₜ::Int,
        Nₗ::Int,
        fib.D,
        fib.γ,
        fib.α,
        fib.L,
        T,
        λ, conf)

end

function configure(
    Nₜ::Int,
    Nₗ::Int,
    D,
    γ,
    α,
    L,
    T,
    λ, conf::SimulationConfig=default_config)
    t = T * 0.5 * range(-1, stop=1, length=Nₜ) # Time vector
    l = L * range(0,stop=1,length=Nₗ)
    Ψ = Matrix{ComplexF32}(zeros((Nₗ, Nₜ)))
    return Simulation(Ψ, L / Nₗ, T / Nₜ, Nₜ, Nₗ, D, γ, α, L, T, λ, -D * λ^2 / (2 * pi * c), t, conf), t, l
end

function configure(
    Nₜ::Int,
    Nₗ::Int,
    fiber::Fiber,
    T,
    λ, conf::SimulationConfig=default_config)
    t = T * 0.5 * range(-1, stop=1, length=Nₜ) # Time vector
    l = fiber.L * range(0,stop=1,length=Nₗ)
    Ψ = Matrix{ComplexF32}(zeros((Nₗ, Nₜ)))
    return Simulation(Ψ, fiber.L / Nₗ, T / Nₜ, Nₜ, Nₗ, fiber.D, fiber.γ, fiber.α, fiber.L, T, λ, - fiber.D * λ^2 / (2 * pi * c), t, conf), t, l
end

function configure(
    fibs::Bundle,
    Nₜ::Int,
    T,
    λ, conf::SimulationConfig=default_config)

    t = T * 0.5 * range(-1, stop=1, length=Nₜ) # Time vector
    sims = Simulation[];

    l=Float64[];
    for (idx,fiber) in enumerate(fibs.elems)
        if typeof(fibs.Ns)==Int

                Nₗ = fibs.Ns;
        else

                Nₗ = fibs.Ns[idx];
        end

                li = fiber.L * range(0,stop=1,length=Nₗ)
                if isempty(l)
                    l=vcat(l,li)
                else
                    l=vcat(l,li.+l[end])
                end
                Ψ = Matrix{ComplexF32}(zeros((Nₗ, Nₜ)))
                sim = Simulation(Ψ, fiber.L / Nₗ, T / Nₜ, Nₜ, Nₗ, fiber.D, fiber.γ, fiber.α, fiber.L, T, λ, - fiber.D * λ^2 / (2 * pi * c), t, conf)
        push!(sims,sim)
    end

    return sims,t,l
end
function inputSignal(sim, ψₒ)
    sim.Ψ[1, :] = ψₒ
end

function edfa(L::Real, G::Real)
    g  = 0.1*log(10)*G/L
    fib = Fiber(-10*ps/nm/km, -g, 1/W/km, L);
    return fib
end

function smf28(L::Real)
    fib = Fiber(-17*ps/nm/km, 0.046/km, 1.1/W/km, L);
    return fib
end

function simulate(sim::Simulation, progress::Bool=false)

    D(sim, ν) = @. -sim.β2 * 0.5im * (2 * pi * ν)^2
    N̂(sim, u) = @. abs(u)^2 * sim.γ * 1im

    # Check wether to show progressbar or not
    if progress
        iter_z = ProgressBar(range(2, sim.Nₗ))
    else
        iter_z = range(2, sim.Nₗ)
    end

    ν = FFTW.fftfreq(sim.Nₜ, 1.0 / (sim.dt))

    dl = sim.dz
    α = sim.α

    D̂ = D(sim, ν)

    for i in iter_z

        ψₗ = sim.Ψ[i-1, :]

        ψₗ = ifft(exp.(0.5 * dl .* (D̂ .- 0.5α)) .* fft(ψₗ))
        ψₗ = exp.(dl * N̂(sim, ψₗ)) .* ψₗ
        ψₗ = ifft(exp.(0.5 * dl .* (D̂ .- 0.5α)) .* fft(ψₗ))

        sim.Ψ[i, :] = ψₗ

    end

end

function simulate(sims::Vector{Simulation}, Ψₒ::Vector{Float64},progress::Bool=false)
    S = []
    for (idx, sim) in enumerate(sims)
        if idx == 1
            inputSignal(sim,Ψₒ)
            simulate(sim,progress)
            S = sim.Ψ
        else
            inputSignal(sim, S[end, :])
            simulate(sim,progress)
            S = vcat(S, sim.Ψ)
        end

    end
    return S
end

function transition(sim1::Simulation, sim2::Simulation) # Use the final state of sim1 as initial state of sim2
    inputSignal(sim2, sim1.Ψ[end, :])
end

function rhs()
    # ! Lire le code de gnlse-python (les mettre dans les crédits)
end

function simulate2(sim::Simulation)
    # ! Integration RK4 de l'ode rhs
end
end


text:
using Distributions, StatsBase
using PyPlot

function qqplot(z; color="black")
    n = length(z)

    grid = [(1 / (n + 1)):(1 / (n + 1)):(1.0 - (1 / (n + 1)));]

    qz = quantile(z, grid)
    qd = quantile.(Ref(Distributions.Normal()), grid)

    lims = 3.290
    x = range(-lims, stop=lims)

    plot(x, x, color="grey", linestyle=":", linewidth=.25)
    scatter(qz, qd, s=.75, color=color)

    xlim([-lims, lims])
    ylim([-lims, lims])

    nothing
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using BenchmarkTools
using CUDA
using Oceananigans
using Benchmarks
using Plots
pyplot()
# Benchmark function

function benchmark_nonhydrostatic_model(Arch, FT, N)
    grid = RectilinearGrid(Arch(), FT, size=(N, N, N), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid=grid)

    time_step!(model, 1) # warmup

    trial = @benchmark begin
        @sync_gpu time_step!($model, 1)
    end samples=10

    return trial
end

function benchmark_hydrostatic_model(Arch, FT, N)
    grid = RectilinearGrid(Arch(), FT, size=(N, N, 10), extent=(1, 1, 1))
    model = HydrostaticFreeSurfaceModel(grid=grid,
                                        tracers = (),
                                        buoyancy = nothing,
                                        free_surface=ImplicitFreeSurface())

    time_step!(model, 0.001) # warmup

    trial = @benchmark begin
        @sync_gpu time_step!($model, 0.001)
    end samples=10

    return trial
end

function benchmark_shallowwater_model(Arch, FT, N)
    grid = RectilinearGrid(Arch(), FT, size=(N, N), extent=(1, 1), topology = (Periodic, Periodic, Flat))
    model = ShallowWaterModel(grid=grid, gravitational_acceleration=1.0)

    time_step!(model, 1) # warmup

    trial = @benchmark begin
        @sync_gpu time_step!($model, 1)
    end samples=10

    return trial
end

# Benchmark parameters

Architectures = has_cuda() ? [CPU, GPU] : [CPU]
Float_types = [Float64]
Ns = [32, 64, 128, 256]

# Run and summarize benchmarks

print_system_info()

for (model, name) in zip((:nonhydrostatic, :hydrostatic, :shallowwater), ("NonhydrostaticModel", "HydrostaticFreeSurfaceModel", "ShallowWaterModel"))

    benchmark_func = Symbol(:benchmark_, model, :_model)
    @eval begin
        suite = run_benchmarks($benchmark_func; Architectures, Float_types, Ns)
    end

    df = benchmarks_dataframe(suite)
    benchmarks_pretty_table(df, title=name * " benchmarks")
end


text: @deprecate NelderMead(a::Real, g::Real, b::Real) NelderMead(initial_simplex = AffineSimplexer(), parameters = FixedParameters(a, g, b, 0.5))


text:
function normalize_factor{R <: Real}(arr::Array{R})
    arr ./ sum(arr)
end

function normpdf(mu, sigma)
    x -> (1/(sigma*sqrt(2*pi)))*e^((-(x-mu)^2)/(2*sigma^2))
end

function group_combinations(groups::Array{Array,1})
    if length(groups) > 0
        elem_type=eltype(groups[1])
        combinations=prod(map(length, groups))
        res=Array(elem_type, combinations, length(groups))

        groups=groups[sortperm(map(length, groups))]

        prev_group_elem_size=combinations

        for i=1:length(groups)

            g=groups[i]
            prev_steps=convert(Int,combinations/prev_group_elem_size)
            group_elem_size=convert(Int,prev_group_elem_size/length(g))

            step_elems=reduce(vcat,map(x -> fill(x,group_elem_size), g))
            steps=reduce(vcat,fill(step_elems,prev_steps))

            res[:,i]=steps
            prev_group_elem_size=group_elem_size
        end

        sets=Array(Set{elem_type},combinations)
        sets[:] = [ Set(res[i,:]) for i=1:combinations ]
    else
        []
    end
end

function single_to_multi_index(i::Integer, dims::(Integer,Integer))
    dim1, dim2 = dims
    convert(Integer,ceil(i/dim1)), i%dim1 > 0 ? i%dim1 : dim1
end

order_perm(A1::Array, A2::Array) = map(x -> findfirst(A1,x), A2)


text:
using Base: front, tail

"""
    ⊗(t1::Tensor{T}, t2::Tensor{S}) where {T,S} -> Tensor{promote_type(T, S)}

Computes the outer product of the tensors `t1` and `t2`, also referred to as the tensor
product of `t1` and `t2`.
Arguments can be of arbitrary size and dimension.

# Arguments
- `t1::Tensor{T}`: the first factor
- `t2::Tensor{S}`: the second factor

# Returns
- `Tensor{promote_type(T, S)}`: the outer product of `t1` and `t2`
"""
function ⊗(t1::Tensor{T}, t2::Tensor{S}) where {T,S}
    res = Array{promote_type(T, S)}(undef, size(t1)..., size(t2)...)

    j, k = 1, 1
    for i in eachindex(res)
        @inbounds res[i] = t1.data[j] * t2.data[k]
        j = j % length(t1) + 1
        k += j == 1 ? 1 : 0
    end

    return Tensor(res)
end

"""
    ⊗(ts...) -> Tensor{promote_type(eltype.(ts)...)}

Computes the outer product of an arbitrary number of tensors `ts...`. The empty outer
product is defined as `Tensor(1)` (i.e., the neutral element of the outer product).
Arguments can be of arbitrary size and dimension.

# Arguments
- `ts...`: the (possibly empty) list of tensors to be multiplied

# Returns
- `Tensor{promote_type(eltype.(ts)...)}`: the outer product of the tensors `ts...`
"""
⊗(ts...) = reduce(⊗, ts; init=Tensor(1))

"""
    ⋅(t1::Tensor{T}, t2::Tensor{S}) where {T,S} -> Tensor{promote_type(T, S)}

Computes the inner product of the tensors `t1` and `t2`, also referred to as the scalar
product of `t1` and `t2`.
Arguments need to be of non-zero dimension and the last dimension of `t1` needs to match the
first dimension of `t2`.

# Arguments
- `t1::Tensor{T}`: the first factor
- `t2::Tensor{S}`: the second factor

# Returns
- `Tensor{promote_type(T, S)}`: the inner product of `t1` and `t2`

# Throws
- `ArgumentError`: if either `t1` or `t2` is 0-dimensional
- `DimensionMismatch`: if the last dimension of `t1` does not match the first dimension of
    `t2`
"""
function ⋅(t1::Tensor{T}, t2::Tensor{S}) where {T,S}
    @argcheck ndims(t1) > 0
    @argcheck ndims(t2) > 0
    @argcheck last(size(t1)) == first(size(t2)) DimensionMismatch

    return pushover(conj(t1), t2, ndims(t1), 1)
end

"""
    contract(t::Tensor, d1::Integer, d2::Integer) -> Tensor

Computes the contraction of the tensor `t` along the dimensions `d1` and `d2`.
`d1` and `d2` need to be distinct, valid indices for the dimensions of `t`. Additionally,
the `d1`-th dimension of `t` has to match its `d2`-th dimension.

# Arguments
- `t::Tensor`: the tensor to be contracted
- `d1::Integer`: the first contraction dimension
- `d2::Integer`: the second contraction dimension

# Returns
- `Tensor`: the contracted tensor

# Throws
- `ArgumentError`: if `d1` and `d2` are not distinct
- `BoundsError`: if `d1` and `d2` do not index valid dimensions of `t`
- `DimensionMismatch`: if the `d1`-th dimension does not match the `d2`-th dimension of `t`
"""
function contract(t::Tensor, d1::Integer, d2::Integer)
    @argcheck d1 != d2
    @argcheck 1 <= d1 <= ndims(t) BoundsError(size(t), d1)
    @argcheck 1 <= d2 <= ndims(t) BoundsError(size(t), d2)
    @argcheck size(t, d1) == size(t, d2) DimensionMismatch

    # allow more efficient access of the needed (column-major) array elements
    a = permutedims(t.data, (d1, d2, Tuple(k for k in 1:ndims(t) if k != d1 && k != d2)...))

    res = Array{eltype(t)}(undef, size(a)[3:ndims(a)])
    inds = CartesianIndices(res)

    # NOTE: replacing the cartesian indexing by on-the-fly linear indexing hurts performance
    for i in eachindex(res)
        @inbounds res[i] = sum(a[j, j, inds[i]] for j in axes(a, 1))
    end

    return Tensor(res)
end

"""
    trace(t::Tensor{T,2}) where {T} -> Tensor{T,0}

Computes the contraction of the tensor `t` along its two dimensions, equivalent to the trace
of the underlying matrix of `t`.
The two dimensions `t` need to match, i.e., `t` needs to be a 'square' tensor.

# Arguments
- `t::Tensor{T,2}`: the tensor to be contracted

# Returns
- `Tensor{T,0}`: the contracted scalar tensor

# Throws
- `DimensionMismatch`: if the first and the second dimension of `t` do not match
"""
trace(t::Tensor{T,2}) where {T} = contract(t, 1, 2)

"""
    pushover(t1::Tensor{T}, t2::Tensor{S}, d1::Integer, d2::Integer) where {T,S} -> Tensor{promote_type(T, S)}

Computes the 'Überschiebung' of the tensors `t1` and `t2` along the dimensions `d1` and `d2`
respectively.
`d1` and `d2` need to be valid indices for the dimensions of `t1` and `t2` respectively.
Additionally, the `d1`-th dimension of `t1` has to match the `d2`-th dimension of `t2`.

# Arguments
- `t1::Tensor{T}`: the first tensor to push over
- `t2::Tensor{S}`: the second tensor to push over
- `d1::Integer`: the first push-over dimension
- `d2::Integer`: the second push-over dimension

# Returns
- `Tensor{promote_type(T, S)}`: the push-over tensor

# Throws
- `BoundsError`: if `d1` and `d2` do not index valid dimensions of `t1` and `t2`
    respectively
- `DimensionMismatch`: if the `d1`-th dimension of `t1` does not match the `d2`-th dimension
    of `t2`
"""
function pushover(t1::Tensor{T}, t2::Tensor{S}, d1::Integer, d2::Integer) where {T,S}
    @argcheck 1 <= d1 <= ndims(t1) BoundsError(size(t1), d1)
    @argcheck 1 <= d2 <= ndims(t2) BoundsError(size(t2), d2)
    @argcheck size(t1, d1) == size(t2, d2) DimensionMismatch

    # allow more efficient access of the needed (column-major) array elements
    a1 = permutedims(t1.data, (d1, Tuple(k for k in 1:ndims(t1) if k != d1)...))
    a2 = permutedims(t2.data, (d2, Tuple(k for k in 1:ndims(t2) if k != d2)...))

    res = Array{promote_type(T, S)}(
        undef, size(a1)[2:ndims(a1)]..., size(a2)[2:ndims(a2)]...
    )

    # NOTE: improved indexing, but still not great, more optimization?
    j, k = 1, 1
    for i in eachindex(res)
        @inbounds res[i] = sum(a1[j:(j + size(a1, 1) - 1)] .* a2[k:(k + size(a1, 1) - 1)])
        j = (j + size(a1, 1)) % length(t1)
        k += j == 1 ? size(a1, 1) : 0
    end

    return Tensor(res)
end

# TODO: Implement epsilon tensor
# TODO: Implement delta tensor
# TODO: Think about sparse tensors for efficiency?
# TODO: Implement cross product (how?)
# TODO: Add examples to all relevant docstrings
# TODO: Maybe rename actual implementations of ⊗ and ⋅ to outer and inner and alias symbols