text:
@testset "ModelKit - Instructions" begin
    @testset "evaluate" begin
        function test_instr_interpreter(instr, vars, expected)
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(instr, vars...))
            I = ModelKit.Interpreter(list, [out]; variables = vars)
            cache =
                ModelKit.InterpreterCache(zeros(Expression, ModelKit.cache_min_length(I)))
            u = Expression[0]
            ModelKit.execute!(u, I, Variable.(vars), nothing, cache)
            @test expand(u[1] - expected) == 0
        end
        function test_instr_pow_interpreter(instr, r, expected)
            @var x
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(instr, :x, r))
            I = ModelKit.Interpreter(list, [out]; variables = [:x])
            cache =
                ModelKit.InterpreterCache(zeros(Expression, ModelKit.cache_min_length(I)))
            u = Expression[0]
            ModelKit.execute!(u, I, [x], nothing, cache)
            @test expand(u[1] - expected) == 0
        end

        @var x y z
        test_instr_interpreter(ModelKit.INSTR_NEG, [:x], -x)
        test_instr_interpreter(ModelKit.INSTR_SQR, [:x], x^2)
        test_instr_interpreter(ModelKit.INSTR_SIN, [:x], sin(x))
        test_instr_interpreter(ModelKit.INSTR_COS, [:x], cos(x))

        test_instr_pow_interpreter(ModelKit.INSTR_POW, 2, x^2)
        test_instr_pow_interpreter(ModelKit.INSTR_POW, 5, x^5)
        test_instr_pow_interpreter(ModelKit.INSTR_POW, 6, x^6)
        test_instr_interpreter(ModelKit.INSTR_ADD, [:x, :y], x + y)
        test_instr_interpreter(ModelKit.INSTR_MUL, [:x, :y], x * y)
        test_instr_interpreter(ModelKit.INSTR_SUB, [:x, :y], x - y)
        test_instr_interpreter(ModelKit.INSTR_DIV, [:x, :y], x / y)

        test_instr_interpreter(ModelKit.INSTR_MULADD, [:x, :y, :z], x * y + z)
        test_instr_interpreter(ModelKit.INSTR_MULSUB, [:x, :y, :z], x * y - z)
        test_instr_interpreter(ModelKit.INSTR_SUBMUL, [:x, :y, :z], z - x * y)
    end

    @testset "evaluate with arb" begin
        function test_instr_interpreter(instr, vars, values, expected)
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(instr, vars...))
            I = ModelKit.Interpreter(list, [out]; variables = vars)
            cache = ModelKit.InterpreterCache(Arblib.Acb, I; prec = 256)
            u = Arblib.AcbVector(1; prec = 256)
            ModelKit.execute!(u, I, values, nothing, cache)
            @test Float64(Arblib.get!(Arblib.Mag(), u[1] - expected)) < 1e-14
        end
        function test_instr_pow_interpreter(instr, r, x, expected)
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(instr, :x, r))
            I = ModelKit.Interpreter(list, [out]; variables = [:x])
            cache = ModelKit.InterpreterCache(Arblib.Acb, I; prec = 256)
            u = Arblib.AcbVector(1; prec = 256)
            ModelKit.execute!(u, I, [x], nothing, cache)
            @test Float64(Arblib.get!(Arblib.Mag(), u[1] - expected)) < 1e-14
        end

        x, y, z = Arblib.Acb.(randn(ComplexF64, 3); prec = 256)
        test_instr_interpreter(ModelKit.INSTR_NEG, [:x], [x], -x)
        test_instr_interpreter(ModelKit.INSTR_SQR, [:x], [x], x^2)
        test_instr_interpreter(ModelKit.INSTR_SIN, [:x], [x], sin(x))
        test_instr_interpreter(ModelKit.INSTR_COS, [:x], [x], cos(x))

        test_instr_pow_interpreter(ModelKit.INSTR_POW, 2, x, x^2)
        test_instr_pow_interpreter(ModelKit.INSTR_POW, 5, x, x^5)
        test_instr_pow_interpreter(ModelKit.INSTR_POW, 6, x, x^6)
        test_instr_interpreter(ModelKit.INSTR_ADD, [:x, :y], [x, y], x + y)
        test_instr_interpreter(ModelKit.INSTR_MUL, [:x, :y], [x, y], x * y)
        test_instr_interpreter(ModelKit.INSTR_SUB, [:x, :y], [x, y], x - y)
        test_instr_interpreter(ModelKit.INSTR_DIV, [:x, :y], [x, y], x / y)

        test_instr_interpreter(ModelKit.INSTR_MULADD, [:x, :y, :z], [x, y, z], x * y + z)
        test_instr_interpreter(ModelKit.INSTR_MULSUB, [:x, :y, :z], [x, y, z], x * y - z)
        test_instr_interpreter(ModelKit.INSTR_SUBMUL, [:x, :y, :z], [x, y, z], z - x * y)
    end

    @testset "gradient" begin
        function test_instr_interpreter_grad(instr, vars, expected)
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(instr, vars...))
            grad_list, ev_out, jac_out = ModelKit.gradient(list, [out], vars)
            I = ModelKit.Interpreter(
                grad_list,
                [ev_out; vec(jac_out)];
                variables = vars,
                nexpressions = 1,
            )
            cache =
                ModelKit.InterpreterCache(zeros(Expression, ModelKit.cache_min_length(I)))
            u = zeros(Expression, 1 + length(vars))
            ModelKit.execute!(u, I, Variable.(vars), nothing, cache)
            @test expand.(u - expected) == zeros(Expression, 1 + length(vars))
        end
        function test_instr_pow_interpreter_grad(instr, r, expected)
            @var x
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(instr, :x, r))
            grad_list, ev_out, jac_out = ModelKit.gradient(list, [out], [:x])
            I = ModelKit.Interpreter(
                grad_list,
                [ev_out; vec(jac_out)];
                variables = [:x],
                nexpressions = 1,
            )
            cache =
                ModelKit.InterpreterCache(zeros(Expression, ModelKit.cache_min_length(I)))
            u = Expression[0, 0]
            ModelKit.execute!(u, I, [x], nothing, cache)
            @test expand.(u - expected) == Expression[0, 0]
        end
        @var x y z

        test_instr_interpreter_grad(ModelKit.INSTR_NEG, [:x], [-x, -1])
        test_instr_interpreter_grad(ModelKit.INSTR_SQR, [:x], [x^2, 2x])
        test_instr_interpreter_grad(ModelKit.INSTR_SIN, [:x], [sin(x), cos(x)])
        test_instr_interpreter_grad(ModelKit.INSTR_COS, [:x], [cos(x), -sin(x)])

        test_instr_pow_interpreter_grad(ModelKit.INSTR_POW, 2, [x^2, 2x])
        test_instr_pow_interpreter_grad(ModelKit.INSTR_POW, 5, [x^5, 5x^4])
        test_instr_pow_interpreter_grad(ModelKit.INSTR_POW, 6, [x^6, 6x^5])
        test_instr_interpreter_grad(ModelKit.INSTR_ADD, [:x, :y], [x + y, 1, 1])
        test_instr_interpreter_grad(ModelKit.INSTR_MUL, [:x, :y], [x * y, y, x])
        test_instr_interpreter_grad(ModelKit.INSTR_SUB, [:x, :y], [x - y, 1, -1])
        test_instr_interpreter_grad(ModelKit.INSTR_DIV, [:x, :y], [x / y, 1 / y, -x / y^2])

        test_instr_interpreter_grad(
            ModelKit.INSTR_MULADD,
            [:x, :y, :z],
            [x * y + z, y, x, 1],
        )
        test_instr_interpreter_grad(
            ModelKit.INSTR_MULSUB,
            [:x, :y, :z],
            [x * y - z, y, x, -1],
        )
        test_instr_interpreter_grad(
            ModelKit.INSTR_SUBMUL,
            [:x, :y, :z],
            [z - x * y, -y, -x, 1],
        )
    end

    @testset "taylor" begin
        function test_instr_taylor1(instr, op, K, M)
            vars = [:x]
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(instr, vars...))
            I = ModelKit.Interpreter(list, [out]; variables = vars)
            tape = ModelKit.TaylorVector{K + 1}(Expression, ModelKit.cache_min_length(I))
            cache = ModelKit.InterpreterCache(tape)
            u = ModelKit.TaylorVector{K + 1}(Expression, 1)
            @var x[0:M]
            X = ModelKit.TaylorVector{M + 1}(Expression, 1)
            X[1] = tuple(x...)
            ModelKit.execute!(u, Val(K), I, X, nothing, cache)

            @var λ
            tx = sum(λ .^ (k - 1) .* xk for (k, xk) in enumerate(x))
            true_value = [differentiate(op(tx), λ, k)(λ => 0) / factorial(k) for k = 0:K]
            @test expand.(u[1] .- true_value) == zeros(Expression, K + 1)
        end

        function test_instr_taylor_pow(r, K, M)
            vars = [:x]
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(ModelKit.INSTR_POW, :x, r))
            I = ModelKit.Interpreter(list, [out]; variables = vars)
            tape = ModelKit.TaylorVector{K + 1}(Expression, ModelKit.cache_min_length(I))
            cache = ModelKit.InterpreterCache(tape)
            u = ModelKit.TaylorVector{K + 1}(Expression, 1)
            @var x[0:M]
            X = ModelKit.TaylorVector{M + 1}(Expression, 1)
            X[1] = tuple(x...)
            ModelKit.execute!(u, Val(K), I, X, nothing, cache)

            @var λ
            tx = sum(λ .^ (k - 1) .* xk for (k, xk) in enumerate(x))
            true_value = [differentiate(tx^r, λ, k)(λ => 0) / factorial(k) for k = 0:K]
            @test expand.(u[1] .- true_value) == zeros(Expression, K + 1)
        end
        function test_instr_taylor2(instr, op, K, M)
            vars = [:x, :y]
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(instr, vars...))
            I = ModelKit.Interpreter(list, [out]; variables = vars)
            tape = ModelKit.TaylorVector{K + 1}(Expression, ModelKit.cache_min_length(I))
            cache = ModelKit.InterpreterCache(tape)
            u = ModelKit.TaylorVector{K + 1}(Expression, 1)
            @var x[0:M] y[0:M]
            X = ModelKit.TaylorVector{M + 1}(Expression, 2)
            X[1] = tuple(x...)
            X[2] = tuple(y...)
            ModelKit.execute!(u, Val(K), I, X, nothing, cache)

            @var λ
            tx = sum(λ .^ (k - 1) .* xk for (k, xk) in enumerate(x))
            ty = sum(λ .^ (k - 1) .* yk for (k, yk) in enumerate(y))
            true_value =
                [differentiate(op(tx, ty), λ, k)(λ => 0) / factorial(k) for k = 0:K]
            @test expand.(u[1] .- true_value) == zeros(Expression, K + 1)
        end

        function test_instr_taylor3(instr, op, K, M)
            vars = [:x, :y, :z]
            list = ModelKit.InstructionList()
            out = push!(list, ModelKit.Instruction(instr, vars...))
            I = ModelKit.Interpreter(list, [out]; variables = vars)
            tape = ModelKit.TaylorVector{K + 1}(Expression, ModelKit.cache_min_length(I))
            cache = ModelKit.InterpreterCache(tape)
            u = ModelKit.TaylorVector{K + 1}(Expression, 1)
            @var x[0:M] y[0:M] z[0:M]
            X = ModelKit.TaylorVector{M + 1}(Expression, 3)
            X[1] = tuple(x...)
            X[2] = tuple(y...)
            X[3] = tuple(z...)
            ModelKit.execute!(u, Val(K), I, X, nothing, cache)

            @var λ
            tx = sum(λ .^ (k - 1) .* xk for (k, xk) in enumerate(x))
            ty = sum(λ .^ (k - 1) .* yk for (k, yk) in enumerate(y))
            tz = sum(λ .^ (k - 1) .* yk for (k, yk) in enumerate(z))
            true_value =
                [differentiate(op(tx, ty, tz), λ, k)(λ => 0) / factorial(k) for k = 0:K]
            @test expand.(u[1] .- true_value) == zeros(Expression, K + 1)
        end

        for K = 0:4, M = 0:K
            test_instr_taylor1(ModelKit.INSTR_NEG, -, K, M)
            test_instr_taylor1(ModelKit.INSTR_SQR, a -> a^2, K, M)
            test_instr_taylor1(ModelKit.INSTR_SIN, sin, K, M)
            test_instr_taylor1(ModelKit.INSTR_COS, cos, K, M)
            test_instr_taylor2(ModelKit.INSTR_ADD, +, K, M)
            test_instr_taylor2(ModelKit.INSTR_SUB, -, K, M)
            test_instr_taylor2(ModelKit.INSTR_MUL, *, K, M)
            test_instr_taylor2(ModelKit.INSTR_DIV, /, K, M)
            test_instr_taylor_pow(2, K, M)
            test_instr_taylor_pow(4, K, M)
            test_instr_taylor_pow(7, K, M)
            test_instr_taylor3(ModelKit.INSTR_MULADD, muladd, K, M)
            test_instr_taylor3(ModelKit.INSTR_MULSUB, (x, y, z) -> x * y - z, K, M)
            test_instr_taylor3(ModelKit.INSTR_SUBMUL, (x, y, z) -> z - x * y, K, M)
        end
    end
end


text:
"""
    IRLS(d,operators,parameters;<keyword arguments>)

Non-quadratic regularization with Iteratively Reweighted Least Squares (IRLS).

# Arguments
- `Niter_external=3`
- 'Niter_internal=10'
- `mu=0`

"""
function IRLS(d,operators,parameters;Niter_external=3,Niter_internal=10,mu=0)

	cost = Float64[]
	weights = ones(Float64,size(d))
	parameters[end][:w] = weights
	m = []
	for iter_external = 1 : Niter_external
		v,cost1 = ConjugateGradients(d,operators,parameters,Niter=Niter_internal,mu=mu)
		append!(cost,cost1)
		m = v .* weights
		weights = abs.(m./maximum(abs.(m[:])))
		parameters[end][:w] = weights
	end

	return m, cost

end

ext: jl


text:
module VTK

export writevtk, writepvtu

include("writemesh.jl")
include("writevtk.jl")
include("writepvtu.jl")

function __init__()
    tictoc()
end

end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

################################################################################
# Common models used in testing
################################################################################
function rand_dims(max_ndim=6)
  tuple(rand(1:10, rand(1:max_ndim))...)
end

function mlp2()
  data = mx.Variable(:data)
  out = mx.FullyConnected(data, name=:fc1, num_hidden=1000)
  out = mx.Activation(out, act_type=:relu)
  out = mx.FullyConnected(out, name=:fc2, num_hidden=10)
  return out
end

function mlpchain()
  mx.@chain mx.Variable(:data) =>
            mx.FullyConnected(name=:fc1, num_hidden=1000) =>
            mx.Activation(act_type=:relu) =>
            mx.FullyConnected(name=:fc2, num_hidden=10)
end

"""
execution helper of SymbolicNode
"""
function exec(x::mx.SymbolicNode; feed...)
  ks, vs = zip(feed...)
  vs′ = mx.NDArray.(vs)

  e = mx.bind(x, context = mx.cpu(), args = Dict(zip(ks, vs′)))
  mx.forward(e)
  e.outputs
end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
# SPDX-License-Identifier: MIT

const N_THREADS_1D = 256
const N_THREADS_2D = (16, 16)

#======================================================================
elmul_svdm!

  Element-wise multiplication of a sparse vector and with each row of a
  dense matrix. The results is a dense matrix.
======================================================================#

function elmul_svdm!(out::AbstractMatrix{T},
                     x::CuSparseVector{T},
                     M::AbstractMatrix{T}) where T
    @boundscheck (size(out) == size(M) || throw(DimensionMistmatch()))
    @boundscheck (length(x) == size(M, 1) || throw(DimensionMistatch()))

    fill!(out, zero(T))
    I = SparseArrays.nonzeroinds(x)
    V = nonzeros(x)
    N = length(x.nzVal)

    if N > 0
        nb = ceil(Int, N/N_THREADS_1D)
        @cuda threads=N_THREADS_1D blocks=nb _cukernel_elmul_svdm!(out, V, I, M)
    end
    return out
end

function _cukernel_elmul_svdm!(out, nzVal, nzInd, dsMat)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i = index:stride:length(nzInd)
        idx = nzInd[i]
        for j in 1:size(dsMat, 2)
            out[idx, j] = nzVal[i] * dsMat[idx,j ]
        end
    end
    return
end

#======================================================================
elmul_svdv!

  Element-wise multiplication of a sparse vector and a dense vector.
======================================================================#

function elmul_svdv!(out::CuVector,
                     x::CuSparseVector{T},
                     y::CuVector{T}) where T
    @boundscheck size(out) == size(y) || throw(DimensionMistmatch())
    @boundscheck length(x) == length(y) || throw(DimensionMistatch())

    fill!(out, zero(T))
    I = SparseArrays.nonzeroinds(x)
    V = nonzeros(x)
    N = length(x.nzVal)

    if N > 0
        nb = ceil(Int, N/N_THREADS_1D)
        @cuda threads=N_THREADS_1D blocks=nb _cukernel_elmul_svdv!(out, V, I, y)
    end
    return out
end

function _cukernel_elmul_svdv!(out, nzVal, nzInd, dsVec)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i = index:stride:length(nzInd)
        idx = nzInd[i]
        out[idx] = nzVal[i] * dsVec[idx]
    end
    return
end

#======================================================================
mul_smdv!

  Product of a sparse matrix and a dense vector.
======================================================================#

function mul_smdv!(out::CuVector{T},
                   M::CuSparseMatrixCSR{T},
                   x::CuVector{T}) where T
    @boundscheck size(out) == size(x) || throw(DimensionMistmatch())
    @boundscheck size(M, 2) == length(x) || throw(DimensionMistatch())

    fill!(out, zero(T))
    rowPtr = M.rowPtr
    colVal = M.colVal
    nzVal = M.nzVal
    N = size(M,1)

    if length(nzVal) > 0
        nb = ceil(Int, N/N_THREADS_1D)
        @cuda threads=N_THREADS_1D blocks=nb _cukernel_mul_smdv!(out, N, rowPtr,
                                                                 colVal, nzVal, x)
    end
    return out
end

function _cukernel_mul_smdv!(out, N, rowPtr, colVal, nzVal, x)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x

    for ti = index:stride:N
        for i = rowPtr[ti]:(rowPtr[ti+1]-1)
            out[ti] += nzVal[i] * x[colVal[i]]
        end
    end
    return
end

#======================================================================
mul_smdm!

  Product of a sparse matrix and a dense matrix.
======================================================================#

function mul_smdm!(out::AbstractMatrix{T},
                   X::CuSparseMatrixCSR{T},
                   Y::AbstractMatrix{T}) where T
    @boundscheck size(out) == (size(X,1), size(Y,2)) || throw(DimensionMistmatch())
    @boundscheck size(X, 2) == size(Y,1) || throw(DimensionMistatch())

    fill!(out, zero(T))
    rowPtr = X.rowPtr
    colVal = X.colVal
    nzVal = X.nzVal
    N = size(X, 1)

    if length(nzVal) > 0
        nb = ceil(Int, N/N_THREADS_1D)
        @cuda threads=N_THREADS_1D blocks=nb _cukernel_mul_smdm!(out, N, rowPtr,
                                                                 colVal, nzVal, Y)
    end
    return out
end

function _cukernel_mul_smdm!(out, N, rowPtr, colVal, nzVal, Y)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x

    for ti = index:stride:N
        for i = rowPtr[ti]:(rowPtr[ti+1]-1)
            for j in 1:size(Y, 2)
                out[ti,j] += nzVal[i] * Y[colVal[i],j]
            end
        end
    end
    return
end

#======================================================================
mul_dmdm!

  Product of a dense matrix and a dense matrix.
======================================================================#

function _cukernel_mul_dmdm!(out, X, Y)
    index_i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride_i = blockDim().x * gridDim().x
    index_j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_j = blockDim().y * gridDim().y

    for i = index_i:stride_i:size(X, 1)
        for j in index_j:stride_j:size(Y, 2)
            for k in 1:size(X, 2)
                out[i,j] += X[i,k] * Y[k,j]
            end
        end
    end
    return
end

function mul_dmdm!(out::AbstractMatrix{T}, X::CuMatrix{T},
                   Y::CuMatrix{T}) where T
    fill!(out, zero(T))
    M, N = size(X, 1), size(Y, 2)

    nb = (
        ceil(Int, M/N_THREADS_2D[1]),
        ceil(Int, N/N_THREADS_2D[2])
    )
    @cuda threads=N_THREADS_2D blocks=nb _cukernel_mul_dmdm!(out, X, Y)
    return out
end

max_stars_repo_name: MartinKocour/MarkovModels

max_issues_repo_name: MartinKocour/MarkovModels

max_forks_repo_name: MartinKocour/MarkovModels


text:
# rotate2d.jl

using MIRT: rotate2d
using Test: @test

@test [rotate2d(1, 1, π/2)...] ≈ [(1.,-1.)...]


text:
struct TestDatum
    input::Vector{Float64};
    output::Vector{Float64};
end

input(t::TestDatum) = t.input;
output(t::TestDatum) = t.output;
inlength(t::TestDatum) = length(t.input);
outlength(t::TestDatum) = length(t.output);


text:
"""
Remove a template signal from each column of an array

#### Arguments

* `signals`: Original signals to be modified  (samples x channels)
* `template`: Template to remove from each signal

#### Returns
Signals with template removed
"""
function remove_template(signals::Array{T, 2}, template::AbstractVector{T}) where T <: AbstractFloat

    @assert size(signals, 1) == size(template, 1)

    for chan = 1:size(signals)[end]
        signals[:, chan] -=  template
    end

    return signals
end

"""
Re reference a signals to specific signal channel by index.

If multiple channels are specififed, their average is used as the reference.

#### Arguments

* `signals`: Original signals to be modified
* `refChan`: Index of channels to be used as reference

#### Returns

Rereferenced signals
"""
function rereference(signals::Array{T, 2}, refChan::Union{Int, Array{Int}}) where T <: AbstractFloat

    @debug("Re referencing $(size(signals)[end]) channels to $(length(refChan)) channels")
    @debug("Reference channels = $refChan")

    reference_signal = signals[:, refChan]

    # If using the average of several channels
    if size(reference_signal, 2) > 1
        reference_signal = vec(mean(reference_signal, dims = 2))
    end

    remove_template(signals, reference_signal)
end

"""
Re-reference a signals to specific signal channel by name.

If multiple channels are specififed, their average is used as the reference.
Or you can specify to use the `average` reference.

#### Arguments

* `signals`: Original signals to be modified
* `refChan`: List of channels to be used as reference or `average`
* `chanNames`: List of channel names associated with signals array

#### Returns

Rereferenced signals
"""
function rereference(signals::Array{T, 2}, refChan::S, chanNames::Vector{S}) where {S <: AbstractString, T <: AbstractFloat}

    @debug("Reference channels = $refChan")

    if refChan == "car" || refChan == "average"
        refChan_Idx = collect(1:size(signals, 2))
    elseif isa(refChan, AbstractString)
        refChan_Idx = something(findfirst(isequal(refChan), chanNames), 0)
    end

    if refChan == 0; error("Requested channel is not in the provided list of channels"); end

    rereference(signals, refChan_Idx)
end

function rereference(signals::Array{T, 2}, refChan::Vector{S}, chanNames::Vector{S}) where {S <: AbstractString, T <: AbstractFloat}

    @debug("Reference channels = $refChan")

    refChan_Idx = [something(findfirst(isequal(i), chanNames), 0)  for i = refChan]

    if refChan == 0; error("Requested channel is not in the provided list of channels"); end

    rereference(signals, refChan_Idx)
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using SparseArrays;
using ViewerGL, LinearAlgebra
GL = ViewerGL

store=[]
V, (VV, EV, FV, CV) = Lar.cuboid([1,1,1], true);
mybox = (V,CV,FV,EV);
str = Lar.Struct([ mybox ])
obj = Lar.struct2lar(str)
push!(store, obj)
V, (VV, EV, FV, CV) = Lar.cuboid([0.5,0.5,0.5], true);
mybox = (V,CV,FV,EV);
transl = Lar.t(0.25,0.25,0.25)
str = Lar.Struct([ transl, mybox ])
obj = Lar.struct2lar(str)
push!(store, obj)
str = Lar.Struct(store);
V,CV,FV,EV = Lar.struct2lar(str);
GL.VIEW([ GL.GLPol(V,CV, GL.COLORS[1]) ]);
cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W=convert(Lar.Points,V')
V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(W::Lar.Points, cop_EV::Lar.ChainOp, cop_FE::Lar.ChainOp);
triangulated_faces = Lar.triangulate2D(V, [copEV, copFE]);
FVs = convert(Array{Lar.Cells}, triangulated_faces);
VVV = convert(Lar.Points, V');
try
    GL.VIEW( GL.GLExplode(VVV,FVs,1,1,1,99,1) );
catch
    e="visualizzazione errata causa triangolazione errata"
    @show FVs
    println("La presenza di array vuoti indica la mal riuscita della triangolazione")
finally
    temp=[]
    for x in FVs
       if !isempty(x)
           push!(temp, x)
       end
    end
    GL.VIEW( GL.GLExplode(VVV,temp,1,1,1,99,1) );
end


text:
include("gdtw_graph_implementation.jl")

w = GDTWWorkspace(20,20)
@test w isa GDTWWorkspace

@testset "GDTW: test against graph implementation" begin
    ts = range(0, stop=4π, length=50)
    rnd = () -> rand(2)
    x = LinearInterpolation(sin.(ts)' .* rnd())
    y = LinearInterpolation(sin.(1.1 .* ts)' .* rnd())

    N = 30
    M = 20

    # choices that change the result
    for metric in (Distances.SqEuclidean(), (x, y) -> norm(x - y, 1)),
        symmetric in (true, false)

        kwargs = (N=N, M=M, verbose=false, metric=metric, symmetric = symmetric)

        data = prepare_gdtw(x, y; kwargs...)
        cost = DynamicAxisWarping.single_gdtw!(data)
        warp = copy(data.warp)

        # Check against graph based implementation
        # Should not change the result.
        for algo in (LightGraphs.dijkstra_shortest_paths,
 LightGraphs.desopo_pape_shortest_paths, LightGraphs.bellman_ford_shortest_paths)
            new_cost = DynamicAxisWarping.single_gdtw!(data, algo)
            @test new_cost ≈ cost
            @test data.warp ≈ warp
        end

    end

end

@testset "GDTW: test that refinements reduce cost" begin
    x = LinearInterpolation(rand(2, 20))
    y = LinearInterpolation(rand(2, 20))
    for symmetric in (true, false)
        costs = [ gdtw(x,y; max_iters = i, symmetric=symmetric)[1] for i = 1:3]
        @test issorted(costs; rev=true)
    end
end

@testset "GDTW: test that warping helps" begin
    x = LinearInterpolation(rand(20))
    y = LinearInterpolation(rand(20))
    for symmetric in (true, false)
        cost, ϕ, ψ = gdtw(x, y; symmetric=symmetric)
        t = range(0, stop = 1, length=20)
        @test norm(x.(t) - y.(t)) >= norm( x.(ϕ.(t)) - y.(ψ.(t)) )
    end
end

@testset "GDTW: Test objective function" begin
    # short time span to make it easy
    ts = range(0, stop=π/12, length=128)
    x = LinearInterpolation(sin.(ts))
    y = LinearInterpolation(sin.(1.1 .* ts))
    metric = (x,y) -> norm(x-y)
    λinst = .1
    λcum = .1
    for symmetric in (false, true)
        # We need high `N` and `M` to have low discretization error
        cost, ϕ, ψ = gdtw(x, y; symmetric=symmetric, N = 400, M = 300, λinst=λinst, λcum=λcum)
        ϕ′ = x -> ForwardDiff.derivative(ϕ, x)
        res = quadgk(t -> metric(x(ϕ(t)),  y(ψ(t))) + λinst*(ϕ′(t)-1)^2 + λcum*(ϕ(t) - t)^2, 0, 1 )[1]
        @test res ≈ cost rtol=1e-2
    end
end

@testset "GDTW: test symmetry" begin
    x = LinearInterpolation(rand(20))
    y = LinearInterpolation(rand(20))
    cost_xy, ϕ_xy, ψ_xy = gdtw(x, y)
    cost_yx, ϕ_yx, ψ_yx = gdtw(y, x)
    @test cost_xy ≈ cost_yx
    t = range(0, stop = 1, length=100)
    @test ϕ_xy.(t) ≈ ψ_yx.(t)
    @test ψ_xy.(t) ≈ ϕ_yx.(t)
end

@testset "GDTW: `iterative_gdtw!`" begin
    x = LinearInterpolation(rand(20))
    y = LinearInterpolation(rand(20))
    t = range(0, stop = 1, length=100)
    for symmetric in (true, false), max_iters in (1, 2, 5)
        cost, ϕ, ψ = gdtw(x, y; symmetric = symmetric, max_iters = max_iters)

        data = prepare_gdtw(x, y; symmetric = symmetric)
        local cost2
        for i = 1:max_iters
            cost2 = iterative_gdtw!(data, max_iters = i)
        end
        ϕ2, ψ2 = gdtw_warpings(data)
        @test cost ≈ cost2
        @test ϕ.(t) ≈ ϕ2.(t)
        @test ψ.(t) ≈ ψ2.(t)

        @info "The following warning is intentional"
        c = iterative_gdtw!(data, max_iters = 1)
        @test c == 0
    end
end


text:
@doc raw"""
    HardCoreGraph{T<:Integer} <: AbstractGraphPointProcess{T}

Concrete type representing a point process on the vertices of a `graph` ``=(V, E)`` parametrized by ``\beta \geq 0`` which characterizes the distribution on the [independent sets](https://en.wikipedia.org/wiki/Independent_set_(graph_theory)) of `graph`, where each vertex is present with marginal probability ``\frac{\beta}{1+\beta}``.

In other words, it can also be viewed as the product distribution ``\operatorname{Bernoulli}(\frac{\beta}{1+\beta})^{\otimes |V|}`` on the vertices of `graph` conditioned on forming an independent set,

```math
    \mathbb{P}\!\left[ \mathcal{X} = X \right]
    \propto
    \prod_{x\in X}
        \frac{\beta}{1+\beta}
     1_{X \text{forms an independent set}}.
```

**See also**

- Section 7.2 of [GuJeLi19](@cite),
- Example 4.1 of [MoKr20](@cite),
- [`PRS.HardCorePointProcess`](@ref), the spatial counterpart of [`PRS.HardCoreGraph`](@ref).

# Example

A realization from a ``5\times 5`` grid graph, with ``\beta = 0.5``.

![assets/hard_core_graph.png](assets/hard_core_graph.png)
"""
struct HardCoreGraph{T<:Integer} <: AbstractGraphPointProcess{T}
    "Graph"
    graph::LG.SimpleGraph{T}
    β::Float64
end

function Base.show(io::IO, pp::HardCoreGraph{T}) where {T}
    print(io, "HardCoreGraph{$T}\n- graph = $(pp.graph)\n- β = $(pp.β)")
end

"""
    HardCoreGraph(
        graph::LG.SimpleGraph{T},
        β::Real
    ) where {T<:Integer}

Construct a [`PRS.HardCoreGraph`](@ref).

```jldoctest; output = true
using PartialRejectionSampling
using LightGraphs; const LG = LightGraphs

g, β = LG.grid([5, 5]), 1
PRS.HardCoreGraph(g, β)

# output

HardCoreGraph{Int64}
- graph = {25, 40} undirected simple Int64 graph
- β = 1.0
```
"""
function HardCoreGraph(
    graph::LG.SimpleGraph{T},
    β::Real
) where {T<:Integer}
    @assert β >= 0
    return HardCoreGraph{T}(graph, β)
end

# Sampling

"""
    generate_sample(
        [rng=Random.AbstractRNG,]
        pp::HardCoreGraph{T}
    )::Vector{T} where {T}

Generate an exact sample from the [`PRS.SinkFreeGraph`](@ref).

Default sampler is [`PRS.generate_sample_prs`](@ref).
"""
function generate_sample(
    rng::Random.AbstractRNG,
    pp::HardCoreGraph{T}
)::Vector{T} where {T}
    return generate_sample_prs(rng, pp)
end

# generate_sample(pp::HardCoreGraph) = generate_sample(Random.default_rng(), pp)

## Partial Rejeciton Sampling (PRS)

"""
    generate_sample_prs(
        [rng=Random.AbstractRNG,]
        pp::HardCoreGraph{T}
    )::Vector{T} where {T}

Sample from [`PRS.HardCoreGraph`](@ref) using Partial Rejection Sampling (PRS), see Section 7.2 of [GuJeLi19](@cite)

**See also**

- Example 4.1 of [MoKr20](@ref).

# Example

An illustration of the procedure on a ``5\\times 5`` grid graph.

- red: variables involved in constraints violation (gray neighboring nodes)
- orange: variables to be resampled (red nodes and their neighborhood)

![assets/hard_core_graph.gif](assets/hard_core_graph.gif)
"""
function generate_sample_prs(
    rng::Random.AbstractRNG,
    pp::HardCoreGraph{T}
)::Vector{T} where {T}

    β_max = 1 / (2 * sqrt(exp(1)) * LG.Δ(pp.graph) - 1)
    if pp.β > β_max
        @warn "The arguments do not satisfy Theorem 35 of Guo, Jerrum and Liu (2019).\nPartial rejection sampling may not be efficient.\nGiven `graph`, consider choosing β ≤ $(β_max)"
    end

    proba = pp.β / (one(pp.β) + pp.β)

    adj = LG.adjacency_matrix(pp.graph)
    occupied = randsubseq(rng, LG.vertices(pp.graph), proba)
    while true
        # Check if occupied vertices form an independent set
        sub_graph = LG.SimpleGraph(adj[occupied, occupied])
        LG.ne(sub_graph) == 0 && break

        independent, resample = T[], Set{T}()
        for cc in LG.connected_components(sub_graph)
            if length(cc) == 1  # Identify current independent vertices
                append!(independent, occupied[cc])
            else  # Construct the resampling set of vertices
                union!(resample, occupied[cc])
                for v in cc
                    union!(resample, LG.neighbors(pp.graph, occupied[v]))
                end
            end
        end
        randsubseq!(rng, occupied, collect(resample), proba)
        append!(occupied, independent)
    end
    return occupied
end

# generate_sample_prs(pp::HardCoreGraph) = generate_sample_prs(Random.default_rng(), pp)


text:
using Parameters, StaticArrays, BenchmarkTools, LinearAlgebra, Interpolations
# using MAT
using TrajectoryOptimization
import TrajectoryOptimization: dynamics, AbstractConstraint, evaluate, state_dim
const TO = TrajectoryOptimization

@with_kw mutable struct BicycleCar{T,P} <: AbstractModel
    path::P = StraightPath(30.0, 0.)
    a::T = 1.2169  # dist to front axle (m)
    b::T = 1.4131  # dist to rear axle (m)
    h::T = 0.5       # height of center of gravity? (m)
    mass::T = 1744.  # mass (kg)
    Iz::T = 2.7635e3  # moment of inertia (kg⋅m²)
    Cαf::T = 66000.  # (N/rad)
    Cαr::T = 60000.  # (N/rad)
    μ::T = 0.9       # coefficient of friction
    Cd0::T = 218.06  # Cd0 (N)
    Cd1::T = 0.0 # (N/mps)
    g::T = 9.81  # gravity (m/s²)
end

labels(::BicycleCar) = ["delta" "fx" "r" "Uy" "Ux" "dpsi" "e" "t"]

Base.size(::BicycleCar) = 8,2

function dynamics(car::BicycleCar, x, u, s)
    # Scale the dynamics
    Sx, Su = scaling(car)
    x = Sx*x
    u = Su*u

    δ_dot  = u[1]  # steering angle rate
    fx_dot = u[2]  # accel rate
    δ  = x[1]  # steering angle
    fx = x[2]  # acceleration
    r  = x[3]  # yaw rate
    Uy = x[4]  # lateral velocity
    Ux = x[5]  # longitudinal velocity
    Δψ = x[6]  # heading error
    e  = x[7]  # lateral error
    t  = x[8]  # time

    # Road curvature
    k = curvature(car.path, s)

    # Tire forces
    α_f = atan(Uy + car.a * r, Ux) - δ
    α_r = atan(Uy - car.b * r, Ux)

    Fxf, Fxr, Fzf, Fzr = FWD_force_model(car, fx)  # long. and downward forces
    Fyf = fiala_tire_model(car.μ, car.Cαf, α_f, Fxf, Fzf)  # front lateral force
    Fyr = fiala_tire_model(car.μ, car.Cαr, α_r, Fxr, Fzr)  # rear  lateral force

    # # Drag Force
    Fx_drag = 0.0
    #
    # # Spatial derivative
    s_dot = (Ux * cos(Δψ) - Uy * sin(Δψ)) / (1 - k*e)
    #
    # # Get tire forces
    # Fxf, Fxr, Fyf, Fyr = logit_lateral_force_model(car, x)
    #
    # # State Derivatives
    r_dot  = (car.a * (Fyf * cos(δ) + Fxf*sin(δ)) - car.b * Fyr) / car.Iz
    Ux_dot =  r * Uy + (Fxf * cos(δ) - Fyf * sin(δ) + Fxr - Fx_drag) / car.mass
    Uy_dot = -r * Ux + (Fyf * cos(δ) + Fxf * sin(δ) + Fyr) / car.mass
    Δψ_dot = r - k * s_dot
    e_dot = Ux * sin(Δψ) + Uy * cos(Δψ)
    t_dot = 1 / s_dot

    # Slip angles
    αf = atan(Uy + car.a*r, Ux) - δ
    αr = atan(Uy - car.b*r, Ux)

    # # Get longitudinal forces
    # Fxf, Fxr, Fzf, Fzr = FWD_force_model(car, fx)

    # # Get lateral forces
    # Fyf = logit_tire_model(αf, car.μ, car.Cαf, Fxf, Fzf)
    # Fyr = logit_tire_model(αr, car.μ, car.Cαr, Fxr, Fzr)
    #
    # # Dynamics from paper
    # s_dot = Ux - Uy*Δψ
    # e_dot = Uy + Ux*Δψ
    #
    # Δψ_dot = r - k*Ux
    # Ux_dot = (Fxf + Fxr)/car.mass + r*Uy
    # Uy_dot = (Fyf + Fyr)/car.mass - r*Ux
    # r_dot = (car.a*Fyf - car.b*Fyr)/car.Iz
    # t_dot = 1/s_dot

    return t_dot * @SVector [δ_dot, fx_dot, r_dot, Uy_dot, Ux_dot, Δψ_dot, e_dot, 1]
end

function scaling(model::BicycleCar)
    Sx = Diagonal(@SVector [1,1e3,1,1,1,1,1,1])
    Su = Diagonal(@SVector [1,1e2])
    return Sx, Su
end

function FWD_force_model(car::BicycleCar, fx)
    l = car.a + car.b
    Fxf = max(fx * car.b / l, fx)  # lower bound
    Fxr = min(fx * car.a / l, 0.0)  # upper bound (can't accelerate)
    # Fxf = fx*car.b / l
    # Fxr = fx*car.a / l
    Fzf = car.mass * car.g * car.b / l
    Fzr = car.mass * car.g * car.a / l
    return Fxf, Fxr, Fzf, Fzr
end

function logit_lateral_force_model(car::BicycleCar, x)
    δ = x[1]
    fx = x[2]
    r  = x[3]
    Uy = x[4]
    Ux = x[5]

    α_f = atan(Uy + car.a * r, Ux) - δ
    α_r = atan(Uy - car.b * r, Ux)

    Fxf, Fxr, Fzf, Fzr = FWD_force_model(car, fx)
    # Fyf = logit_tire_model(car.μ, car.Cαf, α_f, Fxf, Fzf)
    # Fyr = logit_tire_model(car.μ, car.Cαr, α_r, Fxr, Fzr)
    Fyf = fiala_tire_model(car.μ, car.Cαf, α_f, Fxf, Fzf)
    Fyr = fiala_tire_model(car.μ, car.Cαr, α_r, Fxr, Fzr)
    return Fxf, Fxr, Fyf, Fyr
end

function logit_tire_model(α, μ, Cα, Fx, Fz)
    Fy_max = sqrt((μ*Fz)^2 - Fx^2)
    slope = 2*Cα / Fy_max
    Fy = Fy_max * (1-2*exp(slope*α)) / (1 + exp(slope*α))
    return Fy
end

function fiala_tire_model(μ, Cα, α, Fx, Fz)
    Fx = min( 0.99*μ*Fz, Fx)  # try removing this?
    Fx = max(-0.99*μ*Fz, Fx)
    Fy_max = sqrt((μ*Fz)^2 - Fx^2)
    a_slide = atan(3*Fy_max, Cα)
    tan_a = tan(α)
    Fy_unsat = ( -Cα * tan_a + Cα^2/(3*Fy_max)*tan_a*abs(tan_a) - Cα^3/(27*Fy_max^2) * tan_a^3 )
    Fy_sat = -Fy_max*sign(α)
    if abs(α) < a_slide
        return Fy_unsat
    else
        return Fy_sat
    end
end

function load_scenario!(car::BicycleCar, s, k, s_interp)
    itp = CubicSplineInterpolation(s,k)
    k_interp = itp.(s_interp)
    car.k_c = Dict(s_interp .=> k_interp)
    return nothing
end

struct BrakeForceConstraint{L} <: AbstractConstraint{Inequality,State,1}
    car::L
end

state_dim(con::BrakeForceConstraint) = size(con.car)[1]

function evaluate(con::BrakeForceConstraint, x::SVector)
    Sx, = scaling(con.car)
    x = Sx*x
    car = con.car
    l = car.a + car.b

    Fx = x[2]
    r  = x[3]
    Uy = x[4]
    Ux = x[5]

    # Down force on front
    αf = atan(Uy + car.a*r, Ux)
    Fzf = car.mass*car.g*car.b/l

    # Long. force on front
    Fxf = Fx*car.b/l
    SVector{1}(-cos(αf)*car.μ*Fzf - Fxf)
end

struct ContingencyCar{T} <: AbstractModel
    cars::Vector{BicycleCar{T}}
    ix::Vector{SVector{8,Int}}
    iu::Vector{SVector{2,Int}}
end

function ContingencyCar(num_cars=2)
    cars = [BicycleCar() for i = 1:num_cars]
    n,m = 8,2
    ix = [SVector{n}((1:n) .+ (i-1)*n) for i = 1:num_cars]
    iu = [SVector{m}((1:m) .+ (i-1)*m) for i = 1:num_cars]
    ContingencyCar(cars, ix, iu)
end

Base.size(::ContingencyCar) = 16,4

function dynamics(car::ContingencyCar, x, u, s)
    ix, iu = car.ix, car.iu
    cars = car.cars
    x1,u1 = x[ix[1]], u[iu[1]]
    ẋ1 = dynamics(cars[1], x1, u1, s)
    x2,u2 = x[ix[2]], u[iu[2]]
    ẋ2 = dynamics(cars[2], x1, u1, s)
    return [ẋ1; ẋ2]
end

function localToGlobal(path::AbstractPath, Z::Traj)
    e = [z.z[7] for z in Z]
    s = [z.t for z in Z]
    localToGlobal(path, s, e)
end

function Plots.plot(solver::TrajectoryOptimization.AbstractSolver)
    plot(get_model(solver), get_trajectory(solver))
end

function Plots.plot(car::BicycleCar, Z::Traj)
    path = car.path
    plot(path, aspect_ratio=:equal)
    x,y = localToGlobal(path, Z)
    plot!(x,y, linewidth=2)
end


text:
"""
ilqr_solve.jl
"""

"iLQR solve method"
function solve!(solver::iLQRSolver{IR}) where {IR}
    # initialize
    reset!(solver)
    set_verbosity!(solver)
    clear_cache!(solver)
    solver.ρ[1] = solver.opts.bp_reg_initial
    solver.dρ[1] = 0.0
    n, m, N = size(solver)
    X = solver.X
    U = solver.U
    ts = solver.ts

    # initial rollout and cost calculation
    J = J_prev = 0.
    for k = 1:N-1
        J_prev += cost(solver.obj, X, U, k)
        dt = ts[k + 1] - ts[k]
        discrete_dynamics!(X[k + 1], IR, solver.model, X[k], U[k], ts[k], dt)
    end
    J_prev += cost(solver.obj, X, U, N)

    # run iLQR iterations
    for i = 1:solver.opts.ilqr_max_iterations
        # step
	    backwardpass!(solver)
        J, reg_flag = forwardpass!(solver, J_prev)
        # exit if solve succeeded
        if solver.stats.status > SOLVE_SUCCEEDED
            break
        end
        # accept the updated trajectory
        for k = 1:N-1
            solver.X[k] .= solver.X_tmp[k]
            solver.U[k] .= solver.U_tmp[k]
        end
        solver.X[N] .= solver.X_tmp[N]
        # record iteration
        dJ = abs(J - J_prev)
        J_prev = J
        gradient_todorov!(solver)
        record_iteration!(solver, J, dJ)
        if is_verbose(solver)
            print_level(InnerLoop, global_logger())
        end
        # exit if converged
        exit = evaluate_convergence(solver, i, reg_flag)
        exit && break
    end
    terminate!(solver)
    return J
end

"""
$(SIGNATURES)
Simulate the system forward using the optimal feedback gains from the backward pass,
projecting the system on the dynamically feasible subspace. Performs a line search to ensure
adequate progress on the nonlinear problem.
"""
function forwardpass!(solver::iLQRSolver, J_prev)
    ΔV = solver.ΔV
    J = Inf
    α = 1.0
    iter = 0
    z = -1.0
    expected = 0.0
    rollout_flag = false
    reg_flag = false

    while ((z ≤ solver.opts.line_search_lower_bound
            || z > solver.opts.line_search_upper_bound)
           && J >= J_prev)
        # exit and regularize if the maximum number of line search decrements
        # has occured
        if iter > solver.opts.iterations_linesearch
            J = J_prev
            z = 0
            α = 0.0
            expected = 0.0
            regularization_update!(solver, :increase)
            solver.ρ[1] += solver.opts.bp_reg_fp
            reg_flag = true
            break
        end
        # otherwise, rollout a new trajectory for current alpha
        J, rollout_flag = rollout!(solver, α)
        # reduce step size if rollout returns non-finite values (NaN or Inf)
        if rollout_flag
            J = J_prev
            iter += 1
            α /= 2.0
            continue
        end
        # update
        expected = -α * (ΔV[1] + α * ΔV[2])
        if expected > 0.0
            z  = (J_prev - J)/expected
        else
            z = -1.0
        end
        iter += 1
        α /= 2.0
    end

    if J > J_prev
        # error("Error: Cost increased during Forward Pass")
        solver.stats.status = COST_INCREASE
        return NaN
    end

    @logmsg InnerLoop :expected value=expected
    @logmsg InnerLoop :z value=z
    @logmsg InnerLoop :α value=2*α
    @logmsg InnerLoop :ρ value=solver.ρ[1]

    return J, reg_flag
end

"""
Stash iteration statistics
"""
function record_iteration!(solver::iLQRSolver, J, dJ)
    gradient = mean(solver.grad)
    record_iteration!(solver.stats, cost=J, dJ=dJ, gradient=gradient)
    i = solver.stats.iterations::Int

    if dJ ≈ 0
        solver.stats.dJ_zero_counter += 1
    else
        solver.stats.dJ_zero_counter = 0
    end

    @logmsg InnerLoop :iter value=i
    @logmsg InnerLoop :cost value=J
    @logmsg InnerLoop :dJ   value=dJ
    @logmsg InnerLoop :grad value=gradient
    # @logmsg InnerLoop :zero_count value=solver.stats[:dJ_zero_counter][end]
    return nothing
end

"""
$(SIGNATURES)
    Calculate the problem gradient using heuristic from iLQG (Todorov) solver
"""
function gradient_todorov!(solver::iLQRSolver)
    N = solver.N
    tmp = solver.U_tmp[N]
    u = solver.U_tmp[N+1]
    for k in eachindex(solver.d)
	tmp .= abs.(solver.d[k])
	u .= abs.(solver.U[k]) .+ 1
	tmp ./= u
	solver.grad[k] = maximum(tmp)
    end
end

"""
$(SIGNATURES)
Check convergence conditions for iLQR
"""
function evaluate_convergence(solver::iLQRSolver, ilqr_iterations::Int, reg_flag::Bool)
    # Get current iterations
    i = solver.stats.iterations
    grad = solver.stats.gradient[i]
    dJ = solver.stats.dJ[i]
    J = solver.stats.cost[i]

    # If the change in cost is small, the gradient is small, and
    # the last step was not a regularization step, we have converged
    if ((0.0 <= dJ < solver.opts.ilqr_ctol) && (grad < solver.opts.ilqr_gtol)
        && !reg_flag)
        @logmsg InnerLoop "Cost criteria satisfied."
        solver.stats.status = SOLVE_SUCCEEDED
        return true
    end

    # Check total iterations
    if i >= solver.opts.iterations
        @logmsg InnerLoop "Hit max iterations. Terminating."
        solver.stats.status = MAX_ITERATIONS
        return true
    end

    # Exit iLQR and continue solving if the maximum number of iLQR iterations
    # has been reached.
    if ilqr_iterations >= solver.opts.ilqr_max_iterations
        return true
    end

    # Outer loop update if forward pass is repeatedly unsuccessful
    if solver.stats.dJ_zero_counter > solver.opts.dJ_counter_limit
        @logmsg InnerLoop "dJ Counter hit max. Terminating."
        solver.stats.status = NO_PROGRESS
        return true
    end

    # Terminate if max cost value was reached.
    if J > solver.opts.max_cost_value
        @logmsg InnerLoop "Hit max cost value. Terminating."
        solver.stats.status = MAXIMUM_COST
        return true
    end

    return false
end

"""
$(SIGNATURES)
Update the regularzation for the iLQR backward pass
"""
function regularization_update!(solver::iLQRSolver,status::Symbol=:increase)
    # println("reg $(status)")
    if status == :increase # increase regularization
        # @logmsg InnerLoop "Regularization Increased"
        solver.dρ[1] = max(solver.dρ[1] * solver.opts.bp_reg_increase_factor,
                           solver.opts.bp_reg_increase_factor)
        solver.ρ[1] = max(solver.ρ[1] * solver.dρ[1], solver.opts.bp_reg_min)
        # if solver.ρ[1] > solver.opts.bp_reg_max
        #     @warn "Max regularization exceeded"
        # end
    elseif status == :decrease # decrease regularization
        # TODO: Avoid divides by storing the decrease factor (divides are 10x slower)
        solver.dρ[1] = min(solver.dρ[1] / solver.opts.bp_reg_increase_factor,
                           1. / solver.opts.bp_reg_increase_factor)
        solver.ρ[1] = solver.ρ[1] * solver.dρ[1] * (solver.ρ[1] * solver.dρ[1]
                                                    > solver.opts.bp_reg_min)
    end
end


text:
module OptionalUnits

using Unitful
using MacroTools

export @optionalunits

"""
Can be applied to function and struct definitions that use Unitful types to add default units, when the parameters that are passed to the function do not have a unit attached.

In a function definition the type parameter e.g. Unitful.Length is appended with the default unit (e.g. x::Unitful.Length→u"mm")

In a struct definition nothing additional has to be added, the field types with Unitful.Quantity will be detected

# Examples
```julia
julia> using OptionalUnits, Unitful

julia>  @optionalunits function addOneMeter(x::Unitful.Length→u"m")
            return x+1u"m"
        end

julia>  @optionalunits function addOneMeter(x::Vector{Unitful.Length→u"m"})
            return x.+1u"m"
        end

julia> @optionalunits struct Point
            x::typeof(1.0u"m")
            y::typeof(1.0u"m")
        end

julia>  Point(1,2)
```
"""
macro optionalunits(ex)

    @assert isa(ex, Expr)

    function create_unit_check(s::Symbol, u::Unitful.Unitlike, isArray::Bool)
        if isArray
            return quote
                dims = dimension.($s)
                @assert all(dims .== [first(dims)]) "The array-like type $($s) has mixed dimensions which is not supported by @optionalunits"
                if first(dims) == NoDims
                    @warn "Used default unit $($u) for array-like type" $s
                    $s *= $u
                end
            end
        else
            return quote
                if dimension($s) == NoDims
                    @warn "Used default unit $($u) for" $s
                    $s *= $u
                end
            end
        end
    end

    if ex.head == :function
        unitpairs = Vector{Tuple{Symbol,Unitful.Unitlike,Expr,Bool}}()
        ex = MacroTools.postwalk(ex) do x
            if @capture(x, name_::container_{par__})

                par_new = copy(par)
                for i in range(1, stop=length(par))
                    if @capture(par[i], dimension_ → unit_)
                        @assert eval(container) <: AbstractArray "Type $(container) is not supported for the use with optional units, only implementations of the AbstractArray interface are supported"
                        if length(unitpairs) != 0 && unitpairs[end][1] == name
                            error("Found two default units in type parameters for single argument $(name)")
                        end

                        @assert typeof(eval(unit)) <: Unitful.Unitlike "Annotated unit $unit is not a Unitful.Unitlike"
                        @assert (1 * eval(unit) isa eval(dimension)) "Annotated unit $unit does not have the correct dimension $dimension"

                        par[i] = :(<:$dimension)
                        par_new[i] = :(<:Real)
                        push!(unitpairs, (name, eval(unit), dimension, true))
                    end
                end

                return :($name::Union{$container{$(par...)},$container{$(par_new...)}})
            end
            if @capture(x, name_::dimension_ → unit_)
                @assert typeof(eval(unit)) <: Unitful.Unitlike "Annotated unit $unit is not a Unitful.Unitlike"
                @assert (1 * eval(unit) isa eval(dimension)) "Annotated unit $unit does not have the correct dimension $dimension"

                push!(unitpairs, (name, eval(unit), dimension, false))
                return :($name::Union{$dimension,Real})
            end

            return x
        end

        for pairs in unitpairs
            insert!(ex.args[2].args, 1, create_unit_check(pairs[1], pairs[2], pairs[4]))
        end

        return esc(prettify(ex))

    elseif ex.head == :struct
        typename = ex.args[2]
        eval(ex) # evaluate struct in local scope, to be able to execute Base.fieldnames, there has to be a better way to get the fieldnames from the expression without evaluating
        type = eval(typename)
        fieldnames = Base.fieldnames(type)

        out = Expr(:block)
        push!(out.args, ex)

        tmp = :(function $(typename)($(fieldnames...)) end)

        for i in range(1, stop=length(fieldnames))

            fieldtype = Core.fieldtype(type, i)

            if fieldtype <: Unitful.Quantity
                push!(tmp.args[2].args, create_unit_check(fieldnames[i], unit(fieldtype), false))
            elseif fieldtype <: AbstractArray && eltype(fieldtype) <: Unitful.Quantity
                push!(tmp.args[2].args, create_unit_check(fieldnames[i], unit(eltype(fieldtype)), true))
            end

            push!(tmp.args[2].args, :($(fieldnames[i]) = Base.convert(Core.fieldtype($typename, $i), $(fieldnames[i]))))

        end

        push!(tmp.args[2].args, :(return $(typename)($(fieldnames...))))

        push!(out.args, tmp)
        return esc(out)
    end
end

end


text:
export VertexIter, EliminatedVertexIter, eliminated_vertices, iseliminated
export NeighborIter

"""
    VertexIter{GT}<:AbstractArray{Int,1}

Vertex enumerator for a graph.
"""
struct VertexIter{GT}<:AbstractArray{Int,1}
    eg::GT
end
Base.length(vs::VertexIter) = vs.eg.nv
Base.getindex(vs::VertexIter, i::Int) = vs.eg.vertices[end-vs.eg.nv+i]

"""
    EliminatedVertexIter{GT}<:AbstractArray{Int,1}

Eliminated vertex enumerator for a graph.
"""
struct EliminatedVertexIter{GT}<:AbstractArray{Int,1}
    eg::GT
end
Base.length(vs::EliminatedVertexIter) = nv0(vs.eg)-vs.eg.nv
Base.getindex(vs::EliminatedVertexIter, i::Int) = vs.eg.vertices[i]

"""eliminated vertices of a `EliminateGraph`."""
eliminated_vertices(eg::EliminateGraph) = EliminatedVertexIter(eg)

"""
    iseliminated(eg::EliminateGraph, i::Int) -> Bool

Return true if a vertex of a `EliminateGraph` is eliminated.
"""
iseliminated(eg::EliminateGraph, i::Int) = i in eliminated_vertices(eg)

for V in [:VertexIter, :EliminatedVertexIter]
    @eval Base.size(vs::$V) = (length(vs),)
    @eval Base.size(vs::$V, i::Int) = i==1 ? length(vs) : 1
    @eval Base.IteratorEltype(::Type{$V}) = Base.HasEltype()
    @eval Base.IteratorSize(::Type{$V}) = Base.HasLength()
    @eval Base.eltype(vs::$V) = Int
    @eval Base.eltype(vs::Type{$V}) = Int

    VI = V==:VertexIter ? :(eg.vertices[end-eg.nv+state]) : :(eg.vertices[state])
    @eval function Base.iterate(vs::$V, state=1)
        eg = vs.eg
        if state > length(vs)
            return nothing
        else
            @inbounds vi = $VI
            return vi, state+1
        end
    end
end

"""vertices of a graph."""
vertices(eg::EliminateGraph) = VertexIter(eg)

"""
    NeighborIter{SP,NTH,GT}<:AbstractArray{Int,1}

Neighbor enumerator for a graph.
"""
struct NeighborIter{SP,NTH,GT}
    eg::GT
    i::Int
end

NeighborIter{SP,NTH}(eg::GT, i::Int) where {SP,NTH,GT} = NeighborIter{SP,NTH,GT}(eg, i)

Base.IteratorEltype(::Type{NeighborIter}) = Base.HasEltype()
Base.eltype(vs::NeighborIter) = Int
Base.eltype(vs::Type{NeighborIter}) = Int

Base.:(==)(ni::NeighborIter, v) = Set([ni...]) == Set([v...])
Base.:(==)(v, ni::NeighborIter) = Set([v...]) == Set([nj...])
Base.:(==)(ni::NeighborIter, nj::NeighborIter) = Set([ni...]) == Set([nj...])

function Base.show(io::IO, ni::NeighborIter)
    println(io, summary(ni))
    print(io, (ni...,))
end

@generated function Base.iterate(ni::NeighborIter{SP,1}, state=1) where SP
    condition = SP == CLOSED ? :(vj==ni.i || isconnected(eg,vj,ni.i)) : :(isconnected(eg,vj,ni.i))
    quote
        eg = ni.eg
        if state > nv(eg)
            return nothing
        end
        for j = state:nv(eg)
            vj = @inbounds vertices(eg)[j]
            $condition && return (vj, j+1)
        end
        return nothing
    end
end

#neighbors(eg::EliminateGraph, i::Int) = generate_set(eg, NearestNeighbors{OPEN}(i))
neighbors(eg::EliminateGraph, vi::Int) = NeighborIter{OPEN,1}(eg, vi)

#neighborcover(eg::EliminateGraph, i::Int) = generate_set(eg, NearestNeighbors{CLOSED}(i))
neighborcover(eg::EliminateGraph, vi::Int) = NeighborIter{CLOSED,1}(eg, vi)


text:
using Distributed
using TimerOutputs

to = TimerOutput()

if nprocs() == 1
    addprocs(3)
end
# Number of logical CPU cores available in the system.
# println("CPU cores: ", Sys.CPU_CORES)
Sys.cpu_summary()
# Number of available processes.
# println("nprocs: ", nprocs())
# println("nworkers: ", nworkers())
println("nprocs: ", nprocs())
println("nworkers: ", nworkers())

function prepare_data(fib_parameter, n)
    input = rand(1:fib_parameter, n)
    return input
end

# The goal is make suma of fibonaci of every element of input array

@everywhere fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)
# Ordinary for loop.
function forloop(input)
    a = 0
    for i = 1:length(input)
         a += fib(input[i]);
    end
    return a
end

function parallel_forloop(input)
    a = @distributed (+) for i = 1:length(input)
         fib(input[i]);
    end
    return a
end
function parallel_pmap(input)
    pout = pmap(fib, input)
    a = 0
    for i = 1:length(pout)
        a += pout[i]
    end
    return a
end

println("Warming")
input = prepare_data(10, 10)
@timeit to "warming" begin
    @timeit to "forloop" forloop(input)
    @timeit to "parallel forloop" parallel_forloop(input)
    @timeit to "parallel pmap" parallel_pmap(input)
end

println("---")
input = prepare_data(30, 10000)
@timeit to "10000" begin
    @timeit to "forloop" forloop(input)
    @timeit to "parallel forloop" parallel_forloop(input)
    @timeit to "parallel pmap" parallel_pmap(input)
end
input = prepare_data(55, 5)

@timeit to "5" begin
    @timeit to "forloop" forloop(input)
    @timeit to "parallel forloop" parallel_forloop(input)
    @timeit to "parallel pmap" parallel_pmap(input)
end

display(to)
# Parallel for loop
# @time @parallel (+) for i = 1:20000
# @time @distributed (+) for i = 1:n
#     Int(rand(Bool));
# end
#
# # Predefined function.
# @time sum(rand(0:1, n))
#
# # For a not so small amount of work:
# # n = 200000000;
# # n = 20000;
# # @time @parallel (+) for i = 1:n
# @time @distributed (+) for i = 1:n
#     Int(rand(Bool))
# end
#
# DO NOT run sum(rand(0:1, n)) or the ordinary for loop with 'n': You may run
# out of memory.


text:
function ic_eqm(lx::Float64,ly::Float64,nx::Int,ny::Int,Ξ::Float64;jw::Float64=0.05)
    ζ0 = zeros(ComplexF64,2*ny-1,nx)
    ζjet = acoeffs(ly,ny,Ξ,1.0,jw=jw)
    for y in 1:2*ny-1
        ζ0[y,1] = ζjet[y]
    end
    ζ0
end

function ic_rand(lx::Float64,ly::Float64,nx::Int,ny::Int)
    uxy = randn(Float64,2*ny-1,2*nx-1)
    umn = fftshift(fft(uxy))
    umn[ny,nx] = 0.0 + 0.0im
    umn[:,nx:2*nx-1]
end

function ic_pert_eqm(lx::Float64,ly::Float64,nx::Int,ny::Int,Ξ::Float64;jw::Float64=0.05)
    @. ic_eqm(lx,ly,nx,ny,Ξ,jw=jw) + ic_rand(lx,ly,nx,ny)*1e-4
end

function ic_cumulants(nx::Int,ny::Int,Λ::Int,u0::Array{ComplexF64,2})
    u0_low::Array{ComplexF64,2} = u0[:,1:Λ+1]
    for n = 1:ny-1
        u0_low[n,1] = conj(u0_low[2*ny - n,1])
    end
    u0_low[ny,1] = 0.0 + im*0.0

    # u0_high = zeros(ComplexF64,2*ny-1,nx-Λ-1,2*ny-1,nx-Λ-1)
    u0_high::Array{ComplexF64,4} = zeros(ComplexF64,2*ny-1,nx-Λ,2*ny-1,nx-Λ)
    for m1=Λ+1:nx-1
        for n1=-(ny-1):ny-1
            for m2=Λ+1:nx-1
                for n2=-(ny-1):ny-1

                    u0_high[n2+ny,m2-Λ,n1+ny,m1-Λ] = u0[n2+ny,m2+1]*conj(u0[n1+ny,m1+1])

                end
            end
        end
    end
    ArrayPartition(u0_low,u0_high)
end

function ic_cumulants(nx::Int,ny::Int,Λ::Int,σ::Float64,u0::Array{ComplexF64,2})
    u0_low::Array{ComplexF64,2} = u0[:,1:Λ+1]
    for n = 1:ny-1
        u0_low[n,1] = conj(u0_low[2*ny - n,1])
    end
    u0_low[ny,1] = 0.0 + im*0.0

    twopoint = σ*Matrix{ComplexF64}(I,(2*ny-1)*(nx-Λ),(2*ny-1)*(nx-Λ))
    u0_high = reshape(twopoint,2*ny-1,nx-Λ,2*ny-1,nx-Λ)

    ArrayPartition(u0_low,u0_high)
end


text:
# ---
# title: 1227. Airplane Seat Assignment Probability
# id: problem1227
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Dynamic Programming, Brainteaser
# link: <https://leetcode.com/problems/airplane-seat-assignment-probability/description/>
# hidden: true
# ---
#
# `n` passengers board an airplane with exactly `n` seats. The first passenger
# has lost the ticket and picks a seat randomly. But after that, the rest of
# passengers will:
#
#   * Take their own seat if it is still available,
#   * Pick other seats randomly when they find their seat occupied
#
# What is the probability that the n-th person can get his own seat?
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1
#     Output: 1.00000
#     Explanation: The first person can only get the first seat.
#
# **Example 2:**
#
#
#
#     Input: n = 2
#     Output: 0.50000
#     Explanation: The second person has a probability of 0.5 to get the second seat (when first person gets the first seat).
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file was generated, do not modify it. # hide
using Pkg # hide
Pkg.status(["SimilaritySearch", "Plots"]) # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
struct PropaneRefConsts <: EoSParam
    Mw::Float64
    T_c::Float64
    P_c::Float64
    rho_c::Float64
    v_c::Float64
    n::Vector{Float64}
    t::Vector{Float64}
    d::Vector{Int}
    l::Vector{Int}
    beta::Vector{Float64}
    gamma::Vector{Float64}
    eta::Vector{Float64}
    epsilon::Vector{Float64}
    function PropaneRefConsts()
        Mw = 44.09562 #g·mol-1
        T_c = 369.89    #K
        P_c = 4.2512e6 #Pa
        rho_c= 5000.0 # mol·m-3
        v_c = 1/rho_c
        n = [0.042910051, 1.7313671, -2.4516524, 0.34157466, -0.46047898, -0.66847295, 0.20889705, 0.19421381, -0.22917851, -0.60405866, 0.066680654, 0.017534618, 0.33874242, 0.22228777, -0.23219062, -0.09220694, -0.47575718, -0.017486824]
        t = [1.0, 0.33, 0.8, 0.43, 0.9, 2.46, 2.09, 0.88, 1.09, 3.25, 4.62, 0.76, 2.5, 2.75, 3.05, 2.55, 8.4, 6.75]
        d = [4,1,1,2,2,1,3,6,6,2,3,1,1,1,2,2,4,1]
        l = [1,1,1,1,2,2]
        η = [0.963,1.977,1.917,2.307,2.546,3.28,14.6]
        β = [2.33,3.47,3.15,3.19,0.92,18.8,547.8]
        γ = [0.684,0.829,1.419,0.817,1.500,1.426,1.093]
        ε = [1.283,0.6936,0.788,0.473,0.8577,0.271,0.948]
        return new(Mw,T_c,P_c,rho_c,v_c,
        n,t,d,l,β,γ,η,ε)
    end
end

struct PropaneRef <: EmpiricHelmholtzModel
    components::Vector{String}
    consts::PropaneRefConsts
    references::Vector{String}
end

PropaneRef() = PropaneRef(["propane"],PropaneRefConsts(),["1021/je900217v"])

is_splittable(::PropaneRef) = false

function _f0(::PropaneRef,δ,τ)
    δ,τ = promote(δ,τ)
    _1 = one(δ)
    a₁ = -4.970583
    a₂ = 4.29352
    b = (1.062478, 3.344237,5.363757,11.762957)
    v = (3.043,5.874,9.337,7.922)
    α₀ = log(δ) + 3*log(τ) + a₁ + a₂*τ +
     v[1]*log(_1-exp(-b[1]*τ)) +
     v[2]*log(_1-exp(-b[2]*τ)) +
     v[3]*log(_1-exp(-b[3]*τ)) +
     v[4]*log(_1-exp(-b[4]*τ))
     return α₀
end

function _fr1(model::PropaneRef,δ,τ)
    δ,τ = promote(δ,τ)
    n = model.consts.n
    t = model.consts.t
    d = model.consts.d
    l = model.consts.l
    η = model.consts.eta
    β = model.consts.beta
    γ = model.consts.gamma
    ε = model.consts.epsilon
    αᵣ = zero(δ)

    for k in 1:5
        αᵣ += n[k]*(δ^d[k])*(τ^t[k])
    end

    for k in 6:11
        αᵣ += n[k]*(δ^d[k])*(τ^t[k])*exp(-δ^l[k-5])
    end

    for k in 12:18
        _k = k - 11
        mul = exp(-η[_k]*(δ-ε[_k])^2  - β[_k]*(τ-γ[_k])^2)
        αᵣ += n[k]*(δ^d[k])*(τ^t[k])*mul
    end
    return αᵣ
end
#ancillary equations for calculation of P_sat, rhovsat y rholsat
function _propaneref_tsat(T)
    T_c = 369.89
    P_c = 4.2512e6
    T>T_c && return zero(T)/zero(T)
    Tr = T/T_c
    θ = 1.0-Tr
    lnPsatPc = (-6.7722*θ + 1.6938*θ^1.5 -1.3341*θ^2.2 -3.1876*θ^4.8 + 0.94937*θ^6.2)/Tr
    Psat = exp(lnPsatPc)*P_c
    return Psat
end

function _propaneref_rholsat(T)
    T_c = 369.89
    ρ_c = 5000.0
    T>T_c && return zero(T)/zero(T)
    Tr = T/T_c
    θ = 1.0-Tr
    #
    ρ_l = (1.0 + 1.82205*θ^0.345 + 0.65802*θ^0.74 + 0.21109*θ^2.6 + 0.083973*θ^7.2)*ρ_c
    return ρ_l
end

function _propaneref_rhovsat(T)
    T_c = 369.89
    ρ_c = 5000.0
    T>T_c && return zero(T)/zero(T)
    Tr = T/T_c
    θ = 1.0 - Tr
    log_ρ_v_ρ_c = (-2.4887*θ^0.3785 -5.1069*θ^1.07 -12.174*θ^2.7 -30.495*θ^5.5 -52.192*θ^10 -134.89*θ^20)
    ρ_v = exp(log_ρ_v_ρ_c)*ρ_c
    return ρ_v
end

function a_ideal(model::PropaneRef,V,T,z=SA[1.])
    T_c = model.consts.T_c
    rho_c = model.consts.rho_c
    N = only(z)
    rho = (N/V)
    δ = rho/rho_c
    τ = T_c/T
    return  _f0(model,δ,τ)
end

function a_res(model::PropaneRef,V,T,z=SA[1.])
    T_c = model.consts.T_c
    rho_c = model.consts.rho_c
    N = only(z)
    rho = (N/V)
    δ = rho/rho_c
    τ = T_c/T
    return  _fr1(model,δ,τ)
end

function eos(model::PropaneRef, V, T, z=SA[1.0];phase=:unknown)
    R =8.314472
    T_c = model.consts.T_c
    rho_c = model.consts.rho_c
    N = only(z)
    rho = (N/V)
    δ = rho/rho_c
    τ = T_c/T
    return N*R*T*(_f0(model,δ,τ)+_fr1(model,δ,τ))
end

function eos_res(model::PropaneRef,V,T,z=SA[1.0];phase=:unknown)
    R =8.314472
    T_c = model.consts.T_c
    rho_c = model.consts.rho_c
    N = only(z)
    rho = (N/V)
    δ = rho/rho_c
    τ = T_c/T
    return N*R*T*_fr1(model,δ,τ)
end

mw(model::PropaneRef) = SA[model.consts.Mw]
molecular_weight(model::PropaneRef,z = @SVector [1.]) = model.consts.Mw*0.001
T_scale(model::PropaneRef,z=SA[1.0]) = model.consts.T_c
p_scale(model::PropaneRef,z=SA[1.0]) = model.consts.P_c
lb_volume(model::PropaneRef,z=SA[1.0]) = 6.0647250138479785e-5 #calculated at 1000 MPa and 650 K
Base.length(::PropaneRef) = 1
function Base.show(io::IO,mime::MIME"text/plain",model::PropaneRef)
    print(io,"Propane Reference Equation of State")
end
function x0_sat_pure(model::PropaneRef,T,z=SA[1.0])
    log10vv = log10(1.0/_propaneref_rhovsat(T))
    log10vl = log10(1.0/_propaneref_rholsat(T))
    return [log10vl,log10vv]
end

export PropaneRef

#=
    R = 8.314472 #J·mol-1·K-1
    ,Mw = 44.09562 #g·mol-1
    ,T_c = 369.89    #K
    ,P_c = 4.2512e6 #Pa
    ,rho_c= 5000.0 # mol·m-3
    ,T_triple = 85.525 #K
    ,P_Triple = 0.00017 #Pa
    ,rho_v_triple = 2.4e-07 # mol·m-3
    ,rho_l_triple = 16626.0 # mol·m-3
    ,T_0 = 231.036 #K
    ,rho_v_0 = 54.8 # mol·m-3
    ,rho_l_0 = 13173.0 # mol·m-3
    ,ref_T_ideal = 273.15
    ,ref_P_0 = 100000.0
    ,ref_T_ideal_H = 26148.48 ## J·mol-1
    ,ref_T_ideal_S = 157.9105 ##J·mol-1·K-1

=#


text:
module FilteredArrays

using Test
using Gridap.Arrays
using Gridap
using FillArrays

v1 = collect(1:4)
v2 = collect(5:8)
v3 = collect(9:12)
values = [v1,v2,v3]
ptrs = [1,2,3,3,2,2]
a = CompressedArray(values,ptrs)
r = values[ptrs]
test_array(a,r)

filter = [false, true, true, false]
filters = Fill(filter,6)

k = FilterMap()

cache = return_cache(k,filters[1],a[1])

return_type(k,filters[1],a[1])

evaluate!(cache,k,filters[1],a[1])

test_map(a[1][2:3],k,filters[1],a[1])

@test lazy_map(k,filters,a) == [a[i][2:3] for i in 1:6]

#

v1 = [ [2  4  3  3]; [2  4  3  3]]
v2 = [ [5  8  4  3]; [2  4  4  4]]
v3 = [ [2  2  8  8]; [2  9  8  6]]

values = [v1,v2,v3]
ptrs = [1,2,3,3,2,2]
a = CompressedArray(values,ptrs)
r = values[ptrs]
test_array(a,r)

b = Array{Bool,2}(undef,2,4)
b .= true
b[2,3] = b[1,4] = false

filters = Fill(b,6)

r1 = [2, 2, 4, 4, 3, 3]
r2 = [5, 2, 8, 4, 4, 4]
r3 = [2, 2, 2, 9, 8, 6]
res = [r1,r2,r3]
r = CompressedArray(res,ptrs)

@test lazy_map(k,filters,a) == r

end #module


text:
@timedtestset "test rrule for TensorKit.dot" for ix in 1:10
    psi = TensorMap(rand, ComplexF64, ℂ^5*ℂ^2, ℂ^5)
    psi_arr = reshape(psi.data, (5, 2, 5))

    tmp1 = gradient(arr->norm(arr)^2, psi_arr)[1]
    tmp2 = gradient(psi->real(dot(psi, psi)), psi)[1]
    tmp2 = reshape(tmp2.data, (5,2,5))
    @test tmp1 ≈ tmp2
end

@timedtestset "test element-wise multiplication and its rrule" for ix in 1:10
    psi = TensorMap(rand, ComplexF64, ℂ^3*ℂ^2, ℂ^3)
    phi = TensorMap(rand, ComplexF64, ℂ^3*ℂ^2, ℂ^3)
    psi_arr = convert_to_array(psi)
    phi_arr = convert_to_array(phi)

    @test elem_mult(psi, phi) ≈ convert_to_tensormap(psi_arr .* phi_arr, 2)

    function f(arr::Array{ComplexF64, 3})
        a = convert_to_tensormap(arr, 2)
        return norm(elem_mult(a, phi))
    end

    @test f'(psi_arr) ≈ FiniteDifferences.grad(central_fdm(5, 1), f, psi_arr)[1]

end


text:
using Oceananigans.Fields: AbstractField, compute!

fetch_output(output, model, field_slicer) = output(model)

# Needed to support `fetch_output` with `model::Nothing`.
time(model) = model.clock.time
time(::Nothing) = nothing

function fetch_output(field::AbstractField, model, field_slicer)
    compute!(field, time(model))
    return slice_parent(field_slicer, field)
end

convert_output(output, writer) = output
convert_output(output::AbstractArray, writer) = writer.array_type(output)

fetch_and_convert_output(output, model, writer) =
    convert_output(fetch_output(output, model, writer.field_slicer), writer)


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

JLLWrappers.@generate_main_file_header("OpenFOAM")
JLLWrappers.@generate_main_file("OpenFOAM", UUID("89340185-c6b7-517e-a5e0-0a450377a6e4"))
end  # module OpenFOAM_jll


text:
#===============================================================================
tec_to_vtk.jl

This script converts tecplot formatted fehm cont files into the vtk format for
reading into paraview. Cont files must be output in the tecplot format only.
To run, call this script in the folder where the fehmn.files is stored.
The script uses fehmn.files to point to the grid and other outfiles needed for
node locations.

Command to run: julia path/to/directory/tec_to_vtk.jl
Dylan R. Harp (dharp@lanl.gov)
================================================================================

USAGE:

    --sca_diff			Sca file to subtract to create delta variables
    --sca_gdkm_diff		GDKM sca file to subtract to create delta variables
    --var_diff			Variables to create delta variables for
    --output_dir, -o		Folder to create and place vtu2 files
    --mesh_file, -m		Mesh file location
    --run_name, -r		Run name: prefix for visualization file names

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

using WriteVTK
using Glob
using DataFrames
using ArgParse

function get_tec_time(f)
	# Tecplot to dataframe
    fh = open(f)
    # Read in header
    i = 1
    header = true
	tm = ""
    while header
        l = readline(fh)
        if startswith(l,"TITLE")
			title = strip(split(l,"=")[2])
        elseif startswith(l,"VARIABLES")
			variables = split(split(l,"=")[2],"\" \"")
			variables = [replace(v,"\""=>"") for v in variables]
			variables = [strip(replace(v,r"\(.+\)"=>"")) for v in variables]
        elseif startswith(l,"ZONE")
            tm = split(split(l,"\"")[2])[3]
            if occursin("days\"",tm)
				tm = replace(tm,"days\""=>"")
            elseif occursin("days",tm)
				tm = replace(tm,"days"=>"")
			end
        elseif startswith(l,"FILETYPE")
			continue
        elseif startswith(l,"SOLUTION")
			continue
        elseif occursin(r"^[a-z,\s]+$",split(l)[1])
            println(string("ERROR: Unrecognized header: ",split(l)[1],"; currently looks for TITLE, VARIABLES, and ZONE; might need to add it"))
        else
			header = false
		end
        i = i + 1
	end

	return tm
end

function tec_to_df(f,variables=0)
	# Tecplot to dataframe
    fh = open(f)
    # Read in header
    i = 1
    header = true
	tm = ""
    while header
        l = readline(fh)
        if startswith(l,"TITLE")
			title = strip(split(l,"=")[2])
        elseif startswith(l,"VARIABLES")
			variables = split(split(l,"=")[2],"\" \"")
			variables = [replace(v,"\""=>"") for v in variables]
			variables = [strip(replace(v,r"\(.+\)"=>"")) for v in variables]
        elseif startswith(l,"ZONE")
            tm = split(split(l,"\"")[2])
			if length(tm) <= 2
				i = i + 1
				continue
			else # Probably material file
            	tm = tm[3]
			end
            if occursin("days\"",tm)
				tm = replace(tm,"days\""=>"")
			end
        elseif startswith(l,"FILETYPE")
        	i = i + 1
			continue
        elseif startswith(l,"SOLUTION")
        	i = i + 1
			continue
		elseif isreal(parse(Float64,split(l)[1]))
			header = false
        else
            println(string("ERROR: Unrecognized header: ",split(l)[1],"; currently looks for TITLE, VARIABLES, and ZONE; might need to add it"))
		end
        i = i + 1
	end

	seek(fh,0)
	for i in 1:i-2
		readline(fh)
	end
	lns = readlines(fh)

	if occursin("gdkm",f)
		n_lines = n_gdkm
	else
		n_lines = numnodes
	end
	dtemp = Array{Float64,2}(undef,n_lines,length(variables))
	for i in 1:n_lines
		vs = split(strip(lns[i]))
		for j in 1:length(vs)
			dtemp[i,j] = parse(Float64,vs[j])
		end
	end

	d = Array{Float64,2}(undef,numnodes,length(variables))
	d[:,:] .= 0.
	if occursin("gdkm",f)
		d[gdkm_map[:,1],:] .= dtemp[:,:]
	else
		d[:,:] .= dtemp[:,:]
	end

	#d = readdlm(fh)
	close(fh)
	df = convert(DataFrame,d)
	var_sym = Array{Symbol,length(variables)}
	variables = [replace(v,"/"=>"_") for v in variables]
	variables = [replace(v,"-"=>"_") for v in variables]
	variables = [replace(v," "=>"_") for v in variables]
	var_sym = [Symbol(v) for v in variables]
	rename!(df,var_sym)

	if (in(:X_coordinate,names(df)))
		delete!(df,:X_coordinate)
	end
	if (in(:Y_coordinate,names(df)))
		delete!(df,:Y_coordinate)
	end
	if (in(:Z_coordinate,names(df)))
		delete!(df,:Z_coordinate)
	end

	if (in(:X,names(df)))
		delete!(df,:X)
	end
	if (in(:Y,names(df)))
		delete!(df,:Y)
	end
	if (in(:Z,names(df)))
		delete!(df,:Z)
	end
	if (in(:Zone,names(df)))
		delete!(df,:Zone)
	end

	return df
end

function write_vtk(f,coords, elems; diff_df=DataFrame(), variables=0, fi=-1)
	println(f)
	fnlst = split(f,'.')
    root = fnlst[1]
    if occursin("sca_node",fnlst[2])
		lst = split(fnlst[2],"_")
		if (fi==-1)
			fi = lst[1]
		end
        root = string(root,"_sca_node")
    elseif occursin(fnlst[2],"sca_gdkm_node")
		lst = split(fnlst[2],"_")
		if (fi==-1)
			fi = lst[1]
		end
        root = string(root,"_sca_gdkm_node")
    elseif occursin("con_node",fnlst[2])
		lst = split(fnlst[2],"_")
		if (fi==-1)
			fi = lst[1]
		end
        root = string(root,"_con_node")
    elseif occursin("con_gdkm_node",fnlst[2])
		lst = split(fnlst[2],"_")
		if (fi==-1)
			fi = lst[1]
		end
        root = string(root,"_con_gdkm_node")
    elseif occursin("mat_node",fnlst[2])
        root = string(root,"_mat_node")
	end

	df = tec_to_df(f,variables)

	if !isdir(parsed_args["output_dir"])
		mkdir(parsed_args["output_dir"])
	end
	if fi > -1
		fname = string(parsed_args["output_dir"],"/",root,lpad(fi,4,"0"))
	else
		fname = string(parsed_args["output_dir"],"/",root)
	end
	vtkfile = vtk_grid(fname, coords, elems)
	for v in names(df)
		vtk_point_data(vtkfile, df[v], string(v))
		if in(v,names(diff_df))
			if length(var_diff)>0
				if in(string(v),var_diff)
					vtk_point_data(vtkfile, df[v]-diff_df[v], string("Delta ",v))
				end
			else
				vtk_point_data(vtkfile, df[v]-diff_df[v], string("Delta ",v))
			end
		end
	end

	#outfiles = vtk_save(vtkfile)

	return vtkfile

end

s = ArgParseSettings(description = "Convert FEHM gdkm tecplot output to vtu2 files")

@add_arg_table! s begin
    "--sca_diff"
		help = "Sca file to subtract to create delta variables"
		arg_type = String
		default = "none"
    "--sca_gdkm_diff"
		help = "GDKM sca file to subtract to create delta variables"
		arg_type = String
		default = "none"
    "--var_diff"
		help = "Variables to create delta variables for"
		arg_type = String
		default = "none"
    "--output_dir", "-o"         # another option, with short form
		help = "Folder to create and place vtu2 files"
		arg_type = String
		default = "vtu"
    "--mesh_file", "-m"         # another option, with short form
		help = "Mesh file location"
		arg_type = String
		default = "none"
    "--run_name", "-r"         # another option, with short form
		help = "Run name: prefix for visualization file names"
		arg_type = String
		default = "*"
end

parsed_args = parse_args(s) # the result is a Dict{String,Any}
println("Parsed args:")
for (key,val) in parsed_args
    println("  $key  =>  $(repr(val))")
end

if parsed_args["var_diff"] != "none"
	var_diff = split(parsed_args["var_diff"],",")
else
	var_diff = Set()
end

if parsed_args["mesh_file"] == "none"
	fh = open("fehmn.files")
	mesh = ""
	for l in eachline(fh)
		if (startswith(l,"grid"))
			print(l)
		mesh = strip(split(l,":")[2])
			break
		end
	end
	close(fh)
else
	mesh = parsed_args["mesh_file"]
end

fh = open(mesh)
lns = readlines(fh)
close(fh)
numnodes = parse(Int,strip(lns[2]))

# Points
coords = zeros(Float64,(3,numnodes))
for i in 3:numnodes+2
    vs = split(strip(lns[i]))
    coords[:,i-2] = [parse(Float64,vs[2]), parse(Float64,vs[3]), parse(Float64,vs[4])]
end

# Elements
vs = split(strip(lns[numnodes+5]))

numconns = parse(Int,vs[1])
numelems = parse(Int,vs[2])
elems = Array{MeshCell}(undef,numelems)
# Determine cell type
if ( numconns == 3 )
	celltype = VTKCellTypes.VTK_TRIANGLE
elseif ( numconns == 4 )
	if length(unique(coords[1,:]))>1 && length(unique(coords[2,:]))>1 && length(unique(coords[3,:]))>1
		celltype = VTKCellTypes.VTK_TETRA
	else
		celltype = VTKCellTypes.VTK_QUAD
	end
elseif ( numconns == 8 )
	celltype = VTKCellTypes.VTK_HEXAHEDRON
end
j = 1
for i in numnodes+6:numnodes+5+numelems
    vs = split(strip(lns[i]))
	if ( numconns == 3 )
    	elems[j] = MeshCell(celltype, [parse(Int,vs[2]), parse(Int,vs[3]), parse(Int,vs[4])])
	elseif ( numconns == 4 )
		# Check for triangle elements for DFNs
		if parse(Int,vs[5]) == 0
    			elems[j] = MeshCell(VTKCellTypes.VTK_TRIANGLE, [parse(Int,vs[2]), parse(Int,vs[3]), parse(Int,vs[4])])
		else
    			elems[j] = MeshCell(celltype, [parse(Int,vs[2]), parse(Int,vs[3]), parse(Int,vs[4]), parse(Int,vs[5])])
		end
	elseif ( numconns == 8 )
    	elems[j] = MeshCell(celltype, [parse(Int,vs[2]), parse(Int,vs[3]), parse(Int,vs[4]), parse(Int,vs[5]), parse(Int,vs[6]), parse(Int,vs[7]), parse(Int,vs[8]), parse(Int,vs[9])])
	end
	global j = j + 1
end

# Scalar files
fs = Glob.glob(string(parsed_args["run_name"],".","*_sca_node.dat"))
if (length(fs) > 0)
	pvd_sca = paraview_collection(string(parsed_args["output_dir"],"/sca"))
	# Collect variables from first file
	sca_o = tec_to_df(fs[1])
	#print(names(sca_o))
	variables = [string(nm) for nm in names(sca_o)]
	if occursin("none",parsed_args["sca_diff"])
		sca_zero = DataFrame()
	else
		println(parsed_args["sca_diff"])
		sca_zero = tec_to_df(parsed_args["sca_diff"],variables)
	end
	#fi = Array(Int,len(fs))
	for (fi,f) in enumerate(fs)
	#	fi[fi] = fi
		vfile = write_vtk(f,coords,elems,diff_df=sca_zero,variables=variables,fi=fi);
		tm = get_tec_time(f)
		collection_add_timestep(pvd_sca, vfile, parse(Float64,tm))
	end
	#pmap(write_vtk,fs,fi)
	vtk_save(pvd_sca)
end

# Concentration files
fs = Glob.glob(string(parsed_args["run_name"],".","*_con_node.dat"))
if (length(fs) > 0)
	pvd_con = paraview_collection(string(parsed_args["output_dir"],"/con"))
	# Collect variables from first file
	con_zero = tec_to_df(fs[1])
	variables = [string(nm) for nm in names(con_zero)]
	con_zero = DataFrame()
	for (fi,f) in enumerate(fs)
		vfile = write_vtk(f,coords,elems,diff_df=con_zero,variables=variables,fi=fi)
		tm = get_tec_time(f)
		collection_add_timestep(pvd_con, vfile, parse(Float64,tm))
	end
	vtk_save(pvd_con)
end

fs = Glob.glob(string(parsed_args["run_name"],".","mat_node.dat"))
if (length(fs) > 0)
	for (fi,f) in enumerate(fs)
		vfile = write_vtk(f,coords,elems,fi=fi)
		vtk_save(vfile)
	end
end

# Write gdkm.vtu
if isfile("gdkm.dat")
	fh = open("gdkm.dat")
	readline(fh)
	n_gdkm = parse(Int,split(strip(readline(fh)))[2])
	d = Array{Float64,(n_gdkm,2)}
	for i in 1:n_gdkm
		vs = split(strip(readline(fh)))
		d[i,:] = [parse(Float64,vs[2]),parse(Float64,vs[3])]
	end
	readline(fh)
	gdkm_map = Array{Float64,(n_gdkm,2)}
	for i in 1:n_gdkm
		vs = split(strip(readline(fh)))
		gdkm_map[i,:] = [parse(Int,vs[1]),numnodes+parse(Int,vs[4])]
	end
	close(fh)
	gdkm_dict = Dict(zip(gdkm_map[:,2],gdkm_map[:,1]))
	gdkm_dat = ones(Float64,(numnodes,2))*-9999
	gdkm_dat[gdkm_map[:,1],:] = d[:,:]

	if !isdir(parsed_args["output_dir"])
		mkdir(parsed_args["output_dir"])
	end
	vtkfile = vtk_grid(string(parsed_args["output_dir"],"/gdkm"), coords, elems)
	vtk_point_data(vtkfile, gdkm_dat[:,1], "VFrac Primary")
	vtk_point_data(vtkfile, gdkm_dat[:,2], "GDKM X")
	outfiles = vtk_save(vtkfile)
end

# Write perm.vtu
if isfile("perm.dat")
	perm_dat = ones(Float64,(numnodes,3))*-9999
	fh = open("perm.dat")
	readline(fh)
	for l in eachline(fh)
		vs = split(strip(l))
		if isempty(strip(l))
			break
		end
		if parse(Int,vs[1]) > 1
			perm_dat[gdkm_dict[parse(Int,vs[1])],:] = [parse(Float64,vs[4]),parse(Float64,vs[5]),parse(Float64,vs[6])]
		end
	end
	close(fh)

	if !isdir(parsed_args["output_dir"])
		mkdir(parsed_args["output_dir"])
	end
	vtkfile = vtk_grid(string(parsed_args["output_dir"],"/perm"), coords, elems)
	vtk_point_data(vtkfile, perm_dat[:,1], "Fracture Perm X")
	vtk_point_data(vtkfile, perm_dat[:,1], "Fracture Perm Y")
	vtk_point_data(vtkfile, perm_dat[:,1], "Fracture Perm Z")
	outfiles = vtk_save(vtkfile)
end

# Scalar gdkm files
fs = Glob.glob(string(parsed_args["run_name"],".","*_sca_gdkm_node.dat"))
if (length(fs) > 0)
	pvd_sca = paraview_collection(string(parsed_args["output_dir"],"/sca_gdkm"))
	# Collect variables from first file
	sca_o = tec_to_df(fs[1])
	#print(names(sca_o))
	variables = [string(nm) for nm in names(sca_o)]
	if occursin("none",parsed_args["sca_gdkm_diff"])
		sca_zero = DataFrame()
	else
		println(parsed_args["sca_gdkm_diff"])
		sca_zero = tec_to_df(parsed_args["sca_gdkm_diff"])
	end
	#fi = Array(Int,len(fs))
	for (fi,f) in enumerate(fs)
	#	fi[fi] = fi
		vfile = write_vtk(f,coords,elems,diff_df=sca_zero,variables=variables,fi=fi);
		tm = get_tec_time(f)
		collection_add_timestep(pvd_sca, vfile, parse(Float64,tm))
	end
	#pmap(write_vtk,fs,fi)
	vtk_save(pvd_sca)
end

# Concentration files
fs = Glob.glob(string(parsed_args["run_name"],".","*_con_gdkm_node.dat"))
if (length(fs) > 0)
	pvd_con = paraview_collection(string(parsed_args["output_dir"],"/con_gdkm"))
	# Collect variables from first file
	con_zero = tec_to_df(fs[1])
	variables = [string(nm) for nm in names(con_zero)]
	con_zero = DataFrame()
	for (fi,f) in enumerate(fs)
		vfile = write_vtk(f,coords,elems,diff_df=con_zero,variables=variables,fi=fi)
		tm = get_tec_time(f)
		collection_add_timestep(pvd_con, vfile, parse(Float64,tm))
	end
	vtk_save(pvd_con)
end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
"""
_compute_mutliplier: Additional function to estimate Lagrange multiplier of the problems
    (guarantee if LICQ holds)
"""
function _compute_mutliplier(pb    :: AbstractNLPModel,
                             x     :: Iterate,
                             gx    :: Iterate,
                             cx    :: Iterate,
                             Jx    :: AbstractMatrix;
                             active_prec_c :: Float64 = 1e-6,
                			 active_prec_b :: Float64 = 1e-6)

 n  = length(x)
 nc = cx == nothing ? 0 : length(cx)

 #active res_bounds
 Ib = findall(x->(norm(x) <= active_prec_b),
			      min(abs.(x - pb.meta.lvar),
				      abs.(x - pb.meta.uvar)))
 if nc != 0
  #active constraints
  Ic = findall(x->(norm(x) <= active_prec_c),
                   min(abs.(cx-pb.meta.ucon),
                   abs.(cx-pb.meta.lcon)))

  Jc = hcat(Matrix(1.0I, n, n)[:,Ib], Jx'[:,Ic])
 else
  Ic = []
  Jc = hcat(Matrix(1.0I, n, n)[:,Ib])
 end

 mu, lambda = zeros(n), zeros(nc)
 if (Ib != []) || (Ic  != [])
  l = Jc \ (- gx)
  mu[Ib], lambda[Ic] = l[1:length(Ib)], l[length(Ib)+1:length(l)]
 end

 return mu, lambda
end

function _compute_mutliplier(pb    :: AbstractMPCCModel,
                             x     :: AbstractVector,
                             gx    :: AbstractVector,
                             cx    :: AbstractVector,
                             Jx    :: AbstractMatrix,
                             Gx    :: AbstractVector,
                             JGx   :: AbstractMatrix,
                             Hx    :: AbstractVector,
                             JHx   :: AbstractMatrix;
                             active_prec_c  :: Float64 = 1e-6,
                             active_prec_cc :: Float64 = 1e-3,
                			 active_prec_b  :: Float64 = 1e-6)

 n, ncc  = length(x), pb.meta.ncc
 nc = cx == nothing ? 0 : length(cx)

 #active res_bounds
 Ib = findall(x->(norm(x) <= active_prec_b),
			      min(abs.(x - pb.meta.lvar),
				      abs.(x - pb.meta.uvar)))
 if nc != 0
  #active constraints
  Ic = findall(x->(norm(x) <= active_prec_c),
                   min(abs.(cx-pb.meta.ucon),
                   abs.(cx-pb.meta.lcon)))

  Jc = hcat(Matrix(1.0I, n, n)[:,Ib], Jx'[:,Ic])
 else
  Ic = []
  Jc = hcat(Matrix(1.0I, n, n)[:,Ib])
 end

 if ncc != 0
  #active constraints
  IG = findall(x->(norm(x) <= active_prec_cc),abs.(Gx))
  IH = findall(x->(norm(x) <= active_prec_cc),abs.(Hx))
  Jc = (length(Ic)+length(Ib)) != 0 ? hcat(Jc,-JGx'[:,IG],-JHx'[:,IH]) : hcat(-JGx'[:,IG],-JHx'[:,IH])
 else
  IG, IH = [], []
 end

 mu, lambda, etaG, etaH = zeros(n), zeros(nc), zeros(ncc), zeros(ncc)
 if (Ib != []) || (Ic  != []) || (IG != []) || (IH  != [])
  l = Jc \ (- gx)
  mu[Ib], lambda[Ic], etaG[IG], etaH[IH] = l[1:length(Ib)], l[length(Ib)+1:length(Ib)+length(Ic)],
                                                            l[length(Ib)+length(Ic)+1:length(Ib)+length(Ic)+length(IG)],
                                                            l[length(Ib)+length(Ic)+length(IG)+1:length(Ib)+length(Ic)+length(IG)+length(IH)]
 end

 return mu, lambda, etaG, etaH
end


text:
using ChainMap
import DataStructures
import NullableArrays
import DataFrames
using Base.Test
na =  NullableArrays.NullableArray
a = na(["one", "two"], [false, true])
result = @chain begin
    Dict(:b => na([1, 2]), :c => na(["I", "II"]))
    @chain_map begin
        :b
        sum
        get
        string
        *(~a, " ", _, " ", ~:c)
    end
end

@test get(result[1]) == "one 3 I"
@test result.isnull == [false, true]
@test vcat(1) == @chain 1 vcat
@test vcat(2, 1) == @chain 1 vcat(2, _)
@test ( @chain 1 vcat(_, 2) vcat(_, 3) ) ==
    @chain ( @chain 1 vcat(_, 2) ) vcat(_, 3)
@test 1 == @chain 1

chain_block = @chain begin
    1
    vcat(_, 2)
end

@test chain_block == @chain 1 vcat(_, 2)
lambda_function = @lambda vcat(_, 2)
@test lambda_function(1) == vcat(1, 2)
f = :map
e = :(vcat(_, 1))
lambda(f, e)

_ = [1, 2]
@test map(_ -> vcat(_, 1), _) == @lambda map vcat(_, 1)
e = :(_ + 1)
f = :(NullableArrays.map(lift = true))

lambda(f, e)

_ = NullableArrays.NullableArray([1, 2], [false, true])

result = @lambda NullableArrays.broadcast(lift = true) _ + 1

@test result.values[1] == 2
@test result.isnull == [false, true]

# `f` must be a call
@test_throws ErrorException lambda(:(import ChainMap), :(_ + 1) )
_ = [1, 2]
@test map(_ -> vcat(_, 1), _) == @map vcat(_, 1)
merge_test = @chain begin
    collect_arguments(1, a = 2, b = 3)
    merge(_, collect_arguments(4, a = 5, c = 6) )
end

@test merge_test == collect_arguments(1, 4, a = 5, b = 3, c = 6)
merge_test = @chain begin
    collect_arguments([1, 2])
    unshift(_, vcat)
    LazyCall(_, map)
    merge(_, collect_arguments([3, 4]) )
    run
end

@test merge_test == [[1, 3], [2, 4]]
push_test = @chain begin
    collect_arguments(1, a = 2, b = 3)
    push(_, 4, a = 5, c = 6)
end

@test push_test == collect_arguments(1, 4, a = 5, b = 3, c = 6)
push_test = @chain begin
    collect_arguments([1, 2])
    unshift(_, vcat)
    LazyCall(_, map)
    push(_, [3, 4])
    run
end

@test push_test == [[1, 3], [2, 4]]
unshift_test = @chain begin
    collect_arguments(2, a = 3)
    unshift(_, 1)
end

@test unshift_test == collect_arguments(1, 2, a = 3)
unshift_test = @chain begin
    collect_arguments([1, 2], [3, 4])
    LazyCall(_, map)
    unshift(_, vcat)
    run
end

@test unshift_test == [[1, 3], [2, 4]]
a = collect_arguments(1, 2, a = 3, b = 4)
@test a.positional == (1, 2)
@test a.keyword == Dict{Symbol, Any}(:a => 3, :b => 4)
l = collect_call(vcat, [1, 2], [3, 4])
@test l.function_call == vcat
@test l.arguments == collect_arguments([1, 2], [3, 4])
run_test = @chain begin
    collect_arguments([1, 2], [3, 4])
    unshift(_, vcat)
    collect_arguments(_, map)
    run
end

@test run_test == map(vcat, [1, 2], [3, 4])
run_test = @chain begin
    collect_arguments([1, 2], [3, 4])
    unshift(_, vcat)
    LazyCall(_, map)
    run
end

@test run_test == map(vcat, [1, 2], [3, 4])
run_test = @chain begin
    collect_arguments([1, 2], [3, 4])
    unshift(_, vcat)
    run(_, map)
end

@test run_test == map(vcat, [1, 2], [3, 4])
run_test = @chain begin
    collect_arguments([1, 2], [3,4])
    LazyCall(_, vcat)
    run(_, map)
end

@test run_test == map(vcat, [1, 2], [3, 4])
test_function(arguments...; keyword_arguments...) =
    (arguments, keyword_arguments)

@test ( @lazy_call test_function(1, 2, a = 3) ) ==
    collect_call(test_function, 1, 2, a = 3)
binaryfun(a, b, c) = Expr(:call, b, a, c)
chainback(a, b, c) = Expr(:call, c, b, a)

@nonstandard binaryfun chainback

@test vcat(1, 2) == @binaryfun 1 vcat 2
@test vcat(3, 2) == @chainback 2 3 vcat

@test "See documentation of [`binaryfun`](@ref)" ==
    @chain (@doc @binaryfun) string chomp

a = 1
_ = DataFrames.DataFrame(a = 2)

@test DataFrames.DataFrame(b = _[:a] + a, c = :d) ==
   @with DataFrames.DataFrame(b = :a + a, c = ^(:d))
e = Expr(:parameters, Expr(:..., :a) )
first = :parameters
second = :...

@test ChainMap.double_match(e, first, second)
@test !ChainMap.double_match(:b, :parameters, :...)
e = Expr(:parameters, Expr(:..., :a) )
@test ChainMap.replace_key(e, :b) == Expr(:parameters, Expr(:..., :b) )

e = Expr(:..., :a)
@test ChainMap.replace_key(e, :b) == Expr(:..., :b)

@test ChainMap.replace_key(:a, :b) == :b
a = Expr(:..., :a)

@test ChainMap.unparameterize(Expr(:parameters, a) ) == a

@test ChainMap.unparameterize(:b) == :b
e = :a
symbol = :z
d = Dict(:a => :b, :(a + 1) => :(b + 1))

ChainMap.add_key!(d, e, symbol)
@test d[:a] == :b

ChainMap.add_key!(d, :c, symbol)
@test d[:c] == :z

e = Expr(:parameters, Expr(:..., :d) )
@test ChainMap.add_key!(d, e, symbol) == Expr(:..., :z)
@test d[e] == Expr(:parameters, Expr(:..., :z) )
d = Dict()
e = :( 1 + ~(a))
ChainMap.replace_record!(e, d)
@test :a in keys(d)
f = x -> x == :a
@test ChainMap.negate(f)(:b)
e = :(~_ + 1)
ChainMap.split_anonymous(e)
@test 1 == @unweave 1
@test 1 == @unweave +(1)

A = [1, 2]
B = ( [5, 6], [7, 8] )

unweave_test = @chain begin
    @unweave vcat(~A, ~[3, 4], ~(B...) )
    run(_, map)
end

@test unweave_test ==
      map((a, c, b...) -> vcat(a, c, b...), A, [3, 4], B...)

keyword_test(; keyword_arguments...) = keyword_arguments

a = keyword_test(a = 1, b = 2)

unweave_keyword_test = @chain begin
    @unweave keyword_test(c = 3; ~(a...))
    run
end

@test unweave_keyword_test == keyword_test(c = 3; a... )

# Can splat no more than one positional argument
@test_throws ErrorException unweave(:( ~(a...) + ~(b...) ))

# Can splat no more than one keyword argument
@test_throws ErrorException unweave(:( ~(;a...) + ~(;b...) ))
@test 1 == @unweave 1 1
@test 1 == @unweave 1 +(1)

e = :(vcat(~a, ~b) )
f = :broadcast
unweave(f, e)

a = [1, 2]
b = [3, 4]

@test broadcast((a, b) -> vcat(a, b), a, b) ==
    @unweave broadcast vcat(~a, ~b)
@test 1 == @unweave +(1) 1
@test 1 == @unweave +(1) +(1)

e = :(~a + ~b)
f = :(NullableArrays.broadcast(lift = true))

unweave(f, e)

a = NullableArrays.NullableArray([1, 2])
b = NullableArrays.NullableArray([3, 4], [false, true])

result = @unweave broadcast(lift = true) ~a + ~b

@test result.values[1] == 4
@test result.isnull == [false, true]

# `f` must be a call
@test_throws ErrorException unweave(:(import ChainMap), :(~_ + 1) )
a = [1, 2]
b = [3, 4]

@test broadcast((a, b) -> vcat(a, b), a, b) ==
    @broadcast vcat(~a, ~b)
@test bitnot(1) == ~1
na = NullableArrays.NullableArray
a = na(["one", "two"], [false, true])
result = @chain begin
    Dict(:b => na([1, 2]), :c => na(["I", "II"]))
    @chain_map begin
        :b
        sum
        get
        string
        *(~a, " ", _, " ", ~:c)
    end
end

@test get(result[1]) == "one 3 I"
@test result.isnull == [false, true]
along() = "dummy function; could be a fancy view some day"

Base.run(A::AbstractArray,
         map_call::typeof(map), map_function::Function,
         along_call::LazyCall{typeof(along)},
         reduce_call::typeof(reduce), reduce_function::Function) =
    mapreducedim(map_function, reduce_function, A,
                 along_call.arguments.positional[1] )
fancy = @chain begin
    [1, 2, 3, 4]
    reshape(_, 2, 2)
    collect_arguments(
        _,
        map,
        _ -> -(_, 1),
        @lazy_call( along(1) ),
        reduce,
        +)
    run(_)
end

boring = mapreducedim(x -> x - 1, +, reshape([1, 2, 3, 4], 2, 2), 1)

@test fancy == boring


text:
struct TrackedReal{T<:Real} <: Real
  tracker::Tracked{T}
end

TrackedReal(x::Real) = TrackedReal(Tracked(Call(nothing), x, zero(x)))

tracker(x::TrackedReal) = x.tracker

track(f::Call, x::Real) = TrackedReal(Tracked(f, x, zero(x)))

back!(x::TrackedReal) = back!(x, 1)

function Base.show(io::IO, x::TrackedReal)
  show(io, data(x))
  print(io, " (tracked)")
end

Base.decompose(x::TrackedReal) = Base.decompose(data(x))

Base.convert(::Type{TrackedReal{T}}, x::TrackedReal{T}) where T = x

# This cuts derivatives, fix if needed.
# Base.convert(::Type{TrackedReal{T}}, x::TrackedReal) where T =
#   TrackedReal(Tracked(x.tracker.f, convert(T, x.tracker.data)))

Base.convert(::Type{TrackedReal{T}}, x::Real) where T = TrackedReal(convert(T, x))

Base.:(<)(x::TrackedReal, y::TrackedReal) = data(x) < data(y)
Base.:(==)(x::TrackedReal, y::TrackedReal) = data(x) == data(y)

for f in :[isinf, isnan, isfinite].args
  @eval Base.$f(x::TrackedReal) = Base.$f(data(x))
end

Base.Printf.fix_dec(x::TrackedReal, n::Int) = Base.Printf.fix_dec(data(x), n)

Base.promote_rule(::Type{TrackedReal{S}},::Type{T}) where {S,T} =
  TrackedReal{promote_type(S,T)}

using DiffRules, SpecialFunctions, NaNMath

for (M, f, arity) in DiffRules.diffrules()
  arity == 1 || continue
  @eval begin
    $M.$f(a::TrackedReal) = track($M.$f, a)
    back(::typeof($M.$f), Δ::Real, a::TrackedReal) =
      back(a, Δ * $(DiffRules.diffrule(M, f, :(data(a)))))
  end
end

for (M, f, arity) in DiffRules.diffrules()
  arity == 2 || continue
  da, db = DiffRules.diffrule(M, f, :(data(a)), :(data(b)))
  @eval begin
    $M.$f(a::TrackedReal, b::TrackedReal)  = track($M.$f, a, b)
    $M.$f(a::TrackedReal, b::Real) = track($M.$f, a, b)
    $M.$f(a::Real, b::TrackedReal) = track($M.$f, a, b)
    function back(::typeof($M.$f), Δ::Real, a::Real, b::Real)
      @back(a, Δ * $da)
      @back(b, Δ * $db)
    end
  end
end

# Eliminating ambiguity
import Base:^

^(a::TrackedReal, b::Integer) = track(^, a, b)

# Tuples

struct TrackedTuple{T<:Tuple}
  tracker::Tracked{T}
end

tracker(xs::TrackedTuple) = xs.tracker

accum!(x::Tuple, Δ::Tuple) = accum!.(x, Δ)
init_grad(x::Tuple) = init_grad.(x)
zero_grad!(x::Tuple) = zero_grad!.(x)

track(f::Call, xs::Tuple) = TrackedTuple(Tracked(f, xs))

function Base.show(io::IO, xs::TrackedTuple)
  show(io, data(xs))
  print(io, " (tracked)")
end

Base.length(x::TrackedTuple) = length(data(x))

Base.getindex(xs::TrackedTuple, i::Integer) = track(getindex, xs, i)

back(::typeof(getindex), Δ, t, i) =
  back(t, ntuple(j -> i == j ? Δ : 0, length(t)))

# Array collection

function collect(xs)
  xs = Base.collect(xs)
  track(Call(collect, xs), data.(xs))
end

function scan(c::Call{typeof(collect)})
  foreach(scan, c.args[1])
end

function back(::typeof(collect), Δ, xs)
  foreach((x, Δ) -> @back(x, Δ), xs, Δ)
end