text:
# ------------------------------------------------------------------------------------------
# # Strings
#
# Topics:
# 1. How to get a string
# 2. String interpolation
# 3. String concatenation
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## How to get a string
#
# Enclose your characters in " " or """ """!
# ------------------------------------------------------------------------------------------

s1 = "I am a string."

s2 = """I am also a string. """

# ------------------------------------------------------------------------------------------
# There are a couple functional differences between strings enclosed in single and triple
# quotes.

# One difference is that, in the latter case, you can use quotation marks within your
# string.
# ------------------------------------------------------------------------------------------

"Here, we get an "error" because it's ambiguous where this string ends "

"""Look, Mom, no "errors"!!! """

# ------------------------------------------------------------------------------------------
# Note that ' ' define a character, but NOT a string!
# ------------------------------------------------------------------------------------------

typeof('a')

'We will get an error here'

# ------------------------------------------------------------------------------------------
# ## String interpolation
#
# We can use the $ sign to insert existing variables into a string and to evaluate
# expressions within a string.

# Below is an example that contains some highly sensitive personal information.
# ------------------------------------------------------------------------------------------

name = "Jane"
num_fingers = 10
num_toes = 10

println("Hello, my name is $name.")
println("I have $num_fingers fingers and $num_toes toes.")

 println("That is $(num_fingers + num_toes) digits in all!!")

# ------------------------------------------------------------------------------------------
# ## String concatenation
#
# Below are three ways we can concatenate strings!

# The first way is to use the `string()` function.

# `string()` converts non-string inputs to strings.
# ------------------------------------------------------------------------------------------

s3 = "How many cats ";
s4 = "is too many cats?";
😺 = 10

string(s3, s4)

string("I don't know, but ", 😺, " is too few.")

# ------------------------------------------------------------------------------------------
# We can also use `*` for concatenation!
# ------------------------------------------------------------------------------------------

s3*s4

# ------------------------------------------------------------------------------------------
# ### Exercises
#
# #### 2.1
# Create a string that says "hi" 1000 times, first with `repeat` and then with the
# exponentiation operator, which can call `*` under the hood. Assign it the variable `hi`
# below.
# ------------------------------------------------------------------------------------------

@assert hi == "hihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihi"

# ------------------------------------------------------------------------------------------
# #### 2.2
# Declare two variables
#
# ```julia
# a = 3
# b = 4
# ```
# and use them to create two strings:
# ```julia
# "3 + 4"
# "7"
# ```
# and store the results in `c` and `d` respectively
# ------------------------------------------------------------------------------------------

@assert c == "3 + 4"
@assert d == "7"

# ------------------------------------------------------------------------------------------
# Please click on `Validate` on the top, once you are done with the exercises.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# Setup

run(`git clone https://github.com/mohamed82008/TopOpt.jl TopOpt`)
cd(() -> run(`git checkout cuda_error`), "TopOpt")
using Pkg
Pkg.activate("./TopOpt")
Pkg.instantiate()

# Error example

using Revise, TopOpt, CuArrays, LinearAlgebra
CuArrays.allowscalar(false)

function gpu_testcase(s)
	problem = PointLoadCantilever(Val{:Linear}, s, (1.0, 1.0), 1.0, 0.3, 1.0);
	solver = FEASolver(Displacement, CG, MatrixFree, problem);
	cusolver = cu(solver);
	cusolver.vars .= 1;
	cuop = TopOpt.buildoperator(cusolver)
	cux = CuArray(ones(size(cuop, 1)))
	cuy = similar(cux)
	mul!(cuy, cuop, cux)
end

function cpu_testcase(s)
	problem = PointLoadCantilever(Val{:Linear}, s, (1.0, 1.0), 1.0, 0.3, 1.0);
	solver = FEASolver(Displacement, CG, MatrixFree, problem);
	solver.vars .= 1;
	op = TopOpt.buildoperator(solver)
	x = ones(size(op, 1))
	y = similar(x)
	mul!(y, op, x)
end

s = (4, 4)
cpu_testcase(s) # works
gpu_testcase(s) # works

s = (10, 4)
cpu_testcase(s) # works
gpu_testcase(s)
# Some times works and other times errors with:
# ERROR: CUDA error: unspecified launch failure (code #719, ERROR_LAUNCH_FAILED)

s = (100, 10)
gpu_testcase(s)
# Works the first time and then fails the second time with:
# ERROR: CUDA error: an illegal memory access was encountered (code #700, ERROR_ILLEGAL_ADDRESS)

#=
The block-thread config is found using https://github.com/mohamed82008/TopOpt.jl/blob/cuda_error/src/GPUUtils/GPUUtils.jl#L101.

The number of blocks and threads can be overwritten in https://github.com/mohamed82008/TopOpt.jl/blob/cuda_error/src/GPUUtils/GPUUtils.jl#L95.
=#

ext: jl


text:
# MarkovMaps
function inv(m::MarkovMap)
  ~isa(domain(m),AbstractInterval) && error("Non-interval domains not supported")
  nbranches(m) > 1 && error("Cannot invert map with multiple branches")
  domain(m.branches[1]) != domain(m) && error("Cannot invert map that isn't defined everywhere")
  unsafe_inv(m)
end

function unsafe_inv(m::MarkovMap{D,R,B}) where {D,R} where B<:FwdExpandingBranch
  b = m.branches[1]
  b2 = RevExpandingBranch(b.f,b.dfdx,b.rangedomain,b.domain)
  MarkovMap([b2],m.rangedomain,m.domain)
end

function unsafe_inv(m::MarkovMap{D,R,B}) where {D,R} where B<:RevExpandingBranch
  b = m.branches[1]
  b2 = FwdExpandingBranch(b.v,b.dvdx,b.rangedomain,b.domain)
  MarkovMap([b2],m.rangedomain,m.domain)
end

# ComposedMarkovMaps
inv(m::ComposedMarkovMap) = ComposedMarkovMap([inv(mm) for mm in reverse(m.maps)],m.rangedomain,m.domain)

# CircleMaps
function inv(m::AbstractCircleMap)
  ncover(m) > 1 && error("Cannot invert map with covering index > 1")
  unsafe_inv(m)
end
unsafe_inv(m::FwdCircleMap) = RevCircleMap(m.f,m.dfdx,m.rangedomain,m.domain)
unsafe_inv(m::RevCircleMap) = FwdCircleMap(m.v,m.dvdx,m.rangedomain,m.domain)


text:
using MPI

function runmpi(file; ntasks = 1, localhost = false)
    localhostenv = try
        parse(Bool, get(ENV, "CLIMATEMACHINE_TEST_RUNMPI_LOCALHOST", "false"))
    catch
        false
    end
    # Force mpiexec to exec on the localhost node
    # TODO: MicrosoftMPI mpiexec has issues if mpiexec.exe is in a different
    # folder as the MPI script to run so ignore for now.
    if (localhost || localhostenv) && (MPI.MPI_LIBRARY != MPI.MicrosoftMPI)
        localhostonly = `-host localhost`
    else
        localhostonly = ``
    end
    # by default some mpi runtimes will
    # complain if more resources (processes)
    # are requested than available on the node
    if MPI.MPI_LIBRARY == MPI.OpenMPI
        oversubscribe = `--oversubscribe`
    else
        oversubscribe = ``
    end
    @info "Running MPI test..." file ntasks
    # Running this way prevents:
    #   Balance Law Solver | No tests
    # since external tests are not returned as passed/fail
    @time @test MPI.mpiexec() do cmd
        Base.run(
            `$cmd $localhostonly $oversubscribe -np $ntasks $(Base.julia_cmd()) --startup-file=no --project=$(Base.active_project()) $file`;
            wait = true,
        )
        true
    end
end


text:
function random_sampling(A,b,obj,noisyp,Atest,btest,testx,instance)
    I = 100 # I experiments
    J = 30 # J samples for each experiment
    TS = 100 # Number of test samples

    true_x  = zeros(Float64, I, n)
    noisy_x = zeros(Float64, J, n, I)
    b_train = zeros(Float64, I, c)
    A_train = zeros(Float64, c, n, I)

    # Generate training data
    for T = 1:I
        mu_train = zeros(Float64, II, JJ) # Generate random inputs
        D_train = zeros(Float64, JJ, TT)
        for i = 1:II
            for j = 1:JJ
                mu_train[i,j] = (rand(rng, 50:100)/100)*Mumax[i,j]
                for t = 1:TT
                    D_train[j, t] = (rand(rng, 90:110)/100)*Nominal_Demand[j]
                end
            end
        end
        A_train[:,:,T], b_train[T, :] = ModelParameters(II,JJ,TT,mu_train,D_train)
        if T == 1
            A_train[:,:,1] = deepcopy(A)
            b_train[1, :] = deepcopy(b)
        end
        true_x[T, 1:n] = SolveModel(A_train[:, :, T], obj, b_train[T, :], n)
    end

    t_instance_s = time()
    re_solve_counter = 0
    global cone = fill(zeros(c,n), 1)
    Pred_Error_Evolution = zeros(Float64, I, 10)
    hat_x = zeros(Float64, I, n)
    hat_z = zeros(Float64, I, c)

    for T = 1:I
        istight = []

        for j = 1:J
            for k = 1:n
                noisy_x[j, k, T] = true_x[T,k] + 3*randn(rng, Float64)
            end
        end

        hat_x[T, :], r, z = SolveP1(A_train[:,:,T], b_train[T,:], noisy_x[:,:,T], 1, J)
        istight, z = tightind(A_train[:,:,T], b_train[T,:], hat_x[T, :])
        hat_z[T, :] = z

        if T == 1
            for j = 1:length(istight)
                cone[T][j,:] = A_train[istight[j], :, T]/norm(A_train[istight[j], :, T])
            end
        else
            B = zeros(Float64, length(istight), n)
            for j = 1:length(istight)
                B[j,:] = A_train[istight[j], :, T]/norm(A_train[istight[j], :, T])
            end
            push!(cone, B)
        end

        status = SolveFP(cone, T)

        if status == :Infeasible
            re_solve_counter += 1
            z = []
            x, r, z, status = SolveP1m(A_train[:,:,1:T], b_train[1:T,:], noisy_x[:,:,1:T], T, J, instance, hat_x, hat_z)

            if status2 != :Optimal
                break
            end
            hat_x[1:T, :] = x
            hat_z[1:T, :] = z

            global cone = []
            global cone = fill(zeros(c,n), 1)

            for k = 1:T
                istight = []
                istight, z = tightind(A_train[:,:,k], b_train[k,:], hat_x[k, :])
                if k == 1
                    for j = 1:length(istight)
                        cone[k][j,:] = A_train[istight[j], :, k]/norm(A_train[istight[j], :, k])
                    end
                else
                    B = zeros(Float64, length(istight), n)
                    for j = 1:length(istight)
                        B[j,:] = A_train[istight[j], :, k]/norm(A_train[istight[j], :, k])
                    end
                    push!(cone, B)
                end
            end
            status = SolveFP(cone, T)
        end

        for estimate = 1:10
            hatc, loss = SolveP2(noisyp[estimate, :], cone, T)
            x_pred  = zeros(Float16, TS, n)

            for test = 1:TS
                x_pred[test, 1:n] = SolveModel(Atest[:, :, test], hatc, btest[test, :], n)
            end
            Pred_Error_Evolution[T, estimate] = infnorm(testx, x_pred)
        end
    end
    cd(string(root,"Results"))
    df = DataFrame(Pred_Error_Evolution)
    s = @sprintf "PredictionError_dim_%03.d_randomsamp_instance_%03.d.csv" dim instance
    CSV.write(s, df)

    println("Runtime ", time() - t_instance_s)
    println("Total no of times full problem was solved, ", re_solve_counter)
    println("-------------------------------------------------")
end

max_stars_repo_name: ddolab/InvLP

max_issues_repo_name: ddolab/InvLP

max_forks_repo_name: ddolab/InvLP


text: #Nothing here yet


text:
export filt_seis_chans!, get_seis_channels

"""
    get_seis_channels(S::GphysData)

Get an array with the channel numbers of all seismic data channels in `S`.

Assumes each ID in S ends with an alphanumeric three-digit channel code and that
ID fields are separated by periods: for example, two channels with IDs
"XX.YYY.00.EHZ" and "_.YHY" both have an instrument code of 'H'.

Channel codes less than two characters long (e.g. "Z" in "AA.BBB.CC.Z") are ignored.

SEED channel codes of seismic and seismoacoustic data (for which operations
like detrend! and taper! are sane) include D, G, H, J, L, M, N, P, Z.

See also: `inst_codes`
"""
function get_seis_channels(S::GphysData;
  chans::ChanSpec=Int64[])

  if chans == Int64[]
    chans = Int64.(collect(1:S.n))
  elseif typeof(chans) == UnitRange
    chans = Int64.(collect(chans))
  elseif typeof(chans) <: Integer
    chans = [Int64(chans)]
  end
  keep = falses(length(chans))

  @inbounds for (n,i) in enumerate(chans)
    id = S.id[i]
    L = length(id)
    for j = L:-1:1
      if id[j] == '.'
        j > L-2 && break
        if (id[j+2] in seis_inst_codes) && (length(S.x[i]) > 0)
          keep[n] = true
          break
        end
      end
    end
  end

  return chans[keep]
end

@doc """
    filt_seis_chans!(chans::ChanSpec, S::GphysData)
    filt_seis_chans(chans::ChanSpec, S::GphysData)

Filter a channel list `chans` to channels in `S` that contain seismic data.

Assumes each ID in S ends with an alphanumeric three-digit channel code and that
ID fields are separated by periods: for example, if two channels have IDs
"XX.YYY.00.EHZ" and "_.YHY", each has an instrument code of 'H'.

Channel codes less than two characters long (e.g. "Z" in "AA.BBB.CC.Z") are ignored.

SEED channel codes of seismic and seismoacoustic data (for which operations
like detrend! and taper! are sane) include D, G, H, J, L, M, N, P, Z.

See also: `inst_codes`
""" filt_seis_chans!
function filt_seis_chans!(chans::Array{Int64,1}, S::GphysData)

  @inbounds for n = length(chans):-1:1
    id = S.id[chans[n]]
    L = length(id)
    for j = L:-1:1
      if id[j] == '.'
        if (j > L-2)
          deleteat!(chans, n)
          break
        end
        if (id[j+2] in seis_inst_codes) == false
          deleteat!(chans, n)
        end
        break
      elseif j == 1
        deleteat!(chans, n)
      end
    end
  end
  return nothing
end


text:
using Convex
Random.seed!(2)

m = 40;  n = 50
A = randn(m, n); b = randn(m, 1)
x = Variable(n)
problem = minimize(sumsquares(A * x - b), [x >= 0])

ϵ = 1e-8
opt = 12.38418747141913
solver =  LineSearchWrapper(GAP(0.5, 1.0, 1.0, eps=ϵ, verbose=1, checki=1, max_iters=10000), lsinterval=100)
solve!(problem, solver)


text:
#=
Public interface for (incomplete) LDLt factorizations.
=#

using SparseArrays
using LinearAlgebra

include("ldlt.jl")

"""
    SimplicialLDLt{T} <: Factorization{T}

Simplicial, potentially incomplete LDLt factorization of a sparse symmetric
matrix. This is the return type of [`simplicial_ldlt`](@ref),
[`ildlt_lof`](@ref), [`ildlt_tol`](@ref).

The individual components of the factorization can be accessed via `getproperty`:

| Component | Description                                 |
|:---------:|:--------------------------------------------|
| `F.L`     | `L` (unit lower triangular) part of `LDLt`  |
| `F.D`     | `D` (diagonal) part of `LDLt`               |
| `F.Lt`    | `Lt` (unit upper triangular) part of `LDLt` |
| `F.d`     | diagonal values of `D` as a `Vector`        |
"""
struct SimplicialLDLt{T,C} <: LinearAlgebra.Factorization{T}
    colptr::Vector{Int}
    rowval::Vector{Int}
    nzval::Vector{T}
    conj::C
end

Base.size(F::SimplicialLDLt) = (length(F.colptr)-1).*(1,1)
Base.size(F::SimplicialLDLt, i::Integer) = size(F, i)

function Base.getproperty(F::SimplicialLDLt, d::Symbol)
    Fp,Fi,Fv,conj = getfield.(Ref(F), (:colptr,:rowval,:nzval,:conj))
    if d === :d
        return Fv[@view(Fp[1:end-1])]
    elseif d === :D
        return Diagonal(F.d)
    elseif d === :L
        return I+tril(SparseMatrixCSC(size(F)...,Fp,Fi,Fv),-1)
    elseif d === :Lt
        if conj === Base.conj
            return adjoint(F.L)
        elseif conj === identity
            return transpose(F.L)
        else
            return conj.(transpose(F.L))
        end
    else
        return getfield(F, d)
    end
end

"""
    ldlt_structure(A::SparseMatrixCSC) -> P::SparseMatrixCSC{Bool}

Compute the structurally nonzero entries of the LDLt factorization of `A`.
Only the lower-triangular part of `A` is accessed, and only the
lower-triangular part of `P` is returned.
"""
function ldlt_structure(A)
    @assert size(A,1) == size(A,2)
    @assert size(A,1) == length(A.colptr)-1
    Ap,Ai = A.colptr,A.rowval,A.nzval
    Fp,Fi = symbolic_ldlt(Ap,Ai)
    return SparseMatrixCSC(size(A)...,Fp,Fi,ones(Bool,length(Fi)))
end

"""
    level_of_fill(A::SparseMatrixCSC, c = size(A,2)) -> lof::SparseMatrixCSC{Int}

Compute the level-of-fill of each entry in the sparsity pattern of the
incomplete LDLt factorization of `A` where fill-in is restricted to entries
with `level_of_fill <= c`.

Only the lower-triangular part of `A` is accessed, and only the
lower-triangular part of `P` is returned.

Note that both structural zeros and entries with level-of-fill 0 will satisfy
`lof[i,j] = 0`, but the former will be structural zeros of `lof` while the
latter will be structural nonzeros. The sparsity pattern of `lof` can be
conveniently reduced using the `SparseArrays.fkeep!()` function.
"""
function level_of_fill(A::SparseMatrixCSC,c=size(A,2))
    @assert size(A,1) == size(A,2)
    @assert size(A,1) == length(A.colptr)-1
    Ap,Ai = A.colptr,A.rowval,A.nzval
    Fp,Fi,Fl = symbolic_ildlt_lof(Ap,Ai,c)
    return SparseMatrixCSC(size(A)...,Fp,Fi,Fl)
end

"""
    simplicial_ldlt(A::SparseMatrixCSC; conj=conj) -> F

Compute the LDLt factorization of a sparse symmetric matrix `A`. Only the
lower-triangular part of `A` is accessed.

If `conj = conj` (default), then this function computes the Hermitian
LDLt factorization `A = L*D*L'`. If `conj = identity`, then this function
computes the complex symmetric LDLt factorization `A = L*D*transpose(L)`.

The individual components of the factorization can be accessed via `getproperty`:

| Component | Description                                 |
|:---------:|:--------------------------------------------|
| `F.L`     | `L` (unit lower triangular) part of `LDLt`  |
| `F.D`     | `D` (diagonal) part of `LDLt`               |
| `F.Lt`    | `Lt` (unit upper triangular) part of `LDLt` |
| `F.d`     | diagonal values of `D` as a `Vector`        |
"""
function simplicial_ldlt(A::SparseMatrixCSC; conj=conj)
    @assert size(A,1) == size(A,2)
    @assert size(A,1) == length(A.colptr)-1
    Ap,Ai,Av = A.colptr,A.rowval,A.nzval
    Fp,Fi = symbolic_ldlt(Ap,Ai)
    Fv = numeric_ldlt(Ap,Ai,Av,Fp,Fi; conj=conj)
    return SimplicialLDLt(Fp,Fi,Fv,conj)
end

"""
    ildlt(A::SparseMatrixCSC, P::AbstractSparseMatrixCSC; conj=conj) -> F

Compute the incomplete LDLt factorization of a sparse symmetric matrix `A`
where fill-in is restricted to nonzero entries of `P`. Only the
lower-triangular part of `A` and `P` are accessed.

If `conj = conj` (default), then this function computes the Hermitian
LDLt factorization `A = L*D*L'`. If `conj = identity`, then this function
computes the complex symmetric LDLt factorization `A = L*D*transpose(L)`.

The individual components of the factorization can be accessed via `getproperty`:

| Component | Description                                 |
|:---------:|:--------------------------------------------|
| `F.L`     | `L` (unit lower triangular) part of `LDLt`  |
| `F.D`     | `D` (diagonal) part of `LDLt`               |
| `F.Lt`    | `Lt` (unit upper triangular) part of `LDLt` |
| `F.d`     | diagonal values of `D` as a `Vector`        |
"""
function ildlt(A::SparseMatrixCSC, P::SparseArrays.AbstractSparseMatrixCSC; conj=conj)
    @assert size(A,1) == size(A,2)
    @assert size(A,1) == length(A.colptr)-1
    Ap,Ai,Av = A.colptr,A.rowval,A.nzval
    Fp,Fi = P.colptr,P.rowval
    Fv = numeric_ldlt(Ap,Ai,Av,Fp,Fi; conj=conj)
    return SimplicialLDLt(Fp,Fi,Fv,conj)
end

"""
    ildlt_lof(A::SparseMatrixCSC, c; conj=conj) -> F

Compute the incomplete LDLt factorization of a sparse symmetric matrix `A`
where fill-in is restricted to entries with `level_of_fill <= c`. Only the
lower-triangular part of `A` is accessed.

If `conj = conj` (default), then this function computes the Hermitian
LDLt factorization `A = L*D*L'`. If `conj = identity`, then this function
computes the complex symmetric LDLt factorization `A = L*D*transpose(L)`.

The individual components of the factorization can be accessed via `getproperty`:

| Component | Description                                 |
|:---------:|:--------------------------------------------|
| `F.L`     | `L` (unit lower triangular) part of `LDLt`  |
| `F.D`     | `D` (diagonal) part of `LDLt`               |
| `F.Lt`    | `Lt` (unit upper triangular) part of `LDLt` |
| `F.d`     | diagonal values of `D` as a `Vector`        |
"""
function ildlt_lof(A::SparseMatrixCSC, c::Integer; conj=conj)
    @assert size(A,1) == size(A,2)
    @assert size(A,1) == length(A.colptr)-1
    Ap,Ai,Av = A.colptr,A.rowval,A.nzval
    Fp,Fi = symbolic_ildlt_lof(Ap,Ai,c)
    Fv = numeric_ldlt(Ap,Ai,Av,Fp,Fi; conj=conj)
    return SimplicialLDLt(Fp,Fi,Fv,conj)
end

"""
    ildlt_tol(A::SparseMatrixCSC,tol; conj=conj) -> F

Compute the incomplete LDLt factorization of a`A` where fill-in is restricted
to entries with `|F[i,j]| > tol`. Only the lower-triangular part of `A` is
accessed.

If `conj = conj` (default), then this function computes the Hermitian
LDLt factorization `A = L*D*L'`. If `conj = identity`, then this function
computes the complex symmetric LDLt factorization `A = L*D*transpose(L)`.

The individual components of the factorization can be accessed via `getproperty`:

| Component | Description                                 |
|:---------:|:--------------------------------------------|
| `F.L`     | `L` (unit lower triangular) part of `LDLt`  |
| `F.D`     | `D` (diagonal) part of `LDLt`               |
| `F.Lt`    | `Lt` (unit upper triangular) part of `LDLt` |
| `F.d`     | diagonal values of `D` as a `Vector`        |
"""
function ildlt_tol(A::SparseMatrixCSC, tol::Real; conj=conj)
    @assert size(A,1) == size(A,2)
    @assert size(A,1) == length(A.colptr)-1
    Ap,Ai,Av = A.colptr,A.rowval,A.nzval
    Fp,Fi,Fv = ildlt_tol(Ap,Ai,Av,tol; conj=conj)
    return SimplicialLDLt(Fp,Fi,Fv,conj)
end

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text:
using CSV
using DataFrames
using YAML

using SLiDE  # see src/SLiDE.jl

"""
# 3: Edit data
Standardize data read from csv file into DataFrame and manipulated using the
`edit_with` function from the SLiDE module.

This example will most-likely eventually be incorporated into some build stream method.
Relevant functions can be found in the associated files:
- load_from() - src/parse/load_data.jl
- edit_with() - src/parse/edit_data.jl
- read_file() - src/parse/load_data.jl
"""

DATA_DIR = joinpath("tests", "data")

y = read_file(joinpath(DATA_DIR, "test_datastream.yml"))
df = read_file(joinpath(DATA_DIR, "test_datastream.csv"))

# TERNARY OPERATORS are a type of CONDITIONAL EVALUATION and can be used to write
# single-line if-then-else statements:
#   <if> ? <then> : <else>
# && is a type of SHORT-CIRCUIT EVALUATION that can be used to write
# single-line if-then statements:
#   <if> && (<else>)
# See:
# https://discourse.julialang.org/t/style-question-ternary-operator-or-short-circuit-operator-or-if-end/34224/2
# https://docs.julialang.org/en/v1/manual/control-flow/#man-conditional-evaluation-1
# https://docs.julialang.org/en/v1/manual/control-flow/#Short-Circuit-Evaluation-1
"Drop"     in keys(y) && (df = edit_with(df, y["Drop"]))
"Rename"   in keys(y) && (df = edit_with(df, y["Rename"]))
"Group"    in keys(y) && (df = edit_with(df, y["Group"]))
"Match"    in keys(y) && (df = edit_with(df, y["Match"]))
"Melt"     in keys(y) && (df = edit_with(df, y["Melt"]))
"Add"      in keys(y) && (df = edit_with(df, y["Add"]))
"Map"      in keys(y) && (df = edit_with(df, y["Map"]))
"Replace"  in keys(y) && (df = edit_with(df, y["Replace"]))
"Drop"     in keys(y) && (df = edit_with(df, y["Drop"]))
"Operate"  in keys(y) && (df = edit_with(df, y["Operate"]))
"Describe" in keys(y) && (df = edit_with(df, y["Describe"], y["CSVInput"]))
"Order"    in keys(y) && (df = edit_with(df, y["Order"]))

show(df)

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
# Stolen from Flux as well

for n in [2, 20, 200, 2000]
    x = randn(Float32, n, 2000, 100)
    model = FourierLayer(n, n, 2000, 100, 16)
    println("CPU n=$n")
    run_benchmark(model, x, cuda=false)
    println("CUDA n=$n")
    run_benchmark(model, x, cuda=true)
end


text:
function codegen_lowlevel_variable_method(
  f::Function,
  statetype::Union{Symbol, Void}=nothing,
  states::Bool=true,
  returnfields::Vector{Symbol}=Symbol[],
  nkeys::Integer=0,
  vfarg::Bool=false
)
  body::Expr
  fargs::Vector
  rvalues::Expr

  if nkeys == 0
    fargs = [:(_state.value)]
  elseif nkeys > 0
    if vfarg
      fargs = [Expr(:ref, :Any, [:(_states[$i].value) for i in 1:nkeys]...)]
    else
      fargs = [:(_state.value), Expr(:ref, :Any, [:(_states[$i].value) for i in 1:nkeys]...)]
    end
  else
    error("nkeys must be non-negative, got $nkeys")
  end

  nr = length(returnfields)
  if nr == 0
    body = :($(f)($(fargs...)))
  elseif nr == 1
    rvalues = Expr(:., :_state, QuoteNode(returnfields[1]))
    body = :($rvalues = $(f)($(fargs...)))
  elseif nr > 1
    rvalues = Expr(:tuple, [Expr(:., :_state, QuoteNode(returnfields[i])) for i in 1:nr]...)
    body = :($rvalues = $(f)($(fargs...)))
  else
    error("Vector of return symbols must have one or more elements")
  end

  args = [statetype == nothing ? :_state : :(_state::$statetype)]
  if states
    push!(args, :(_states::VariableStateVector))
  end

  @gensym lowlevel_variable_method

  quote
    function $lowlevel_variable_method($(args...))
      $body
    end
  end
end


text:
using Test
using PETSc, MPI
using ForwardDiff, SparseArrays
MPI.Initialized() || MPI.Init()

@testset "DMStagCreate1d" begin

    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt    = PETSc.inttype(petsclib)

        # Create 1D DMStag
        dm = PETSc.DMStagCreate1d(
                petsclib,
                comm,
                PETSc.DM_BOUNDARY_PERIODIC,
                20,
                2,
                2,
                PETSc.DMSTAG_STENCIL_BOX,
                2)
        @test PETSc.DMStagGetBoundaryTypes(dm)==PETSc.DM_BOUNDARY_PERIODIC
        PETSc.destroy(dm)

        # Create 1D DMStag with array of local @ of points
        dm = PETSc.DMStagCreate1d(petsclib,comm,PETSc.DM_BOUNDARY_NONE,20,2,2,PETSc.DMSTAG_STENCIL_BOX,2,[20])

        # Test get size
        @test PETSc.DMStagGetGlobalSizes(dm) == 20
        @test PETSc.DMStagGetLocalSizes(dm) == 20

        # Test
        @test PETSc.gettype(dm) == "stag"
        @test PETSc.getdimension(dm) == 1

        # Info about ranks
        @test PETSc.DMStagGetIsFirstRank(dm) == (true,false,false)
        @test PETSc.DMStagGetIsLastRank(dm) == (true,false,false)

        # Boundary
        @test PETSc.DMStagGetBoundaryTypes(dm)==PETSc.DM_BOUNDARY_NONE

        # Corners
        corners         = PETSc.getcorners(dm)
        ghost_corners   = PETSc.getghostcorners(dm)

        @test corners.lower[1] == 1
        @test corners.upper[1] == 20
        @test corners.size[1]  == 20
        @test corners.extra[1] == 1

        @test ghost_corners.lower[1] == 1
        @test ghost_corners.upper[1] == 21
        @test ghost_corners.size[1]  == 21
        @test ghost_corners.extra[1]  == 0

        # DOF
        @test PETSc.DMStagGetDOF(dm) == (2,2)
        PETSc.destroy(dm)

        # Create new struct and pass keyword arguments
        dm_1D = PETSc.DMStagCreate1d(petsclib,comm,PETSc.DM_BOUNDARY_NONE,200,2,2; stag_grid_x=10);
        @test PETSc.DMStagGetGlobalSizes(dm_1D) == 10
        @test PETSc.DMStagGetEntriesPerElement(dm_1D)==4

        # Stencil width & type
        @test  PETSc.DMStagGetStencilWidth(dm_1D)==2
        @test  PETSc.DMStagGetBoundaryTypes(dm_1D) == PETSc.DM_BOUNDARY_NONE

        PETSc.destroy(dm_1D)

        # test ghosted array set using keywords
        dm_ghosted = PETSc.DMStagCreate1d(petsclib,comm,PETSc.DM_BOUNDARY_GHOSTED,200,2,2; stag_grid_x=10);
        @test  PETSc.DMStagGetStencilWidth(dm_ghosted)==2
        corners         = PETSc.getcorners(dm_ghosted)

        @test corners.size[1]==10       # keyword overrides the specified value
        @test  PETSc.DMStagGetBoundaryTypes(dm_ghosted) == PETSc.DM_BOUNDARY_GHOSTED

        ind = PETSc.DMStagGetIndices(dm_ghosted);
        @test ind.center.x[3] == 5

        # simple test to retrieve the KSP object
        # NOTE: need to implement a similar SNES routine
        ksp = PETSc.KSP(dm_ghosted)
        @test PETSc.gettype(ksp)=="gmres"

        PETSc.destroy(dm_ghosted)

        PETSc.finalize(petsclib)
    end
end

@testset "DMStagCreate2d" begin

    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt    = PETSc.inttype(petsclib)

        # Create 2D DMStag
        dm_2D = PETSc.DMStagCreate2d(petsclib, comm,
                PETSc.DM_BOUNDARY_NONE,
                PETSc.DM_BOUNDARY_NONE,
                20,
                21,
                1,
                1,
                1,
                1,
                1,
                PETSc.DMSTAG_STENCIL_BOX,
                2)

        @test PETSc.DMStagGetGlobalSizes(dm_2D) == (20,21)
        corners = PETSc.getcorners(dm_2D)
        @test corners.size  == [20,21,0]
        @test corners.extra == [1, 1, 0]
        @test corners.lower == [1, 1, 1]
        @test corners.upper == [20,21,0]

        PETSc.finalize(petsclib)
    end
end

@testset "DMStagCreate3d" begin

    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt    = PETSc.inttype(petsclib)

        # Create 3D DMStag
        dm_3D = PETSc.DMStagCreate3d(petsclib,comm,
                PETSc.DM_BOUNDARY_NONE,
                PETSc.DM_BOUNDARY_NONE,
                PETSc.DM_BOUNDARY_NONE,
                20,
                21,
                22,
                1,
                1,
                1,
                2,
                2,
                2,
                2,
                PETSc.DMSTAG_STENCIL_BOX,
                1,
                [],
                [],
                [])
        @test PETSc.DMStagGetGlobalSizes(dm_3D) == (20, 21, 22)

        # copy struct
        dmnew = PETSc.DMStagCreateCompatibleDMStag(dm_3D,1,1,2,2)
        @test PETSc.DMStagGetGlobalSizes(dmnew) == (20, 21, 22)

        @test PETSc.DMStagGetGlobalSizes(dm_3D) == (20,21,22)
        corners = PETSc.getcorners(dm_3D)
        @test corners.size  == [20,21,22]
        @test corners.extra == [1, 1, 1]
        @test corners.lower == [1, 1, 1]
        @test corners.upper == [20,21,22]

        PETSc.finalize(petsclib)
    end
end

@testset "DMStag Vectors and Coordinates" begin
    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt    = PETSc.inttype(petsclib)

        # Create 1D DMStag
        dm_1D = PETSc.DMStagCreate1d(petsclib,comm,PETSc.DM_BOUNDARY_NONE,200,2,2; stag_grid_x=10);
        @test PETSc.DMStagGetGlobalSizes(dm_1D) == 10

        # Set coordinates
        PETSc.setuniformcoordinates!(dm_1D, (0,), (10,))

        DMcoord = PETSc.getcoordinateDM(dm_1D)
        @test PETSc.gettype(DMcoord)=="stag"

        # Retrieve array with staggered coordinates
        coord_vec   = PETSc.getcoordinateslocal(dm_1D)
        X_coord     = PETSc.DMStagVecGetArray(DMcoord, coord_vec);
        @test  X_coord[1,2] == 0.5

        # Set coordinates using product (1D) arrays
        dm_1D = PETSc.DMStagCreate1d(petsclib,comm,PETSc.DM_BOUNDARY_NONE,200,2,2; stag_grid_x=10);
        PETSc.setuniformcoordinatesproduct!(dm_1D, (0,), (10,))

        # retrieve DM with coordinates
        DMcoord = PETSc.getcoordinateDM(dm_1D)
        @test PETSc.gettype(DMcoord)=="product"
        coord_vec   = PETSc.getcoordinateslocal(dm_1D)

        # Note: retrieving 1D coordinate vectors using the "product" type appears broken
        #  This is something to be looked at later
        #x_coord,y_coord,z_coord = PETSc.DMStagGetProductCoordinateArraysRead(dm_1D); # BROKEN

        PETSc.DMStagGetLocationSlot(dm_1D, PETSc.DMSTAG_RIGHT, 0) ==4
        @test PETSc.DMStagGetProductCoordinateLocationSlot(dm_1D, PETSc.DMSTAG_RIGHT) == 2

        global_vec      = PETSc.createglobalvector(dm_1D)
        local_vec       = PETSc.createlocalvector(dm_1D)

        # Fill everything with some data
        fill!(local_vec, mpisize)
        fill!(global_vec, mpisize)
        @test global_vec[3] == 1.0

        # Add the local values to the global values
        PETSc.update!(global_vec, local_vec, PETSc.ADD_VALUES)
        @test global_vec[3] == 2.0

        # PETSc.destroy(DMcoord);

        # Do 2D tests
        dm_2D = PETSc.DMStagCreate2d(petsclib, comm,
            PETSc.DM_BOUNDARY_NONE,
            PETSc.DM_BOUNDARY_NONE,
            3,
            4,
            1,
            1,
            1,
            1,
            1,
            PETSc.DMSTAG_STENCIL_BOX,
            2)

        PETSc.setuniformcoordinates!(dm_2D, (1,3), (10,11))
        coord_vec   = PETSc.getcoordinateslocal(dm_2D)

        # Retrieve array with staggered coordinates
        DMcoord_2D = PETSc.getcoordinateDM(dm_2D)
        X_coord_2D = PETSc.DMStagVecGetArray(DMcoord_2D, coord_vec);

        @test X_coord_2D[3,3,1] ≈ 7.0
        @test X_coord_2D[3,4,2] ≈ 9.0
        @test X_coord_2D[3,3,3] ≈ 8.5
        @test X_coord_2D[3,3,4] ≈ 7.0

        vec_test_2D     = PETSc.createlocalvector(dm_2D)
        X               = PETSc.DMStagVecGetArray(dm_2D,vec_test_2D);
        X[end,end,end] = 111;                   # modify 3D array @ some point and DOF
        @test vec_test_2D[end]==111.0           # verify that this modified the vector as well
        Base.finalize(X)                        # release from memory

        #test stencil locations
        pos1 = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_LEFT,1,0,0,1)
        @test pos1.c == 1
        pos2 = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_RIGHT,4,0,0,0)
        pos  = [pos1, pos2]
        @test pos2.loc == PETSc.DMSTAG_RIGHT
        @test pos2.i == 4

        # Retrieve value from stencil
        vec_test       = PETSc.createlocalvector(dm_1D)
        vec_test      .= 1:length(vec_test)                 # point wise copy of data to PetscVec
        val             = PETSc.DMStagVecGetValuesStencil(dm_1D, vec_test, pos1) # this gets a single value
        @test val ==6
        vals            = PETSc.DMStagVecGetValuesStencil(dm_1D, vec_test, 2, pos)     # this gets an array of values
        @test vals[1] == 6

        X_1D            = PETSc.DMStagVecGetArray(dm_1D,vec_test);
        @test X_1D[2,3] == 7.0

        # Set values using stencils
        vec_test_global = PETSc.createglobalvector(dm_1D)
        val1 = PetscScalar.([2222.2, 3.2]);
        PETSc.DMStagVecSetValuesStencil(dm_1D, vec_test_global,    pos1, val1[1],   PETSc.INSERT_VALUES)
        @test vec_test_global[6] ≈ 2222.2
        PETSc.DMStagVecSetValuesStencil(dm_1D, vec_test_global, 2, pos,  val1,      PETSc.INSERT_VALUES)
        @test vec_test_global[21] ≈ 3.2

        pos3 = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_LEFT,1,0,0,1)
        val = PETSc.DMStagVecGetValuesStencil(dm_1D, vec_test, 2, [pos3; pos3])
        @test val[2] == 6.0
        PETSc.destroy(dm_1D);

        PETSc.finalize(petsclib)
    end
end

@testset "DMStag create matrixes" begin
    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt    = PETSc.inttype(petsclib)

        dm_1D = PETSc.DMStagCreate1d(petsclib,comm,PETSc.DM_BOUNDARY_NONE,200,2,2; stag_grid_x=10);
        PETSc.setuniformcoordinatesproduct!(dm_1D, (0,), (10,))

        A = PETSc.creatematrix(dm_1D);  #
        @test size(A) == (42,42)
        PETSc.assembled(A)

        # set some values using normal indices:
        A[1,1]  = 1.0
        A[1,10] = 1.0

        pos1 = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_LEFT,1,0,0,1)
        pos2 = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_RIGHT,4,0,0,0)
        pos  = [pos1, pos2]
        val1 = PetscScalar.([2222.2, 3.2]);
        PETSc.DMStagMatSetValuesStencil(dm_1D, A, pos1, pos1, 11.1, PETSc.INSERT_VALUES)
        PETSc.DMStagMatSetValuesStencil(dm_1D, A, 1, [pos2], 2, pos, val1, PETSc.INSERT_VALUES)

        PETSc.assemble(A)
        @test A[1,10] == 1.0

        # Reads a value from the matrix, using the stencil structure
        @test PETSc.DMStagMatGetValuesStencil(dm_1D, A, pos1, pos1)== PetscScalar(11.1)
        @test PETSc.DMStagMatGetValuesStencil(dm_1D, A, 1, [pos2], 2, pos)==val1

        PETSc.destroy(dm_1D);

        dofCenter       =   1;
        dofEdge         =   1;
        dofVertex       =   0
        stencilWidth    =   1;
        dm_2D = PETSc.DMStagCreate2d(petsclib,comm,
                                        PETSc.DM_BOUNDARY_GHOSTED,
                                        PETSc.DM_BOUNDARY_GHOSTED,
                                        10,11,
                                        PETSc.PETSC_DECIDE,PETSc.PETSC_DECIDE,
                                        dofVertex,dofEdge,dofCenter,
                                        PETSc.DMSTAG_STENCIL_BOX,stencilWidth)

        vec_test_2D_global      =   PETSc.createglobalvector(dm_2D)
        vec_test_2D_local       =   PETSc.createlocalvector(dm_2D)

        corners                 =   PETSc.getcorners(dm_2D)
        ghost_corners           =   PETSc.getghostcorners(dm_2D)

        for ix=corners.lower[1]:corners.upper[1]
            for iy=corners.lower[2]:corners.upper[2]
                local dof
                # DOF at the center point
                dof     = 0;
                posA    = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_DOWN,ix,iy,0,dof)
                value   = PetscScalar(ix+10);
                PETSc.DMStagVecSetValuesStencil(dm_2D, vec_test_2D_global, posA, value, PETSc.INSERT_VALUES)

                dof     = 0;
                posB    = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_LEFT,ix,iy,0,dof)
                value   = PetscScalar(33);
                PETSc.DMStagVecSetValuesStencil(dm_2D, vec_test_2D_global, posB, value, PETSc.INSERT_VALUES)

                dof     = 0;
                posC    = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_ELEMENT,ix,iy,0,dof)
                value   = PetscScalar(44);
                PETSc.DMStagVecSetValuesStencil(dm_2D, vec_test_2D_global, posC, value, PETSc.INSERT_VALUES)

            end
        end
        PETSc.assemble(vec_test_2D_global) # assemble global vector

        # Add the global values to the local values
        PETSc.update!(vec_test_2D_local, vec_test_2D_global,PETSc.INSERT_VALUES)

        # retrieve value back from the local array and check that it agrees with global one
        dof     = 0;
        pos     = PETSc.DMStagStencil{PetscInt}(PETSc.DMSTAG_DOWN,2,2,0,dof)
        @test PETSc.DMStagVecGetValuesStencil(dm_2D, vec_test_2D_local, pos) == 12.0

        # Extract an array that holds all DOF's
        X2D_dofs  = PETSc.DMStagVecGetArray(dm_2D,vec_test_2D_local)           # extract arrays with all DOF (mostly for visualizing)
        @test X2D_dofs[4,4,1] ≈ PetscScalar(12.0)
        @test X2D_dofs[4,4,2] ≈ PetscScalar(33.0)
        @test X2D_dofs[4,4,3] ≈ PetscScalar(44.0)

        # Extract an array of a specific DOF (here a face velocity @ the left)
        Xarray = PETSc.DMStagGetGhostArrayLocationSlot(dm_2D,vec_test_2D_local, PETSc.DMSTAG_LEFT, 0)
        @test sum(X2D_dofs[:,:,2]-Xarray)==0                # check if the local array is identical to the full array

        Xarray .= 111.                                      # Set a value @ a specific location
        @test vec_test_2D_local[2] ≈ PetscScalar(111)       # verify that this is changed in the PETSc Vec

        PETSc.finalize(petsclib)
    end
end

# -----------------
# Example of DMStag & SNES with AD jacobian
@testset "DMStag: 1D SNES AD" begin

    # Tell AD that it can handle Complex as scalars
    ForwardDiff.can_dual(::Type{ComplexF64}) = true
    ForwardDiff.can_dual(::Type{ComplexF32}) = true

    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)

    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt    = PETSc.inttype(petsclib)

        # Define a struct that holds data we need in the local SNES routines below
        mutable struct Data_1{PetscScalar,PetscInt}
            dm
            x_l
            f_l
        end

        user_ctx = Data_1{PetscScalar,PetscInt}(nothing, nothing, nothing);  # holds data we need in the local

        function FormRes!(ptr_fx_g, ptr_x_g, user_ctx)

            # Note that in PETSc, ptr_x_g and ptr_fx_g are pointers to global vectors.
            # Copy global to local vectors that are stored in user_ctx
            PETSc.update!(user_ctx.x_l, ptr_x_g,   PETSc.INSERT_VALUES)
            PETSc.update!(user_ctx.f_l, ptr_fx_g,  PETSc.INSERT_VALUES)

            # Retrieve arrays from the local vectors
            ArrayLocal_x     =   PETSc.DMStagVecGetArrayRead(user_ctx.dm, user_ctx.x_l);  # array with all local x-data
            ArrayLocal_f     =   PETSc.DMStagVecGetArray(user_ctx.dm, user_ctx.f_l);      # array with all local residual

            # Compute local residual
            ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx)

            # Finalize local arrays
            Base.finalize(ArrayLocal_x)
            Base.finalize(ArrayLocal_f)

            # Copy local into global residual vector
            PETSc.update!(ptr_fx_g, user_ctx.f_l,   PETSc.INSERT_VALUES)

        end

        function ComputeLocalResidual(dm, ArrayLocal_x, ArrayLocal_f, user_ctx)
            # Compute the local residual. The vectors include ghost points

            T              =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x, PETSc.DMSTAG_LEFT,    0);
            fT             =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f, PETSc.DMSTAG_LEFT,    0);

            P              =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x, PETSc.DMSTAG_ELEMENT, 0);
            fP             =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f, PETSc.DMSTAG_ELEMENT, 0);

            # compute the FD stencil
            indices         =     PETSc.DMStagGetIndices(dm);      # indices of (center/element) points, not including ghost values.
            gc              =     PETSc.getghostcorners(user_ctx.dm);   # start and end of loop including ghost points
            c               =     PETSc.getcorners(user_ctx.dm);        # start and end of loop including ghost points

            nT             =     length(T);                             # array length
            dx             =     1.0/(c.size[1]-1);
            xp             =     (gc.lower[1]:gc.upper[1]).*dx;         # coordinates including ghost points (to define source term)
            F              =     6.0.*xp .+ (xp .+1.e-12).^6.0;         # define source term function

            # Nonlinear equation @ nodal points
            ind            =     indices.vertex.x;                         #  There is one more "vertex" point
            i              =     ind[2:end-1]
            fT[ind[1]]     =     T[ind[1]  ]-0.5;                       # left BC
            fT[ind[end]]   =     T[ind[end]]-2.0;                       # right BC
            fT[i]          =     (T[i .+ 1] - 2*T[i] + T[i .- 1])/dx^2  + T[i].*T[i] - F[i] # NL diffusion with source term

            # second, non-coupled, equation @ center points
            ind            =     indices.center.x;                             #  There is one more "vertex" point
            i              =     ind[2:end-1];
            fP[ind[1]]     =     P[ind[1]]-30.;                             # left BC
            fP[ind[end]]   =     P[ind[end]]-20.;                           # right BC
            fP[i]          =     (P[i .+ 1] - 2*P[i] + P[i .- 1])/dx^2      # steady state diffusion

        end

        function  ForwardDiff_res(x, user_ctx)

            f   = zero(x)               # vector of zeros, of same type as x (local vector)

            ArrayLocal_x     =   PETSc.DMStagVecGetArray(user_ctx.dm, x);        # array with all local x-data
            ArrayLocal_f     =   PETSc.DMStagVecGetArray(user_ctx.dm, f);        # array with all local residual

            ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx);
            # As the residual vector f is linked with ArrayLocal_f, we don't need to pass ArrayLocal_f back

            return f;
        end

        function FormJacobian!(ptr_x_g, J, P, user_ctx)

            # This requires several steps:
            #
            #   1) Extract local vector from global solution (x) vector
            #   2) Compute local jacobian from the residual routine (note that
            #       this routine requires julia vectors as input)

            # Extract the local vector
            PETSc.update!(user_ctx.x_l, ptr_x_g,  PETSc.INSERT_VALUES)
            x               =   PETSc.unsafe_localarray(PetscScalar, user_ctx.x_l.ptr;  write=false, read=true)
            f_Residual      =   (x -> ForwardDiff_res(x, user_ctx));        # pass additional arguments into the routine
            J_julia         =   ForwardDiff.jacobian(f_Residual,x);

            # Note: since x is the LOCAL vector, J_julia also ends up having the same size.
            ind             =   PETSc.LocalInGlobalIndices(user_ctx.dm);
            if PETSc.assembled(J) == false
                J           =   PETSc.MatSeqAIJ(sparse(J_julia[ind,ind]));
            else
                J           .=   sparse(J_julia[ind,ind]);
            end

           Base.finalize(x_g)
           return sparse(J_julia[ind,ind]), ind
        end

        # Main part

        # Construct a 1D test case for a coupled P-T diffusion solver, with 1 DOF @ the center & 1 DOF @ faces
        nx              =   21;
        user_ctx.dm     =   PETSc.DMStagCreate1d(petsclib,comm,
                                PETSc.DM_BOUNDARY_GHOSTED,
                                nx,
                                1,                              # DOF @ vertex
                                1,                              # DOF @ center
                                PETSc.DMSTAG_STENCIL_BOX,
                                1);                             # Stencil width

        x_g             =   PETSc.createglobalvector(user_ctx.dm)
        f_g             =   PETSc.createglobalvector(user_ctx.dm)
        user_ctx.x_l    =   PETSc.createlocalvector(user_ctx.dm)
        user_ctx.f_l    =   PETSc.createlocalvector(user_ctx.dm)

        PJ           =      PETSc.creatematrix(user_ctx.dm);                  # extract (global) matrix from DMStag
        J_julia, ind =      FormJacobian!(x_g, PJ, PJ, user_ctx)
        PJ           =      PETSc.MatSeqAIJ(J_julia)                # assemble non-zero structure

        S = PETSc.SNES{PetscScalar}(petsclib, comm;
                snes_rtol=1e-12,
                snes_view=false,
                snes_monitor=true,
                ksp_view=true,
                # pc_type="none",
                snes_monitor_true_residual=false,
                snes_converged_reason=false);
        S.user_ctx  =       user_ctx;

        PETSc.setfunction!(S, FormRes!, f_g)
        PETSc.setjacobian!(S, FormJacobian!, PJ, PJ)

        # Solve
        PETSc.solve!(x_g, S);
       # @show x_g

        # check
        @test x_g[4] ≈ 29.5
        @test x_g[11] ≈ 0.63797 rtol=1e-4

        PETSc.finalize(petsclib)

    end
end

@testset "DMStag: 2D SNES AD"  begin

    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)

    # Tell AD that it can handle Complex as scalars
    ForwardDiff.can_dual(::Type{ComplexF64}) = true
    ForwardDiff.can_dual(::Type{ComplexF32}) = true

    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt    = PETSc.inttype(petsclib)

        mutable struct Data_2D{PetscScalar,PetscInt}
            dm
            x_l
            f_l
        end
        user_ctx = Data_2D{PetscScalar,PetscInt}(nothing, nothing, nothing);  # holds data we need in the local

        function FormRes!(ptr_fx_g, ptr_x_g, user_ctx)
            # Note that in PETSc, ptr_x_g and ptr_fx_g are pointers to global vectors.

            # Copy global to local vectors
            PETSc.update!(user_ctx.x_l, ptr_x_g,  PETSc.INSERT_VALUES)
            PETSc.update!(user_ctx.f_l, ptr_fx_g, PETSc.INSERT_VALUES)

            # Retrieve arrays from the local vectors
            ArrayLocal_x     =   PETSc.DMStagVecGetArrayRead(user_ctx.dm, user_ctx.x_l);  # array with all local x-data
            ArrayLocal_f     =   PETSc.DMStagVecGetArray(user_ctx.dm, user_ctx.f_l);      # array with all local residual

            # Compute local residual
            ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx)

            # Finalize local arrays
            Base.finalize(ArrayLocal_x)
            Base.finalize(ArrayLocal_f)

            # Copy local into global residual vector
            PETSc.update!(ptr_fx_g, user_ctx.f_l, PETSc.INSERT_VALUES)

        end

        function  ForwardDiff_res(x, user_ctx)
            f   = zero(x)               # vector of zeros, of same type as x (local vector)

            ArrayLocal_x     =   PETSc.DMStagVecGetArray(user_ctx.dm, x);        # array with all ocal x-data
            ArrayLocal_f     =   PETSc.DMStagVecGetArray(user_ctx.dm, f);        # array with all ocal residual

            ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx);

            # As the residual vector f is linked with ArrayLocal_f, we don't need to pass ArrayLocal_f back

            return f;
        end

        function ComputeLocalResidual(dm, ArrayLocal_x, ArrayLocal_f, user_ctx)
            # Compute the local residual. The vectors include ghost points

            # Important! Make sure you retrieve the values from the correct locations. In this example we have a
            T              =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x, PETSc.DMSTAG_ELEMENT, 0);
            fT             =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f, PETSc.DMSTAG_ELEMENT, 0);

            # compute the FD stencil
            indices         =     PETSc.DMStagGetIndices(dm);      # indices of (center/element) points, not including ghost values.

            sz             =     size(user_ctx.dm);                                 # array length
            dx             =     1.0/(sz[1]-1);
            dz             =     1.0/(sz[2]-1);

            # set ghost points for BC'S
            bnd            =    PETSc.DMStagGetBoundaryTypes(user_ctx.dm)
            if bnd[1] == PETSc.DM_BOUNDARY_GHOSTED
                T[1,:]     =    T[2,:];        # zero flux; dT/dx=0
                T[end,:]   =    T[end-1,:];    # zero flux
            end

            # Diffusion @ center points
            indx        = indices.center.x;
            indz        = indices.center.y;

            ix          =     indx[1:end]                             # use ghost points in x  (required GHOSTED x-boundary)
            iz          =     indz[2:end-1]                           # center points

            # upper and lower BC (including corners)
            fT[indx,indz[1]]   =    T[indx,indz[1]] .- 0.5;                             # bottom BC
            fT[indx,indz[end]] =    T[indx,indz[end]] .- 2.0;                             # top BC

            fT[ix,iz]       =    (T[ix .+ 1,iz] - 2*T[ix,iz] + T[ix .- 1,iz])/dx^2   +
                                 (T[ix,iz .+ 1] - 2*T[ix,iz] + T[ix,iz .- 1])/dz^2

        end

        function FormJacobian!(ptr_x_g, J, P, user_ctx)
            # This requires several steps:
            #
            #   1) Extract local vector from global solution (x) vector
            #   2) Compute local jacobian from the residual routine (note that
            #       this routine requires julia vectors as input)

            # Extract the local vector
            PETSc.update!(user_ctx.x_l, ptr_x_g,  PETSc.INSERT_VALUES)
            x               =   PETSc.unsafe_localarray(PetscScalar, user_ctx.x_l.ptr;  write=false, read=true)

            f_Residual      =   (x -> ForwardDiff_res(x, user_ctx));        # pass additional rguments into the routine
            J_julia         =   ForwardDiff.jacobian(f_Residual,x);

            # Note: since x is the LOCAL vector, J_julia also ends up having the same size.
            ind             =   PETSc.LocalInGlobalIndices(user_ctx.dm);
            if PETSc.assembled(J) == false
                J           =   PETSc.MatSeqAIJ(sparse(J_julia[ind,ind]));
            else
                J           .=   sparse(J_julia[ind,ind]);
            end

            return sparse(J_julia[ind,ind]), ind, sparse(J_julia)
        end

        # Main routine starts here ----

        dofVertex   =   0
        dofEdge     =   0
        dofCenter   =   1
        nx,nz       =   6,25
        user_ctx.dm =   PETSc.DMStagCreate2d(petsclib,comm,
                    PETSc.DM_BOUNDARY_GHOSTED,
                    PETSc.DM_BOUNDARY_NONE,
                    nx,
                    nz,
                    1,
                    1,
                    dofVertex,
                    dofEdge,
                    dofCenter,
                    PETSc.DMSTAG_STENCIL_BOX,
                    1)
        PJ           =      PETSc.creatematrix(user_ctx.dm);                  # extract global) matrix from DMStag
        x_g             =   PETSc.createglobalvector(user_ctx.dm)
        f_g             =   PETSc.createglobalvector(user_ctx.dm)
        user_ctx.x_l    =   PETSc.createlocalvector(user_ctx.dm)
        user_ctx.f_l    =   PETSc.createlocalvector(user_ctx.dm)

        S = PETSc.SNES{PetscScalar}(petsclib, comm;
                snes_rtol=1e-12,
                snes_monitor=true,
                pc_type="none",
                snes_monitor_true_residual=true,
                snes_converged_reason=false);
        S.user_ctx  =       user_ctx;

        J_julia, ind, J_full =      FormJacobian!(x_g, PJ, PJ, user_ctx)
        PJ           =      PETSc.MatSeqAIJ(J_julia)                # assemble non-zero structure

        PETSc.setfunction!(S, FormRes!, f_g)
        PETSc.setjacobian!(S, FormJacobian!, PJ, PJ)

        # Solve 2D system
        sol = PETSc.solve!(x_g, S);

        PETSc.update!(user_ctx.x_l,sol, PETSc.INSERT_VALUES);   # copy global solution -> local vector
        T2d =   PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,user_ctx.x_l, PETSc.DMSTAG_ELEMENT,    0);

        @test T2d[5,5] ≈ 0.75 rtol=1e-3
        #
        # -----------------

        PETSc.finalize(petsclib)

    end
end


text:
###
# UDY Oracle Tests
###
#doesn't use any bounds
facts("portTest UDY") do
	srand(8675309); data = randn(500, 2)
	w = UDYOracle(data, .1, .2)
	portTest(w, -2.450295808139833, [0.5372084714704519, 0.4627915285295481], TOL=1e-6,
			unc_lower=[-10, -10], unc_upper=[10, 10])
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ConnectionMonitor <: SwaggerModel
    location::Any # spec type: Union{ Nothing, String } # spec name: location
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags
    properties::Any # spec type: Union{ Nothing, ConnectionMonitorParameters } # spec name: properties

    function ConnectionMonitor(;location=nothing, tags=nothing, properties=nothing)
        o = new()
        validate_property(ConnectionMonitor, Symbol("location"), location)
        setfield!(o, Symbol("location"), location)
        validate_property(ConnectionMonitor, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        validate_property(ConnectionMonitor, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type ConnectionMonitor

const _property_map_ConnectionMonitor = Dict{Symbol,Symbol}(Symbol("location")=>Symbol("location"), Symbol("tags")=>Symbol("tags"), Symbol("properties")=>Symbol("properties"))
const _property_types_ConnectionMonitor = Dict{Symbol,String}(Symbol("location")=>"String", Symbol("tags")=>"Dict{String, String}", Symbol("properties")=>"ConnectionMonitorParameters")
Base.propertynames(::Type{ ConnectionMonitor }) = collect(keys(_property_map_ConnectionMonitor))
Swagger.property_type(::Type{ ConnectionMonitor }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitor[name]))}
Swagger.field_name(::Type{ ConnectionMonitor }, property_name::Symbol) =  _property_map_ConnectionMonitor[property_name]

function check_required(o::ConnectionMonitor)
    (getproperty(o, Symbol("properties")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ConnectionMonitor }, name::Symbol, val)
end


text:
cfg = get_config("test.yaml")

function test_individual(ind::NEATInd)
    inputs = rand(cfg["n_in"])
    set_inputs(ind, inputs)
    @test ind.neurons[1].output == inputs[1]

    outputs = get_outputs(ind)
    @test length(outputs) == cfg["n_out"]
    @test all(outputs .== 0.0)

    inputs = zeros(cfg["n_in"])
    outputs = process(ind, inputs)
    @test length(outputs) == cfg["n_out"]
    @test any(outputs .!= 0.0)
    @test all(outputs .>= 0.0)
    @test all(outputs .<= 1.0)

    inputs = 10 .* rand(cfg["n_in"])
    outputs = process(ind, inputs)
    @test any(outputs .!= 0.0)
    @test all(outputs .>= 0.0)
    @test all(outputs .<= 1.0)
end

@testset "Processing individuals" begin
    ind = NEATInd(cfg)
    test_individual(ind)

    # TODO: add after mutation is done
    # ind2 = mutate_add_node(cfg, ind)
    # test_individual(ind2)

    # ind3 = mutate_add_link(cfg, ind)
    # test_individual(ind3)
end


text:
using LaTeXStrings
using Test

using MathTeXEngine

import MathTeXEngine: manual_texexpr, texparse
import MathTeXEngine: TeXParseError

import MathTeXEngine: tex_layout, generate_tex_elements
import MathTeXEngine: Space, TeXElement
import MathTeXEngine: _default_fonts, FontFamily, load_font
import MathTeXEngine: inkheight, inkwidth

include("parser.jl")
include("fonts.jl")
include("layout.jl")


text:
mutable struct RODEProblem{uType,tType,isinplace,P,NP,F,C,ND} <: AbstractRODEProblem{uType,tType,isinplace,ND}
  f::F
  u0::uType
  tspan::tType
  p::P
  noise::NP
  callback::C
  rand_prototype::ND
  seed::UInt64
  @add_kwonly function RODEProblem{iip}(f::RODEFunction{iip},u0,tspan,p=nothing;
                       rand_prototype = nothing,
                       noise= nothing, seed = UInt64(0),
                       callback=nothing) where {iip}
  _tspan = promote_tspan(tspan)
  new{typeof(u0),typeof(_tspan),
              isinplace(f),typeof(p),
              typeof(noise),typeof(f),typeof(callback),
              typeof(rand_prototype)}(
              f,u0,_tspan,p,noise,callback,
              rand_prototype,seed)
  end
  function RODEProblem{iip}(f,u0,tspan,p=nothing;kwargs...) where {iip}
    RODEProblem(convert(RODEFunction{iip},f),u0,tspan,p;kwargs...)
  end
end

function RODEProblem(f::RODEFunction,u0,tspan,p=nothing;kwargs...)
  RODEProblem{isinplace(f)}(f,u0,tspan,p;kwargs...)
end

function RODEProblem(f,u0,tspan,p=nothing;kwargs...)
  RODEProblem(convert(RODEFunction,f),u0,tspan,p;kwargs...)
end


text:
module BetaVQE
using LinearAlgebra
using Yao, Yao.EasyBuild, Yao.BitBasis

using Zygote
import Optimisers
import ChainRulesCore: @non_differentiable, NoTangent, Tangent, rrule
using StatsBase
include("VAN/VAN.jl")
using .VAN

export qaoa_circuit, tns_circuit
export free_energy, free_energy_local, energy, entropy
export TFIM, hamiltonian
export train, exact

include("hamiltonian.jl")
include("circuit.jl")
include("free_energy.jl")
include("zygote_patch.jl")

# exact value
function exact(β, h)
    H = mat(h)
    w, _ = eigen(Matrix(H))
    Z = sum(e->exp(-β*e), w)
    F = -log(Z)/β
    E = sum(e->exp(-β*e)*e, w)/Z
    E2 = sum(e->exp(-β*e)*e*e, w)/Z
    γ =  sum(e->exp(-2*β*e), w)/Z^2
    S = (E-F)*β
    Cv = (E2-E^2)*β^2
    F, E, S, Cv, γ
end

function loss(β, H, sampler, circuit, nbatch::Int)
    samples = gen_samples(sampler, nbatch)
    free_energy(β, H, sampler, circuit, samples)
end

function train(β::Real, H::AbstractBlock, sampler::AbstractSampler, circuit::AbstractBlock, logfile=nothing;
                optimizer=Optimisers.ADAM(0.1), nbatch::Int=1000, niter::Int=100)
    ϕ = model_parameters(sampler)
    θ = parameters(circuit)
    opt = Optimisers.setup(optimizer, (ϕ=ϕ,θ=θ))
    for i = 1:niter
        _, _, gϕ, gθ, _ = Zygote.gradient(loss, β, H, sampler, circuit, nbatch)
        Optimisers.update!(opt, (ϕ=ϕ, θ=θ), (ϕ=collect_gradients(sampler, gϕ), θ=gθ))
        model_dispatch!(sampler, ϕ)
        dispatch!(circuit, θ)
        message = "$i $(loss(β, H, sampler, circuit, nbatch))"
        println(message)
        flush(stdout)
        if !(logfile === nothing)
            write(logfile, message*"\n")
            flush(logfile)
        end
    end
    return ϕ, θ
end

# collect parameters into a tuple
collect_gradients(model::AutoRegressiveModel, g) = (g.W..., g.b...)
collect_gradients(model::PSAModel, g) = (g.w, )

end # module


text:
abstract type AbstractEpiModel end

mutable struct SIR <: AbstractEpiModel
    tspan::Tuple
    parameter:: NamedTuple
    saveat
end


text:
__precompile__()

module MujocoSim

import GLFW
using MuJoCo
using StaticArrays

##################################################### globals
const fontscale = mj.FONTSCALE_200 # can be 100, 150, 200
const maxgeom       = 5000 # preallocated geom array in mjvScene
const syncmisalign  = 0.1  # maximum time mis-alignment before re-sync
const refreshfactor = 0.7  # fraction of refresh available for simulation

mutable struct mjSim
   # visual interaction controls
   lastx::Float64
   lasty::Float64
   button_left::Bool
   button_middle::Bool
   button_right::Bool

   lastbutton::GLFW.MouseButton
   lastclicktm::Float64

   refreshrate::Int

   showhelp::Int
   showoption::Bool
   showdepth::Bool
   showfullscreen::Bool
   showsensor::Bool
   slowmotion::Bool

   showinfo::Bool
   paused::Bool
   keyreset::Int
   record::Any
   vidbuff::Vector{UInt8}

   framecount::Float64
   framenum::Int
   lastframenum::Int

   # MuJoCo things
   scn::Ref{mjvScene}
   cam::Ref{mjvCamera}
   vopt::Ref{mjvOption}
   pert::Ref{mjvPerturb}
   con::Ref{mjrContext}
   figsensor::Ref{mjvFigure}
   m::jlModel
   d::jlData

   # GLFW handle
   window::GLFW.Window
   vmode::GLFW.VidMode

   #uistate::mjuiState
   #ui0::mjUI
   #ui1::mjUI

   function mjSim(m::jlModel, d::jlData, name::String;
                  width=0, height=0)
      vmode = GLFW.GetVideoMode(GLFW.GetPrimaryMonitor())
      w = width > 0 ? width : floor(2*vmode.width / 3)
      h = height > 0 ? height : floor(2*vmode.height / 3)
      new(0.0, 0.0, false, false, false, GLFW.MOUSE_BUTTON_1, 0.0,
          vmode.refreshrate,
          0, false, false, false, false, false, true, true, 0,
          nothing, #open("/tmp/test.bin","w"),
          Vector{UInt8}(undef, 5120 * 2880 * 3),
      0.0, 0, 0,
      Ref(mjvScene()),
      Ref(mjvCamera()),
      Ref(mjvOption()),
      Ref(mjvPerturb()),
      Ref(mjrContext()),
      Ref(mjvFigure()),
      m, d,
      GLFW.CreateWindow(w, h, "Simulate"),
      vmode
     )
   end
end

export mjSim

const keycmds = Dict{GLFW.Key, Function}(GLFW.KEY_F1=>(s)->begin  # help
                                            s.showhelp += 1
                                            if s.showhelp > 2 s.showhelp = 0 end
                                         end,
                                         GLFW.KEY_F2=>(s)->begin  # option
                                            s.showoption = !s.showoption;
                                         end,
                                         GLFW.KEY_F3=>(s)->begin  # info
                                            s.showinfo = !s.showinfo;
                                         end,
                                         GLFW.KEY_F4=>(s)->begin  # depth
                                            s.showdepth = !s.showdepth;
                                         end,
                                         GLFW.KEY_F5=>(s)->begin  # toggle full screen
                                            s.showfullscreen = !s.showfullscreen;
                                            s.showfullscreen ? GLFW.MaximizeWindow(s.window) : GLFW.RestoreWindow(s.window)
                                         end,
                                         #GLFW.KEY_F6=>(s)->begin  # stereo
                                         #   s.stereo = s.scn.stereo == mj.mjSTEREO_NONE ? mjSTEREO_QUADBUFFERED : mj.mjSTEREO_NONE
                                         #   s.scn[].stereo
                                         #end,
                                         GLFW.KEY_F7=>(s)->begin  # sensor figure
                                            s.showsensor = !s.showsensor;
                                         end,
                                         GLFW.KEY_F8=>(s)->begin  # profiler
                                            s.showprofiler = !s.showprofiler;
                                         end,
                                         GLFW.KEY_ENTER=>(s)->begin  # slow motion
                                            s.slowmotion = !s.slowmotion;
                                            s.slowmotion ? println("Slow Motion Mode!") : println("Normal Speed Mode!")
                                         end,
                                         GLFW.KEY_SPACE=>(s)->begin  # pause
                                            s.paused = !s.paused
                                            s.paused ? println("Paused") : println("Running")
                                         end,
                                         GLFW.KEY_V=>(s)->begin
                                            if s.record == nothing
                                               println("Saving video to /tmp/video.mp4")
                                               #wi, hi = GLFW.GetFramebufferSize(w)
                                               s.record = open(`ffmpeg -y
                                                               -f rawvideo -pixel_format rgb24
                                                               -video_size 1280x800 -framerate 60
                                                               -i pipe:0
                                                               -preset fast -threads 0
                                                               -vf "vflip" /tmp/video.mp4`, "w")
                                            else
                                               println("Closing /tmp/video.mp4")
                                               close(s.record)
                                               s.record = nothing
                                            end
                                         end,
                                         GLFW.KEY_PAGE_UP=>(s)->begin    # previous keyreset
                                            s.keyreset = min(s.m.m[].nkey - 1, s.keyreset + 1)
                                         end,
                                         GLFW.KEY_PAGE_DOWN=>(s)->begin  # next keyreset
                                            s.keyreset = max(-1, s.keyreset - 1)
                                         end,
                                         # continue with reset
                                         GLFW.KEY_BACKSPACE=>(s)->begin  # reset
                                            mj_resetData(s.m.m, s.d.d)
                                            if s.keyreset >= 0 && s.keyreset < s.m.m[].nkey
                                               s.d[].time = s.m.key_time[s.keyreset+1]
                                               s.d.qpos[:] = s.m.key_qpos[:,s.keyreset+1]
                                               s.d.qvel[:] = s.m.key_qvel[:,s.keyreset+1]
                                               s.d.act[:]  = s.m.key_act[:,s.keyreset+1]
                                            end
                                            mj_forward(s.m, s.d)
                                            #profilerupdate()
                                            sensorupdate(s)
                                         end,
                                         GLFW.KEY_RIGHT=>(s)->begin  # step forward
                                            if s.paused
                                               mj_step(s.m, s.d)
                                               #profilerupdate()
                                               sensorupdate(s)
                                            end
                                         end,
                                         GLFW.KEY_LEFT=>(s)->begin  # step back
                                            if s.paused
                                               dt = s.m.m[].opt.timestep
                                               s.m.m[].opt.timestep = -dt
                                               #cleartimers(s.d);
                                               mj_step(s.m, s.d);
                                               s.m.m[].opt.timestep = dt
                                               #profilerupdate()
                                               sensorupdate(s)
                                            end
                                         end,
                                         GLFW.KEY_DOWN=>(s)->begin  # step forward 100
                                            if s.paused
                                               #cleartimers(d);
                                               for n=1:100 mj_step(s.m, s.d) end
                                               #profilerupdate();
                                               sensorupdate(s)
                                            end
                                         end,
                                         GLFW.KEY_UP=>(s)->begin  # step back 100
                                            if s.paused
                                               dt = s.m.m[].opt.timestep
                                               s.m.m[].opt.timestep = -dt
                                               #cleartimers(d)
                                               for n=1:100 mj_step(s.m, s.d) end
                                               s.m.m[].opt.timestep = dt
                                               #profilerupdate();
                                               sensorupdate(s)
                                            end
                                         end,
                                         GLFW.KEY_ESCAPE=>(s)->begin  # free camera
                                            s.cam[]._type = Int(mj.CAMERA_FREE)
                                         end,
                                         GLFW.KEY_EQUAL=>(s)->begin  # bigger font
                                            if fontscale < 200
                                               fontscale += 50
                                               mjr_makeContext(s.m.m, s.con, fontscale)
                                            end
                                         end,
                                         GLFW.KEY_MINUS=>(s)->begin  # smaller font
                                            if fontscale > 100
                                               fontscale -= 50;
                                               mjr_makeContext(s.m.m, s.con, fontscale);
                                            end
                                         end,
                                         GLFW.KEY_LEFT_BRACKET=>(s)->begin  # '[' previous fixed camera or free
                                            fixedcam = s.cam._type
                                            if s.m.m[].ncam > 0 && fixedcam == Int(mj.CAMERA_FIXED)
                                               fixedcamid = s.cam.fixedcamid
                                               if (fixedcamid  > 0)
                                                  s.cam[].fixedcamid = fixedcamid-1
                                               else
                                                  s.cam[]._type = Int(mj.CAMERA_FREE)
                                               end
                                            end
                                         end,
                                         GLFW.KEY_RIGHT_BRACKET=>(s)->begin  # ']' next fixed camera
                                            if s.m.m[].ncam > 0
                                               fixedcam = s.cam.type
                                                  fixedcamid = s.cam.fixedcamid
                                                  if fixedcam != Int(mj.CAMERA_FIXED)
                                                     s.cam[]._type = Int(mj.CAMERA_FIXED)
                                                  elseif fixedcamid < s.m.m[].ncam - 1
                                                     s.cam[].fixedcamid = fixedcamid+1
                                                  end
                                               end
                                            end,
                                            GLFW.KEY_SEMICOLON=>(s)->begin  # cycle over frame rendering modes
                                               frame = s.vopt.frame
                                               s.vopt[].frame = max(0, frame - 1)
                                            end,
                                            GLFW.KEY_APOSTROPHE=>(s)->begin  # cycle over frame rendering modes
                                               frame = s.vopt.frame
                                               s.vopt[].frame = min(Int(mj.NFRAME)-1, frame+1)
                                            end,
                                            GLFW.KEY_PERIOD=>(s)->begin  # cycle over label rendering modes
                                               label = s.vopt.label
                                               s.vopt[].label = max(0, label-1)
                                            end,
                                            GLFW.KEY_SLASH=>(s)->begin  # cycle over label rendering modes
                                               label = s.vopt.label
                                               s.vopt[].label = min(Int(mj.NLABEL)-1, label+1)
                                            end)

##################################################### functions
function alignscale(s::mjSim)
   s.cam[].lookat = s.m.m[].stat.center
   s.cam[].distance = 1.5*s.m.m[].stat.extent

   # set to free camera
   s.cam[]._type = Cint(mj.CAMERA_FREE)
end

function str2vec(s::String, len::Int)
   str = zeros(UInt8, len)
   str[1:length(s)] = codeunits(s)
   return str
end

# init sensor figure
function sensorinit(s::mjSim)
   # set figure to default
   mjv_defaultFigure(s.figsensor)

   # set flags
   s.figsensor[].flg_extend = Cint(1)
   s.figsensor[].flg_barplot = Cint(1)

   s.figsensor[].title = str2vec("Sensor data", length(s.figsensor[].title))

   # y-tick nubmer format
   s.figsensor[].yformat = str2vec("%.0f", length(s.figsensor[].yformat))

   # grid size
   s.figsensor[].gridsize = [2, 3]

   # minimum range
   s.figsensor[].range = [[0 1],[-1 1]]
end

# update sensor figure
function sensorupdate(s::mjSim)
   maxline = 10

   for i=1:maxline # clear linepnt
      mj.set(s.figsensor, :linepnt, Cint(0), i)
   end

   lineid = 1 # start with line 0
   m = s.m
   d = s.d

   # loop over sensors
   for n=1:m.m[].nsensor
      # go to next line if type is different
      if (n > 1 && m.sensor_type[n] != m.sensor_type[n - 1])
         lineid = min(lineid+1, maxline)
      end

      # get info about this sensor
      cutoff = m.sensor_cutoff[n] > 0 ? m.sensor_cutoff[n] : 1.0
      adr = m.sensor_adr[n]
      dim = m.sensor_dim[n]

      # data pointer in line
      p = mj.get(s.figsensor, :linepnt, lineid)

      # fill in data for this sensor
      for i=0:(dim-1)
         # check size
         if ((p + 2i) >= Int(mj.MAXLINEPNT) / 2) break end

         x1 = 2p + 4i + 1
         x2 = 2p + 4i + 3
         mj.set(s.figsensor, :linedata, adr+i, lineid, x1)
         mj.set(s.figsensor, :linedata, adr+i, lineid, x2)

         y1 = 2p + 4i + 2
         y2 = 2p + 4i + 4
         se = d.sensordata[adr+i+1]/cutoff
         mj.set(s.figsensor, :linedata,  0, lineid, y1)
         mj.set(s.figsensor, :linedata, se, lineid, y2)
      end

      # update linepnt
      mj.set(s.figsensor, :linepnt,
             min(Int(mj.MAXLINEPNT)-1, p+2dim),
             lineid)
   end
end

# show sensor figure
function sensorshow(s::mjSim, rect::mjrRect)
   # render figure on the right
   viewport = mjrRect(rect.width - rect.width / 4,
                      rect.bottom,
                      rect.width / 4,
                      rect.height / 3)
   mjr_figure(viewport, s.figsensor, s.con)
end

##################################################### callbacks

function mykeyboard(s::mjSim, window::GLFW.Window,
                    key::GLFW.Key, scancode::Int32, act::GLFW.Action, mods::Int32)
   # do not act on release
   if act == GLFW.RELEASE return end

   try
      keycmds[key](s) # call anon function in Dict with s struct passed in
   catch
      # control keys
      if mods & GLFW.MOD_CONTROL > 0
         if key == GLFW.KEY_A
            alignscale(s)
         elseif key == GLFW.KEY_P
            println(s.d.qpos)
            #elseif key == GLFW.KEY_L && lastfile[0]
            #   loadmodel(window, s.)
         elseif key == GLFW.KEY_Q
            if s.record != nothing
               close(s.record)
               s.record = nothing
            end
            GLFW.SetWindowShouldClose(window, true)
         end
      end

      # toggle visualiztion flag
      for i=1:Int(mj.NVISFLAG)
         if Int(key) == Int(mj.VISSTRING[i,3][1])
            flags = MVector(s.vopt[].flags)
            flags[i] = flags[i] == 0 ? 1 : 0
            s.vopt[].flags = flags
            return
         end
      end
      # toggle rendering flag
      for i=1:Int(mj.NRNDFLAG)
         if Int(key) == Int(mj.RNDSTRING[i,3][1])
            flags = MVector(s.scn[].flags)
            flags[i] = flags[i] == 0 ? 1 : 0
            s.scn[].flags = flags
            return
         end
      end
      # toggle geom/site group
      for i=1:Int(mj.NGROUP)
         if Int(key) == i + Int('0')
            if mods & GLFW.MOD_SHIFT == true
               sitegroup = MVector(s.vopt[].sitegroup)
               sitegroup[i] = sitegroup[i] > 0 ? 0 : 1
               s.vopt[].sitegroup[i] = sitegroup
               return
            else
               geomgroup = MVector(s.vopt[].geomgroup)
               geomgroup[i] = geomgroup[i] > 0 ? 0 : 1
               s.vopt[].geomgroup = geomgroup
               return
            end
         end
      end
   end
end

function mouse_move(s::mjSim, window::GLFW.Window,
                    xpos::Float64, ypos::Float64)
   # no buttons down: nothing to do
   if !s.button_left && !s.button_middle && !s.button_right
      return
   end

   # compute mouse displacement, save
   dx = xpos - s.lastx
   dy = ypos - s.lasty
   s.lastx = xpos
   s.lasty = ypos

   width, height = GLFW.GetWindowSize(window)

   mod_shift = GLFW.GetKey(window, GLFW.KEY_LEFT_SHIFT) || GLFW.GetKey(window, GLFW.KEY_RIGHT_SHIFT)

   # determine action based on mouse button
   if s.button_right
      action = mod_shift ? Int(mj.MOUSE_MOVE_H) : Int(mj.MOUSE_MOVE_V)
   elseif s.button_left
      action = mod_shift ? Int(mj.MOUSE_ROTATE_H) : Int(mj.MOUSE_ROTATE_V)
   else
      action = Int(mj.MOUSE_ZOOM)
   end

   # move perturb or camera
   if s.pert[].active != 0
      mjv_movePerturb(s.m.m, s.d.d, action,
                      dx / height, dy / height,
                      s.scn, s.pert);
   else
      mjv_moveCamera(s.m.m, action,
                     dx / height, dy / height,
                     s.scn, s.cam)
   end
end

# past data for double-click detection
function mouse_button(s::mjSim, window::GLFW.Window,
                      button::GLFW.MouseButton, act::GLFW.Action, mods::Int32)
   # update button state
   s.button_left = GLFW.GetMouseButton(window, GLFW.MOUSE_BUTTON_LEFT)
   s.button_middle = GLFW.GetMouseButton(window, GLFW.MOUSE_BUTTON_MIDDLE)
   s.button_right = GLFW.GetMouseButton(window, GLFW.MOUSE_BUTTON_RIGHT)

   # Alt: swap left and right
   if mods == GLFW.MOD_ALT
      tmp = s.button_left
      s.button_left = s.button_right
      s.button_right = tmp

      if button == GLFW.MOUSE_BUTTON_LEFT
         button = GLFW.MOUSE_BUTTON_RIGHT;
      elseif button == GLFW.MOUSE_BUTTON_RIGHT
         button = GLFW.MOUSE_BUTTON_LEFT;
      end
   end

   # update mouse position
   x, y = GLFW.GetCursorPos(window)
   s.lastx = x
   s.lasty = y

   # set perturbation
   newperturb = 0;
   if act == GLFW.PRESS && mods == GLFW.MOD_CONTROL && s.pert[].select > 0
      # right: translate;  left: rotate
      if s.button_right
         newperturb = Int(mj.PERT_TRANSLATE)
      elseif s.button_left
         newperturb = Int(mj.PERT_ROTATE)
      end
      # perturbation onset: reset reference
      if newperturb>0 && s.pert[].active==0
         mjv_initPerturb(s.m.m, s.d.d, s.scn, s.pert)
      end
   end
   s.pert[].active = newperturb

   # detect double-click (250 msec)
   if act == GLFW.PRESS && (time() - s.lastclicktm < 0.25) && (button == s.lastbutton)
      # determine selection mode
      if button == GLFW.MOUSE_BUTTON_LEFT
         selmode = 1;
      elseif mods == GLFW.MOD_CONTROL
         selmode = 3; # CTRL + Right Click
      else
         selmode = 2; # Right Click
      end
      # get current window size
      width, height = GLFW.GetWindowSize(window)

      # find geom and 3D click point, get corresponding body
      selpnt = zeros(3)
      selgeom, selskin = Int32(0), Int32(0)
      selbody = mjv_select(s.m.m, s.d.d, s.vopt,
                           width / height, x / width,
                           (height - y) / height,
                           s.scn, selpnt, selgeom, selskin)

      # set lookat point, start tracking is requested
      if selmode == 2 || selmode == 3
         # copy selpnt if geom clicked
         if selbody >= 0
            s.cam[].lookat = SVector{3,Float64}(selpnt...)
         end

         # switch to tracking camera
         if selmode == 3 && selbody >= 0
            s.cam[]._type = Int(mj.CAMERA_TRACKING)
            s.cam[].trackbodyid = selbody
            s.cam[].fixedcamid = -1
         end
      else # set body selection
         if selbody >= 0
            # compute localpos
            tmp = selpnt - s.d.xpos[:,selbody+1]
            res = reshape(s.d.xmat[:,selbody+1], 3, 3)' * tmp
            s.pert[].localpos = SVector{3}(res)

            # record selection
            s.pert[].select = selbody
            s.pert[].skinselect = selskin
         else
            s.pert[].select = 0
            s.pert[].skinselect = -1
         end
      end

      # stop perturbation on select
      s.pert[].active = 0
   end
   # save info
   if act == GLFW.PRESS
      s.lastbutton = button
      s.lastclicktm = time()
   end
end

function scroll(s::mjSim, window::GLFW.Window,
                xoffset::Float64, yoffset::Float64)
   # scroll: emulate vertical mouse motion = 5% of window height
   mjv_moveCamera(s.m.m, Int(mj.MOUSE_ZOOM),
                  0.0, -0.05 * yoffset, s.scn, s.cam);

end

function drop(window::GLFW.Window,
              count::Int, paths::String)
end

function start(mm::jlModel, dd::jlData, width=1200, height=900) # TODO named args for callbacks
   GLFW.WindowHint(GLFW.SAMPLES, 4)
   GLFW.WindowHint(GLFW.VISIBLE, 1)

   s = mjSim(mm, dd, "Simulate"; width=width, height=height)

   @info("Refresh Rate: $(s.refreshrate)")
   @info("Resolution: $(width)x$(height)")

   # Make the window's context current
   GLFW.MakeContextCurrent(s.window)
   GLFW.SwapInterval(1)

   # init abstract visualization
   mjv_defaultCamera(s.cam)
   mjv_defaultOption(s.vopt)
   #profilerinit();
   sensorinit(s)

   # make empty scene
   mjv_defaultScene(s.scn)
   mjv_makeScene(s.m, s.scn, maxgeom)

   # mujoco setup
   mjv_defaultPerturb(s.pert)
   mjr_defaultContext(s.con)
   mjr_makeContext(s.m, s.con, Int(fontscale)) # model specific setup

   alignscale(s)
   mjv_updateScene(s.m, s.d,
                   s.vopt, s.pert, s.cam, Int(mj.CAT_ALL), s.scn)

   GLFW.SetKeyCallback(s.window, (w,k,sc,a,m)->mykeyboard(s,w,k,sc,a,m))

   GLFW.SetCursorPosCallback(s.window, (w,x,y)->mouse_move(s,w,x,y))
   GLFW.SetMouseButtonCallback(s.window, (w,b,a,m)->mouse_button(s,w,b,a,m))
   GLFW.SetScrollCallback(s.window, (w,x,y)->scroll(s,w,x,y))
   ##GLFW.SetDropCallback(s.window, drop)

   return s
end

#### To customize what is rendered, change the following functions

function render(s::MujocoSim.mjSim, w::GLFW.Window)
   wi, hi = GLFW.GetFramebufferSize(w)
   rect = mjrRect(Cint(0), Cint(0), Cint(wi), Cint(hi))
   smallrect = mjrRect(Cint(0), Cint(0), Cint(wi), Cint(hi))

   # update scene
   mjv_updateScene(s.m, s.d,
                   s.vopt, s.pert, s.cam, Int(mj.CAT_ALL), s.scn)
   # render
   mjr_render(rect, s.scn, s.con)

   if s.record != nothing
      mjr_readPixels(s.vidbuff, C_NULL, rect, s.con);
      write(s.record, s.vidbuff[1:3*rect.width*rect.height]);
   end

   # Swap front and back buffers
   GLFW.SwapBuffers(w)
end

function loadmodel(modelfile, width, height)
   ptr_m = mj_loadXML(modelfile, C_NULL)
   ptr_d = mj_makeData(ptr_m)
   m, d = mj.mapmujoco(ptr_m, ptr_d)
   mj_forward(m, d)
   s = MujocoSim.start(m, d, width, height)
   @info("Model file: $modelfile")

   GLFW.SetWindowRefreshCallback(s.window, (w)->render(s,w))

   return s
end

function simstep(s::mjSim)
   d = s.d
   m = s.m
   if s.paused
      if s.pert[].active > 0
         mjv_applyPerturbPose(m, d, s.pert, 1)  # move mocap and dynamic bodies
         mj_forward(m, d)
      end
   else
      #slow motion factor: 10x
      factor = s.slowmotion ? 10 : 1

      # advance effective simulation time by 1/refreshrate
      startsimtm = d.d[].time
      starttm = time()
      refreshtm = 1.0/(factor*s.refreshrate)
      updates = refreshtm / m.m[].opt.timestep

      steps = round(Int, round(s.framecount+updates)-s.framecount)
      s.framecount += updates

      for i=1:steps
         # clear old perturbations, apply new
         d.xfrc_applied .= 0.0
         if s.pert[].select > 0
            mjv_applyPerturbPose(m, d, s.pert, 0) # move mocap bodies only
            mjv_applyPerturbForce(m, d, s.pert)
         end

         mj_step(s.m, s.d)

         # break on reset
         (d.d[].time < startsimtm) && break
      end
   end
end

function simulate(file)
   s = loadmodel(file, 1200, 900)
   # Loop until the user closes the window
   MujocoSim.alignscale(s)
   while !GLFW.WindowShouldClose(s.window)

      simstep(s)

      render(s, s.window)
      GLFW.PollEvents()
   end
   GLFW.DestroyWindow(s.window)
end

end


text:
using Random

include("get_div_grad.jl")

function test_cr()
  cr_tol = 1.0e-6

  # Cubic spline matrix _ case: ‖x*‖ > Δ
  n = 10
  A = spdiagm(-1 => ones(n-1), 0 => 4 * ones(n), 1 => ones(n-1))
  b = A * [1:n;]

  (x, stats) = cr(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("CR: Relative residual: %8.1e\n", resid)
  @test(resid <= cr_tol)
  @test(stats.solved)

  # Code coverage
  (x, stats) = cr(Matrix(A), b)
  show(stats)

  radius = 0.75 * norm(x)
  (x, stats) = cr(A, b, radius=radius)
  show(stats)
  @test(stats.solved)
  @test abs(norm(x) - radius) ≤ cr_tol * radius

  # Sparse Laplacian
  A = get_div_grad(16, 16, 16)
  b = randn(size(A, 1))
  itmax = 0
  # case: ‖x*‖ > Δ
  radius = 10.
  (x, stats) = cr(A, b, radius=radius)
  xNorm = norm(x)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("CR: Relative residual: %8.1e\n", resid)
  @test abs(xNorm - radius) ≤ cr_tol * radius
  @test(stats.solved)
  # case: ‖x*‖ < Δ
  radius = 30.
  (x, stats) = cr(A, b, radius=radius)
  xNorm = norm(x)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("CR: Relative residual: %8.1e\n", resid)
  @test(resid <= cr_tol)
  @test(stats.solved)

  radius = 0.75 * xNorm
  (x, stats) = cr(A, b, radius=radius)
  show(stats)
  @test(stats.solved)
  @test(abs(radius - norm(x)) <= cr_tol * radius)

  opA = LinearOperator(A)
  (xop, statsop) = cr(opA, b, radius=radius)
  @test(abs(radius - norm(xop)) <= cr_tol * radius)

  n = 100
  itmax = 2 * n
  B = LBFGSOperator(n)
  Random.seed!(0)
  for i = 1:5
    push!(B, rand(n), rand(n))
  end
  b = B * ones(n)
  (x, stats) = cr(B, b)
  @test norm(x - ones(n)) ≤ cr_tol * norm(x)
  @test stats.solved

  # Test b == 0
  (x, stats) = cr(A, zeros(size(A,1)))
  @test x == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = [4 -1 0; -1 4 -1; 0 -1 4]
  b = [7; 2; -1]
  (x, stats) = cr(A, b)
  @test stats.solved

  # Test with Jacobi (or diagonal) preconditioner
  A = ones(10,10) + 9 * I
  b = 10 * [1:10;]
  M = 1/10 * opEye(10)
  (x, stats) = cr(A, b, M=M, itmax=10);
  show(stats)
  r = b - A * x;
  resid = norm(r) / norm(b);
  @printf("CR: Relative residual: %8.1e\n", resid);
  @test(resid <= cr_tol);
  # @test(stats.solved);
end

test_cr()


text:
# SPDX-License-Identifier: X11
# 2020-08-29
# Multiplication 1 (100pt)

println(*(parse.(Int, split(readline(), " "))...))

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming