text:
@with_kw struct CSGGeneticBuildParameters{R<:Real} @deftype R
    # scoring
    ϵ_d = 0.01
    α = deg2rad(10)
    # mutation and crossover
    χ = 0.3
    μ = 0.3
    μ_0 = 0.3
    # tree & genetic iteration parameters
    maxdepth::Int = 10
    itermax::Int = 3000
    populationsize::Int = 150
    keepbestn::Int = 2
    # tournament selection
    tournamentsize::Int = 30
    selectionprob = 0.5
    rungc::Bool = false
end

function crossover(creatures, params)
    @unpack χ, maxdepth = params
    if rand() >= χ
        #@debug "rand() >= χ"
        return creatures
    end
    creat_copy = deepcopy(creatures)
    podfs = PreOrderDFS.(creat_copy)
    nofnodes = numberofnodes.(podfs)
    #@debug "nofnodes: $nofnodes"
    (isempty(1:nofnodes[1]) || isempty(1:nofnodes[2])) && return creatures
    randnodes = map(x->rand(1:x), nofnodes)
    #@debug "randnodes: $randnodes"
    changenodes = reverse(deepcopy(selectfirstchildnode.(podfs, randnodes)))
    for i in 1:2
        n = 0
        for j in podfs[i]
            isempty(j.children) && continue
            n += 1
            if n == randnodes[i]
                j.children[1] = changenodes[i]
                break
            end
        end
    end
    for c in creat_copy
        if depth(c) > maxdepth
            #@debug "depth(c) > maxdepth"
            return creatures
        end
    end
    return creat_copy
end

function mutatecached(creature, nodes, params)
    @unpack μ, μ_0, maxdepth = params

    if rand() < μ
        # return mutated
        if rand() < μ_0
            # return new tree
            return randomcachedtree(nodes, maxdepth)
        else
            # mutate tree
            creat_copy = deepcopy(creature)
            podfs = PreOrderDFS(creat_copy)
            nofnodes = numberofnodes(podfs)
            isempty(1:nofnodes) && return creature
            randnode = rand(1:nofnodes)
            selected_node = selectfirstchildnode(podfs, randnode)
            nd = depth(selected_node)
            newdepth = nd > 0 ? nd : 1
            newnode = randomcachedtree(nodes, newdepth)
            n = 0
            for j in podfs
                isempty(j.children) && continue
                n += 1
                if n == randnode
                    j.children[1] = newnode
                    break
                end
            end
            if depth(creat_copy) > maxdepth
                return creature
            end
            return creat_copy
        end
    else
        # don't mutate tree
        return creature
    end
end

function tournament(params)
    # asserts that the population is sorted in descending order
    @unpack tournamentsize, selectionprob, populationsize = params
    tourn = [rand(1:populationsize) for _ in 1:tournamentsize]
    sort!(tourn)
    for ind in tourn
        rand() < selectionprob && return ind
    end
    return tourn[1]
end

function projectnormal(res::CachedResult, surfaces, n)
    pl = surfaces[res.index]
    cf = pl.coordframe
    n_ = SVector{2,Float64}(dot(cf[1],n), dot(cf[2],n))
    resn = res.signint*n_
    return resn
end

function rawscorefunc(tree, surfaces, cpoints, cnormals, normals, params)
    @unpack ϵ_d, α = params
    # λ should be unpacked?
    λ = log(size(cpoints, 1))
    score = 0.
    f = tree2func(tree)
    for i in eachindex(cpoints)
        res = Base.invokelatest(f, cpoints, i)::CachedResult
        v = value(res)
        n = normal(res, cnormals, i)
        d_i = v/ϵ_d
        #TODO: this is not really safe :D
        if size(n,1) == 2
            n_comp = projectnormal(res, surfaces, normals[i])
            ddot = dot(n, n_comp)
        else
            ddot = dot(n, normals[i])
        end
        ddot = clamp(ddot, -1, 1)
        θ_i = acos(ddot)/α
        score += exp(-d_i^2) + exp(-θ_i^2)
    end
    return score - λ*treesize(tree)
end

function rankcachedpopulationfunc(population, surfaces, cpoints, cnormals, normals, params)
    sum = 0.
    score = Array{Float64,1}(undef, size(population))
    normed = similar(score)

    #smphr = Semaphore(1)
    #@threads for i in eachindex(population)
    for i in eachindex(population)
        score[i] = rawscorefunc(population[i], surfaces, cpoints, cnormals, normals, params)
    end
    #=
    for i in eachindex(score)
        normed[i] = 1/(1+score[i])
        sum += normed[i]
    end
    for i in eachindex(score)
        score[i] = normed[i]/sum
    end
    # this is the rescaled score
    =#
    p = sortperm(score, rev=true)
    bestscore = score[p][1]
    return population[p], bestscore
end

function cachedfuncgeneticbuildtree(surfaces, points, normals, params)
    @unpack itermax, maxdepth, populationsize = params
    @unpack keepbestn, rungc = params
    notifit = itermax > 10 ? div(itermax,10) : 1
    cnodes, cvalues, cnormals = cachenodes(surfaces, points)
    population = [randomcachedtree(cnodes, maxdepth) for _ in 1:populationsize]
    npopulation = similar(population)
    scores = zeros(itermax)
    start_time = time_ns()
    @info "Iteration in progress..."
    for i in 1:itermax
        population, nsc = rankcachedpopulationfunc(population, surfaces, cvalues, cnormals, normals, params)
        if i%notifit == 0
            @info "$i-th iteration - best score: $nsc"
            if rungc
                GC.gc()
            end
        end
        scores[i] = nsc
        # save the best
        npopulation[1:keepbestn] = population[1:keepbestn]
        n = keepbestn+1
        while true
            # select 2 creatures by tournament
            t1 = deepcopy(population[tournament(params)])
            t2 = deepcopy(population[tournament(params)])
            # crossover them
            newc = crossover([t1, t2], params)
            # mutate first
            cv = mutatecached(newc[1], cnodes, params)
            npopulation[n] = cv
            n += 1
            # finish if we have enough
            n > populationsize && break
            # mutate second
            cv = mutatecached(newc[2], cnodes, params)
            npopulation[n] = cv
            n += 1
            n > populationsize && break
        end
        # npopulation is populated with a new set of creature
        population = npopulation
    end
    fint = trunc((time_ns() - start_time)/1_000_000_000, digits=2)
    @info "Finished: $itermax iteration in $fint seconds."
    return population, cached2normaltree(population[1], surfaces), fint, scores
end


text:
# =======================================================
#         CycleGraph struct
# =======================================================

"""
    CycleGraph <: AbstractGraph

A structure representing an undirected cycle graph.

A `CycleGraph` with one vertex is a single vertex without any edges (no self-loops)
and a `CycleGraph` with two vertices is a single edge.

See also: [`LightGraphs.cycle_graph`](@ref)
"""
struct CycleGraph{T<:Integer} <: AbstractGraph{T}
    nv::T

    function CycleGraph{T}(nv) where {T}

        nv = convert(T, nv)
        nv >= zero(T) || throw(ArgumentError("nv must be >= 0"))

        return new{T}(nv)
    end
end

CycleGraph(nv::T) where {T<: Integer} = CycleGraph{T}(nv)

# =======================================================
#         traits
# =======================================================

LG.is_directed(::Type{<:CycleGraph}) = false

#
# =======================================================
#        vertices
# =======================================================

Base.eltype(::Type{CycleGraph{T}}) where {T} = T
Base.eltype(g::CycleGraph) = eltype(typeof(g))

LG.edgetype(::Type{CycleGraph{T}}) where {T} = Edge{T}

LG.nv(g::CycleGraph) = g.nv

LG.vertices(g::CycleGraph{T}) where {T} = Base.OneTo(nv(g))

LG.has_vertex(g::CycleGraph, v) = v in vertices(g)

# =======================================================
#        edges
# =======================================================

LG.edgetype(g::CycleGraph{T}) where {T} = edgetype(typeof(g))

function LG.ne(g::CycleGraph)
    nvg = Int(nv(g))

    nvg >= 3 && return nvg
    nvg == 2 && return 1

    return 0
end

function LG.has_edge(g::CycleGraph{T}, u, v) where {T}

    u, v = minmax(u, v)
    nvg = nv(g)
    oneT = one(T)
    isinbounds = (oneT <= u) & (v <= nvg)
    isedge = ((v - u == oneT) | (v - u == nvg - oneT)) & (u != v)
    return isinbounds & isedge
end

# ---- edges vector -----------------------------------

LG.edges(g::CycleGraph) = SimpleEdgeVector(g)

function Base.size(edgevec::SimpleEdgeVector{V, G}) where {V, G <: CycleGraph}

    g = edgevec.graph

    return (ne(g), )
end

@inline function Base.getindex(edgevec::SimpleEdgeVector{V, G}, i::Int) where {V, G <: CycleGraph}

    @boundscheck i ∈ Base.OneTo(length(edgevec)) || throw(BoundsError(edgevec, i))

    g = edgevec.graph
    T = eltype(g)
    nvg::T = nv(g)

    i == 1 && return Edge(one(T), T(2))
    i == 2 && return Edge(one(T), nvg)
    return Edge(T(i - 1), T(i))
end

Base.IndexStyle(::Type{<:SimpleEdgeVector{V, G}}) where {V, G <: CycleGraph} = IndexLinear()

# =======================================================
#        neighbors
# =======================================================

@inline function LG.outneighbors(g::CycleGraph, v::Integer)

    @boundscheck v ∈ vertices(g) || throw(BoundError(g, v))

    return OutNeighborVector(g, eltype(g)(v))
end

LG.inneighbors(g::CycleGraph, v::Integer) = outneighbors(g, v)

# ---- neighbors vector -----------------------------------

function Base.size(nbs::OutNeighborVector{V, G}) where {V, G <: CycleGraph}

    g = nbs.graph
    T = eltype(g)
    nvg::T = nv(g)

    nvg <= one(T) && return (0,)
    nvg == T(2) && return (1,)
    return (2,)
end

@inline function Base.getindex(nbs::OutNeighborVector{V, G}, i::Int) where {V, G <: CycleGraph}

    @boundscheck i ∈ Base.OneTo(length(nbs)) || throw(BoundsError(nbs, i))

    g = nbs.graph
    T = eltype(g)
    v::T = nbs.vertex
    nvg::T = nv(g)

    if i == 1
        v == one(T) && return T(2)
        v == nvg && return one(T)
        return v - one(T)
    end
    # i == 2
    v == one(T) && return nvg
    v == nvg && return nvg - T(1)
    return v + one(T)
end

Base.IndexStyle(::Type{<:OutNeighborVector{V, G}}) where {V, G <: CycleGraph} = IndexLinear()

# =======================================================
#         converting
# =======================================================

Base.convert(::Type{SimpleGraph}, g::CycleGraph{T}) where {T} = cycle_graph(nv(g))
Base.convert(::Type{SimpleGraph{T}}, g::CycleGraph) where {T} = cycle_graph(T(nv(g)))

# =======================================================
#         overrides
# =======================================================

# we use this check so that we have the same convention as in LightGraphs
LG.is_connected(g::CycleGraph) = nv(g) > 0

function LG.connected_components(g::CycleGraph)

    nv(g) == 0 && return typeof(vertices(g))[]
    return [vertices(g)]
end

# has_self_loops is defined in terms of this
LG.num_self_loops(g::CycleGraph) = 0

LG.is_bipartite(g::CycleGraph) = (nv(g) == 1) | iseven(nv(g))

function LG.squash(g::CycleGraph)
    nvg = nv(g)
    for T ∈ (UInt8, UInt16, UInt32, UInt64)
        nv(g) < typemax(T) && return CycleGraph{T}(nvg)
    end
end

# ---- degree -----------------------------------

function LG.Δ(g::CycleGraph)
    nvg = nv(g)
    return typemin(Int) * (nvg == 0) + (nvg >= 2) + (nvg >= 3)
end

LG.Δout(g::CycleGraph) = Δ(g)
LG.Δin(g::CycleGraph) = Δ(g)

function LG.δ(g::CycleGraph)
    nvg = nv(g)
    return typemax(Int) * (nvg == 0) + (nvg >= 2) + (nvg >= 3)
end

LG.δout(g::CycleGraph) = δ(g)
LG.δin(g::CycleGraph) = δ(g)

LG.degree(g::CycleGraph, vs::AbstractVector=vertices(g)) = Fill(Δ(g), length(vs))

LG.indegree(g::CycleGraph, vs::AbstractVector=vertices(g)) = degree(g, vs)
LG.outdegree(g::CycleGraph, vs::AbstractVector=vertices(g)) = degree(g, vs)

ext: jl


text:
s1 = "Yo soy una cadena."

s2 = """Yo también soy una cadena. """

"Aquí me sale un "error" porque es ambiguo dónde se acaba la cadena. "

"""Mira Mamá, sin "errors"!!! """

typeof('a')

'Esto va a dar un error'

nombre = "Jane"
num_dedos = 10
num_dedos_de_los_pies = 10

println("Hola, me llamo $nombre.")
println("Yo tengo $num_dedos dedos y $num_dedos_de_los_pies. ¡Esos son $(num_dedos + num_dedos_de_los_pies) dígitos en total!!")

string("¿Cuántos gatos ", "son demasiados?")

string("No lo sé, pero ", 10, " son muy pocos.")

s3 = "¿Cuántos gatos ";
s4 = "son demasiados?";

s3*s4

"$s3$s4"

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
using Knet: @diff, Param, value, grad, params
using Knet: sigm, tanh, softmax

mutable struct Layer
    wai::Param
    was::Param
    ba::Param
    wri::Param
    wrs::Param
    br::Param
    wsi::Param
    bs::Param
    state
end

Layer(in_size, out_size) =
begin
    sq = sqrt(2/(in_size+out_size))
    wai = Param(2*sq .* randn(Float32, in_size, out_size) .- sq)
    was = Param(2*sq .* randn(Float32, out_size, out_size) .- sq)
    ba  = Param(zeros(Float32, 1, out_size))
    wri = Param(2*sq .* randn(Float32, in_size, out_size) .- sq)
    wrs = Param(2*sq .* randn(Float32, out_size, out_size) .- sq)
    br  = Param(zeros(Float32, 1, out_size))
    wsi = Param(2*sq .* randn(Float32, in_size, out_size) .- sq)
    bs  = Param(zeros(Float32, 1, out_size))
    state = zeros(Float32, 1, out_size)
    layer = Layer(wai, was, ba, wri, wrs, br, wsi, bs, state)
layer
end

(layer::Layer)(in) =
begin
    focus  = sigm.(in * layer.wai + layer.state * layer.was + layer.ba)
    keep   = sigm.(in * layer.wri + layer.state * layer.wrs + layer.br)
    interm = tanh.(in * layer.wsi + layer.state .* focus + layer.bs)
    layer.state = keep .* interm + (1 .- keep) .* layer.state
end

mk_def(in_size, layers, out_size) =
begin
    hm_layers = length(layers)+1

    model_def = "struct Model\n"
    for i in 1:hm_layers
        model_def *= "l$i::Layer\n"
    end
    model_def *= "end"

    eval(Meta.parse(model_def))

    model_call = "(model::Model)(io) =\n"
    model_call *= "begin\n"
    for i in 1:hm_layers
        model_call *= "io = model.l$i(io)\n"
    end
    model_call *= "end"

    eval(Meta.parse(model_call))

end

mk_layers(in_size, layers, out_size) =
begin
    hm_layers = length(layers)+1

    hm_layers = length(layers) +1
    internal = Layer[]
    for i in 1:hm_layers
        if     i == 1         layer = Layer(in_size, layers[i])
        elseif i == hm_layers layer = Layer(layers[end], out_size)
        else                  layer = Layer(layers[i-1], layers[i])
        end
        push!(internal, layer)
    end
internal
end


text:
using RobertoMD
using Test

@testset "monoliquid" begin

    control=Dict("dt"=>0.01,
                "zero velocity"=>true,
                "steps"=>1000,
                "velocity verlet"=>true,
                "BerendsenNVT"=>Dict("tau"=>1.0),
                "Canonical field"=>Dict(
                "χ"=>Dict("1"=>[0.0]),
                "κ"=>0.1,
                "uniform mesh"=>true,
                "update"=>1,
                "Lcell"=>1.0),
                "thermo information"=>Dict( "freq"=>100,
                                            "energy"=>true,
                                        "momentum"=>true))

    Waters=Dict()
    for i in 1:10
        Waters[string(i)]=Dict("atoms"=>Dict("1"=>Dict("type"=>1,"mass"=>1,"coords"=>rand(Float64, 3)*4)))
    end

    configs=Dict(
        "box"=>[5.0,5.0,5.0],
        "molecules"=>Waters
        )

    Simulate(control,configs)
end


text:
struct Sum{H} end
struct LocalOperator{T, N, O, MT} end
struct SubHam{TT, OP} end
struct Vertex end

const Alias{T, MT} = Sum{SubHam{NTuple{1, Vertex}, LocalOperator{T, 1, 2, MT}}}

Alias{Float64}

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    l = Exo2Binion(Exo2Real(1,2,3,4))
    r = Exo2Binion(1,0,2,0,3,0,4)
    l == r
end

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(Exo2Binion(a))
end

@test begin
    a = 1
    !isreal(Exo2Binion(a, a, a, a))
end

@test begin
    z = random(Exo2Binion{Int})
    z == +(z)
end

@test begin
    a = CayleyDickson.randomBigFloat()
    real(Exo2Binion(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, Exo2Binion(1.0, 2, 3, 4, 5, 6, 7, 8))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, a: 3.0, ia: 4.0, b: 5.0, ib: 6.0, ab: 7.0, (ia)b: 8.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo2Binion(1.0, 2, 3, 4, 5, 6, 7))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, a: 3.0, ia: 4.0, b: 5.0, ib: 6.0, ab: 7.0, (ia)b: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo2Binion(1.0, 2, 3, 4, 5, 6))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, a: 3.0, ia: 4.0, b: 5.0, ib: 6.0, ab: 0.0, (ia)b: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo2Binion(1.0, 2, 3, 4, 5))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, a: 3.0, ia: 4.0, b: 5.0, ib: 0.0, ab: 0.0, (ia)b: 0.0]"
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo2Binion(1.0, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1.0, i: 2.0, a: 3.0, ia: 4.0, b: 0.0, ib: 0.0, ab: 0.0, (ia)b: 0.0]"
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    d = CayleyDickson.randomBigFloat()
    l = Exo2Binion(Exo1Binion(a, b, c, d))
    r = Exo2Binion(a, b, c, d)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    l = Exo2Binion(Binion(a, b))
    r = Exo2Binion(a, b)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    b = CayleyDickson.randomBigFloat()
    c = CayleyDickson.randomBigInt()
    l = Exo2Binion(a, b, c)
    r = Exo2Binion(a, b, c, 0)
    l == r
end

@test begin
    x = random(Exo2Binion{BigFloat})
    y = random(Exo2Binion{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo2Binion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = random(Exo1Binion{BigInt})
    z = random(Exo2Binion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo2Binion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    a = random(Exo1Binion{BigInt})
    z = random(Exo2Binion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = x * y
    r = y * x
    l != r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    z = random(Exo2Binion{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    z = random(Exo2Binion{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l != r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = alternatorL(x, y)
    r = zero(Exo2Binion{BigInt})
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = alternatorR(x, y)
    r = zero(Exo2Binion{BigInt})
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(Exo2Binion{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(Exo2Binion{Rational{BigInt}})
    y = random(Exo2Binion{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(Exo2Binion{Rational{BigInt}})
    y = random(Exo2Binion{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(Exo2Binion{Rational{BigInt}})
    y = random(Exo2Binion{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    z = random(Exo2Binion{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo2Binion{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    a = random(Exo1Binion{BigInt})
    x = random(Exo2Binion{BigInt})
    l = commutator(a, x)
    r = zero(Exo2Binion{BigInt})
    l != r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    z = random(Exo2Binion{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    z = random(Exo2Binion{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    z = random(Exo2Binion{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Exo2Binion{BigInt})
    quadrance(x) > big(0) || quadrance(x) == big(0)
end

@test begin
    x = random(Exo2Binion{BigInt})
    y = random(Exo2Binion{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(Exo2Binion(0, 0.0))
end

@test_throws ErrorException begin
    random(Exo2Binion{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Exo2Binion{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo2Binion{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo2Binion{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Exo2Binion{BigInt})
    z == z + zero(Exo2Binion{BigInt})
end

@test begin
    z = random(Exo2Binion{BigInt})
    z == z * one(Exo2Binion{BigInt})
end

@test begin
    z = random(Exo2Binion{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Exo2Binion{BigInt})
    z == z * one(z)
end


text:
function calculate_initial_conditions!(sim::Simulation, inputs::SimulationInputs)
    @debug "Start state intialization routine"
    sys = get_system(inputs)
    initial_guess = sim.x0_init
    res = PSY.solve_powerflow!(sys)
    if !res
        error("Power Flow fail")
    end
    var_count = get_variable_count(inputs)

    @debug "Setting up initilization indexing"
    bus_size = get_bus_count(inputs)
    bus_vars_count = 2 * bus_size
    bus_range = 1:bus_vars_count
    injection_start = get_injection_pointer(inputs)
    injection_count = 1
    branches_start = get_branches_pointer(inputs)
    branches_count = 1

    @debug "Updating Voltage guess"
    for bus in PSY.get_components(PSY.Bus, sys)
        @debug PSY.get_name(bus)
        @debug "V_r" PSY.get_magnitude(bus) * cos(PSY.get_angle(bus))
        @debug "V_i" PSY.get_magnitude(bus) * sin(PSY.get_angle(bus))
        bus_n = PSY.get_number(bus)
        bus_ix = get_lookup(inputs)[bus_n]
        initial_guess[bus_ix] = PSY.get_magnitude(bus) * cos(PSY.get_angle(bus))
        initial_guess[bus_ix + bus_size] = PSY.get_magnitude(bus) * sin(PSY.get_angle(bus))
    end

    @debug "Updating Source internal voltages"
    sources = PSY.get_components(PSY.Source, sys)
    if !isempty(sources)
        for s in sources
            initialize_device!(s)
        end
    end

    @debug "Updating Dynamic Injection Component Initial Guess"
    for d in get_injectors_data(inputs)
        dynamic_device = PSY.get_dynamic_injector(d)
        @debug PSY.get_name(d) typeof(d) typeof(dynamic_device)
        bus = PSY.get_bus(d)
        bus_n = PSY.get_number(PSY.get_bus(d))
        bus_ix = get_lookup(inputs)[bus_n]
        n_states = PSY.get_n_states(dynamic_device)
        ix_range = range(injection_start, length = n_states)
        injection_start = injection_start + n_states
        x0_device = initialize_device!(d)
        @assert length(x0_device) == n_states
        initial_guess[ix_range] = x0_device
    end

    @debug "Updating Component Initial Guess"
    dyn_branches = PSY.get_components(PSY.DynamicBranch, sys)
    if !isempty(dyn_branches)
        for br in dyn_branches
            @debug PSY.get_name(br) typeof(br)
            arc = PSY.get_arc(br)
            n_states = PSY.get_n_states(br)
            from_bus_number = PSY.get_number(arc.from)
            to_bus_number = PSY.get_number(arc.to)
            bus_ix_from = get_lookup(inputs)[from_bus_number]
            bus_ix_to = get_lookup(inputs)[to_bus_number]
            ix_range = range(branches_start, length = n_states)
            branches_start = branches_start + n_states
            x0_branch = initialize_device!(br)
            @assert length(x0_branch) == n_states
            initial_guess[ix_range] = x0_branch
        end
    else
        @debug "No Dynamic Branches in the system"
    end

    dx0 = zeros(var_count) #Define a vector of zeros for the derivative
    inif! = (out, x) -> system!(
        out,    #output of the function
        dx0,    #derivatives equal to zero
        x,      #states
        inputs,    #Parameters
        0.0,    #time equals to zero.
    )
    #Refine initial solution
    @debug "Start NLSolve Run"
    sys_solve = NLsolve.nlsolve(
        inif!,
        initial_guess,
        xtol = STRICT_NL_SOLVE_TOLERANCE,
        ftol = STRICT_NL_SOLVE_TOLERANCE,
        method = :trust_region,
    ) #Solve using initial guess x0
    if !NLsolve.converged(sys_solve)
        @warn("Initialization failed, initial conditions do not meet conditions for an stable equilibrium.\nTrying to solve again reducing numeric tolerance from $(STRICT_NL_SOLVE_TOLERANCE):")
        sys_solve = NLsolve.nlsolve(
            inif!,
            initial_guess,
            xtol = RELAXED_NL_SOLVE_TOLERANCE,
            ftol = RELAXED_NL_SOLVE_TOLERANCE,
            method = :trust_region,
        ) #Solve using initial guess x0
        if NLsolve.converged(sys_solve)
            @info("Initialization succeeded with a relaxed tolerance of $(RELAXED_NL_SOLVE_TOLERANCE). Saving solution")
        else
            @warn("Initialization failed again. Initial conditions do not meet conditions for an stable equilibrium\nSaving best result.")
        end
    end
    @debug "Write result to initial guess vector"
    initial_guess .= sys_solve.zero
    return NLsolve.converged(sys_solve)
end

"""
Returns a Dictionary with the resulting initial conditions of the simulation
"""
function get_initial_conditions(sim::Simulation)
    system = get_system(sim.simulation_inputs)
    bus_size = get_bus_count(sim.simulation_inputs)
    V_R = Dict{Int64, Float64}()
    V_I = Dict{Int64, Float64}()
    Vm = Dict{Int64, Float64}()
    θ = Dict{Int64, Float64}()
    for bus in PSY.get_components(PSY.Bus, system)
        bus_n = PSY.get_number(bus)
        bus_ix = get_lookup(sim.simulation_inputs)[bus_n]
        V_R[bus_n] = sim.x0_init[bus_ix]
        V_I[bus_n] = sim.x0_init[bus_ix + bus_size]
        Vm[bus_n] = sqrt(sim.x0_init[bus_ix]^2 + sim.x0_init[bus_ix + bus_size]^2)
        θ[bus_n] = angle(sim.x0_init[bus_ix] + sim.x0_init[bus_ix + bus_size] * 1im)
    end
    results = Dict{String, Any}("V_R" => V_R, "V_I" => V_I, "Vm" => Vm, "θ" => θ)
    for device in PSY.get_components(PSY.DynamicInjection, system)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        global_index = get_global_index(sim.simulation_inputs)[name]
        x0_device = Dict{Symbol, Float64}()
        for (i, s) in enumerate(states)
            x0_device[s] = sim.x0_init[global_index[s]]
        end
        results[name] = x0_device
    end
    dyn_branches = PSY.get_components(PSY.DynamicBranch, system)
    if !isempty(dyn_branches)
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            global_index = get_global_index(sim.simulation_inputs)[name]
            x0_br = Dict{Symbol, Float64}()
            for (i, s) in enumerate(states)
                x0_br[s] = sim.x0_init[global_index[s]]
            end
            printed_name = "Line " * name
            results[printed_name] = x0_br
        end
    end
    return results
end


text:
##############################################################################
#
# Basic Corpus type
#
##############################################################################

# TODO: Make this a parametric type?

type Corpus
    documents::Vector{GenericDocument}
    total_terms::Int
    lexicon::Dict{Compat.UTF8String, Int}
    inverse_index::Dict{Compat.UTF8String, Vector{Int}}
    h::TextHashFunction
end

function Corpus(docs::Vector{GenericDocument})
    Corpus(
        docs,
        0,
        Dict{Compat.UTF8String, Int}(),
        Dict{Compat.UTF8String, Vector{Int}}(),
        TextHashFunction()
    )
end

Corpus(docs::Vector{Any}) = Corpus(convert(Array{GenericDocument,1}, docs))

##############################################################################
#
# Construct a Corpus from a directory of text files
#
##############################################################################

function DirectoryCorpus(dirname::AbstractString)
    # Recursive descent of directory
    # Add all non-hidden files to Corpus

    docs = GenericDocument[]

    function add_files(dirname::AbstractString)
        if !isdir(dirname)
            error("DirectoryCorpus() can only be called on directories")
        end

        starting_dir = pwd()

        cd(dirname)
        for filename in readdir(".")
            if isfile(filename) && !ismatch(r"^\.", filename)
                push!(docs, FileDocument(abspath(filename)))
            end
            if isdir(filename) && !islink(filename)
                add_files(filename)
            end
        end
        cd(starting_dir)
    end

    add_files(dirname)

    return Corpus(docs)
end

##############################################################################
#
# Basic Corpus properties
#
##############################################################################

documents(c::Corpus) = c.documents
Base.length(crps::Corpus) = length(crps.documents)

##############################################################################
#
# Convert a Corpus to a DataFrame
#
##############################################################################

function Base.convert(::Type{DataFrame}, crps::Corpus)
    df = DataFrame()
    n = length(crps)
    df["Language"] = DataArray(Compat.UTF8String, n)
    df["Name"] = DataArray(Compat.UTF8String, n)
    df["Author"] = DataArray(Compat.UTF8String, n)
    df["TimeStamp"] = DataArray(Compat.UTF8String, n)
    df["Length"] = DataArray(Int, n)
    df["Text"] = DataArray(Compat.UTF8String, n)
    for i in 1:n
        d = crps.documents[i]
        df[i, "Language"] = string(language(d))
        df[i, "Name"] = name(d)
        df[i, "Author"] = author(d)
        df[i, "TimeStamp"] = timestamp(d)
        df[i, "Length"] = length(d)
        df[i, "Text"] = text(d)
    end
    return df
end

##############################################################################
#
# Treat a Corpus as an iterable
#
##############################################################################

Base.start(crps::Corpus) = 1
Base.next(crps::Corpus, ind::Int) = (crps.documents[ind], ind + 1)
Base.done(crps::Corpus, ind::Int) = ind > length(crps.documents)

##############################################################################
#
# Treat a Corpus as a container
#
##############################################################################

Base.push!(crps::Corpus, d::AbstractDocument) = push!(crps.documents, d)
Base.pop!(crps::Corpus) = pop!(crps.documents)

Base.unshift!(crps::Corpus, d::AbstractDocument) = unshift!(crps.documents, d)
Base.shift!(crps::Corpus) = shift!(crps.documents)

function Base.insert!(crps::Corpus, index::Int, d::AbstractDocument)
    insert!(crps.documents, index, d)
end
Base.delete!(crps::Corpus, index::Integer) = delete!(crps.documents, index)

##############################################################################
#
# Indexing into a Corpus
#
# (a) Numeric indexing just provides the n-th document
# (b) String indexing is effectively a trivial search engine
#
##############################################################################

Base.getindex(crps::Corpus, ind::Real) = crps.documents[ind]
Base.getindex{T <: Real}(crps::Corpus, inds::Vector{T}) = crps.documents[inds]
Base.getindex(crps::Corpus, r::Range) = crps.documents[r]
Base.getindex(crps::Corpus, term::AbstractString) = get(crps.inverse_index, term, Int[])

##############################################################################
#
# Assignment into a Corpus
#
##############################################################################

function Base.setindex!(crps::Corpus, d::AbstractDocument, ind::Real)
    crps.documents[ind] = d
    return d
end

##############################################################################
#
# Basic Corpus properties
#
# TODO: Offer progressive update that only changes based on current document
#
##############################################################################

lexicon(crps::Corpus) = crps.lexicon

function update_lexicon!(crps::Corpus, doc::AbstractDocument)
    ngs = ngrams(doc)
    for (ngram, counts) in ngs
        crps.total_terms += counts
        crps.lexicon[ngram] = get(crps.lexicon, ngram, 0) + counts
    end
end

function update_lexicon!(crps::Corpus)
    crps.total_terms = 0
    crps.lexicon = Dict{Compat.UTF8String,Int}()
    for doc in crps
        update_lexicon!(crps, doc)
    end
end

lexicon_size(crps::Corpus) = length(keys(crps.lexicon))
lexical_frequency(crps::Corpus, term::AbstractString) =
    (get(crps.lexicon, term, 0) / crps.total_terms)

##############################################################################
#
# Work with the Corpus's inverse index
#
# TODO: offer progressive update that only changes based on current document
#
##############################################################################

inverse_index(crps::Corpus) = crps.inverse_index

function update_inverse_index!(crps::Corpus)
    idx = Dict{Compat.UTF8String, Array{Int, 1}}()
    for i in 1:length(crps)
        doc = crps.documents[i]
        ngram_arr = isa(doc, NGramDocument) ? collect(keys(ngrams(doc))) : tokens(doc)
        ngram_arr = convert(Array{Compat.UTF8String,1}, ngram_arr)
        for ngram in ngram_arr
            if haskey(idx, ngram)
                push!(idx[ngram], i)
            else
                idx[ngram] = [i]
            end
        end
    end
    for key in keys(idx)
        idx[key] = unique(idx[key])
    end
    crps.inverse_index = idx
    nothing
end

index_size(crps::Corpus) = length(crps.inverse_index)

##############################################################################
#
# Every Corpus prespecifies a hash function for hash trick analysis
#
##############################################################################

hash_function(crps::Corpus) = crps.h
hash_function!(crps::Corpus, f::TextHashFunction) = (crps.h = f; nothing)

##############################################################################
#
# Standardize the documents in a Corpus to a common type
#
##############################################################################

function standardize!{T <: AbstractDocument}(crps::Corpus, ::Type{T})
    for i in 1:length(crps)
        crps.documents[i] = convert(T, crps.documents[i])
    end
end


text:
import SDDGeometry: AbstractLine, AbstractRay, AbstractLineSegment,
    AbstractCircle, AbstractArc, basepointR2, initialpointR2, finalpointR2,
    angle, centerR2, radius

_drawingkind = :png

"""
Set the current kind of drawing
- `:pixels` or `png`: Drawing pixels. Support depends on backend.
- `:vectorial` or `svg`: Drawing using vectorial shapes. Support depends on backend.
"""
function drawingkind(dk::Symbol)
    if dk == :svg || dk == :vectorial
        global _drawingkind = :svg
    else
        global _drawingkind = :png
    end
end # function

"Returns the current drawing kind symbol."
drawingkind() = _drawingkind

"Create a new drawing, using drawing kind and canvas data."
newdrawing() = _backend.newdrawing()

_axes = false

"Set true to draw coordinated axes."
axes(b::Bool) = global _axes = b

"Drawing coordinated axes?"
axes() = _axes

"Set the point size."
pointsize(ps::Real) = _backend.pointsize(ps)

"Get the point size."
pointsize() = _backend.pointsize()

"Set the line width."
linewidth(lw::Real) = _backend.linewidth(lw)

"Get the line width."
linewidth() = _backend.linewidth()

"""
Set the drawing style for shapes. Options:
- `:stroke`.
- `:fill`.
- `:fillstroke`.
- `:fillgrad`: Using a gradient coloring fill.
"""
style(s::Symbol) = _backend.style(s)

"Get the drawing style."
style() = _backend.style()

"Draw the (i,j) pixel."
drawpixel(i::Integer, j::Integer) = _backend.drawpixel(i,j)

"Draw the point (x,y)."
drawpoint(x::Real, y::Real) = _backend.drawpoint(x,y)

"Draw the point z."
drawpoint(z::Complex) = _backend.drawpoint(z)

"Draw the point x."
drawpoint(x::Real) = _backend.drawpoint(x)

drawlinesegment(x1::Real, y1::Real, x2::Real, y2::Real) =
_backend.drawlinesegment(x1,y1,x2,y2)
drawlinesegment(z1::Number, z2::Number) = _backend.drawlinesegment(z1,z2)
drawlinesegment(l::AbstractLineSegment) = _backend.drawlinesegment(l)

drawline(x::Real, y::Real, θ::Real) = _backend.drawline(x,y,θ)
drawline(z::Number, θ::Real) = _backend.drawline(z,θ)
drawline(l::AbstractLine) = _backend.drawline(l)

drawcircle(x::Real, y::Real, r::Real) = _backend.drawcircle(x,y,r)
drawcircle(z::Number, r::Real) = _backend.drawcircle(z,r)
drawcircle(c::AbstractCircle) = _backend.drawcircle(c)

drawarc(x1::Real,y1::Real,x2::Real,y2::Real,x3::Real,y3::Real) =
_backend.drawarc(x1,y1,x2,y2,x3,y3)
drawarc(z1::Number,z2::Number,z3::Number) = _backend.drawarc(z1,z2,z3)
drawarc(a::AbstractArc) = _backend.drawarc(a)
drawarc(x1::Real,y1::Real,x2::Real,y2::Real) =
_backend.drawarc(x1,y1,x2,y2)
drawarc(z1::Number,z2::Number) = _backend.drawarc(z1,z2)

"Returns the current created drawing."
drawing() = _backend.drawing()


text:
# dispatch on matrix inner product to handle sparse JuMP container
import LinearAlgebra.tr
function tr(X::SparseMatrixCSC, Y::JuMP.Containers.SparseAxisArray)
    return sum(v*Y[j,k] for (j,k,v) in zip(findnz(X)...))
end

# A + S ∈ PSDCone(), where S ⪰ 0
"""
    build_constraints_lmi!(m::JuMP.Model, A::SparseMatrixCSC{JuMP.AffExpr, Int}; verbose=false)

Adds the constraint `A + S ∈ PSDCone()`, where `S ⪰ 0` to JuMP model `m`
"""
function build_constraints_lmi!(m::JuMP.Model, A::SparseMatrixCSC{JuMP.AffExpr, Int}; verbose=false)
    perm, iperm, Tls = get_selectors(A; verbose=verbose, ret_cliques=false)
    num_cliques = length(Tls)
    A = A[perm, perm]
    S = Vector{AbstractMatrix}(undef, num_cliques)
    for p=1:num_cliques
            len_p = size(Tls[p], 1)
            S[p] = sparse(@variable(m, [1:len_p, 1:len_p] in PSDCone(), base_name="Cℓ$p"))
            for (j, k, v) ∈ zip(findnz(-Tls[p]'*S[p]*Tls[p])...)
                add_to_expression!(A[j,k], v)
            end
            # @info "Build constraint $p of $num_cliques"
    end
    for p = 1:num_cliques
            @constraint(m, Tls[p]*A*Tls[p]' .== 0)
    end
end

"""
    build_constraints_lmi!(m::JuMP.Model, A::SparseMatrixCSC{JuMP.AffExpr, Int}; verbose=false)

Adds the constraint `F_1 y_1 + F_2 y_2 + ... + F_n y_n + G ∈ PSDCone()` to JuMP model `m`
"""
function build_constraints_lmi!(
        m::JuMP.Model,
        y::Vector{VariableRef},
        F::AbstractVector{SparseMatrixCSC{T, S}},
        G::SparseMatrixCSC{T, S};
        verbose=false
) where {T <: Number, S <: Integer}
    build_constraints_primal!(m, build_A(y, F, G); verbose=verbose)
end

# A = F_1 y_1 + F_2 y_2 + ... + F_n y_n + G
"""
    build_A(y::Vector{JuMP.VariableRef}, F::AbstractVector{SparseMatrixCSC{T, S}}, G::SparseMatrixCSC{T, S}) where {T <: Number, S <: Integer}

Builds a JuMP.GenericAffExpr `A = F_1 y_1 + F_2 y_2 + ... + F_n y_n + G`.
"""
function build_A(y::Vector{JuMP.VariableRef}, F::AbstractVector{SparseMatrixCSC{T, S}}, G::SparseMatrixCSC{T, S}) where {T <: Number, S <: Integer}
    n = length(y)
    A = spzeros(GenericAffExpr{Float64, VariableRef}, n, n)

    # This is a hack that seems to work...
    for (j, k, v) ∈ zip(findnz(sparsity_pattern([F..., G]))...)
        A[j,k] = 1
    end

    for i in 1:n
        for (j, k, v) ∈ zip(findnz(F[i])...)
            add_to_expression!(A[j,k], F[i][j,k], y[i])
        end
    end
    for (j, k, v) ∈ zip(findnz(G)...)
        add_to_expression!(A[j,k], G[j,k])
    end

    # hack part II
    for (j, k, v) ∈ zip(findnz(sparsity_pattern([F..., G]))...)
        A[j,k] -= 1
    end
    return A
end

"""
    get_sparsity_pattern_from_cliques(cliques)

Returns a sparse matrix `sp` corresponding to the graph with maximal cliques `cliques`
"""
function get_sparsity_pattern_from_cliques(cliques)
    n = maximum(maximum.(cliques))

    sp = spzeros(Float64, n, n)
    for c in cliques
        for i in c, j in c
            sp[i,j] = 1.0
        end
    end

    return sp
end

"""
    reconstruct_from_sparse_varref(Zref, n)

Gets the value of JuMP.Containers.SparseAxisArray `Zref` and returns the result
as a SparseMatrixCSC
"""
function reconstruct_from_sparse_varref(Zref, n)
    Zcv = value.(Zref)
    Z_uncomp = spzeros(Float64, n, n)
    for ind in eachindex(Zcv)
        Z_uncomp[ind[1], ind[2]] = Zcv[ind]
    end
    return Z_uncomp
end

#FIXME
function build_constraints_standard!(
        model::JuMP.Model,
        X,
        F::AbstractVector{SparseMatrixCSC{T, S}},
        G::SparseMatrixCSC{T, S},
        c::Vector{T};
        verbose=false
) where {T <: Number, S <: Integer}
    perm, iperm, Tls, cliques = get_selectors(sparsity_pattern([F..., G]); verbose=verbose, ret_cliques=true)
    num_cliques = length(Tls)

    G .= G[perm, perm]
    for i in 1:length(F)
        F[i] .= F[i][perm,perm]
    end

    for i in 1:length(c)
        @constraint(model, tr(F[i], X) == c[i])
    end

    for (p, cl) in enumerate(cliques)
        len_cl = length(cl)
        Xp = @variable(model, [1:len_cl, 1:len_cl], base_name="X$p")
        for i in 1:len_cl, j in 1:len_cl
            Xp[i,j] = Z[cl[i], cl[j]]
        end
        @constraint(model, Xp in PSDCone())
    end
end


text:
# functions to make simpleODE work

import PDESolver.evalResidual

@doc """
### SimpleODEMod.evalResidual

This function evaluates the simple ODE equation.

** Inputs **

*  `mesh` : Abstract mesh object
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Simple ODE equation object
*  `opts` : Options dictionary
*  `t`    :

**Outputs**

*  None

"""->
function evalResidual(mesh::AbstractMesh{Tmsh},
sbp::AbstractSBP, eqn::SimpleODEData{Tsol, Tres, Tdim},
opts::Dict, t = 0.0) where {Tmsh, Tsol, Tres, Tdim}

  myrank = mesh.myrank
  params = eqn.params

  # TODO: difference between eqn.params.time and eqn.t?
  eqn.t = t

  eqn.res = fill!(eqn.res, 0.0)         # zero out the residual for the next function evaluation

  # start communication right away
  if opts["parallel_type"] == 1
    params.time.t_send += @elapsed if mesh.commsize > 1
      sendParallelData(mesh, sbp, eqn, opts)
    end
  end

  params.time.t_volume += @elapsed evalEquation(mesh, sbp, eqn, opts, t)
#   params.time.t_face += @elapsed if mesh.isDG
#     evalFaceTerm(mesh, sbp, eqn, opts)
#   end
#   params.time.t_source += @elapsed evalSRCTerm(mesh, sbp, eqn, opts)
#   params.time.t_bndry += @elapsed evalBndry(mesh, sbp, eqn)
#   params.time.t_sharedface += @elapsed if mesh.commsize > 1
#     evalSharedFaceIntegrals(mesh, sbp, eqn, opts)
#   end

  # !!!!!!!!!!!
  # Simple ODE does not need use_Minv. Forcing it to be false.
  opts["use_Minv"] = false

  if opts["use_Minv"]
    applyMassMatrixInv3D(mesh, sbp, eqn, opts, eqn.res)
  end

end   # end of function evalResidual

function evalEquation(mesh::AbstractMesh{Tmsh},
                   sbp::AbstractSBP,
                   eqn::SimpleODEData{Tsol, Tres, Tdim}, opts,
                   t) where {Tmsh, Tsol, Tres, Tdim}

  dxidx = mesh.dxidx
  q = eqn.q

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      for k=1:Tdim  # loop over parametric dimensions

        # Here, the expression for du/dt needs to be stored in eqn.res
        # For example:
        #   if u = x^2 + t^4   =>   then res = 4*t^3
        #       do this with:  eqn.res[1, j, i] = calc_4t3(mesh.coords[:, j, i], eqn.params, t)
        #   if u = x^2 + t^3   =>   then res = 3*t^2
        #       do this with:  eqn.res[1, j, i] = calc_3t2(mesh.coords[:, j, i], eqn.params, t)
        #   if u = x^2 + t^2   =>   then res = 2t
        #       do this with:  eqn.res[1, j, i] = calc_2t(mesh.coords[:, j, i], eqn.params, t)
        #   if u = x^2 + t     =>   then res = 1
        #   if u = t           =>   then res = 1
        #       do this with:  eqn.res[1, j, i] = 1.0

        if opts["simpleODE_eqn"] == 1
          eqn.res[1, j, i] = calc_4t3(mesh.coords[:, j, i], eqn.params, t)
        elseif opts["simpleODE_eqn"] == 2
          eqn.res[1, j, i] = calc_3t2(mesh.coords[:, j, i], eqn.params, t)
        elseif opts["simpleODE_eqn"] == 3
          eqn.res[1, j, i] = calc_2t(mesh.coords[:, j, i], eqn.params, t)
        elseif opts["simpleODE_eqn"] == 4 || opts["simpleODE_eqn"] == 5
          eqn.res[1, j, i] = 1.0
        elseif opts["simpleODE_eqn"] == 6
          # for du/dt = u + 2*t
          eqn.res[1, j, i] = q[1, j, i] + calc_2t(mesh.coords[:, j, i], eqn.params, t)
        elseif opts["simpleODE_eqn"] == 7
          # for du/dt = u
          eqn.res[1, j, i] = q[1, j, i]
        else
          throw(ErrorException("Incorrect or no option simpleODE_eqn specified."))
        end # end of simpleODE_eqn handling

      end # end of loop over parametric dimensions
    end # end of loop over nodes on the element
  end # end of loop over elements

  return nothing

end

@doc """
### SimpleODEMod.init

  This function initializes the mesh and equation objects with any module
  specific data, such as boundary condition and source term functors.

  Inputs:
    mesh
    sbp
    eqn
    opts

  Outputs: none

  Aliasing restrictions: none
"""->
function init(mesh::AbstractMesh{Tmsh}, sbp::AbstractSBP,
eqn::SimpleODEData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}

  println("Entering Simple ODE Module")
  return nothing

end

@doc """
### SimpleODEMod.majorIterationCallback

  This function gets called by the NonlinearSolvers during every major
  iteration.  This provides the opportunity to do output or monitoring.

  Inputs:
    itr: major iteration number
    mesh
    sbp
    eqn
    opts

    Outputs: none

    Aliasing restrictions: none

"""->
function majorIterationCallback(itr::Integer, mesh::AbstractMesh, sbp::AbstractSBP,
                                eqn::AbstractSimpleODEData, opts, f::IO)

  return nothing

end

@doc """
### NonlinearSolvers.ode_pre_func

  The pre-function for solving ordinary differential equations with a physics
  module.  The only operation it performs is disassembling eqn.q_vec into
  eqn.q

  Compare with NonlinearSolvers.pde_pre_func in NonlinearSolvers/rk4.jl

  Inputs:
    mesh
    sbp
    eqn
    opts
"""->
function ode_pre_func(mesh, sbp, eqn, opts)

  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
end

@doc """
### NonlinearSolvers.pde_ode_func

  The post-function for solving ordinary differential equations with a physics
  module.  This function DOES NOT multiply by A0inv.
  It assembles eqn.res into eqn.res_vec.
  Then it DOES NOT multiply by the inverse mass matrix, then it calculates
  the SBP approximation to the integral L2 norm.

  Compare with NonlinearSolvers.pde_post_func in NonlinearSolvers/rk4.jl

  Inputs:
    mesh
    sbp
    eqn
    opts

"""->
function ode_post_func(mesh, sbp, eqn, opts; calc_norm=true)
  # IN pde_post_func:  eqn.multiplyA0inv(mesh, sbp, eqn, opts, eqn.res)
  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
  # IN pde_post_func:  for j=1:length(eqn.res_vec) eqn.res_vec[j] = eqn.Minv[j]*eqn.res_vec[j] end
  if calc_norm

    # since simpleODE doesn't use SBP, need to use standard 2-norm instead of SBP norm
    # IN pde_post_func:    local_norm = calcNorm(eqn, eqn.res_vec)
    local_norm = norm(eqn.res_vec, 2)/mesh.numDof
    eqn.params.time.t_allreduce += @elapsed global_norm = MPI.Allreduce(local_norm*local_norm, MPI.SUM, mesh.comm)
    return sqrt(global_norm)
  end

  return nothing
end

@doc """
### SimpleODeMod.applyMassMatrixInv3D

  This function applies the 3D inverse mass matrix to an array.
    The array passed in should always be eqn.res

  Inputs:
    mesh: mesh object, needed for numEl and numDofPerNode fields
    sbp: sbp object, needed for numnodes field
    eqn: equation object, needed for Minv3D field
    opts
    arr: the 3D array to have the 3D mass matrix inverse applied to it

"""->
function applyMassMatrixInv3D(mesh, sbp, eqn, opts, arr)

  for i = 1:mesh.numEl
    for j = 1:sbp.numnodes
      for k = 1:mesh.numDofPerNode
        arr[k, j, i] = eqn.Minv3D[k, j, i] * arr[k, j, i]
      end
    end
  end

  return arr
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule DASKR_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("DASKR")
JLLWrappers.@generate_main_file("DASKR", UUID("655fdf9c-4522-5e7e-b537-e46a52300ea6"))
end  # module DASKR_jll


text:
"""
   Solid(ρ, cₚ, k)

Propriedades termodinâmicas de sólidos
"""
struct Solid <: AbstractSubstance
    "Massa específica em kg/m³"
    ρ::Float64
    "Calor específico em J/kg.K"
    cₚ::Float64
    "Condutividade térmica W/m.K"
    k::Float64
end

"Óxido de ferro. Usado nos termistores NTC (?)"
const Fe2O3 = Solid(3e3, 650.0, 0.58)

"""
    thermaldiff(x::AbstractSubstance)

Difusividade térmica em m²/s
"""
thermaldiff(x::AbstractSubstance) = x.k / (x.ρ*x.cₚ)

"""
    densidade(x::AbstractSubstance)

Densidade em kg/m³
"""
density(x::AbstractSubstance) = x.ρ

"""
    specheat(x::AbstractSubstance)

Calor específico em J/kg.K
"""
specheat(x::AbstractSubstance) = x.cₚ

"""
    thermalcond(x::AbstractSubstance)

Condutividade térmica em W/m.K
"""
thermalcond(x::AbstractSubstance) = x.k

"""
    thermalcond(x::AbstractSubstance)

Condutividade térmica em W/m.K
"""
heatcond(x::AbstractSubstance) = x.k

"Óxido de ferro. Usado nos termistores NTC (?)"
const Fe₂O₃ = Solid(5242.0, 650.6, 0.58)


text:
"""
        moora_ref(decisionMat, weights, fns)

Apply MOORA (Multi-Objective Optimization By Ratio Analysis) method for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n candidate (or strategy) and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of function that are either maximum or minimum.

# Description
moora() applies the MOORA method to rank n strategies subject to m criteria which are supposed to be
either maximized or minimized. Note that this is the reference version of the MOORA method. For the
ratio method, look at `moora_ratio`.

# Output
- `::MooraResult`: MooraResult object that holds multiple outputs including scores and best index.

# Examples
```julia-repl
julia> w =  [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077];

julia> Amat = [
             100 92 10 2 80 70 95 80 ;
             80  70 8  4 100 80 80 90 ;
             90 85 5 0 75 95 70 70 ;
             70 88 20 18 60 90 95 85
           ];

julia> dmat = makeDecisionMatrix(Amat)
4×8 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4     Crt5     Crt6     Crt7     Crt8
     │ Float64  Float64  Float64  Float64  Float64  Float64  Float64  Float64
─────┼────────────────────────────────────────────────────────────────────────
   1 │   100.0     92.0     10.0      2.0     80.0     70.0     95.0     80.0
   2 │    80.0     70.0      8.0      4.0    100.0     80.0     80.0     90.0
   3 │    90.0     85.0      5.0      0.0     75.0     95.0     70.0     70.0
   4 │    70.0     88.0     20.0     18.0     60.0     90.0     95.0     85.0

julia> fns = makeminmax([maximum, maximum, maximum, maximum, maximum, maximum, maximum, maximum]);

julia> result = moora(dmat, w, fns)

julia> result.scores
4-element Array{Float64,1}:
 0.3315938731541169
 0.2901446390098523
 0.3730431072983815
 0.019265256092245782

julia> result.bestIndex
4
```

# References
Celikbilek Yakup, Cok Kriterli Karar Verme Yontemleri, Aciklamali ve Karsilastirmali
Saglik Bilimleri Uygulamalari ile. Editor: Muhlis Ozdemir, Nobel Kitabevi, Ankara, 2018

İşletmeciler, Mühendisler ve Yöneticiler için Operasyonel, Yönetsel ve Stratejik Problemlerin
Çözümünde Çok Kriterli Karar verme Yöntemleri, Editörler: Bahadır Fatih Yıldırım ve Emrah Önder,
Dora, 2. Basım, 2015, ISBN: 978-605-9929-44-8
"""
function moora_ref(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::MooraResult

    w = unitize(weights)

    nalternatives, ncriteria = size(decisionMat)

    normalizedMat = normalize(decisionMat)
    weightednormalizedMat = w * normalizedMat

    # cmaxs = colmaxs(weightednormalizedMat)
    cmaxs = apply_columns(fns, weightednormalizedMat)
    cmins = apply_columns(reverseminmax(fns), weightednormalizedMat)

    refmat = similar(weightednormalizedMat)

    #=
    # This implementation is buggy and will be removed
    # in later releases.
    # note: the problem is fns[rowind] is non-sense
    # because fns are due to columns not rows.
    for rowind in 1:nalternatives
        if fns[rowind] == maximum
            refmat[rowind, :] .= cmaxs - weightednormalizedMat[rowind, :]
        elseif fns[rowind] == minimum
            refmat[rowind, :] .= weightednormalizedMat[rowind, :] - cmins
        else
            @warn fns[rowind]
            error("Function must be either maximize or minimize")
        end
    end
    =#

    for rowind in 1:nalternatives
        for colind in 1:ncriteria
            if fns[colind] == maximum
                refmat[rowind, colind] = cmaxs[colind] - weightednormalizedMat[rowind, colind]
            elseif fns[colind] == minimum
                refmat[rowind, colind] = weightednormalizedMat[rowind, colind] - cmins[colind]
            else
                @warn fns[colind]
                error("Function must be either maximize or minimize")
            end
        end
    end

    scores = rmaxs = rowmaxs(refmat)

    bestIndex = sortperm(rmaxs) |> first

    result = MooraResult(
       :reference,
       decisionMat,
       w,
       weightednormalizedMat,
       refmat,
       scores,
       bestIndex
   )

    return result
end

"""
        moora_ratio(decisionMat, weights, fns)

Apply MOORA (Multi-Objective Optimization By Ratio Analysis) method for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n candidate (or strategy) and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of function that are either maximum or minimum.

# Description
moora() applies the MOORA method to rank n strategies subject to m criteria which are supposed to be
either maximized or minimized. Note that this is the ratio version of the MOORA method. For the
reference method, look at `moora_ref`.

# Output
- `::MooraResult`: MooraResult object that holds multiple outputs including scores and best index.

# References
KUNDAKCI, Nilsen. "Combined multi-criteria decision making approach based on MACBETH
and MULTI-MOORA methods." Alphanumeric Journal 4.1 (2016): 17-26.
"""
function moora_ratio(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::MooraResult
    w = unitize(weights)

    nalternatives, ncriteria = size(decisionMat)

    mat = Matrix{Float64}(decisionMat)
    normalizedMatrix = similar(mat)
    weightednormalizedMat = similar(mat)

    for i in 1:ncriteria
        normalizedMatrix[:, i] = mat[:, i] ./ sqrt(sum(mat[:, i] .^ 2.0))
        weightednormalizedMat[:, i] = normalizedMatrix[:, i] .* w[i]
    end

    scores = zeros(nalternatives)
    for i in 1:nalternatives
        for j in 1:ncriteria
            if fns[j] == maximum
                scores[i] += weightednormalizedMat[i, j]
            elseif fns[j] == minimum
                scores[i] -= weightednormalizedMat[i, j]
            else
                error("In Moora, direction of optimization must be either minimum or maximum.")
            end
        end
    end

    bestIndex = scores |> sortperm |> last

    return MooraResult(
        :ratio,
       decisionMat,
       w,
       DataFrame(weightednormalizedMat, :auto),
       nothing, # refmat
       scores,
       bestIndex
    )
end

"""
        moora_ratio(decisionMat, weights, fns; method = :reference)

Apply MOORA (Multi-Objective Optimization By Ratio Analysis) method for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n candidate (or strategy) and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of function that are either maximum or minimum.
 - `method::Symbol`: Either `:reference` or `:ratio`. By default, it is `:reference`.

# Description
moora() applies the MOORA method to rank n strategies subject to m criteria which are supposed to be
either maximized or minimized. This method has two different versions. The method parameter determines the method used. It is `:reference` by default. For the other version, it can be set to `:ratio`.

# Output
- `::MooraResult`: MooraResult object that holds multiple outputs including scores and best index.

# References
KUNDAKCI, Nilsen. "Combined multi-criteria decision making approach based on MACBETH
and MULTI-MOORA methods." Alphanumeric Journal 4.1 (2016): 17-26.

Celikbilek Yakup, Cok Kriterli Karar Verme Yontemleri, Aciklamali ve Karsilastirmali
Saglik Bilimleri Uygulamalari ile. Editor: Muhlis Ozdemir, Nobel Kitabevi, Ankara, 2018

İşletmeciler, Mühendisler ve Yöneticiler için Operasyonel, Yönetsel ve Stratejik Problemlerin
Çözümünde Çok Kriterli Karar verme Yöntemleri, Editörler: Bahadır Fatih Yıldırım ve Emrah Önder,
Dora, 2. Basım, 2015, ISBN: 978-605-9929-44-8
"""
function moora(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1}; method::Symbol = :reference)::MooraResult
    if method == :reference
        return moora_ref(decisionMat, weights, fns)
    elseif method == :ratio
        return moora_ratio(decisionMat, weights, fns)
    else
        @error "Method not found: " method
        @error "Moora is defined for methods :reference and :ratio"
        error("Terminating.")
    end
end

"""
        moora(setting; method = :reference)

Apply MOORA (Multi-Objective Optimization By Ratio Analysis) method for a given matrix and weights.

# Arguments:
 - `setting::MCDMSetting`: MCDMSetting object.
 - `method::Symbol`: Either `:reference` or `:ratio`. By default, it is `:reference`.

# Description
moora() applies the MOORA method to rank n strategies subject to m criteria which are supposed to be either maximized or minimized.

# Output
- `::MooraResult`: MooraResult object that holds multiple outputs including scores and best index.
"""
function moora(setting::MCDMSetting; method::Symbol = :reference)::MooraResult
    moora(
        setting.df,
        setting.weights,
        setting.fns,
        method = method
    )
end

max_stars_repo_name: jbytecode/jMCDM

max_issues_repo_name: jbytecode/jMCDM

max_forks_repo_name: jbytecode/jMCDM


text:
"""
$(TYPEDEF)

The solution to a `DataDrivenProblem` derived via a certain algorithm.
The solution is represented via an `AbstractBasis`, which makes it callable.

# Fields
$(FIELDS)
"""
struct DataDrivenSolution{R <: Union{AbstractBasis,Nothing} , S, P , A, O,M} <: AbstractDataDrivenSolution
    """Result"""
    res::R # The result
    """Returncode"""
    retcode::Symbol
    """Parameters used to derive the equations"""
    ps::S
    """Algorithm used for solution"""
    alg::A # Solution algorithm
    """Original Output"""
    out::O # E.g. coefficients
    """Original Input"""
    inp::P
    """Error metrics"""
    metrics::M # Named tuples
end

# Make it callable
(r::DataDrivenSolution)(args...) = r.res(args...)

is_implicit(r::DataDrivenSolution) = isa(r.alg, Optimize.AbstractSubspaceOptimizer)

"""
$(SIGNATURES)

Returns the result of in form of an `AbstractBasis`.
"""
result(r::DataDrivenSolution) = r.res

"""
$(SIGNATURES)

Returns the estimated parameters in form of an `Vector`.
"""
ModelingToolkit.parameters(r::DataDrivenSolution) = r.ps

"""
$(SIGNATURES)

Generate a mapping of the parameter values and symbolic representation useable
to `solve` and `ODESystem`.
"""
function parameter_map(r::DataDrivenSolution)
    return [
        ps_ => p_ for (ps_, p_) in zip(parameters(r.res), r.ps)
    ]
end

"""
$(SIGNATURES)

Returns the metrics of the result in form of a `NamedTuple`.
"""
metrics(r::DataDrivenSolution) = r.metrics

"""
$(SIGNATURES)

Returns the original output of the algorithm, e.g. an `AbstractArray` of coefficients for sparse regression.
"""
output(r::DataDrivenSolution) = r.out

"""
$(SIGNATURES)

Returns the algorithm used to derive the solution.
"""
algorithm(r::DataDrivenSolution) = r.alg

"""
$(SIGNATURES)

Returns the original inputs, most commonly the `DataDrivenProblem` and the `Basis` used to derive the solution.
"""
inputs(r::DataDrivenSolution) = r.inp

function Base.summary(io::IO, r::DataDrivenSolution)
    is_implicit(r) ? println(io,"Implicit Result") : println(io,"Explicit Result")
    println(io, "Solution with $(length(r.res)) equations and $(length(r.ps)) parameters.")
    println(io, "Returncode: $(r.retcode)")
    haskey(r.metrics, :Sparsity) && println(io, "Sparsity: $(r.metrics.Sparsity)")
    haskey(r.metrics, :Error) && println(io, "L2 Norm Error: $(r.metrics.Error)")
    haskey(r.metrics, :AICC) && println(io, "AICC: $(r.metrics.AICC)")
    return
end

function Base.print(io::IO, r::DataDrivenSolution, fullview::DataType)

    fullview != Val{true} && return summary(io, r)

    is_implicit(r) ? println(io,"Implicit Result") : println(io,"Explicit Result")
    println(io, "Solution with $(length(r.res)) equations and $(length(r.ps)) parameters.")
    println(io, "Returncode: $(r.retcode)")
    haskey(r.metrics, :Sparsity) && println(io, "Sparsity: $(r.metrics.Sparsity)")
    haskey(r.metrics, :Error) && println(io, "L2 Norm Error: $(r.metrics.Error)")
    haskey(r.metrics, :AICC) && println(io, "AICC: $(r.metrics.AICC)")
    println(io, "")
    print(io, r.res)
    println(io, "")
    if length(r.res.ps) > 0
        x = parameter_map(r)
        println(io, "Parameters:")
        for v in x
            println(io, "   $(v[1]) : $(v[2])")
        end
    end

    return
end

Base.print(io::IO, r::DataDrivenSolution) = summary(io, r)
Base.show(io::IO, r::DataDrivenSolution) = is_implicit(r) ? show(io,"Implicit Result") : show(io,"Explicit Result")

function build_parametrized_eqs(X::AbstractMatrix, b::Basis)
    # Create additional variables
    sp = Int(norm(X, 0))
    sps = norm.(eachcol(X), 0)
    inds = sps .> zero(eltype(X))
    pl = length(parameters(b))

    p = [Symbolics.variable(:p, i) for i in (pl+1):(pl+sp)]
    p = collect(p)
    ps = zeros(eltype(X), sp)

    eqs = zeros(Num, sum(inds))
    eqs_ = [e.rhs for e in equations(b)]
    cnt = 1
    for j in 1:size(X, 2)
        if sps[j] == zero(eltype(X))
            continue
        end
        for i in 1:size(X, 1)
            if iszero(X[i,j])
                continue
            end
            ps[cnt] = X[i,j]
            eqs[j] += p[cnt]*eqs_[i]
            cnt += 1
        end
    end
    return eqs, ps, p
end

# Explicit sindy
function build_solution(prob::DataDrivenProblem, Ξ::AbstractMatrix, opt::Optimize.AbstractOptimizer, b::Basis;
    eval_expression = false)
    if all(iszero.(Ξ))
        @warn "Sparse regression failed! All coefficients are zero."
        return DataDrivenSolution(
        nothing , :failed, nothing, opt, Ξ, (Problem = prob, Basis = b, nothing),
    )
    end

    eqs, ps, p_ = build_parametrized_eqs(Ξ, b)

    # Build the lhs
    if length(eqs) == length(states(b))
        xs = states(b)
        d = Differential(get_iv(b))
        eqs = [d(xs[i]) ~ eq for (i,eq) in enumerate(eqs)]
    end

    # Build a basis
    res_ = Basis(
        eqs, states(b),
        parameters = [parameters(b); p_], iv = get_iv(b),
        controls = controls(b), observed = observed(b),
        name = gensym(:Basis),
        eval_expression = eval_expression
    )

    sparsity = norm(Ξ, 0)
    sparsities = map(i->norm(i, 0), eachcol(Ξ))

    retcode = size(Ξ, 2) == size(prob.DX, 1) ? :success : :incomplete
    pnew = !isempty(parameters(b)) ? [prob.p; ps] : ps
    X = get_target(prob)
    x, _, t, c = get_oop_args(prob)
    Y = res_(x, pnew, t, c)

    # Build the metrics
    sparsity = norm(Ξ, 0)
    sparsities = map(i->norm(i, 0), eachcol(Ξ))
    inds = sparsities .> zero(eltype(X))
    error = norm(X[inds, :]-Y, 2)
    k = free_parameters(res_)
    aic = AICC(k, X[inds, :], Y)
    errors = zeros(eltype(X), sum(inds))
    aiccs = zeros(eltype(X), sum(inds))
    j = 1
    for i in 1:size(X,1)
        if inds[i]
            errors[i] = norm(X[i,:].-Y[j,:],2)
            aiccs[i] = AICC(k, X[i:i, :], Y[j:j,:])
            j += 1
        end
    end

    metrics = (
        Sparsity = sparsity,
        Error = error,
        AICC = aic,
        Sparsities = sparsities,
        Errors = errors,
        AICCs = aiccs,
    )

    inputs = (
        Problem = prob,
        Basis = b,
    )

    return DataDrivenSolution(
        res_, retcode, pnew, opt, Ξ, inputs, metrics
    )
end

function build_solution(prob::DataDrivenProblem, Ξ::AbstractMatrix, opt::Optimize.AbstractSubspaceOptimizer,
    b::Basis, implicits::Vector{Num}; eval_expression = false)

    if all(iszero(Ξ))
        @warn "Sparse regression failed! All coefficients are zero."
        return DataDrivenSolution(
        nothing , :failed, nothing, opt, Ξ, (Problem = prob, Basis = b, nothing),
    )
    end

    eqs, ps, p_ = build_parametrized_eqs(Ξ, b)
    eqs = [0 .~ eq for eq in eqs]

    # Build a basis
    res_ = Basis(
        collect(eqs), states(b),
        parameters = [parameters(b); p_], iv = get_iv(b),
        controls = controls(b), observed = observed(b),
        name = gensym(:Basis),
        eval_expression = eval_expression
    )

    sparsity = norm(Ξ, 0)
    sparsities = map(i->norm(i, 0), eachcol(Ξ))

    retcode = size(Ξ, 2) == size(prob.DX, 1) ? :success : :incomplete
    pnew = !isempty(parameters(b)) ? [prob.p; ps] : ps
    Y = res_([prob.X; prob.DX], pnew, prob.t, prob.U)

    # Build the metrics
    sparsity = norm(Ξ, 0)
    sparsities = map(i->norm(i, 0), eachcol(Ξ))
    inds = sparsities .> zero(eltype(Y))
    error = norm(Y, 2)
    k = free_parameters(res_)
    aic = AICC(k, Y, zero(Y))
    errors = zeros(eltype(Y), sum(inds))
    aiccs = zeros(eltype(Y), sum(inds))
    j = 1
    for i in 1:size(Y,1)
        if inds[i]
            errors[i] = norm(Y[j,:],2)
            aiccs[i] = AICC(k, Y[j:j,:], zero(Y[j:j,:]))
            j += 1
        end
    end

    metrics = (
        Sparsity = sparsity,
        Error = error,
        AICC = aic,
        Sparsities = sparsities,
        Errors = errors,
        AICCs = aiccs,
    )

    inputs = (
        Problem = prob,
        Basis = b,
        implicit_variables = implicits,
    )

    return DataDrivenSolution(
        res_, retcode, pnew, opt, Ξ, inputs, metrics
    )
end

function _round!(x::AbstractArray{T, N}, digits::Int) where {T, N}
    for i in eachindex(x)
        x[i] = round(x[i], digits = digits)
    end
    return x
end

function build_solution(prob::DataDrivenProblem, k, C, B, Q, P, inds, b::AbstractBasis, alg::AbstractKoopmanAlgorithm; digits::Int = 10, eval_expression = false)

    # Build parameterized equations, inds indicate the location of basis elements containing an input
    Ξ = zeros(eltype(B), size(C,2), length(b))

    Ξ[:, inds] .= real.(Matrix(k))
    if !isempty(B)
        Ξ[:, .! inds] .= B
    end

    # Transpose because of the nature of build_parametrized_eqs
    if !eval_expression
        eqs, ps, p_ = build_parametrized_eqs(_round!(C*Ξ, digits)', b)
    else
        K̃ = _round!(C*Ξ, digits)
        eqs = K̃*Num[states(b); controls(b)]
        p_ = []
        ps = [K̃...]
    end

    # Build the lhs
    if length(eqs) == length(states(b))
        xs = states(b)
        d = Differential(get_iv(b))
        eqs = [d(xs[i]) ~ eq for (i,eq) in enumerate(eqs)]
    end

    res_ = Koopman(eqs, states(b),
        parameters = [parameters(b); p_],
        controls = controls(b), iv = get_iv(b),
        K = k, C = C, Q = Q, P = P, lift = b.f,
        is_discrete = is_discrete(prob),
        eval_expression = eval_expression)

    retcode = :success
    X = get_target(prob)
    X_, p_, t, U = get_oop_args(prob)

    pnew = !isempty(parameters(b)) ? [p_; ps] : ps

    if !eval_expression
        # Equation space
        Y = res_(X_, pnew, t, U)
    else
        Y = K̃*b(X_, p_, t, U)
    end

    # Build the metrics
    error = norm(X-Y, 2)
    k = free_parameters(res_)
    aic = AICC(k, X, Y)
    errors = zeros(eltype(X), sum(inds))
    aiccs = zeros(eltype(X), sum(inds))
    j = 1
    for i in 1:size(X,1)
        errors[i] = norm(X[i,:].-Y[i,:],2)
        aiccs[i] = AICC(k, X[i:i, :], Y[i:i,:])
    end

    metrics = (
        Error = error,
        AICC = aic,
        Errors = errors,
        AICCs = aiccs,
    )

    inputs = (
        Problem = prob,
        Basis = b,
    )

    return DataDrivenSolution(
        res_, retcode, pnew, alg, Ξ, inputs, metrics
    )
end


text:
is_circ(e, set) = false
function is_circ(e::Expr, set)
    if e.head === :call
        if length(e.args) == 1 || ((length(e.args) == 2) && (e.args[2] ∈ set))
            # f() or f(x)
            return true
        elseif length(e.args) == 2
            # f(g(...))
            return is_circ(e.args[2], set)
        end
    end
    return false
end

function make_circ(e::Expr, set)
    if e.head === :call
        if length(e.args) == 1 || ((length(e.args) == 2) && (e.args[2] ∈ set))
            # f() or f(x)
            return e.args[1]
        elseif length(e.args) == 2
            # f(g(...))
            return Expr(:call, :(Base.:∘), e.args[1], make_circ(e.args[2], set))
        end
    end
end


text:
using Termbox

with_term() do
    for N in 1:4
        for row in 1:N
            for col in 1:N
                tb_change_cell(col, row, '#')
            end
        end
        tb_present()
        sleep(1e-2)
    end
end


text:
@with_kw struct ImageListRequest
    limit::Int = 50; @assert limit > 0
    offset::Int = 0; @assert offset >= 0
    sort::Symbol = :lowerName
    sortdir::Int = 1; @assert sortdir == 1 || sortdir == -1
    datasetId::String
end

ImageListRequest(ds::Union{ListEntry,DatasetMetadata}; kw...) = ImageListRequest(datasetId = ds.id; kw...)

function Base.get(req::ImageListRequest)
    query = "https://isic-archive.com:443/api/v1/image?limit=$(req.limit)&offset=$(req.offset)&sort=$(req.sort)&sortdir=$(req.sortdir)&datasetId=$(req.datasetId)"
    [ListEntry(o) for o in clean_json(get(query))]
end

# ====================================================================

mutable struct ImageMetadata
    id::String
    modelType::String
    created::DateTime
    #creatorId::String
    #description::String
    name::String
    updated::DateTime
    class::Symbol
    meta::Dict{String,Any}
end

function ImageMetadata(o::Dict)
    ImageMetadata(
        o["_id"],
        o["_modelType"],
        parse_datetime(o["created"]),
        #o["creatorId"],
        #o["description"],
        o["name"],
        parse_datetime(o["updated"]),
        extract_class(o["meta"]),
        o["meta"])
end

function JSON.json(o::ImageMetadata)
    dict = Dict{String, Any}()
    dict["_id"] = o.id
    dict["_modelType"] = o.modelType
    dict["created"] = tostring(o.created)
    #dict["creatorId"] = o.creatorId
    #dict["description"] = o.description
    dict["name"] = o.name
    dict["updated"] = tostring(o.updated)
    dict["meta"] = o.meta
    json(dict)
end

save(file, o::ImageMetadata) = writeall(file, json(o))
load_imagemeta(file) = ImageMetadata(JSON.parse(readall(file)))

for op = (:<, :>, :(==), :(!=), :(<=), :(>=))
  @eval function Base.$op(i1::ImageMetadata, i2::ImageMetadata)
      Base.$op(i1.updated, i2.updated)
  end
end

function extract_class(o::Dict)
    if haskey(o, "clinical")
        c = o["clinical"]
        if haskey(c, "benign_malignant")
            return Symbol(c["benign_malignant"])
        elseif haskey(c, "ben_mal")
            return Symbol(c["ben_mal"])
        end
    end
    :none
end

function Base.show(io::IO, o::ImageMetadata)
    print_with_color(:white, io, string(typeof(o)), "\n")
    print(io,   "  .name: ")
    print_with_color(:blue, io, o.name, "\n")
    println(io, "  .id: ", o.id)
    println(io, "  .modelType: ", o.modelType)
    #println(io, "  .creatorId: ", o.creatorId)
    println(io, "  .created: ", o.created)
    println(io, "  .updated: ", o.updated)
    println(io, "  .class: ", o.class)
    #println(io, "  .description: ", o.description)
    println(io, "")
    print(io,   "  .meta: parsed JSON as ")
    Base.showdict(io, o.meta)
end

# ====================================================================

@with_kw struct ImageMetadataRequest
    id::String
end

ImageMetadataRequest(le::Union{ListEntry,ImageMetadata}) = ImageMetadataRequest(id = le.id)

function Base.get(req::ImageMetadataRequest)
    query = "https://isic-archive.com:443/api/v1/image/$(req.id)"
    ImageMetadata(clean_json(get(query)))
end

# ====================================================================

@with_kw struct ImageDownloadRequest
    id::String
end

ImageDownloadRequest(le::Union{ListEntry,ImageMetadata}) = ImageDownloadRequest(id = le.id)

function Base.get(req::ImageDownloadRequest)
    query = "https://isic-archive.com:443/api/v1/image/$(req.id)/download"
    ImageMagick.readblob(get(query).data)
end

# ====================================================================

@with_kw struct ImageThumbnailRequest
    id::String
end

ImageThumbnailRequest(le::Union{ListEntry,ImageMetadata}) = ImageThumbnailRequest(id = le.id)

function Base.get(req::ImageThumbnailRequest)
    query = "https://isic-archive.com:443/api/v1/image/$(req.id)/thumbnail"
    ImageMagick.readblob(get(query).data)
end


text:
""" BFGS algorithm using a line search supporting stopping
    The initial B₀ may be a matrix, an identity or a operator.
    Using the operator implements L-BFGS
    kwargs may specify options for the line search.
"""
function bfgs_StopLS(nlp       :: AbstractNLPModel;
                     x         :: Vector{T}=copy(nlp.meta.x0),
                     stp       :: NLPStopping = NLPStopping(nlp,
                                                      NLPAtX(nlp.meta.x0)),
                     scaling   :: Bool = true,
                     LS_algo   :: Function = bracket,
                     LS_logger :: AbstractLogger = Logging.NullLogger(),
                     B₀        :: Union{AbstractLinearOperator,
                                        AbstractMatrix,
                                        UniformScaling{T},
                                        Nothing}              = nothing,
                     kwargs...      # eventually options for the line search
                     ) where T

    @info log_header([:iter, :f, :dual, :step, :slope], [Int, T, T, T, T],
                     hdr_override=Dict(:f=>"f(x)", :dual=>"‖∇f‖", :slope=>"∇fᵀd"))

    f = obj(nlp,x)
    ∇f = grad(nlp, x)

    xt = similar(x)
    ∇ft = similar(∇f)

    B = one(eltype(x))*I
    if B₀ != nothing
        B = B₀
    end

    τ₀ = 0.0005
    τ₁ = 0.9999

    ϕ, ϕstp = prepare_LS(stp, x, ∇f, τ₀, f, ∇f)

    OK = update_and_start!(stp, x = x, fx = f, gx = ∇f)
    update_and_stop!(stp,  x = x, fx = f, gx = ∇f)
    @info log_row(Any[0, f, norm(∇f)])

    while !OK
        d = - B*∇f

        # Simple line search call
        # returns  t, xt = x + t*d,  ft=f(xt), ∇ft = ∇f(xt)

        t, xt, ft, ∇ft = linesearch(ϕ, ϕstp, x, d, f, ∇f, τ₀, τ₁, logger = LS_logger, algo = LS_algo; kwargs...)
        fail_sub_pb = ~ϕstp.meta.optimal

        # Update BFGS approximation.
        B = push!(B, t * d, ∇ft - ∇f)
        #B = push!(B, t * d, ∇ft - ∇f, scaling)

        #move on
        x .= xt
        f = ft
        ∇f .= ∇ft

        if fail_sub_pb
            OK = true
            stp.meta.fail_sub_pb = true
        else
            OK = update_and_stop!(stp, x = x, gx = ∇f, fx = f)
        end

        @info log_row(Any[stp.meta.nb_of_stop, f, norm(∇f), t, ϕ.counters.neval_obj])
    end

    if !stp.meta.optimal
        @warn status(stp,list=true)
        @debug x
    end

    return tuple(stp, B)
end

""" L-BFGS wrapper of the gereral BFGS implementation.
"""
function L_bfgs_StopLS(nlp :: AbstractNLPModel;
                       x :: Vector{T}=copy(nlp.meta.x0),
                       stp :: NLPStopping = NLPStopping(nlp,
                                                        NLPAtX(nlp.meta.x0)),
                       mem :: Int = 5,
                       scaling :: Bool = true,
                       LS_algo   :: Function = bracket,
                       LS_logger :: AbstractLogger = Logging.NullLogger(),
                       kwargs...      # eventually options for the line search
                       ) where T

    @debug "U_Solver = L_bfgs_StopLS"
    n = nlp.meta.nvar
    B₀ =  InverseLBFGSOperator(Float64, n, mem=mem, scaling=scaling)

    return bfgs_StopLS(nlp; x=x, stp=stp, scaling=scaling, LS_logger=LS_logger, LS_algo=LS_algo, B₀=B₀, kwargs...)
end

max_stars_repo_name: tmigot/LSDescent

max_issues_repo_name: tmigot/LSDescent

max_forks_repo_name: tmigot/LSDescent


text:
using ArrayLayouts, LinearAlgebra, FillArrays, Base64, Test
import ArrayLayouts: sub_materialize, MemoryLayout, ColumnNorm, RowMaximum

struct MyMatrix <: LayoutMatrix{Float64}
    A::Matrix{Float64}
end

Base.getindex(A::MyMatrix, k::Int, j::Int) = A.A[k,j]
Base.setindex!(A::MyMatrix, v, k::Int, j::Int) = setindex!(A.A, v, k, j)
Base.size(A::MyMatrix) = size(A.A)
Base.strides(A::MyMatrix) = strides(A.A)
Base.unsafe_convert(::Type{Ptr{T}}, A::MyMatrix) where T = Base.unsafe_convert(Ptr{T}, A.A)
MemoryLayout(::Type{MyMatrix}) = DenseColumnMajor()
Base.copy(A::MyMatrix) = MyMatrix(copy(A.A))

struct MyVector <: LayoutVector{Float64}
    A::Vector{Float64}
end

Base.getindex(A::MyVector, k::Int) = A.A[k]
Base.setindex!(A::MyVector, v, k::Int) = setindex!(A.A, v, k)
Base.size(A::MyVector) = size(A.A)
Base.strides(A::MyVector) = strides(A.A)
Base.unsafe_convert(::Type{Ptr{T}}, A::MyVector) where T = Base.unsafe_convert(Ptr{T}, A.A)
MemoryLayout(::Type{MyVector}) = DenseColumnMajor()

# These need to test dispatch reduces to ArrayLayouts.mul, etc.
@testset "LayoutArray" begin
    @testset "LayoutVector" begin
        a = MyVector([1.,2,3])
        B = randn(3,3)
        b = randn(3)

        @test a == a.A == Vector(a)
        @test a[1:3] == a.A[1:3]
        @test a[:] == a
        @test (a')[1,:] == (a')[1,1:3] == a
        @test stringmime("text/plain", a) == "3-element MyVector:\n 1.0\n 2.0\n 3.0"
        @test B*a ≈ B*a.A
        @test B'*a ≈ B'*a.A
        @test transpose(B)*a ≈ transpose(B)*a.A
        @test b'a ≈ transpose(b)a ≈ a'b ≈ transpose(a)b ≈ b'a.A
        @test qr(B).Q*a ≈ qr(B).Q*a.A

        @test a'a == transpose(a)a == dot(a,a) == dot(a,a.A) == dot(a.A,a) == 14
        v = view(a,1:3)
        @test copy(v) == sub_materialize(v) == a[1:3]
        @test dot(v,a) == dot(v,a.A) == dot(a,v) == dot(a.A,v) == dot(v,v) == 14
        @test norm(v) == norm(a) == norm([1,2,3])

        V = view(a',:,1:3)
        @test copy(V) == sub_materialize(V) == (a')[:,1:3]

        s = SparseVector(3, [1], [2])
        @test a's == s'a == dot(a,s) == dot(s,a) == dot(s,a.A)
    end

    @testset "LayoutMatrix" begin
        A = MyMatrix(randn(5,5))
        for (kr,jr) in ((1:2,2:3), (:,:), (:,1:2), (2:3,:), ([1,2],3:4), (:,[1,2]), ([2,3],:),
                        (2,:), (:,2), (2,1:3), (1:3,2))
            @test A[kr,jr] == A.A[kr,jr]
        end
        b = randn(5)
        for Tri in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular)
            @test ldiv!(Tri(A), copy(b)) ≈ ldiv!(Tri(A.A), copy(b))
            @test lmul!(Tri(A), copy(b)) ≈ lmul!(Tri(A.A), copy(b))
        end

        @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), A) == A
        @test copyto!(Array{Float64}(undef,5,5), A) == A
        @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), A.A) == A
        @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),1:3,1:3), view(A,1:3,1:3)) == A[1:3,1:3]
        @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), A) == A
        @test copyto!(MyMatrix(Array{Float64}(undef,3,3)), view(A,1:3,1:3)) == A[1:3,1:3]
        @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), A.A) == A
        @test copyto!(Array{Float64}(undef,3,3), view(A,1:3,1:3)) == A[1:3,1:3]
        @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), A') == A'
        @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), view(A',:,:)) == A'
        @test copyto!(Array{Float64}(undef,5,5), A') == A'
        @test copyto!(Array{Float64}(undef,5,5), view(A',:,:)) == A'
        @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), A') == A'

        @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), view(A',:,:)) == A'

        @testset "factorizations" begin
            @test qr(A).factors ≈ qr(A.A).factors
            @test qr(A,ColumnNorm()).factors ≈ qr(A.A,ColumnNorm()).factors
            @test lu(A).factors ≈ lu(A.A).factors
            @test lu(A,RowMaximum()).factors ≈ lu(A.A,RowMaximum()).factors
            @test_throws ErrorException qr!(A)
            @test lu!(copy(A)).factors ≈ lu(A.A).factors
            b = randn(5)
            @test all(A \ b .≡ A.A \ b)
            @test all(lu(A).L .≡ lu(A.A).L)
            @test all(lu(A).U .≡ lu(A.A).U)
            @test lu(A).p == lu(A.A).p
            @test lu(A).P == lu(A.A).P

            @test qr(A) isa LinearAlgebra.QRCompactWY
            @test inv(A) ≈ inv(A.A)
            if VERSION ≥ v"1.7-"
                @test qr(A, Val(true)) == qr(A, ColumnNorm())
            end

            S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I))
            @test cholesky(S).U ≈ @inferred(cholesky!(deepcopy(S))).U
            @test cholesky(S,Val(true)).U ≈ cholesky(Matrix(S),Val(true)).U

            S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I),:L)
            @test cholesky(S).U ≈ @inferred(cholesky!(deepcopy(S))).U
            @test cholesky(S,Val(true)).U ≈ cholesky(Matrix(S),Val(true)).U
        end
        Bin = randn(5,5)
        B = MyMatrix(copy(Bin))
        muladd!(1.0, A, A, 2.0, B)
        @test all(B .=== A.A^2 + 2Bin)

        @testset "tiled_blasmul!" begin
            B = MyMatrix(copy(Bin))
            muladd!(1.0, Ones(5,5), A, 2.0, B)
        end

        @testset "mul!" begin
            @test mul!(B, A, A) ≈ A*A
            @test mul!(B, A', A) ≈ A'*A
            @test mul!(B, A, A') ≈ A*A'
            @test mul!(B, A', A') ≈ A'*A'
            @test mul!(B, A, Bin) ≈ A*Bin

            @test mul!(copy(B), A, A, 2, 3) ≈ 2A*A + 3B
            @test mul!(copy(B), A', A, 2, 3) ≈ 2A'*A + 3B
            @test mul!(copy(B), A, A', 2, 3) ≈ 2A*A' + 3B
            @test mul!(copy(B), A', A', 2, 3) ≈ 2A'*A' + 3B
            @test mul!(copy(B), Bin, A, 2, 3) ≈ 2Bin*A + 3B
            @test mul!(copy(B), Bin, A', 2, 3) ≈ 2Bin*A' + 3B
            @test mul!(copy(B), A, Bin, 2, 3) ≈ 2A*Bin + 3B
            @test mul!(copy(B), A', Bin, 2, 3) ≈ 2A'*Bin + 3B
            @test mul!(copy(B), A, Bin, 2, 3) ≈ 2A*Bin + 3B
            @test mul!(copy(B), A, Bin', 2, 3) ≈ 2A*Bin' + 3B
            @test mul!(copy(B), Bin', A, 2, 3) ≈ 2Bin'*A + 3B
        end

        @testset "generic_blasmul!" begin
            A = BigFloat.(randn(5,5))
            Bin = BigFloat.(randn(5,5))
            B = copy(Bin)
            muladd!(1.0, Ones(5,5), A, 2.0, B)
            @test B == Ones(5,5)*A + 2.0Bin
        end

        C = MyMatrix([1 2; 3 4])
        @test stringmime("text/plain", C) == "2×2 MyMatrix:\n 1.0  2.0\n 3.0  4.0"

        @testset "layoutldiv" begin
            A = MyMatrix(randn(5,5))
            x = randn(5)
            X = randn(5,5)
            t = view(randn(10),[1,3,4,6,7])
            T = view(randn(10,5),[1,3,4,6,7],:)
            t̃ = copy(t)
            T̃ = copy(T)
            B = Bidiagonal(randn(5),randn(4),:U)
            @test ldiv!(A, copy(x)) ≈ A\x
            @test A\t ≈ A\t̃
            # QR is not general enough
            @test_broken ldiv!(A, t) ≈ A\t
            @test ldiv!(A, copy(X)) ≈ A\X
            @test A\T ≈ A\T̃
            @test_broken A/T ≈ A/T̃
            @test_broken ldiv!(A, T) ≈ A\T
            @test B\A ≈ B\Matrix(A)
            @test transpose(B)\A ≈ transpose(B)\Matrix(A) ≈ Transpose(B)\A ≈ Adjoint(B)\A
            @test B'\A ≈ B'\Matrix(A)
            @test A\A ≈ I
            @test_broken A/A ≈ I
            @test A\MyVector(x) ≈ A\x
            @test A\MyMatrix(X) ≈ A\X
        end

        @testset "dot" begin
            A = MyMatrix(randn(5,5))
            b = randn(5)
            @test dot(b, A, b) ≈ b'*(A*b) ≈ b'A*b
            @test dot(b, A, b) ≈ transpose(b)*(A*b) ≈ transpose(b)A*b
        end

        @testset "dual vector * symmetric (#40)" begin
            A = randn(3,3)
            x = rand(3)
            @test x' * Symmetric(MyMatrix(A)) ≈ x'Symmetric(A)
            @test transpose(x) * Symmetric(MyMatrix(A)) ≈ transpose(x)Symmetric(A)
        end

        @testset "map(copy, ::Diagonal)" begin
            # this is needed in BlockArrays
            D = Diagonal([MyMatrix(randn(2,2)), MyMatrix(randn(2,2))])
            @test map(copy, D) == D
        end

        @testset "permutedims(::Diagonal)" begin
            D = Diagonal(MyVector(randn(5)))
            @test permutedims(D) ≡ D
        end
    end

    @testset "l/rmul!" begin
        b = MyVector(randn(5))
        A = MyMatrix(randn(5,5))
        @test lmul!(2, deepcopy(b)) == rmul!(deepcopy(b), 2) == 2b
        @test lmul!(2, deepcopy(A)) == rmul!(deepcopy(A), 2) == 2A
        @test lmul!(2, deepcopy(A)') == rmul!(deepcopy(A)', 2) == 2A'
        @test lmul!(2, transpose(deepcopy(A))) == rmul!(transpose(deepcopy(A)), 2) == 2transpose(A)
        @test lmul!(2, Symmetric(deepcopy(A))) == rmul!(Symmetric(deepcopy(A)), 2) == 2Symmetric(A)
        @test lmul!(2, Hermitian(deepcopy(A))) == rmul!(Hermitian(deepcopy(A)), 2) == 2Hermitian(A)

        C = randn(ComplexF64,5,5)
        @test ArrayLayouts.lmul!(2, Hermitian(copy(C))) == ArrayLayouts.rmul!(Hermitian(copy(C)), 2) == 2Hermitian(C)

        @test ldiv!(2, deepcopy(b)) == rdiv!(deepcopy(b), 2) == 2\b
        @test ldiv!(2, deepcopy(A)) == rdiv!(deepcopy(A), 2) == 2\A
        @test ldiv!(2, deepcopy(A)') == rdiv!(deepcopy(A)', 2) == 2\A'
        @test ldiv!(2, transpose(deepcopy(A))) == rdiv!(transpose(deepcopy(A)), 2) == 2\transpose(A)
        @test ldiv!(2, Symmetric(deepcopy(A))) == rdiv!(Symmetric(deepcopy(A)), 2) == 2\Symmetric(A)
        @test ldiv!(2, Hermitian(deepcopy(A))) == rdiv!(Hermitian(deepcopy(A)), 2) == 2\Hermitian(A)
        @test ArrayLayouts.ldiv!(2, Hermitian(copy(C))) == ArrayLayouts.rdiv!(Hermitian(copy(C)), 2) == 2\Hermitian(C)
    end

    @testset "pow/I" begin
        A = randn(2,2)
        B = MyMatrix(A)
        @test B^2 ≈ A^2
        @test B^2.3 ≈ A^2.3
        @test B^(-1) ≈ inv(A)
        @test B + I ≈ I + B ≈ A + I
        @test B - I ≈ A - I
        @test I - B ≈ I - B
    end

    @testset "Diagonal" begin
        D = Diagonal(MyVector(randn(5)))
        D̃ = Diagonal(Vector(D.diag))
        B = randn(5,5)
        B̃ = MyMatrix(B)
        @test D*D ≈ Matrix(D)^2
        if VERSION ≥ v"1.7-"
            @test D^2 ≈ D*D
        end
        @test D*B ≈ Matrix(D)*B
        @test B*D ≈ B*Matrix(D)
        @test D*B̃ ≈ Matrix(D)*B̃
        @test B̃*D ≈ B̃*Matrix(D)
        @test D*D̃ ≈ D̃*D

        @test D\D ≈ I
        @test D\B ≈ Matrix(D)\B
        @test B\D ≈ B\Matrix(D)
        @test D\B̃ ≈ Matrix(D)\B̃
        @test B̃\D ≈ B̃\Matrix(D)
        @test D\D̃ ≈ D̃\D
    end

    @testset "Adj/Trans" begin
        A = MyMatrix(randn(5,5))
        T = UpperTriangular(randn(5,5))
        D = Diagonal(MyVector(randn(5)))

        @test D * A' ≈ D * A.A'
        @test A' * D ≈ A.A' * D

        @test A * Adjoint(T) ≈ A.A * Adjoint(T)
        @test A * Transpose(T) ≈ A.A * Transpose(T)
        @test Adjoint(T) * A ≈ Adjoint(T) * A.A
        @test Transpose(T) * A ≈ Transpose(T) * A.A
        @test Transpose(T)A' ≈ Adjoint(T)A' ≈ Adjoint(T)Transpose(A) ≈ Transpose(T)Transpose(A)
        @test Transpose(A)Adjoint(T) ≈ A'Adjoint(T) ≈ A'Transpose(T) ≈ Transpose(A)Transpose(T)

        @test Zeros(5)' * A ≡ Zeros(5)'
        @test transpose(Zeros(5)) * A ≡ transpose(Zeros(5))

        @test A' * Zeros(5) ≡ Zeros(5)
        @test Zeros(3,5) * A ≡ Zeros(3,5)
        @test A * Zeros(5,3) ≡ Zeros(5,3)
        @test A' * Zeros(5,3) ≡ Zeros(5,3)
        @test transpose(A) * Zeros(5,3) ≡ Zeros(5,3)
        @test A' * Zeros(5) ≡ Zeros(5)
        @test transpose(A) * Zeros(5) ≡ Zeros(5)

        b = MyVector(randn(5))
        @test A' * b ≈ A' * b.A
    end

    @testset "AbstractQ" begin
        A = MyMatrix(randn(5,5))
        Q = qr(randn(5,5)).Q
        @test Q'*A ≈ Q'*A.A
        @test Q*A ≈ Q*A.A
        @test A*Q ≈ A.A*Q
        @test A*Q' ≈ A.A*Q'
    end

    @testset "concat" begin
        a = MyVector(randn(5))
        A = MyMatrix(randn(5,6))

        @test [a; a] == [Array(a); a] == [a; Array(a)]
        @test [A; A] == [Array(A); A] == [A; Array(A)]

        @test [Array(a) A A] == [Array(a) Array(A) Array(A)]
        @test [a A A] == [Array(a) Array(A) Array(A)]
        @test [a Array(A) A] == [Array(a) Array(A) Array(A)]
    end

    @testset "dot" begin
        a = MyVector(randn(5))
        @test dot(a, Zeros(5)) ≡ dot(Zeros(5), a) ≡ 0.0
    end
end

struct MyUpperTriangular{T} <: AbstractMatrix{T}
    A::UpperTriangular{T,Matrix{T}}
end

MyUpperTriangular{T}(::UndefInitializer, n::Int, m::Int) where T = MyUpperTriangular{T}(UpperTriangular(Array{T}(undef, n, m)))
MyUpperTriangular(A::AbstractMatrix{T}) where T = MyUpperTriangular{T}(UpperTriangular(Matrix{T}(A)))
Base.convert(::Type{MyUpperTriangular{T}}, A::MyUpperTriangular{T}) where T = A
Base.convert(::Type{MyUpperTriangular{T}}, A::MyUpperTriangular) where T = MyUpperTriangular(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyUpperTriangular}, A::MyUpperTriangular)= A
Base.convert(::Type{AbstractArray{T}}, A::MyUpperTriangular) where T = MyUpperTriangular(convert(AbstractArray{T}, A.A))
Base.convert(::Type{AbstractMatrix{T}}, A::MyUpperTriangular) where T = MyUpperTriangular(convert(AbstractArray{T}, A.A))
Base.convert(::Type{MyUpperTriangular{T}}, A::AbstractArray{T}) where T = MyUpperTriangular{T}(A)
Base.convert(::Type{MyUpperTriangular{T}}, A::AbstractArray) where T = MyUpperTriangular{T}(convert(AbstractArray{T}, A))
Base.convert(::Type{MyUpperTriangular}, A::AbstractArray{T}) where T = MyUpperTriangular{T}(A)
Base.getindex(A::MyUpperTriangular, kj...) = A.A[kj...]
Base.getindex(A::MyUpperTriangular, ::Colon, j::AbstractVector) = MyUpperTriangular(A.A[:,j])
Base.setindex!(A::MyUpperTriangular, v, kj...) = setindex!(A.A, v, kj...)
Base.size(A::MyUpperTriangular) = size(A.A)
Base.similar(::Type{MyUpperTriangular{T}}, m::Int, n::Int) where T = MyUpperTriangular{T}(undef, m, n)
Base.similar(::MyUpperTriangular{T}, m::Int, n::Int) where T = MyUpperTriangular{T}(undef, m, n)
Base.similar(::MyUpperTriangular, ::Type{T}, m::Int, n::Int) where T = MyUpperTriangular{T}(undef, m, n)
LinearAlgebra.factorize(A::MyUpperTriangular) = factorize(A.A)

MemoryLayout(::Type{MyUpperTriangular{T}}) where T = MemoryLayout(UpperTriangular{T,Matrix{T}})
triangulardata(A::MyUpperTriangular) = triangulardata(A.A)

@_layoutlmul MyUpperTriangular

@testset "MyUpperTriangular" begin
    A = randn(5,5)
    B = randn(5,5)
    x = randn(5)
    U = MyUpperTriangular(A)

    @test lmul!(U, copy(x)) ≈ U * x
    @test lmul!(U, copy(B)) ≈ U * B

    @test_skip lmul!(U,view(copy(B),collect(1:5),1:5)) ≈ U * B
end


text:
using UtilitiesForMRI, Random

# Cartesian domain
n = (128,128,128)
h = (1f0, 1f0, 1f0)
X = spatial_sampling(Float32, n; h=h)

# Cartesian sampling (dense)
phase_encoding = (1,2)
K = kspace_Cartesian_sampling(X; phase_encoding=phase_encoding)
# using PyPlot
# figure()
# plot3D(K[1][:,1],  K[1][:,2],  K[1][:,3], ".")
# plot3D(K[10][:,1], K[10][:,2], K[10][:,3], ".")
# plot3D(K[20][:,1], K[20][:,2], K[20][:,3], ".")
# plot3D(K[128*(10-1)+1][:,1], K[128*(10-1)+1][:,2], K[128*(10-1)+1][:,3], ".")
# plot3D(K[128*(20-1)+1][:,1], K[128*(20-1)+1][:,2], K[128*(20-1)+1][:,3], ".")
# plot3D(K[128^2][:,1], K[128^2][:,2], K[128^2][:,3], ".")

# Cartesian sampling (randomized)
subsampling = (1:128^2)[randperm(128^2)][1:32]
K = kspace_Cartesian_sampling(X; phase_encoding=(1,3), subsampling=subsampling)
# figure()
# for i = 1:32
#     plot3D(K[i][:,1], K[i][:,2], K[i][:,3], ".")
# end

# for i = 1:50:128^2
#     plot3D(K[i][:,1], K[i][:,2], K[i][:,3], ".")
# end

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
"""
    lambdipole(U, R, grid::TwoDGrid; center=(mean(grid.x), mean(grid.y))

Return the two-dimensional vorticity field of the Lamb dipole with strength `U` and radius `R`,
centered on `center=(xc, yc)` and on the `grid`. Default value for `center` is the center of
the domain.
"""
function lambdipole(U, R, grid::TwoDGrid{T, A}; center=(mean(grid.x), mean(grid.y))) where {T, A}
  firstzero = 3.8317059702075123156
  k = firstzero / R # dipole wavenumber for radius R in terms of first zero of besselj
  q₀ = -2U * k / besselj(0, k * R) # dipole amplitude for strength U and radius R
  x, y = gridpoints(grid)
  xc, yc = center
  r = @. sqrt((x - xc)^2 + (y - yc)^2)
  besseljorder1 = CUDA.@allowscalar A([besselj1(k * r[i, j]) for i=1:grid.nx, j=1:grid.ny])
  q = @. q₀ * besseljorder1 * (y - yc) / r
  @. q = ifelse(r >= R, 0, q)

  return q
end

"""
    peakedisotropicspectrum(grid, kpeak, E₀; mask=mask, allones=false)

Generate a random two-dimensional relative vorticity field ``q(x, y)`` with Fourier spectrum
peaked around a central non-dimensional wavenumber `kpeak` and normalized so that its total
kinetic energy is `E₀`.
"""
function peakedisotropicspectrum(grid::TwoDGrid{T, A}, kpeak::Real, E₀::Real; mask=ones(size(grid.Krsq)), allones=false) where {T, A}
  if grid.Lx !== grid.Ly
    error("the domain is not square")
  else
    k₀ = kpeak * 2π / grid.Lx
    modk = sqrt.(grid.Krsq)
    modψ = A(zeros(T, (grid.nk, grid.nl)))
    modψ = @. (modk^2 * (1 + (modk / k₀)^4))^(-0.5)
    CUDA.@allowscalar modψ[1, 1] = 0.0

    phases = randn(Complex{T}, size(grid.Krsq))
    phases_real, phases_imag = real.(phases), imag.(phases)
    phases = A(phases_real) + im * A(phases_imag)
    ψh = @. phases * modψ
    if allones; ψh = modψ; end
    ψh = ψh .* A(mask)
    energy_initial = sum(grid.Krsq .* abs2.(ψh)) / (grid.nx * grid.ny)^2
    ψh *= sqrt(E₀ / energy_initial)
    q = A(irfft(-grid.Krsq .* ψh, grid.nx))
  end

  return q
end