text:
# # Custom workaround OrdinaryDiffEq: https://github.com/SciML/OrdinaryDiffEq.jl/blob/54fb35870fa402fc95d665cd5f9502e2759ea436/src/derivative_utils.jl#L606-L608
# # as proposed in: https://github.com/CliMA/ClimaCore.jl/blob/af6b01bc8f945ada0794147d8f2c7ca5303adf62/examples/hybrid/inertial_gravity_wave_implicit.jl#L411-L414
# struct CustomWFact
#     W
#     γref
# end
# Base.similar(W::CustomWFact) = deepcopy(W)
# function linsolve!(::Type{Val{:init}}, f, u0; kwargs...)
#     function _linsolve!(x, W, b, update_matrix = false; kwargs...)
#         W isa CustomWFact || error("Linsolve specialized for CustomWFact")
#         ldiv!(x, W.W, b)
#     end
# end

"""
    solve_btpde(model, matrices, experiment)

Solve the Bloch-Torrey partial differential equation using P1 finite elements.
"""
function solve_btpde(model::Model, matrices, experiment::Experiment)

    # Measure function evalutation time
    starttime = Base.time()

    # Extract input parameters
    @unpack mesh, D, T₂, ρ = model
    @unpack M, S, R, Mx, Q, M_cmpts = matrices
    @unpack directions, sequences, values, values_type = experiment.gradient
    @unpack reltol, abstol, nsave = experiment.btpde

    # Deduce sizes
    ncompartment = length(ρ)
    npoint_cmpts = size.(mesh.points, 2)
    inds_cmpts = cumsum([0; npoint_cmpts])
    ndirection = size(directions, 2)
    nsequence = length(sequences)
    namplitude = length(values)

    qvalues, bvalues = get_values(experiment.gradient)

    # Create initial conditions (enforce complex values)
    ρ = vcat(fill.(complex.(ρ), npoint_cmpts)...)

    # Allocate output arrays
    signal = zeros(ComplexF64, ncompartment, namplitude, nsequence, ndirection)
    signal_allcmpts = zeros(ComplexF64, namplitude, nsequence, ndirection)
    magnetization =
        Array{Matrix{ComplexF64},4}(undef, ncompartment, namplitude, nsequence, ndirection)
    time = Array{Vector{Float64},3}(undef, namplitude, nsequence, ndirection)
    itertimes = zeros(namplitude, nsequence, ndirection)

    # Time dependent ODE function
    function Mdξ!(dξ, ξ, p, t)
        # @show t
        @unpack J, S, Q, A, R, q, f = p
        @. J = -(S + Q + R + im * f(t) * q * A)
        mul!(dξ, J, ξ)
        nothing
    end

    # Time independent ODE function, given Jacobian
    function Mdξ_constant!(dξ, ξ, p, t)
        # @show t
        J = p.J
        mul!(dξ, J, ξ)
        nothing
    end

    # Time dependent Jacobian of ODE function with respect to the state `ξ`
    function Jac!(J, ξ, p, t)
        # println("Jac at $t")
        @unpack S, Q, A, R, q, f = p
        @. J = -(S + Q + R + im * f(t) * q * A)
        nothing
    end

    # # TODO: propagate types
    # # TODO: look for caching packages
    # # TODO: detect change in f and recompute J and W
    # function Wfact(u, p, γ, t)
    #     println("Wfact at $t")
    #     @unpack J = p
    #     CustomWFact(lu(M .- γ .* J), Ref(0.0))
    # end
    # function Wfact!(W, u, p, γ, t)
    #     # println("Wfact! at $t")
    #     @unpack J, γ_cache, factorization_cache = p
    #     ind = findfirst(isapprox(γ), γ_cache)
    #     # if isnothing(ind)
    #     if !(γ ≈ W.γref[])
    #         println("Refact at γ = $γ, t = $t")
    #         # W = factorize(M .- γ .* J)
    #         lu!(W.W, M .- γ .* J)
    #         W.γref[] = γ
    #     else
    #         println("No fact at γ = $γ, t = $t")
    #     end
    #     #     push!(γ_cache, γ)
    #     #     push!(factorization_cache, W)
    #     # else
    #     #     W = factorization_cache[ind]
    #     # end
    #     W
    # end

    # Jacobian sparsity pattern
    jac_prototype = -(S + Q + im * sum(Mx))

    # Iterate over gradient amplitudes, time profiles and directions
    for idir = 1:ndirection, iseq = 1:nsequence, iamp = 1:namplitude

        # Gradient amplitude
        q = qvalues[iamp, iseq]
        b = bvalues[iamp, iseq]

        # Time profile
        f = sequences[iseq]
        interval = (0, echotime(f))
        if nsave == 1
            saveat = [interval[2]]
            # saveat = [interval[1], interval[2]]
        else
            saveat = LinRange(interval..., nsave)
        end

        # Gradient direction
        dir = directions[:, idir]

        # Display state of iterations
        @printf "Solving BTPDE with size %d\n" (sum(npoint_cmpts))
        @printf "  Direction %d of %d: g = [%.2f, %.2f, %.2f]\n" idir ndirection dir...
        @printf "  Sequence  %d of %d: f = %s\n" iseq nsequence f
        @printf "  Amplitude %d of %d: q = %g, b = %g\n" iamp namplitude q b

        # Gradient direction dependent finite element matrix
        A = dir' * Mx

        # ODE problem
        J = -(S + Q + im * q * A)
        # γ_cache = Float64[]
        # factorization_cache = LinearAlgebra.Factorization[]
        # p = (; J, S, Q, A, R, q, f, γ_cache, factorization_cache)
        p = (; J, S, Q, A, R, q, f)

        # Gather ODE function
        if all(constant_intervals(p.f)) # is_constant(p.f, t)
            func = Mdξ_constant!
            Jac!(p.J, ρ, p, 0)
            jac = (J, _, p, t) -> (J .= p.J; nothing)
            # jac = (J, _, p, t) -> (println("Jac called at $t"); J .= p.J; nothing)
        else
            func = Mdξ!
            jac = Jac!
        end
        odefunction = ODEFunction(
            func;
            mass_matrix = M,
            jac,
            jac_prototype, # = Wfact(ρ, p, 0.1, 0.0),
            # Wfact = Wfact!,
        )
        odeproblem = ODEProblem(odefunction, ρ, interval, p, progress = false)#, dtmax = 50)

        tstops = intervals(f)[2:end-1]

        # Callback for updating problem
        function affect!(integrator)
            @unpack u, p, t = integrator
            # print("    t = $t: constant time profile, f = $(p.f(t))")
            Jac!(p.J, u, p, t)
        end
        callback = PresetTimeCallback(tstops, affect!)

        # Measure solving time
        itertime = Base.time()

        # Solve ODE problem
        sol = solve(
            odeproblem,
            QNDF(); # , linsolve = linsolve!);
            saveat,
            reltol,
            abstol,
            callback,
        )

        itertimes[iamp, iseq, idir] = Base.time() - itertime

        # Extract solution
        time[iamp, iseq, idir] = sol.t
        ξ = hcat(sol.u...)

        if nsave == 1
            time[iamp, iseq, idir] = time[iamp, iseq, idir][[end]]
            ξ = ξ[:, [end]]
        end

        # Split solution into compartments
        for icmpt = 1:ncompartment
            inds = inds_cmpts[icmpt]+1:inds_cmpts[icmpt+1]

            # Store magnetization in compartment
            magnetization[icmpt, iamp, iseq, idir] = ξ[inds, :]

            # Integrate final magnetization over compartment
            signal[icmpt, iamp, iseq, idir] =
                sum(M_cmpts[icmpt] * ξ[inds, end], dims = 1)[1]

        end
    end

    signal_allcmpts = sum(signal, dims = 1)[1, :, :, :]

    totaltime = Base.time() - starttime

    (; magnetization, signal, signal_allcmpts, time, itertimes, totaltime)
end


text:
function wordcloud(; x::AbstractVector =  String[], wordAngles::AbstractVector = [-45, 0, 45], minThreshold::Int = 0)

    v = VegaVisualization(height = 400,
                      width = 800,
                      padding = "auto",
                      name = "wordcloud",
                      data = Array{VegaData}(1),
                      scales = Array{VegaScale}(1),
                      marks = Array{VegaMark}(1))

    #Make unique table name
    table = string("table_", lowercase(randstring(5)))

    v.data[1] = VegaData(name = table,
                         values = x,
                         transform = Array{VegaTransform}(6))

    v.data[1].transform[1] = VegaTransform(Dict{Any, Any}("type" => "countpattern",
                                                        "field" => "data",
                                                        "case" => "upper",
                                                        "pattern" => "[\\w']{3,}",
                                                        "stopwords" => "(i|me|my|myself|we|us|our|ours|ourselves|you|your|yours|yourself|yourselves|he|him|his|himself|she|her|hers|herself|it|its|itself|they|them|their|theirs|themselves|what|which|who|whom|whose|this|that|these|those|am|is|are|was|were|be|been|being|have|has|had|having|do|does|did|doing|will|would|should|can|could|ought|i'm|you're|he's|she's|it's|we're|they're|i've|you've|we've|they've|i'd|you'd|he'd|she'd|we'd|they'd|i'll|you'll|he'll|she'll|we'll|they'll|isn't|aren't|wasn't|weren't|hasn't|haven't|hadn't|doesn't|don't|didn't|won't|wouldn't|shan't|shouldn't|can't|cannot|couldn't|mustn't|let's|that's|who's|what's|here's|there's|when's|where's|why's|how's|a|an|the|and|but|if|or|because|as|until|while|of|at|by|for|with|about|against|between|into|through|during|before|after|above|below|to|from|up|upon|down|in|out|on|off|over|under|again|further|then|once|here|there|when|where|why|how|all|any|both|each|few|more|most|other|some|such|no|nor|not|only|own|same|so|than|too|very|say|says|said|shall)"
                                                        ))

    #Calculate angles
    v.data[1].transform[3] = VegaTransform(Dict{Any, Any}("type" => "formula",
                                                            "field" => "angle",
                                                            "expr" => "$(wordAngles)[~~(random() * $(length(wordAngles)))]"
                                                            ))

    #Specify wordcloud properties
    v.data[1].transform[4] = VegaTransform(Dict{Any, Any}("type" => "wordcloud",
                                                            "size" => [800,400],
                                                            "text" => Dict{Any, Any}("field" => "text"),
                                                            "rotate" => Dict{Any, Any}("field" => "angle"),
                                                            "font" => Dict{Any, Any}("value" => "Helvetica Neue"),
                                                            "fontSize" => Dict{Any, Any}("field" => "count"),
                                                            "fontWeight" => Dict{Any, Any}("field" => "weight"),
                                                            "fontScale" => [12,56]
                                                            ))

    #Set weights
    v.data[1].transform[2] = VegaTransform(Dict{Any, Any}("type" => "formula",
    "field" => "weight",
    "expr" => "if(datum.text=='VEGA', 600, 300)"

                                                         ))

    #Set min word count
    v.data[1].transform[5] = VegaTransform(Dict{Any, Any}("type" => "filter",
    "test" => "datum.count > $(minThreshold)"

                                                         ))

    #Create field "group" for consistency
    v.data[1].transform[6] = VegaTransform(Dict{Any, Any}("type" => "formula",
    "field" => "group",
    "expr" => "datum.text"

                                                         ))

    #Define color palette
    v.scales[1] = VegaScale(name = "group", _type = "ordinal", range = ["#d5a928", "#652c90", "#939597"])

    #Draw words
    v.marks[1] = VegaMark(_type = "text",
                          from = VegaMarkFrom(data = table),
                          properties = VegaMarkProperties(enter = VegaMarkPropertySet(
                                                                                x = VegaValueRef(field = "layout_x"),
                                                                                y = VegaValueRef(field = "layout_y"),
                                                                                angle = VegaValueRef(field = "layout_rotate"),
                                                                                font = VegaValueRef(field = "layout_font"),
                                                                                fontSize = VegaValueRef(field = "layout_fontSize"),
                                                                                fontStyle = VegaValueRef(field = "layout_fontStyle"),
                                                                                fontWeight = VegaValueRef(field = "layout_fontWeight"),
                                                                                text = VegaValueRef(field = "text"),
                                                                                align = VegaValueRef(value = "center"),
                                                                                baseline = VegaValueRef(value = "alphabetic"),
                                                                                fill = VegaValueRef(scale = "group", field = "text")
                                                                                )

                                                          )
                         )

    #Default to Paired color scale, 12
    colorscheme!(v; palette = ("Paired", 12))
    return v
end

ext: jl


text:
using MAT
cd(dirname(Base.source_path()))
include("../src/utils.jl")
include("../src/2d/cpab2dIntegration.jl")
include("../src/2d/imageWarp.jl")

file = matopen("tri_R2toR2_MINS_0_0_MAXS_512_512_nc_4_4_Xbdr_0_0_ext_0.mat")
As = read(file, "As")
pts_src = read(file, "pts_src")
tw_args = read(file, "tw_args")
nCols = tw_args["nCols"]
nRows = tw_args["nRows"]
close(file)

nC1 = 4
nC2 = 4
incX = nCols/nC1
incY = nRows/nC2
nPts = size(pts_src,1)
t=0.001
nSteps = 100
Nsteps = 10
Δt = t/Nsteps
ΔtAs = Δt*As
solutionSingleThread = Array(Float64, nPts, 2)

@time solver2d(pts_src,solutionSingleThread,t,As,nCols,nRows,incX,incY,Nsteps,Δt,nSteps,ΔtAs)

using PyPlot
PyPlot.figure("CPAB 2D Integration",figsize=(12,6))
PyPlot.subplot(1,2,1)
PyPlot.title("Original Image")
img = load("lena.jpg")
imr = reinterpret(UInt8, separate(img))
PyPlot.imshow(imr.data)

PyPlot.subplot(1,2,2)
PyPlot.title("Destination Image")
newimg1 = reinterpret(UInt8, separate(imwarp(solutionSingleThread, img, nCols, nRows)))
PyPlot.imshow(newimg1.data)


text:
#=
Note that the following styles used values from the matplotlib style library
(https://github.com/matplotlib/matplotlib/tree/master/lib/matplotlib/mpl-data/stylelib):

- ggplot
- fivethirtyeight
- seaborn

=#

immutable Style
    color_cycle::Vector
    layout::Layout
    global_trace::PlotlyAttribute
    trace::Dict{Symbol,PlotlyAttribute}
end

function Style(;color_cycle=[], layout=Layout(), global_trace=attr(),
                trace=Dict{Symbol,PlotlyAttribute}())
    Style(color_cycle, layout, global_trace, trace)
end

function Style(ps1::Style, ps2::Style)
    cs = isempty(ps2.color_cycle) ? ps1.color_cycle: ps2.color_cycle

    la = deepcopy(ps1.layout)
    for (k, v) in ps2.layout.fields
        la[k] = v
    end

    gta = merge(ps1.global_trace, ps2.global_trace)

    ta = deepcopy(ps1.trace)
    for (k, v) in ps2.trace
        ta_k = get(ta, k, attr())
        merge!(ta_k, v)
        ta[k] = ta_k
    end

    Style(color_cycle=cs, layout=la, global_trace=gta, trace=ta)
end

Style(pss::Style...) = foldl(Style, pss[1], pss[2:end])

function Style(base::Style; color_cycle=[], layout=Layout(),
               global_trace=attr(), trace=Dict{Symbol,PlotlyAttribute}())
    new_style = Style(color_cycle, layout, global_trace, trace)
    Style(base, new_style)
end

@compat function Base.show(io::IO, ::MIME"text/plain", s::Style)
    println(io, "Style with:")

    if !isempty(s.color_cycle)
        println(io, "  - color_cycle: ", s.color_cycle)
    end

    if !isempty(s.layout)
        print(io, "  - "); show(io, MIME"text/plain"(), s.layout)
    end

    if !isempty(s.global_trace)
        print(io, "  - global_trace: ")
        show(io, MIME"text/plain"(), s.global_trace)
    end

    if !isempty(s.trace)
        println(io, "  - trace: ")
        for (k, v) in s.trace
            print(io, "    - ", k, ": "); show(io, MIME"text/plain"(), v)
        end
    end
end

const DEFAULT_STYLE = [Style()]
const STYLES = [:default, :ggplot, :fivethirtyeight, :seaborn]

function ggplot_style()
    axis = attr(showgrid=true, gridcolor="white", linewidth=1.0,
                linecolor="white", titlefont_color="#555555",
                titlefont_size=14, ticks="outside",
                tickcolor="#555555"
                )
    layout = Layout(plot_bgcolor="#E5E5E5",
                    paper_bgcolor="white",
                    font_size=10,
                    xaxis=axis,
                    yaxis=axis,
                    titlefont_size=14)

    gta = attr(marker_line_width=0.5, marker_line_color="#348ABD")

    colors = ["#E24A33", "#348ABD", "#988ED5", "#777777", "#FBC15E",
              "#8EBA42", "#FFB5B8"]
    Style(layout=layout, color_cycle=colors, global_trace=gta)
end

function fivethirtyeight_style()
    ta = Dict(:scatter=>attr(line_width=4))
    axis = attr(showgrid=true, gridcolor="#cbcbcb", linewidth=1.0,
                      ticklen=0.0,
                      linecolor="#f0f0f0", titlefont_color="#555555",
                      titlefont_size=12, ticks="outside", showgrid=true,
                      tickcolor="#555555"
                      )
    layout = Layout(plot_bgcolor="#f0f0f0",
                    paper_bgcolor="#f0f0f0",
                    font_size=14,
                    xaxis=axis,
                    yaxis=axis,
                    legend=attr(borderwidth=1.0,
                                bgcolor="f0f0f0", bordercolor="f0f0f0"),
                    titlefont_size=14)
    colors = ["#008fd5", "#fc4f30", "#e5ae38", "#6d904f",
              "#8b8b8b", "#810f7c"]
    Style(layout=layout, color_cycle=colors, trace=ta)
end

function seaborn_style()
    ta = Dict(:heatmap=>attr(colorscale="Greys"),
              :scatter=>attr(marker=attr(size=7, line_width=0.0),
                             line_width=1.75))
    axis = attr(showgrid=true, gridcolor="white", linewidth=1.0,
                      ticklen=0.0, showgrid=true,
                      linecolor="white", titlefont_color="#555555",
                      titlefont_size=12, ticks="outside",
                      tickfont_size=10,
                      tickcolor="#555555"
                      )
    # TODO: no concept of major vs minor ticks...
    layout = Layout(plot_bgcolor="EAEAF2",
                    paper_bgcolor="white",
                    width=800,
                    height=550, # TODO: what does font_color=0.15 mean??
                    font=attr(family="Arial", size=14, color=0.15),
                    xaxis=axis,
                    yaxis=axis,
                    legend=attr(font_size=10,
                                bgcolor="white", bordercolor="white"),
                    titlefont_size=14)
    colors = ["#4C72B0", "#55A868", "#C44E52", "#8172B2", "#CCB974", "#64B5CD"]
    Style(color_cycle=colors, trace=ta, layout=layout)
end

reset_style!() = DEFAULT_STYLE[1] = Style()
use_style!(sty::Symbol) = DEFAULT_STYLE[1] = style(sty)
use_style!(s::Style) = DEFAULT_STYLE[1] = s

function style(sty::Symbol)
    sty == :ggplot ? ggplot_style() :
    sty == :fivethirtyeight ? fivethirtyeight_style() :
    sty == :seaborn ? seaborn_style() :
    sty == :default ? Style() :
    error("Uknown style $sty")
end


text:
# GapFindFill Interface

#=
Set up MILP model for Gap-Find
Input:
  isRev: true if corresponding reaction is reversible;
  isCyt: true if corresponding metabolite is in cytosol;
  isExt: true if corresponding metabolite is in extracellular compartment;
  stoiMatrix: stoichiometric matrix, #compounds * #reactions
  fluxLB: flux lower bound;
  fluxUB: flux upper bound;
  solver: CPLEX, Gurobi, GLPK;
  epsilon: minimum amount to be considered active;
  bigM: constant used in MILP model;
  nonZero: minimum stoichiometric coefficient to be considered non-zero.
Output:
  m: the JuMP model;
  objVal: objective value, i.e. # non-blocked compounds;
  status: termination status;
  binX: 1 if corresponding compound is non-blocked;
=#
function find_gaps(isRev::Array{Bool}, isCyt::Array{Bool}, isExt::Array{Bool},
    stoiMatrix::Array{Float64}, fluxLB::Array{Float64}, fluxUB::Array{Float64};
    epsilon::Float64 =0.001, bigM::Float64 =1000.0, nonZero::Float64 =1e-8,
    solver::Module=Gurobi)
  t1 = time()
  noCompound, noRxn = size(stoiMatrix)
  m = Model(with_optimizer(solver.Optimizer))
  # define variables using @variable
  @variable(m, fluxLB[i] <= v[i=1:noRxn] <= fluxUB[i])  # fluxes
  @variable(m, x[1:noCompound], Bin)  # 1 if metabolite i could be produced
  @variable(m, w[1:noCompound, 1:noRxn], Bin)  # 1 if metabolite i is produced in reaction j
  # define objective
  @objective(m, Max, sum(x))
  # define constraints
  # Mass Balance
  for i = 1:noCompound
    if (isCyt[i] && !isExt[i])  # cytosl but not extracellular
      @constraint(m, sum(stoiMatrix[i,j]*v[j] for j=1:noRxn) >= 0)
    elseif (!isCyt[i] && !isExt[i])  # other compartments but not extracellular
      @constraint(m, sum(stoiMatrix[i,j]*v[j] for j=1:noRxn) == 0)
    end
  end
  # Production Constrains for (ir)reversible reactions
  for i = 1:noCompound
    for j = 1:noRxn
      if (abs(stoiMatrix[i,j]) > nonZero)  # non-zero coefficient
        @constraint(m, stoiMatrix[i,j]*v[j] >= epsilon - bigM*(1-w[i,j]))
        @constraint(m, stoiMatrix[i,j]*v[j] <= bigM*w[i,j])
      end
    end
  end
  # Binary Constraints
  for i = 1:noCompound
    @constraint(m, sum(w[i,j] for j=1:noRxn if (
      (!isRev[j] && (stoiMatrix[i,j] > nonZero)) ||
      (isRev[j] && (abs(stoiMatrix[i,j]) > nonZero)))) >= x[i])
  end
  # Calling solver
  JuMP.optimize!(m)
  t2 = time()
  # get some results
  objVal = objective_value(m)
  status = termination_status(m)
  solveTime = t2 - t1
  binX = value.(x)
  return m, objVal, status, solveTime, binX
end

#=
Set up MILP model for Gap-Fill
Input:
  isMd: true if corresponding reaction is in the model;
  isDb: true if corresponding reaction is in the database,
        i.e. not in the model;
  isRev: true if corresponding reaction is reversible;
  isCyt: true if corresponding metabolite is in cytosol;
  isExt: true if corresponding metabolite is in extracellular compartment;
  noProdID: indices of no-production-metabolites
  stoiMatrix: stoichiometric matrix, #compounds * #reactions
  fluxLB: flux lower bound;
  fluxUB: flux upper bound;
  epsilon: minimum amount to be considered active;
  bigM: constant used in MILP model;
  nonZero: minimum stoichiometric coefficient to be considered non-zero.
Output:
  results: a dictionary, each key is an index of no-production-mebanolites,
    each value is a Set containing all possible gap filling ways with minimum
    number of added reactions.
=#
function fill_gaps_min(isMd::Array{Bool}, isDb::Array{Bool}, isRev::Array{Bool},
    isCyt::Array{Bool}, isExt::Array{Bool}, noProdID::Array{Integer},
    stoiMatrix::Array{Float64}, fluxLB::Array{Float64}, fluxUB::Array{Float64};
    epsilon::Float64 =0.001, bigM::Float64 =1000.0, nonZero::Float64 =1e-9)
  solver=Gurobi
  results = Dict{Int,Any}()
  for i in noProdID
    m, objVal, status, solveTime, yresult = fix_metabolite(
        isMd, isDb, isRev, isCyt, isExt, i,
        stoiMatrix, fluxLB, fluxUB,
        solver, epsilon, bigM, nonZero)
    println(objVal)
    println(yresult)
    # store result anyway (maybe no result)
    if haskey(results, i)
      push!(results[i], yresult)
    else
      results[i] = Set()
      push!(results[i], yresult)
    end
    if objVal > 0.5
      # apply integer cuts
      Z = objVal
      while (objVal <= Z) && (objVal > 0.5)
        # println("=============current Z value $Z")
        m, objVal, status, solveTime, yresult = fix_metabolite(
            isMd, isDb, isRev, isCyt, isExt, i,
            stoiMatrix, fluxLB, fluxUB,
            solver, epsilon, bigM, nonZero, Z, results[i])
        println(objVal)
        println(yresult)
        # store valid result only
        if (objVal <= Z) && (objVal > 0.5)
          push!(results[i], yresult)
        end
      end
    end
  end
  return results
end

function fix_metabolite(isMd, isDb, isRev, isCyt, isExt, ID,
    stoiMatrix, fluxLB, fluxUB, solver, epsilon, bigM, nonZero, intCuts...)
  t1 = time()
  println("---tackling no-production metabolite \"$ID\"..")
  noCompound, noRxn = size(stoiMatrix)
  m = Model(with_optimizer(solver.Optimizer))
  # define variables using @variables blocks
  @variable(m, v[1:noRxn])  # fluxes
  # # 1 if the reversibility of rnx i in model is relaxed
  # @variable(m, x[i=1:noRxn; (isMd[i] && !isRev[i])], Bin)
  # 1 if rxn i is relaxed or added
  @variable(m, y[i=1:noRxn; (isDb[i] || (isMd[i] && !isRev[i]))], Bin)
  # 1 if metabolite ID is produced in reaction j
  @variable(m, w[1:noRxn], Bin)
  # @variable(m, w[1:noCompound, 1:noRxn], Bin)

  # define objective
  @objective(m, Min, sum(y))
  # define constraints
  # Mass Balance (Eq. 6.9 & 6.10 on P124)
  for i = 1:noCompound
    if (isCyt[i] && !isExt[i])
      # cytosol but not extracellular Eq. 6.9
      @constraint(m, sum(stoiMatrix[i,j]*v[j] for j=1:noRxn) >= 0)
    elseif (!isCyt[i] && !isExt[i])
      # other compartments but not extracellular, Eq. 6.10
      @constraint(m, sum(stoiMatrix[i,j]*v[j] for j=1:noRxn) == 0)
    end
  end
  # Production Constrains for (ir)reversible reactions (Eq. 6.11 on P124)
  for j in 1:noRxn
    # keep this pre-condition to reduce # of constraints
    if (abs(stoiMatrix[ID,j]) > nonZero)  # non-zero coefficient
      @constraint(m, stoiMatrix[ID,j]*v[j] >= epsilon - bigM*(1-w[j]))
      @constraint(m, stoiMatrix[ID,j]*v[j] <= bigM*w[j])
    end
  end
  # Binary Constraints (Eq. 6.12 on P124)
  @constraint(m, sum(w[j] for j=1:noRxn if
    (abs(stoiMatrix[ID,j]) > nonZero)) >= 1)
  # Bound constraints (Eq. 6.13-15 on P124)
  for j in 1:noRxn
    if isMd[j] # Eq. 6.13-14
      @constraint(m, v[j] <= fluxUB[j])
      if isRev[j]  # Eq. 6.13
        @constraint(m, v[j] >= fluxLB[j])
      else # Eq. 6.14
        @constraint(m, v[j] >= -bigM*y[j])
      end
    else  # Eq. 6.15
      @constraint(m, v[j] >= fluxLB[j]*y[j])
      @constraint(m, v[j] <= fluxUB[j]*y[j])
    end
  end
  # add integer cuts
  if length(intCuts) != 0
    for cut in intCuts[2]
      @constraint(m, sum(y[j] for j in cut) <= intCuts[1]-1)
    end
  end
  # Calling solver
  JuMP.optimize!(m)
  t2 = time()
  # get some results
  status = termination_status(m)
  solveTime = t2 - t1
  yresult = Set{Int64}()
  if (status == MOI.OPTIMAL) || (status == MOI.TIME_LIMIT && has_values(m))
    objVal = objective_value(m)
    # extract y
    for i in 1:noRxn
      if isDb[i] || (isMd[i] && !isRev[i])
        if value(y[i]) == 1
          push!(yresult, i)
        end
      end
    end
  else
    println("The model was not solved correctly.")
    objVal = -1
  end
  return m, objVal, status, solveTime, yresult
end


text:
struct UNIQUACParam <: EoSParam
    a::PairParam{Float64}
    r::SingleParam{Float64}
    q::SingleParam{Float64}
    q_p::SingleParam{Float64}
    Mw::SingleParam{Float64}
end

abstract type UNIQUACModel <: ActivityModel end

struct UNIQUAC{c<:EoSModel} <: UNIQUACModel
    components::Array{String,1}
    icomponents::UnitRange{Int}
    params::UNIQUACParam
    puremodel::Vector{c}
    absolutetolerance::Float64
    references::Array{String,1}
end
@registermodel UNIQUAC
export UNIQUAC

function UNIQUAC(components::Vector{String}; puremodel=PR,
    userlocations=String[],
     verbose=false)
    params = getparams(components, ["Activity/UNIQUAC/UNIQUAC_like.csv", "properties/molarmass.csv","Activity/UNIQUAC/UNIQUAC_unlike.csv"]; userlocations=userlocations, asymmetricparams=["a"], ignore_missing_singleparams=["a"], verbose=verbose)
    a  = params["a"]
    r  = params["r"]
    q  = params["q"]
    q_p = params["q_p"]
    Mw  = params["Mw"]
    icomponents = 1:length(components)

    init_puremodel = [puremodel([components[i]]) for i in icomponents]
    packagedparams = UNIQUACParam(a,r,q,q_p,Mw)
    references = String[]
    model = UNIQUAC(components,icomponents,packagedparams,init_puremodel,1e-12,references)
    return model
end

function activity_coefficient(model::UNIQUACModel,p,T,z)
    a = model.params.a.values
    q = model.params.q.values
    q_p = model.params.q_p.values
    r = model.params.r.values

    x = z ./ sum(z)

    Φ =  x.*r/sum(x[i]*r[i] for i ∈ @comps)
    θ = x.*q/sum(x[i]*q[i] for i ∈ @comps)
    θ_p = x.*q_p/sum(x[i]*q_p[i] for i ∈ @comps)
    τ = @. exp(-a/T)
    lnγ_comb = @. log(Φ/x)+(1-Φ/x)-5*q*(log(Φ/θ)+(1-Φ/θ))
    lnγ_res  = q_p.*(1 .-log.(sum(θ_p[i]*τ[i,:] for i ∈ @comps)) .-sum(θ_p[i]*τ[:,i]/sum(θ_p[j]*τ[j,i] for j ∈ @comps) for i ∈ @comps))
    return exp.(lnγ_comb+lnγ_res)
end


text:
# This file contains specification and implementation
# Of the homoiconic DSL for defining rules and theories
using MatchCore

function rule_sym_map(ex)
    @smatch ex begin
        :($a => $b) => RewriteRule
        :($a |> $b) => DynamicRule
        :($a == $b) => EqualityRule
        :($a != $b) => UnequalRule
        :($a ≠ $b) => UnequalRule
        _ => error("Cannot parse rule from $ex")
    end
end

interp_dollar(x, mod::Module) = x
function interp_dollar(ex::Expr, mod::Module)
    if Meta.isexpr(ex, :$)
        mod.eval(ex.args[1])
    else
        Expr(ex.head, map(x -> interp_dollar(x, mod), ex.args)...)
    end
end

"""
Construct an `AbstractRule` from a quoted expression.
You can also use the [`@rule`] macro to
create a `Rule`.
"""
function Rule(e::Expr, mod::Module=@__MODULE__, resolve_fun=false)
    op = gethead(e)
    RuleType = rule_sym_map(e)
    l, r = e.args[Meta.isexpr(e, :call) ? (2:3) : (1:2)]

    l = interp_dollar(l, mod)

    if RuleType !== DynamicRule
        r = interp_dollar(r, mod)
    end

    lhs = Pattern(l, mod, resolve_fun)
    rhs = r

    if RuleType <: SymbolicRule
        rhs = Pattern(rhs, mod, resolve_fun)
    end

    return RuleType(lhs, rhs)
end

# fallback when defining theories and there's already a rule
function Rule(r::AbstractRule, mod::Module=@__MODULE__, resolve_fun=false)
    r
end

macro rule(e)
    e = macroexpand(__module__, e)
    e = rmlines(copy(e))
    Rule(e, __module__, false)
end

macro methodrule(e)
    e = macroexpand(__module__, e)
    e = rmlines(copy(e))
    Rule(e, __module__, true)
end

# Theories can just be vectors of rules!

macro theory(e)
    e = macroexpand(__module__, e)
    e = rmlines(e)
    # e = interp_dollar(e, __module__)
    if Meta.isexpr(e, :block)
        Vector{AbstractRule}(e.args .|> x -> Rule(x, __module__, false))
    else
        error("theory is not in form begin a => b; ... end")
    end
end

# TODO document this puts the function as pattern head instead of symbols
macro methodtheory(e)
    e = macroexpand(__module__, e)
    e = rmlines(e)
    # e = interp_dollar(e, __module__)
    if Meta.isexpr(e, :block)
        Vector{AbstractRule}(e.args .|> x -> Rule(x, __module__, true))
    else
        error("theory is not in form begin a => b; ... end")
    end
end

"""
A Theory is either a vector of [`Rule`](@ref) or
a compiled, callable function.
"""
const Theory = Union{Vector{<:AbstractRule}, Function}


text: tests_are_sorted = false


text:
function jaccard_index_basic(a::Vector{Int},
                             b::Vector{Int},
                             v::Int)
  agree = 0.0
  disagree = 0.0

  for j in 1:(v - 1), i in (j + 1):v
    if (a[i] == a[j]) && (b[i] == b[j])
      agree += 1
    elseif ((a[i] == a[j]) && (b[i] != b[j])) ||
           ((a[i] != a[j]) && (b[i] == b[j]))
      disagree += 1
    end
  end

  (agree > 0 || disagree > 0) ? agree / (agree + disagree) : 1.0
end

function jaccard_index_table(a::Vector{Int},
                             ka::Int,
                             b::Vector{Int},
                             kb::Int,
                             v::Int)
  agree = 0.0
  disagree = 0.0

  count = zeros(Float64, ka, kb)
  for i in 1:v
    count[a[i], b[i]] += 1.0
  end

  za = zeros(Float64, ka)
  zb = 0.0
  zab = 0.0

  for j in 1:kb
    zb = 0.0

    for i in 1:ka
      za[i] += count[i, j]
      zb += count[i, j]
      zab += count[i, j]^2

      agree += count[i, j] * (count[i, j] - 1) / 2
    end

    disagree += zb^2
  end

  for i in 1:ka
    disagree += za[i]^2
  end

  disagree /= 2
  disagree -= zab

  (agree > 0 || disagree > 0) ? agree / (agree + disagree) : 1.0
end

function jaccard(a::Vector{Int},
                 ka::Int,
                 b::Vector{Int},
                 kb::Int,
                 v::Int)
  # total number of for iterations for jaccard_index_table is v + ka + ka * kb
  # v iterations to compute counts
  # ka * kb iterations to compute various quantities
  # ka final iterations to compute row sums

  # total number of for iterations for jaccard_index_basic is v * (v - 1) / 2
  # which is equivalent to the total number of distinct pairs
  index = if (v + ka * (1 + kb)) < (v * (v - 1) / 2)
    jaccard_index_table(a, ka, b, kb, v)
  else
    jaccard_index_basic(a, b, v)
  end

  1.0 - index
end

function jaccardlower(n::Vector{T},
                      k::Int,
                      v::Real) where T <: Real
  agree = 0.0
  disagree = 0.0

  zab = 0.0
  for g in 1:k
    zab += n[g]^2

    agree += n[g] * (n[g] - 1) / 2
    disagree += n[g]^2
  end

  disagree = (disagree + v^2) / 2 - zab

  (agree > 0 || disagree > 0) ? 1.0 - agree / (agree + disagree) : 0.0
end


text:
"Parse command line arguments"
function parse_commandline()
    setting = ArgParseSettings()

    @add_arg_table! setting begin
        "--size", "-s"
            dest_name="size"
            help="size of the Cayley tree"
            arg_type=Int
            default=5
        "--num", "-n"
            dest_name="num"
            help="number of trees to generate"
            arg_type=Int
            default=1
        "--log-size", "-l"
            dest_name="log_size"
            help="log base 10 of the size of the tree"
            arg_type=Int
        "--height", "-g"
            dest_name="tall"
            help="height of the tree (only takes effect for full d-ary trees)"
            arg_type=Int
            default=0
        "--tree", "-t"
            dest_name="tree_type"
            help="type of tree: " * join(keys(TREE_DICT), ", ")
            arg_type=String
            default="Cayley"
        "--seed", "-e"
            dest_name="seed"
            help="seed for random number generators"
            arg_type=Int
        "-d"
            dest_name="d"
            help="d-ary tree"
            arg_type=Int
        "kcut"
            action = :command
            help="simulation of the k-cut number"
        "log-prod"
            action = :command
            help="simulation of log product of subtree sizes"
        "size-power"
            action = :command
            help="simulation of sum of power of subtree sizes"
        "leaf"
            action = :command
            help="simulation of the number of leaves coutning"
        "height"
            action = :command
            help="simulation of the height"
        "total-path"
            action = :command
            help="simulation of the total path length"
    end

    @add_arg_table! setting["kcut"] begin
        "-k"
            dest_name="k"
            help="k-records"
            arg_type=Int
            default=2
        "--segment", "-g"
            dest_name="segment"
            help="number of segments to track"
            arg_type=Int
            default=0
    end

    @add_arg_table! setting["log-prod"] begin
        "--pow", "-p"
            dest_name="pow"
            help="log(subtree size)^pow"
            arg_type=Int
            default=1
    end

    @add_arg_table! setting["size-power"] begin
        "--pow", "-p"
            dest_name="pow"
            help="subtree size^pow"
            arg_type=Real
            default=0.5
    end

    return parse_args(setting)
end


text:
using HSARuntime
using HSARuntime.HSA
using Test

using Pkg
Pkg.add(PackageSpec(;name="GPUCompiler",rev="master"))
Pkg.add(PackageSpec(;name="AMDGPUnative",rev="master"))
using AMDGPUnative

if HSARuntime.configured
    @testset "HSA Status Error" begin
        errorcode = HSARuntime.HSAError(HSA.STATUS_SUCCESS)
        @test HSARuntime.description(errorcode) == "HSA_STATUS_SUCCESS: The function has been executed successfully."
    end
    agents = get_agents()
    if length(agents) > 0
        include("agent.jl")
        include("array.jl")
        include("memory.jl")

        if AMDGPUnative.configured
            include("vadd.jl")
            include("global.jl")
        else
            @warn "AMDGPUnative.jl has not been configured; skipping on-device tests"
        end
    else
        @warn "No devices detected; skipping runtime tests"
    end
else
    error("HSARuntime.jl has not been configured; skipping all tests")
end


text:
"""
   solving ||Ax-b|| via CGLS, This code is a translated version of SOL's matlab code
   the resNE (normalized relative error is ||A'b-A'Ax||/||A'b|| )
"""
function cgls(A, b::Vector{Tv}; params=Dict(),
              x0=0.0, shift=0.0, tol=1e-6, maxIter=30, print_flag=false) where {Tv<:Number}

    # determine the type of A
    if isa(A, Matrix)
       explicitA = true
    elseif isa(A, Function)
       explicitA = false
    else
       error("A must be a matrix or a function")
    end

    # initial guess of the unknowns
    if explicitA

       (m, n) = size(A)

       # provided initial guess
       if x0 != 0.0
          x = copy(x0)
          r = b - A * x
          s = A' * r - shift * x
       # zero vector as initial guess
       else
          x = zeros(eltype(b), n)
          r = copy(b)
          s = A' * r
       end

    else

       m = length(b)

       # provided initial guess
       if x0 != 0.0
          x = copy(x0)
          r = b - A(x, 1; params...)
          s = A(r, 2; params...) - shift * x
          n = length(s)
       # zero vector as initial guess
       else
          r = copy(b)
          s = A(b, 2; params...)
          n = length(s)
          x = zeros(eltype(b), n)
       end

    end

    data_fitting = dot(r,r)
    constraint   = 0.0
    cost0 = data_fitting + constraint

    # initialize some intermidiate vectors
    p = copy(s)
    norms0= vecnorm(s)
    gamma = norms0^2
    normx = vecnorm(x)
    xmax  = normx     # keep the initial one
    resNE = 1.0

    gamma0= copy(gamma)
    delta = 0.0

    # iteration counter and stop condition
    k = 0
    run_flag = true

    if print_flag
       header = "  k         data_fitting           constraint             normx                resNE"
       println(""); println(header);
       @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
    end

    while (k < maxIter) && run_flag

          k = k + 1

          if explicitA
             q = A * p
          else
             q = A(p, 1; params...)
          end

          delta = (vecnorm(q))^2 + shift * (vecnorm(p))^2
          indefinite = delta <= 0 ? true : false
          delta      = delta == 0.? eps(): delta

          alpha = gamma / delta

          x = x + alpha * p
          r = r - alpha * q

          data_fitting = dot(r,r)
          constraint   = shift * dot(x, x)
          cost = data_fitting + constraint

          if explicitA
             s = A' * r - shift * x
          else
             s = A(r, 2; params...) - shift * x
          end

          norms  = norm(s)
          gamma0 = copy(gamma)
          gamma  = norms^2
          beta   = gamma / gamma0

          p = s + beta * p

          # check the stopping crietia
          normx = norm(x)
          xmax  = normx > xmax ? normx : xmax
          if norms <= norms0 * tol || normx * tol >= 1.0
             run_flag = false
          end

          # print information
          # resNE = norms / norms0
          resNE = cost  / cost0
          if print_flag
             @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
          end

    end

    return x, resNE

end

# ==============================================================================
#                 CGLS based on hard drive
# ==============================================================================
"""
   solving ||Ax-b|| via CGLS
"""
function cgls(A::Function, b::Ts; x0="NULL", path="NULL", params=[],
         shift=0.0, tol=1e-6, maxIter=50, print_flag=false) where {Ts<:Union{String, Vector{String}}}

    if path == "NULL"
       path = pwd()
    end

    # create folder to save files of model size
    path_model = joinpath(path, "model")
    if !isdir(path_model)
       mkdir(path_model)
       if !isdir(path_model)
          error("could not create directory")
       end
    end
    x = join([path_model "/x.bin"])
    s = join([path_model "/s.bin"])
    p = join([path_model "/p.bin"])

    # create folders to save residue
    path_residue   = joinpath(path, "residue")
    if !isdir(path_residue)
       mkdir(path_residue)
       if !isdir(path_residue)
          error("could not create directory")
       end
    end

    # create folders to save synthetic data with current model
    path_forward = joinpath(path, "forward")
    if !isdir(path_forward)
       mkdir(path_forward)
       if !isdir(path_forward)
          error("could not create directory")
       end
    end

    # intermidiate variables
    if typeof(b) == String
       r = join([path_residue "/res.bin"])
       q = join([path_forward "/fwd.bin"])
    elseif typeof(b) == Vector{String}
       n = length(b)
       r = Vector{String}(n)
       q = Vector{String}(n)
       for i = 1 : n
           r[i] = join([path_residue "/res" "_" "$i" ".bin"])
           q[i] = join([path_forward "/fwd" "_" "$i" ".bin"])
       end
    end

    # provided initial guess
    if x0 != "NULL"

       cp(x0, x, remove_destination=true)
       dcal = A(x, 1; params=params)
       x_plus_alpha_q!(r, b, -1.0, dcal)

       if typeof(params) <: Dict{Symbol, Any}
          params[:path_m] = s
       elseif typeof(params) == Vector{Dict}
         for i = 1 : length(params)
             params[i][:path_m] = s
         end
       end
       s = A(r, 2; params=params)
       x_plus_alpha_q!(s, -shift, x)

    # zero vector as initial guess
    else

       cp(b, r, remove_destination=true)
       if typeof(params) <: Dict{Symbol, Any}
          params[:path_m] = s
       elseif typeof(params) == Vector{Dict}
          for i = 1 : length(params)
              params[i][:path_m] = s
          end
       end
       s = A(r, 2; params=params)

       # create x
       (hdr, ds) = read_USdata(s)
       tmp = zeros(ds)
       write_USdata(x, hdr, tmp)
    end

    # compute residue
    data_fitting = (norm(r))^2
    constraint   = 0.0
    cost0 = data_fitting + constraint
    convergence = Float64[]; push!(convergence, 1.0);

    # initialize some intermidiate vectors
    cp(s, p, remove_destination=true)

    norms0= norm(s)
    gamma = norms0^2
    normx = norm(x)
    xmax  = normx     # keep the initial one
    resNE = 1.0

    gamma0= copy(gamma)
    delta = 0.0

    # iteration counter and stop condition
    k = 0
    run_flag = true

    if print_flag
       header = "  k         data_fitting           constraint             normx                resNE"
       println(""); println(header);
       @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
    end

    if typeof(params) <: Dict{Symbol, Any}
       params[:path_fwd] = q
    elseif typeof(params) == Vector{Dict}
       for i = 1 : length(params)
           params[i][:path_fwd] = q[i]
       end
    end

    while (k < maxIter) && run_flag

          k = k + 1

          q = A(p, 1; params=params)

          delta = (norm(q))^2 + shift * (norm(p))^2
          indefinite = delta <= 0 ? true : false
          delta      = delta == 0.? eps(): delta

          alpha = gamma / delta

          x_plus_alpha_q!(x,  alpha, p)
          x_plus_alpha_q!(r, -alpha, q)

          data_fitting = (norm(r))^2
          constraint   = shift * (norm(x))^2
          cost = data_fitting + constraint

          # save the intermidiate result
          path_iter = join([path_model "/iteration" "_" "$k" ".bin"])
          cp(x, path_iter, remove_destination=true)

          s = A(r, 2; params=params)
          x_plus_alpha_q!(s, -shift, x)

          norms  = norm(s)
          gamma0 = gamma
          gamma  = norms^2
          beta   = gamma / gamma0

          (hdr, ds) = read_USdata(s)
          (hdr, dp) = read_USdata(p)
          dp = ds + beta * dp
          write_USdata(p, hdr, dp)

          # check the stopping crietia
          normx = norm(x)
          xmax  = normx > xmax ? normx : xmax
          if norms <= norms0 * tol || normx * tol >= 1.0
             run_flag = false
          end

          # print information
          resNE = cost / cost0
          if print_flag
             @printf("%3.0f %20.10e %20.10e %20.10e %20.10e\n", k, data_fitting, constraint, normx, resNE);
          end
          push!(convergence, resNE);
    end

    return x, convergence

end


text:
using Base.Test, Convex, SCS, SchattenNorms

set_default_solver(SCSSolver(verbose=0))

x,y,z = [0 1; 1 0], [0 -1im; 1im 0], [1 0; 0 -1]

liou(x) = kron(conj(x),x)
liou(x,y) = kron(transpose(y),x)

function rand_unitary(d)
    rm = randn(d,d)+1im*randn(d,d)
    u,_,_ = svd(rm)
    return u
end

cnot = [1. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 0. 1.; 0. 0. 1. 0.]
cnotL = liou(cnot)

#Check the diamond norm of randomly chosen Pauli channels
p1 = rand(4); p1 = p1/sum(p1)
p2 = rand(4); p2 = p2/sum(p2)
pauli_dnorm = sum(abs(p1-p2))
pauli = Matrix{Complex128}[ [1 0; 0 1],
                            [0 1; 1 0],
                            [0 -im; im 0],
                            [1 0; 0 -1] ]

u = liou(rand_unitary(2))
v = liou(rand_unitary(2))

#ee = rand_cp_map(2)
#ff = rand_cp_map(2)

@test_approx_eq_eps 2.0 dnorm(cnotL-eye(16)) 1e-5
@test_approx_eq_eps 2.0 dnorm(eye(4)-liou(x,x)) 1e-5
@test_approx_eq_eps 2.0 dnorm(eye(4)-liou(z,z)) 1e-5
@test_approx_eq_eps 2.0 dnorm(eye(4)-liou(y,y)) 1e-5
@test_approx_eq_eps 2.0 dnorm(liou(x,x)-liou(y,y)) 1e-5
@test_approx_eq_eps 2.0 dnorm(liou(x,x)-liou(z,z)) 1e-5
@test_approx_eq_eps 2.0 dnorm(liou(z,z)-liou(y,y)) 1e-5

@test_approx_eq_eps 2.0 dnorm(u*eye(4)-u*liou(x,x)) 1e-5
@test_approx_eq_eps 2.0 dnorm(u*eye(4)-u*liou(z,z)) 1e-5
@test_approx_eq_eps 2.0 dnorm(u*eye(4)-u*liou(y,y)) 1e-5
@test_approx_eq_eps 2.0 dnorm(u*liou(x,x)-u*liou(y,y)) 1e-5
@test_approx_eq_eps 2.0 dnorm(u*liou(x,x)-u*liou(z,z)) 1e-5
@test_approx_eq_eps 2.0 dnorm(u*liou(z,z)-u*liou(y,y)) 1e-5

pc1 = sum([ p1[ii+1]*liou(pauli[ii+1]) for ii in 0:3])
pc2 = sum([ p2[ii+1]*liou(pauli[ii+1]) for ii in 0:3])
calc_dnorm1 = dnorm(pc1-pc2)

@test_approx_eq_eps pauli_dnorm calc_dnorm1  1e-5

duv  = dnorm(u-v)
@time begin
    duv2 = dnorm(u*v'-eye(4))
    duv3 = dnorm(v'u-eye(4))
    duv4 = dnorm(eye(4)-u'*v)
    duv5 = dnorm(eye(4)-v*u')
end

@test_approx_eq_eps duv duv2 1e-5
@test_approx_eq_eps duv duv3 1e-5
@test_approx_eq_eps duv duv4 1e-5
@test_approx_eq_eps duv duv5 1e-5

#def = dnorm(ee-ff)
#@time begin
#    def2 = dnorm(u*(ee-ff)*v)
#    def3 = dnorm(u*(ee-ff))
#    def4 = dnorm((ee-ff)*v)
#    def5 = dnorm(u'*(ee-ff))
#end

#println(def)
#println(def2)
#println(def3)
#println(def4)
#println(def5)


text:
"""
    convert_to_FlowFields(U::Array{T,2},V::Array{T,2},t1::T) where T

Convert a pair of U,V arrays (staggered C-grid velocity field in 2D) to
a `𝐹_MeshArray2D` struct ready for integration of individual displacements
from time `t0=0` to time `t1`.

```
_,u,v=random_flow_field()
𝐹=convert_to_FlowFields(u,v,10.0)
```
"""
function convert_to_FlowFields(U::Array{T,2},V::Array{T,2},t1::T) where T
    np,nq=size(U)
    Γ=simple_periodic_domain(np,nq)

    g=Γ.XC.grid
    u=MeshArray(g,[U])
    v=MeshArray(g,[V])
    (u,v)=exchange(u,v,1)
    func=(u -> IndividualDisplacements.update_location_dpdo!(u,g))

    𝐹_MeshArray2D{eltype(u)}(u,u,v,v,[0,t1],func)
end

"""
    postprocess_MeshArray(sol,𝑃::FlowFields; id=missing, 𝑇=missing)

Copy `sol` to a `DataFrame` & map position to lon,lat coordinates
using "exchanged" 𝐷.XC, 𝐷.YC via `add_lonlat!`
"""
function postprocess_MeshArray(sol::ODESolution,𝑃::FlowFields; id=missing, 𝑇=missing)
    ismissing(id) ? id=collect(1:size(sol,2)) : nothing
    ismissing(𝑇) ? 𝑇=𝑃.𝑇 : nothing

    nd=length(size(sol))
    nt=size(sol,nd)
    nf=size(sol,nd-1)
    id=id*ones(1,size(sol,nd))
    if (size(sol,1)>1)&&(nd>2)
        x=sol[1,:,:]
        y=sol[2,:,:]
        fIndex=sol[end,:,:]
    elseif (nd>2)
        x=[sol[1,i,j][1] for i in 1:nf, j in 1:nt]
        y=[sol[1,i,j][2] for i in 1:nf, j in 1:nt]
        fIndex=[sol[1,i,j][end] for i in 1:nf, j in 1:nt]
    else
        x=sol[1,:]
        y=sol[2,:]
        fIndex=sol[end,:]
        nf=1
    end

    𝑃.u0.grid.nFaces==1 ? fIndex=ones(size(x)) : nothing

    t=[ceil(i/nf)-1 for i in 1:nt*nf]
    t=𝑇[1] .+ (𝑇[2]-𝑇[1])/t[end].*t

    df = DataFrame(ID=Int.(id[:]), x=x[:], y=y[:], fid=Int.(fIndex[:]), t=t[:])
    return df
end

"""
    add_lonlat!(df::DataFrame,XC,YC)

Add lon & lat to dataframe using "exchanged" XC, YC
"""
function add_lonlat!(df::DataFrame,XC,YC)
    x = cosd.(YC) * cosd.(XC)
    y = cosd.(YC) * sind.(XC)
    z = sind.(YC)

    df.lon=0*df.x
    df.lat=0*df.x
    for i in eachindex(df.lon)
        xx=interp_to_xy(df.x[i],df.y[i],df.fid[i],x)
        yy=interp_to_xy(df.x[i],df.y[i],df.fid[i],y)
        zz=interp_to_xy(df.x[i],df.y[i],df.fid[i],z)
        df.lat[i] = asind(zz/sqrt(xx^2+yy^2+zz^2))
        df.lon[i] = atand(yy, xx)
    end

    return df
end

"""
    add_lonlat!(df::DataFrame,XC,YC,func::Function)

Add lon & lat to dataframe using "exchanged" XC, YC after updating
subdomain indices (via func) if needed (location_is_out)
"""
function add_lonlat!(df::DataFrame,XC,YC,func::Function)
    g=XC.grid
    u=zeros(3)

    for i in eachindex(df.lon)
        u[:]=[df.x[i];df.y[i];df.fid[i]]
        while location_is_out(u,g)
            func(u)
            df.x[i]=u[1]
            df.y[i]=u[2]
            df.fid[i]=u[3]
        end
    end

    add_lonlat!(df,XC,YC)

    return df
end

"""
    postprocess_xy()

Copy `sol` to a `DataFrame` & map position to x,y coordinates,
and define time axis for a simple doubly periodic domain
"""
function postprocess_xy(sol,𝑃::FlowFields; id=missing, 𝑇=missing)
    ismissing(id) ? id=collect(1:size(sol,2)) : nothing
    ismissing(𝑇) ? 𝑇=𝑃.𝑇 : nothing

    nf=size(sol,2)
    nt=size(sol,3)

    isa(𝑃.u0,MeshArray) ? (nx,ny)=𝑃.u0.grid.ioSize[1:2] : (nx,ny)=size(𝑃.u0)[1:2]
    nd=length(size(sol))

    id=id*ones(1,size(sol,nd))
    if (size(sol,1)>1)&&(nd>2)
        x=mod.(sol[1,:,:],Ref(nx))
        y=mod.(sol[2,:,:],Ref(ny))
    elseif (nd>2)
        x=[mod(sol[1,i,j][1],nx) for i in 1:nf, j in 1:nt]
        y=[mod(sol[1,i,j][2],ny) for i in 1:nf, j in 1:nt]
    else
        x=mod.(sol[1,:],Ref(nx))
        y=mod.(sol[2,:],Ref(ny))
    end
    t=[ceil(i/nf)-1 for i in 1:nt*nf]
    #size(𝐷.XC,1)>1 ? fIndex=sol[3,:,:] : fIndex=fill(1.0,size(x))
    t=𝑇[1] .+ (𝑇[2]-𝑇[1])/t[end].*t

    return DataFrame(ID=Int.(id[:]), t=t[:], x=x[:], y=y[:])
end

"""
    randn_lonlat(nn=1,seed=missing)

Randomly distributed longitude, latitude positions on the sphere.
"""
function randn_lonlat(nn=1;seed=missing)
    !ismissing(seed) ? rng = MersenneTwister(1234) : rng = MersenneTwister()
    tmp = randn(rng, Float64, (nn, 3))
    tmpn = tmp ./ sqrt.(sum(tmp.^2, dims=2))
    lon = rad2deg.(atan.(tmpn[:,2], tmpn[:,1]))
    lat = 90.0 .- rad2deg.(acos.(tmpn[:,3]))
    return lon, lat
end

"""
    nearest_to_xy(α::MeshArray,x,y,f)

Value of α at eachindex of the grid cell center nearest to `x,y` on subdomain array / facet `f`
"""
nearest_to_xy(α::MeshArray,x,y,f) = [α[Int(f[i]),1][ Int(round(x[i] .+ 0.5)), Int(round(y[i] .+ 0.5)) ] for i in eachindex(x)]

"""
    nearest_to_xy(α::Array,x,y)

Value of α at eachindex of the grid cell center nearest to `x,y`
"""
nearest_to_xy(α::Array,x,y) = [α[ Int(round(x[i] .+ 0.5)), Int(round(y[i] .+ 0.5)) ] for i in eachindex(x)]

"""
    interp_to_lonlat

Use MeshArrays.Interpolate() to interpolate to e.g. a regular grid (e.g. maps for plotting purposes).

```jldoctest
using IndividualDisplacements
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/jupyter/helper_functions.jl"))
𝑃,𝐷=global_ocean_circulation(k=1,ny=2);

lon=[i for i=20.:20.0:380., j=-70.:10.0:70.]
lat=[j for i=20.:20.0:380., j=-70.:10.0:70.]
(f,i,j,w,_,_,_)=InterpolationFactors(𝐷.Γ,vec(lon),vec(lat))
IntFac=(lon=lon,lat=lat,f=f,i=i,j=j,w=w)

tmp1=interp_to_lonlat(𝐷.Γ.Depth,𝐷.Γ,lon,lat)
tmp2=interp_to_lonlat(𝐷.Γ.Depth,IntFac)

ref=[5896. 5896.]
prod(isapprox.([maximum(tmp1) maximum(tmp2)],ref,atol=1.0))

# output

true
```
"""
function interp_to_lonlat(X::MeshArray,Γ::NamedTuple,lon,lat)
    (f,i,j,w,_,_,_)=InterpolationFactors(Γ,vec(lon),vec(lat))
    return reshape(Interpolate(X,f,i,j,w),size(lon))
end

function interp_to_lonlat(X::MeshArray,IntFac::NamedTuple)
    @unpack f,i,j,w,lon,lat = IntFac
    return reshape(Interpolate(X,f,i,j,w),size(lon))
end

"""
    interp_to_xy(df::DataFrame,Zin::MeshArray)

Interpolate "exchanged" / "hallo-included" Zin to df[!,:x], df[!,:y] on df[!,:fid]
"""
function interp_to_xy(df::DataFrame,Zin::MeshArray)
    x=df[!,:x];
    y=df[!,:y];
    f=Int.(df[!,:fid]);
    dx,dy=(x - floor.(x) .+ 0.5,y - floor.(y) .+ 0.5);
    i_c = Int32.(floor.(x)) .+ 1;
    j_c = Int32.(floor.(y)) .+ 1;

    Z=zeros(length(x),4)
    [Z[k,:]=Zin[f[k]][i_c[k]:i_c[k]+1,j_c[k]:j_c[k]+1][:] for k in 1:length(i_c)]

    return (1.0 .-dx).*(1.0 .-dy).*Z[:,1]+dx.*(1.0 .-dy).*Z[:,2] +
           (1.0 .-dx).*dy.*Z[:,3]+dx.*dy.*Z[:,4]
end

function interp_to_xy(x,y,f,Zin::MeshArray)
    dx,dy=(x - floor(x) + 0.5,y - floor(y) + 0.5)
    i_c = Int(floor(x)) + 1
    j_c = Int(floor(y)) + 1
    ff=Int(f)

    return (1.0-dx)*(1.0-dy)*Zin[ff][i_c,j_c] +
    dx*(1.0-dy)*Zin[ff][i_c+1,j_c] +
    (1.0-dx)*dy*Zin[ff][i_c,j_c+1] +
    dx*dy*Zin[ff][i_c+1,j_c+1]
end

"""
    stproj(XC,YC;XC0=0.0,YC0=90.0)

Apply to XC,YC (longitude, latitude) the stereographic projection
which would put XC0,YC0 (longitude, latitude) at x,y=0,0
"""
function stproj(XC,YC,XC0=0.0,YC0=90.0)

# compute spherical coordinates:
phi=pi/180*XC; theta=pi/180*(90-YC);
phi0=pi/180*XC0; theta0=pi/180*(90-YC0);

# compute cartesian coordinates:
x=sin(theta)*cos(phi)
y=sin(theta)*sin(phi)
z=cos(theta)

# bring chosen longitude to x>0,y=0:
xx=cos(phi0)*x+sin(phi0)*y
yy=-sin(phi0)*x+cos(phi0)*y
zz=z
# bring chosen point to South Pole:
x=cos(theta0)*xx-sin(theta0)*zz
y=yy
z=sin(theta0)*xx+cos(theta0)*zz

# stereographic projection from South Pole:
xx=x/(1+z)
yy=y/(1+z)

return xx,yy
end

"""
    stproj_inv(xx,yy;XC0=0.0,YC0=90.0)

Apply to xx,yy (stereographic projection coordinates) the reverse
of the stereographic projection which puts XC0,YC0 (longitude,
latitude) at x,y=0,0
"""
function stproj_inv(xx,yy,XC0=0.0,YC0=90.0)
    phi0=pi/180*XC0; theta0=pi/180*(90-YC0);

    # Reverse stereographic projection from North Pole:
    (x,y,z)=2*[xx yy (1-xx^2-yy^2)/2]./(1+xx^2+yy^2)

    # bring chosen point back from South Pole:
    xx=cos(theta0)*x+sin(theta0)*z
    yy=y
    zz=-sin(theta0)*x+cos(theta0)*z

    # bring chosen longitude back from x>0,y=0:
    x=cos(phi0)*xx-sin(phi0)*yy
    y=sin(phi0)*xx+cos(phi0)*yy
    z=zz

    # compute spherical coordinates:
    theta=atan(sqrt(x^2+y^2)/z)
    phi=atan(y, x)

    XC=180/pi*phi
    theta>=0 ? YC=90-180/pi*theta : YC=-90-180/pi*theta

    return XC,YC
end


text:
module test_test_broken

using Test

@test_broken 1 == 2
@test_skip   1 == 2

end # module test_test_broken


text:
@with_kw struct ContinuousWave <: AbstractWaveform
    fc::Float64
    t0::Float64 = 6 / fc
end

function (s::ContinuousWave)(t)
    g(t) = sin(π / 2 * (t / s.t0))
    0 ≤ t ≤ s.t0 ? g(t) * sin(2π * s.fc * (t - s.t0)) : sin(2π * s.fc * (t - s.t0))
end

@recipe function f(s::ContinuousWave; nsamples = 30, spectrum = true)
    linewidth := 2
    label := "Sinusoidal"
    δt = 1 / s.fc / nsamples
    t_samples = range(0, s.t0 + 10 * 1 / s.fc, step = δt)
    if spectrum
        return fourier(t_samples, s.(t_samples))
    else
        return (t_samples, s.(t_samples))
    end
end


text:
# total_cpi_rebase.jl - Función de inflación total con cambio de base sintético.

# El parámetro de maxchanges permite realizar hasta un máximo de cambios de base
# sintéticos en cada VarCPIBase. El valor cero indica que se utilizarán todos
# los cambios de base sintéticos posibles dentro de cada base
Base.@kwdef struct InflationTotalRebaseCPI <: InflationFunction
    period::Int = 60 # períodos para realizar el cambio de base sintético
    maxchanges::Int = 0
end

# Constructor de convenienica para número máximo de cambios de base en cada
# VarCPIBase
InflationTotalRebaseCPI(period::Int) = InflationTotalRebaseCPI(period, 0)

# Nombre de la medida
measure_name(inflfn::InflationTotalRebaseCPI) =
    "Variación interanual IPC con cambios de base sintéticos ($(inflfn.period), $(inflfn.maxchanges))"

# Parámetros
params(totalrebasefn::InflationTotalRebaseCPI) = (totalrebasefn.period, )

# Computar variación intermensual resumen de medida de inflación aplicando
# metodología de cambio de base sintético
function (totalrebasefn::InflationTotalRebaseCPI)(base::VarCPIBase)

    # Número de períodos para realizar cambio de base
    period = totalrebasefn.period

    # Obtener períodos de la base
    T = periods(base)

    # Obtener vector de índices
    startidxs, endidxs = _getranges(T, period, totalrebasefn.maxchanges)

    # Mapear cada rango de índices en el resumen intermensual obtenido con
    # fórmula del IPC
    vinterm = mapreduce(vcat, startidxs, endidxs) do startidx, endidx
        cpi = capitalize(view(base.v, startidx:endidx, :)) * base.w / 100
        varinterm(cpi)
    end

    # Resumen intermensual es la concatenación vertical del resumen intermensual
    # en las distintas bases sintéticas
    vinterm

end

# Función de ayuda para obtener vector de rangos
function _getranges(T, period, max_changes=0)

    # Si hay menos observaciones que el período, devolver el único rango
    T <= period && return 1, T

    # max_changes == 0 => todos los cambios de base posibles
    if max_changes == 0
        blocks = cld(T, period)
    else
        blocks = min(cld(T, period), max_changes+1)
    end

    startidxs = 1:period:T # índices iniciales
    finalidxs = Vector{Int}(undef, blocks) # vector de índices finales
    c = 0 # contador de cambios de base

    # Completar los índices finales
    for j in 1:blocks
        finalidxs[j] = startidxs[j] + period - 1
        c += 1
        if finalidxs[j] >= T || (max_changes != 0 && c > max_changes)
            finalidxs[j] = T
            break
        end
    end

    startidxs[1:blocks], finalidxs
end


text:
###############################################################################
#
# The Stopping structure eases the implementation of algorithms and the
# stopping criterion.
#
# The following examples illustrate solver for linear algebra:
#  Ax = b with A an m x n matrix.
#
# This tutorial illustrates the different step in preparing the resolution of a
# new problem.
# - we create a LinearAlgebraPRoblem (that stores A, b)
# - we use the GennericState storing x and the current_time
# - we create a LinearAlgebraStopping
# - the optimality function linear_system_check!
#
# A more sophisticated version can be found in LAStopping.
#
###############################################################################
using LinearAlgebra, Stopping, Test

m, n = 400, 200 #size of A: m x n
A    = 100 * rand(m, n)
xref = 100 * rand(n)
b    = A * xref

#Our initial guess
x0 = zeros(n)

mutable struct LinearAlgebraProblem
    A :: Any #matrix type
    b :: Vector
end

la_pb = LinearAlgebraProblem(A, b)
la_state = GenericState(xref)

@test norm(la_pb.A * xref - la_pb.b) <= 1e-6

mutable struct LinearAlgebraStopping <: AbstractStopping

        # problem
        pb :: LinearAlgebraProblem

        # Common parameters
        meta :: AbstractStoppingMeta

        # current state of the problem
        current_state :: AbstractState

        # Stopping of the main problem, or nothing
        main_stp :: Union{AbstractStopping, Nothing}

        function LinearAlgebraStopping(pb               :: LinearAlgebraProblem,
                                       current_state    :: AbstractState; kwargs...)
         return new(pb, linear_system_check!, StoppingMeta(;optimality_check = linear_system_check, kwargs...), la_state, nothing)
        end
end

function linear_system_check(pb    :: LinearAlgebraProblem,
                             state :: AbstractState; kwargs...)
 return norm(pb.A * state.x - pb.b)
end

@test linear_system_check(la_pb, la_state) == 0.0
update!(la_state, x = x0)
@test linear_system_check(la_pb, la_state) != 0.0

la_stop = LinearAlgebraStopping(la_pb, la_state,
                                max_iter = 150000, rtol = 1e-6, optimality_check = linear_system_check)

"""
Randomized block Kaczmarz
"""
function RandomizedBlockKaczmarz(stp :: AbstractStopping; kwargs...)

    #A,b = stp.current_state.Jx, stp.current_state.cx
    A,b = stp.pb.A, stp.pb.b
    x0  = stp.current_state.x

    m,n = size(A)
    xk  = x0

    OK = start!(stp)

    while !OK

     i  = Int(floor(rand() * m)+1) #rand a number between 1 and m
     Ai = A[i,:]
     xk  = Ai == 0 ? x0 : x0 - (dot(Ai,x0)-b[i])/dot(Ai,Ai) * Ai

     OK = update_and_stop!(stp, x = xk)
     x0  = xk

    end

 return stp
end

RandomizedBlockKaczmarz(la_stop)
@test status(la_stop) == :Optimal


text:
using LinearAlgebra

function sum_dim1(A::CuSparseMatrixCSR)
    function kernel(Tnorm, out, dA)
        idx = (blockIdx().x-1) * blockDim().x + threadIdx().x
        idx < length(dA.rowPtr) || return
        s = zero(Tnorm)
        for k in dA.rowPtr[idx]:dA.rowPtr[idx+1]-1
            s += abs(dA.nzVal[k])
        end
        out[idx] = s
        return
    end

    m, n = size(A)
    Tnorm = typeof(float(real(zero(eltype(A)))))
    Tsum = promote_type(Float64,Tnorm)
    rowsum = CUDA.CuArray{Tsum}(undef, m)
    kernel_f = @cuda launch=false kernel(Tnorm, rowsum, A)

    config = launch_configuration(kernel_f.fun)
    threads = min(n, config.threads)
    blocks = cld(n, threads)
    kernel_f(Tnorm, rowsum, A; threads, blocks)
    return rowsum
end

function LinearAlgebra.opnorm(A::CuSparseMatrixCSR, p::Real=2)
    if p == Inf
        return maximum(sum_dim1(A))
    else
        error("p=$p is not supported")
    end
end

LinearAlgebra.opnorm(A::CuSparseMatrixCSC, p::Real=2) = opnorm(CuSparseMatrixCSR(A), p)


text:
# # Solve the Stokes equation of colliding flow: Q1-Q0 element, general formulation

# Synopsis: Compute the solution of the Stokes equation of two-dimensional
# incompressible viscous flow for a manufactured problem of colliding flow.
# Continuous velocity/discontinuous pressure quadrilateral elements are used.

# The "manufactured" colliding flow example from Elman et al 2014. The
# Continuous velocity/discontinuous pressure formulation with linear
# quadrilaterals. These elements suffer from pressure oscillation instability.
# The results of this simulation demonstrate it.

# The analytical results are shown here: pressure is shown  with contours,
# and the velocities visualized with arrows at random points.
# ![Pressure and velocity](colliding.png)

# The formulation is the general elasticity-like scheme with
# strain-rate/velocity matrices. It can be manipulated into the one derived
# in Reddy, Introduction to the finite element method, 1993. Page 486 ff.

# The complete code is in the file [`tut_stokes_q1_q0_gen.jl`](tut_stokes_q1_q0_gen.jl).

# The solution will be defined  within a module in order to eliminate conflicts
# with data or functions defined elsewhere.

module tut_stokes_q1_q0_gen

# We'll need some functionality from linear algebra, static arrays, and the mesh
# libraries. Some plotting will be produced to visualize structure of the
# stiffness matrix. Finally we will need the `Elfel` functionality.
using LinearAlgebra
using StaticArrays
using MeshCore.Exports
using MeshSteward.Exports
using Elfel.Exports
using UnicodePlots

# The boundary value problem is expressed in this weak form
# ```math
#  \int_{V}{\underline{\varepsilon}}(\underline{\delta v})^T\;
#  \underline{\underline{D}}\; {\underline{\varepsilon}}(\underline{u})\; \mathrm{d} V
# - \int_{V} \mathrm{div}(\underline{\delta v})\; p\; \mathrm{d} V = 0,\quad \forall \underline{\delta v}
# ```
# ```math
#  - \int_{V} \delta q\; \mathrm{div}(\underline{u}) \; \mathrm{d} V = 0,\quad \forall \delta q
# ```
# Here ``\underline{\delta v}`` are the test functions in the velocity space,
# and ``\delta q`` are the pressure test functions. Further ``\underline
# {u}`` is the trial velocity, and ``p`` is the trial pressure.

function run()
    mu = 1.0 # dynamic viscosity
    # This is the material-property matrix ``\underline{\underline{D}}``:
    D = SMatrix{3, 3}(
        [2*mu  0   0
          0  2*mu  0
          0    0   mu])
    A = 1.0 # half of the length of the side of the square
    N = 100 # number of element edges per side of the square
    # These three functions define the true velocity components and the true
    # pressure.
    trueux = (x, y) -> 20 * x * y ^ 3
    trueuy = (x, y) -> 5 * x ^ 4 - 5 * y ^ 4
    truep = (x, y) -> 60 * x ^ 2 * y - 20 * y ^ 3

    # Construct the two meshes for the mixed method. They need to support the
    # velocity and pressure spaces.
    mesh = genmesh(A, N)

    # Construct the velocity space: it is a vector space with two components. The
    # degrees of freedom are real numbers (`Float64`). The velocity mesh
    # carries the finite elements of  the continuity ``H ^1``, i. e. both the
    # function values and the derivatives are square integrable. Each node
    # carries 2 degrees of freedom, hence there are two velocity components per
    # node.
    Uh = FESpace(Float64, mesh, FEH1_Q4(), 2)

    # Now we apply the boundary conditions at the nodes around the
    # circumference.
    locs = geometry(mesh)
    # We use searching based on the presence of the node within a box. The
    # entire boundary will be captured within these four boxes, provided we
    # inflate those boxes with a little tolerance (we can't rely on those
    # nodes to be precisely at the coordinates given, we need to introduce some
    # tolerance).
    boxes = [[-A A -A -A], [-A -A -A A], [A A -A A], [-A A A A]]
    inflate = A / N / 100
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            # Remember that all  components of the velocity are known at the
            # boundary.
            setebc!(Uh, 0, i, 1, trueux(locs[i]...))
            setebc!(Uh, 0, i, 2, trueuy(locs[i]...))
        end
    end

    # No we construct the pressure space. It is a discontinuous, piecewise
    # constant space supported on a mesh of 4-node quadrilaterals
    # (``L ^2``-continuity elements). Each quadrilateral carries a single
    # degree of freedom.
    Ph = FESpace(Float64, mesh, FEL2_Q4(), 1)

    # The pressure in this "enclosed" flow example is only known up to a
    # constant. By setting  pressure degree of freedom at one degree of freedom
    # to zero will make the solution unique.
    atcenter = vselect(geometry(mesh); nearestto = [0.0, 0.0])
    setebc!(Ph, 2, atcenter[1], 1, 0.0)

    # Number the degrees of freedom. First all the free degrees of freedom are
    # numbered, both velocities and pressures. Next all the data degrees of
    # freedom are numbered, again both for the velocities and for the
    # pressures.
    numberdofs!([Uh, Ph])
    # The total number of degrees of freedom is now calculated.
    tndof = ndofs(Uh) + ndofs(Ph)
    # As is the total number of unknowns.
    tnunk = nunknowns(Uh) + nunknowns(Ph)

    # Assemble the coefficient matrix.
    K = assembleK(Uh, Ph, tndof, D)

    # Display the structure of the indefinite stiffness matrix. Note that this
    # is the complete matrix, including rows and columns for all the degrees of
    # freedom, unknown and known.
    p = spy(K, canvas = DotCanvas)
    display(p)

    # Solve the linear algebraic system. First construct system vector of all
    # the degrees of freedom, in the first `tnunk` rows that corresponds to the
    # unknowns, and the subsequent rows are for the data degrees of freedom.
    U = fill(0.0, tndof)
    gathersysvec!(U, [Uh, Ph])
    # Note that  the vector `U` consists of nonzero numbers in rows are for the
    # data degrees of freedom. Multiplying the stiffness matrix with this
    # vector will generate a load vector  on the right-hand side. Otherwise
    # there is no loading, hence the vector `F` consists of all zeros.
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    # Once we have solved the system of linear equations, we can distribute the
    # solution from the vector `U` into the finite element spaces.
    scattersysvec!([Uh, Ph], U)

    # Given that the solution is manufactured, i. e. exactly known, we can
    # calculate the true errors.
    @show ep = evaluate_pressure_error(Ph, truep)
    @show ev = evaluate_velocity_error(Uh, trueux, trueuy)

    # Postprocessing. First we make attributes, scalar nodal attributes,
    # associated with the meshes for the pressures and the velocity.
    makeattribute(Ph, "p", 1)
    makeattribute(Uh, "ux", 1)
    makeattribute(Uh, "uy", 2)
    # The pressure and the velocity components are then written out into two VTK
    # files. The pressure is of the "cell-data" type, the velocity is of
    # the "point-data" type.
    vtkwrite("tut_stokes_q1_q0_gen-p", baseincrel(mesh), [(name = "p",), ])
    vtkwrite("tut_stokes_q1_q0_gen-v", baseincrel(mesh), [(name = "ux",), (name = "uy",)])

    return true
end

function genmesh(A, N)
    # Linear quadrilateral mesh is used for both the velocity space and the
    # pressure space.
    mesh = attach!(Mesh(), Q4block(2 * A, 2 * A, N, N), "velocity")
    # Now translate so that the center of the square is at the origin of the
    # coordinates.
    ir = baseincrel(mesh)
    transform(ir, x -> x .- A)
    # The pressure degree of freedom is associated with the element itself. The
    # mesh will therefore need be equipped with the incidence relation ``
    # (2, 2)``. The finite element space for the velocity will  have degrees of
    # freedom associated with the vertices. The finite element space for the
    # pressure will have degrees of freedom associated with  the faces
    # (elements themselves). The finite element space does that by associating
    # the pressure field with the ``(2, 2)`` incidence relation.
    eidir = ir_identity(ir)
    attach!(mesh, eidir)
    return mesh
end

function assembleK(Uh, Ph, tndof, D)
    function integrate!(ass, elits, qpits, D)
        # Consider the elementwise definition of the test strain rate, ``
        # {\underline{\varepsilon}}(\underline{\delta v})``. It is calculated
        # from the elementwise degrees of freedom and the associated basis
        # functions  as
        # ```math
        # {\underline{\varepsilon}}(\underline{\delta v}) =
        #  \sum_i{\delta V}_i {\underline{B}_{c(i)}(N_i)}
        # ```
        # where ``i = 1, \ldots, n_{du}``, and ``n_{du}`` is the number of
        # velocity degrees of freedom per element, ``c(i)`` is the number of
        # the component corresponding to the degree of freedom ``i``. This is
        # either 1 when degree of freedom ``i`` is the ``x``-component of the
        # velocity, 2 otherwise(for the ``y``-component of the velocity).
        # Analogously for the trial strain rate.

        # The strain-rate/velocity matrices are defined as
        # ```math
        # {\underline{B}_{1}(N_i)} =
        # \left[\begin{array}{c}
        #      \partial{N_i}/\partial{x}  \\
        #      0 \\
        #      \partial{N_i}/\partial{y}
        # \end{array}\right],
        # ```
        # and
        # ```math
        # {\underline{B}_{2}(N_i)} =
        # \left[\begin{array}{c}
        #      0 \\
        #      \partial{N_i}/\partial{y}  \\
        #      \partial{N_i}/\partial{x}
        # \end{array}\right].
        # ```
        # This tiny function evaluates the strain-rate/velocity matrices defined above
        # from the gradient of a basis function and the given number of the
        # component corresponding to the current degree of freedom.
        B = (g, k) -> (k == 1 ?
            SVector{3}((g[1], 0, g[2])) :
            SVector{3}((0, g[2], g[1])))
        # This array defines the components for the element degrees of freedom,
        # as defined above as ``c(i)``.
        c = edofcompnt(Uh)
        # These are the totals of the velocity and pressure degrees of freedom
        # per element.
        n_du, n_dp = ndofsperel.((Uh, Ph))
        # The local matrix assemblers are used as if they were ordinary
        # elementwise dense matrices. Here they are defined.
        kuu = LocalMatrixAssembler(n_du, n_du, 0.0)
        kup = LocalMatrixAssembler(n_du, n_dp, 0.0)
        for el in zip(elits...)
            uel, pel = el
            # The local matrix assemblers are initialized with zeros for the
            # values, and with the element degree of freedom vectors to be used
            # in the assembly. The assembler `kuu` is used for the velocity
            # degrees of freedom, and the assembler `kup` collect the coupling
            # coefficients between the velocity and the pressure. The function
            # `eldofs` collects the global numbers of the degrees of freedom
            # either for the velocity space, or for the pressure space
            # (`eldofs(pel)`).
            init!(kuu, eldofs(uel), eldofs(uel))
            init!(kup, eldofs(uel), eldofs(pel))
            for qp in zip(qpits...)
                uqp, pqp = qp
                # The integration is performed using the velocity quadrature points.
                Jac, J = jacjac(uel, uqp)
                JxW = J * weight(uqp)
                gradNu = bfungrad(uqp, Jac) # gradients of the velocity basis functions
                Np = bfun(pqp) # pressure basis functions
                # This double loop corresponds precisely to the integrals of the
                # weak form. This is the  matrix in the upper left corner.
                for i in 1:n_du
                    DBi = D * B(gradNu[i], c[i])
                    for j in 1:n_du
                        Bj = B(gradNu[j], c[j])
                        kuu[j, i] += dot(Bj, DBi) * (JxW)
                    end
                end
                # And this is the coupling matrix in the top right corner.
                for i in 1:n_dp, j in 1:n_du
                    kup[j, i] += gradNu[j][c[j]] * (-JxW * Np[i])
                end
            end
            # Assemble the matrices. The submatrix off the diagonal is assembled
            # twice, once as itself, and once as its transpose.
            assemble!(ass, kuu)
            assemble!(ass, kup) # top right corner
            assemble!(ass, transpose(kup)) # bottom left corner
        end
        return ass # return the updated assembler of the global matrix
    end

    # In the `assembleK` function we first we create the element iterators. We
    # can go through all the elements, both in the velocity finite element
    # space and in the pressure finite element space, that define the domain of
    # integration using this iterator. Each time a new element is accessed,
    # some data are precomputed such as the element degrees of freedom,
    # components of the degree of freedom, etc. Note that we need to iterate
    # two finite element spaces, hence we create a tuple of iterators.
    elits = (FEIterator(Uh), FEIterator(Ph))
    # These are the quadrature point iterators. We know that the elements are
    # triangular. We choose the three-point rule, to capture the quadratic
    # component in the velocity space. Quadrature-point iterators provide
    # access to basis function values and gradients, the Jacobian matrix and
    # the Jacobian determinant, the location of the quadrature point and so
    # on. Note that we need to iterate the quadrature rules of
    # two finite element spaces, hence we create a tuple of iterators.
    qargs = (kind = :Gauss, order = 2)
    qpits = (QPIterator(Uh, qargs), QPIterator(Ph, qargs))
    # The matrix will be assembled into this assembler. Which is initialized
    # with the total number of degrees of freedom (dimension of the coefficient
    # matrix before partitioning into unknowns and data degrees of freedom).
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    # The integration is carried out, and then...
    integrate!(ass, elits, qpits, D)
    # ...we materialize the sparse stiffness matrix and return it.
    return finish!(ass)
end

# The linear algebraic system is solved by partitioning. The vector `U` is
# initially all zero, except in the degrees of freedom which are prescribed as
# nonzero. Therefore the product of the stiffness matrix and the vector `U`
# are the loads due to nonzero essential boundary conditions.  The
# submatrix of the stiffness conduction matrix corresponding to the free degrees of
# freedom (unknowns), `K[1:nu, 1:nu]` is then used to solve for the unknowns `U
# [1:nu]`.
function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

# The function `evaluate_pressure_error` evaluates the true ``L^2``  error of
# the pressure. It does that by integrating the square of the difference
# between the approximate pressure and the  true pressure, the true pressure
# being provided by the `truep` function.
function evaluate_pressure_error(Ph, truep)
    function integrate!(elit, qpit, truep)
        n_dp = ndofsperel(elit)
        E = 0.0
        for el in elit
            dofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Np = bfun(qp)
                pt = truep(location(el, qp)...)
                pa = 0.0
                for j in 1:n_dp
                    pa += (dofvals[j] * Np[j])
                end
                E += (JxW) * (pa - pt)^2
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(Ph)
    qargs = (kind = :Gauss, order = 2)
    qpit = QPIterator(Ph, qargs)
    return integrate!(elit, qpit, truep)
end

# The function `evaluate_velocity_error` evaluates the true ``L^2``  error of
# the velocity. It does that by integrating the square of the difference
# between the approximate pressure and the  true velocity, the true velocity
# being provided by the `trueux`, `trueuy` functions.
function evaluate_velocity_error(Uh, trueux, trueuy)
    function integrate!(elit, qpit, trueux, trueuy)
        n_du = ndofsperel(elit)
        uedofcomp = edofcompnt(Uh)
        E = 0.0
        for el in elit
            udofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Nu = bfun(qp)
                uxt = trueux(location(el, qp)...)
                uyt = trueuy(location(el, qp)...)
                uxa = 0.0
                uya = 0.0
                for j in 1:n_du
                    (uedofcomp[j] == 1) && (uxa += (udofvals[j] * Nu[j]))
                    (uedofcomp[j] == 2) && (uya += (udofvals[j] * Nu[j]))
                end
                E += (JxW) * ((uxa - uxt)^2 + (uya - uyt)^2)
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(Uh)
    qargs = (kind = :Gauss, order = 2)
    qpit = QPIterator(Uh, qargs)
    return integrate!(elit, qpit, trueux, trueuy)
end

end

# To run the example, evaluate this file which will  compile the module
# `.tut_stokes_q1_q0_gen`.
using .tut_stokes_q1_q0_gen
tut_stokes_q1_q0_gen.run()


text:
############################################################
## joAbstractOperator - overloaded Base functions

# eltype(jo)
eltype(A::joAbstractOperator) = jo_method_error(A,"eltype(jo) not implemented")

# deltype(jo)
deltype(A::joAbstractOperator) = jo_method_error(A,"deltype(jo) not implemented")

# reltype(jo)
reltype(A::joAbstractOperator) = jo_method_error(A,"reltype(jo) not implemented")

# show(jo)
show(A::joAbstractOperator) = jo_method_error(A,"show(jo) not implemented")

# display(jo)
display(A::joAbstractOperator) = jo_method_error(A,"display(jo) not implemented")

# size(jo)
size(A::joAbstractOperator) = jo_method_error(A,"size(jo) not implemented")

# size(jo,1/2)
size(A::joAbstractOperator,ind) = jo_method_error(A,"size(jo,1/2) not implemented")

# length(jo)
length(A::joAbstractOperator) = jo_method_error(A,"length(jo) not implemented")

# jo_full(jo)
jo_full(A::joAbstractOperator) = jo_method_error(A,"jo_full(jo) not implemented")

# norm(jo)
norm(A::joAbstractOperator,p) = jo_method_error(A,"norm(jo) not implemented")

# real(jo)
real(A::joAbstractOperator) = jo_method_error(A,"real(jo) not implemented")

# imag(jo)
imag(A::joAbstractOperator) = jo_method_error(A,"imag(jo) not implemented")

# conj(jo)
conj(A::joAbstractOperator) = jo_method_error(A,"conj(jo) not implemented")

# transpose(jo)
transpose(A::joAbstractOperator) = jo_method_error(A,"transpose(jo) not implemented")

# adjoint(jo)
adjoint(A::joAbstractOperator) = jo_method_error(A,"adjoint(jo) not implemented")

# isreal(jo)
isreal(A::joAbstractOperator) = jo_method_error(A,"isreal(jo) not implemented")

# issymmetric(jo)
issymmetric(A::joAbstractOperator) = jo_method_error(A,"issymmetric(jo) not implemented")

# ishermitian(jo)
ishermitian(A::joAbstractOperator) = jo_method_error(A,"ishermitian(jo) not implemented")

# isequal(jo,jo)
isequal(A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,"isequal(jo,jo) not implemented or type mismatch")

# isapprox(jo,jo)
isapprox(A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,"isapprox(jo,jo) not implemented or type mismatch")

# getindex(jo,...)
getindex(A::joAbstractOperator,args...) = jo_method_error(A,B,"getindex(jo,...) not implemented or type mismatch")

############################################################
## overloaded Base *(...jo...)

# *(jo,jo)
*(A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,"*(jo,jo) not implemented or type mismatch")

# *(jo,mvec)
*(A::joAbstractOperator,mv::AbstractMatrix) = jo_method_error(A,mv,"*(jo,mvec) not implemented or type mismatch")

# *(mvec,jo)
*(mv::AbstractMatrix,A::joAbstractOperator) = jo_method_error(mv,A,"*(mvec,jo) not implemented or type mismatch")

# *(jo,vec)
*(A::joAbstractOperator,v::AbstractVector) = jo_method_error(A,v,"*(jo,vec) not implemented or type mismatch")

# *(vec,jo)
*(v::AbstractVector,A::joAbstractOperator) = jo_method_error(v,A,"*(vec,jo) not implemented or type mismatch")

# *(num,jo)
*(a,A::joAbstractOperator) = jo_method_error(a,A,"*(any,jo) not implemented or type mismatch")

# *(jo,num)
*(A::joAbstractOperator,a) = jo_method_error(A,a,"*(jo,any) not implemented or type mismatch")

############################################################
## overloaded Base \(...jo...)

# \(jo,jo)
\(A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,"\\(jo,jo) not implemented or type mismatch")

# \(jo,mvec)
\(A::joAbstractOperator,mv::AbstractMatrix) = jo_method_error(A,mv,"\\(jo,mvec) not implemented or type mismatch")

# \(mvec,jo)
\(mv::AbstractMatrix,A::joAbstractOperator) = jo_method_error(mv,A,"\\(mvec,jo) not implemented or type mismatch")

# \(jo,vec)
\(A::joAbstractOperator,v::AbstractVector) = jo_method_error(A,v,"\\(jo,vec) not implemented or type mismatch")

# \(vec,jo)
\(v::AbstractVector,A::joAbstractOperator) = jo_method_error(v,A,"\\(vec,jo) not implemented or type mismatch")

# \(num,jo)
\(a,A::joAbstractOperator) = jo_method_error(a,A,"\\(any,jo) not implemented or type mismatch")

# \(jo,num)
\(A::joAbstractOperator,a) = jo_method_error(A,a,"\\(jo,any) not implemented or type mismatch")

############################################################
## overloaded Base +(...jo...)

# +(jo)
+(A::joAbstractOperator) = jo_method_error(A,"+(jo) not implemented")

# +(jo,jo)
+(A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,"+(jo,jo) not implemented or type mismatch")

# +(jo,mvec)
+(A::joAbstractOperator,mv::AbstractMatrix) = jo_method_error(A,mv,"+(jo,mvec) not implemented or type mismatch")

# +(mvec,jo)
+(mv::AbstractMatrix,A::joAbstractOperator) = jo_method_error(mv,A,"+(mvec,jo) not implemented or type mismatch")

# +(jo,vec)
+(A::joAbstractOperator,v::AbstractVector) = jo_method_error(A,v,"+(jo,vec) not implemented or type mismatch")

# +(vec,jo)
+(v::AbstractVector,A::joAbstractOperator) = jo_method_error(v,A,"+(vec,jo) not implemented or type mismatch")

# +(jo,num)
+(A::joAbstractOperator,b) = jo_method_error(A,b,"+(jo,any) not implemented or type mismatch")

# +(num,jo)
+(b,A::joAbstractOperator) = jo_method_error(b,A,"+(any,jo) not implemented or type mismatch")

############################################################
## overloaded Base -(...jo...)

# -(jo)
-(A::joAbstractOperator) = jo_method_error(A,"-(jo) not implemented")

# -(jo,jo)
-(A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,"-(jo,jo) not implemented or type mismatch")

# -(jo,mvec)
-(A::joAbstractOperator,mv::AbstractMatrix) = jo_method_error(A,mv,"-(jo,mvec) not implemented or type mismatch")

# -(mvec,jo)
-(mv::AbstractMatrix,A::joAbstractOperator) = jo_method_error(mv,A,"-(mvec,jo) not implemented or type mismatch")

# -(jo,vec)
-(A::joAbstractOperator,v::AbstractVector) = jo_method_error(A,v,"-(jo,vec) not implemented or type mismatch")

# -(vec,jo)
-(v::AbstractVector,A::joAbstractOperator) = jo_method_error(v,A,"-(vec,jo) not implemented or type mismatch")

# -(jo,num)
-(A::joAbstractOperator,b) = jo_method_error(A,b,"-(jo,any) not implemented or type mismatch")

# -(num,jo)
-(b,A::joAbstractOperator) = jo_method_error(b,A,"-(any,jo) not implemented or type mismatch")

############################################################
## overloaded Base .*(...jo...)
## function Base.Broadcast.broadcasted(::typeof(*), ...)

# .*(jo,jo)
Base.Broadcast.broadcasted(::typeof(*),A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,".*(jo,jo) not implemented")

# .*(jo,mvec)
Base.Broadcast.broadcasted(::typeof(*),A::joAbstractOperator,mv::AbstractMatrix) = jo_method_error(A,mv,".*(jo,mvec) not implemented")

# .*(mvec,jo)
Base.Broadcast.broadcasted(::typeof(*),mv::AbstractMatrix,A::joAbstractOperator) = jo_method_error(mv,A,".*(mvec,jo) not implemented")

# .*(jo,vec)
Base.Broadcast.broadcasted(::typeof(*),A::joAbstractOperator,v::AbstractVector) = jo_method_error(A,v,".*(jo,vec) not implemented")

# .*(vec,jo)
Base.Broadcast.broadcasted(::typeof(*),v::AbstractVector,A::joAbstractOperator) = jo_method_error(v,A,".*(vec,jo) not implemented")

# .*(num,jo)
Base.Broadcast.broadcasted(::typeof(*),a,A::joAbstractOperator) = jo_method_error(a,A,".*(any,jo) not implemented")

# .*(jo,num)
Base.Broadcast.broadcasted(::typeof(*),A::joAbstractOperator,a) = jo_method_error(A,a,".*(jo,any) not implemented")

############################################################
## overloaded Base .\(...jo...)
## function Base.Broadcast.broadcasted(::typeof(\), ...)

# .\(jo,jo)
Base.Broadcast.broadcasted(::typeof(\),A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,".(jo,jo) not implemented")

# .\(jo,mvec)
Base.Broadcast.broadcasted(::typeof(\),A::joAbstractOperator,mv::AbstractMatrix) = jo_method_error(A,mv,".(jo,mvec) not implemented")

# .\(mvec,jo)
Base.Broadcast.broadcasted(::typeof(\),mv::AbstractMatrix,A::joAbstractOperator) = jo_method_error(mv,A,".(mvec,jo) not implemented")

# .\(jo,vec)
Base.Broadcast.broadcasted(::typeof(\),A::joAbstractOperator,v::AbstractVector) = jo_method_error(A,v,".(jo,vec) not implemented")

# .\(vec,jo)
Base.Broadcast.broadcasted(::typeof(\),v::AbstractVector,A::joAbstractOperator) = jo_method_error(v,A,".(vec,jo) not implemented")

# .\(num,jo)
Base.Broadcast.broadcasted(::typeof(\),a,A::joAbstractOperator) = jo_method_error(a,A,".(any,jo) not implemented")

# .\(jo,num)
Base.Broadcast.broadcasted(::typeof(\),A::joAbstractOperator,a) = jo_method_error(A,a,".(jo,any) not implemented")

############################################################
## overloaded Base .+(...jo...)
## function Base.Broadcast.broadcasted(::typeof(+), ...)

# .+(jo,jo)
Base.Broadcast.broadcasted(::typeof(+),A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,".+(jo,jo) not implemented")

# .+(jo,mvec)
Base.Broadcast.broadcasted(::typeof(+),A::joAbstractOperator,mv::AbstractMatrix) = jo_method_error(A,mv,".+(jo,mvec) not implemented")

# .+(mvec,jo)
Base.Broadcast.broadcasted(::typeof(+),mv::AbstractMatrix,A::joAbstractOperator) = jo_method_error(mv,A,".+(mvec,jo) not implemented")

# .+(jo,vec)
Base.Broadcast.broadcasted(::typeof(+),A::joAbstractOperator,v::AbstractVector) = jo_method_error(A,v,".+(jo,vec) not implemented")

# .+(vec,jo)
Base.Broadcast.broadcasted(::typeof(+),v::AbstractVector,A::joAbstractOperator) = jo_method_error(v,A,".+(vec,jo) not implemented")

# .+(jo,num)
Base.Broadcast.broadcasted(::typeof(+),A::joAbstractOperator,b) = jo_method_error(A,b,".+(jo,any) not implemented")

# .+(num,jo)
Base.Broadcast.broadcasted(::typeof(+),b,A::joAbstractOperator) = jo_method_error(b,A,".+(any,jo) not implemented")

############################################################
## overloaded Base .-(...jo...)
## function Base.Broadcast.broadcasted(::typeof(-), ...)

# .-(jo,jo)
Base.Broadcast.broadcasted(::typeof(-),A::joAbstractOperator,B::joAbstractOperator) = jo_method_error(A,B,".-(jo,jo) not implemented")

# .-(jo,mvec)
Base.Broadcast.broadcasted(::typeof(-),A::joAbstractOperator,mv::AbstractMatrix) = jo_method_error(A,mv,".-(jo,mvec) not implemented")

# .-(mvec,jo)
Base.Broadcast.broadcasted(::typeof(-),mv::AbstractMatrix,A::joAbstractOperator) = jo_method_error(mv,A,".-(mvec,jo) not implemented")

# .-(jo,vec)
Base.Broadcast.broadcasted(::typeof(-),A::joAbstractOperator,v::AbstractVector) = jo_method_error(A,v,".-(jo,vec) not implemented")

# .-(vec,jo)
Base.Broadcast.broadcasted(::typeof(-),v::AbstractVector,A::joAbstractOperator) = jo_method_error(v,A,".-(vec,jo) not implemented")

# .-(jo,num)
Base.Broadcast.broadcasted(::typeof(-),A::joAbstractOperator,b) = jo_method_error(A,b,".-(jo,any) not implemented")

# .-(num,jo)
Base.Broadcast.broadcasted(::typeof(-),b,A::joAbstractOperator) = jo_method_error(b,A,".-(any,jo) not implemented")

############################################################
## overloaded Base block methods

# hcat(...jo...)
hcat(ops::joAbstractOperator...) = jo_method_error(ops[1],"hcat(jo...) not implemented")

# vcat(...jo...)
vcat(ops::joAbstractOperator...) = jo_method_error(ops[1],"vcat(jo...) not implemented")

# hvcat(...jo...)
hvcat(rows::Tuple{Vararg{Int}}, ops::joAbstractOperator...) = jo_method_error(ops[1],"hvcat(jo...) not implemented")

############################################################
## overloaded LinearAlgebra functions

# mul!(...,jo,...)
mul!(y::AbstractVector,A::joAbstractOperator,x::AbstractVector)=jo_method_error(y,A,x,"mul!(vec,jo,vec) not implemented or type mismatch")
mul!(y::AbstractMatrix,A::joAbstractOperator,x::AbstractMatrix)=jo_method_error(y,A,x,"mul!(mat,jo,mat) not implemented or type mismatch")

# ldiv!(...,jo,...)
ldiv!(y::AbstractVector,A::joAbstractOperator,x::AbstractVector)=jo_method_error(y,A,x,"ldiv!(vec,jo,vec) not implemented or type mismatch")
ldiv!(y::AbstractMatrix,A::joAbstractOperator,x::AbstractMatrix)=jo_method_error(y,A,x,"ldiv!(mat,jo,mat) not implemented or type mismatch")


text:
using LinearAlgebra;
using SparseArrays;
using IterativeSolvers;
using Krylov;
using LinearOperators;

function LaLdlt!(vX, mP, mA, vV, ρ = 1e6, σ = 1e-6;
    numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500)

    # Direct Solver using Base LDLT

    numElementsX            = size(vX, 1);
    numRowsP, numColsP      = size(mP);
    numRowsA, numColsA      = size(mA);

    ρ¹ = 1 / ρ;

    hDL = ldlt([mP + (σ * sparse(I, numElementsX, numElementsX)) mA'; mA -ρ¹ * sparse(I, numRowsA, numRowsA)]);

    for ii in 1:(numFactor - 1)
        hDL = ldlt([mP + (σ * sparse(I, numElementsX, numElementsX)) mA'; mA -ρ¹ * sparse(I, numRowsA, numRowsA)]);
    end

    for ii in 1:numIterations
        vV .= hDL \ vV;
        vX .= vV[1:numElements];
    end

end

function ItrSolCg!(vX, mP, mA, vV, ρ = 1e6, σ = 1e-6;
    numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500)

    # Iterative Solver using IterativeSolver cg!()

    numElementsX            = size(vX, 1);
    numRowsP, numColsP      = size(mP);
    numRowsA, numColsA      = size(mA);

    ρ¹ = 1 / ρ;

    mL = mP + (σ * sparse(I, numElementsX, numElementsX)) + (ρ * (mA' * mA));

    for ii in 1:(numFactor - 1)
        mL = mP + (σ * sparse(I, numElementsX, numElementsX)) + (ρ * (mA' * mA));
    end

    for ii in 1:numIterations
        IterativeSolvers.cg!(vX, mL, vX, abstol = ϵSolver, maxiter = numItrSolver); #<! Also compatible with LinearOperators
    end

end

function KrylovCg!(vX, mP, mA, vV, ρ = 1e6, σ = 1e-6;
    numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500)

    # Iterative Solver using IterativeSolver cg!()

    numElementsX            = size(vX, 1);
    numRowsP, numColsP      = size(mP);
    numRowsA, numColsA      = size(mA);

    ρ¹ = 1 / ρ;

    mL = mP + (σ * sparse(I, numElementsX, numElementsX)) + (ρ * (mA' * mA));

    for ii in 1:(numFactor - 1)
        mL = mP + (σ * sparse(I, numElementsX, numElementsX)) + (ρ * (mA' * mA));
    end

    sCgSolver = Krylov.CgSolver(mL, vX);

    for ii in 1:numIterations
        # vT, _ = Krylov.cg(mL, vX, atol = ϵSolver, itmax = numItrSolver);
        vT, _ = Krylov.cg!(sCgSolver, mL, vX, atol = ϵSolver, itmax = numItrSolver);
        vX .= vT;
    end

end

function KrylovCr!(vX, mP, mA, vV, ρ = 1e6, σ = 1e-6;
    numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500)

    # Iterative Solver using IterativeSolver cg!()

    numElementsX            = size(vX, 1);
    numRowsP, numColsP      = size(mP);
    numRowsA, numColsA      = size(mA);

    ρ¹ = 1 / ρ;

    mL = mP + (σ * sparse(I, numElementsX, numElementsX)) + (ρ * (mA' * mA));

    for ii in 1:(numFactor - 1)
        mL = mP + (σ * sparse(I, numElementsX, numElementsX)) + (ρ * (mA' * mA));
    end

    for ii in 1:numIterations
        vT, _ = Krylov.cr(mL, vX, atol = ϵSolver, itmax = numItrSolver);
        vX .= vT;
    end

end

function KrylovCgLanczos!(vX, mP, mA, vV, ρ = 1e6, σ = 1e-6;
    numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500)

    # Iterative Solver using IterativeSolver cg!()

    numElementsX            = size(vX, 1);
    numRowsP, numColsP      = size(mP);
    numRowsA, numColsA      = size(mA);

    ρ¹ = 1 / ρ;

    mL = mP + (σ * sparse(I, numElementsX, numElementsX)) + (ρ * (mA' * mA));

    for ii in 1:(numFactor - 1)
        mL = mP + (σ * sparse(I, numElementsX, numElementsX)) + (ρ * (mA' * mA));
    end

    for ii in 1:numIterations
        vT, _ = Krylov.cg_lanczos(mL, vX, atol = ϵSolver, itmax = numItrSolver);
        vX .= vT;
    end

end

function LinOpCg!(vX, mP, mA, vV, ρ = 1e6, σ = 1e-6;
    numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500)

    # Iterative Solver using IterativeSolver cg!()

    numElementsX            = size(vX, 1);
    numRowsP, numColsP      = size(mP);
    numRowsA, numColsA      = size(mA);

    ρ¹ = 1 / ρ;

    # Define linear operator that models u <- (P + AᵀA) w.
    # We need to allocate one temporary vector.
    vT = zeros(numRowsA); #<! Buffer
    opL = LinearOperator(Float64, numElements, numElements, true, true, (vU, vW, α, β) -> begin
    mul!(vT, mA, vW);
    mul!(vU, mA', vT);
    mul!(vU, mP, vW, one(Float64), ρ);
    vU .= vU .+ (σ .* vW);
end);

for ii in 1:(numFactor - 1)
    opL = LinearOperator(Float64, numElements, numElements, true, true, (vU, vW, α, β) -> begin
        mul!(vT, mA, vW);
        mul!(vU, mA', vT);
        mul!(vU, mP, vW, one(Float64), ρ);
        vU .= vU .+ (σ .* vW);
    end);
end

for ii in 1:numIterations
    vT, _ = Krylov.cg(opL, vX, atol = ϵSolver, itmax = numItrSolver);
    vX .= vT;
end

end

function KrylovTriCg!(vX, mP, mA, vV, ρ = 1e6, σ = 1e-6;
    numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500)

    # Iterative Solver using IterativeSolver cg!()

    numElementsX            = size(vX, 1);
    numRowsP, numColsP      = size(mP);
    numRowsA, numColsA      = size(mA);

    ρ¹ = 1 / ρ;

    vVV = rand(numRowsA);

    mPI = mP + (σ * sparse(I, numElementsX, numElementsX));
    # mM = cholesky(mPI);
    mM = opCholesky(mPI);

    mN = ρ * sparse(I, numRowsA, numRowsA);

    for ii in 1:(numFactor - 1)
        mN = ρ * sparse(I, numRowsA, numRowsA);
    end

    for ii in 1:numIterations
        vXX, vVV, _ = tricg(mA', vX, vVV; M = mM, N = mN, atol = ϵSolver, itmax = numItrSolver);
        vX .= vXX;
    end

end

function KrylovTriMr!(vX, mP, mA, vV, ρ = 1e6, σ = 1e-6;
    numFactor = 3, numIterations = 50, ϵSolver = 1e-6, numItrSolver = 500)

    # Iterative Solver using IterativeSolver cg!()

    numElementsX            = size(vX, 1);
    numRowsP, numColsP      = size(mP);
    numRowsA, numColsA      = size(mA);

    ρ¹ = 1 / ρ;

    vVV = rand(numRowsA);

    mPI = mP + (σ * sparse(I, numElementsX, numElementsX));
    # mM = cholesky(mPI);
    mM = opCholesky(mPI);

    mN = ρ * sparse(I, numRowsA, numRowsA);

    for ii in 1:(numFactor - 1)
        mN = ρ * sparse(I, numRowsA, numRowsA);
    end

    for ii in 1:numIterations
        vXX, vVV, _ = trimr(mA', vX, vVV; M = mM, N = mN, atol = ϵSolver, itmax = numItrSolver);
        vX .= vXX;
    end

end

max_stars_repo_name: RoyiAvital/QuadraticProgramSolver

max_issues_repo_name: RoyiAvital/QuadraticProgramSolver

max_forks_repo_name: RoyiAvital/QuadraticProgramSolver