text:
#=2013 LLVM Developers’ Meeting: “Julia: An LLVM-based approach to scientific computing”
by Mr. Keno Fischer - Harvard College/MIT CSAIL "Computer Scientist"

slides
https://llvm.org/devmtg/2013-11/slides/Fischer-Julia.html#/7

Video: https://www.youtube.com/watch?v=UsdQcbpVWdM&ab_channel=LLVM
2013 -
Benchmarks = readtable("http")

=#

#=Agenda

- Part1: Why LLVM
- Part2:
    2.1. introduction
    2.2 High-level Overview of how we go from JuliaCode to LLVM
    2.3 How we use LLVM as a Part of Julia
=#

#=Is It any fast?

    When introducing a scientific COmputing Language,
    the first Question is:
        Q. Is it any fast?
       #benchmark with a Log-Scale plot
       deleted #REASON INVALID URI

    -short answer: Yes, we Are , however, being fast Is Not What We All About !

    -being fast:
    1. Emerges from the way we DO Things

    2. Something you need: for Scientific Computing
    But, if we All wanted to do is (to) Be Fast,
    We wouldn't need a new Programming Language

    3. It Is Just a Side Effect   that we will pass , but:
        "It's not the main thing, that we're about "

    with that, let's get into a Brief Overview of what Julia is About

     =#

# the equivalent of Hello world
username = run(`git config --get user.name`)

#=Dynamic Multiple Dispatch Everywhere

    (one of the) the main Paradigms in Julia is
    this Idea "Multiple Dispatch":

    Definition: every function depends on ALL the Types
    of the Arguments you pass in

##functions

1.Convert
 - probably (was) the most overloaded function , now it's got only 10 methods
 - it has in the base system, just withouot loading

 2. the Plus +

 -Highly Polymorphic  with 196 methods

-for any argument (to a function):
    write the call & compiler will tell which method is selected

    if we want to convert 1 as int64 to Float64
        it would choose the (proper) convert method

    the plus + (between float & int )

    - is a more generic function

    - that Dispatches based on the number Type
        (which is a common SuperType of both of them  )
        that's called promotion()   - later
    =#

convert  # 10 method generic function

+ #196 method generic function (still the largest)

@which convert(Float64, 1) #no unique method found #checking uniqueness of method (if method has n= 1 passed with it )

@which 1.0 + 2 # can add 2 numbers of different types i.e +( float, Int64)
#=+(x::Number, y::Number) via promotion

=#

#= Basic Dispatch

the first function, its first argument: Anyy, unconstrained (can be anything ) -accepts different variables & types

the second function has with both Arguments of Abstract type Number
            we have `Type Heirarchy`
      we do `Type Annotations` by using  ::Any
        A return type is `implicit` (when no type is specified)
        but, if want to have more specific, use  Annotations ::

       In Julia,  we don't have (Structural) Inheritance,an idea of sth we may do (in the future) -haven't missed it (so far that much)

       in Julia, we  define a Type Heirarchy, that  we have
       1.1. Abstract Type :  like Number & Integers (inherits from Number), which everything inherits from it (i.e. floating points Floats, Ints (Int64,32,..),Rationals
       1.2. Concrete Type: which are leaf types (i.e implementable), &  have Fields (for specidic application)

=#
#--- Basic Dispatch

f(a::Any, b) = "fallback" # a generic fallback
f(a::Number, b::Number) = "a and b are both numbers" # both of Type Number
f(a::Number, b) = "a is a number" # second argument isa ::Any
f(a, b::Number) = "b is a number"
f(a::Integer, b::Integer) = "a and b are both integers" #we could be more specific (constrained)

f(1, 5.2) # "a and b are both numbers"
f(1, "bar") #= "a is a number" =#
f(1, 2) #= "a and b are both integers" =#

f("foo", [1, 2]) # : "fallback" -- it compares all functions, none has array explicitly specified, hence  it picks  the `fallback``  (because its 2nd  argument is dynamic ) # Infer 1st: f(a::Any, b) = "fallback"

methods(f) #returns all different methods (& its different types function has ) #convenient
#=julia> methods(f)
# 5 methods for generic function "f":
[1] f(a::Integer, b::Integer) in Main at c:\Users\adamus\.git\DeepLearner\DeepLearner\Tutorial\LLVM\llvm.jl:32
[2] f(a::Number, b::Number) in Main at c:\Users\adamus\.git\DeepLearner\DeepLearner\Tutorial\LLVM\llvm.jl:29
[3] f(a, b::Number) in Main at c:\Users\adamus\.git\DeepLearner\DeepLearner\Tutorial\LLVM\llvm.jl:31[4] f(a::Number, b) in Main at c:\Users\adamus\.git\DeepLearner\DeepLearner\Tutorial\LLVM\llvm.jl:30[5] f(a, b) in Main at c:\Users\adamus\.git\DeepLearner\DeepLearner\Tutorial\LLVM\llvm.jl:28
=#
#=--- Diagonal Dispatch/ Parametric Dispatch =#

f{T<:Number}(a::T, b::T) = "a and b are both $(T)s"
f{T<:Number}(a::Vector{T}, b::Vector{T}) = "a and b are both vectors of $(T)s"

#=Crossreference  Dispatch on tuples of types https://github.com/JuliaLang/julia/issues/10947
#10947
(Int,) is not a type at all; it's a tuple that contains a type. So you can't use <: on it.

Tuples now "infer" their parameters from constructor arguments just like other parametric types do:
=#
typeof(Int => String) #Pair{DataType,DataType}
typeof((Int, String)) # Tuple{DataType,DataType}
[Int] #1-element Array{DataType,1}: #Int64 # Core.Int64

#=
The difference is that tuples don't provide a way to request parameter types. If you want, you can construct Pair{Any,Any}(1,2), but you can't do the same thing with tuples. Tuple "parameters" are always inferred.

This is normally what you'd expect, but the corner case is that it's impossible to construct a Tuple{Type{Int}}.

To do this dispatch, I think you'll have to pass tuple types --- dispatch on Type{Tuple{T,T}}. Unfortunately the input syntax Tuple{A,B} is much less pleasant than (A,B). Another example in favor of using {A,B} for tuple types.

The difference is that tuples don't provide a way to request parameter types. If you want, you can construct Pair{Any,Any}(1,2), but you can't do the same thing with tuples. Tuple "parameters" are always inferred.

This is normally what you'd expect, but the corner case is that it's impossible to construct a Tuple{Type{Int}}.

To do this dispatch, I think you'll have to pass tuple types --- dispatch on Type{Tuple{T,T}}. Unfortunately the input syntax Tuple{A,B} is much less pleasant than (A,B)
Another example in favor of using {A,B} for tuple types.
=#

isa(Int, Any)
#true
isa(Int, Type{Int})
#true
isa(Int, DataType) #true
isa(("",), Tuple{Any}) #true
isa(("",), Tuple{String}) #true
isa(("",), Tuple{ASCIIString}) #true #2022: Error UndefVarErro (is this normal?)
isa((Int,), Tuple{Any}) #true
isa((Int,), Tuple{Type{Int}}) # false
isa((Int,), Tuple{DataType}) #true
#Edit: Actually, upon further reflection, I understand now. typeof((Int,)) is a Tuple{DataType}. And this is behaving like dispatch:

f(::Tuple{Type{Int}}) = 1 #f (generic function with 1 method) #works

f((Int,)) #ERROR: MethodError: `f` has no method matching f(::Tuple{DataType}) # Closest candidates are: f(::Tuple{Type{Int64}}) #me : so the ....
#ERROR: MethodError: `f` has no method matching f(::Tuple{DataType}) Closest candidates are: f(::Tuple{Type{Int64}})

#= Warning f(::Tuple{Type{Int64}}) #f(::Tuple{DataType}) #OPEN! #BUGGY ERROR Detected (for more review: src\UI\Gtk\issues\open\isa_issue.jl for more examples of buggy code)
#=Remark:in 2013, pointing to a hidden bug, with a one word f((Int,)) -that breaks Julia itself- was later discovered by chance in 2015; till today the issue is open (because it is as deep as LLVM)
the presenter is the least, but the fellow 'scientists' that he resemble are, hereby, spooky=#

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


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

# Wrapping

function ansi_length(s)
    replace(s, r"\e\[[0-9]+m" => "") |> length
end

words(s) = split(s, " ")
lines(s) = split(s, "\n")

function wrapped_lines!(lines, io::IO, s::AbstractString, width, i)
    ws = words(s)
    for word in ws
        word_length = ansi_length(word)
        word_length == 0 && continue
        if isempty(lines) || i + word_length + 1 > width
            i = word_length
            push!(lines, word)
        else
            i += word_length + 1
            lines[end] *= " " * word   # this could be more efficient
        end
    end
    return i
end

function wrapped_lines(io::IO, s::AbstractString; width = 80, i = 0)
    lines = AbstractString[]
    if occursin(r"\n", s)
        for ss in split(s, "\n")
            i = wrapped_lines!(lines, io, ss, width, i)
        end
    else
        wrapped_lines!(lines, io, s, width, i)
    end
    return lines
end

wrapped_lines(io::IO, f::Function, args...; width = 80, i = 0) =
    wrapped_lines(io, sprint(f, args...; context=io), width = width, i = 0)

function print_wrapped(io::IO, s...; width = 80, pre = "", i = 0)
    lines = wrapped_lines(io, s..., width = width, i = i)
    isempty(lines) && return 0, 0
    print(io, lines[1])
    for line in lines[2:end]
        print(io, '\n', pre, line)
    end
    length(lines), length(pre) + ansi_length(lines[end])
end

print_wrapped(f::Function, io::IO, args...; kws...) = print_wrapped(io, f, args...; kws...)

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
###################################################
# METODO DE PONTOS INTERIORES PRIMAL DUAL AFIM
# ESCALA PARA PL's NA FORMA PADRÃO
#
# min c^t x
# s.a. Ax = b
#       x >= 0
#
# Autor: Leonardo D. Secchin
# Data : 24/06/2021
# leonardosecchin.github.io/topicospo/
#
# Entrada:
#    -- Obrigatórios --
#    A,b,c:   dados do problema (esparsos ou não)
#
#    -- Opcionais --
#    maxit : número máximo de iterações (padrão = 100)
#    eps   : precisao para parada com sucesso (padrão = 1e-8)
#    tau   : ajuste de passos (padrão = 0.99995)
#    epsz  : epsilon para z inicial (padrão = 30)
#    saidas: imprimi saídas na tela? true(padrão) ou false
#
# Saida:
#    x:       solução ou último iterando
#    f:       valor final da FO
#    iter:    número de iterações executadas
#    status:
#       0 = sucesso
#       1 = problema ilimitado ou inviável
#       2 = maxit atingido
#       3 = A tem posto numericamente deficiente
#       4 = erro em Cholesky de ADA' (ADA' numericamente singular?)
#       9 = outro erro
#
# Para testar o método nos problemas pré-selecionados, execute
#   testar()
# ou
#   testar(saidas=true)
# Esta função fará o download dos problemas automaticamente.
#
# Exemplos de uso da função principal:
#   pdafimescala(A,b,c)
#   x, f, iter, status = ptoint(A,b,c, maxit=100, eps=1e-8, tau=0.99995, epsz=30, saidas=true)
###################################################

# pacotes para manipulação matricial
using LinearAlgebra, SparseArrays

# pacotes impressão e gravação de resultados
using Printf, DataFrames

# pacote para download de matrizes de https://sparse.tamu.edu
using MatrixDepot

################################
# MÉTODO PRIMAL DUAL AFIM ESCALA
################################

function pdafimescala(A,b,c; maxit=100, eps=1e-8, tau=0.99995, epsz=30, saidas=true)

    # Lê dimensões
    m,n = size(A)

    # Verifica se as dimensões dos dados são compatíveis
    if (length(b) != m) | (length(c) != n) | (m > n)
        error("Dados com dimensões incompatíveis!")
    end

    # Certifica-se que b e c são matrizes coluna
    if size(b)[1] == 1; b = Vector(b'[:,1]); end
    if size(c)[1] == 1; c = Vector(c'[:,1]); end

    #####################
    # Inicialização
    #####################

    status = 9
    iter   = 0
    viab   = Inf
    KKT    = Inf

    # Realoca A em estrutura esparsa, caso não esteja
    A = sparse(A)

    # c e b são vetores cheios
    c = Vector(c[:,1])
    b = Vector(b[:,1])

    # Guarda normas 1 de b e c para teste de convergência
    normb = norm(b,1)
    normc = norm(c,1)

    # Tenta fatorar A*A' por Cholesky.
    # Caso falhe, é porque A não tem posto completo. Neste caso retorna erro.
    # Caso tenha sucesso, capturamos a permutação p obtida por AMD.
    # "Symmetric" garante que internamente a matriz seja tratada como simétrica perfeita.
    # check=true faz cholesky retornar erro caso G*G' != P'(A*A')P
    C = []

    try
        C  = cholesky(Symmetric(sparse(A*A')), check=true, perm=nothing)
    catch
        println("A tem posto numericamente deficiente. Não se pode continuar!")
        status = 3

        return NaN,NaN,iter,viab,KKT,status
    end

    pAMD = C.p  # permutação de AMD

    # A permutação usada internamente por cholesky é gravada em C.p.
    # Assim, temos
    #   G*G' = (A*A')[C.p,C.p]
    # Dado um vetor v,
    #   v[C.p]
    # retorna P'*v. Para retornar P*v, precisamos da permutação
    # inversa de C.p:
    invpAMD = Vector(1:m)
    invpAMD[pAMD] = Vector(1:m)

    #####################
    # Ponto inicial
    #####################

    # Calcula A'(AA')^-1 b
    xbar = A'*((A*A')\b)

    # x, y, z iniciais
    x = max.( xbar , maximum( [ -minimum(xbar); 100; 0.01*normb/norm(A,1) ] ) )
    y = zeros(m,1)
    z = abs.(c) + epsz * (c .>= 0)
    z[(c .> -epsz) .& (c .< 0)] .= epsz

    # x, y, z, resíduos e direções são vetores cheios ao longo do processo
    x = Vector(x[:,1])
    y = Vector(y[:,1])
    z = Vector(z[:,1])

    # avalia FO primal f(x)=c^t x no ponto inicial
    f = dot(c,x)

    # Limite para critério de inviabilidade
    M = 1e+6 * norm([x;z],1)

    # inicia passos = Inf para correta impressão da 1a iteração
    alphap = Inf
    alphad = Inf

    # LOOP PRINCIPAL
    while iter < maxit

        # Calcula resíduos
        rp = b - A*x
        rd = c - A'*y - z
        rc = -x.*z  # =-XZe=-(x1z1, ... ,xnzn)

        #####################
        # Critério de parada
        #####################

        # Teste de convergência KKT
        KKTp = norm(rp,1)/(1 + normb)  # viabilidade primal relativa
        KKTd = norm(rd,1)/(1 + normc)  # viabilidade dual relativa
        fd   = dot(b,y)
        KKTc = abs(f - fd)/(1 + abs(f) + abs(fd))  # brecha de dualidade relativa

        KKT  = maximum([KKTp KKTd KKTc])  # medida de otimalidade (viab primal + dual + dualidade)

        viab = norm(rp,1)  # viabilidade primal absoluta

        # imprimi dados da iteração na tela
        if saidas
            if mod(iter,10) == 0
                @printf("\nit     \tf         viab      alphamax  KKT")
                @printf("\n===========================================================\n")
            end
            @printf("%d\t%8.2e  %8.2e  %8.2e  %8.2e\n",
                    iter, f, viab, max(alphap,alphad), KKT)
        end

        # Parar com sucesso?
        if KKT < eps
            println("\n--- SOLUÇÃO ENCONTRADA! ---")
            status = 0

            # retorno da função principal
            return x,f,iter,viab,KKT,status
        end

        # Teste de ilimitabilidade ou inviabilidade: parar com insucesso?
        # TEOREMA: se P é inviável ou ilimitado então |x|->Inf ou |z|->Inf
        # Assim, um critério para "suspeitar" que P não tem solução é |(x,z)| >> 1
        if norm([x,z],1) > M
            println("\n--- Problema provavelmente inviável ou ilimitado ---")
            status = 1
            f = NaN

            # retorno da função principal
            return x,f,iter,viab,KKT,status
        end

        # D = Z^{-1} X  (D será um vetor aqui, pois só precisamos armazenar a diagonal)
        D = x./z

        #####################
        # Direção Newtoniana
        #####################

        # Decomposição de Cholesky de ADA'.
        # spdiagm(0 => D) é a matriz diagonal esparsa. "0=>D" indica diagonal principal = D.
        AD = sparse(A*spdiagm(0 => D))

        # Usa a permutação pAMD calculada no inicio do método
        try
            C  = cholesky(Symmetric(sparse(AD*A')), check=true, perm=pAMD)
        catch
            println("Cholesky de ADA' falhou!")
            status = 4

            # retorno da função principal
            return x,f,iter,viab,KKT,status
        end

        # G será o fator de Cholesky triangular INFERIOR
        # Portanto temos GG' = P'(ADA')P.
        G  = sparse(C.L)

        # (ADA') dy = rp + AD(rd - X^-1 rc)
        #        dx = D(A'dy - rd + X^-1 rc)
        #        dz = X^-1 (rc - Zdx)
        # Temos GG' = P'(ADA')P. Resolvemos então
        #  (i) GG'h = P'(ADA')(Ph) = P'(rp + AD(rd - X^-1 rc))
        # (ii)   dy = Ph
        rcx = rc./x                    # rcx = X^-1 rc
        r   = rd - rcx                 # r  = rd - X^-1 rc
        s   = G\(rp + AD*r)[pAMD]      # Gs = P'(rp + AD(rd - X^-1 tc))
        tdy = G'\s                     # G' tdy = s
        dy  = tdy[invpAMD]             # dy = P tdy
        dx  = AD'*dy - D.*r
        dz  = rcx - (dx.*z)./x

        #####################
        # Tamanhos de passo (<= 1)
        #####################
        xdx = x./dx  # =(x1/dx1, ... ,xn/dxn)
        zdz = z./dz  # =(z1/dz1, ... ,zn/dzn)
        alphap = minimum( [ 1 ; -tau * xdx[dx .< 0] ] )  # .< testa < em cada coordenada
        alphad = minimum( [ 1 ; -tau * zdz[dz .< 0] ] )

        #####################
        # Atualização do ponto corrente
        #####################
        x = x + alphap*dx
        y = y + alphad*dy
        z = z + alphad*dz

        # Recalcula FO primal
        f = dot(c,x)

        iter += 1

    end

    # Neste ponto, maxit foi atingido... retorna último x obtido
    println("\n--- Número máximo de iterações atingido ---")
    status = 2

    return x,f,iter,viab,KKT,status
end
################################
################################

# Função para testes com problemas da NETLIB
# http://www.netlib.org/lp/data/index.html
function testar(; saidas=false)

    # PARAMETROS SUGERIDOS
    maxit = 100
    eps   = 1e-8
    tau   = 0.99995
    epsz  = 30

    # Problemas na forma padrão selecionados da NETLIB
    #   Fonte dos arquivos: grupo "LPnetlib" - https://sparse.tamu.edu/LPnetlib
    netlib = [
        "LPnetlib/lp_adlittle",
        "LPnetlib/lp_afiro",
        "LPnetlib/lp_agg",
        "LPnetlib/lp_agg2",
        "LPnetlib/lp_agg3",
        "LPnetlib/lp_bandm",
        "LPnetlib/lp_beaconfd",
        "LPnetlib/lp_blend",
        "LPnetlib/lp_bnl2",
        "LPnetlib/lp_d2q06c",
        "LPnetlib/lp_e226",
        "LPnetlib/lp_fffff800",
        "LPnetlib/lp_osa_60",
        "LPnetlib/lp_israel",
        "LPnetlib/lp_lotfi",
        "LPnetlib/lp_maros_r7",
        "LPnetlib/lp_osa_07",
        "LPnetlib/lp_osa_14",
        "LPnetlib/lp_osa_30",
        "LPnetlib/lp_sc50a",
        "LPnetlib/lp_sc105",
        "LPnetlib/lp_sc205",
        "LPnetlib/lp_sc50b",
        "LPnetlib/lp_scagr25",
        "LPnetlib/lp_scagr7",
        "LPnetlib/lp_scfxm1",
        "LPnetlib/lp_scfxm2",
        "LPnetlib/lp_scfxm3",
        "LPnetlib/lp_scrs8",
        "LPnetlib/lp_scsd1",
        "LPnetlib/lp_scsd6",
        "LPnetlib/lp_scsd8",
        "LPnetlib/lp_sctap1",
        "LPnetlib/lp_sctap2",
        "LPnetlib/lp_sctap3",
        "LPnetlib/lp_share1b",
        "LPnetlib/lp_share2b",
        "LPnetlib/lp_woodw",
        "LPnetlib/lp_stocfor1",
        "LPnetlib/lp_stocfor2",
        "LPnetlib/lp_stocfor3",
        "LPnetlib/lp_truss",
        "LPnetlib/lpi_itest2", #inviável
        "LPnetlib/lpi_itest6", #inviável
        "LPnetlib/lpi_bgindy", #inviável
        "LPnetlib/lpi_bgprtr", #inviável
        "LPnetlib/lpi_klein1", #inviável
        "LPnetlib/lpi_klein2", #inviável
        "LPnetlib/lpi_klein3"  #inviável
    ];

    # cabeçalho tabela com resultados
    resultados = DataFrame(Prob=[],n=[],m=[],st=[],it=[],f=[],viab=[],KKT=[],tempo=[])

    # executa pontos interiores em cada problema
    for p in netlib

        # lê problema (e baixa caso não tenha sido baixado)
        PL = mdopen(p);
        nomeprob = p[10:end]

        # PL é uma estutura com as seguintes propriedades:
        # PL.A, PL.b : sistema Ax=b
        # PL.c       : coeficiente função objetivo
        # PL.n, PL,m : número de variáveis e restrições
        # PL.nnz     : número de não zeros em A

        PLstatus = -1
        PLiter   = -1
        PLviab   = -1
        PLKKT    = -1
        PLf      = NaN
        PLtempo  = -1

        try
            println("\n==================")
            println("Prob: "*p)
            println("n="*string(PL.n)*", m="*string(PL.m))
            println("nnz em A: "*string(PL.nnz))
            println("==================")

            # aplica pontos interiores com reordenamento AMD
            PLtempo = @elapsed ~, PLf, PLiter, PLviab, PLKKT, PLstatus =
                pdafimescala(PL.A, PL.b, PL.c, maxit=maxit, eps=eps, tau=tau, epsz=epsz, saidas=saidas);

            print("Tempo (seg): ")
            println(PLtempo)

            # grava linha em resultados
            push!(resultados, (nomeprob, PL.n, PL.m, PLstatus, PLiter, PLf, PLviab, PLKKT, PLtempo))

        catch
            println("ERRO ao resolver $(nomeprob)!")
        end

    end

    # ordena resultados por nome do problema (coluna 1)
    sort!(resultados, 1)

    # salva resultados em TXT
    txt = open("resultados.txt", "w")
    write(txt, @sprintf("%s", resultados))
    close(txt)

    # Você pode salvar o resultado em planilhas.
    # Veja https://leonardosecchin.github.io/juliaopt_ex11/
end

max_stars_repo_name: leonardosecchin/temp

max_issues_repo_name: leonardosecchin/temp

max_forks_repo_name: leonardosecchin/temp


text:
# SPD-License-Identifier: MIT

using ArgParse
using JLD2
using MarkovModels

function main(args)
    hmmsdata = load(args["hmms"])
    lexicon = load(args["lexiconfsm"])
    hmms = hmmsdata["units"]
    pdfid_mapping = hmmsdata["pdfid_mapping"]
    sil = args["sil-sym"]
    bprob = args["between-silprob"]
    eprob = args["edge-silprob"]

    SF = LogSemifield{Float32}
    jldopen(args["alifsms"], "w") do f
        open(args["ali"], "r") do rf

            for line in eachline(rf)
                tokens = split(line)
                uttid, sentence = tokens[1], tokens[2:end]

                fsm = VectorFSM{SF}()

                # Initial silence state.
                prev = addstate!(fsm, sil;
                                 initweight = SF(log(eprob)),
                                 finalweight = zero(SF))

                for (i, word) in enumerate(sentence)
                    if word ∉ keys(lexicon)
                        @info "no pronunciation found for $word"
                        word = "<UNK>"
                    end

                    initweight = i == 1 ? SF(log(1-eprob)) : zero(SF)
                    finalweight = i == length(sentence) ? SF(log(1-eprob)) : zero(SF)
                    s = addstate!(fsm, word; initweight, finalweight)
                    addarc!(fsm, prev, s)

                    if 1 < i
                        silstate = addstate!(fsm, sil)
                        addarc!(fsm, prev, silstate, SF(log(bprob)))
                        addarc!(fsm, silstate, s, SF(log(bprob)))
                    end
                    prev = s
                end

                # Final silence state.
                silstate = addstate!(fsm, sil;
                                     initweight = zero(SF),
                                     finalweight = SF(log(eprob)))
                addarc!(fsm, prev, silstate)

                # G: Grammar L: Lexicon H: HMM
                GL = HierarchicalFSM(fsm |> renormalize, lexicon)
                GLH = HierarchicalFSM(GL, hmms)
                f["$uttid"] = MatrixFSM(GLH, pdfid_mapping, t -> (t[2], t[3]))
            end
        end
    end
end

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--sil-sym", "-s"
            default = "<SIL>"
            help = "optional silence symbol"
        "--edge-silprob", "-e"
            arg_type = Float64
            default = 0.8
            help = "probability of optional silence at the edges of the utterance"
        "--between-silprob", "-b"
            arg_type = Float64
            default = 0.2
            help = "probability of optional silence between words"
        "hmms"
            required = true
            help = "per-unit hmms"
        "lexiconfsm"
            required = true
            help = "input lexicon fsm"
        "ali"
            required = true
            help = "alignment file"
        "alifsms"
            required = true
            help = "output alignment fsms"
    end
    s.description = """
    Build alignment fsms.
    """
    parse_args(s)
end

args = parse_commandline()
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
##### MEAN-VARIANCE OPTIMIZATION
#This is not supposed to produce great results, since mean and covariance
#is supposed to be very exact before good weights are found
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
include("CEOSupportFunctions.jl")
println("Leeeeroooy Jenkins")
#Esben's path
#path = "$(homedir())/Documents/GitHub/Thesis/Data/IndexData"

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"
path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/IndexDataDiff/"

trainingSize = 240
possibilities = 5
industries = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
industriesTotal = length(industries)

modelMatrix = zeros(industriesTotal, possibilities)
#noDurModel = [1 0 1 1 1]
#noDurModel = [1 0 1 0 1]
noDurModel = [1 0 0 1 1]
testModel = [1 0 0 1 1]
modelMatrix[1, :] = noDurModel
for i=2:industriesTotal
    modelMatrix[i, :] = noDurModel
end

XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)

XArrays, YArrays = generateXandYs(industries, modelMatrix)

nGammas = 5
standY = YArrays[1]
SSTO = sum((standY[i]-mean(standY[:]))^2 for i=1:length(standY))
lambdaValues = log.(logspace(100, SSTO/2, nGammas))
nRows = size(standY)[1]
amountOfModels = nGammas^4

modelConfig = zeros(amountOfModels, 4)
counter = 1
for l1 = 1:nGammas
    for l2 = 1:nGammas
        for l3 = 1:nGammas
            for l4 = 1:nGammas
                modelConfig[counter,:] = [lambdaValues[l1] lambdaValues[l2] lambdaValues[l3] lambdaValues[l4]]
                counter += 1
            end
        end
    end
end
modelConfig

#Initialization of parameters
w1N = repeat([0.1], outer = 10) #1/N weights
gamma = 2.4 #risk aversion
validationPeriod = 5
PMatrix = zeros(nRows-trainingSize, amountOfModels)
return1NMatrix = zeros(nRows-trainingSize)
returnCEOMatrix = zeros(nRows-trainingSize, amountOfModels)
returnPerfectMatrix = zeros(nRows-trainingSize)

bestModelAmount = 5
bestModelConfigs = zeros(bestModelAmount, 4)
bestModelIndexes = zeros(bestModelAmount)

weightsPerfect = zeros(nRows-trainingSize, 10)
weightsCEO     = zeros(nRows-trainingSize, 10)

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEO/"
weightsPerfect = Array{Float64}(CSV.read(path*"weightsPerfect.csv",header=false, datarow=1, nullable=false))
returnPerfectMatrix = Array{Float64}(CSV.read(path*"returnPerfectMatrix.csv", header=false, datarow=1, nullable=false))

for i = 1:5
    for j = 0:24
        #path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
        path = "/zhome/9f/d/88706/SpecialeCode/Results/CEO/"
        PMatrixTemp = Array{Float64}(CSV.read(path*string(i)*"_"*string(j)*"_PMatrix.csv",header=false, datarow=1, nullable=false, types=[Float64]))
        return1NMatrixTemp = Array{Float64}(CSV.read(path*string(i)*"_"*string(j)*"_return1NMatrix.csv",header=false, datarow=1, nullable=false, types=[Float64]))
        returnCEOMatrixTemp = Array{Float64}(CSV.read(path*string(i)*"_"*string(j)*"_returnCEOMatrix.csv",header=false, datarow=1, nullable=false, types=[Float64]))

        PMatrix[i,:] += PMatrixTemp
        return1NMatrix[i] += return1NMatrixTemp[i]
        returnCEOMatrix[i,:] += returnCEOMatrixTemp
    end
end

#Load weightsPerfect and return perfect matrix
#Load PMatrix

modelMeans = mean(PMatrix, 1)
for i = 1:bestModelAmount
    bestModelIndex = indmax(modelMeans)
    bestModelConfigs[i,:] = modelConfig[bestModelIndex, :]
    bestModelIndexes[i] = bestModelIndex
    #place 0 at index place now and take new max as "next best model"
    modelMeans[bestModelIndex] = 0
end

inputArg1 = 83 #83 breaks the trainingSize
inputArg1 = parse(Int64,ARGS[1])
if 10+validationPeriod+(10*inputArg1) <= nRows-trainingSize-2
    println("Starting CEO Validation loop")
    for t=validationPeriod+1+(10*inputArg1):10+validationPeriod+(10*inputArg1)#(nRows-trainingSize-2)
        println("Time $t/50")
        trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)

        for m = 1:bestModelAmount
            betaArray, U = @time(runCEO(trainingXArrays, trainingYArrays, bestModelConfigs[m, :], gamma))
            expectedReturns = generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)

            #Need to send OOSRow to mean-variance optimization to get "perfect information" since validationY is the values in OOSRow[1:10]
            valY = zeros(10)
            for i = 1:10
                valY[i] = validationY[i][1]
            end
            return1N, returnCEO, wStar = performMVOptimization(expectedReturns, U, gamma, valY, valY)
            weightsCEO[t, 1:10]     = wStar
            return1NMatrix[t]      = return1N
            returnCEOMatrix[t, Int64(bestModelIndexes[m])]     = returnCEO
            returnPerfect = returnPerfectMatrix[t]
            println("1N returns is $return1N, returnPerfect is $returnPerfect and returnCEO is $returnCEO")
            PMatrix[t, Int64(bestModelIndexes[m])] = calculatePvalue(return1N, returnPerfect, returnCEO)
        end
    end
elseif validationPeriod+(10*inputArg1) <= nRows-trainingSize-2
    for t=validationPeriod+1+(10*inputArg1):(nRows-trainingSize-2)
        println("Time $t/50")
        trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)

        for m = 1:bestModelAmount
            betaArray, U = @time(runCEO(trainingXArrays, trainingYArrays, bestModelConfigs[m, :], gamma))
            expectedReturns = generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)

            #Need to send OOSRow to mean-variance optimization to get "perfect information" since validationY is the values in OOSRow[1:10]
            valY = zeros(10)
            for i = 1:10
                valY[i] = validationY[i][1]
            end
            return1N, returnCEO, wStar = performMVOptimization(expectedReturns, U, gamma, valY, valY)
            weightsCEO[t, 1:10]     = wStar
            return1NMatrix[t]      = return1N
            returnCEOMatrix[t, Int64(bestModelIndexes[m])]     = returnCEO
            returnPerfect = returnPerfectMatrix[t]
            println("1N returns is $return1N, returnPerfect is $returnPerfect and returnCEO is $returnCEO")
            PMatrix[t, Int64(bestModelIndexes[m])] = calculatePvalue(return1N, returnPerfect, returnCEO)
        end
    end
end

combinedPortfolios = hcat(returnPerfectMatrix[1:nRows-trainingSize-2, 1], return1NMatrix[1:nRows-trainingSize-2, 1],
    returnCEOMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)], PMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)])
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEO/"
writedlm(path*string(inputArg1)*"_returnPvalueOutcome1to200.csv", combinedPortfolios, ",")
println("Finished Everything")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
@testset "Testing sampler.jl" begin

    n = 4
    mu = collect(range(0.2, stop=0.6, length=n))
    @testset "check positive definite" begin
        Sigma = [3.0 1.0 1.0 1.0;
                 1.0 2.0 1.0 1.0;
                 1.0 1.0 2.0 1.0;
                 1.0 1.0 1.0 1.0]
        mvns = MVNSampler(mu, Sigma)
        @test isapprox(mvns.Q * mvns.Q', mvns.Sigma)
    end

    @testset "check positive semi-definite zeros" begin
        mvns = MVNSampler(mu, zeros(n, n))
        @test rand(mvns) == mu
    end

    @testset "check positive semi-definite ones" begin
        mvns = MVNSampler(mu, ones(n, n))
        c = rand(mvns)-mvns.mu
        @test all(broadcast(isapprox,c[1],c))
    end

    @testset "check positive semi-definite 1 and -1/(n-1)" begin
        Sigma = -1/(n-1)*ones(n, n) + n/(n-1)*Matrix(Matrix(I, n, n))
        mvns = MVNSampler(mu, Sigma)
        @test isapprox(sum(rand(mvns)) , sum(mu), atol=1e-4, rtol=1e-4)
    end

    @testset "check non-positive definite" begin
        Sigma = [2.0 1.0 3.0 1.0;
                 1.0 2.0 1.0 1.0;
                 3.0 1.0 2.0 1.0;
                 1.0 1.0 1.0 1.0]
        @test_throws ArgumentError MVNSampler(mu, Sigma)
    end

    @testset "check availability of rank deficient matrix" begin
        A = randn(n,n)
        for r=1:n-2
            r=2
            for i = 1:r
                A[:, end+1-i] = sum(A[:, 1:end-r], dims = 2)
            end
            Sigma = A * A'
            @test typeof(MVNSampler(mu,Sigma)) <: MVNSampler
        end
    end

    @testset "test covariance matrices of Int and Rational" begin
        n = 2
        mu = zeros(2)
        for T in [Int, Rational{Int}]
            Sigma = Matrix{T}(I, n, n)
            @test typeof(MVNSampler(mu, Sigma)) <: MVNSampler
        end
    end

end  # @testset


text:
#=
  Original Author: Anastasia Yendiki

  Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH"

  Terms and conditions for use, reproduction, distribution and contribution
  are found in the 'FreeSurfer Software License Agreement' contained
  in the file 'LICENSE' found in the FreeSurfer distribution, and here:

  https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense

  Reporting: freesurfer@nmr.mgh.harvard.edu
=#

using LinearAlgebra, StaticArrays, Statistics

export DTI, adc_fit, dti_fit, dti_write

"Container for outputs of a DTI fit"
struct DTI
  s0::MRI
  eigval1::MRI
  eigval2::MRI
  eigval3::MRI
  eigvec1::MRI
  eigvec2::MRI
  eigvec3::MRI
  rd::MRI
  md::MRI
  fa::MRI
end

"""
    ADCwork{T}

Pre-allocated workspace for ADC fit computations

- `T::DataType`                : Data type for computations (default: `Float32`)
- `nvol::Int`                  : Number of volumes in DWI series
- `ib0::Vector{Bool}`          : Indicator for b=0 volumes [`nvol`]
- `ipos::Vector{Vector{Bool}}` : Indicator for volumes with positive DWI signal [`nvol`]
- `logs::Vector{Vector{T}}`    : Logarithm of DWI signal [`nvol`]
- `d::Vector{Vector{T}}`       : Linear system solution vector [2]
- `A::Matrix{T}`               : System matrix [nvol x 2]
- `pA::Matrix{T}`              : Pseudo-inverse of system matrix [2 x nvol]
"""
struct ADCwork{T}
  nvol::Int
  ib0::Vector{Bool}
  ipos::Vector{Vector{Bool}}
  logs::Vector{Vector{T}}
  d::Vector{Vector{T}}
  A::Matrix{T}
  pA::Matrix{T}

  function ADCwork(bval::Vector{Float32}, T::DataType=Float32)

    # Number of volumes in DWI series
    nvol = length(bval)

    # Indicator for b=0 volumes
    ib0  = (bval .== minimum(bval))

    # Indicator for volumes with positive DWI signal
    ipos = [Vector{Bool}(undef, nvol) for tid in 1:Threads.nthreads()]

    # Log of DWI signal
    logs = [Vector{T}(undef, nvol) for tid in 1:Threads.nthreads()]

    # Linear system solution vector
    d    = [Vector{T}(undef, 2) for tid in 1:Threads.nthreads()]

    # System matrix
    A = Matrix{T}(undef, nvol, 2)

    A[:,1] = -bval
    A[:,2] = ones(T, nvol)

    # Pseudo-inverse of system matrix
    pA = pinv(A)

    new{T}(
      nvol,
      ib0,
      ipos,
      logs,
      d,
      A,
      pA
    )
  end
end

"""
    DTIwork{T}

Pre-allocated workspace for DTI fit computations

- `T::DataType`                : Data type for computations (default: `Float32`)
- `nvol::Int`                  : Number of volumes in DWI series
- `ib0::Vector{Bool}`          : Indicator for b=0 volumes [`nvol`]
- `ipos::Vector{Vector{Bool}}` : Indicator for volumes with positive DWI signal [`nvol`]
- `logs::Vector{Vector{T}}`    : Logarithm of DWI signal [`nvol`]
- `d::Vector{Vector{T}}`       : Linear system solution vector [7]
- `A::Matrix{T}`               : System matrix [nvol x 7]
- `pA::Matrix{T}`              : Pseudo-inverse of system matrix [7 x nvol]
"""
struct DTIwork{T}
  nvol::Int
  ib0::Vector{Bool}
  ipos::Vector{Vector{Bool}}
  logs::Vector{Vector{T}}
  d::Vector{Vector{T}}
  A::Matrix{T}
  pA::Matrix{T}

  function DTIwork(bval::Vector{Float32},
                   bvec::Matrix{Float32}, T::DataType=Float32)

    # Number of volumes in DWI series
    nvol = length(bval)

    # Indicator for b=0 volumes
    ib0  = (bval .== minimum(bval))

    # Indicator for volumes with positive DWI signal
    ipos = [Vector{Bool}(undef, nvol) for tid in 1:Threads.nthreads()]

    # Log of DWI signal
    logs = [Vector{T}(undef, nvol) for tid in 1:Threads.nthreads()]

    # Linear system solution vector
    d    = [Vector{T}(undef, 7) for tid in 1:Threads.nthreads()]

    # System matrix
    A = Matrix{T}(undef, nvol, 7)

    A[:,1] = bvec[:,1].^2
    A[:,2] = 2*bvec[:,1].*bvec[:,2]
    A[:,3] = 2*bvec[:,1].*bvec[:,3]
    A[:,4] = bvec[:,2].^2
    A[:,5] = 2*bvec[:,2].*bvec[:,3]
    A[:,6] = bvec[:,3].^2

    A[:,1:6] .*= -bval

    A[:,7] = ones(T, nvol)

    # Pseudo-inverse of system matrix
    pA = pinv(A)

    new{T}(
      nvol,
      ib0,
      ipos,
      logs,
      d,
      A,
      pA
    )
  end
end

"""
    adc_fit(dwi::MRI, mask::MRI)

Fit the apparent diffusion coefficient (ADC) to DWIs and return it as an
`MRI` structure.
"""
function adc_fit(dwi::MRI, mask::MRI)

  if isempty(dwi.bval)
    error("Missing b-value table from input DWI structure")
  end

  W = ADCwork(dwi.bval)

  adc = MRI(mask, 1, Float32)
  s0  = MRI(mask, 1, Float32)

  Threads.@threads for iz in 1:size(dwi.vol, 3)
    for iy in 1:size(dwi.vol, 2)
      for ix in 1:size(dwi.vol, 1)
        mask.vol[ix, iy, iz] == 0 && continue

        adc.vol[ix, iy, iz],
        s0.vol[ix, iy, iz] = adc_fit(dwi.vol[ix, iy, iz, :], W)
      end
    end
  end

  return adc, s0
end

"""
    adc_fit(dwi::Vector{T}, W::ADCwork{T}) where T<:AbstractFloat

Fit the ADC for a single voxel
"""
function adc_fit(dwi::Vector{T}, W::ADCwork{T}) where T<:AbstractFloat

  tid = Threads.threadid()

  # Only use positive DWI values to fit the model
  W.ipos[tid] .= dwi .> 0
  npos = sum(W.ipos[tid])

  if npos == W.nvol
    W.logs[tid] .= log.(dwi)
    mul!(W.d[tid], W.pA, W.logs[tid])
  elseif npos > 6 && any(W.ipos[tid][W.ib0])
    mul!(W.d[tid], pinv(W.A[W.ipos[tid], :]), log.(dwi[W.ipos[tid]]))
  else
    return zero(T), zero(T)
  end

  return W.d[tid][1], exp(W.d[tid][2])
end

"""
    dti_fit(dwi::MRI, mask::MRI)

Fit tensors to DWIs and return a `DTI` structure.
"""
function dti_fit(dwi::MRI, mask::MRI)

  if isempty(dwi.bval)
    error("Missing b-value table from input DWI structure")
  end

  if isempty(dwi.bvec)
    error("Missing gradient table from input DWI structure")
  end

  dti_fit_ls(dwi::MRI, mask::MRI)
end

"""
    dti_fit_ls(dwi::MRI, mask::MRI)

Perform least-squares fitting of tensors from DWIs and return a `DTI` structure.

If you use this method, please cite:
Peter Basser et al. (1994). Estimation of the effective self-diffusion tensor from the NMR spin echo. Journal of Magnetic Resonance Series B, 103(3), 247–254. https://doi.org/10.1006/jmrb.1994.1037
"""
function dti_fit_ls(dwi::MRI, mask::MRI)

  W = DTIwork(dwi.bval, dwi.bvec)

  S0    = MRI(mask, 1, Float32)
  Eval1 = MRI(mask, 1, Float32)
  Eval2 = MRI(mask, 1, Float32)
  Eval3 = MRI(mask, 1, Float32)
  Evec1 = MRI(mask, 3, Float32)
  Evec2 = MRI(mask, 3, Float32)
  Evec3 = MRI(mask, 3, Float32)
  RD    = MRI(mask, 1, Float32)
  MD    = MRI(mask, 1, Float32)
  FA    = MRI(mask, 1, Float32)

  Threads.@threads for iz in 1:size(dwi.vol, 3)
    for iy in 1:size(dwi.vol, 2)
      for ix in 1:size(dwi.vol, 1)
        mask.vol[ix, iy, iz] == 0 && continue

        S0.vol[ix, iy, iz],
        Eval1.vol[ix, iy, iz],
        Eval2.vol[ix, iy, iz],
        Eval3.vol[ix, iy, iz],
        Evec1.vol[ix, iy, iz, :],
        Evec2.vol[ix, iy, iz, :],
        Evec3.vol[ix, iy, iz, :],
        RD.vol[ix, iy, iz],
        MD.vol[ix, iy, iz],
        FA.vol[ix, iy, iz] = dti_fit_ls(dwi.vol[ix, iy, iz, :], W)
      end
    end
  end

  return DTI(S0, Eval1, Eval2, Eval3, Evec1, Evec2, Evec3, RD, MD, FA)
end

"""
    dti_fit_ls(dwi::Vector{T}, W::DTIwork{T}) where T<:AbstractFloat

Perform least-squares fitting of tensor for a single voxel
"""
function dti_fit_ls(dwi::Vector{T}, W::DTIwork{T}) where T<:AbstractFloat

  tid = Threads.threadid()

  # Only use positive DWI values to fit the model
  W.ipos[tid] .= dwi .> 0
  npos = sum(W.ipos[tid])

  if npos == W.nvol
    W.logs[tid] .= log.(dwi)
    mul!(W.d[tid], W.pA, W.logs[tid])
  elseif npos > 6 && any(W.ipos[tid][W.ib0])
    mul!(W.d[tid], pinv(W.A[W.ipos[tid], :]), log.(dwi[W.ipos[tid]]))
  else
    return zero(T), zero(T), zero(T), zero(T),
                    zeros(T, 3), zeros(T, 3), zeros(T, 3),
                    zero(T), zero(T), zero(T)
  end

  s0 = exp(W.d[tid][7])

  D = @SMatrix [ W.d[tid][1] 0           0;
                 W.d[tid][2] W.d[tid][4] 0;
                 W.d[tid][3] W.d[tid][5] W.d[tid][6] ]

  E = eigen(Symmetric(D, :L))

  return s0, E.values[3], E.values[2], E.values[1],
             E.vectors[:, 3], E.vectors[:, 2], E.vectors[:, 1],
             dti_maps(E.values[3], E.values[2], E.values[1])...
end

"""
    dti_maps(eigval1::T, eigval2::T, eigval3::T) where T<:AbstractFloat

Return the radial diffusivity (RD), mean diffusivity (MD), and fractional
anisotropy (FA) given the 3 eigenvalues the diffusion tensor.
"""
function dti_maps(eigval1::T, eigval2::T, eigval3::T) where T<:AbstractFloat

  rd = eigval2 + eigval3
  md = (eigval1 + rd) / 3
  rd /= 2

  fa = sqrt(((eigval1 - md)^2 + (eigval2 - md)^2 + (eigval3 - md)^2) /
            (eigval1^2 + eigval2^2 + eigval3^2) * T(1.5))

  return rd, md, fa
end

"""
    dti_write(dti::DTI, basename::String)

Write the volumes from a `DTI` structure that was created by `dti_fit()`
to files whose names start with the specified base name.
"""
function dti_write(dti::DTI, basename::String)

  for var in fieldnames(DTI)
    mri_write(getfield(dti, var), basename * "_" * string(var) * ".nii.gz")
  end
end

max_stars_repo_name: freesurfer/jupysurf

max_issues_repo_name: freesurfer/jupysurf

max_forks_repo_name: freesurfer/jupysurf


text:
mutable struct ShowSettings
    tableformat::PrettyTableFormat
    floatformat::String
    border::Crayon
    header::Crayon
    subheader::Crayon
    fixed::Crayon
    error::Crayon
    section::Crayon
    fixedpars::Bool
    ShowSettings() = new(unicode_rounded, "%10.4g",
                         crayon"light_blue", crayon"light_blue bold",
                         crayon"dark_gray bold", crayon"dark_gray",
                         crayon"light_red blink", crayon"yellow bold",
                         false)
end

const showsettings = ShowSettings()

printtable(args...; kw...) = pretty_table(args..., showsettings.tableformat; kw...,
                                          border_crayon=showsettings.border,
                                          header_crayon=showsettings.header,
                                          subheader_crayon=showsettings.subheader)

section(io, args...) = println(io, showsettings.section, args...)

function left(s::String, maxlen::Int)
    (length(s) <= maxlen)  &&  (return s)
    return s[1:maxlen]
end

function show(io::IO, dom::AbstractCartesianDomain)
    section(io, "Cartesian domain (ndims: ", ndims(dom), ", length: ", length(dom), ")")
    table = Matrix{Union{Int,Float64}}(undef, ndims(dom), 6)
    for i in 1:ndims(dom)
        a = dom[i]
        b = 0
        if length(a) > 1
            b = a .- circshift(a, 1)
            b = b[2:end]
        end
        table[i, 1] = i
        table[i, 2] = length(a)
        table[i, 3:6] = [minimum(a), maximum(a), minimum(b), maximum(b)]
    end
    printtable(io, table, ["Dim", "Size", "Min val", "Max val", "Min step", "Max step"],
               formatter=ft_printf(showsettings.floatformat, 3:6))
end

function show(io::IO, dom::AbstractLinearDomain)
    section(io, "Linear domain (ndims: ", ndims(dom), ", length: ", length(dom), ")")
    table = Matrix{Union{Int,Float64}}(undef, ndims(dom), 3)
    for i in 1:ndims(dom)
        table[i, 1] = i
        table[i, 2] = getaxismin(dom, i)
        table[i, 3] = getaxismax(dom, i)
    end
    printtable(io, table, ["Dim", "Min val", "Max val"],
               formatter=ft_printf(showsettings.floatformat, 2:3))
end

# Special case for Domain_1D: treat it as a Cartesian domain, despite it is a Linear one.
function show(io::IO, dom::Domain_1D)
    section(io, "Domain (ndims: ", ndims(dom), ", length: ", length(dom), ")")
    table = Matrix{Union{Int,Float64}}(undef, ndims(dom), 6)
    for i in 1:ndims(dom)
        a = dom[i]
        b = 0
        if length(a) > 1
            b = a .- circshift(a, 1)
            b = b[2:end]
        end
        table[i, 1] = i
        table[i, 2] = length(a)
        table[i, 3:6] = [minimum(a), maximum(a), minimum(b), maximum(b)]
    end
    printtable(io, table, ["Dim", "Size", "Min val", "Max val", "Min step", "Max step"],
               formatter=ft_printf(showsettings.floatformat, 3:6))
end

function show(io::IO, data::AbstractData)
    section(io, typeof(data), "   length: ", (length(data.val)))
    table = Matrix{Union{String,Float64}}(undef, 0, 7)

    names = fieldnames(typeof(data))
    for name in names
        a = getfield(data, name)
        nan = length(findall(isnan.(a))) + length(findall(isinf.(a)))
        a = a[findall(isfinite.(a))]
        table = vcat(table, [(nan > 0  ?  "⚠"  :  "") string(name) minimum(a) maximum(a) mean(a) median(a) std(a)])
    end
    printtable(io, table, ["⚠", "", "Min", "Max", "Mean", "Median", "Std. dev."],
               formatter=ft_printf(showsettings.floatformat, 3:7))
end

function preparetable(wcomp::WComponent)
    table = Matrix{Union{String,Float64}}(undef, 0, 7)
    fixed = Vector{Bool}()
    error = Vector{Bool}()
    comp = wcomp.comp
    cname = string(wcomp.cname)
    for (pname, param) in getparams(wcomp)
        (!showsettings.fixedpars)  &&  (wcomp.fixed  ||  param.fixed)  &&  continue
        ss = string(param._private.pname) * (param._private.index >= 1   ?   "["*string(param._private.index)*"]"  :  "")
        range = (param.fixed  ?  "FIXED"  :  strip(@sprintf("%7.2g:%-7.2g", param.low, param.high)))
        (range == "-Inf:Inf")  &&  (range = "")
        log = (param.log  ?  "LOG"  :  "")
        table = vcat(table, [cname ss param.val range log param.expr description(comp, param._private.pname)])
        push!(fixed, param.fixed)
        push!(error, !(param.low <= param.val <= param.high))
        cname = ""
    end
    return (table, fixed, error)
end

function preparetable(dict::OrderedDict{Symbol, WComponent})
    table = Matrix{Union{String,Float64}}(undef, 0, 7)
    fixed = Vector{Bool}()
    error = Vector{Bool}()
    hrule = Vector{Int}()
    for (cname, wcomp) in dict
        (t, f, e) = preparetable(wcomp)
        table = vcat(table, t)
        append!(fixed, f)
        append!(error, e)
        push!(hrule, length(error))
    end
    return (table, fixed, error, hrule)
end

show(io::IO, w::UI{WComponent}) = show(io, wrappee(w))
show(io::IO, wcomp::WComponent) = show(OrderedDict{Symbol, WComponent}(wcomp.cname => wcomp))

function show(io::IO, dict::OrderedDict{Symbol, WComponent})
    (table, fixed, error, hrule) = preparetable(dict)
    (length(table) > 0)  ||  (return nothing)
    printtable(io, table , ["Component" "Param." "Value" "Range" "Log" "Expr" "Notes"], alignment=:l,
               hlines=hrule, formatter=ft_printf(showsettings.floatformat, [3]),
               highlighters=(Highlighter((data,i,j) -> fixed[i], showsettings.fixed),
                             Highlighter((data,i,j) -> (error[i] && j==5), showsettings.error)))
end

show(io::IO, mime::MIME"text/plain", model::Model) = show(io, model)
show(io::IO, w::UI{Model}) = show(io, wrappee(w))
function show(io::IO, model::Model)
    _evaluate!(model)
    section(io, "Model components:")
    length(model.comp) != 0  || (return nothing)

    table = Matrix{Union{String,Float64}}(undef, 0, 4)
    fixed = Vector{Bool}()
    for (cname, wcomp) in model.comp
        ctype = split(string(typeof(wcomp.comp)), ".")
        (ctype[1] == "DataFitting")  &&   (ctype = ctype[2:end])
        ctype = join(ctype, ".")
        table = vcat(table, [string(cname) (wcomp.fixed  ?  "F"  :  "") ctype description(wcomp.comp)])
        push!(fixed, wcomp.fixed)
    end
    printtable(io, table, ["Component" "F" "Type" "Description"], alignment=:l,
               highlighters=(Highlighter((data,i,j) -> fixed[i], showsettings.fixed),))

    println(io)
    section(io, "Parameters:")
    show(io, model.comp)

    if length(model.instruments) == 0
        println(io)
        section(io, "Instrument(s): 0")
        return nothing
    end

    for i in 1:length(model.instruments)
        println(io)
        section(io, "Instrument #$i ")
        show(io, model.instruments[i])
    end

     tmp = length(model.index1d) - 1
     (tmp < 0)  &&  (tmp = 0)
     section(io, "Instrument(s): ", length(model.instruments),
               ".  Dataset(s) required for fitting: ", tmp)
end

function show(io::IO, instr::Instrument)
#     if ((isa(io, Base.TTY)  ||  isa(io, IOContext))  &&  (displaysize(io)[2] >= 80))
    show(io, instr.domain)
    println(io)

    (length(instr.compevals) == 0)  &&  (return nothing)
    section(io, "Evaluated expression(s):")
    table = Matrix{Union{String,Int,Float64}}(undef, length(instr.compevals)+length(instr.exprs), 8)
    fixed = Vector{Bool}()
    error = Vector{Bool}()
    hrule = Vector{Int}()

    for i in 1:length(instr.compevals)
        cname = instr.compnames[i]
        ceval = instr.compevals[i]
        result = ceval.result
        v = view(result, findall(isfinite.(result)))
        (length(v) == 0)  &&  (v = [NaN])
        nan = length(findall(isnan.(result)))
        inf = length(findall(isinf.(result)))
        table[i, 1] = ""
        table[i, 2] = string(cname)
        table[i, 3] = ceval.counter
        table[i, 4:6] = [minimum(v), maximum(v), mean(v)]
        table[i, 7] = (nan+inf > 0 ? "⚠" : "")
        table[i, 8] = ""
        push!(fixed, ceval.fixed)
        push!(error, (nan+inf > 0))
    end
    i0 = length(instr.compevals)

    for i in 1:length(instr.exprs)
        result = instr.exprevals[i]
        v = view(result, findall(isfinite.(result)))
        nan = length(findall(isnan.(result)))
        inf = length(findall(isinf.(result)))
        table[i0+i, 1] = (instr.exprcmp[i]  ?  "⇒"  :  "")
        table[i0+i, 2] = string(instr.exprnames[i])
        table[i0+i, 3] = instr.counter
        table[i0+i, 4:6] = [minimum(v), maximum(v), mean(v)]
        table[i0+i, 7] = (nan+inf > 0 ? "⚠" : "")
        table[i0+i, 8] = left(string(instr.exprs[i]), 30)
        push!(fixed, false)
        push!(error, (nan+inf > 0))
    end

    printtable(io, table, ["", "Label", "Counter", "Min", "Max", "Mean", "⚠", "Expr"], alignment=:l,
               hlines=[i0], formatter=ft_printf(showsettings.floatformat, 4:6),
               highlighters=(Highlighter((data,i,j) -> fixed[i], showsettings.fixed),
                             Highlighter((data,i,j) -> (error[i] && j==5), showsettings.error)))
    println(io)
end

function show(io::IO, par::Parameter)
    if par.fixed
        println(io, "Value : ", par.val, "   (FIXED)")
    else
        println(io, "Value : ", par.val, "  [", par.low , " : ", par.high, "]")
        if par.expr != ""
            println(io, "Expr : ", par.expr)
        end
    end
    if par.log
        println(io, "(use logarithmic value in fit)")
    end
end

show(io::IO, par::FitParam) = println(io, par.val, " ± ", par.unc)

show(io::IO, w::UI{FitComp}) = show(io, wrappee(w))

function preparetable(comp::FitComp)
    table = Matrix{Union{String,Float64}}(undef, 0, 4)
    fixed = Vector{Bool}()
    error = Vector{Bool}()

    for (pname, params) in comp.params
        if typeof(params) == Vector{FitParam}
            for ii in 1:length(params)
                par = params[ii]
                (!showsettings.fixedpars)  &&  (par.fixed)  &&  continue
                spname = string(pname) * "[" * string(ii) * "]"
                table = vcat(table, ["" spname par.val par.unc])
                push!(fixed, par.fixed)
                push!(error, !isfinite(par.unc))
            end
        else
            par = params
            (!showsettings.fixedpars)  &&  (par.fixed)  &&  continue
            spname = string(pname)
            table = vcat(table, ["" spname par.val par.unc])
            push!(fixed, par.fixed)
            push!(error, !isfinite(par.unc))
        end
    end
    return (table, fixed, error)
end

function show(io::IO, w::UI{FitResult})
    res = wrappee(w)
    section(io, "Best Fit results:")

    table = Matrix{Union{String,Float64}}(undef, 0, 4)
    fixed = Vector{Bool}()
    error = Vector{Bool}()
    hrule = Vector{Int}()
    for (cname, comp) in res.bestfit
        if length(comp.params) > 0
            (t, f, e) = preparetable(comp)
            (length(t) > 0)  ||  continue
            t[1,1] = string(cname)
            table = vcat(table, t)
            append!(fixed, f)
            append!(error, e)
            push!(hrule, length(error))
        end
    end
    printtable(io, table , ["Component" "Param." "Value" "Uncert."], alignment=:l,
               hlines=hrule, formatter=ft_printf(showsettings.floatformat, [3,4]),
               highlighters=(Highlighter((data,i,j) -> fixed[i], showsettings.fixed),
                             Highlighter((data,i,j) -> (error[i]  &&  (!fixed[i])  &&  (j==4)), showsettings.error)))

    println(io)
    println(io, @sprintf("    #Data  : %10d              Cost : %-10.5g", res.ndata, res.cost))
    println(io, @sprintf("    #Param : %10d              Red. : %-10.4g", res.ndata-res.dof, res.cost / res.dof))
    print(  io, @sprintf("    DOF    : %10d              ", res.dof))
    if res.log10testprob < -3
        println(io, @sprintf("Prob.: 10^%-10.4g", res.log10testprob))
    else
        println(io, @sprintf("Prob.: %10.4g", 10^res.log10testprob))
    end
    printstyled(io, "    Status :  ")
    if res.status == :Optimal
        printstyled(color=:green, io, @sprintf("%-15s", "Optimal"))
    elseif res.status == :NonOptimal
        printstyled(color=printcolorerr(), io, @sprintf("%-15s", "Non Optimal"))
    elseif res.status == :Warn
         printstyled(color=printcolorerr(), io, @sprintf("%-15s", "Warning"))
    elseif res.status == :Error
        printstyled(color=printcolorerr(), io, @sprintf("%-15s", "Error"))
    else
        printstyled(color=printcolorerr(), io, @sprintf("%-15s", "Unknown (" * string(res.status) * "), see fitter output"))
    end
    println(io, @sprintf("        Elapsed: %-10.4g s", res.elapsed))
end


text:
using Test

include("all-your-base.jl")

@testset "single bit one to decimal" begin
    @test all_your_base([1], 2, 10) == [1]
end

@testset "binary to single decimal" begin
    @test all_your_base([1, 0, 1], 2, 10) == [5]
end

@testset "single decimal to binary" begin
    @test all_your_base([5], 10, 2) == [1, 0, 1]
end

@testset "binary to multiple decimal" begin
    @test all_your_base([1, 0, 1, 0, 1, 0], 2, 10) == [4, 2]
end

@testset "decimal to binary" begin
    @test all_your_base([4, 2], 10, 2) == [1, 0, 1, 0, 1, 0]
end

@testset "trinary to hexadecimal" begin
    @test all_your_base([1, 1, 2, 0], 3, 16) == [2, 10]
end

@testset "hexadecimal to trinary" begin
    @test all_your_base([2, 10], 16, 3) == [1, 1, 2, 0]
end

@testset "15-bit integer" begin
    @test all_your_base([3, 46, 60], 97, 73) == [6, 10, 45]
end

@testset "empty list" begin
    @test all_your_base([], 2, 10) == [0]
end

@testset "single zero" begin
    @test all_your_base([0], 10, 2) == [0]
end

@testset "multiple zeros" begin
    @test all_your_base([0, 0, 0], 10, 2) == [0]
end

@testset "leading zeros" begin
    @test all_your_base([0, 6, 0], 7, 10) == [4, 2]
end

@testset "input base is one" begin
    @test_throws DomainError all_your_base([0], 1, 10)
end

@testset "input base is zero" begin
    @test_throws DomainError all_your_base([], 0, 10)
end

@testset "input base is negative" begin
    @test_throws DomainError all_your_base([1], -2, 10)
end

@testset "negative digit" begin
    @test_throws DomainError all_your_base([1, -1, 1, 0, 1, 0], 2, 10)
end

@testset "invalid positive digit" begin
    @test_throws DomainError all_your_base([1, 2, 1, 0, 1, 0], 2, 10)
end

@testset "output base is one" begin
    @test_throws DomainError all_your_base([1, 0, 1, 0, 1, 0], 2, 1)
end

@testset "output base is zero" begin
    @test_throws DomainError all_your_base([7], 10, 0)
end

@testset "output base is negative" begin
    @test_throws DomainError all_your_base([1], 2, -7)
end

@testset "both bases are negative" begin
    @test_throws DomainError all_your_base([1], -2, -7)
end


text:
using DFTK
using LinearAlgebra
using Test

@testset "Anyons" begin
    # Test that the magnetic potential satisfies ∇∧A = 2π ρref, ∇⋅A = 0
    x = 1.23
    y = -1.8
    ε = 1e-8
    M = 2.31
    σ = 1.81
    dAdx = (  DFTK.magnetic_potential_produced_by_ρref(x + ε, y, M, σ)
            - DFTK.magnetic_potential_produced_by_ρref(x,     y, M, σ)) / ε
    dAdy = (  DFTK.magnetic_potential_produced_by_ρref(x, y + ε, M, σ)
            - DFTK.magnetic_potential_produced_by_ρref(x, y,     M, σ)) / ε
    curlA = dAdx[2] - dAdy[1]
    divA = dAdx[1] + dAdy[2]
    @test norm(curlA - 2π*DFTK.ρref_real_2D(x, y, M, σ)) < 1e-4
    @test abs(divA) < 1e-6
end


text:
show_fix2(io::IO, f::Base.Fix2) = print(io, f.f, "(", f.x, ")")
Base.print(io::IO, f::Base.Fix2) = show_fix2(io, f)
Base.show(io::IO, ::MIME"text/plain", f::Base.Fix2) = show_fix2(io, f)


text:
using Test
using BOMoD
using BOMoD: filter_constrain, Combination, Compose_Construct_Constrains
using Random
using StatsBase

@testset "BOMod" begin
    include("Mod_test.jl")
    include("Filter_test.jl")
    include("Space_test.jl")
    include("Constrain_test.jl")
    include("Combination_test.jl")
    include("Constructs_test.jl")
end


text:
using AcuteML:
    Document,  # Not used, but required
    Node,  # Not used, but required
    parsexml,
    setnodename!,
    @aml,
    @att_str,
    @doc_str,
    @empty_str,
    @txt_str
using Pseudopotentials: Pseudopotential

export UPF, getdata

istrue(str) = occursin(r"t(rue)?"i, str)

function parsevec(str)
    vec = Float64[]
    for line in split(str, r"\R"; keepempty = false)
        append!(vec, map(x -> parse(Float64, x), split(line, r"[ \t]"; keepempty = false)))
    end
    return vec
end

@aml struct Info "PP_INFO"
    text::String, txt""
    inputfile::UN{String}, "PP_INPUTFILE"
end

@aml struct Header empty"PP_HEADER"
    generated::UN{String}, att"generated"
    author::UN{String}, att"author"
    date::UN{String}, att"date"
    comment::UN{String}, att"comment"
    element::String, att"element"
    pseudo_type::String, att"pseudo_type"
    relativistic::String, att"relativistic"
    is_ultrasoft::String, att"is_ultrasoft"
    is_paw::String, att"is_paw"
    is_coulomb::UN{String} = ".false.", att"is_coulomb"
    has_so::UN{String} = ".false.", att"has_so"
    has_wfc::String, att"has_wfc"
    has_gipaw::UN{String} = ".false.", att"has_gipaw"
    paw_as_gipaw::UN{String}, att"paw_as_gipaw"  # Suggested "required" but can be missing
    core_correction::String, att"core_correction"
    functional::String, att"functional"
    z_valence::Float64, att"z_valence"
    total_psenergy::UN{Float64} = 0, att"total_psenergy"
    wfc_cutoff::UN{Float64} = 0, att"wfc_cutoff"
    rho_cutoff::UN{Float64} = 0, att"rho_cutoff"
    l_max::Float64, att"l_max"
    l_max_rho::UN{Float64}, att"l_max_rho"  # Suggested "required" but can be missing
    l_local::UN{Int}, att"l_local"
    mesh_size::UInt, att"mesh_size"
    number_of_wfc::UInt, att"number_of_wfc"
    number_of_proj::UInt, att"number_of_proj"
end

@aml struct R "PP_R"
    text::String, txt""
end

@aml struct Rab "PP_RAB"
    text::String, txt""
end

@aml struct Mesh "PP_MESH"
    dx::UN{Float64}, att"dx"
    mesh::UN{Int}, att"mesh"
    xmin::UN{Float64}, att"xmin"
    rmax::UN{Float64}, att"rmax"  # Suggested "required" but can be missing
    zmesh::UN{Float64}, att"zmesh"
    r::R, "PP_R"
    rab::Rab, "PP_RAB"
end

@aml struct Chi "PP_CHI"
    n::UN{UInt}, att"n"
    l::UInt, att"l"
    index::UN{Int}, att"index"
    label::UN{String}, att"label"
    occupation::Float64, att"occupation"
    pseudo_energy::UN{Float64}, att"pseudo_energy"
    cutoff_radius::UN{Float64}, att"cutoff_radius"
    ultrasoft_cutoff_radius::UN{Float64}, att"ultrasoft_cutoff_radius"
    text::String, txt""
end

@aml struct Pswfc "PP_PSWFC"
    chi::Vector{Chi}, "PP_CHI"
end

@aml struct Vnl "PP_VNL"
    l::String, att"L"
    j::String, att"J"
    text::String, txt""
end

@aml struct Semilocal "PP_SEMILOCAL"
    chi::Vector{Vnl}, "PP_VNL"
end

@aml struct Nlcc "PP_NLCC"
    text::String, txt""
end

@aml struct Local "PP_LOCAL"
    text::String, txt""
end

@aml struct Beta "PP_BETA"
    angular_momentum::Int, att"angular_momentum"
    index::UN{Int}, att"index"
    label::UN{String}, att"label"
    cutoff_radius::UN{Float64}, att"cutoff_radius"
    cutoff_radius_index::UN{Int}, att"cutoff_radius_index"
    norm_conserving_radius::UN{Float64}, att"norm_conserving_radius"
    ultrasoft_cutoff_radius::UN{Float64}, att"ultrasoft_cutoff_radius"  # Suggested "required" but can be missing
    text::String, txt""
end

@aml struct Dij "PP_DIJ"
    text::String, txt""
end

@aml struct Q "PP_Q"
    text::String, txt""
end

@aml struct Multipoles "PP_MULTIPOLES"
    text::String, txt""
end

@aml struct Qfcoeff "PP_QFCOEFF"
    text::String, txt""
end

@aml struct Rinner "PP_RINNER"
    text::String, txt""
end

@aml struct Qij "PP_QIJ"
    first_index::UN{Int}, att"first_index"
    second_index::UN{Int}, att"second_index"
    composite_index::UN{Int}, att"composite_index"
    is_null::UN{String}, att"is_null"
    text::String, txt""
end

@aml struct Qijl "PP_QIJL"
    angular_momentum::Int, att"angular_momentum"
    first_index::UN{Int}, att"first_index"
    second_index::UN{Int}, att"second_index"
    composite_index::UN{Int}, att"composite_index"
    is_null::UN{String}, att"is_null"
    text::String, txt""
end

@aml struct Augmentation "PP_AUGMENTATION"
    q_with_l::String, att"q_with_l"
    nqf::UN{Int}, att"nqf"
    nqlc::UN{Int}, att"nqlc"
    shape::UN{String}, att"shape"
    iraug::UN{Int}, att"iraug"
    raug::UN{Float64}, att"raug"
    augmentation_epsilon::UN{Float64}, att"augmentation_epsilon"
    cutoff_r::UN{Float64}, att"cutoff_r"
    cutoff_r_index::UN{Int}, att"cutoff_r_index"
    l_max_aug::UN{Int}, att"l_max_aug"
    q::Q, "PP_Q"
    multipoles::UN{Multipoles}, "PP_MULTIPOLES"
    qfcoeff::UN{Qfcoeff}, "PP_QFCOEFF"
    rinner::UN{Rinner}, "PP_RINNER"
    qijl::Vector{Qijl}, "PP_QIJL"
end

@aml struct Nonlocal "PP_NONLOCAL"
    beta::Vector{Beta}, "PP_BETA"
    dij::Dij, "PP_DIJ"
    augmentation::UN{Augmentation}, "PP_AUGMENTATION"
end

@aml struct Rhoatom "PP_RHOATOM"
    text::String, txt""
end

@aml struct Aewfc "PP_AEWFC"
    l::UInt, att"l"
    index::UN{Int}, att"index"
    label::UN{String}, att"label"
    text::String, txt""
end

@aml struct FullWfc "PP_FULL_WFC"
    aewfc::Vector{Aewfc}, "PP_AEWFC"
end

@aml struct UPF <: Pseudopotential doc"UPF"
    version::VersionNumber, att"version"
    info::Info, "PP_INFO"
    header::Header, "PP_HEADER"
    mesh::Mesh, "PP_MESH", validate
    nlcc::UN{Nlcc}, "PP_NLCC"
    loc::Local, "PP_LOCAL"
    nonlocal::Nonlocal, "PP_NONLOCAL"
    semilocal::UN{Semilocal}, "PP_SEMILOCAL"
    pswfc::UN{Pswfc}, "PP_PSWFC"
    full_wfc::UN{FullWfc}, "PP_FULL_WFC"
    rhoatom::Rhoatom, "PP_RHOATOM"
    # paw::UN{Paw}, "PP_PAW"
end

function validate(x::Mesh)
    r, rab = getdata.((x.r, x.rab))
    return x.mesh === nothing ? size(r) == size(rab) : x.mesh == length(r) == length(rab)
end

function fixenumeration!(doc, name)
    children = findall("//*[contains(name(), '$name')]", doc)  # See https://stackoverflow.com/a/40124534/3260253
    if !isempty(children)  # Need to change something
        for child in children
            setnodename!(child, name)
        end
    end
    return doc
end

function Base.parse(::Type{UPF}, str)
    doc = parsexml(str)
    fixenumeration!(doc, "PP_CHI")
    fixenumeration!(doc, "PP_BETA")
    fixenumeration!(doc, "PP_AEWFC")
    fixenumeration!(doc, "PP_QIJL")
    return UPF(doc)
end

# Idea from https://github.com/JuliaData/CSV.jl/blob/c3af297/src/CSV.jl#L64-L69
function Base.read(io::IO, ::Type{UPF})
    str = read(io, String)
    return parse(UPF, str)
end
function Base.read(filename::AbstractString, ::Type{UPF})
    str = read(filename, String)
    return parse(UPF, str)
end

getdata(x::Union{Rhoatom,Nlcc,Local,R,Rab,Chi,Beta,Dij,Q,Multipoles,Qijl}) =
    parsevec(x.text)

function Base.getproperty(x::Header, name::Symbol)
    if name in (
        :is_ultrasoft,
        :is_paw,
        :is_coulomb,
        :has_so,
        :has_wfc,
        :has_gipaw,
        :paw_as_gipaw,
        :core_correction,
    )
        return istrue(getfield(x, name))
    else
        return getfield(x, name)
    end
end

# From https://github.com/mauro3/Parameters.jl/blob/ecbf8df/src/Parameters.jl#L554-L561
function Base.show(io::IO, x::UPF)
    if get(io, :compact, false) || get(io, :typeinfo, nothing) == typeof(x)
        Base.show_default(IOContext(io, :limit => true), x)
    else
        # just dumping seems to give ok output, in particular for big data-sets:
        dump(IOContext(io, :limit => true), x, maxdepth = 1)
    end
end


text:
mutable struct Easy
    handle   :: Ptr{Cvoid}
    input    :: Union{Vector{UInt8},Nothing}
    ready    :: Threads.Event
    seeker   :: Union{Function,Nothing}
    output   :: Channel{Vector{UInt8}}
    progress :: Channel{NTuple{4,Int}}
    req_hdrs :: Ptr{curl_slist_t}
    res_hdrs :: Vector{String}
    code     :: CURLcode
    errbuf   :: Vector{UInt8}
end

const EMPTY_BYTE_VECTOR = UInt8[]

function Easy()
    easy = Easy(
        curl_easy_init(),
        EMPTY_BYTE_VECTOR,
        Threads.Event(),
        nothing,
        Channel{Vector{UInt8}}(Inf),
        Channel{NTuple{4,Int}}(Inf),
        C_NULL,
        String[],
        typemax(CURLcode),
        zeros(UInt8, CURL_ERROR_SIZE),
    )
    finalizer(done!, easy)
    add_callbacks(easy)
    set_defaults(easy)
    return easy
end

function done!(easy::Easy)
    easy.handle == C_NULL && return
    curl_easy_cleanup(easy.handle)
    curl_slist_free_all(easy.req_hdrs)
    easy.handle = C_NULL
    return
end

# request options

function set_defaults(easy::Easy)
    # curl options
    setopt(easy, CURLOPT_NOSIGNAL, true)
    setopt(easy, CURLOPT_FOLLOWLOCATION, true)
    setopt(easy, CURLOPT_MAXREDIRS, 50)
    setopt(easy, CURLOPT_POSTREDIR, CURL_REDIR_POST_ALL)
    setopt(easy, CURLOPT_USERAGENT, USER_AGENT)
    # ssh-related options
    setopt(easy, CURLOPT_SSH_PRIVATE_KEYFILE, ssh_key_path())
    setopt(easy, CURLOPT_SSH_PUBLIC_KEYFILE, ssh_pub_key_path())
    key_pass = something(ssh_key_pass(), C_NULL)
    setopt(easy, CURLOPT_KEYPASSWD, ssh_pub_key_path())
end

function set_ca_roots_path(easy::Easy, path::AbstractString)
    Base.unsafe_convert(Cstring, path) # error checking
    opt = isdir(path) ? CURLOPT_CAPATH : CURLOPT_CAINFO
    setopt(easy, opt, path)
end

function set_url(easy::Easy, url::Union{String, SubString{String}})
    # TODO: ideally, Clang would generate Cstring signatures
    Base.unsafe_convert(Cstring, url) # error checking
    setopt(easy, CURLOPT_URL, url)
    set_ssl_verify(easy, verify_host(url, "ssl"))
    set_ssh_verify(easy, verify_host(url, "ssh"))
end
set_url(easy::Easy, url::AbstractString) = set_url(easy, String(url))

function set_ssl_verify(easy::Easy, verify::Bool)
    setopt(easy, CURLOPT_SSL_VERIFYPEER, verify)
    setopt(easy, CURLOPT_SSL_VERIFYHOST, verify*2)
end

function set_ssh_verify(easy::Easy, verify::Bool)
    if !verify
        setopt(easy, CURLOPT_SSH_KNOWNHOSTS, C_NULL)
    else
        file = ssh_known_hosts_file()
        setopt(easy, CURLOPT_SSH_KNOWNHOSTS, file)
    end
end

function set_method(easy::Easy, method::Union{String, SubString{String}})
    # TODO: ideally, Clang would generate Cstring signatures
    Base.unsafe_convert(Cstring, method) # error checking
    setopt(easy, CURLOPT_CUSTOMREQUEST, method)
end
set_method(easy::Easy, method::AbstractString) = set_method(easy, String(method))

function set_verbose(easy::Easy, verbose::Bool)
    setopt(easy, CURLOPT_VERBOSE, verbose)
end

function set_body(easy::Easy, body::Bool)
    setopt(easy, CURLOPT_NOBODY, !body)
end

function set_upload_size(easy::Easy, size::Integer)
    opt = Sys.WORD_SIZE ≥ 64 ? CURLOPT_INFILESIZE_LARGE : CURLOPT_INFILESIZE
    setopt(easy, opt, size)
end

function set_seeker(seeker::Function, easy::Easy)
    add_seek_callbacks(easy)
    easy.seeker = seeker
end

function set_timeout(easy::Easy, timeout::Real)
    timeout > 0 ||
        throw(ArgumentError("timeout must be positive, got $timeout"))
    if timeout ≤ typemax(Clong) ÷ 1000
        timeout_ms = round(Clong, timeout * 1000)
        setopt(easy, CURLOPT_TIMEOUT_MS, timeout_ms)
    else
        timeout = timeout ≤ typemax(Clong) ? round(Clong, timeout) : Clong(0)
        setopt(easy, CURLOPT_TIMEOUT, timeout)
    end
end

function add_header(easy::Easy, hdr::Union{String, SubString{String}})
    # TODO: ideally, Clang would generate Cstring signatures
    Base.unsafe_convert(Cstring, hdr) # error checking
    easy.req_hdrs = curl_slist_append(easy.req_hdrs, hdr)
    setopt(easy, CURLOPT_HTTPHEADER, easy.req_hdrs)
end

add_header(easy::Easy, hdr::AbstractString) = add_header(easy, string(hdr)::String)
add_header(easy::Easy, key::AbstractString, val::AbstractString) =
    add_header(easy, isempty(val) ? "$key;" : "$key: $val")
add_header(easy::Easy, key::AbstractString, val::Nothing) =
    add_header(easy, "$key:")
add_header(easy::Easy, pair::Pair) = add_header(easy, pair...)

function add_headers(easy::Easy, headers::Union{AbstractVector, AbstractDict})
    for hdr in headers
        hdr isa Pair{<:AbstractString, <:Union{AbstractString, Nothing}} ||
            throw(ArgumentError("invalid header: $(repr(hdr))"))
        add_header(easy, hdr)
    end
end

function enable_progress(easy::Easy, on::Bool=true)
    setopt(easy, CURLOPT_NOPROGRESS, !on)
end

function enable_upload(easy::Easy)
    add_upload_callbacks(easy::Easy)
    setopt(easy, CURLOPT_UPLOAD, true)
end

# response info

function get_protocol(easy::Easy)
    proto_ref = Ref{Clong}()
    r = @check curl_easy_getinfo(easy.handle, CURLINFO_PROTOCOL, proto_ref)
    r == CURLE_UNKNOWN_OPTION && error("The `libcurl` version you are using is too old and does not include the `CURLINFO_PROTOCOL` feature. Please upgrade or us a Julia build that uses its own `libcurl` library.")
    proto = proto_ref[]
    proto == CURLPROTO_DICT   && return "dict"
    proto == CURLPROTO_FILE   && return "file"
    proto == CURLPROTO_FTP    && return "ftp"
    proto == CURLPROTO_FTPS   && return "ftps"
    proto == CURLPROTO_GOPHER && return "gopher"
    proto == CURLPROTO_HTTP   && return "http"
    proto == CURLPROTO_HTTPS  && return "https"
    proto == CURLPROTO_IMAP   && return "imap"
    proto == CURLPROTO_IMAPS  && return "imaps"
    proto == CURLPROTO_LDAP   && return "ldap"
    proto == CURLPROTO_LDAPS  && return "ldaps"
    proto == CURLPROTO_POP3   && return "pop3"
    proto == CURLPROTO_POP3S  && return "pop3s"
    proto == CURLPROTO_RTMP   && return "rtmp"
    proto == CURLPROTO_RTMPE  && return "rtmpe"
    proto == CURLPROTO_RTMPS  && return "rtmps"
    proto == CURLPROTO_RTMPT  && return "rtmpt"
    proto == CURLPROTO_RTMPTE && return "rtmpte"
    proto == CURLPROTO_RTMPTS && return "rtmpts"
    proto == CURLPROTO_RTSP   && return "rtsp"
    proto == CURLPROTO_SCP    && return "scp"
    proto == CURLPROTO_SFTP   && return "sftp"
    proto == CURLPROTO_SMB    && return "smb"
    proto == CURLPROTO_SMBS   && return "smbs"
    proto == CURLPROTO_SMTP   && return "smtp"
    proto == CURLPROTO_SMTPS  && return "smtps"
    proto == CURLPROTO_TELNET && return "telnet"
    proto == CURLPROTO_TFTP   && return "tftp"
    return nothing
end

status_2xx_ok(status::Integer) = 200 ≤ status < 300
status_zero_ok(status::Integer) = status == 0

const PROTOCOL_STATUS = Dict{String,Function}(
    "dict"   => status_2xx_ok,
    "file"   => status_zero_ok,
    "ftp"    => status_2xx_ok,
    "ftps"   => status_2xx_ok,
    "http"   => status_2xx_ok,
    "https"  => status_2xx_ok,
    "ldap"   => status_zero_ok,
    "ldaps"  => status_zero_ok,
    "pop3"   => status_2xx_ok,
    "pop3s"  => status_2xx_ok,
    "rtsp"   => status_2xx_ok,
    "scp"    => status_zero_ok,
    "sftp"   => status_2xx_ok,
    "smtp"   => status_2xx_ok,
    "smtps"  => status_2xx_ok,
)

function status_ok(proto::AbstractString, status::Integer)
    test = get(PROTOCOL_STATUS, proto, nothing)
    test !== nothing && return test(status)::Bool
    error("Downloads.jl doesn't know the correct request success criterion for $proto: you can use `request` and check the `status` field yourself or open an issue with Downloads with details an example URL that you are trying to download.")
end
status_ok(proto::Nothing, status::Integer) = false

function get_effective_url(easy::Easy)
    url_ref = Ref{Ptr{Cchar}}()
    @check curl_easy_getinfo(easy.handle, CURLINFO_EFFECTIVE_URL, url_ref)
    return unsafe_string(url_ref[])
end

function get_response_status(easy::Easy)
    code_ref = Ref{Clong}()
    @check curl_easy_getinfo(easy.handle, CURLINFO_RESPONSE_CODE, code_ref)
    return Int(code_ref[])
end

function get_response_info(easy::Easy)
    proto = get_protocol(easy)
    url = get_effective_url(easy)
    status = get_response_status(easy)
    message = ""
    headers = Pair{String,String}[]
    if proto in ("http", "https")
        message = isempty(easy.res_hdrs) ? "" : easy.res_hdrs[1]
        for hdr in easy.res_hdrs
            if contains(hdr, r"^\s*$")
                # ignore
            elseif (m = match(r"^(HTTP/\d+(?:.\d+)?\s+\d+\b.*?)\s*$", hdr)) !== nothing
                message = m.captures[1]
                empty!(headers)
            elseif (m = match(r"^(\S[^:]*?)\s*:\s*(.*?)\s*$", hdr)) !== nothing
                push!(headers, lowercase(m.captures[1]) => m.captures[2])
            else
                @warn "malformed HTTP header" url status header=hdr
            end
        end
    elseif proto in ("ftp", "ftps", "sftp")
        message = isempty(easy.res_hdrs) ? "" : easy.res_hdrs[end]
    else
        # TODO: parse headers of other protocols
    end
    message = chomp(message)
    endswith(message, '.') && (message = chop(message))
    return proto, url, status, message, headers
end

function get_curl_errstr(easy::Easy)
    easy.code == Curl.CURLE_OK && return ""
    errstr = easy.errbuf[1] == 0 ?
        unsafe_string(Curl.curl_easy_strerror(easy.code)) :
        GC.@preserve easy unsafe_string(pointer(easy.errbuf))
    return chomp(errstr)
end

# callbacks

function header_callback(
    data   :: Ptr{Cchar},
    size   :: Csize_t,
    count  :: Csize_t,
    easy_p :: Ptr{Cvoid},
)::Csize_t
    easy = unsafe_pointer_to_objref(easy_p)::Easy
    n = size * count
    hdr = unsafe_string(data, n)
    push!(easy.res_hdrs, hdr)
    return n
end

# feed data to read_callback
function upload_data(easy::Easy, input::IO)
    while true
        data = eof(input) ? nothing : readavailable(input)
        easy.input === nothing && break
        easy.input = data
        curl_easy_pause(easy.handle, Curl.CURLPAUSE_CONT)
        wait(easy.ready)
        easy.input === nothing && break
        easy.ready = Threads.Event()
    end
end

function read_callback(
    data   :: Ptr{Cchar},
    size   :: Csize_t,
    count  :: Csize_t,
    easy_p :: Ptr{Cvoid},
)::Csize_t
    easy = unsafe_pointer_to_objref(easy_p)::Easy
    buf = easy.input
    if buf === nothing
        notify(easy.ready)
        return 0 # done uploading
    end
    if isempty(buf)
        notify(easy.ready)
        return CURL_READFUNC_PAUSE # wait for more data
    end
    n = min(size * count, length(buf))
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), data, buf, n)
    deleteat!(buf, 1:n)
    return n
end

function seek_callback(
    easy_p :: Ptr{Cvoid},
    offset :: curl_off_t,
    origin :: Cint,
)::Cint
    if origin != 0
        @async @error("seek_callback: unsupported seek origin", origin)
        return CURL_SEEKFUNC_CANTSEEK
    end
    easy = unsafe_pointer_to_objref(easy_p)::Easy
    easy.seeker === nothing && return CURL_SEEKFUNC_CANTSEEK
    try easy.seeker(offset)
    catch err
        @async @error("seek_callback: seeker failed", err)
        return CURL_SEEKFUNC_FAIL
    end
    return CURL_SEEKFUNC_OK
end

function write_callback(
    data   :: Ptr{Cchar},
    size   :: Csize_t,
    count  :: Csize_t,
    easy_p :: Ptr{Cvoid},
)::Csize_t
    easy = unsafe_pointer_to_objref(easy_p)::Easy
    n = size * count
    buf = Array{UInt8}(undef, n)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), buf, data, n)
    put!(easy.output, buf)
    return n
end

function progress_callback(
    easy_p   :: Ptr{Cvoid},
    dl_total :: curl_off_t,
    dl_now   :: curl_off_t,
    ul_total :: curl_off_t,
    ul_now   :: curl_off_t,
)::Cint
    easy = unsafe_pointer_to_objref(easy_p)::Easy
    put!(easy.progress, (dl_total, dl_now, ul_total, ul_now))
    return 0
end

function add_callbacks(easy::Easy)
    # pointer to easy object
    easy_p = pointer_from_objref(easy)
    setopt(easy, CURLOPT_PRIVATE, easy_p)

    # pointer to error buffer
    errbuf_p = pointer(easy.errbuf)
    setopt(easy, CURLOPT_ERRORBUFFER, errbuf_p)

    # set header callback
    header_cb = @cfunction(header_callback,
        Csize_t, (Ptr{Cchar}, Csize_t, Csize_t, Ptr{Cvoid}))
    setopt(easy, CURLOPT_HEADERFUNCTION, header_cb)
    setopt(easy, CURLOPT_HEADERDATA, easy_p)

    # set write callback
    write_cb = @cfunction(write_callback,
        Csize_t, (Ptr{Cchar}, Csize_t, Csize_t, Ptr{Cvoid}))
    setopt(easy, CURLOPT_WRITEFUNCTION, write_cb)
    setopt(easy, CURLOPT_WRITEDATA, easy_p)

    # set progress callback
    progress_cb = @cfunction(progress_callback,
        Cint, (Ptr{Cvoid}, curl_off_t, curl_off_t, curl_off_t, curl_off_t))
    setopt(easy, CURLOPT_XFERINFOFUNCTION, progress_cb)
    setopt(easy, CURLOPT_XFERINFODATA, easy_p)
end

function add_upload_callbacks(easy::Easy)
    # pointer to easy object
    easy_p = pointer_from_objref(easy)

    # set read callback
    read_cb = @cfunction(read_callback,
        Csize_t, (Ptr{Cchar}, Csize_t, Csize_t, Ptr{Cvoid}))
    setopt(easy, CURLOPT_READFUNCTION, read_cb)
    setopt(easy, CURLOPT_READDATA, easy_p)
end

function add_seek_callbacks(easy::Easy)
    # pointer to easy object
    easy_p = pointer_from_objref(easy)

    # set seek callback
    seek_cb = @cfunction(seek_callback,
        Cint, (Ptr{Cvoid}, curl_off_t, Cint))
    setopt(easy, CURLOPT_SEEKFUNCTION, seek_cb)
    setopt(easy, CURLOPT_SEEKDATA, easy_p)
end


text:
using Test
using MLStyle
using Expronicon

@testset "JLFunction pattern" begin
    f = @λ begin
        JLFunction(;head=:function, args) => args
        JLFunction(;head=:(=), name) => name
        _ => nothing
    end

    @test f(:(function foo(x) x end)) == [:x]
    @test f(:(goo(x) = x)) === :goo
    @test f(:(1 + 1)) === nothing
end

@testset "JLStruct/JLKwStruct pattern" begin
    f = @λ begin
            JLStruct(;name) => (name, )
            JLKwStruct(;name, typealias) => (name, typealias)
            _ => nothing
        end

    ex = @expr struct KwStructP
        x::Int = 1
    end

    @test f(ex) == (:KwStructP, nothing)

    ex = @expr struct StructP
        x::Int
    end

    @test f(ex) == (:StructP, )

    ex = :(foo(x) = x)
    @test f(ex) === nothing
end

@testset "JLIfElse pattern" begin
    f = @λ begin
        JLIfElse(;conds, stmts, otherwise) => (conds, stmts, otherwise)
        _ => nothing
    end

    ex = @expr if x > 1
        x + 1
    else
        nothing
    end
    conds, stmts, otherwise = f(ex)
    @test conds == Any[:(x > 1)]
    @test prettify.(stmts) == Any[:(x + 1)]

    @test_expr otherwise == quote
        nothing
    end

    ex = :(foo(x) = x)
    @test f(ex) === nothing
end

@testset "JLFor pattern" begin
    f = @λ begin
        JLFor(;vars) => vars
        _ => nothing
    end

    ex = @expr for i in 1:10, j in 1:l
        M[i, j] += 1
    end

    @test f(ex) == [:i, :j]
    @test f(:(1 + 1)) === nothing
end

@testset "patch" begin
    f = @λ begin
        Symbol(x) => x
        _ => nothing
    end

    @test f(:x) == "x"
    @test f(1) === nothing

    f = @λ begin
        GlobalRef(m, s) => (m, s)
        _ => nothing
    end

    @test f(GlobalRef(Main, :sin)) == (Main, :sin)
    @test f(:(1 + 1)) === nothing
end

@testset "codegen_match" begin
    ex = codegen_match(:x) do
        quote
            1 => true
            2 => false
            _ => nothing
        end
    end

    eval(codegen_ast(JLFunction(;name=:test_match, args=[:x], body=ex)))

    @test test_match(1) == true
    @test test_match(2) == false
    @test test_match(3) === nothing
end

@testset "JLMatch" begin
    jl = JLMatch(:x)
    jl[1] = true
    jl[2] = :(sin(x))
    jl[3] = quote
        1 + 2
        3 + 4
    end
    println(jl)
    ex = codegen_ast(jl)
    jl = JLFunction(;name=:test_match, args=[:x], body=ex)
    println(jl)
    eval(codegen_ast(jl))

    @test test_match(1) == true
    @test test_match(2) == sin(2)
    @test test_match("abc") === nothing
end


text:
struct Acoustic <: Equation end
@declare_dofs Acoustic [:u, :v, :pressure, :rho, :lambda]
struct GaussianWave <: Scenario
end

function is_periodic_boundary(equation::Acoustic, scenario::GaussianWave)
    # GaussianWave scenario doesn't require periodic boundary conditions
    false
end

function evaluate_boundary(eq::Acoustic, scenario::GaussianWave, face, normalidx, dofsface, dofsfaceneigh)
# dofsface and dofsfaceneigh have shape (num_2d_quadpoints, dofs)
# you need to set dofsfaceneigh
    dofsfaceneigh .= dofsface
    if normalidx == 1
        dofsfaceneigh[:, 1] = -dofsface[:, 1]
    else
        dofsfaceneigh[:, 2] = -dofsface[:, 2]
    end
end

function get_initial_values(eq::Acoustic, scenario::GaussianWave, global_position; t=0.0)
    x, y = global_position
    if is_analytical_solution(eq, scenario)
        println("No analytical solution implemented")
    else
        [ 0, 0, exp(-100*(x - 0.5)^2 - 100*(y - 0.5)^2), 1.0, if x<= 0.5 1/5 else 1 end]
    end
end

function is_analytical_solution(equation::Acoustic, scenario::GaussianWave)
    false
end

function evaluate_flux(eq::Acoustic, celldofs, cellflux)
    # size(cellflux)=(2*order*order, ndof)=(2,5)
    # size(celldofs)=(order*order, ndof)=(1,5)

    order_sq, ndof = size(celldofs)
    for j=1:order_sq
        # explicit names of dofs for clarity
        vx = celldofs[j, 1]
        vy = celldofs[j, 2]
        pressure = celldofs[j, 3]
        density = celldofs[j, 4]
        bulk = celldofs[j, 5]

        # here flux is computed elementwise
        cellflux[j,1]= pressure / density
        cellflux[j + order_sq, 1] = 0

        cellflux[j,2] = 0
        cellflux[j + order_sq, 2] = pressure / density

        cellflux[j,3] = bulk*vx
        cellflux[j + order_sq, 3] = bulk*vy

        cellflux[j,4] = 0
        cellflux[j + order_sq, 4] = 0

        cellflux[j,5] = 0
        cellflux[j + order_sq, 5] = 0

    end

end

function max_eigenval(eq::Acoustic, celldata, normalidx)
    # sqrt( K_0 / rho_0 )
    maxeigenval = 0.
    for i in 1:size(celldata)[1]
        vx, vy, pressure, density, bulk = celldata[i, :]
        maxeigenval_new = sqrt(bulk / density)
        if maxeigenval_new > maxeigenval
            maxeigenval = maxeigenval_new
        end
    end
    maxeigenval
end


text:
# Use `perl ../deps/imports.pl base.jl` to generate the next line
import Base: !=, !==, *, +, -, /, <, <=, ==, >, >=, \, ^, abs, abs2, all, any, axes, big, broadcast, ceil, checkbounds, copy, count, div, eachindex, eltype, eps, float, floor, identity, isassigned, isempty, isequal, isfinite, isinf, isinteger, isless, isnan, lastindex, length, literal_pow, maximum, minimum, ndims, oftype, one, ones, permutedims, prod, rem, reshape, round, sign, signbit, similar, size, stride, strides, sum, trunc, typemax, typemin, unsafe_trunc, values, vec, widemul, zero
import Base.Broadcast: broadcasted

# The following list copied from relevant portions of julia/base/exports.jl

### Operators
# !   Bool function
@zerograd !=(x1,x2)
# ≠   Same as !=
@zerograd !==(x1,x2)
# ≡   Same as ===, Core builtin function, cannot add methods
# ≢   Same as !==
# xor Int function
# ⊻   Same as xor
# %   Same as rem
# ÷   Same as div
# &   Int function

@primitive +(x),dy  dy
@primitive +(x1,x2),dy  unbroadcast(x1,dy)  unbroadcast(x2,dy)
@primitive -(x),dy  -dy
@primitive -(x1,x2),dy  unbroadcast(x1,dy)  unbroadcast(x2,-dy)

@primitive *(x),dy                      dy
@primitive1 *(x1,x2),dy                 (dy*x2')  (x1'*dy)
@primitive2 *(x1,x2),dy                 unbroadcast(x1,dy.*x2)  unbroadcast(x2,x1.*dy)
@primitive *(x1::Number,x2::Number),dy  unbroadcast(x1,dy.*x2)  unbroadcast(x2,x1.*dy)
@primitive *(x1::Number,x2),dy          unbroadcast(x1,dy.*x2)  unbroadcast(x2,x1.*dy)
@primitive *(x1,x2::Number),dy          unbroadcast(x1,dy.*x2)  unbroadcast(x2,x1.*dy)

@primitive1 /(x1,x2),dy                 error("A/B grad undefined") error("A/B grad undefined") #TODO
@primitive2 /(x1,x2),dy                 unbroadcast(x1,dy./x2)  unbroadcast(x2,-dy.*x1./abs2.(x2))
@primitive /(x1::Number,x2::Number),dy  unbroadcast(x1,dy./x2)  unbroadcast(x2,-dy.*x1./abs2.(x2))
@primitive /(x1::Number,x2),dy          unbroadcast(x1,dy./x2)  unbroadcast(x2,-dy.*x1./abs2.(x2))
@primitive /(x1,x2::Number),dy          unbroadcast(x1,dy./x2)  unbroadcast(x2,-dy.*x1./abs2.(x2))

@primitive1 \(x1,x2),dy                 error("A\\B grad undefined") error("A\\B grad undefined") #TODO
@primitive2 \(x1,x2),dy                 unbroadcast(x1,-dy.*x2./abs2.(x1))  unbroadcast(x2,dy./x1)
@primitive \(x1::Number,x2::Number),dy  unbroadcast(x1,-dy.*x2./abs2.(x1))  unbroadcast(x2,dy./x1)
@primitive \(x1::Number,x2),dy          unbroadcast(x1,-dy.*x2./abs2.(x1))  unbroadcast(x2,dy./x1)
@primitive \(x1,x2::Number),dy          unbroadcast(x1,-dy.*x2./abs2.(x1))  unbroadcast(x2,dy./x1)

@primitive1 ^(x1,x2),dy,y                  error("A^B grad undefined") error("A^B grad undefined") #TODO: Arr^Int (Arr^Num can be imaginary)
@primitive2 ^(x1,x2),dy,y                  unbroadcast(x1,dxndx(x1,x2,dy))  unbroadcast(x2,dy.*y.*log.(x1))
@primitive ^(x1::Number,x2::Number),dy,y   dxndx(x1,x2,dy)  dy*y*log(x1)
literal_pow(::typeof(^), x::Value, ::Val{N}) where N = forw(^,x,N) # x^p for any literal integer p is lowered to Base.literal_pow(^, x, Val(p))
broadcasted(::typeof(literal_pow), ::typeof(^), x::Value, ::Val{N}) where N = forw(broadcasted,^,x,N)
dxndx(x1,x2,dy)=(if x2==0; zero(dy); elseif x2==1; dy; elseif x2==2; 2 .* x1 .* dy; else; dy.*x2.*x1.^(x2 .- 1); end) # optimize common cases

# //  Int function
@zerograd <(x1,x2)
# <:  Type function
# <<  Int function
@zerograd <=(x1,x2)
# ≤   Same as <=
@zerograd ==(x1,x2)
@zerograd >(x1,x2)
# >:  Type function
@zerograd >=(x1,x2)
# ≥   Same as >=
# >>  Int function
# >>> Int function

# |   Int function
# |>  Function chaining
# ~   Int function
# :   Range operator
# =>  Pair constructor
# ∘   Function composition

### scalar math (some of this is in math.jl)
@primitive abs(x),dy (dy.*sign.(x))
@primitive abs2(x),dy (dy.*(2x))
# angle(z): Compute the phase angle in radians of a complex number z.
@primitive big(x),dy oftype(x,dy)
# binomial(n,k): Number of ways to choose k out of n items.
# bswap(n): Reverse the byte order of n.
@zerograd ceil(x)
# cis(z): Return \exp(iz).
# cld(x, y): Smallest integer larger than or equal to x/y.
# cmp(x,y): Return -1, 0, or 1 depending on whether x is less than, equal to, or greater than y, respectively.
# complex(r, [i]): Convert real numbers or arrays to complex. i defaults to zero.
# conj(z): Compute the complex conjugate of a complex number z.
# copysign(x, y) -> z: Return z which has the magnitude of x and the same sign as y. (TODO)
# count_ones(x::Integer): Number of ones in the binary representation of x.
# count_zeros(x::Integer): Number of zeros in the binary representation of x.
# denominator(x): Denominator of the rational representation of x.
@zerograd div(x,y)
# divrem(x,y): The quotient and remainder from Euclidean division. Equivalent to (div(x,y), rem(x,y)) or (x÷y, x%y).
@zerograd eps(x)
# factorial,
# fld,
# fld1,
# fldmod,
# fldmod1,
# flipsign,
@primitive float(x),dy dy
# tryparse,
@zerograd floor(x)
# fma,
# gcd,
# gcdx,
# imag,
# inv,
# invmod,
# isapprox,
# iseven,
@zerograd isfinite(x)
@zerograd isinf(x)
@zerograd isinteger(x)
@zerograd isnan(x)
# isodd,
# ispow2,
# isqrt,
# isreal,
# issubnormal,
# iszero,
# isone,
# lcm,
# leading_ones,
# leading_zeros,
# maxintfloat,
# mod,
# mod1,
# muladd,
# nextfloat,
# nextpow,
# nextpow2,
# nextprod,
# numerator,
@zerograd one(x)
# oneunit,
# powermod,
# prevfloat,
# prevpow,
# prevpow2,
# rationalize,
# real,
# realmax,
# realmin,
# reim,
# reinterpret,
@primitive rem(x1,x2),dy,y  unbroadcast(x1,dy)  unbroadcast(x2,-dy.*div.(x1,x2))
@zerograd round(x)
@zerograd sign(x)
@zerograd signbit(x)
# signed: Int function
# trailing_ones:  Int function
# trailing_zeros: Int function
@zerograd trunc(x)
@zerograd unsafe_trunc(x)
@zerograd typemax(x)
@zerograd typemin(x)
# unsigned: Int function
@primitive widemul(x1,x2),dy  unbroadcast(x1,dy.*x2)  unbroadcast(x2,x1.*dy)
@zerograd zero(x)
# √  Same as sqrt
# ∛  Same as cbrt
# ≈  Same as isapprox
# ≉  This is equivalent to !isapprox(x,y) (see isapprox).

### arrays
@zerograd axes(x,i...)
# broadcast!,
# broadcast  Handled in broadcast.jl
# cat  Handled in cat.jl
@zerograd checkbounds(x,i...)
# checkindex,
# circcopy!,
# circshift,
# circshift!,
# clamp!,
# conj!,
# copy!,
# copyto!,
# diff,
# cumprod,
# cumprod!,
# cumsum,
# cumsum!,
# accumulate,
# accumulate!,
@zerograd eachindex(x,i...)
# extrema,
# fill!,
# fill,
# first,
# hcat  Handled in cat.jl
# hvcat,
# indexin,
# argmax,
# argmin,
# invperm,
# invpermute!,
@zerograd isassigned(x,i...)
# isperm,
# issorted,
# last,
# mapslices,
# maximum!,
@primitive maximum(x;d...),dy,y  (dy.*(y.==x))
@primitive maximum(f::typeof(abs),x;d...),dy,y  nothing  (dy.*(y.==abs.(x)).*sign.(x))
# minimum!,
@primitive minimum(x;d...),dy,y  (dy.*(y.==x))
@primitive minimum(f::typeof(abs),x;d...),dy,y  nothing  (dy.*(y.==abs.(x)).*sign.(x))
# minmax,
@zerograd ndims(x)
@zerograd ones(x)
# parent,
# parentindices,
# partialsort,
# partialsort!,
# partialsortperm,
# partialsortperm!,
# permute!,
@primitive permutedims(x,d),dy  permutedims(dy,invperm(d))
@primitive permutedims(x),dy  reshape(permutedims(dy),size(x)) # need reshape for vectors
# permutedims!,
# prod!,
@primitive prod(x;d...),dy,y  (dy.*(y./x))  # TODO: prod with abs, abs2
# promote_shape,
# range,
@primitive reshape(x,i...),dy  reshape(dy,size(x))
# reverse!,
# reverse,
# rot180,
# rotl90,
# rotr90,
# shuffle,
# shuffle!,
@zerograd size(x,i...)
# selectdim, # implemented in getindex.jl
# sort!,
# sort,
# sortcols,
# sortperm,
# sortperm!,
# sortrows,
# squeeze,
# step,
@zerograd stride(x,i...)
@zerograd strides(x)
# sum!,
@primitive sum(x;d...),dy  (dy.*one.(x))
@primitive sum(f::typeof(abs),x;d...),dy   nothing  (dy.*sign.(x))
@primitive sum(f::typeof(abs2),x;d...),dy  nothing  (dy.*(2x))
# to_indices,
# vcat  Handled in cat.jl
@primitive vec(x),dy  reshape(dy,size(x))
# view, # implemented in getindex.jl
# zeros: no longer called with array arguments, using zero instead.

# collections
# all!,
@zerograd all(a;dims=:)
@zerograd all(f,a;dims=:)
# allunique,
# any!,
@zerograd any(a;dims=:)
@zerograd any(f,a;dims=:)
# firstindex,
# collect,
@zerograd count(a;dims=:)
@zerograd count(f,a;dims=:)
# delete!,
# deleteat!,
@zerograd eltype(x)
# empty!,
# empty,
@zerograd lastindex(x,d...)
# filter!,
# filter,
# foldl,
# foldr,
# foreach,
# get  Handled in getindex.jl
# get!,
# getindex  Handled in getindex.jl
# getkey,
# haskey,
# in,
# intersect!,
# intersect,
@zerograd isempty(x)
# issubset,
# issetequal,
# keys,
# keytype,
@zerograd length(x)
# map!,
# map,
# mapfoldl,
# mapfoldr,
# mapreduce,
# merge!,
# merge,
# pairs,
# reduce,
# setdiff!,
# setdiff,
# setindex!,
@zerograd similar(x,i...)
similar(x::Value, dims::Base.DimOrInd...) = similar(value(x), dims...) # 0.7 ambiguity fix
# sizehint!,
# splice!,
# symdiff!,
# symdiff,
# union!,
# union,
# unique!,
# unique,
values(a::Value{T})  where {T<:AbstractDict} = (a[k] for k in keys(value(a)))
# valtype,
# ∈,
# ∉,
# ∋,
# ∌,
# ⊆,
# ⊈,
# ⊊,
# ⊇,
# ⊉,
# ⊋,
# ∩,
# ∪,

### iteration
# done  deprecated
# next  deprecated
# start deprecated
# iterate    Handled in iterate.jl
# enumerate  Implemented with iterate
# zip        Implemented with iterate

### object identity and equality
@primitive copy(x),dy dy
# deepcopy,
# hash,
@primitive identity(x),dy dy
# isbits,
@zerograd isequal(x,y)
# isimmutable,
@zerograd isless(x,y)
# ifelse,
# objectid,
# sizeof,

### types
# convert,
# getproperty,
# setproperty!,
# fieldoffset,
# fieldname,
# fieldnames,
# fieldcount,
# propertynames,
# isabstracttype,
# isbitstype,
# isprimitivetype,
# isstructtype,
# isconcretetype,
# isdispatchtuple,
@primitive oftype(t,x),dy nothing oftype(x,dy)
# promote,
# promote_rule,
# promote_type,
# instances,
# supertype,
# typeintersect,
# typejoin,
# widen,


text:
function ueberpruefeAbleitung(Ableitung,A2,Messwerte)
Arows,Acols = size(A2);
f_nah = zeros(Acols,1).+1;
a =
if round.(Ableitung(f_nah);digits=2) ==  round.(2*transpose(A2)*A2*f_nah-2*transpose(A2)*Messwerte;digits=2)
    println("Eure Ableitung ist korrekt!");
else
    println("Eure Ableitung stimmt noch nicht. Stellt sicher, dass ihr nirgends das Transponieren vergessen habt.");
end

end


text:
using WaveSimulator
using ComputationalResources
using Base.Test

if Base.Threads.nthreads() > 1
    addresource(CPUThreads)
end
if CUDAdrv.vendor() != "none"
    addresource(CUDALibs)
end
info("Testing resources: $(ComputationalResources.resources)")

# Generate resource instances
test_resources = AbstractResource[]
for resource_type in ComputationalResources.resources
    if resource_type == CPU1
        push!(test_resources, CPU1())
    elseif resource_type == CPUThreads
        push!(test_resources, CPUThreads((64,64,1)))
    elseif resource_type == CUDALibs
        push!(test_resources, CUDALibs())
    end
end

@testset "Backend Simulation State" begin
    test_results = Dict()
    for resource in test_resources
        # Setup
        wave = UniformWave{3}(fmax=2e3)
        sim = Simulator(wave, resource=resource, duration=0.01)
        domain = BoxDomain(6,8,4, gamma=0.05)
        f0 = WaveSimulator.gauss(domain)
        backend = WaveSimulator.backend_init(sim.resource, domain, sim)
        state = WaveSimulator.state_init(f0, backend, domain, sim)
        # Simulation
        for i in 1:5
            simulate!(state, backend, sim)
        end
        test_results[typeof(resource)] = state.current
    end

    # Test that all considered backends return similar results
    cpu1_result = first(values(filter((res,_)->res<:CPU1, test_results)))
    filter!((res,result)->!(res isa CPU1), test_results)
    for resource in keys(test_results)
        info("Comparing CPU1 to $(resource)")
        @test isapprox(cpu1_result, Array(test_results[resource]))
    end
end