text:
"""
    qrdelete!(Q, R, k)

Delete the left-most column of F = Q[:, 1:k] * R[1:k, 1:k] by updating Q and R.
Only Q[:, 1:(k-1)] and R[1:(k-1), 1:(k-1)] are valid on exit.
"""
function qrdelete!(Q::AbstractMatrix, R::AbstractMatrix, k::Int)
  n, m = size(Q)
  m == LinearAlgebra.checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  # apply Givens rotations
  for i in 2:k
      g = first(givens(R, i - 1, i, i))
      lmul!(g, R)
      rmul!(Q, g')
  end

  # move columns of R
  @inbounds for j in 1:(k-1)
    for i in 1:(k-1)
      R[i, j] = R[i, j + 1]
    end
  end

  Q, R
end

"""
    qradd!(Q, R, v, k)

Replace the right-most column of F = Q[:, 1:k] * R[1:k, 1:k] with v by updating Q and R.
This implementation modifies vector v as well. Only Q[:, 1:k] and R[1:k, 1:k] are valid on
exit.
"""
function qradd!(Q::AbstractMatrix, R::AbstractMatrix, v::AbstractVector, k::Int)
  n, m = size(Q)
  n == length(v) || throw(DimensionMismatch())
  m == LinearAlgebra.checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  @inbounds for i in 1:(k-1)
    q = view(Q, :, i)
    r = dot(q, v)

    R[i, k] = r
    axpy!(-r, q, v)
  end

  @inbounds begin
    d = norm(v)
    R[k, k] = d
    @.. @view(Q[:, k]) = v / d
  end

  Q, R
end

function qradd!(Q::AbstractMatrix, R::AbstractMatrix, v::Number, k::Int)
  1 == LinearAlgebra.checksquare(Q) == LinearAlgebra.checksquare(R) ||
    throw(DimensionMismatch())
  k == 1 || throw(ArgumentError())

  R[1, 1] = abs(v)
  Q[1, 1] = one(v)

  Q, R
end

get_status(nlsolver::NLSolver) = nlsolver.status
nlsolvefail(nlsolver::NLSolver) = nlsolvefail(get_status(nlsolver))
nlsolvefail(status::NLStatus) = Int8(status) < 0

isnewton(nlsolver::NLSolver) = isnewton(nlsolver.cache)
isnewton(nlcache::Union{NLNewtonCache,NLNewtonConstantCache}) = true
isnewton(nlcache::AbstractNLSolverCache) = false

set_new_W!(nlsolver::NLSolver, val::Bool)::Bool = set_new_W!(nlsolver.cache, val)
set_new_W!(nlcache::NLNewtonCache, val::Bool)::Bool = nlcache.new_W = val
set_new_W!(nlcache::AbstractNLSolverCache, val::Bool)::Bool = val

get_W(nlsolver::NLSolver) = get_W(nlsolver.cache)
get_W(nlcache::Union{NLNewtonCache,NLNewtonConstantCache}) = nlcache.W

set_W!(nlsolver::NLSolver, W) = set_W!(nlsolver.cache, W)
set_W!(nlcache::Union{NLNewtonCache,NLNewtonConstantCache}, W) = (nlcache.W = W; W)

set_W_dt!(nlsolver::NLSolver, W_dt) = set_W_dt!(nlsolver.cache, W_dt)
set_W_dt!(nlcache::NLNewtonCache, W_dt) = (nlcache.W_dt = W_dt; W_dt)
set_W_dt!(nlcache::NLNewtonConstantCache, W_dt) = W_dt

function nlsolve_f end
function iip_get_uf end
function oop_get_uf end
function build_jac_config end
function resize_jac_config! end

DiffEqBase.@def iipnlsolve begin
  @unpack κ, fast_convergence_cutoff = alg.nlsolve

  # define additional fields of cache of non-linear solver
  z = similar(u); dz = similar(u); tmp = similar(u); b = similar(u)
  k = zero(rate_prototype)

  uTolType = real(uBottomEltypeNoUnits)

  # create cache of non-linear solver
  if alg.nlsolve isa NLNewton
    nf = nlsolve_f(f, alg)

    if islinear(f)
      # get the operator
      J = nf.f
      W = WOperator(f.mass_matrix, dt, J, true)
    else
      if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
        W = WOperator(f, dt, true)
        J = nothing # is J = W.J better?
      else
        J = false .* vec(u) .* vec(u)'
        W = similar(J)
      end
    end

    nlcache = DiffEqBase.NLNewtonCache(true,W,J,dt,alg.nlsolve.new_W_dt_cutoff)
  elseif alg.nlsolve isa NLFunctional
    z₊ = similar(z)

    nlcache = DiffEqBase.NLFunctionalCache(z₊)
  elseif alg.nlsolve isa NLAnderson
    z₊ = similar(z)

    max_history = min(alg.nlsolve.max_history, alg.nlsolve.max_iter, length(z))
    Δz₊s = [zero(z) for i in 1:max_history]
    Q = Matrix{uEltypeNoUnits}(undef, length(z), max_history)
    R = Matrix{uEltypeNoUnits}(undef, max_history, max_history)
    γs = Vector{uEltypeNoUnits}(undef, max_history)
    dzold = zero(z)
    z₊old = zero(z)

    nlcache = DiffEqBase.NLAndersonCache(z₊,dzold,z₊old,Δz₊s,Q,R,γs,alg.nlsolve.aa_start,alg.nlsolve.droptol)
  end

  # define additional fields of cache
  fsalfirst = zero(rate_prototype)
  if alg.nlsolve isa NLNewton
    if islinear(f)
      du1 = rate_prototype
      uf = nothing
      jac_config = nothing
      linsolve = alg.linsolve(Val{:init},nf,u)
    else
      du1 = zero(rate_prototype)
      # if the algorithm specializes on split problems the use `nf`
      uf = UJacobianWrapper(nf,t,p)
      jac_config = build_jac_config(alg,nf,uf,du1,uprev,u,tmp,dz)
      linsolve = alg.linsolve(Val{:init},uf,u)
    end
    # TODO: check if the solver is iterative
    weight = similar(u)
  else
    J = nothing
    W = nothing
    du1 = rate_prototype
    uf = nothing
    jac_config = nothing
    linsolve = nothing
    weight = z
  end

  # create non-linear solver
  nlsolver = NLSolver{true,typeof(z),typeof(k),uTolType,typeof(κ),typeof(γ),typeof(c),typeof(du1),typeof(uf),typeof(jac_config),typeof(linsolve),typeof(fast_convergence_cutoff),typeof(nlcache)}(z,dz,tmp,b,k,one(uTolType),κ,γ,c,alg.nlsolve.max_iter,10000,Convergence,fast_convergence_cutoff,du1,uf,jac_config,linsolve,weight,nlcache)
end

DiffEqBase.@def oopnlsolve begin
  @unpack κ, fast_convergence_cutoff = alg.nlsolve

  # define additional fields of cache of non-linear solver (all aliased)
  z = uprev; dz = z; tmp = z; b = z; k = rate_prototype

  uTolType = real(uBottomEltypeNoUnits)

  # create cache of non-linear solver
  if alg.nlsolve isa NLNewton
    nf = nlsolve_f(f, alg)

    # only use `nf` if the algorithm specializes on split eqs
    uf = UDerivativeWrapper(nf,t,p)

    if islinear(f) || DiffEqBase.has_jac(f)
      # get the operator
      J = islinear(f) ? nf.f : f.jac(uprev, p, t)
      if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
        J = DiffEqArrayOperator(J)
      end
      W = WOperator(f.mass_matrix, dt, J, false)
    else
      # https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/672
      if u isa StaticArray
        # get a "fake" `J`
        J = if u isa AbstractMatrix && size(u, 1) > 1 # `u` is already a matrix
          u
        elseif size(u, 1) == 1 # `u` is a row vector
          vcat(u, u)
        else # `u` is a column vector
          hcat(u, u)
        end
        W = lu(J)
      else
        W = u isa Number ? u : LU{LinearAlgebra.lutype(uEltypeNoUnits)}(Matrix{uEltypeNoUnits}(undef, 0, 0),
                                                                        Vector{LinearAlgebra.BlasInt}(undef, 0),
                                                                        zero(LinearAlgebra.BlasInt))
      end
    end

    nlcache = DiffEqBase.NLNewtonConstantCache(W,J,alg.nlsolve.new_W_dt_cutoff)
  elseif alg.nlsolve isa NLFunctional
    uf = nothing

    nlcache = DiffEqBase.NLFunctionalConstantCache()
  elseif alg.nlsolve isa NLAnderson
    uf = nothing

    max_history = min(alg.nlsolve.max_history, alg.nlsolve.max_iter, length(z))
    Δz₊s = Vector{typeof(z)}(undef, max_history)
    Q = Matrix{uEltypeNoUnits}(undef, length(z), max_history)
    R = Matrix{uEltypeNoUnits}(undef, max_history, max_history)
    γs = Vector{uEltypeNoUnits}(undef, max_history)

    nlcache = DiffEqBase.NLAndersonConstantCache(Δz₊s,Q,R,γs,alg.nlsolve.aa_start,alg.nlsolve.droptol)
  end

  # create non-linear solver
  nlsolver = NLSolver{false,typeof(z),typeof(k),uTolType,typeof(κ),typeof(γ),typeof(c),Nothing,typeof(uf),Nothing,Nothing,typeof(fast_convergence_cutoff),typeof(nlcache)}(z,dz,tmp,b,k,one(uTolType),κ,γ,c,alg.nlsolve.max_iter,10000,Convergence,fast_convergence_cutoff,nothing,uf,nothing,nothing,z,nlcache)
end

# No J version
function iipnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  iipnlsolve(alg,u,uprev,p,t,dt,f,W,nothing,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
end

function iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @unpack κ, fast_convergence_cutoff = alg.nlsolve

  # define additional fields of cache of non-linear solver
  z = similar(u); dz = similar(u); tmp = similar(u); b = similar(u)
  k = zero(rate_prototype)

  uTolType = real(uBottomEltypeNoUnits)

  # create cache of non-linear solver
  if alg.nlsolve isa NLNewton
    nlcache = DiffEqBase.NLNewtonCache(true,W,J,dt,alg.nlsolve.new_W_dt_cutoff)
  elseif alg.nlsolve isa NLFunctional
    z₊ = similar(z)

    nlcache = DiffEqBase.NLFunctionalCache(z₊)
  elseif alg.nlsolve isa NLAnderson
    z₊ = similar(z)

    max_history = min(alg.nlsolve.max_history, alg.nlsolve.max_iter, length(z))
    Δz₊s = [zero(z) for i in 1:max_history]
    Q = Matrix{uEltypeNoUnits}(undef, length(z), max_history)
    R = Matrix{uEltypeNoUnits}(undef, max_history, max_history)
    γs = Vector{uEltypeNoUnits}(undef, max_history)
    dzold = zero(z)
    z₊old = zero(z)

    nlcache = DiffEqBase.NLAndersonCache(z₊,dzold,z₊old,Δz₊s,Q,R,γs,alg.nlsolve.aa_start,alg.nlsolve.droptol)
  end

  # define additional fields of cache
  if alg.nlsolve isa NLNewton
    nf = nlsolve_f(f, alg)

    if islinear(f)
      du1 = rate_prototype
      uf = nothing
      jac_config = nothing
      linsolve = alg.linsolve(Val{:init},nf,u)
    else
      du1 = zero(rate_prototype)
      # if the algorithm specializes on split problems the use `nf`
      # we pass this `alg` here just for identification purpose, because get_uf would be overloaded in different repos
      uf = iip_get_uf(alg,nf,t,p)
      jac_config = build_jac_config(alg,nf,uf,du1,uprev,u,tmp,dz)
      linsolve = alg.linsolve(Val{:init},uf,u)
    end
    # TODO: check if the solver is iterative
    weight = similar(u)
  else
    du1 = rate_prototype
    uf = nothing
    jac_config = nothing
    linsolve = nothing
    weight = z
  end

  # create non-linear solver
  nlsolver = NLSolver{true,typeof(z),typeof(k),uTolType,typeof(κ),typeof(γ),typeof(c),typeof(du1),typeof(uf),typeof(jac_config),typeof(linsolve),typeof(fast_convergence_cutoff),typeof(nlcache)}(z,dz,tmp,b,k,one(uTolType),κ,γ,c,alg.nlsolve.max_iter,10000,Convergence,fast_convergence_cutoff,du1,uf,jac_config,linsolve,weight,nlcache)
end

DiffEqBase.@def getiipnlsolvefields begin
  @unpack z,dz,tmp,k,uf,du1,jac_config,linsolve = nlsolver
  b = nlsolver.ztmp
  fsalfirst = zero(rate_prototype)
end

# No J version
function oopnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  oopnlsolve(alg,u,uprev,p,t,dt,f,W,nothing,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
end

function oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @unpack κ, fast_convergence_cutoff = alg.nlsolve

  # define additional fields of cache of non-linear solver (all aliased)
  z = uprev; dz = z; tmp = z; b = z; k = rate_prototype

  uTolType = real(uBottomEltypeNoUnits)

  # create cache of non-linear solver
  if alg.nlsolve isa NLNewton
    nf = nlsolve_f(f, alg)
    # only use `nf` if the algorithm specializes on split eqs
    uf = oop_get_uf(alg,nf,t,p)
    nlcache = DiffEqBase.NLNewtonConstantCache(W,J,alg.nlsolve.new_W_dt_cutoff)
  elseif alg.nlsolve isa NLFunctional
    uf = nothing
    nlcache = DiffEqBase.NLFunctionalConstantCache()
  elseif alg.nlsolve isa NLAnderson
    uf = nothing
    max_history = min(alg.nlsolve.max_history, alg.nlsolve.max_iter, length(z))
    Δz₊s = Vector{typeof(z)}(undef, max_history)
    Q = Matrix{uEltypeNoUnits}(undef, length(z), max_history)
    R = Matrix{uEltypeNoUnits}(undef, max_history, max_history)
    γs = Vector{uEltypeNoUnits}(undef, max_history)

    nlcache = DiffEqBase.NLAndersonConstantCache(Δz₊s,Q,R,γs,alg.nlsolve.aa_start,alg.nlsolve.droptol)
  end

  # create non-linear solver
  nlsolver = NLSolver{false,typeof(z),typeof(k),uTolType,typeof(κ),typeof(γ),typeof(c),Nothing,typeof(uf),Nothing,Nothing,typeof(fast_convergence_cutoff),typeof(nlcache)}(z,dz,tmp,b,k,one(uTolType),κ,γ,c,alg.nlsolve.max_iter,10000,Convergence,fast_convergence_cutoff,nothing,uf,nothing,nothing,z,nlcache)
end

DiffEqBase.@def getoopnlsolvefields begin
  uf = nlsolver.uf
end

function nlsolve_resize!(integrator::DEIntegrator, i::Int)
  if !isdefined(integrator.cache, :nlsolver)
    return nothing
  end
  alg = integrator.alg; nlsolver = integrator.cache.nlsolver
  if nlsolver isa AbstractArray
    for idx in eachindex(nlsolver) # looping because we may have multiple nlsolver for threaded case
      _nlsolver = nlsolver[idx]
      @unpack z,dz,tmp,ztmp,k,du1,uf,jac_config,linsolve,weight,cache = _nlsolver
      # doubt: if these fields are always going to be in alg cache too, then we shouldnt do this here.
      # double resize doesn't do any bad I think though
      resize!(z,i)
      resize!(dz,i)
      resize!(tmp,i)
      resize!(ztmp,i)
      resize!(k,i)
      resize!(du1,i)
      if jac_config !== nothing
        _nlsolver.jac_config = resize_jac_config!(jac_config, i)
      end
      resize!(weight, i)
      nlsolve_cache_resize!(cache,alg,i)
    end
  else
    @unpack z,dz,tmp,ztmp,k,du1,uf,jac_config,linsolve,weight,cache = nlsolver
    resize!(z,i)
    resize!(dz,i)
    resize!(tmp,i)
    resize!(ztmp,i)
    resize!(k,i)
    resize!(du1,i)
    if jac_config !== nothing
      nlsolver.jac_config = resize_jac_config!(jac_config,i)
    end
    resize!(weight, i)
    nlsolve_cache_resize!(cache,alg,i)
  end
  nothing
end

function nlsolve_cache_resize!(cache::NLNewtonCache, alg, i::Int)
  nothing
end

function nlsolve_cache_resize!(cache::NLNewtonConstantCache, alg, i::Int)
  nothing
end

function nlsolve_cache_resize!(cache::NLAndersonCache, alg, i::Int)
  resize!(cache.z₊, i)
  resize!(cache.dzold, i)
  resize!(cache.z₊old, i)
  max_history = min(alg.nlsolve.max_history, alg.nlsolve.max_iter, i)
  prev_max_history = length(cache.Δz₊s)
  resize!(cache.γs, max_history)
  resize!(cache.Δz₊s, max_history)
  if max_history > prev_max_history
    for i in (max_history - prev_max_history):max_history
      cache.Δz₊s[i] = zero(z₊)
    end
  end
  cache.Q = typeof(cache.Q)(undef, i, max_history)
  cache.R = typeof(cache.R)(undef, max_history, max_history)
  nothing
end

function nlsolve_cache_resize!(cache::NLAndersonConstantCache, alg, i::Int)
  max_history = min(alg.nlsolve.max_history, alg.nlsolve.max_iter, i)
  resize!(cache.Δz₊s, max_history)
  cache.Q = typeof(cache.Q)(undef, i, max_history)
  cache.R = typeof(cache.R)(undef, max_history, max_history)
  resize!(cache.γs, max_history)
  nothing
end

function nlsolve_cache_resize!(cache::NLFunctionalCache, alg, i::Int)
  resize!(cache.z₊, i)
  nothing
end

function nlsolve_cache_resize!(cache::NLFunctionalConstantCache, alg, i::Int)
  nothing
end


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ ffb77e5e-3003-11eb-336e-955303c13b26
md"A notebook on Gases is here [here](/open?path=./chemistry/gases.jl)"

# ╔═╡ 4d847a06-3004-11eb-334a-dfdb8ee22cff
md"A notebook on kinetics is [here](/open?path=./chemistry/kinetics.jl)"

# ╔═╡ 016d1054-34d9-11eb-18a5-113cea681f28
md"A notebook on equilibrium is [here](/open?path=./chemistry/Equilibrium.jl)"

# ╔═╡ c73e418c-6f25-11eb-26d7-8972ef3a7920
md"A notebook on Balancing Chemical Equations is [here](/open?path=./chemistry/chemEquations.jl)"

# ╔═╡ c395c30a-2486-11eb-16ec-9968033a9d5a
md"### Acknowledgements:

##### [Julia chemistry resources](https://github.com/svaksha/Julia.jl/blob/master/Chemistry.md)

##### Domain specific language (DSL) for high performance simulation and modeling of chemical reaction networks [Catalyst.jl] (https://catalyst.sciml.ai/stable/)

##### [Julia packages tagged as chemistry] (https://juliaobserver.com/categories/Chemistry)
"

# ╔═╡ Cell order:
# ╟─ffb77e5e-3003-11eb-336e-955303c13b26
# ╟─4d847a06-3004-11eb-334a-dfdb8ee22cff
# ╟─016d1054-34d9-11eb-18a5-113cea681f28
# ╟─c73e418c-6f25-11eb-26d7-8972ef3a7920
# ╟─c395c30a-2486-11eb-16ec-9968033a9d5a

ext: jl


text:
# saxpy.jl
using LIKWID
using CUDA
using LinearAlgebra

@assert CUDA.functional()

N = 100_000_000
a = 3.141f0
z = zeros(Float32, N)
x = rand(Float32, N)
y = rand(Float32, N)

z_gpu = CUDA.zeros(Float32, N)
x_gpu = CUDA.rand(Float32, N)
y_gpu = CUDA.rand(Float32, N)

function saxpy_cpu!(z,a,x,y)
    z .= a .* x .+ y
end

function saxpy_gpu!(z,a,x,y)
    CUDA.@sync z .= a .* x .+ y
end

Marker.init()
GPUMarker.init()

saxpy_cpu!(z,a,x,y)
@region "saxpy_cpu" saxpy_cpu!(z,a,x,y)

saxpy_gpu!(z_gpu,a,x_gpu,y_gpu)
@gpuregion "saxpy_gpu" saxpy_gpu!(z_gpu,a,x_gpu,y_gpu)

Marker.close()
GPUMarker.close()


text:
module RediSearch

using Base: String, Number
using Jedis

include("utilities.jl")
include("document.jl")
include("query.jl")
include("command.jl")
include("fields.jl")
include("client.jl")
include("result.jl")
include("index_definition.jl")
include("search.jl")

end # module


text:
export bissect
function bissect(h :: AbstractNLPModel;
                t1 :: Float64 = h.meta.x0[1],
                t2 :: Float64 = 100.0,
                tol :: Float64 = 1e-7,
                maxiter :: Int = 50,
                verbose :: Bool = true)

    (length(h.meta.x0) > 1) && warn("Not a 1-D problem ")
    (tₐ, tᵦ) = trouve_intervalle(h, t1, 1.0)
    gₐ = grad(h, [tₐ])[1]
    gᵦ = grad(h, [tᵦ])[1]
    iter = 0

    verbose &&
        @printf(" iter        tₐ        tᵦ         gₐ        gᵦ        \n")
    verbose &&
        @printf(" %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n", iter, tₐ, tᵦ, gₐ, gᵦ)
    while ((max(abs(gᵦ), abs(gₐ))) > tol) && (iter < maxiter) || iter == 0
      tₚ = (tₐ + tᵦ) / 2
      gₚ = grad(h, [tₚ])[1]
      if gₚ <= 0.0
        tₐ = tₚ
        gₐ = gₚ
      else
        tᵦ = tₚ
        gᵦ = gₚ
      end

      iter += 1
      verbose &&
        @printf(" %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n", iter, tₐ, tᵦ, gₐ, gᵦ)
    end

    topt = (tₐ + tᵦ) / 2.0
    gopt = grad(h, [topt])[1]
    status = :NotSolved
    (abs(gopt) < tol) && (status = :Optimal)
    (iter >= maxiter) && (status = :Tired)
    tired = iter > maxiter
    optimal = abs(gopt) < tol
    return (topt, obj(h, [topt])[1], norm(gopt, Inf), iter, optimal, tired, status, h.counters.neval_obj, h.counters.neval_grad, h.counters.neval_hess)
end


text:
"""
Performs loop distribution for any top-level for loops in `statements`
"""
function distribute_batch_loops(statements)
    new_body = []
    for statement in statements
        if isa(statement, Expr) && statement.head == :for
            for expr in statement.args[2].args
                loop = deepcopy(statement)
                if isa(expr, LineNumberNode) || isa(expr, LabelNode) || expr.head == :line
                    continue
                end
                loop.args[2].args = [expr]
                push!(new_body, loop)
            end
        else
            push!(new_body, statement)
        end
    end
    new_body
end


text:
"""

```
partialsort_cutoff!(x,cutoff; by = x -> x)
```

Function that reorders x vector by putting in the first positions the
elements with values smaller than cutoff

"""
function partialsort_cutoff!(x, cutoff; by = x -> x)
    iswap = 1
    for i = 1:length(x)
        if by(x[i]) < cutoff
            if iswap == i
                iswap = iswap + 1
            else
                xtmp = x[iswap]
                x[iswap] = x[i]
                x[i] = xtmp
                iswap = iswap + 1
            end
        end
    end
    nothing
end


text:
@testset "Capsule" begin
    @testset "constructors" begin
        o, extr, r = Point2f0(1, 2), Point2f0(3, 4), 5f0
        s = Capsule(o, extr, r)
        @test typeof(s) == Capsule{2,Float32}
        @test typeof(s) == Capsule2{Float32}
        @test origin(s) == o
        @test extremity(s) == extr
        @test radius(s) == r
        #@test abs(height(s)- norm([1,2]-[3,4]))<1e-5
        h = norm(o - extr)
        @test isapprox(height(s), h)
        #@test norm(direction(s) - Point{2,Float32}([2,2]./norm([1,2]-[3,4])))<1e-5
        @test isapprox(direction(s), Point2f0(2, 2) ./ h)
        v1 = rand(Point{3, Float64}); v2 = rand(Point{3, Float64}); R = rand()
        s = Capsule(v1, v2, R)
        @test typeof(s) == Capsule{3, Float64}
        @test typeof(s) == Capsule3{Float64}
        @test origin(s) == v1
        @test extremity(s) == v2
        @test radius(s) == R
        @test height(s) == norm(v2 - v1)
        #@test norm(direction(s) - Point{3,Float64}((v2-v1)./norm(v2-v1)))<1e-10
        @test isapprox(direction(s), (v2-v1) ./ norm(v2 .- v1))
      end

    @testset "decompose" begin
        o, extr, r = Point2f0(1, 2), Point2f0(3, 4), 5f0
        s = Capsule(o, extr, r)
        positions = Point{3, Float32}[
          (4.535534, -1.5355341, 0.0),
          (2.9134173, -2.619398, 0.0),
          (1.0000002, -3.0, 0.0),
          (-0.913417, -2.6193976, 0.0),
          (-2.5355341, -1.5355337, 0.0),
          (-3.6193976, 0.08658302, 0.0),
          (-4.0, 2.0000002, 0.0),
          (-3.619398, 3.9134173, 0.0),
          (-2.535534, 5.5355344, 0.0),
          (-0.53553426, 7.535534, 0.0),
          (1.0865824, 8.619398, 0.0),
          (2.9999998, 9.0, 0.0),
          (4.913417, 8.619397, 0.0),
          (6.535534, 7.535534, 0.0),
          (7.619397, 5.913417, 0.0),
          (8.0, 3.9999998, 0.0),
          (7.619398, 2.0865827, 0.0),
          (6.535534, 0.4644655, 0.0)
        ]
        @test decompose(Point3f0, s, 4) ≈ positions
        FT = Face{3, Int}
        faces = FT[
          (1, 3, 2),
          (1, 4, 3),
          (1, 5, 4),
          (1, 6, 5),
          (1, 7, 6),
          (1, 8, 7),
          (1, 9, 8),
          (1, 10, 9),
          (1, 11, 10),
          (1, 12, 11),
          (1, 13, 12),
          (1, 14, 13),
          (1, 15, 14),
          (1, 16, 15),
          (1, 17, 16),
          (1, 18, 17)
        ]
        @test faces == decompose(FT, s, 4)

        v1 = Point{3, Float64}(1,2,3); v2 = Point{3, Float64}(4,5,6); R = 5.0
        s = Capsule(v1, v2, R)
        positions = Point{3,Float64}[
            (4.535533905932738, -1.5355339059327373, 3.0),
            (3.1616954987389683, -2.3711193256429137, 1.8952880865480273),
            (1.4587585476806852, -2.5412414523193148, 0.9587585476806852),
            (-0.31402023360539255, -2.020000734336378, 0.332989516029115),
            (7.535533905932738, 1.4644660940672627, 6.0),
            (8.371119325642914, 2.8383045012610317, 7.104711913451973),
            (8.541241452319316, 4.541241452319315, 8.041241452319316),
            (8.020000734336378, 6.314020233605392, 8.667010483970884),
            (3.0412414523193148, 4.041241452319315, -1.0824829046386295),
            (1.7811492852594561, 2.7811492852594557, -1.8764343108748296),
            (0.4021342206547498, 1.4021342206547494, -1.927992798267443),
            (-0.8858611987114287, 0.11413880128857112, -1.2293090544897973),
            (6.041241452319315, 7.041241452319315, 1.9175170953613705),
            (6.990573112163402, 7.990573112163402, 3.3329895160291154),
            (7.484617125293379, 8.48461712529338, 5.1544901063711865),
            (7.448159769230341, 8.44815976923034, 7.104711913451973),
            (-2.535533905932737, 5.535533905932738, 2.9999999999999996),
            (-3.3711193256429137, 4.161695498738968, 1.8952880865480268),
            (-3.5412414523193148, 2.4587585476806852, 0.9587585476806848),
            (-3.020000734336378, 0.6859797663946074, 0.332989516029115),
            (0.46446609406726314, 8.535533905932738, 6.0),
            (1.8383045012610322, 9.371119325642914, 7.104711913451972),
            (3.5412414523193148, 9.541241452319316, 8.041241452319316),
            (5.314020233605392, 9.020000734336378, 8.667010483970884),
            (-1.0412414523193152, -0.04124145231931431, 7.0824829046386295),
            (-1.990573112163402, -0.9905731121634007, 5.667010483970884),
            (-2.4846171252933793, -1.4846171252933784, 3.845509893628814),
            (-2.4481597692303416, -1.4481597692303412, 1.8952880865480275),
            (1.9587585476806848, 2.9587585476806857, 10.08248290463863),
            (3.2188507147405434, 4.218850714740545, 10.87643431087483),
            (4.59786577934525, 5.597865779345251, 10.927992798267443),
            (5.885861198711429, 6.885861198711429, 10.229309054489796),
            (-1.8867513459481287, -0.8867513459481287, 0.11324865405187134),
            (6.886751345948129, 7.886751345948129, 8.886751345948129)
        ]
        @test decompose(Point3{Float64},s,(4, 4)) ≈ positions
        faces = Face[
          (13, 9, 1),
          (5, 13, 1),
          (4, 12, 33),
          (16, 8, 34),
          (9, 10, 1),
          (10, 2, 1),
          (14, 13, 5),
          (14, 5, 6),
          (10, 11, 2),
          (11, 3, 2),
          (15, 14, 6),
          (15, 6, 7),
          (11, 12, 3),
          (12, 4, 3),
          (16, 15, 7),
          (16, 7, 8),
          (21, 17, 9),
          (13, 21, 9),
          (12, 20, 33),
          (24, 16, 34),
          (17, 18, 9),
          (18, 10, 9),
          (22, 21, 13),
          (22, 13, 14),
          (18, 19, 10),
          (19, 11, 10),
          (23, 22, 14),
          (23, 14, 15),
          (19, 20, 11),
          (20, 12, 11),
          (24, 23, 15),
          (24, 15, 16),
          (29, 25, 17),
          (21, 29, 17),
          (20, 28, 33),
          (32, 24, 34),
          (25, 26, 17),
          (26, 18, 17),
          (30, 29, 21),
          (30, 21, 22),
          (26, 27, 18),
          (27, 19, 18),
          (31, 30, 22),
          (31, 22, 23),
          (27, 28, 19),
          (28, 20, 19),
          (32, 31, 23),
          (32, 23, 24),
          (5, 1, 25),
          (29, 5, 25),
          (28, 4, 33),
          (8, 32, 34),
          (1, 2, 25),
          (2, 26, 25),
          (6, 5, 29),
          (6, 29, 30),
          (2, 3, 26),
          (3, 27, 26),
          (7, 6, 30),
          (7, 30, 31),
          (3, 4, 27),
          (4, 28, 27),
          (8, 7, 31),
          (8, 31, 32)
        ]
        @test faces == decompose(Face{3, Int}, s, (4, 4))
        m = GLPlainMesh(s, (4, 4))
        @test m.faces == faces
        @test m.vertices ≈ positions
        m = GLNormalMesh(s)# just test that it works without explicit resolution parameter
        @test m isa GLNormalMesh

    end
end


text:
"""
    SemidirectProductOperation(action::AbstractGroupAction)

Group operation of a semidirect product group. The operation consists of the operation
`opN` on a normal subgroup `N`, the operation `opH` on a subgroup `H`, and an automorphism
`action` of elements of `H` on `N`. Only the action is stored.
"""
struct SemidirectProductOperation{A<:AbstractGroupAction} <: AbstractGroupOperation
    action::A
    function SemidirectProductOperation{A}(action::A) where {A<:AbstractGroupAction}
        return new(action)
    end
end
function SemidirectProductOperation(action::A) where {A<:AbstractGroupAction}
    return SemidirectProductOperation{A}(action)
end

function Base.show(io::IO, op::SemidirectProductOperation)
    return print(io, "SemidirectProductOperation($(op.action))")
end

const SemidirectProductGroup{𝔽,N,H,A} =
    GroupManifold{𝔽,ProductManifold{𝔽,Tuple{N,H}},SemidirectProductOperation{A}}

@doc raw"""
    SemidirectProductGroup(N::GroupManifold, H::GroupManifold, A::AbstractGroupAction)

A group that is the semidirect product of a normal group ``\mathcal{N}`` and a subgroup
``\mathcal{H}``, written ``\mathcal{G} = \mathcal{N} ⋊_θ \mathcal{H}``, where
``θ: \mathcal{H} × \mathcal{N} → \mathcal{N}`` is an automorphism action of ``\mathcal{H}`` on
``\mathcal{N}``. The group ``\mathcal{G}`` has the composition rule

````math
g \circ g' = (n, h) \circ (n', h') = (n \circ θ_h(n'), h \circ h')
````

and the inverse

````math
g^{-1} = (n, h)^{-1} = (θ_{h^{-1}}(n^{-1}), h^{-1}).
````
"""
function SemidirectProductGroup(
    N::GroupManifold{𝔽},
    H::GroupManifold{𝔽},
    A::AbstractGroupAction,
) where {𝔽}
    N === g_manifold(A) || error("Subgroup $(N) must be the G-manifold of action $(A)")
    H === base_group(A) || error("Subgroup $(H) must be the base group of action $(A)")
    op = SemidirectProductOperation(A)
    M = ProductManifold(N, H)
    return GroupManifold(M, op)
end

function allocate_result(G::SemidirectProductGroup, ::typeof(identity_element))
    M = base_manifold(G)
    N, H = M.manifolds
    np = allocate_result(N, identity_element)
    hp = allocate_result(H, identity_element)
    return ProductRepr(np, hp)
end

"""
    identity_element(G::SemidirectProductGroup)

Get the identity element of [`SemidirectProductGroup`](@ref) `G`. Uses [`ProductRepr`](@ref)
to represent the point.
"""
identity_element(G::SemidirectProductGroup)

function identity_element!(G::SemidirectProductGroup, q)
    M = base_manifold(G)
    N, H = M.manifolds
    nq, hq = submanifold_components(G, q)
    identity_element!(N, nq)
    identity_element!(H, hq)
    @inbounds _padpoint!(G, q)
    return q
end

function is_identity(G::SemidirectProductGroup, p; kwargs...)
    M = base_manifold(G)
    N, H = M.manifolds
    nq, hq = submanifold_components(G, p)
    return is_identity(N, nq; kwargs...) && is_identity(H, hq; kwargs...)
end
function is_identity(G::SemidirectProductGroup, e::Identity; kwargs...)
    return invoke(is_identity, Tuple{AbstractGroupManifold,typeof(e)}, G, e; kwargs...)
end

function Base.show(io::IO, G::SemidirectProductGroup)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    return print(io, "SemidirectProductGroup($(N), $(H), $(A))")
end

submanifold(G::SemidirectProductGroup, i) = submanifold(base_manifold(G), i)

_padpoint!(::SemidirectProductGroup, q) = q

_padvector!(::SemidirectProductGroup, X) = X

function inv!(G::SemidirectProductGroup, q, p)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    inv!(H, hq, hp)
    npinv = inv(N, np)
    apply!(A, nq, hq, npinv)
    @inbounds _padpoint!(G, q)
    return q
end
function inv!(G::SemidirectProductGroup, q, ::Identity{<:SemidirectProductOperation})
    return identity_element!(G, q)
end
function inv!(
    G::SemidirectProductGroup,
    q::Identity{<:SemidirectProductOperation},
    ::Identity{<:SemidirectProductOperation},
)
    return q
end

function _compose!(G::SemidirectProductGroup, x, p, q)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    nx, hx = submanifold_components(G, x)
    compose!(H, hx, hp, hq)
    nxtmp = apply(A, hp, nq)
    compose!(N, nx, np, nxtmp)
    @inbounds _padpoint!(G, x)
    return x
end

@doc raw"""
    translate_diff(G::SemidirectProductGroup, p, q, X, conX::LeftAction)

Perform differential of the left translation on the semidirect product group `G`.

Since the left translation is defined as (cf. [`SemidirectProductGroup`](@ref)):

````math
L_{(n', h')} (n, h) = ( L_{n'} θ_{h'}(n), L_{h'} h)
````

then its differential can be computed as

````math
\mathrm{d}L_{(n', h')}(X_n, X_h) = ( \mathrm{d}L_{n'} (\mathrm{d}θ_{h'}(X_n)), \mathrm{d}L_{h'} X_h).
````
"""
translate_diff(G::SemidirectProductGroup, p, q, X, conX::LeftAction)

function translate_diff!(G::SemidirectProductGroup, Y, p, q, X, conX::LeftAction)
    M = base_manifold(G)
    N, H = M.manifolds
    A = G.op.action
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    nX, hX = submanifold_components(G, X)
    nY, hY = submanifold_components(G, Y)
    translate_diff!(H, hY, hp, hq, hX, conX)
    nZ = apply_diff(A, hp, nq, nX)
    nr = apply(A, hp, nq)
    translate_diff!(N, nY, np, nr, nZ, conX)
    @inbounds _padvector!(G, Y)
    return Y
end

# We need to prevent decorator unwrapping so that the correct `get_vector!` gets called
# and applies proper padding to the result if `X` happens to be a matrix.
# Otherwise rare random bugs happen where the padding is not applied.
function get_vector(G::SemidirectProductGroup, p, X, B::VeeOrthogonalBasis)
    Y = allocate_result(G, get_vector, p, X)
    return get_vector!(G, Y, p, X, B)
end

function get_vector!(G::SemidirectProductGroup, Y, p, X, B::VeeOrthogonalBasis)
    M = base_manifold(G)
    N, H = M.manifolds
    dimN = manifold_dimension(N)
    dimH = manifold_dimension(H)
    @assert length(X) == dimN + dimH
    np, hp = submanifold_components(G, p)
    nY, hY = submanifold_components(G, Y)
    get_vector!(N, nY, np, view(X, 1:dimN), B)
    get_vector!(H, hY, hp, view(X, (dimN + 1):(dimN + dimH)), B)
    @inbounds _padvector!(G, Y)
    return Y
end

function get_coordinates!(G::SemidirectProductGroup, Y, p, X, B::VeeOrthogonalBasis)
    M = base_manifold(G)
    N, H = M.manifolds
    dimN = manifold_dimension(N)
    dimH = manifold_dimension(H)
    @assert length(Y) == dimN + dimH
    np, hp = submanifold_components(G, p)
    nY, hY = submanifold_components(G, X)
    get_coordinates!(N, view(Y, 1:dimN), np, nY, B)
    get_coordinates!(H, view(Y, (dimN + 1):(dimN + dimH)), hp, hY, B)
    return Y
end

function zero_vector(G::SemidirectProductGroup, p)
    X = allocate_result(G, zero_vector, p)
    zero_vector!(G, X, p)
    return X
end

function zero_vector!(G::SemidirectProductGroup, X, p)
    M = base_manifold(G)
    N, H = M.manifolds
    np, hp = submanifold_components(G, p)
    nX, hX = submanifold_components(G, X)
    zero_vector!(N, nX, np)
    zero_vector!(H, hX, hp)
    return X
end

function Base.isapprox(
    G::SemidirectProductGroup,
    p::AbstractArray,
    q::AbstractArray;
    kwargs...,
)
    M = base_manifold(G)
    N, H = M.manifolds
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    return isapprox(N, np, nq; kwargs...) && isapprox(H, hp, hq; kwargs...)
end

function Base.isapprox(
    G::SemidirectProductGroup,
    p,
    X::AbstractMatrix,
    Y::AbstractMatrix;
    kwargs...,
)
    M = base_manifold(G)
    N, H = M.manifolds
    np, hp = submanifold_components(G, p)
    nX, hX = submanifold_components(G, X)
    nY, hY = submanifold_components(G, Y)
    return isapprox(N, np, nX, nY; kwargs...) && isapprox(H, hp, hX, hY; kwargs...)
end
function isapprox(
    G::SemidirectProductGroup{𝔽,N,H,A},
    ::Identity{SemidirectProductOperation{A}},
    X::AbstractMatrix,
    Y::AbstractMatrix;
    kwargs...,
) where {𝔽,N<:AbstractManifold,H<:AbstractManifold,A<:AbstractGroupAction}
    return isapprox(G, identity_element(G), X, Y; kwargs...)
end


text:
using BenchmarkTools

function compute_px1(siz,px,n,hx,var)
   for k=1:siz[3], j=1:siz[2], i=1:n-2
         px[i+1,j,k] = (var[i+2,j,k,1] - var[i,j,k,1])/(hx[i+2] - hx[i])
   end
end
function compute_px2(siz,px,n,hx,var)
   for i=1:n-2, j=1:siz[2], k=1:siz[3]
         px[i+1,j,k] = (var[i+2,j,k,1] - var[i,j,k,1])/(hx[i+2] - hx[i])
   end
end

function compute_px3(siz,px,n,hx,var)
   @inbounds for k=1:siz[3], j=1:siz[2], i=1:n-2
         px[i+1,j,k] = (var[i+2,j,k,1] - var[i,j,k,1])/(hx[i+2] - hx[i])
   end
end
function compute_px4(siz,px,n,hx,var)
   @inbounds for i=1:n-2, j=1:siz[2], k=1:siz[3]
         px[i+1,j,k] = (var[i+2,j,k,1] - var[i,j,k,1])/(hx[i+2] - hx[i])
   end
end
compute_px5(siz,px,n,hx,var) = @. px[2:n-1,:,:] = (var[3:n,:,:,1] - var[1:n-2,:,:,1])/(hx[3:n] - hx[1:n-2])
compute_px6(siz,px,n,hx,var) = @inbounds @. px[2:n-1,:,:] = (var[3:n,:,:,1] - var[1:n-2,:,:,1])/(hx[3:n] - hx[1:n-2])
compute_px7(siz,px,n,hx,var) = @views @. px[2:n-1,:,:] = (var[3:n,:,:,1] - var[1:n-2,:,:,1])/(hx[3:n] - hx[1:n-2])

function foo(hx::StepRangeLen,hy::StepRangeLen,hz::StepRangeLen,
   var::Array{Float64,4})::Array{Float64,3}

   siz = size(var)::NTuple{4,Int64}

   px = zeros(Float64,siz[1:3]...)
   qy = zeros(Float64,siz[1:3]...)
   rz = zeros(Float64,siz[1:3]...)

   n = size(hx,1)

   if n > 2
      @btime compute_px1($siz,$px,$n,$hx,$var)
      @btime compute_px2($siz,$px,$n,$hx,$var)
      @btime compute_px3($siz,$px,$n,$hx,$var)
      @btime compute_px4($siz,$px,$n,$hx,$var)
      @btime compute_px5($siz,$px,$n,$hx,$var)
      @btime compute_px6($siz,$px,$n,$hx,$var)
      @btime compute_px7($siz,$px,$n,$hx,$var)
      @time compute_px1(siz,px,n,hx,var)
      @time compute_px2(siz,px,n,hx,var)
      @time compute_px3(siz,px,n,hx,var)
      @time compute_px4(siz,px,n,hx,var)
      @time compute_px5(siz,px,n,hx,var)
      @time compute_px6(siz,px,n,hx,var)
      @time compute_px7(siz,px,n,hx,var)
   end

   px
end

function foo_stable(hx::StepRangeLen,hy::StepRangeLen,hz::StepRangeLen,
   var::Array{Float64,4})::Array{Float64,3}

   siz = size(var)::NTuple{4,Int64}

   px = zeros(Float64,Base.front(siz))
   qy = zeros(Float64,Base.front(siz))
   rz = zeros(Float64,Base.front(siz))

   n = size(hx,1)

   if n > 2
      @btime compute_px1($siz,$px,$n,$hx,$var)
      @btime compute_px2($siz,$px,$n,$hx,$var)
      @btime compute_px3($siz,$px,$n,$hx,$var)
      @btime compute_px4($siz,$px,$n,$hx,$var)
      @btime compute_px5($siz,$px,$n,$hx,$var)
      @btime compute_px6($siz,$px,$n,$hx,$var)
      @btime compute_px7($siz,$px,$n,$hx,$var)
      @time compute_px1(siz,px,n,hx,var)
      @time compute_px2(siz,px,n,hx,var)
      @time compute_px3(siz,px,n,hx,var)
      @time compute_px4(siz,px,n,hx,var)
      @time compute_px5(siz,px,n,hx,var)
      @time compute_px6(siz,px,n,hx,var)
      @time compute_px7(siz,px,n,hx,var)
   end

   px
end

function test()
   x = range(0.0, 1.0, length=10000)
   y = 0.5:0.5:1.5
   z = 0.5:0.5:1.5
   var = zeros(10000,3,3,3)
   px = foo(x,y,z,var);
   px = foo_stable(x,y,z,var);
end
test()


text:
module BFREPL

import REPL
import REPL.LineEdit

include(joinpath(@__DIR__(), "BFParser.jl"))

do_cmd(repl, input) = do_cmd(repl, String(input))
function do_cmd(repl, input::String)
    terminal = REPL.outstream(repl)
    bfi = BFParser.BFInterpreter(terminal, terminal)
    try
        bfs = BFParser.parseBF(BFParser.BF, input)
        BFParser.run(bfi, bfs)
    catch ex
        err_stream = isa(terminal, REPL.Terminals.TTYTerminal) ? terminal.err_stream : stderr
        if isa(ex, BFParser.ParseError)
            Base.display_error(err_stream, ErrorException(ex.msg), ex.stacktrace)
        elseif isa(ex, BFParser.ErrorWrapper)
            Base.display_error(err_stream, ex.err, ex.backtraces)
        elseif !isa(ex, InterruptException)
            rethrow()
        end
    end
end

function create_mode(repl, main_mode)
    bf_mode = LineEdit.Prompt("bf> ";
        prompt_prefix = Base.text_colors[:magenta],
        prompt_suffix = "",
        # complete = nothing,
        sticky = false)
    bf_mode.repl = repl
    hp = main_mode.hist
    hp.mode_mapping[:bf] = bf_mode
    bf_mode.hist = hp

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    prefix_prompt, prefix_keymap = LineEdit.setup_prefix_keymap(hp, bf_mode)

    bf_mode.on_done = (s, buf, ok) -> begin
        ok || return REPL.transition(s, :abort)
        input = String(take!(buf))
        REPL.reset(repl)
        do_cmd(repl, input)
        REPL.prepare_next(repl)
        REPL.reset_state(s)
        s.current_mode.sticky || REPL.transition(s, main_mode)
    end

    mk = REPL.mode_keymap(main_mode)

    # repl_keymap = Dict()

    b = Dict{Any, Any}[
        skeymap, #=repl_keymap,=# mk, prefix_keymap, LineEdit.history_keymap,
        LineEdit.default_keymap, LineEdit.escape_defaults
    ]
    bf_mode.keymap_dict = LineEdit.keymap(b)
    return bf_mode
end

function repl_init(repl)
    main_mode = repl.interface.modes[1]
    bf_mode = create_mode(repl, main_mode)
    push!(repl.interface.modes, bf_mode)
    org_action = main_mode.keymap_dict['\x02']
    keymap = Dict{Any, Any}(
        # "^B" => function (s, args...)
        '\x02' => function (s, args...)
            # println("^B")
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, bf_mode) do
                    LineEdit.state(s, bf_mode).input_buffer = buf
                end
            else
                # LineEdit.edit_insert(s, ']')
                org_action(s, args...)
            end
        end
    )
    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, keymap)
    return
end

function __init__()
    if isdefined(Base, :active_repl)
        repl_init(Base.active_repl)
    else
        atreplinit() do repl
            if isinteractive() && repl isa REPL.LineEditREPL
                isdefined(repl, :interface) || (repl.interface = REPL.setup_interface(repl))
                repl_init(repl)
            end
        end
    end
end

end # module


text:
using SparseArrays
using LinearAlgebra

struct Equation{T}
    A :: SparseMatrixCSC{Scalar,Label}
    x :: Array{T,1}
    b :: Array{T,1}
end

Ac(eq::Equation{T}) where {T} = Array(diag(eq.A));
H(eq::Equation{T}) where {T} =  -eq.b .+ Ac(eq) .* eq.x .-  eq.A * eq.x;

function *(A::SparseMatrixCSC{Scalar, Label}, x::VectorList)
    n,m = size(A)
    @assert length(x) == m

    y = zeros(Vector, n)
    rows = rowvals(A)
    vals = nonzeros(A)
    for c=1:m
        for k ∈ nzrange(A,c)
            r = rows[k]
            a = vals[k]
            y[r] += a*x[c]
        end
    end
    return y
end

# Scitani a odcitani rovnic
for op ∈ [:+, :-]
    @eval begin
        function ($op)(e1::Equation{T}, e2::Equation{T}) where {T}
            @assert e1.x == e2.x
            Equation{T}(($op)(e1.A, e2.A), e1.x, $(op)(e1.b, e2.b))
        end
    end
end

# Pricteni/odecteni pole o spravne velikosti
for op ∈ [:+, :-]
    @eval begin
        function ($op)(e1::Equation{T}, f::Array{T,1}) where {T}
            @assert size(e1.x) == size(f)
            Equation{T}(copy(e1.A), e1.x, $(op)(e1.b, f))
        end
    end
end

# Vynasobeni rovnice konstantou
*(a::Scalar, eq::Equation{T}) where {T} = Equation{T}(a*eq.A, eq.x, a*eq.b);

# Unarni minus
-(eq::Equation{T}) where {T} = Equation{T}(-eq.A, eq.x, -eq.b);

function solve!(eqn::Equation{Scalar})
    x = eqn.A \ eqn.b
    for i in eachindex(eqn.x)
        eqn.x[i] = -x[i]
    end
end;

function solve!(eqn::Equation{Vector})
    rhs = zeros(length(eqn.b),2)
    for i in eachindex(eqn.b)
        rhs[i,1], rhs[i,2] = eqn.b[i][1], eqn.b[i][2]
    end
    x = eqn.A \ rhs
    for i in eachindex(eqn.x)
        eqn.x[i] = -Vector(x[i,1], x[i,2])
    end
end;

function relax!(eqn::Equation{T}, α) where {T}
    D = diag(eqn.A)
    for i=1:length(D)
        eqn.A[i,i] /= α
    end
    eqn.b[:] -= (1-α)/α * D .* eqn.x
end;

# Aproximace casove derivace
function ddt(f::Field{T}, Δt) where {T}
    n = length(f.values)
    Equation{T}(sparse(I/Δt,n,n), f.values, -f.values/Δt)
end;

# Laplaceuv operator
function Δ(μ::Array{Float64,1}, u::Field{T}) where {T}
    mesh = u.mesh

    A = spzeros(length(u.values),length(u.values))
    b = zeros(T,length(u.values))

    for f in internal_faces(mesh)
        o, n = mesh.owner[f], mesh.neighbor[f]
        xo, xn = mesh.centre[o], mesh.centre[n]
        S = mesh.surface[f]

        μf = (μ[o] + μ[n]) / 2
        g = μf * norm(S) / norm(xn - xo)

        A[o,o] -= g / mesh.volume[o]
        A[o,n] += g / mesh.volume[o]

        A[n,o] += g / mesh.volume[n]
        A[n,n] -= g / mesh.volume[n]
    end

    for p in boundary_patches(mesh)
        name = mesh.patch[p].name
        bc = u.boundaries[name]

        for f in patch_faces(mesh, p)
            o = mesh.owner[f]
            S = mesh.surface[f]
            xo, xf = mesh.centre[o], mesh.facecentre[f]
            c1, c2 = boundary_coeffs(bc, f)

            # normal derivative at boundary dudn = (ub - uin) / delta
            # with ub = c[1] + c[2]*uin, i.e.
            # dudn = c[1]/delta + (c[2]-1)/delta*uin

            δ = dot(xf - xo, S)/norm(S)
            μf = μ[o]

            A[o,o] += μf * (c2-1)/δ * norm(S) / mesh.volume[o]
            b[o] += μf * c1/δ * norm(S) / mesh.volume[o]
        end
    end

    return Equation{T}(A, u.values, b)

end

Δ(μ::Float64, u::Field{T}) where {T} = Δ(μ*ones(length(u.values)), u::Field{T})

function ∂(u::Field{T},dir) where {T}
    A = spzeros(length(u.values),length(u.values))
    b = zeros(T, length(u.values))

    mesh = u.mesh

    for f in internal_faces(mesh)
        o = mesh.owner[f]
        n = mesh.neighbor[f]
        uo = u[o]
        un = u[n]

        ϕ = mesh.surface[f][dir]
        α = ϕ / 2
        β = ϕ / 2

        A[o,o] += α / mesh.volume[o]
        A[o,n] += β / mesh.volume[o]

        A[n,o] -= α / mesh.volume[o]
        A[n,n] -= β / mesh.volume[o]
    end

    for p in boundary_patches(mesh)
        name = mesh.patch[p].name
        bc = u.boundaries[name]

        for f in patch_faces(mesh, p)
            o = mesh.owner[f]
            ϕ = mesh.surface[f][dir]
            α = ϕ / 2
            β = ϕ / 2

            c1, c2 = boundary_coeffs(bc, f)  # tj. ub = c[1] + c[2]*uin
            A[o,o] += (α + β*c2) / mesh.volume[o]
            b[o]   += β*c1 / mesh.volume[o]
        end
    end

    return Equation{T}(A, u.values, b)
end

∂x(u::Field{T}) where {T} = ∂(u::Field{T},1)
∂y(u::Field{T}) where {T} = ∂(u::Field{T},2)

# Explicitni vypocet divergence
function div(u::VectorField)
    mesh = u.mesh
    r = zeros(length(cells(mesh)))

    for f in internal_faces(mesh)
        o, n = mesh.owner[f], mesh.neighbor[f]
        uf = (u[o] + u[n]) / 2.0
        flux = dot(uf, mesh.surface[f])
        r[o] += flux
        r[n] -= flux
    end

    for p in boundary_patches(mesh)
        name = mesh.patch[p].name
        bc = u.boundaries[name]

        for f in patch_faces(mesh, p)
            ub = boundary_value(bc, f)
            flux = dot(ub, mesh.surface[f])
            r[mesh.owner[f]] += flux
        end
    end

    for c in cells(mesh)
        r[c] /= mesh.volume[c]
    end

    return r
end

# Explicitni vypocet gradientu
function grad(u::ScalarField)
    mesh = u.mesh
    r = zeros(Vector,length(cells(mesh)))

    for f in internal_faces(mesh)
        o, n = mesh.owner[f], mesh.neighbor[f]
        uf = (u[o] + u[n]) / 2.0
        flux = uf * mesh.surface[f]
        r[o] += flux
        r[n] -= flux
    end

    for p in boundary_patches(mesh)
        name = mesh.patch[p].name
        bc = u.boundaries[name]

        for f in patch_faces(mesh, p)
            o = mesh.owner[f]
            ub = boundary_value(bc, f)
            flux = ub * mesh.surface[f]
            r[o] += flux
        end
    end

    for c in cells(mesh)
        r[c] /= mesh.volume[c]
    end

    return r
end

max_stars_repo_name: furstj/NSP

max_issues_repo_name: furstj/NSP

max_forks_repo_name: furstj/NSP


text:
module BasisModule
export GaussianBasis, computeBasis, computeBasisShells
import Base.display
import ..IntegralsModule.normalize!

using ..BaseModule
using ..BasisFunctionsModule
using ..BasisSetModule
using ..GeometryModule
using ..AtomModule
using ..IntegralsModule
using ..ShellModule

abstract Basis

"""
A GaussianBasis is a Basis constructed purely from Gaussian type functions (see PrimitiveGaussianBasisFunction).
The list of basis functions defines the basis. It is most easily constructed from a BasisSet and a Geometry by computeBasis().
"""
type GaussianBasis <: Basis
  contractedBFs::Array{ContractedGaussianBasisFunction,1}
end

function computeBasis(basSet::BasisSet,geo::Geometry)
  bas=GaussianBasis([])
  for atom in geo.atoms
		for contractedDefinition in basSet.definitions[atom.element]
			for mqn in MQuantumNumbers(contractedDefinition.lQuantumNumber)
				contractedBF = ContractedGaussianBasisFunction([],[])
				for primitiveDefinition in contractedDefinition.primitives
					exponent=primitiveDefinition.exponent
					primitiveBF = PrimitiveGaussianBasisFunction(atom.position,exponent,mqn)
					norm = IntegralsModule.computeIntegralOverlap(primitiveBF,primitiveBF)
					push!(contractedBF.coefficients,primitiveDefinition.prefactor/sqrt(norm))
					push!(contractedBF.primitiveBFs,primitiveBF)
				end
				push!(bas.contractedBFs,contractedBF)
			end
    end
  end
  return bas
end

function normalize!(basis::GaussianBasis)
  normalize!(basis.contractedBFs)
end

function display(basis::GaussianBasis)
  dump(typeof(basis))
  indent = "  "

  println(indent * "contractedBFs: ")
  indent = indent * "  "

  for cgbf in basis.contractedBFs
    display(cgbf,indent)
  end
end

function computeBasisShells(basSet::BasisSet,geo::Geometry)
  shells = Shell[]
  for atom in geo.atoms
    for contractedDefinition in basSet.definitions[atom.element]
      exponents = [prim.exponent for prim in contractedDefinition.primitives]
      coefficients = [prim.prefactor for prim in contractedDefinition.primitives]
      sh = Shell(atom.position,contractedDefinition.lQuantumNumber,exponents,coefficients)
      push!(shells,sh)
    end
  end
  return shells
end

end # module


text:
# Load packages

using Vann
using PyPlot
using DataFrames
using JLD

# Settings

path_inputs = "C:/Work/Studies/vann/data/norway"

epot_choice = epot_monthly
snow_choice = TinBasic
hydro_choice = Gr4j

tstep = 24.0

calib_start = DateTime(2000,09,01)
calib_stop = DateTime(2014,12,31)

valid_start = DateTime(1985,09,01)
valid_stop = DateTime(2000,08,31)

# Select watersheds

dir_all = readdir(path_inputs)

dir_cur = dir_all[143]

print("Calibration for data in $(dir_cur)\n")

# Run for calibration period

# Load data

date, tair, prec, q_obs, frac = load_data("$path_inputs/$dir_cur")

# Crop data

date, tair, prec, q_obs = crop_data(date, tair, prec, q_obs, calib_start, calib_stop)

# Compute potential evapotranspiration

epot = eval(Expr(:call, epot_choice, date))

# Precipitation correction step

if ~all(isnan(q_obs))

  ikeep = ~isnan(q_obs)

  prec_tmp = mean(prec,1)

  prec_sum = sum(prec_tmp[ikeep])
  q_sum = sum(q_obs[ikeep])
  epot_sum = sum(epot)

  pcorr = (q_sum + 0.5 * epot_sum) / prec_sum

  prec = pcorr * prec

end

# Initilize model

st_snow = eval(Expr(:call, snow_choice, tstep, date[1], frac))
st_hydro = eval(Expr(:call, hydro_choice, tstep, date[1]))

# Run calibration

param_snow, param_hydro = run_model_calib(st_snow, st_hydro, tair, prec, epot, q_obs)

println(param_snow)
println(param_hydro)

# Reinitilize model

st_snow = eval(Expr(:call, snow_choice, tstep, date[1], param_snow, frac))
st_hydro = eval(Expr(:call, hydro_choice, tstep, date[1], param_hydro))

#=init_states(st_snow)
init_states(st_hydro)=#

# Run model with best parameter set

q_sim, st_snow, st_hydro = run_model(st_snow, st_hydro, tair, prec, epot; return_all = true)

# Plot results

plot_sim(st_hydro; q_obs = q_obs)

plot_sim(st_snow)

# Run for validation period

# Load data

date, tair, prec, q_obs, frac = load_data("$path_inputs/$dir_cur")

# Crop data

date, tair, prec, q_obs = crop_data(date, tair, prec, q_obs, valid_start, valid_stop)

# Compute potential evapotranspiration

epot = eval(Expr(:call, epot_choice, date))

# Reinitilize model

st_snow = eval(Expr(:call, snow_choice, tstep, date[1], param_snow, frac))
st_hydro = eval(Expr(:call, hydro_choice, tstep, date[1], param_hydro))

# Precipitation correction step

if isdefined(:pcorr)
  prec = pcorr * prec
end

# Run model with best parameter set

q_sim, st_snow, st_hydro = run_model(st_snow, st_hydro, tair, prec, epot; return_all = true)

# Plot results

plot_sim(st_hydro; q_obs = q_obs)

plot_sim(st_snow)

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
function example(name)
    fn = joinpath(dirname(@__FILE__),"..","..","examples",name*".jl")
    include(fn)
end


text:
# Regularizer

abstract Regularizer

## generic methods

value_and_grad{T<:AbstractFloat}(f::Regularizer, θ::StridedArray{T}) = value_and_addgrad!(f, zero(T), similar(θ), one(T), θ)

prox!(f::Regularizer, θ::StridedArray) = prox!(f, θ, θ)

prox(f::Regularizer, θ::StridedArray) = prox!(f, similar(θ), θ, 1.0)
prox(f::Regularizer, θ::StridedArray, λ::Real) = prox!(f, similar(θ), θ, λ)

function value_and_grad!{T<:AbstractFloat,N}(rm::RiskModel, reg::Regularizer, g::StridedArray{T,N}, θ::StridedArray{T,N}, X, y)
    v_risk, _ = value_and_addgrad!(rm, zero(T), g, one(T), θ, X, y)
    v_regr, _ = value_and_addgrad!(reg, one(T), g, one(T), θ)
    return (v_risk + v_regr, g)
end

## ZeroReg: 0

immutable ZeroReg <: Regularizer
end

value{T<:Real}(::ZeroReg, θ::StridedArray{T}) = zero(T)

value_and_addgrad!{T<:Real,N}(::ZeroReg, β::T, g::StridedArray{T,N}, α::T, θ::StridedArray{T,N}) =
    (zero(T), g)

prox!{T<:Real}(::ZeroReg, r::StridedArray{T}, θ::StridedArray{T}, λ::Real) =
    copy!(r, θ)

## SqrL2Reg: (c/2) * ||θ||_2^2

immutable SqrL2Reg{T<:AbstractFloat} <: Regularizer
    c::T
end

SqrL2Reg{T<:AbstractFloat}(c::T) = SqrL2Reg{T}(c)

value{T<:BlasReal}(f::SqrL2Reg{T}, θ::StridedArray{T}) = half(f.c * sumabs2(θ))

function value_and_addgrad!{T<:BlasReal,N}(f::SqrL2Reg{T}, β::T, g::StridedArray{T,N}, α::T, θ::StridedArray{T,N})
    @_checkdims length(g) == length(θ)
    c = f.c * α
    v = zero(T)
    if β == zero(T)
        @inbounds for i in eachindex(θ)
            θ_i = θ[i]
            v += abs2(θ_i)
            g[i] = c * θ_i
        end
    else
        @inbounds for i in eachindex(θ)
            θ_i = θ[i]
            v += abs2(θ_i)
            g[i] = β * g[i] + c * θ_i
        end
    end
    (half(c * v), g)
end

function prox!{T<:BlasReal}(f::SqrL2Reg{T}, r::StridedArray{T}, θ::StridedArray{T}, λ::Real)
    c = convert(T, λ) * f.c
    scale!(r, θ, one(T) / (one(T) + c))
end

## L1Reg: c * ||θ||_1

immutable L1Reg{T<:AbstractFloat} <: Regularizer
    c::T
end

L1Reg{T<:AbstractFloat}(c::T) = L1Reg{T}(c)

value{T<:BlasReal}(f::L1Reg{T}, θ::StridedArray{T}) = f.c * sumabs(θ)

function value_and_addgrad!{T<:BlasReal,N}(f::L1Reg{T}, β::T, g::StridedArray{T,N}, α::T, θ::StridedArray{T,N})
    @_checkdims length(g) == length(θ)
    c = f.c * α
    v = zero(T)
    if β == zero(T)
        for i in eachindex(θ)
            θ_i = θ[i]
            v += abs(θ_i)
            @inbounds g[i] = c * sign(θ_i)
        end
    else
        for i in eachindex(θ)
            θ_i = θ[i]
            v += abs(θ_i)
            @inbounds g[i] = β * g[i] + c * sign(θ_i)
        end
    end
    (c * v, g)
end

function prox!{T<:BlasReal}(f::L1Reg{T}, r::StridedArray{T}, θ::StridedArray{T}, λ::Real)
    @_checkdims size(r) == size(θ)
    c = convert(T, λ) * f.c
    for i in eachindex(θ)
        r[i] = shrink(θ[i], c)
    end
    r
end

## ElasticNet: c1 * ||θ||_1 + c2 * ||θ||_2^2

immutable ElasticReg{T<:AbstractFloat} <: Regularizer
    c1::T
    c2::T
end

ElasticReg{T<:AbstractFloat}(c1::T, c2::T) = ElasticReg{T}(c1, c2)

function value{T<:BlasReal}(f::ElasticReg{T}, θ::StridedArray{T})
    s = zero(T)
    c1 = f.c1
    c2_h = half(f.c2)
    @inbounds for i in eachindex(θ)
        θ_i = θ[i]
        s += c1 * abs(θ_i) + c2_h * abs2(θ_i)
    end
    s
end

function value_and_addgrad!{T<:BlasReal,N}(f::ElasticReg{T}, β::T, g::StridedArray{T,N}, α::T, θ::StridedArray{T,N})
    @_checkdims length(g) == length(θ)
    c1 = f.c1 * α
    c2 = f.c2 * α

    v1 = zero(T)
    v2 = zero(T)
    if β == zero(T)
        @inbounds for i in eachindex(θ)
            θ_i = θ[i]
            v1 += abs(θ_i)
            v2 += abs2(θ_i)
            g[i] = c1 * sign(θ_i) + c2 * θ_i
        end
    else
        @inbounds for i in eachindex(θ)
            θ_i = θ[i]
            v1 += abs(θ_i)
            v2 += abs2(θ_i)
            g[i] = β * g[i] + (c1 * sign(θ_i) + c2 * θ_i)
        end
    end
    (c1 * v1 + half(c2 * v2), g)
end

function prox!{T<:BlasReal}(f::ElasticReg{T}, r::StridedArray{T}, θ::StridedArray{T}, λ::Real)
    @_checkdims size(r) == size(θ)
    c1 = convert(T, λ) * f.c1
    c2 = convert(T, λ) * f.c2
    c = one(T) / (one(T) + c2)
    t = c1 / (one(T) + c2)
    @inbounds for i in eachindex(θ)
        r[i] = shrink(c * θ[i], t)
    end
    r
end

## NonNegReg: θ[i] >= 0 ∀ i
immutable NonNegReg <: Regularizer
end

function value{T<:Real}(::NonNegReg, θ::StridedArray{T})
    s = zero(T)
    @inbounds for i in eachindex(θ)
        θ[i] < s && return Inf      # value = + Inf if outside domain
    end
    s   # value = 0. otherwise
end

value_and_addgrad!{T<:Real,N}(::NonNegReg, β::T, g::StridedArray{T,N}, α::T, θ::StridedArray{T,N}) =
    error("Gradient undefined for the NonNegReg regularizer")

#  NOTE : λ is ignored
function prox!{T<:Real}(::NonNegReg, r::StridedArray{T}, θ::StridedArray{T}, λ::Real)
    @_checkdims size(r) == size(θ)
    c = convert(T, λ)
    s = zero(T)
    for i in eachindex(θ)
        r[i] = θ[i] < s ? s : θ[i]
    end
    r
end

## Simplex: Σ θ[i] = c, θ[i] >= 0 ∀ i
immutable SimplexReg{T<:AbstractFloat} <: Regularizer
    c::T
end

SimplexReg{T<:AbstractFloat}(c::T=1.0) = SimplexReg{T}(c)

function value{T<:BlasReal}(f::SimplexReg{T}, θ::StridedArray{T})
    z = zero(T)
    s = zero(T)
    @inbounds for i in eachindex(θ)
        θ[i] < z && return Inf       # value = + Inf if outside domain
        s += θ[i]
    end

    return abs(s - f.c) <= 1.0e-12 ? z : Inf
end

value_and_addgrad!{T<:BlasReal,N}(f::SimplexReg{T}, β::T, g::StridedArray{T,N}, α::T, θ::StridedArray{T,N}) =
    error("Gradient undefined for the SimplexReg regularizer")

# Source = Efficient Projections on the l1-Ball for Learning in High Dimensions - J.Duchi
# Note : λ is ignored
function prox!{T<:BlasReal}(f::SimplexReg{T}, r::StridedArray{T}, θ::StridedArray{T}, λ::Real)
    @_checkdims size(r) == size(θ)

    z  = zero(T)
    ρ, s = 0, z
    U  = [1:length(θ);]
    cU = length(θ)

    # pivot search by divide and conquer
    while cU > 0
        k  = U[rand(1:cU)]  # pick an element at random
        vk = θ[k]
        dρ, ds = 0, z
        @inbounds for i in 1:cU
            j = U[i]
            θ[j] >= vk || continue
            dρ += 1
            ds += θ[j]
        end

        if (s+ds) - (ρ+dρ)*vk < f.c
            s += ds
            ρ += dρ
            nCu = 0
            @inbounds for i in 1:cU
                j = U[i]
                θ[j] < vk || continue
                nCu += 1
                U[nCu] = U[i]
            end
            cU = nCu
        else
            nCu = 0
            @inbounds for i in 1:cU
                j = U[i]
                θ[j] >= vk || continue
                j != k     || continue
                nCu += 1
                U[nCu] = U[i]
            end
            cU = nCu
        end
    end

    # apply pivot
    @inbounds for i in eachindex(θ)
        r[i]  = max( θ[i] - (s - f.c)/ρ, z )
    end
    r

    r
end

## L1Ball: Σ |θ[i]| = c
immutable L1BallReg{T<:AbstractFloat} <: Regularizer
    c::T
end

L1BallReg{T<:AbstractFloat}(c::T) = L1BallReg{T}(c)

function value{T<:BlasReal}(f::L1BallReg{T}, θ::StridedArray{T})
    z = zero(T)
    s = zero(T)
    @inbounds for i in eachindex(θ)
        s += abs(θ[i])
        s > f.c && return Inf
    end

    return z
end

value_and_addgrad!{T<:BlasReal,N}(f::L1BallReg{T}, β::T, g::StridedArray{T,N}, α::T, θ::StridedArray{T,N}) =
    error("Gradient undefined for the L1BallReg regularizer")

# Source = Efficient Projections on the l1-Ball for Learning in High Dimensions - J.Duchi
# Note : λ is ignored
function prox!{T<:BlasReal}(f::L1BallReg{T}, r::StridedArray{T}, θ::StridedArray{T}, λ::Real)
    @_checkdims size(r) == size(θ)

    z  = zero(T)
    ρ, s = 0, z
    U  = [1:length(θ);]
    cU = length(θ)

    # early exit if within limits
    if sumabs(θ) <= f.c
        copy!(r, θ)
        return r
    end

    # pivot search by divide and conquer
    while cU > 0
        k  = U[rand(1:cU)]
        vk = abs(θ[k])
        dρ, ds = 0, z
        @inbounds for i in 1:cU
            j = U[i]
            abs(θ[j]) >= vk || continue
            dρ += 1
            ds += abs(θ[j])
        end

        if (s+ds) - (ρ+dρ)*vk < f.c
            s += ds
            ρ += dρ
            nCu = 0
            @inbounds for i in 1:cU
                j = U[i]
                abs(θ[j]) < vk || continue
                nCu += 1
                U[nCu] = U[i]
            end
            cU = nCu
        else
            nCu = 0
            @inbounds for i in 1:cU
                j = U[i]
                abs(θ[j]) >= vk || continue
                j != k || continue
                nCu += 1
                U[nCu] = U[i]
            end
            cU = nCu
        end
    end

    # apply pivot
    @inbounds for i in eachindex(θ)
        r[i]  = sign(θ[i]) * max( abs(θ[i]) - (s - f.c)/ρ, z )
    end
    r
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# Auxiliary, non-exported functions are declared here.

# convert marker string description to gnuplot's expected number
function pointtype(x::String)
    x == "+" && return "1"
    x == "x" && return "2"
    x == "*" && return "3"
    x == "esquare" && return "4"
    x == "fsquare" && return "5"
    x == "ecircle" && return "6"
    x == "fcircle" && return "7"
    x == "etrianup" && return "8"
    x == "ftrianup" && return "9"
    x == "etriandn" && return "10"
    x == "ftriandn" && return "11"
    x == "edmd" && return "12"
    x == "fdmd" && return "13"
    return "1"
end

# create a Z-coordinate matrix from x, y coordinates and a function
function meshgrid(x,y,f)
    Z = zeros(length(x),length(y))
    for k = 1:length(x)
        Z[k,:] = [ f(i,j) for i=x[k], j=y ]
    end
    return Z
end

# create x,y coordinates for a histogram, from a sample vector, using a number
# of bins
function hist(s,bins)
    # When adding an element s to a bin, we use an iequality m < s <= M.
    # In order to account for elements s==m, we need to special-case
    # the computation for the first bin
    ms = minimum(s)
    Ms = maximum(s)
    bins = max(bins, 1)
    if Ms == ms
        # compute a "natural" scale
        g = (10.0^floor(log10(abs(ms)+eps()))) / 2
        ms, Ms = ms - g, ms + g
    end
    delta = (Ms-ms)/bins
    x = ms:delta:Ms
    y = zeros(bins)
    # this is special-cased because we want to include the minimum in the
    # first bin
    y[1] = sum(ms .<= s .<= x[2])
    for i in 2:length(x)-2
        y[i] = sum(x[i] .< s .<= x[i+1])
    end
    # this is special-cased because there is no guarantee that x[end] == Ms
    # (because of how ranges work)
    if length(y) > 1 y[end] = sum(x[end-1] .< s .<= Ms) end

    if bins != 1
        # We want the left bin to start at ms and the right bin to end at Ms
        x = (ms+delta/2):delta:Ms
    else
        # add two empty bins on the sides to provide a scale to gnuplot
        x = (ms-delta/2):delta:(ms+delta/2)
        y = [0.0, y[1], 0.0]
    end
    return x,y
end

function Base.show(io::IO, ::MIME"text/plain", x::Figure)
    isempty(x) && return nothing
    IsJupyterOrJuno && return nothing
    llplot(x)
    if config[:mode] != "null"
        if (config[:term][:terminal] ∈ term_text) || (config[:mode] == "ijulia")
            write(io, x.svg)
        end
    end
    return nothing
end

function Base.show(io::IO, ::MIME"image/svg+xml", x::Figure)
    llplot(x)
    write(io,x.svg)
    return nothing
end

# return configuration string for a single plot
function linestr_single(conf::CurveConf)
    s = ""
    conf.legend != "" && (s *= " title '"*conf.legend*"' ")
    conf.plotstyle != "" && (s *= " with "*conf.plotstyle*" ")
    conf.linecolor != "" && (s *= "lc rgb '"*conf.linecolor*"' ")
    conf.linewidth != "" && (s *= "lw "*conf.linewidth*" ")
    conf.linestyle != "" && (s *= "dt '"*conf.linestyle*"' ")
    conf.fillstyle != "" && (s *= "fs "*conf.fillstyle*" ")
    conf.fillcolor != "" && (s *= "fc \""*conf.fillcolor*"\" ")
    # some plotstyles don't allow point specifiers
    if conf.plotstyle ∈ ps_sup_points
        if conf.pointtype != ""
            if conf.pointtype ∈ supported_pointtypes
                s = s*"pt "*pointtype(conf.pointtype)*" "
            else
                s = s*"pt \""*conf.pointtype*"\" "
            end
            conf.pointsize != "" && (s = s*"ps "*conf.pointsize*" ")
        end
    end
    return s
end

# build a string with plot commands according to configuration
function linestr(curves::Vector{Curve}, cmd, file)
    # We have to insert "," between plot commands. One easy way to do this
    # is create the first plot command, then the rest
    # We also need to keep track of the current index (starts at zero)
    index = 0
    s = cmd*" '"*file*"' "*" i 0 "*linestr_single(curves[1].conf)
    if length(curves) > 1
        for i in curves[2:end]
            index += 1
            s = s*", '"*file*"' "*" i "*string(index)*" "*linestr_single(i.conf)
        end
    end
    return s
end

# Build a "set term" string appropriate for the terminal type
function termstring(f::Figure,print=false)
    global gnuplot_state

    ac = f.axes
    tc = f.term
    pc = f.print

    term = print ? pc.print_term : config[:term][:terminal]

    if term != ""
        # determine font, size, and background
        font = print ? pc.print_font : tc.font
        size = print ? pc.print_size : tc.size
        background = print ? pc.print_background : tc.background

        # build term string
        ts = "set term $term "
        term ∈ term_window && (ts = ts*string(gnuplot_state.current)*" ")

        isempty(font) ? s = "" : s = " font \""*font*"\" "
        term ∈ term_sup_font && (ts *= s)

        isempty(size) ? s = "" : s = " size "*size*" "
        term ∈ term_sup_size && (ts *= s)

        isempty(background) ? s = "" : s = " background \""*background*"\" "
        term ∈ term_sup_bkgnd && (ts *= s)

        # terminal options
        print || (ts *= config[:term][:termopts]*" ")
        print && (ts *= config[:print][:print_termopts]*" ")

        # set output file
        if term ∈ term_file
            s = ""
            isempty(ac.output) || (s = "\nset output \"$(ac.output)\" ")
            ts = ts*s
        end
    end
    return ts
end

# send gnuplot the current figure's configuration
function gnuplot_send_fig_config(config)
    config.title != "" && gnuplot_send("set title '"*config.title*"' ")
    config.fillstyle != "" && gnuplot_send("set style fill "*config.fillstyle)
    if config.grid != ""
        if config.grid == "on"
            gnuplot_send("set grid")
        else
            gnuplot_send("set grid "*config.grid)
        end
    end
    config.keyoptions != "" && gnuplot_send("set key "*config.keyoptions)
    config.boxwidth != "" && gnuplot_send("set boxwidth "*config.boxwidth)
    if config.axis != ""
        config.axis == "semilogx" && gnuplot_send("set logscale x")
        config.axis == "semilogy" && gnuplot_send("set logscale y")
        config.axis == "semilogz" && gnuplot_send("set logscale z")
        config.axis == "loglog" && gnuplot_send("set logscale xyz")
    end
    config.xlabel != "" && gnuplot_send("set xlabel '"*config.xlabel*"' ")
    config.ylabel != "" && gnuplot_send("set ylabel '"*config.ylabel*"' ")
    config.zlabel != "" && gnuplot_send("set zlabel '"*config.zlabel*"' ")
    config.xrange != "" && gnuplot_send("set xrange "*config.xrange)
    config.yrange != "" && gnuplot_send("set yrange "*config.yrange)
    config.zrange != "" && gnuplot_send("set zrange "*config.zrange)
    if config.xzeroaxis != ""
        if config.xzeroaxis == "on"
            gnuplot_send("set xzeroaxis")
        else
            gnuplot_send("set xzeroaxis "*config.xzeroaxis)
        end
    end
    if config.yzeroaxis != ""
        if config.yzeroaxis == "on"
            gnuplot_send("set yzeroaxis")
        else
            gnuplot_send("set yzeroaxis "*config.yzeroaxis)
        end
    end
    if config.zzeroaxis != ""
        if config.zzeroaxis == "on"
            gnuplot_send("set zzeroaxis")
        else
            gnuplot_send("set zzeroaxis "*config.zzeroaxis)
        end
    end
    config.palette != "" && gnuplot_send("set palette "*config.palette)
end

# write commands to gnuplot's pipe
function gnuplot_send(s)
    config[:debug] && println(s)  # print gnuplot commands if debug enabled
    w = write(P.gstdin, string(s,"\n"))
    # check that data was accepted by the pipe
    if !(w > 0)
        println("Something went wrong writing to gnuplot STDIN.")
        return
    end
    flush(P.gstdin)
end


text:
#simple frequency domain FWI to illustrate PARSDMM
#include("/home/slim/bpeters/.julia/v0.6/WAVEFORM/src/Waveform.jl")
using Waveform

#This example used the Waveform package, which is in julia 0.6 and uses the v06 branch of JOLI (https://github.com/slimgroup/JOLI.jl)
#untill everyting is updated to current julia version, you need to got to (on linux) /.julia/v0.6/JOLI and use: git checkout v06

using JOLI
ENV["MPLBACKEND"]="qt5agg"
using PyPlot
#using OptimPackNextGen
#using opesciSLIM.SLIM_optim
#using Waveform

@everywhere using SetIntersectionProjection

#include("SPGslim.jl")
include(string(Pkg.dir("JUDI"),"/src/Optimization/SPGSlim.jl"))
include(string(Pkg.dir("JUDI"),"/src/Optimization/OptimizationFunctions.jl"))

#data directory for loading and writing results
#save_dir = "/data/slim/bpeters/SetIntersection_data_results"
save_dir = pwd()#"/Volumes/Users/bpeters/Downloads"

# Set up the model geometry
# In 2D, we organize our coordinates as (z,x) for legacy + visualization purposes

# Number of gridpoints in z-, x-
n = [151;101];

# Spacing of gridpoints in z-, x-
d = 10.0*[1;1];

# Origin point of domain in z-, x- (typically (0.,0.))
o = 0.0*[1;1];

# Constant velocity background
vel_background = 2500;

# Source wavelet time shift (in seconds)
t0 = 0.0;

# Peak frequency of Ricker Wavelet (in Hz)
f0 = 10.0;

# Units of the model
unit = "m/s";

# Model domain range
L = n.*d;

# x-coordinates of the sources
xsrc = 0.0:50.0:L[2];

# y-coordinates of the sources (irrelevant in 2D)
ysrc = [0.0];

# z-coordinates of the sources
zsrc = [10.0];

# x-coordinates of the receivers
xrec = [950.0];

# y-coordinates of the receivers (irrelevant in 2D)
yrec = [0.0];

# z-coordinates of the receivers (this is a transmission experiment)
#zrec = 0.0:50.0:L[2];
zrec = 0.0:50.0:L[1];

# Frequencies to compute
freqs = [3.0 ; 5.0 ; 7.0]# ; 8.0]#[5.0;10.0;15.0];

# Maximum wavelength
λ = vel_background/maximum(freqs);

# Model type
model = Model{Int64,Float64}(n,d,o,t0,f0,unit,freqs,xsrc,ysrc,zsrc,xrec,yrec,zrec);
(z,x) = odn_to_grid(o,d,n);
nsrc,nrec,nfreq = length(xsrc),length(xrec),length(freqs)

# Set up computational parameters
## Size/spacing of computational domain
comp_n = n;
comp_d = d;
comp_o = o;

# Number of PML points
npml = round(Int,λ/minimum(comp_d))*[1 1; 1 1];

# PDE scheme
scheme = Waveform.helm2d_chen9p;

# Whether to remove elements in the PML when restricting to the model domain (if false, will use stacking)
cut_pml = true;

# If true, use an implicit matrix representation (false for 2D, the matrices are small enough to form explicitly)
implicit_matrix = false;

# Binary mask of the sources and frequencies to compute
srcfreqmask = trues(nsrc,nfreq);

# Misfit function for the objective
misfit = least_squares;

# Linear solver options
lsopts = LinSolveOpts(solver=:lufact);

# Computational options type
opts = PDEopts{Int64,Float64}(scheme,comp_n,comp_d,comp_o,cut_pml,implicit_matrix,npml,misfit,srcfreqmask,lsopts);

# Create velocity model
v0 = vel_background*ones(n...);
v = copy(v0);
#v[div(n[1],3):2*div(n[1],3),div(n[2],3):2*div(n[2],3)] = 1.25*vel_background;
v[61:80,50:70]=2400
v = vec(v);
v0 = vec(v0);

function plot_velocity(vs,title_str,model,keyword,rec_x,rec_z,src_x,src_z)
    fig = figure(figsize=(2.1, 1.6))
    FS  = 4
    LS  = 3
    TML = 1
    PD  = 1
    #ax1 = axes()
    #vplot = imshow(reshape(vs,model.n...),vmin=minimum(v)-50,vmax=maximum(v)+50,axes=ax1,cmap="jet",extent=[0.0 , model.d[2]*model.n[2] , model.d[1]*model.n[1] , 0.0])
    vplot = imshow(reshape(vs,model.n...),vmin=minimum(v)-50,vmax=maximum(v)+50,cmap="jet",extent=[0.0 , model.d[2]*model.n[2] , model.d[1]*model.n[1] , 0.0])
    tick_params(labelsize=LS,length=TML,pad=PD)
    title(title_str,FontSize=FS)
    xlabel("x [m]",FontSize=FS)
    ylabel("z [m]",FontSize=FS)
    zt = 0:20:n[1]
    xt = 0:20:n[2]
    #xticks(xt,round.(Int,xt*model.d[2]))
    #yticks(zt,round.(Int,zt*model.d[1]))
    cbar = colorbar()
    cbar[:ax][:tick_params](labelsize=LS)
    figure;plot(rec_x,rec_z,linewidth=1.0, marker="o",linestyle="",markersize=2)
    figure;plot(src_x,src_z,linewidth=1.0, marker="x",linestyle="","k",markersize=2)
    savefig(joinpath(save_dir,string("CFWI_simple_freq_m_est_",keyword,".eps")),bbox_inches="tight",dpi=300)
    return nothing
end
plot_velocity(v,"a) True velocity",model,"true",repmat(xrec,length(zrec),1),zrec,xsrc,repmat(zsrc,length(xsrc),1))
plot_velocity(v0,"b) Initial velocity",model,"initial",repmat(xrec,length(zrec),1),zrec,xsrc,repmat(zsrc,length(xsrc),1))

# Source weight matrix
Q = eye(nsrc);

# Generate data
D = forw_model(v,Q,model,opts);

imshow(real(D[:,1:nsrc]),aspect="auto");
title("Frequency slice at $(model.freq[1]) Hz")

# Frequency continuation
max_func_evals = 10

# Frequency continuation #original
# size_freq_batch = 2
# overlap = 1
# freq_partition = partition(nfreq,size_freq_batch,overlap)
# vest = v0
# opts1 = deepcopy(opts);
# proj! = (xproj,x)->project_bounds!(x,minimum(v),maximum(v),xproj);
#
# for j in 1:size(freq_partition,2)
#     fbatch = freq_partition[:,j]
#     srcfreqmask = falses(nsrc,nfreq)
#     srcfreqmask[:,fbatch] = true
#     opts1.srcfreqmask = srcfreqmask
#     obj! = construct_pde_misfit(v,Q,D,model,opts1,batch_mode=false)
#     vest = spg(obj!,proj!,vest,3,maxfc=max_func_evals)
# end

#Frequency continuation my own
size_freq_batch = 3
overlap = 1
freq_partition = partition(nfreq,size_freq_batch,overlap)
vest = v0
opts1 = deepcopy(opts);
# proj! = (xproj,x)->project_bounds!(x,minimum(v)-50,maximum(v)+50,xproj);
#
#
# for j in 1:size(freq_partition,2)
#     fbatch = freq_partition[:,j]
#     srcfreqmask = falses(nsrc,nfreq)
#     srcfreqmask[:,fbatch] = true
#     opts1.srcfreqmask = srcfreqmask
#     obj! = construct_pde_misfit(v,Q,D,model,opts1,batch_mode=false)
#     vest = spg(obj!,proj!,vest,3,maxfc=max_func_evals)
# end
#
# plot_velocity(vest,"Inverted velocity - freq continuation",model)

# with constraints:
options=PARSDMM_options()
options.FL=Float64
options=default_PARSDMM_options(options,options.FL)
options.adjust_gamma           = true
options.adjust_rho             = true
options.adjust_feasibility_rho = true
options.Blas_active            = true
options.maxit                  = 1000
options.feas_tol= 0.001
options.obj_tol=0.001
options.evol_rel_tol = 0.00001

options.rho_ini=[1.0f0]

set_zero_subnormals(true)
BLAS.set_num_threads(2)
FFTW.set_num_threads(2)
options.parallel=false
options.feasibility_only = false
options.zero_ini_guess=true

type compgrid
  d :: Tuple
  n :: Tuple
end
comp_grid=compgrid((model.d[1],model.d[2]),(model.n[1],model.n[2]))
#function run_CFWI(freq_partition,nsrc,nfreq,v,Q,D,model,opts1,proj!,vest,max_func_evals)

options_SPG = spg_options(verbose=3, maxIter=max_func_evals, memory=5,suffDec=1f-6)
options_SPG.progTol=1e-6
options_SPG.testOpt=false
options_SPG.interp=2

function run_CFWI(freq_partition,nsrc,nfreq,v,Q,D,model,opts1,options_SPG,proj!,vest,max_func_evals)
  for j in 1:size(freq_partition,2)
      fbatch = freq_partition[:,j]
      srcfreqmask = falses(nsrc,nfreq)
      srcfreqmask[:,fbatch] = true
      opts1.srcfreqmask = srcfreqmask
      obj! = construct_pde_misfit(v,Q,D,model,opts1,batch_mode=false)
      function objective_function(in,obj!)
        f=0.0
        g=zeros(length(in))
        f=obj!(in,g);
        return f,g
      end
      f_obj = in -> objective_function(in,obj!)
      vest, fsave, funEvals= minConf_SPG(f_obj, vec(vest), proj!, options_SPG)
      #vest = spg(obj!,proj!,vest,3,maxfc=max_func_evals,eps1 = 1e-8,eps2 = 1e-8)
  end
  return vest
end

CFWI(ini_model,projector) = run_CFWI(freq_partition,nsrc,nfreq,v,Q,D,model,opts1,options_SPG,projector,ini_model,max_func_evals)

constraint_strategy_list=[1 2 3 4 5]#
for j in constraint_strategy_list
  if j==1
    keyword="bounds_only"
    title_str="c) bounds only"

    constraint = Vector{SetIntersectionProjection.set_definitions}()

    #bounds:
    m_min     = minimum(v)-50.0
    m_max     = maximum(v)+50.0
    set_type  = "bounds"
    TD_OP     = "identity"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

elseif j==2 #various types of cardinality and rank
		keyword="cardmat_cardcol_rank_bounds"
    title_str="c)"# fiber, matrix grad. card. & rank & bounds"

    constraint = Vector{SetIntersectionProjection.set_definitions}()

    #bounds:
    m_min     = minimum(v)-50.0
    m_max     = maximum(v)+50.0
    set_type  = "bounds"
    TD_OP     = "identity"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

		#cardinality on derivatives (column and row wise)
    m_min     = 0
    m_max     = 2
    set_type  = "cardinality"
    TD_OP     = "D_x"
    app_mode  = ("fiber","x")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
		# constraint["use_TD_card_fiber_x"]			= true
		# constraint["card_fiber_x"] 						= 2
	  # constraint["TD_card_fiber_x_operator"]= "D_x"

    m_min     = 0
    m_max     = 2
    set_type  = "cardinality"
    TD_OP     = "D_z"
    app_mode  = ("fiber","z")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

		#cardinality on derivatives (matrix based)
    m_min     = 0
    m_max     = round(Integer,3*0.33*n[2])
    set_type  = "cardinality"
    TD_OP     = "D_x"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
		# constraint["use_TD_card_1"]=true
		# constraint["card_1"]=round(Integer,3*0.33*n[2])
		# constraint["TD_card_operator_1"]="D_x"

    m_min     = 0
    m_max     = round(Integer,3*0.33*n[1])
    set_type  = "cardinality"
    TD_OP     = "D_z"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
		# constraint["use_TD_card_2"]=true
		# constraint["card_2"]=round(Integer,3*0.33*n[1])
		# constraint["TD_card_operator_2"]="D_z"

		#rank constraint
    m_min     = 0
    m_max     = 3
    set_type  = "rank"
    TD_OP     = "identity"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))
		# constraint["use_TD_rank_1"]=true
	  # constraint["TD_max_rank_1"]=3
    # constraint["TD_rank_operator_1"]="identity"

  elseif j==3 #true tv and bounds
      keyword="trueTV_bounds"
      title_str="d) true TV & bounds"

      constraint = Vector{SetIntersectionProjection.set_definitions}()

      #bounds:
      m_min     = minimum(v)-50.0
      m_max     = maximum(v)+50.0
      set_type  = "bounds"
      TD_OP     = "identity"
      app_mode  = ("matrix","")
      custom_TD_OP = ([],false)
      push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

      #TV
      (TV,dummy1,dummy2,dummy3)=get_TD_operator(comp_grid,"TV",options.FL)
      m_min     = 0.0
      m_max     = norm(TV*vec(v),1)
      set_type  = "l1"
      TD_OP     = "TV"
      app_mode  = ("matrix","")
      custom_TD_OP = ([],false)
      push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    elseif j==4 # various types of cardinality
    			keyword="cardmat_cardcol_bounds"
          title_str="b)"#" fiber, matrix grad. card. & bounds"

    			constraint = Vector{SetIntersectionProjection.set_definitions}()

          #bounds:
          m_min     = minimum(v)-50.0
          m_max     = maximum(v)+50.0
          set_type  = "bounds"
          TD_OP     = "identity"
          app_mode  = ("matrix","")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    			#cardinality on derivatives (column and row wise)
          m_min     = 0
          m_max     = 2
          set_type  = "cardinality"
          TD_OP     = "D_x"
          app_mode  = ("fiber","x")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

          m_min     = 0
          m_max     = 2
          set_type  = "cardinality"
          TD_OP     = "D_z"
          app_mode  = ("fiber","z")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    			#cardinality on derivatives (matrix based)
          m_min     = 0
          m_max     = round(Integer,3*0.33*n[2])
          set_type  = "cardinality"
          TD_OP     = "D_x"
          app_mode  = ("matrix","")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

          m_min     = 0
          m_max     = round(Integer,3*0.33*n[1])
          set_type  = "cardinality"
          TD_OP     = "D_z"
          app_mode  = ("matrix","")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

        elseif j==5
    				keyword="cardmat_bounds"
            title_str="a)"# matrix grad. card. & bounds & rank"

    				constraint = Vector{SetIntersectionProjection.set_definitions}()

            #bounds:
            m_min     = minimum(v)-50.0
            m_max     = maximum(v)+50.0
            set_type  = "bounds"
            TD_OP     = "identity"
            app_mode  = ("matrix","")
            custom_TD_OP = ([],false)
            push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    				#cardinality on derivatives (matrix based)
            m_min     = 0
            m_max     = round(Integer,3*0.33*n[2])
            set_type  = "cardinality"
            TD_OP     = "D_x"
            app_mode  = ("matrix","")
            custom_TD_OP = ([],false)
            push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

            m_min     = 0
            m_max     = round(Integer,3*0.33*n[1])
            set_type  = "cardinality"
            TD_OP     = "D_z"
            app_mode  = ("matrix","")
            custom_TD_OP = ([],false)
            push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

            #rank constraint
            m_min     = 0
            m_max     = 3
            set_type  = "rank"
            TD_OP     = "identity"
            app_mode  = ("matrix","")
            custom_TD_OP = ([],false)
            push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

end #end if block for constraint setup choices

  #set up constraints, precompute some things and define projector
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)
  options.rho_ini        = ones(length(TD_OP))*10.0
  # for i=1:length(options.rho_ini)
  #   if set_Prop.ncvx[i]==true
  #     options.rho_ini[i]=1.0
  #   end
  # end
  proj_intersection = x-> PARSDMM(x,AtA,TD_OP,set_Prop,P_sub,comp_grid,options)
  # function prj!{T}(dst::T, src::T)
  #   (dst[:],dummy1,dummy2,dymmy3) = proj_intersection(src)
  #   return dst
  # end
   function prj!(input)
     (x,dummy1,dummy2,dymmy3) = proj_intersection(input)
     return x
   end
  v_ini=deepcopy(v0)
  x=CFWI(v_ini,prj!)
  plot_velocity(x,title_str,model,keyword,repmat(xrec,length(zrec),1),zrec,xsrc,repmat(zsrc,length(xsrc),1))

end #end loop over FWI with various constraint combinations

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module FLoopsBase

using ContextVariablesX
@contextvar EXTRA_STATE_VARIABLES::Union{Vector{Union{Symbol,Expr}},Nothing} = nothing

function with_extra_state_variables(f, variables)
    current_vars = copy(something(EXTRA_STATE_VARIABLES[], Union{Symbol,Expr}[]))
    new_vars = append!(current_vars, variables)
    with_context(f, EXTRA_STATE_VARIABLES => new_vars)
end

abstract type AbstractScratchSpace end

end


text:
## From Juno

import Base: iterate

# An iterator for the parse function: parsit(source) will iterate over the
# expressiosn in a string.
mutable struct ParseIt
    value::AbstractString
end

function parseit(value::AbstractString)
    ParseIt(value)
end

function Base.iterate(it::ParseIt)
    pos = 1
    iterate(it, pos)
#    (ex,newpos) = Meta.parse(it.value, 1)
#    ((it.value[1:(newpos-1)], ex), newpos)
end

function Base.iterate(it::ParseIt, pos)
    if pos > length(it.value)
        nothing
    else
        (ex,newpos) = Meta.parse(it.value, pos)
        ((it.value[pos:(newpos-1)], ex), newpos)
    end
end

# function start(it::ParseIt)
#     1
# end

# function next(it::ParseIt, pos)
#     (ex,newpos) = Meta.parse(it.value, pos)
#     ((it.value[pos:(newpos-1)], ex), newpos)
# end

# function done(it::ParseIt, pos)
#     pos > length(it.value)
# end

# A special dummy module in which a documents code is executed.
module WeaveSandbox
end

"""

Make a module of a given name.

"""
function make_module(nm=randstring())
    nm = "Z"*uppercase(nm)
    eval(Meta.parse("module " * nm * " end"))
    eval(Meta.parse(nm))
end

struct DisplayError
    x
end
Base.show(io::IO, ::MIME"text/plain", e::DisplayError) = println(io, e.x)

# Evaluate an expression and return its result and a string.
safeeval(m, ex::Nothing) = nothing
function safeeval(m, ex::Union{Number,Symbol, Expr})
    try
        res = Core.eval(m, ex)

    catch e
        printstyled("Error with evaluating $ex: $(string(e))\n", color=:red)
        DisplayError(string(e))
    end
end

function process_block(text, m = WeaveSandbox)
    result = ""
    for (cmd, ex) in parseit(strip(text))
        result = safeeval(m, ex)
    end
    result
end


text:
function test_battery()

    T = 2
    bat = DR.DER.Battery(
        index=0, T=T,
        soc_min=zeros(T), soc_max=ones(T),
        charge_pwr_max=1.0,discharge_pwr_max=1.0,
        charge_eff=1.0,
        discharge_eff=1.0
    )

    # GLPK does not support -in-Interval constraints, so we bridge
    model = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)

    lmin = [-2.0, -2.0]
    lmax = [ 2.0,  2.0]
    price = [1.0, 2.0]

    # Instantiate initial model
    h = DR.DER.House(
        0, T,
        lmin, lmax, price,
        [bat]
    )

    model, var2idx, con2idx = DR.DER.build_model!(h, model)

    MOI.optimize!(model)

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ -1.0

    # Check net load
    # Battery should charge then discharge
    pnet1 = var2idx[(:house, h.index, :pnet, 1)]
    pnet2 = var2idx[(:house, h.index, :pnet, 2)]
    @test MOI.get(model, MOI.VariablePrimal(), pnet1) ≈ 1.0
    @test MOI.get(model, MOI.VariablePrimal(), pnet2) ≈ -1.0

    return nothing
end

@testset "Battery" begin
    test_battery()
end


text:
export ConcatCube, concatenateCubes
export mergeAxes
import ..Cubes: YAXArray, caxes, iscompressed, cubechunks, chunkoffset
using DiskArrayTools: diskstack, DiskArrayTools

function Base.map(op, incubes::YAXArray...)
  axlist=copy(caxes(incubes[1]))
  all(i->caxes(i)==axlist,incubes) || error("All axes must match")
  props=merge(getattributes.(incubes)...)
  YAXArray(axlist,broadcast(op,map(c->c.data,incubes)...),props,mapreduce(i->i.cleaner,append!,incubes))
end

"""
    function concatenateCubes(cubelist, cataxis::CategoricalAxis)

Concatenates a vector of datacubes that have identical axes to a new single cube along the new
axis `cataxis`
"""
function concatenatecubes(cl,cataxis::CubeAxis)
  length(cataxis.values)==length(cl) || error("cataxis must have same length as cube list")
  axlist=axcopy.(caxes(cl[1]))
  T=eltype(cl[1])
  N=ndims(cl[1])
  cleaners = CleanMe[]
  append!(cleaners,cl[1].cleaner)
  for i=2:length(cl)
    all(caxes(cl[i]).==axlist) || error("All cubes must have the same axes, cube number $i does not match")
    eltype(cl[i])==T || error("All cubes must have the same element type, cube number $i does not match")
    ndims(cl[i])==N || error("All cubes must have the same dimension")
    append!(cleaners,cl[i].cleaner)
  end
  props=mapreduce(getattributes,merge,cl,init=getattributes(cl[1]))
  YAXArray([axlist...,cataxis],diskstack([c.data for c in cl]),props, cleaners)
end
function concatenatecubes(;kwargs...)
  cubenames = String[]
  for (n,c) in kwargs
    push!(cubenames,string(n))
  end
  cubes = map(i->i[2],collect(kwargs))
  findAxis("Variable",cubes[1]) === nothing || error("Input cubes must not contain a variable kwarg concatenation")
  concatenateCubes(cubes, CategoricalAxis("Variable",cubenames))
end


text:
function rr_bubble_pressure_refine(model,x,y,vl,vv,T)
    pl = pressure(model,vl,T,x)
    pv = pressure(model,vv,T,y)
    p = (pl+pv)/2
    vl = volume(model,p,T,x,phase=:l)
    vv = volume(model,p,T,y,phase=:v)
    μ_l = vt_chemical_potential(model,v_l,T,x)
    #@show "c"
    μ_v = vt_chemical_potential(model,v_v,T,y)
    K = μ_v/μ_l
    y = K.*x
    @show vl,vv,y
    return model,x,y,vl,vv,T
end


text:
@inline incrementp(A::AbstractStridedPointer{T,3} where T, a::Ptr) = VectorizationBase.increment_ptr(A, a, (Zero(), Zero(), One()))
@inline increment2(B::AbstractStridedPointer{T,2} where T, b::Ptr, ::StaticInt{nᵣ}) where {nᵣ} = VectorizationBase.increment_ptr(B, b, (Zero(), StaticInt{nᵣ}()))
@inline increment1(C::AbstractStridedPointer{T,2} where T, c::Ptr, ::StaticInt{mᵣW}) where {mᵣW} = VectorizationBase.increment_ptr(C, c, (StaticInt{mᵣW}(), Zero()))
macro kernel(pack::Bool, ex::Expr)
  ex.head === :for || throw(ArgumentError("Must be a matmul for loop."))
  mincrements = Expr[:(c = increment1(C, c, mᵣW)), :(ãₚ = incrementp(Ãₚ, ãₚ)), :(m = vsub_nsw(m, mᵣW))]
  # massumes = Expr[:(assume(m < mᵣW)),
  #                 :(assume(VectorizationBase.vgt(ãₚ, VectorizationBase.increment_ptr($(esc(:Ãₚ)), ãₚ, (vsub_nsw($(esc(:M)), mᵣW), LoopVectorization.Zero())), $(esc(:Ãₚ))))),
  #                 :(assume(VectorizationBase.vgt(c, VectorizationBase.increment_ptr($(esc(:C)), c, (vsub_nsw($(esc(:M)), mᵣW), LoopVectorization.Zero())), $(esc(:C)))))]
  offsetprecalc = GlobalRef(VectorizationBase,:offsetprecalc)
  preheader = quote
    mᵣ, nᵣ = matmul_params(Val($(esc(:T))))
    mᵣW = pick_vector_width($(esc(:T))) * mᵣ
    m = $(esc(:M)) % Int32
    n = $(esc(:N)) % Int32
    Ãₚ = $(esc(:Ãₚ))
    B = $offsetprecalc($(esc(:B)), Val{(9,9)}())
    C = $offsetprecalc($(esc(:C)), Val{(9,9)}())
    b = pointer(B); c = pointer(C); ãₚ = pointer(Ãₚ)
  end
  if pack
    push!(mincrements, :(a = increment1(A, a, mᵣW)))
    push!(preheader.args, :(A = $(esc(:A))), :(a = pointer(A)))
    areconstruct = Expr[:($(esc(:A)) = VectorizationBase.reconstruct_ptr(A, a))]
    # push!(massumes, :(assume(VectorizationBase.vgt(a, VectorizationBase.increment_ptr($(esc(:A)), a, (vsub_nsw($(esc(:M)), mᵣW), LoopVectorization.Zero())), $(esc(:A))))))
  else
    Ainit = areconstruct = Expr[]
  end
  lvkern = esc(:(@turbo inline=true $ex))

  loopnest = quote
    let ãₚ = ãₚ, c = c, $(esc(:B)) = VectorizationBase.reconstruct_ptr(B, b), m = m
      while m ≥ mᵣW#VectorizationBase.vle(a, amax, A)
        let $(esc(:M)) = mᵣW, $(esc(:N)) = nᵣ, $(esc(:Ãₚ)) = VectorizationBase.reconstruct_ptr(droplastdim(Ãₚ), ãₚ), $(esc(:C)) = VectorizationBase.reconstruct_ptr(C, c), $(areconstruct...)
          $lvkern
          $(mincrements...)
        end
      end
      if m > zero(Int32)#vne(a, amax, A)
        let $(esc(:M)) = UpperBoundedInteger((m%UInt)%Int, mᵣW - One()), $(esc(:N)) = nᵣ, $(esc(:Ãₚ)) = VectorizationBase.reconstruct_ptr(droplastdim(Ãₚ), ãₚ), $(esc(:C)) = VectorizationBase.reconstruct_ptr(C, c), $(areconstruct...)
          # $(massumes...)
          $lvkern
        end
      end
    end
  end

  if !pack
    loopnest = quote
      while n ≥ nᵣ#VectorizationBase.vle(c, cmax, C)
        $loopnest
        c = increment2(C, c, nᵣ)
        b = increment2(B, b, nᵣ)
        n = vsub_nsw(n, nᵣ)
      end
      if n > zero(Int32)#vne(c, cmax, C)
        let $(esc(:B)) = VectorizationBase.reconstruct_ptr(B, b), m = m
          while m ≥ mᵣW#VectorizationBase.vle(a, amax, A)
            let $(esc(:M)) = mᵣW, $(esc(:N)) = UpperBoundedInteger((n%UInt)%Int, nᵣ - One()), $(esc(:Ãₚ)) = VectorizationBase.reconstruct_ptr(droplastdim(Ãₚ), ãₚ), $(esc(:C)) = VectorizationBase.reconstruct_ptr(C, c), $(areconstruct...)
              # assume(n < nᵣ)
              # assume((VectorizationBase.vgt(c, VectorizationBase.increment_ptr($(esc(:C)), c, (Zero(), vsub_nsw($(esc(:N)), nᵣ))), $(esc(:C)))))
              # assume((VectorizationBase.vgt(b, VectorizationBase.increment_ptr($(esc(:B)), b, (Zero(), vsub_nsw($(esc(:N)), nᵣ))), $(esc(:B)))))
              $lvkern
              $(mincrements...)
            end
          end
          if m > zero(Int32)#vne(a, amax, A)
            let $(esc(:M)) = UpperBoundedInteger((m%UInt)%Int, mᵣW - One()), $(esc(:N)) = UpperBoundedInteger((n%UInt)%Int, nᵣ - One()), $(esc(:Ãₚ)) = VectorizationBase.reconstruct_ptr(droplastdim(Ãₚ), ãₚ), $(esc(:C)) = VectorizationBase.reconstruct_ptr(C, c), $(areconstruct...)
              # $(massumes...)
              # assume(n < nᵣ)
              # assume((VectorizationBase.vgt(c, VectorizationBase.increment_ptr($(esc(:C)), c, (Zero(), vsub_nsw($(esc(:N)), nᵣ))), $(esc(:C)))))
              # assume((VectorizationBase.vgt(b, VectorizationBase.increment_ptr($(esc(:B)), b, (Zero(), vsub_nsw($(esc(:N)), nᵣ))), $(esc(:B)))))
              $lvkern
            end
          end
        end
      end
    end
  end
  Expr(:block, preheader, loopnest)
end
@inline function loopmul!(C, A, B, α, β, M, K, N)
    @turbo for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
        Cₘₙ = zero(eltype(C))
        for k ∈ CloseOpen(K)
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = α * Cₘₙ + β * C[m,n]
    end
    nothing
end
@inline function ploopmul!(C::AbstractStridedPointer{T}, Ãₚ, B, α, β, M, K, N) where {T}
  @kernel false for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
    Cₘₙ = zero(eltype(C))
    for k ∈ CloseOpen(K)
      Cₘₙ += Ãₚ[m,k] * B[k,n]
    end
    C[m,n] = α * Cₘₙ + β * C[m,n]
  end
  nothing
end
@inline function packamul!(
    C::AbstractStridedPointer{T}, Ãₚ, A, B,
    α, β, M, K, N
  ) where {T}
  @kernel true for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
    Cₘₙ = zero(eltype(C))
    for k ∈ CloseOpen(K)
      Aₘₖ = A[m,k]
      Cₘₙ += Aₘₖ * B[k,n]
      Ãₚ[m,k] = Aₘₖ
    end
    C[m,n] = α * Cₘₙ + β * C[m,n]
  end
end
@inline function alloc_a_pack(K, ::Val{T}) where {T}
  buffer = first_cache_buffer(Val(T))
  alloc_a_pack(K, buffer), buffer
end
@inline alloc_a_pack(K, buffer::MemoryBuffer) = alloc_a_pack(K, align(pointer(buffer)))
@inline function alloc_a_pack(K, bufferptr::Ptr{T}) where {T}
  mᵣ, nᵣ = matmul_params(Val(T))
  mᵣW = mᵣ * pick_vector_width(T)
  Apack = default_zerobased_stridedpointer(bufferptr, (One(), mᵣW, mᵣW * K)) # mᵣW x K x cld(M, mᵣW)
  Apack
end
function packaloopmul!(
    C::AbstractStridedPointer{T},
    A::AbstractStridedPointer,
    B::AbstractStridedPointer,
    α, β, M, K, N
) where {T}
  Ãₚ, buffer = alloc_a_pack(K, Val(T))
  GC.@preserve buffer begin
    Mᵣ, Nᵣ = matmul_params(Val(T))
    packamul!(C, Ãₚ, A, B, α, β, M, K, Nᵣ)
    ploopmul!(gesp(C, (Zero(), Nᵣ)), Ãₚ, gesp(B, (Zero(), Nᵣ)), α, β, M, K, N - Nᵣ)
  end
  nothing
end

@inline function inlineloopmul!(C, A, B, α, β, M, K, N)
    @turbo inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
        Cₘₙ = zero(eltype(C))
        for k ∈ CloseOpen(K)
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n]  = α * Cₘₙ + β * C[m,n]
    end
    C
end


text:
# A class of Matern isotropic functions including the Matrern 3/2 and 5/2, where d= 3 or 5. Also the exponential function, where d=1
abstract type MaternIso <: Isotropic{Euclidean} end
abstract type MaternARD <: StationaryARD{WeightedEuclidean} end

@inline function dKij_dθp(mat::MaternARD, X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, p::Int, dim::Int)
    r=distij(metric(mat),X1, X2,i,j,dim)
    if p <= dim
        wdiffp=dist2ijk(metric(mat),X1, X2,i,j,p)
        if wdiffp > 0
            return dk_dll(mat,r,wdiffp)
        else
            return zero(dk_dll(mat,r,wdiffp))
        end
    elseif p==dim+1
        return dk_dlσ(mat, r)
    else
        return NaN
    end
end
@inline function dKij_dθp(mat::MaternARD, X1::AbstractMatrix, X2::AbstractMatrix, data::StationaryARDData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(mat,X1,X2,i,j,p,dim)
end

@inline function dk_dθp(mat::MaternIso, r::Real, p::Int)
    if p==1
        return r == 0.0 ? 0.0 : dk_dll(mat, r)
    elseif p==2
        return dk_dlσ(mat, r)
    else
        return NaN
    end
end

include("mat12_iso.jl")
include("mat12_ard.jl")

include("mat32_iso.jl")
include("mat32_ard.jl")

include("mat52_iso.jl")
include("mat52_ard.jl")

"""
    Matern(ν::Real, ll::Union{Real,Vector{Real}}, lσ::Real)

Create Matérn kernel of type `ν` (i.e. `ν = 1/2`, `ν = 3/2`, or `ν = 5/2`) with length scale
`exp.(ll)` and signal standard deviation `exp(σ)`.

See also [`Mat12Iso`](@ref), [`Mat12Ard`](@ref), [`Mat32Iso`](@ref), [`Mat32Ard`](@ref),
[`Mat52Iso`](@ref), and [`Mat52Ard`](@ref).
"""
function Matern(ν::Real, ll::Real, lσ::Real)
    if ν==1/2
        kern = Mat12Iso(ll, lσ)
    elseif ν==3/2
        kern = Mat32Iso(ll, lσ)
    elseif ν==5/2
        kern = Mat52Iso(ll, lσ)
    else throw(ArgumentError("Only Matern 1/2, 3/2 and 5/2 are implementable"))
    end
    return kern
end

function Matern(ν::Real, ll::Vector{<:Real}, lσ::Real)
    if ν==1/2
        kern = Mat12Ard(ll, lσ)
    elseif ν==3/2
        kern = Mat32Ard(ll, lσ)
    elseif ν==5/2
        kern = Mat52Ard(ll, lσ)
    else throw(ArgumentError("Only Matern 1/2, 3/2 and 5/2 are implementable"))
    end
    return kern
end

@deprecate Mat Matern


text:
struct LinearizedStep{T}
	z::Vector{T}
	θ::Vector{T}
	κ::Vector{T}
	r::Vector{T}
	rz::Matrix{T}
	rθ::Matrix{T}
end

function LinearizedStep(s::Simulation, z::Vector{T}, θ::Vector{T}, κ::T) where T
	model = s.model
	env = s.env

	nz = num_var(model, env)
	nθ = num_data(model)

	z0 = copy(z)
	θ0 = copy(θ)
	κ0 = [κ]
	r0 = zeros(nz)
	rz0 = zeros(nz, nz)
	rθ0 = zeros(nz, nθ)

	s.res.r!(r0, z0, θ0, κ0)
	s.res.rz!(rz0, z0, θ0)
	s.res.rθ!(rθ0, z0, θ0)

	return LinearizedStep(z0, θ0, κ0, r0, rz0, rθ0)
end

"""
	Create LinearizedStep.
"""
function LinearizedStep(s::Simulation)
	model = s.model
	env = s.env

	nz = num_var(model, env)
	nθ = num_data(model)
	z0 = zeros(nz)
	θ0 = zeros(nθ)

	κ0 = 0.0

	return LinearizedStep(s, z0, θ0, κ0)
end

function update!(lin::LinearizedStep{T}, s::Simulation{T}, z::Vector{T}, θ::Vector{T}) where T
	lin.z .= z
	lin.θ .= θ
	s.res.r!(lin.r, z, θ, lin.κ)
	s.res.rz!(lin.rz, z, θ)
	s.res.rθ!(lin.rθ, z, θ)
	return nothing
end