text:
"""
    shift_constant(set::MOI.AbstractScalarSet, offset)

Returns a new scalar set `new_set` such that `func`-in-`set` is equivalent to
`func + offset`-in-`new_set`.

Only define this function if it makes sense to!

Use [`supports_shift_constant`](@ref) to check if the set supports shifting:
```Julia
if supports_shift_constant(typeof(old_set))
    new_set = shift_constant(old_set, offset)
    f.constant = 0
    add_constraint(model, f, new_set)
else
    add_constraint(model, f, old_set)
end
```

See also [`supports_shift_constant`](@ref).

## Examples

The call `shift_constant(MOI.Interval(-2, 3), 1)` is equal to
`MOI.Interval(-1, 4)`.
"""
function shift_constant end

"""
    supports_shift_constant(::Type{S}) where {S<:MOI.AbstractSet}

Return `true` if [`shift_constant`](@ref) is defined for set `S`.

See also [`shift_constant`](@ref).
"""
supports_shift_constant(::Type{S}) where {S<:MOI.AbstractSet} = false

function shift_constant(
    set::Union{MOI.LessThan{T},MOI.GreaterThan{T},MOI.EqualTo{T}},
    offset::T,
) where {T}
    return typeof(set)(MOI.constant(set) + offset)
end
supports_shift_constant(::Type{<:MOI.LessThan}) = true
supports_shift_constant(::Type{<:MOI.GreaterThan}) = true
supports_shift_constant(::Type{<:MOI.EqualTo}) = true

function shift_constant(set::MOI.Interval, offset)
    return MOI.Interval(set.lower + offset, set.upper + offset)
end
supports_shift_constant(::Type{<:MOI.Interval}) = true

const ScalarLinearSet{T} =
    Union{MOI.EqualTo{T},MOI.LessThan{T},MOI.GreaterThan{T}}
const VectorLinearSet = Union{MOI.Zeros,MOI.Nonnegatives,MOI.Nonpositives}

vector_set_type(::Type{<:MOI.EqualTo}) = MOI.Zeros
vector_set_type(::Type{<:MOI.LessThan}) = MOI.Nonpositives
vector_set_type(::Type{<:MOI.GreaterThan}) = MOI.Nonnegatives

scalar_set_type(::Type{<:MOI.Zeros}, T::Type) = MOI.EqualTo{T}
scalar_set_type(::Type{<:MOI.Nonpositives}, T::Type) = MOI.LessThan{T}
scalar_set_type(::Type{<:MOI.Nonnegatives}, T::Type) = MOI.GreaterThan{T}

"""
    is_diagonal_vectorized_index(index::Base.Integer)

Return whether `index` is the index of a diagonal element in a
[`MOI.AbstractSymmetricMatrixSetTriangle`](@ref) set.
"""
function is_diagonal_vectorized_index(index::Base.Integer)
    # See https://en.wikipedia.org/wiki/Triangular_number#Triangular_roots_and_tests_for_triangular_numbers
    perfect_square = 1 + 8index
    return isqrt(perfect_square)^2 == perfect_square
end

"""
    side_dimension_for_vectorized_dimension(n::Integer)

Return the dimension `d` such that
`MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(d))` is `n`.
"""
function side_dimension_for_vectorized_dimension(n::Base.Integer)
    return div(isqrt(1 + 8n), 2)
end


text:
using JuMP
using Clp
using DiffOpt
using Test

# This script creates the JuMP problem for a small unit commitment instance

ATOL=1e-4
RTOL=1e-4

## Problem data
unit_codes = [1, 2] # Generator identifiers
load_names = ["Load1", "Load2"] # Load identifiers
n_periods = 4 # Number of time periods
Pmin = Dict(1 => fill(0.5, n_periods), 2 => fill(0.5, n_periods)) # Minimum power output (pu)
Pmax = Dict(1 => fill(3.0, n_periods), 2 => fill(3.0, n_periods)) # Maximum power output (pu)
RR = Dict(1 => 0.25, 2 => 0.25) # Ramp rates (pu/min)
P0 = Dict(1 => 0.0, 2 => 0.0) # Initial power output (pu)
D = Dict("Load1" => [1.0, 1.2, 1.4, 1.6], "Load2" => [1.0, 1.2, 1.4, 1.6]) # Demand
Cp = Dict(1 => 1000.0, 2 => 1500.0) # Generation cost coefficient ($/pu)
Cnl = Dict(1 => 500.0, 2 => 1000.0) # No-load cost ($)

model = Model(() -> diff_optimizer(Clp.Optimizer))

## Variables
@variable(model, 0 <= u[g in unit_codes, t in 1:n_periods] <= 1) # Commitment
@variable(model, p[g in unit_codes, t in 1:n_periods] >= 0) # Power output

## Constraints

# Energy balance
@constraint(
    model,
    energy_balance_cons[t in 1:n_periods],
    sum(p[g, t] for g in unit_codes) == sum(D[l][t] for l in load_names)
)

# Generation limits
@constraint(model, [g in unit_codes, t in 1:n_periods], Pmin[g][t] * u[g, t] <= p[g, t])
@constraint(model, [g in unit_codes, t in 1:n_periods], p[g, t] <= Pmax[g][t] * u[g, t])

# Ramp rates
@constraint(model, [g in unit_codes, t in 2:n_periods], p[g, t] - p[g, t - 1] <= 60 * RR[g])
@constraint(model, [g in unit_codes], p[g, 1] - P0[g] <= 60 * RR[g])
@constraint(model, [g in unit_codes, t in 2:n_periods], p[g, t - 1] - p[g, t] <= 60 * RR[g])
@constraint(model, [g in unit_codes], P0[g] - p[g, 1] <= 60 * RR[g])

# Objective
@objective(
    model,
    Min,
    sum((Cp[g] * p[g, t]) + (Cnl[g] * u[g, t]) for g in unit_codes, t in 1:n_periods),
)

optimize!(model)

diff_opt = backend(model).optimizer.model
v = MOI.get(model, MOI.ListOfVariableIndices())
nv = length(v)

MOI.set.(diff_opt, DiffOpt.BackwardInVariablePrimal(), v, ones(nv))

DiffOpt.backward(diff_opt)

# sensitivity wrt linear objective
grad = MOI.get(diff_opt, DiffOpt.BackwardOutObjective())
@test MOI.Utilities.isapprox_zero(grad, ATOL)

# sensitivity wrt RHS of constraints
for (idx, econs) in enumerate(energy_balance_cons)
    grad_energy_balance = JuMP.constant(MOI.get(model, DiffOpt.BackwardOutConstraint(), econs))
    @test !≈(grad_energy_balance, 0)
end

ext: jl


text:
# chord diagrams

using LinearAlgebra
using SparseArrays
using GraphRecipes
using Plots

adjmat = Symmetric(sparse(rand(0:1,8,8)))

plot(
    graphplot(adjmat,
              method=:chorddiagram,
              names=[text(string(i), 8) for i in 1:8],
              linecolor=:black,
              fillcolor=:lightgray),
    )

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
function build_linear_formula(factors::Array{Symbol, 1})
    linear_formula = Expr(:call)
    linear_formula.args = vcat(:+, factors)
    return Formula(:y, linear_formula)
end

function get_model_variables(formula::Formula)
    variables = Array{Any, 1}()
    push!(variables, Terms(formula).terms...)
    return variables
end

"""
    scale_orthogonal!(design::Array{Float64, 2}, factors::Array{T, 1}) where T <: Any

Orthogonally scale and center factors of a design using design and factor
limits.

# Examples

```jldoctest
julia> using ExperimentalDesign, DataStructures, StatsModels

julia> A = [5. 2. -1.; 5. 3. 0.; -5. 1. -2.]
3×3 Array{Float64,2}:
  5.0  2.0  -1.0
  5.0  3.0   0.0
 -5.0  1.0  -2.0

julia> factors = OrderedDict([(:f1, [-5., 0., 5.]), (:f2, [1., 2., 3.]), (:f3, [-2., -1., 0.])])
DataStructures.OrderedDict{Symbol,Array{Float64,1}} with 3 entries:
  :f1 => [-5.0, 0.0, 5.0]
  :f2 => [1.0, 2.0, 3.0]
  :f3 => [-2.0, -1.0, 0.0]

julia> scale_orthogonal!(A, collect(values(factors)))
3×3 Array{Float64,2}:
  1.0   0.0   0.0
  1.0   1.0   1.0
 -1.0  -1.0  -1.0

julia> A
3×3 Array{Float64,2}:
  1.0   0.0   0.0
  1.0   1.0   1.0
 -1.0  -1.0  -1.0

```

# Formula

For a design ``D_{n,p}`` with ``n`` experiments or rows and ``p`` factors
or columns, scale each factor ``\\mathbf{x}_i`` to ``\\mathbf{x}_{i}^{s}``
according to:

```math
\\mathbf{x}_{i}^{s} = \\dfrac{\\mathbf{x}_i - \\bar{M}}{\\bar{M}_{def}}
```

Where ``\\mathbf{x}_{def}`` is the factor defined in the `factors` parameter and:

```math
\\bar{M} = (max(\\mathbf{x}_i) + min(\\mathbf{x}_i)) / 2
```

```math
\\bar{M}_{def} = (max(\\mathbf{x}_{def}) - min(\\mathbf{x}_{def})) / 2
```
"""
function scale_orthogonal!(design::Array{Float64, 2},
                           factors::Array{T, 1}) where T <: Any
    for i = 1:size(design, 2)
        design_range = (max(design[:, i]...) + min(design[:, i]...)) / 2
        factor_range = (max(factors[i]...) - min(factors[i]...)) / 2
        design[:, i] = (design[:, i] .- design_range) ./ factor_range
    end

    return design
end

"""
    scale_boxdraper_encoding!(design::Array{Float64, 2},
                              factors::Array{T, 1};
                              scale_denominator = true) where T <: Any

Scale factors of a design using the Box and Draper's coding convention from
"Response Surfaces, Mixtures, and Ridge Analyses".

# Examples

With denominator scaling:

```jldoctest
julia> using ExperimentalDesign, DataStructures, StatsModels

julia> A = float(plackett_burman(4))
4×3 Array{Float64,2}:
  1.0   1.0   1.0
  1.0  -1.0  -1.0
 -1.0  -1.0   1.0
 -1.0   1.0  -1.0

julia> factors = OrderedDict([(:f1, [-1., 1.]), (:f2, [-1., 1.]), (:f3, [-1., 1.])])
DataStructures.OrderedDict{Symbol,Array{Float64,1}} with 3 entries:
  :f1 => [-1.0, 1.0]
  :f2 => [-1.0, 1.0]
  :f3 => [-1.0, 1.0]

julia> scale_boxdraper_encoding!(A, collect(values(factors)))
4×3 Array{Float64,2}:
  1.0   1.0   1.0
  1.0  -1.0  -1.0
 -1.0  -1.0   1.0
 -1.0   1.0  -1.0

julia> all([isapprox(1.0, sqrt(sum(A[:, i] .^ 2.0) / size(A, 1))) for i = 1:size(A, 2)])
true

```

Without denominator scaling:

```jldoctest
julia> using ExperimentalDesign, DataStructures, StatsModels

julia> A = float(plackett_burman(4))
4×3 Array{Float64,2}:
  1.0   1.0   1.0
  1.0  -1.0  -1.0
 -1.0  -1.0   1.0
 -1.0   1.0  -1.0

julia> factors = OrderedDict([(:f1, [-1., 1.]), (:f2, [-1., 1.]), (:f3, [-1., 1.])])
DataStructures.OrderedDict{Symbol,Array{Float64,1}} with 3 entries:
  :f1 => [-1.0, 1.0]
  :f2 => [-1.0, 1.0]
  :f3 => [-1.0, 1.0]

julia> scale_boxdraper_encoding!(A, collect(values(factors)), scale_denominator = false)
4×3 Array{Float64,2}:
  0.5   0.5   0.5
  0.5  -0.5  -0.5
 -0.5  -0.5   0.5
 -0.5   0.5  -0.5

julia> all([isapprox(1.0, sqrt(sum(A[:, i] .^ 2.0))) for i = 1:size(A, 2)])
true

```

# Formula

For a design ``D_{n,p}`` with ``n`` experiments or rows and ``p`` factors
or columns, scale each factor ``\\mathbf{x}_i`` to ``\\mathbf{x}_{i}^{s}``
according to:

```math
\\mathbf{x}_{i}^{s} = \\dfrac{\\mathbf{x}_i - \\bar{\\mathbf{x}}_i}{S_i}
```

Where ``\\bar{\\mathbf{x}}_i`` is mean of factor definition values in the
`factors` parameter and:

```math
S_{i}^{2} = \\dfrac{1}{n} \\sum\\limits_{j = 1}^{n}{(x_{ij} - \\bar{\\mathbf{x}}_i)^{2}}
```

If we pass `scale_denominator = false`, ``S_i`` becomes:

```math
S_{i}^{2} = \\sum\\limits_{j = 1}^{n}{(x_{ij} - \\bar{\\mathbf{x}}_i)^{2}}
```
"""
function scale_boxdraper_encoding!(design::Array{Float64, 2},
                                   factors::Array{T, 1};
                                   scale_denominator = true) where T <: Any
    for i = 1:size(design, 2)
        factor_mean = mean(factors[i])
        denominator = sum((design[:, i] .- factor_mean) .^ 2.0)

        if scale_denominator
            denominator /= size(design, 1)
        end

        denominator = sqrt(denominator)

        numerator = design[:, i] .- factor_mean

        if !iszero(denominator) && !iszero(numerator)
            numerator ./= denominator
        end

        design[:, i] = numerator
    end

    return design
end

"""
    get_expanded_values(factors::OrderedDict)

Return a mapping from an original level of a factor to a column
in an expanded design. Does not expand factors.

# Examples

```jldoctest
julia> using ExperimentalDesign, DataStructures, StatsModels

julia> A = OrderedDict(:A => [1., 2.], :B => [1, 2, 3, 4], :C => ["A", "B"], :D => [1.2, 2.2, 3.4], :E => [:a, :b, :c, :d])
DataStructures.OrderedDict{Symbol,Any} with 5 entries:
  :A => [1.0, 2.0]
  :B => [1, 2, 3, 4]
  :C => String["A", "B"]
  :D => [1.2, 2.2, 3.4]
  :E => Symbol[:a, :b, :c, :d]

julia> get_expanded_values(A)
DataStructures.OrderedDict{Any,Any} with 2 entries:
  :C => DataStructures.OrderedDict{Any,Symbol}("B"=>:C_1)
  :E => DataStructures.OrderedDict{Any,Symbol}(:b=>:E_1,:c=>:E_2,:d=>:E_3)

```
"""
function get_expanded_values(factors::OrderedDict)
    expanded_values = OrderedDict()
    for factor in factors
        if !(eltype(factor[2]) <: Real)
            expanded_values[factor[1]] = OrderedDict{Any, Symbol}()
            for l = 2:length(factor[2])
                expanded_values[factor[1]][factor[2][l]] = Symbol(factor[1], "_$(l - 1)")
            end
        end
    end
    return expanded_values
end

function expand_design(design::Array{T, 2}, factors::OrderedDict) where T <: Any
    small_design = DataFrame(design)
    rename!(small_design, OrderedDict(zip(names(small_design), keys(factors))))

    exp_values = get_expanded_values(factors)
    expanded_design = DataFrame()

    for name in names(small_design)
        if name in keys(exp_values)
            for key in keys(exp_values[name])
                expanded_design[exp_values[name][key]] = [small_design[i, name] == key ? 1.0 : 0. for i = 1:length(small_design[name])]
            end
        else
            expanded_design[name] = small_design[name]
        end
    end

    return expanded_design
end

"""
    expand_factors(factors::OrderedDict)

Return an `OrderedDict` with true factors expanded to `[0., 1.]` 2-level numerical factors.

For a true factor with ``n`` levels, creates ``n - 1`` 2-level factors. Only
one of these new factors can be at level `1.` at each row of a design, encoding
each level of the original factor. The first level of the original factor is
encoded by all new factors being at level `0.`.

# Examples
```jldoctest
julia> using ExperimentalDesign, DataStructures, StatsModels

julia> A = OrderedDict(:A => [1., 2.], :B => [1, 2, 3, 4], :C => ["A", "B"], :D => [1.2, 2.2, 3.4], :E => [:a, :b, :c, :d])
DataStructures.OrderedDict{Symbol,Any} with 5 entries:
  :A => [1.0, 2.0]
  :B => [1, 2, 3, 4]
  :C => String["A", "B"]
  :D => [1.2, 2.2, 3.4]
  :E => Symbol[:a, :b, :c, :d]

julia> expand_factors(A)
DataStructures.OrderedDict{Symbol,Any} with 7 entries:
  :A   => [1.0, 2.0]
  :B   => [1, 2, 3, 4]
  :C_1 => [0.0, 1.0]
  :D   => [1.2, 2.2, 3.4]
  :E_1 => [0.0, 1.0]
  :E_2 => [0.0, 1.0]
  :E_3 => [0.0, 1.0]

```
"""
function expand_factors(factors::OrderedDict)
    expanded_factors = OrderedDict{Symbol, Any}()
    for factor in factors
        if !(eltype(factor[2]) <: Real)
            for l = 1:(length(factor[2]) - 1)
                expanded_factors[Symbol(factor[1], "_$l")] = [0., 1.]
            end
        else
            expanded_factors[factor[1]] = factor[2]
        end
    end
    return expanded_factors
end

"""
    generate_model_matrix(formula::Formula,
                          design::Array{Float64, 2},
                          factors::OrderedDict;
                          scale::Function = scale_boxdraper_encoding!)

Generate a `DataFrame` with a scaled model matrix for a given formula, design and factors.

Assumes that `formula` is a linear relationship between all the factors in `factors`.

# Examples

```jldoctest
julia> using ExperimentalDesign, DataStructures, StatsModels

julia> A = plackett_burman(4)
4×3 Array{Int64,2}:
  1   1   1
  1  -1  -1
 -1  -1   1
 -1   1  -1

julia> factors = OrderedDict([(:f1, [-1., 1.]), (:f2, [-1., 1.]), (:f3, [-1., 1.])])
DataStructures.OrderedDict{Symbol,Array{Float64,1}} with 3 entries:
  :f1 => [-1.0, 1.0]
  :f2 => [-1.0, 1.0]
  :f3 => [-1.0, 1.0]

julia> m = generate_model_matrix(@formula(y ~ f1 + f2 + f3), float(A), factors)
4×4 Array{Float64,2}:
 1.0   1.0   1.0   1.0
 1.0   1.0  -1.0  -1.0
 1.0  -1.0  -1.0   1.0
 1.0  -1.0   1.0  -1.0

```
"""
function generate_model_matrix(formula::Formula,
                               input_design::DataFrame,
                               factors::OrderedDict;
                               scale::Function = scale_boxdraper_encoding!)
    design = deepcopy(input_design)
    contrasts = []

    for factor in keys(factors)
        if !(eltype(factors[factor]) <: Real)
            push!(contrasts,
                  (factor, ContrastsMatrix(DummyCoding(), factors[factor])))
        end
    end

    design[:y] = ones(size(design, 1))

    model_frame = ModelFrame(formula, design, contrasts = Dict(contrasts))

    return ModelMatrix(model_frame).m
end


text:
using Weave
set_chunk_defaults( Dict{Symbol, Any}(:eval=>false) )

for ch = filter(s->occursin(r"chapter..",s),readdir())
    for nb = filter(s->occursin(r".ipynb$",s),readdir(ch))
        weave(joinpath(ch,nb),
        out_path=joinpath("markdown",ch),
        doctype="pandoc")
    end
end

max_stars_repo_name: snowdj/fnc-extras

max_issues_repo_name: snowdj/fnc-extras

max_forks_repo_name: snowdj/fnc-extras


text:
abstract type AbstractWald <: SequentialSamplingModel
end

"""
# Wald Constructor
- `υ`: drift rate
- `α`: decision threshold
- `θ`: a encoding-response offset
## Usage
````julia
using SequentialSamplingModels
dist = Wald(υ=3.0, α=.5, θ=.130)
rt = rand(dist, 10)
like = pdf.(dist, rt)
loglike = logpdf.(dist, rt)
````
"""
struct Wald{T1,T2,T3} <: AbstractWald
    ν::T1
    α::T2
    θ::T3
end

Wald(;ν, α, θ) = Wald(ν, α, θ)

function pdf(d::AbstractWald, t::Float64)
    return pdf(InverseGaussian(d.α/d.ν, d.α^2), t - d.θ)
end

function logpdf(d::AbstractWald, t::Float64)
    return logpdf(InverseGaussian(d.α/d.ν, d.α^2), t - d.θ)
end

function logccdf(d::Wald, t::Float64)
    return logccdf(InverseGaussian(d.α/d.ν, d.α^2), t - d.θ)
end

function cdf(d::Wald, t::Float64)
    return cdf(InverseGaussian(d.α/d.ν, d.α^2), t - d.θ)
end

rand(d::AbstractWald) = rand(InverseGaussian(d.α/d.ν, d.α^2)) + d.θ

function rand(d::AbstractWald, n::Int)
    return rand(InverseGaussian(d.α/d.ν, d.α^2), n) .+ d.θ
end

mean(d::AbstractWald) = mean(InverseGaussian(d.α/d.ν, d.α^2)) + d.θ
std(d::AbstractWald) = std(InverseGaussian(d.α/d.ν, d.α^2))

"""
# WaldMixture Constructor
- `υ`: drift rate
- `σ`: standard deviation of drift rate
- `α`: decision threshold
- `θ`: a encoding-response offset
## Usage
````julia
using SequentialSamplingModels
dist = WaldMixture(υ=3.0, σ=.2, α=.5, θ=.130)
rt = rand(dist, 10)
like = pdf.(dist, rt)
loglike = logpdf.(dist, rt)
````
## References
Steingroever, H., Wabersich, D., & Wagenmakers, E. J. (2020).
Modeling across-trial variability in the Wald drift rate parameter.
Behavior Research Methods, 1-17.
"""
struct WaldMixture{T1,T2,T3,T4} <: AbstractWald
    ν::T1
    σ::T2
    α::T3
    θ::T4
end

WaldMixture(;ν, σ, α, θ) = WaldMixture(ν, σ, α, θ)

function pdf(d::WaldMixture, t::Float64)
    @unpack ν, σ, α ,θ = d
    c1 = α/√(2*π*(t - θ)^3)
    c2 = 1/cdf(Normal(0,1), ν/σ)
    c3 = exp(-(ν*(t - θ) - α)^2/(2*(t - θ)*((t - θ)*σ^2 + 1)))
    c4 = (α*σ^2 + ν)/√(σ^2*((t - θ)*σ^2 + 1))
    return c1 * c2 * c3 * cdf(Normal(0,1), c4)
end

function logpdf(d::WaldMixture, t::Float64)
    @unpack ν, σ, α ,θ = d
    c1 = log(α) - log(√(2*π*(t - θ)^3))
    c2 = log(1) - logcdf(Normal(0,1), ν/σ)
    c3 = -(ν*(t - θ) - α)^2/(2*(t - θ)*((t - θ)*σ^2 + 1))
    c4 = (α*σ^2 + ν)/√(σ^2*((t - θ)*σ^2 + 1))
    return c1 + c2 + c3 + logcdf(Normal(0,1), c4)
end

function rand(d::WaldMixture)
    x = rand(truncated(Normal(d.ν, d.σ), 0, Inf))
    return rand(InverseGaussian(d.α/x, d.α^2)) + d.θ
end

function rand(d::WaldMixture, n::Int)
    return map(x->rand(d), 1:n)
end


text:
function nlp_006_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - user defined functions

    model = Model(optimizer)

    function user_function_1d(x)
        if x >= 0
            return max(0.25, exp(x) - 1)
        else
            return log(-x + 3)
        end
    end
    JuMP.register(
        model,
        :user_function_1d,
        1,
        user_function_1d,
        autodiff = true,
    )

    function user_function_2d(x, y)
        return x^2 + y^3
    end
    JuMP.register(
        model,
        :user_function_2d,
        2,
        user_function_2d,
        autodiff = true,
    )

    @variable(model, x)
    @variable(model, y >= 0)

    @objective(model, Max, y + x)
    @NLconstraint(model, y >= user_function_1d(x))
    @NLconstraint(model, user_function_2d(x, y) <= 2)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 1.8813786425753092, tol = objective_tol)
    return check_solution(
        [x, y],
        [0.7546057578960682, 1.1267728846792409],
        tol = primal_tol,
    )
end


text:
function init_simulation_euler_gravity()
  mpi_isparallel() && error("coupled simulations are not yet tested for parallel runs") # TODO parallel

  # Print startup message
  print_startup_message()

  # Get number of dimensions
  ndims_ = parameter("ndims")::Int

  # Check if this is a restart from a previous result or a new simulation
  restart = parameter("restart", false)
  if restart
    error("restarting not yet supported for coupled Euler-gravity simulations")
  end

  # Initialize mesh
  begin
    print("Creating mesh... ")
    @timeit timer() "mesh creation" mesh = generate_mesh()
    mesh.current_filename = save_mesh_file(mesh)
    mesh.unsaved_changes = false
    println("done")
  end

  # Initialize system of equations
  print("Initializing system of equations... ")
  equations_name = parameter("equations")
  @assert equations_name == "euler_gravity" "This only works with 'euler_gravity' as equations type"
  equations_euler = make_equations("CompressibleEulerEquations", ndims_)
  equations_gravity = make_equations("HyperbolicDiffusionEquations", ndims_)
  println("done")

  # Initialize solver
  print("Initializing solver... ")
  solver_name = parameter("solver", valid=["dg"])
  solver = make_solver(solver_name, equations_euler, mesh)
  solver_gravity = make_solver(solver_name, equations_gravity, mesh,
                               surface_flux_function=flux_lax_friedrichs, volume_flux_function=flux_central)
  # `solver` = Euler solver -> this is to keep differences to original method to a minimum
  println("done")

  # Sanity checks
  # If DG volume integral type is weak form, volume flux type must be flux_central,
  # as everything else does not make sense
  if solver.volume_integral_type === Val(:weak_form) && solver.volume_flux_function !== flux_central
    error("using the weak formulation with a volume flux other than 'flux_central' does not make sense")
  end

  # Initialize solution
  amr_interval = parameter("amr_interval", 0)
  adapt_initial_conditions = parameter("adapt_initial_conditions", true)
  adapt_initial_conditions_only_refine = parameter("adapt_initial_conditions_only_refine", true)
  begin
    print("Applying initial conditions... ")
    t_start = parameter("t_start")
    time = t_start
    step = 0
    set_initial_conditions!(solver, time)
    set_initial_conditions!(solver_gravity, time)
    println("done")

    # If AMR is enabled, adapt mesh and re-apply ICs
    if amr_interval > 0 && adapt_initial_conditions
      @timeit timer() "initial condition AMR" has_changed = adapt!(mesh, solver, time,
          only_refine=adapt_initial_conditions_only_refine, passive_solvers=(solver_gravity,))

      # Iterate until mesh does not change anymore
      while has_changed
        set_initial_conditions!(solver, time)
        set_initial_conditions!(solver_gravity, time)
        @timeit timer() "initial condition AMR" has_changed = adapt!(mesh, solver, time,
            only_refine=adapt_initial_conditions_only_refine, passive_solvers=(solver_gravity,))
      end

      # Save mesh file
      mesh.current_filename = save_mesh_file(mesh)
      mesh.unsaved_changes = false
    end
  end
  t_end = parameter("t_end")

  # Init time integration
  time_integration_scheme = Symbol(parameter("time_integration_scheme", "timestep_carpenter_kennedy_erk54_2N!"))
  time_integration_function = eval(time_integration_scheme)

  # Print setup information
  solution_interval = parameter("solution_interval", 0)
  restart_interval = parameter("restart_interval", 0)
  polydeg = parameter("polydeg") # FIXME: This is currently the only DG-specific code in here
  n_steps_max = parameter("n_steps_max")
  cfl = parameter("cfl")
  sources = parameter("sources", "none")
  n_leaf_cells = count_leaf_cells(mesh.tree)
  min_level = minimum_level(mesh.tree)
  max_level = maximum_level(mesh.tree)
  domain_center = mesh.tree.center_level_0
  domain_length = mesh.tree.length_level_0
  min_dx = domain_length / 2^max_level
  max_dx = domain_length / 2^min_level
  s = ""
  s *= """| Simulation setup (Euler + Gravity)
          | ----------------
          | working directory:  $(pwd())
          | parameters file:    $(parameter("parameters_file"))
          | equations:          $(equations_name)
          | | Euler:
          | | | #variables:     $(nvariables(equations_euler))
          | | | variable names: $(join(varnames_cons(equations_euler), ", "))
          | | | sources:        $sources
          | | Gravity:
          | | | #variables:     $(nvariables(equations_gravity))
          | | | variable names: $(join(varnames_cons(equations_gravity), ", "))
          | | | sources:        $sources
          | restart:            $(restart ? "yes" : "no")
          """
  begin
    s *= "| initial conditions: $(get_name(solver.initial_conditions))\n"
    s *= "| t_start:            $t_start\n"
  end
  s *= """| t_end:              $t_end
          | AMR:                $(amr_interval > 0 ? "yes" : "no")
          """
  if amr_interval > 0
    s *= "| | AMR interval:     $amr_interval\n"
    s *= "| | adapt ICs:        $(adapt_initial_conditions ? "yes" : "no")\n"
  end
  s *= """| n_steps_max:        $n_steps_max
          | time integration:   $(get_name(time_integration_function))
          | restart interval:   $restart_interval
          | solution interval:  $solution_interval
          | #parallel threads:  $(Threads.nthreads())
          |
          | Solver
          | | solver:           $solver_name
          | | polydeg:          $polydeg
          | | CFL:              $cfl
          | | Euler solver:
          | | | volume integral:  $(get_name(solver.volume_integral_type))
          | | | volume flux:      $(get_name(solver.volume_flux_function))
          | | | surface flux:     $(get_name(solver.surface_flux_function))
          | | Gravity solver:
          | | | volume integral:  $(get_name(solver_gravity.volume_integral_type))
          | | | volume flux:      $(get_name(solver_gravity.volume_flux_function))
          | | | surface flux:     $(get_name(solver_gravity.surface_flux_function))
          | | #elements:        $(solver.n_elements)
          | | #interfaces:      $(solver.n_interfaces)
          | | #boundaries:      $(solver.n_boundaries)
          | | #l2mortars:       $(solver.n_l2mortars)
          | | #DOFs:            $(ndofs(solver)) + $(ndofs(solver_gravity))
          |
          | Mesh
          | | #cells:           $(length(mesh.tree))
          | | #leaf cells:      $n_leaf_cells
          | | minimum level:    $min_level
          | | maximum level:    $max_level
          | | domain center:    $(join(domain_center, ", "))
          | | domain length:    $domain_length
          | | minimum dx:       $min_dx
          | | maximum dx:       $max_dx
          """
  println()
  println(s)

  # Set up main loop
  save_final_solution = parameter("save_final_solution", true)
  save_final_restart = parameter("save_final_restart", true)
  analysis_interval = parameter("analysis_interval", 0)
  if analysis_interval > 0
    alive_interval = parameter("alive_interval", div(analysis_interval, 10))
  else
    alive_interval = 0
  end

  # Save initial conditions if desired
  if !restart && parameter("save_initial_solution", true)
    # we need to make sure, that derived quantities, such as e.g. blending
    # factor is already computed for the initial condition
    @notimeit timer() rhs!(solver, time)
    save_solution_file(solver, mesh, time, 0, step, "euler")

    @notimeit timer() rhs!(solver_gravity, time)
    save_solution_file(solver_gravity, mesh, time, 0, step, "gravity")
  end
  # Print initial solution analysis and initialize solution analysis
  if analysis_interval > 0
    if get_name(solver.initial_conditions) == "initial_conditions_eoc_test_coupled_euler_gravity"
      analyze_solution(solver, mesh, time, 0, step, 0, 0, solver_gravity=solver_gravity)
      # comment out for anything other than coupling convergence test
      println()
      analyze_solution(solver_gravity, mesh, time, 0, step, 0, 0)
    else
      analyze_solution(solver, mesh, time, 0, step, 0, 0, solver_gravity=solver_gravity)
    end
  end

  time_parameters = (time=time, step=step, t_end=t_end, cfl=cfl,
                    n_steps_max=n_steps_max,
                    save_final_solution=save_final_solution,
                    save_final_restart=save_final_restart,
                    analysis_interval=analysis_interval,
                    alive_interval=alive_interval,
                    solution_interval=solution_interval,
                    amr_interval=amr_interval,
                    restart_interval=restart_interval)
  return mesh, (solver, solver_gravity), time_parameters, time_integration_function
end

function run_simulation_euler_gravity(mesh, solvers, time_parameters, time_integration_function)
  mpi_isparallel() && error("coupled simulations are not yet tested for parallel runs") # TODO parallel

  @unpack time, step, t_end, cfl, n_steps_max,
          save_final_solution, save_final_restart,
          analysis_interval, alive_interval,
          solution_interval, amr_interval,
          restart_interval = time_parameters

  solver_euler, solver_gravity = solvers
  solver = solver_euler # To keep differences to `run_simulation` to a minimum

  loop_start_time = time_ns()
  analysis_start_time = time_ns()
  output_time = 0.0
  n_analysis_timesteps = 0

  # Declare variables to return error norms calculated in main loop
  # (needed such that the values are accessible outside of the main loop)
  local l2_error, linf_error

  # Start main loop (loop until final time step is reached)
  globals[:gravity_subcycles] = 0
  finalstep = false
  first_loop_iteration = true
  @timeit timer() "main loop" while !finalstep
    # Calculate time step size
    @timeit timer() "calc_dt" dt = calc_dt(solver, cfl)

    # Abort if time step size is NaN
    if isnan(dt)
      error("time step size `dt` is NaN")
    end

    # If the next iteration would push the simulation beyond the end time, set dt accordingly
    if time + dt > t_end
      dt = t_end - time
      finalstep = true
    end

    # Evolve solution by one time step
    @timeit timer() "timestep_euler_gravity!" begin
      timestep_euler_gravity!(solver_euler, solver_gravity, time, dt, time_parameters)
    end
    step += 1
    time += dt
    n_analysis_timesteps += 1

    # Check if we reached the maximum number of time steps
    if step == n_steps_max
      finalstep = true
    end

    # Analyze solution errors
    if analysis_interval > 0 && (step % analysis_interval == 0 || finalstep)
      # Calculate absolute and relative runtime
      runtime_absolute = (time_ns() - loop_start_time) / 10^9
      runtime_relative = ((time_ns() - analysis_start_time - output_time) / 10^9 /
                          (n_analysis_timesteps * ndofs(solver)))

      # Analyze solution
      if get_name(solver_euler.initial_conditions) == "initial_conditions_eoc_test_coupled_euler_gravity"
        l2_euler, linf_euler = @timeit timer() "analyze solution" analyze_solution(
            solver, mesh, time, dt, step, runtime_absolute, runtime_relative, solver_gravity=solver_gravity)
        # Pull gravity solver information from file
        timestep_gravity = eval(Symbol(parameter("time_integration_scheme_gravity")))
        cfl_gravity = parameter("cfl_gravity")::Float64
        rho0 = parameter("rho0")::Float64
        G = parameter("G")::Float64
        gravity_parameters = (; timestep_gravity, cfl_gravity, rho0, G)
        update_gravity!(solver_gravity, solver_euler.elements.u, gravity_parameters)
        l2_hypdiff, linf_hypdiff = @timeit timer() "analyze solution" analyze_solution(
            solver_gravity, mesh, time, dt, step, runtime_absolute, runtime_relative)
        l2_error   = vcat(l2_euler  , l2_hypdiff)
        linf_error = vcat(linf_euler, linf_hypdiff)
      else
        l2_error, linf_error = @timeit timer() "analyze solution" analyze_solution(
            solver, mesh, time, dt, step, runtime_absolute, runtime_relative, solver_gravity=solver_gravity)
      end

      # Reset time and counters
      analysis_start_time = time_ns()
      output_time = 0.0
      n_analysis_timesteps = 0
      if finalstep
        println("-"^80)
        println("Trixi simulation run finished.    Final time: $time    Time steps: $step")
        println("-"^80)
        println()
      end
    elseif alive_interval > 0 && step % alive_interval == 0
      runtime_absolute = (time_ns() - loop_start_time) / 10^9
      @printf("#t/s: %6d | dt: %.4e | Sim. time: %.4e | Run time: %.4e s\n",
              step, dt, time, runtime_absolute)
    end

    # Write solution file
    if solution_interval > 0 && (
        step % solution_interval == 0 || (finalstep && save_final_solution))
      output_start_time = time_ns()
      @timeit timer() "I/O" begin
        # Compute current AMR indicator values such that it can be written to
        # the solution file for the current number of elements
        if amr_interval > 0
          calc_amr_indicator(solver, mesh, time)
        end

        # If mesh has changed, write a new mesh file name
        if mesh.unsaved_changes
          mesh.current_filename = save_mesh_file(mesh, step)
          mesh.unsaved_changes = false
        end

        # Then write solution file
        save_solution_file(solver, mesh, time, dt, step, "euler")
        save_solution_file(solver_gravity, mesh, time, dt, step, "gravity")
      end
      output_time += time_ns() - output_start_time
    end

    # Perform adaptive mesh refinement
    if amr_interval > 0 && (step % amr_interval == 0) && !finalstep
      @timeit timer() "AMR" has_changed = adapt!(mesh, solver, time,
                                                 passive_solvers=(solver_gravity,))

      # Store if mesh has changed to write changed mesh file before next restart/solution output
      if has_changed
        mesh.unsaved_changes = true
      end
    end

    # The following call ensures that when doing memory allocation
    # measurements, the memory allocations for JIT compilation are discarded
    # (since virtually all relevant methods have already been called by now)
    if first_loop_iteration
      clear_malloc_data()
      first_loop_iteration = false
    end
  end

  # Print timer information
  print_timer(timer(), title="Trixi.jl", allocations=true, linechars=:ascii, compact=false)
  println()

  # Return error norms for EOC calculation
  println("Number of gravity subcycles: ", globals[:gravity_subcycles])
  if get_name(solver_euler.initial_conditions) == "initial_conditions_eoc_test_coupled_euler_gravity"
    return l2_error, linf_error, vcat(varnames_cons(solver.equations),
                                      varnames_cons(solver_gravity.equations))
  else
    return l2_error, linf_error, varnames_cons(solver.equations)
  end
end


text:
using DeformableRegistration.Transformation
 using DeformableRegistration.ImageProcessing
 using Base.Test

@testset "Transformation" begin
@testset "grid creation" begin
    # test getCellCenteredGrid
    spatialDomain = [-0.1,5.2,-1.1,3.3]
    gridSize=(5,8)
    voxelsize = [(spatialDomain[2]-spatialDomain[1])/gridSize[1],(spatialDomain[4]-spatialDomain[3])/gridSize[2]]
    shift = spatialDomain[1:2:3]
    centeredGrid = getCellCenteredGrid(voxelsize, shift ,gridSize)

    @test length(centeredGrid.data) == prod(gridSize)*2

    @test centeredGrid.data[1] == (shift[1] + voxelsize[1]/2)
    @test centeredGrid.data[prod(gridSize)+1] == (shift[2]+voxelsize[2]/2)

    nodalGrid = getNodalGrid(voxelsize,shift,gridSize)
    @test length(nodalGrid.data) == (gridSize[1]+1)*(gridSize[2]+1)*2

end
end

# @testset "stg <-> cen conversion" begin
#     gridSize=(5,8)
#     voxelsize = [1.06,0.55]
#     shift = [-0.1,-1.1]
#     centeredGrid = getCellCenteredGrid(voxelsize, shift ,gridSize)
#     staggeredGrid = getStaggeredGrid(voxelsize,shift,gridSize)
#     centeredGridFromStaggered = stg2cen(staggeredGrid,gridSize)
#     @test_approx_eq centeredGrid centeredGridFromStaggered
#     # test cen2stg
#     staggeredGridFromCenteredGrid = cen2stg(centeredGrid,gridSize)
# end

# visualize cell centered grid
#using DeformableRegistration.Visualization
#using PyPlot
#pygui(true); close("all"); PyPlot.svg(true)
#imgdata = rand(gridSize[1],gridSize[2])
#img = createImage(imgdata,spatialDomain=spatialDomain)
#centeredGrid = getCellCenteredGrid(img)
#figure()
#showImage(img)
#showGrid(spatialDomain,gridSize,centeredGrid,showPoints=true,showIndices=true)


text:
#
# Definitions for conic sections within
# R2BP dynamics.
#

"""
An abstract type for all conic sections.
"""
abstract type ConicSection end

"""
An abstract type representing the Circular
conic section.
"""
abstract type Circular <: ConicSection end

"""
An abstract type representing the Elliptical
conic section.
"""
abstract type Elliptical <: ConicSection end

"""
An abstract type representing the Parabolic
conic section.
"""
abstract type Parabolic <: ConicSection end

"""
An abstract type representing the Hyperbolic
conic section.
"""
abstract type Hyperbolic <: ConicSection end


text:
module TestArray

using Base.Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType
using Compat

typealias String Compat.ASCIIString

for isordered in (false, true)
    for R in (CategoricalArrays.DefaultRefType, UInt8, UInt, Int8, Int)
        # Vector
        a = ["b", "a", "b"]
        x = CategoricalVector{String, R}(a, ordered=isordered)

        @test x == a
        @test ordered(x) === isordered
        @test levels(x) == sort(unique(a))
        @test size(x) === (3,)
        @test length(x) === 3

        @test convert(CategoricalArray, x) === x
        @test convert(CategoricalArray{String}, x) === x
        @test convert(CategoricalArray{String, 1}, x) === x
        @test convert(CategoricalArray{String, 1, R}, x) === x
        @test convert(CategoricalArray{String, 1, DefaultRefType}, x) == x
        @test convert(CategoricalArray{String, 1, UInt8}, x) == x

        @test convert(CategoricalVector, x) === x
        @test convert(CategoricalVector{String}, x) === x
        @test convert(CategoricalVector{String, R}, x) === x
        @test convert(CategoricalVector{String, DefaultRefType}, x) == x
        @test convert(CategoricalVector{String, UInt8}, x) == x

        for y in (CategoricalArray(x, ordered=isordered),
                  CategoricalArray{String}(x, ordered=isordered),
                  CategoricalArray{String, 1}(x, ordered=isordered),
                  CategoricalArray{String, 1, R}(x, ordered=isordered),
                  CategoricalArray{String, 1, DefaultRefType}(x, ordered=isordered),
                  CategoricalArray{String, 1, UInt8}(x, ordered=isordered),
                  CategoricalVector(x, ordered=isordered),
                  CategoricalVector{String}(x, ordered=isordered),
                  CategoricalVector{String, R}(x, ordered=isordered),
                  CategoricalVector{String, DefaultRefType}(x, ordered=isordered),
                  CategoricalVector{String, UInt8}(x, ordered=isordered),
                  categorical(x, ordered=isordered),
                  categorical(x, false, ordered=isordered),
                  categorical(x, true, ordered=isordered))
            @test isa(y, CategoricalVector{String})
            @test ordered(y) === isordered
            @test y == x
            @test y !== x
            @test y.refs !== x.refs
            @test y.pool !== x.pool
        end

        @test convert(CategoricalArray, a) == x
        @test convert(CategoricalArray{String}, a) == x
        @test convert(CategoricalArray{String, 1}, a) == x
        @test convert(CategoricalArray{String, 1, R}, a) == x
        @test convert(CategoricalArray{String, 1, DefaultRefType}, a) == x
        @test convert(CategoricalArray{String, 1, UInt8}, a) == x

        @test convert(CategoricalVector, a) == x
        @test convert(CategoricalVector{String}, a) == x
        @test convert(CategoricalVector{String, R}, a) == x
        @test convert(CategoricalVector{String, DefaultRefType}, a) == x
        @test convert(CategoricalVector{String, UInt8}, a) == x

        @test CategoricalArray{String}(a, ordered=isordered) == x
        @test CategoricalArray{String, 1}(a, ordered=isordered) == x
        @test CategoricalArray{String, 1, R}(a, ordered=isordered) == x
        @test CategoricalArray{String, 1, DefaultRefType}(a, ordered=isordered) == x
        @test CategoricalArray{String, 1, UInt8}(a, ordered=isordered) == x

        @test CategoricalVector(a, ordered=isordered) == x
        @test CategoricalVector{String}(a, ordered=isordered) == x
        @test CategoricalVector{String, R}(a, ordered=isordered) == x
        @test CategoricalVector{String, DefaultRefType}(a, ordered=isordered) == x
        @test CategoricalVector{String, UInt8}(a, ordered=isordered) == x

        for (y, R1, R2, compact) in ((a, DefaultRefType, UInt8, true),
                                     (a, DefaultRefType, DefaultRefType, false),
                                     (x, R, UInt8, true),
                                     (x, R, R, false))
            x2 = categorical(y, ordered=isordered)
            @test x2 == x
            @test isa(x2, CategoricalVector{String, R1})
            @test ordered(x2) === isordered

            x2 = categorical(y, compact, ordered=isordered)
            @test x2 == x
            @test isa(x2, CategoricalVector{String, R2})
            @test ordered(x2) === isordered
        end

        x2 = compact(x)
        @test x2 == x
        @test isa(x2, CategoricalVector{String, UInt8})
        @test ordered(x2) === ordered(x)
        @test levels(x2) == levels(x)

        x2 = copy(x)
        @test x2 == x
        @test typeof(x2) === typeof(x)
        @test ordered(x2) === ordered(x)
        @test levels(x2) == levels(x)

        if !ordered(x)
            @test ordered!(x, true) === true
        end
        @test x[1] > x[2]
        @test x[3] > x[2]

        @test ordered!(x, false) === false
        @test ordered(x) === false
        @test_throws Exception x[1] > x[2]
        @test_throws Exception x[3] > x[2]

        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[1]
        @test_throws BoundsError x[4]

        @test x[1:2] == ["b", "a"]
        @test typeof(x[1:2]) === typeof(x)

        x[1] = x[2]
        @test x[1] === x.pool.valindex[2]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[1]

        x[3] = "c"
        @test x[1] === x.pool.valindex[2]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[3]
        @test levels(x) == ["a", "b", "c"]

        x[2:3] = "b"
        @test x[1] === x.pool.valindex[2]
        @test x[2] === x.pool.valindex[1]
        @test x[3] === x.pool.valindex[1]
        @test levels(x) == ["a", "b", "c"]

        droplevels!(x) == ["a", "b"]
        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[2]
        @test levels(x) == ["a", "b"]

        levels!(x, ["b", "a"]) == ["b", "a"]
        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[2]
        @test levels(x) == ["b", "a"]

        @test_throws ArgumentError levels!(x, ["a"])
        @test_throws ArgumentError levels!(x, ["e", "b"])
        @test_throws ArgumentError levels!(x, ["e", "a", "b", "a"])

        @test levels!(x, ["e", "a", "b"]) == ["e", "a", "b"]
        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[2]
        @test levels(x) == ["e", "a", "b"]

        x[1] = "c"
        @test x[1] === x.pool.valindex[4]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[2]
        @test levels(x) == ["e", "a", "b", "c"]

        push!(x, "a")
        @test length(x) == 4
        @test x[end] == "a"
        @test levels(x) == ["e", "a", "b", "c"]

        push!(x, "zz")
        @test length(x) == 5
        @test x[end] == "zz"
        @test levels(x) == ["e", "a", "b", "c", "zz"]

        push!(x, x[1])
        @test length(x) == 6
        @test x[1] == x[end]
        @test levels(x) == ["e", "a", "b", "c", "zz"]

        append!(x, x)
        @test length(x) == 12
        @test x == ["c", "b", "b", "a", "zz", "c", "c", "b", "b", "a", "zz", "c"]
        @test ordered(x) === false
        @test levels(x) == ["e", "a", "b", "c", "zz"]

        b = ["z","y","x"]
        y = CategoricalVector{String, R}(b)
        append!(x, y)
        @test length(x) == 15
        @test x == ["c", "b", "b", "a", "zz", "c", "c", "b", "b", "a", "zz", "c", "z", "y", "x"]
        @test levels(x) == ["e", "a", "b", "c", "zz", "x", "y", "z"]

        empty!(x)
        @test length(x) == 0
        @test levels(x) == ["e", "a", "b", "c", "zz", "x", "y", "z"]

        # Vector created from range (i.e. non-Array AbstractArray),
        # direct conversion to a vector with different eltype
        a = 0.0:0.5:1.5
        x = CategoricalVector{Float64, R}(a, ordered=isordered)

        @test x == collect(a)
        @test ordered(x) === isordered
        @test levels(x) == unique(a)
        @test size(x) === (4,)
        @test length(x) === 4

        @test convert(CategoricalArray, x) === x
        @test convert(CategoricalArray{Float64}, x) === x
        @test convert(CategoricalArray{Float64, 1}, x) === x
        @test convert(CategoricalArray{Float64, 1, R}, x) === x
        @test convert(CategoricalArray{Float64, 1, DefaultRefType}, x) == x
        @test convert(CategoricalArray{Float64, 1, UInt8}, x) == x

        @test convert(CategoricalVector, x) === x
        @test convert(CategoricalVector{Float64}, x) === x
        @test convert(CategoricalVector{Float64, R}, x) === x
        @test convert(CategoricalVector{Float64, DefaultRefType}, x) == x
        @test convert(CategoricalVector{Float64, UInt8}, x) == x

        for y in (CategoricalArray(x, ordered=isordered),
                  CategoricalArray{Float64}(x, ordered=isordered),
                  CategoricalArray{Float64, 1}(x, ordered=isordered),
                  CategoricalArray{Float64, 1, R}(x, ordered=isordered),
                  CategoricalArray{Float64, 1, DefaultRefType}(x, ordered=isordered),
                  CategoricalArray{Float64, 1, UInt8}(x, ordered=isordered),
                  CategoricalVector(x, ordered=isordered),
                  CategoricalVector{Float64}(x, ordered=isordered),
                  CategoricalVector{Float64, R}(x, ordered=isordered),
                  CategoricalVector{Float64, DefaultRefType}(x, ordered=isordered),
                  CategoricalVector{Float64, UInt8}(x, ordered=isordered),
                  categorical(x, ordered=isordered),
                  categorical(x, false, ordered=isordered),
                  categorical(x, true, ordered=isordered))
            @test isa(y, CategoricalVector{Float64})
            @test ordered(y) === isordered
            @test y == x
            @test y !== x
            @test y.refs !== x.refs
            @test y.pool !== x.pool
        end

        @test convert(CategoricalArray, a) == x
        @test convert(CategoricalArray{Float64}, a) == x
        @test convert(CategoricalArray{Float32}, a) == x
        @test convert(CategoricalArray{Float64, 1}, a) == x
        @test convert(CategoricalArray{Float32, 1}, a) == x
        @test convert(CategoricalArray{Float64, 1, R}, a) == x
        @test convert(CategoricalArray{Float32, 1, R}, a) == x
        @test convert(CategoricalArray{Float64, 1, DefaultRefType}, a) == x
        @test convert(CategoricalArray{Float32, 1, DefaultRefType}, a) == x
        @test convert(CategoricalArray{Float64, 1, UInt8}, a) == x
        @test convert(CategoricalArray{Float32, 1, UInt8}, a) == x

        @test convert(CategoricalVector, a) == x
        @test convert(CategoricalVector{Float64}, a) == x
        @test convert(CategoricalVector{Float32}, a) == x
        @test convert(CategoricalVector{Float64, R}, a) == x
        @test convert(CategoricalVector{Float32, R}, a) == x
        @test convert(CategoricalVector{Float64, DefaultRefType}, a) == x
        @test convert(CategoricalVector{Float32, DefaultRefType}, a) == x
        @test convert(CategoricalVector{Float64, UInt8}, a) == x
        @test convert(CategoricalVector{Float32, UInt8}, a) == x

        @test CategoricalArray{Float64}(a, ordered=isordered) == x
        @test CategoricalArray{Float32}(a, ordered=isordered) == x
        @test CategoricalArray{Float64, 1}(a, ordered=isordered) == x
        @test CategoricalArray{Float32, 1}(a, ordered=isordered) == x
        @test CategoricalArray{Float64, 1, R}(a, ordered=isordered) == x
        @test CategoricalArray{Float32, 1, R}(a, ordered=isordered) == x
        @test CategoricalArray{Float64, 1, DefaultRefType}(a, ordered=isordered) == x
        @test CategoricalArray{Float32, 1, DefaultRefType}(a, ordered=isordered) == x

        @test CategoricalVector(a, ordered=isordered) == x
        @test CategoricalVector{Float64}(a, ordered=isordered) == x
        @test CategoricalVector{Float32}(a, ordered=isordered) == x
        @test CategoricalVector{Float64, R}(a, ordered=isordered) == x
        @test CategoricalVector{Float32, R}(a, ordered=isordered) == x
        @test CategoricalVector{Float64, DefaultRefType}(a, ordered=isordered) == x
        @test CategoricalVector{Float32, DefaultRefType}(a, ordered=isordered) == x

        for (y, R1, R2, compact) in ((a, DefaultRefType, UInt8, true),
                                     (a, DefaultRefType, DefaultRefType, false),
                                     (x, R, UInt8, true),
                                     (x, R, R, false))
            x2 = categorical(y, ordered=isordered)
            @test x2 == x
            @test isa(x2, CategoricalVector{Float64, R1})
            @test ordered(x2) === isordered

            x2 = categorical(y, compact, ordered=isordered)
            @test x2 == x
            @test isa(x2, CategoricalVector{Float64, R2})
            @test ordered(x2) === isordered
        end

        x2 = copy(x)
        @test x2 == x
        @test typeof(x2) === typeof(x)
        @test ordered(x2) === ordered(x)
        @test levels(x2) == levels(x)

        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[3]
        @test x[4] === x.pool.valindex[4]
        @test_throws BoundsError x[5]

        @test x[1:2] == [0.0, 0.5]
        @test typeof(x[1:2]) === typeof(x)

        x[2] = 1
        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[3]
        @test x[3] === x.pool.valindex[3]
        @test x[4] === x.pool.valindex[4]
        @test levels(x) == unique(a)

        x[1:2] = -1
        @test x[1] === x.pool.valindex[5]
        @test x[2] === x.pool.valindex[5]
        @test x[3] === x.pool.valindex[3]
        @test x[4] === x.pool.valindex[4]
        @test levels(x) == vcat(unique(a), -1)

        push!(x, 2.0)
        @test length(x) == 5
        @test x[end] == 2.0
        @test ordered(x) === isordered
        @test levels(x) == [0.0,  0.5,  1.0,  1.5, -1.0,  2.0]

        push!(x, x[1])
        @test length(x) == 6
        @test x[1] == x[end]
        @test ordered(x) === isordered
        @test levels(x) == [0.0,  0.5,  1.0,  1.5, -1.0,  2.0]

        append!(x, x)
        @test length(x) == 12
        @test x == [-1.0, -1.0, 1.0, 1.5, 2.0, -1.0, -1.0, -1.0, 1.0, 1.5, 2.0, -1.0]
        @test ordered(x) === isordered
        @test levels(x) == [0.0,  0.5,  1.0,  1.5, -1.0,  2.0]

        b = [2.5, 3.0, -3.5]
        y = CategoricalVector{Float64, R}(b, ordered=isordered)
        append!(x, y)
        @test length(x) == 15
        @test x == [-1.0, -1.0, 1.0, 1.5, 2.0, -1.0, -1.0, -1.0, 1.0, 1.5, 2.0, -1.0, 2.5, 3.0, -3.5]
        @test ordered(x) === isordered
        @test levels(x) == [0.0, 0.5, 1.0, 1.5, -1.0, 2.0, -3.5, 2.5, 3.0]

        empty!(x)
        @test length(x) == 0
        @test ordered(x) === isordered
        @test levels(x) == [0.0, 0.5, 1.0, 1.5, -1.0, 2.0, -3.5, 2.5, 3.0]

        # Matrix
        a = ["a" "b" "c"; "b" "a" "c"]
        x = CategoricalMatrix{String, R}(a, ordered=isordered)

        @test x == a
        @test ordered(x) === isordered
        @test levels(x) == unique(a)
        @test size(x) === (2, 3)
        @test length(x) === 6

        @test convert(CategoricalArray, x) === x
        @test convert(CategoricalArray{String}, x) === x
        @test convert(CategoricalArray{String, 2}, x) === x
        @test convert(CategoricalArray{String, 2, R}, x) === x
        @test convert(CategoricalArray{String, 2, DefaultRefType}, x) == x
        @test convert(CategoricalArray{String, 2, UInt8}, x) == x

        @test convert(CategoricalMatrix, x) === x
        @test convert(CategoricalMatrix{String}, x) === x
        @test convert(CategoricalMatrix{String, R}, x) === x
        @test convert(CategoricalMatrix{String, DefaultRefType}, x) == x
        @test convert(CategoricalMatrix{String, UInt8}, x) == x

        for y in (CategoricalArray(x, ordered=isordered),
                  CategoricalArray{String}(x, ordered=isordered),
                  CategoricalArray{String, 2}(x, ordered=isordered),
                  CategoricalArray{String, 2, R}(x, ordered=isordered),
                  CategoricalArray{String, 2, DefaultRefType}(x, ordered=isordered),
                  CategoricalArray{String, 2, UInt8}(x, ordered=isordered),
                  CategoricalMatrix(x, ordered=isordered),
                  CategoricalMatrix{String}(x, ordered=isordered),
                  CategoricalMatrix{String, R}(x, ordered=isordered),
                  CategoricalMatrix{String, DefaultRefType}(x, ordered=isordered),
                  CategoricalMatrix{String, UInt8}(x, ordered=isordered),
                  categorical(x, ordered=isordered),
                  categorical(x, false, ordered=isordered),
                  categorical(x, true, ordered=isordered))
            @test isa(y, CategoricalMatrix{String})
            @test ordered(y) === isordered
            @test y == x
            @test y !== x
            @test y.refs !== x.refs
            @test y.pool !== x.pool
        end

        @test convert(CategoricalArray, a) == x
        @test convert(CategoricalArray{String}, a) == x
        @test convert(CategoricalArray{String, 2, R}, a) == x
        @test convert(CategoricalArray{String, 2, DefaultRefType}, a) == x
        @test convert(CategoricalArray{String, 2, UInt8}, a) == x

        @test convert(CategoricalMatrix, a) == x
        @test convert(CategoricalMatrix{String}, a) == x
        @test convert(CategoricalMatrix{String, R}, a) == x
        @test convert(CategoricalMatrix{String, DefaultRefType}, a) == x
        @test convert(CategoricalMatrix{String, UInt8}, a) == x

        @test CategoricalArray{String}(a, ordered=isordered) == x
        @test CategoricalArray{String, 2}(a, ordered=isordered) == x
        @test CategoricalArray{String, 2}(a, ordered=isordered) == x
        @test CategoricalArray{String, 2, R}(a, ordered=isordered) == x
        @test CategoricalArray{String, 2, DefaultRefType}(a, ordered=isordered) == x
        @test CategoricalArray{String, 2, UInt8}(a, ordered=isordered) == x

        @test CategoricalMatrix(a, ordered=isordered) == x
        @test CategoricalMatrix{String}(a, ordered=isordered) == x
        @test CategoricalMatrix{String, R}(a, ordered=isordered) == x
        @test CategoricalMatrix{String, DefaultRefType}(a, ordered=isordered) == x
        @test CategoricalMatrix{String, UInt8}(a, ordered=isordered) == x

        for (y, R1, R2, compact) in ((a, DefaultRefType, UInt8, true),
                                     (a, DefaultRefType, DefaultRefType, false),
                                     (x, R, UInt8, true),
                                     (x, R, R, false))
            x2 = categorical(y, ordered=isordered)
            @test x2 == x
            @test isa(x2, CategoricalMatrix{String, R1})
            @test ordered(x2) === isordered

            x2 = categorical(y, compact, ordered=isordered)
            @test x2 == x
            @test isa(x2, CategoricalMatrix{String, R2})
            @test ordered(x2) === isordered
        end

        x2 = compact(x)
        @test x2 == x
        @test isa(x2, CategoricalMatrix{String, UInt8})
        @test ordered(x2) === ordered(x)
        @test levels(x2) == levels(x)

        x2 = copy(x)
        @test x2 == x
        @test typeof(x2) === typeof(x)
        @test ordered(x2) === ordered(x)
        @test levels(x2) == levels(x)

        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[2]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[3]
        @test x[6] === x.pool.valindex[3]
        @test_throws BoundsError x[7]

        @test x[1,1] === x.pool.valindex[1]
        @test x[2,1] === x.pool.valindex[2]
        @test x[1,2] === x.pool.valindex[2]
        @test x[2,2] === x.pool.valindex[1]
        @test x[1,3] === x.pool.valindex[3]
        @test x[2,3] === x.pool.valindex[3]
        @test_throws BoundsError x[1,4]
        @test_throws BoundsError x[4,1]
        @test_throws BoundsError x[4,4]

        @test x[1:2,:] == x
        @test typeof(x[1:2,:]) === typeof(x)
        @test x[1:2,1] == ["a", "b"]
        @test typeof(x[1:2,1]) === CategoricalVector{String, R}

        x[1] = "z"
        @test x[1] === x.pool.valindex[4]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[2]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[3]
        @test x[6] === x.pool.valindex[3]
        @test levels(x) == ["a", "b", "c", "z"]

        x[1,:] = "a"
        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[1]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[1]
        @test x[6] === x.pool.valindex[3]
        @test levels(x) == ["a", "b", "c", "z"]

        x[1,1:2] = "z"
        @test x[1] === x.pool.valindex[4]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[4]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[1]
        @test x[6] === x.pool.valindex[3]
        @test levels(x) == ["a", "b", "c", "z"]

        # Uninitialized array
        v = Any[CategoricalArray(2, ordered=isordered),
                CategoricalArray(String, 2, ordered=isordered),
                CategoricalArray{String}(2, ordered=isordered),
                CategoricalArray{String, 1}(2, ordered=isordered),
                CategoricalArray{String, 1, R}(2, ordered=isordered),
                CategoricalVector{String}(2, ordered=isordered),
                CategoricalVector{String, R}(2, ordered=isordered),
                CategoricalArray(2, 3, ordered=isordered),
                CategoricalArray(String, 2, 3, ordered=isordered),
                CategoricalArray{String}(2, 3, ordered=isordered),
                CategoricalArray{String, 2}(2, 3, ordered=isordered),
                CategoricalArray{String, 2, R}(2, 3, ordered=isordered),
                CategoricalMatrix{String}(2, 3, ordered=isordered),
                CategoricalMatrix{String, R}(2, 3, ordered=isordered)]

        # See conditional definition of constructors in array.jl
        if VERSION >= v"0.5.0-dev"
            push!(v, CategoricalVector(2, ordered=isordered),
                     CategoricalVector(String, 2, ordered=isordered),
                     CategoricalMatrix(2, 3, ordered=isordered),
                     CategoricalMatrix(String, 2, 3, ordered=isordered))
        end

        for x in v
            @test !isassigned(x, 1) && isdefined(x, 1)
            @test !isassigned(x, 2) && isdefined(x, 2)
            @test_throws UndefRefError x[1]
            @test_throws UndefRefError x[2]
            @test ordered(x) === isordered
            @test levels(x) == []

            x2 = compact(x)
            if VERSION >= v"0.5.0-dev"
                @test isa(x2, CategoricalArray{String, ndims(x), UInt8})
            else
                @test isa(x2, CategoricalArray{typeof(x).parameters[1], ndims(x), UInt8})
            end
            @test !isassigned(x2, 1) && isdefined(x2, 1)
            @test !isassigned(x2, 2) && isdefined(x2, 2)
            @test_throws UndefRefError x2[1]
            @test_throws UndefRefError x2[2]
            @test levels(x2) == []

            x[1] = "c"
            @test x[1] === x.pool.valindex[1]
            @test !isassigned(x, 2) && isdefined(x, 2)
            @test_throws UndefRefError x[2]
            @test levels(x) == ["c"]

            x[1] = "a"
            @test x[1] === x.pool.valindex[2]
            @test !isassigned(x, 2) && isdefined(x, 2)
            @test_throws UndefRefError x[2]
            @test ordered(x) === isordered
            @test levels(x) == ["c", "a"]

            x[2] = "c"
            @test x[1] === x.pool.valindex[2]
            @test x[2] === x.pool.valindex[1]
            @test levels(x) == ["c", "a"]

            x[1] = "b"
            @test x[1] === x.pool.valindex[3]
            @test x[2] === x.pool.valindex[1]
            @test levels(x) == ["c", "a", "b"]
        end
    end
end

end


text:
#RSDeltaSigmaPort: Generate plots of internal state
#-------------------------------------------------------------------------------

#==Plot state maxima
===============================================================================#
"""`plotStateMaxima(u, ABCD, test_sig=nothing; nlev::Int=2, N::Int=10^5)`

# Inputs
 - `N`: Number of points used to generate a test signal
   (only used if `test_sig=nothing`).
"""
function plotStateMaxima(u, ABCD, test_sig=nothing; nlev::Int=2, N::Int=10^5)
	title = "Simulated State Maxima"; color = :blue
	u = collect(u)
	if isnothing(test_sig)
		test_sig = ones(1,N)
	end
	nu, nq, order = get_nu_nq_order(test_sig, ABCD, nlev)
	maxima = zeros(order,length(u))

	for i = 1:length(u)
		ui = u[i]

		simresult = simulateDSM(ui*test_sig, ABCD, nlev=nlev, trackmax=true)
		maxima[:,i] = simresult.xmax[:]
		if any(simresult.xmax .> 1e2)
			@warn("plotStateMaxima(): umax was too high.")
			umax = ui
			u = u[1:i]
			maxima = maxima[:,1:i]
			break
		end
	end

	plot = cons(:plot, linlog, title = title, legend=true,
#		xyaxes=set(xmin=0, xmax=0.6, ymin=1e-4, ymax=10),
		xyaxes=set(ymin=1e-4, ymax=10),
		labels=set(xaxis="DC input", yaxis="Peak value"),
	)
	colors = range(HSV(0,1,1), stop=HSV(-360,1,1), length=order)
	for i in 1:order
		id = "state $i"; c = colors[i]
		simglyph = cons(:a, glyph=set(shape=:o, size=1, color=c, fillcolor=c))

		_maxima = waveform(u, maxima[i,:])
		push!(plot,
		cons(:wfrm, _maxima, simglyph, line=set(style=:dash, color=c, width=2), label=id),
		)
	end
	return plot
end

#Last line


text:
# wrapper to convert Nonlinear solver into LinearQuadratic solver

# To enable LPQP support from a Nonlinear solver, define, e.g.,
# LinearQuadraticModel(s::IpoptSolver) = NonlinearToLPQPBridge(NonlinearModel(s))

mutable struct QuadConstr
    linearidx::Vector{Int}
    linearval::Vector{Float64}
    quadrowidx::Vector{Int}
    quadcolidx::Vector{Int}
    quadval::Vector{Float64}
    sense::Char
    rhs::Float64
end
getdata(q::QuadConstr) = q.linearidx, q.linearval, q.quadrowidx, q.quadcolidx, q.quadval

mutable struct NonlinearToLPQPBridge <: AbstractLinearQuadraticModel
    nlpmodel::AbstractNonlinearModel
    A::Union{SparseMatrixCSC{Float64,Int64}, <: AbstractMatrix{Float64}}
    l::Vector{Float64}
    u::Vector{Float64}
    c::Vector{Float64}
    lb::Vector{Float64}
    ub::Vector{Float64}
    sense::Symbol
    Qobj::Tuple{Vector{Int},Vector{Int},Vector{Float64}}
    Qconstr::Vector{QuadConstr}
    numvar::Int
    numconstr::Int
    warmstart::Vector{Float64}
    vartypes::Vector{Symbol}
end

function NonlinearToLPQPBridge(m::AbstractNonlinearModel)
    return NonlinearToLPQPBridge(m,Matrix{Float64}(undef,0,0),Vector{Float64}(),Vector{Float64}(),Vector{Float64}(),Vector{Float64}(),Vector{Float64}(),:Min,(Int[],Int[],Float64[]),QuadConstr[],0,0,Float64[],Symbol[])
end

export NonlinearToLPQPBridge

function loadproblem!(model::NonlinearToLPQPBridge, A, l, u, c, lb, ub, sense)
    model.A = A
    model.l = l
    model.u = u
    model.c = c
    model.lb = lb
    model.ub = ub
    model.sense = sense
    model.numvar = size(A,2)
    model.numconstr = size(A,1)
end

function getvarLB(model::NonlinearToLPQPBridge)
    return model.l
end

function setvarLB!(model::NonlinearToLPQPBridge, l::Vector)
    if length(l) != model.numvar
        error("Invalid size of l")
    end
    model.l .= l
    return nothing
end

function getvarUB(model::NonlinearToLPQPBridge)
    return model.u
end

function setvarUB!(model::NonlinearToLPQPBridge, u::Vector)
    if length(u) != model.numvar
        error("Invalid size of u")
    end
    model.u .= u
    return nothing
end

function getconstrLB(model::NonlinearToLPQPBridge)
    return model.lb
end

function setconstrLB!(model::NonlinearToLPQPBridge, lb::Vector)
    if length(lb) != model.numconstr
        error("Invalid size of lb")
    end
    model.lb .= lb
    return nothing
end

function getconstrUB(model::NonlinearToLPQPBridge)
    return model.ub
end

function setconstrUB!(model::NonlinearToLPQPBridge, ub::Vector)
    if length(ub) != model.numconstr
        error("Invalid size of ub")
    end
    model.ub .= ub
    return nothing
end

function getobj(model::NonlinearToLPQPBridge)
    return model.c
end

function setobj!(model::NonlinearToLPQPBridge, c::Vector)
    if length(c) != model.numvar
        error("Invalid size of c")
    end
    model.c .= c
    return nothing
end

function getconstrmatrix(model::NonlinearToLPQPBridge)
    return model.A
end

function addvar!(model::NonlinearToLPQPBridge, constridx::Vector, coeffs::Vector, l::Real, u::Real, objcoeff::Real)
    if length(constridx) != length(coeffs)
        error("constridx and coeffs have different legths")
    end
    model.numvar += 1
    push!(model.l, l)
    push!(model.u, u)
    push!(model.c, objcoeff)
    _addvar!(model, model.A, constridx, coeffs)
    return nothing
end

function _addvar!(model::NonlinearToLPQPBridge, A::SparseMatrixCSC, constridx::Vector, coeffs::Vector)
    I,J,V = findnz(A)
    for (idx,coeff) in zip(constridx, coeffs)
        push!(I, idx)
        push!(J, model.numvar)
        push!(V, coeff)
    end
    model.A = sparse(I,J,V,model.numconstr,model.numvar)
    return nothing
end

function _addvar!(model::NonlinearToLPQPBridge, A::Matrix, constridx::Vector, coeffs::Vector)
    newcol = zeros(model.numconstr,1)
    newcol[constridx] = coeffs
    model.A = hcat(A,newcol)
    return nothing
end

function delvars!(model::NonlinearToLPQPBridge, rmidx::Vector)
    sort!(rmidx)
    rmidx = unique(rmidx)
    deleteat!(model.l, rmidx)
    deleteat!(model.u, rmidx)
    deleteat!(model.c, rmidx)
    delvars!(model, model.A, rmidx)
    return nothing
end

function delvars!(model::NonlinearToLPQPBridge, A::SparseMatrixCSC, rmidx::Vector)
    I = Vector{Int}()
    J = Vector{Int}()
    V = Vector{Float64}()
    rows = rowvals(A)
    vals = nonzeros(A)
    m, n = size(A)
    for j = 1:n
        for r in nzrange(A, j)
            i = rows[r]
            v = vals[r]
            if j ∉ rmidx
                push!(I, i)
                push!(J, j)
                push!(V, v)
            end
        end
    end
    n -= length(rmidx)
    model.numvar = n
    model.A = sparse(I,J,V,m,n)
    return nothing
end

function delvars!(model::NonlinearToLPQPBridge, A::Matrix, rmidx::Vector)
    model.A = A[:, setdiff(1:size(A,2), rmidx)]
    model.numvar -= length(rmidx)
    return nothing
end

function addconstr!(model::NonlinearToLPQPBridge, colidx::Vector, coeffs::Vector, lb::Real, ub::Real)
    if length(colidx) != length(coeffs)
        error("colidx and coeffs have different legths")
    end
    push!(model.lb, lb)
    push!(model.ub, ub)
    model.numconstr += 1
    addconstr!(model, model.A, colidx, coeffs)
    return nothing
end

function addconstr!(model::NonlinearToLPQPBridge, A::SparseMatrixCSC, colidx::Vector, coeffs::Vector)
    I,J,V = findnz(A)
    for (idx,coeff) in zip(colidx, coeffs)
        push!(I, model.numconstr)
        push!(J, idx)
        push!(V, coeff)
    end
    model.A = sparse(I,J,V,model.numconstr,model.numvar)
    return nothing
end

function addconstr!(model::NonlinearToLPQPBridge, A::Matrix, colidx::Vector, coeffs::Vector)
    isempty(colidx) && return nothing
    newrow = zeros(model.numvar)
    newrow[colidx] = coeffs
    model.A = vcat(A,newrow')
    return nothing
end

function delconstrs!(model::NonlinearToLPQPBridge, rmidx::Vector)
    sort!(rmidx)
    rmidx = unique(rmidx)
    deleteat!(model.lb, rmidx)
    deleteat!(model.ub, rmidx)
    delconstrs!(model, model.A, rmidx)
    return nothing
end

function delconstrs!(model::NonlinearToLPQPBridge, A::SparseMatrixCSC, rmidx::Vector)
    I = Vector{Int}()
    J = Vector{Int}()
    V = Vector{Float64}()
    m, n = size(A)
    rows = rowvals(A)
    vals = nonzeros(A)
    for j = 1:n
        for r in nzrange(model.A, j)
            i = rows[r]
            v = vals[r]
            pos = findlast(idx -> i >= idx, rmidx)
            if pos == nothing
                push!(I, i)
            elseif i == rmidx[pos]
                continue
            else
                push!(I, i-pos)
            end
            push!(J, j)
            push!(V, v)
        end
    end
    m -= length(rmidx)
    model.A = sparse(I,J,V,m,n)
    model.numconstr = m
end

function delconstrs!(model::NonlinearToLPQPBridge, A::Matrix, rmidx::Vector)
    model.A = A[setdiff(1:size(A,1), rmidx), :]
    model.numconstr -= length(rmidx)
    return nothing
end

function setsense!(model::NonlinearToLPQPBridge, sense)
    model.sense = sense
end

function setquadobj!(model::NonlinearToLPQPBridge, rowidx, colidx, quadval)
    model.Qobj = (rowidx, colidx, quadval)
end

function addquadconstr!(model::NonlinearToLPQPBridge, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
    push!(model.Qconstr, QuadConstr(linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs))
    model.numconstr += 1
end

function optimize!(model::NonlinearToLPQPBridge)
    lb = copy(model.lb)
    ub = copy(model.ub)
    for q in model.Qconstr
        if q.sense == '<'
            push!(lb,-Inf)
            push!(ub,q.rhs)
        elseif q.sense == '>'
            push!(lb,q.rhs)
            push!(ub,Inf)
        else
            @assert q.sense == '='
            push!(lb,q.rhs)
            push!(ub,q.rhs)
        end
    end

    loadproblem!(model.nlpmodel, model.numvar, model.numconstr, model.l, model.u, lb, ub, model.sense, LPQPEvaluator(model))
    if any(model.vartypes .!= :Cont)
        setvartype!(model.nlpmodel, model.vartypes)
    end
    if length(model.warmstart) > 0
        setwarmstart!(model.nlpmodel, model.warmstart)
    end
    optimize!(model.nlpmodel)
end

for f in [:status, :getsolution, :getobjval, :getobjbound, :getreducedcosts]
    @eval $f(model::NonlinearToLPQPBridge) = $f(model.nlpmodel)
end

getsense(model::NonlinearToLPQPBridge) = model.sense
numvar(model::NonlinearToLPQPBridge) = model.numvar
numconstr(model::NonlinearToLPQPBridge) = model.numconstr

numlinconstr(model::NonlinearToLPQPBridge) = size(model.A,1)
numquadconstr(model::NonlinearToLPQPBridge) = length(model.Qconstr)

setvartype!(model::NonlinearToLPQPBridge,vartypes::Vector{Symbol}) = (model.vartypes = vartypes)
setwarmstart!(model::NonlinearToLPQPBridge,v) = (model.warmstart = v)

function getconstrsolution(model::NonlinearToLPQPBridge)
    x = getsolution(model)
    return model.A*x
end

function getconstrduals(model::NonlinearToLPQPBridge)
    du = getconstrduals(model.nlpmodel)
    return du[1:size(model.A,1)]
end

function getquadconstrduals(model::NonlinearToLPQPBridge)
    du = getconstrduals(model.nlpmodel)
    return du[(numlinconstr(model)+1):end]
end

mutable struct LPQPEvaluator <: AbstractNLPEvaluator
    A::SparseMatrixCSC{Float64,Int}
    c::Vector{Float64}
    Qi::Vector{Int} # quadratic objective terms
    Qj::Vector{Int}
    Qv::Vector{Float64}
    Qconstr::Vector{QuadConstr}
end

function LPQPEvaluator(model::NonlinearToLPQPBridge)
    A = model.A
    c = model.c
    Qi,Qj,Qv = model.Qobj
    return LPQPEvaluator(A,c,Qi,Qj,Qv,model.Qconstr)
end

features_available(d::LPQPEvaluator) = [:Grad, :Jac, :Hess]

function initialize(d::LPQPEvaluator, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in features_available(d))
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
end

function eval_f(d::LPQPEvaluator, x)
    obj = dot(d.c,x)
    Qi = d.Qi
    Qj = d.Qj
    Qv = d.Qv
    for k in 1:length(Qi)
        if Qi[k] == Qj[k]
            obj += Qv[k]*x[Qi[k]]*x[Qi[k]]/2
        else
            obj += Qv[k]*x[Qi[k]]*x[Qj[k]]
        end
    end
    return obj
end

function eval_grad_f(d::LPQPEvaluator, grad_f, x)
    Qi = d.Qi
    Qj = d.Qj
    Qv = d.Qv
    for j = 1:length(d.c)
        grad_f[j] = d.c[j]
    end
    for k in 1:length(Qi)
        if Qi[k] == Qj[k]
            grad_f[Qi[k]] += Qv[k]*x[Qj[k]]
        else
            grad_f[Qi[k]] += Qv[k]*x[Qj[k]]
            grad_f[Qj[k]] += Qv[k]*x[Qi[k]]
        end
    end
end

function eval_g(d::LPQPEvaluator, g, x)
    g_val = d.A*x
    g[1:size(d.A,1)] = g_val
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        val = 0.0
        for k in 1:length(linearidx)
            val += linearval[k]*x[linearidx[k]]
        end
        for k in 1:length(quadrowidx)
            val += quadval[k]*x[quadrowidx[k]]*x[quadcolidx[k]]
        end
        g[i+size(d.A,1)] = val
    end
end

function jac_structure(d::LPQPEvaluator)
    jac_nnz = nnz(d.A)
    for i in 1:length(d.Qconstr)
        jac_nnz += length(d.Qconstr[i].linearidx)
        jac_nnz += 2*length(d.Qconstr[i].quadrowidx)
    end

    I = Vector{Int}(undef, jac_nnz)
    J = Vector{Int}(undef, jac_nnz)
    idx = 1
    for col = 1:size(d.A,2)
        for pos = d.A.colptr[col]:(d.A.colptr[col+1]-1)
            I[idx] = d.A.rowval[pos]
            J[idx] = col
            idx += 1
        end
    end
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        for k in 1:length(linearidx)
            I[idx] = i+size(d.A,1)
            J[idx] = linearidx[k]
            idx += 1
        end
        for k in 1:length(quadrowidx)
            I[idx] = i+size(d.A,1)
            I[idx+1] = i+size(d.A,1)
            J[idx] = quadrowidx[k]
            J[idx+1] = quadcolidx[k]
            idx += 2
        end
    end
    return I,J
end

function eval_jac_g(d::LPQPEvaluator, J, x)
    idx = 1
    for col = 1:size(d.A,2)
        for pos = d.A.colptr[col]:(d.A.colptr[col+1]-1)
            J[idx] = d.A.nzval[pos]
            idx += 1
        end
    end
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        for k in 1:length(linearidx)
            J[idx] = linearval[k]
            idx += 1
        end
        for k in 1:length(quadrowidx)
            J[idx] = quadval[k]*x[quadcolidx[k]]
            J[idx+1] = quadval[k]*x[quadrowidx[k]]
            idx += 2
        end
    end
end

function hesslag_structure(d::LPQPEvaluator)
    hess_nnz = length(d.Qi)
    for i in 1:length(d.Qconstr)
        hess_nnz += length(d.Qconstr[i].quadrowidx)
    end

    I = Vector{Int}(undef, hess_nnz)
    J = Vector{Int}(undef, hess_nnz)

    for k in 1:length(d.Qi)
        I[k] = d.Qi[k]
        J[k] = d.Qj[k]
    end
    idx = length(d.Qi)+1
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        for k in 1:length(quadrowidx)
            qidx1 = quadrowidx[k]
            qidx2 = quadcolidx[k]
            if qidx2 > qidx1
                qidx1, qidx2 = qidx2, qidx1
            end
            I[idx] = qidx1
            J[idx] = qidx2
            idx += 1
        end
    end
    return I,J
end

function eval_hesslag(d::LPQPEvaluator, H, x, σ, μ)
    for k in 1:length(d.Qi)
        H[k] = σ*d.Qv[k]
    end
    idx = length(d.Qi) + 1
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        for k in 1:length(quadrowidx)
            l = μ[size(d.A,1)+i]
            if quadrowidx[k] == quadcolidx[k]
                H[idx] = l*2*quadval[k]
            else
                H[idx] = l*quadval[k]
            end
            idx += 1
        end
    end
end

isobjlinear(d::LPQPEvaluator) = (length(d.Qi) == 0)
isobjquadratic(d::LPQPEvaluator) = true
isconstrlinear(d::LPQPEvaluator,i::Integer) = (i <= size(d.A,1))


text:
let
    while true
        @parameters t σ ρ β
        @variables x(t) y(t) z(t)
        D = Differential(t)

        eqs = [D(D(x)) ~ σ*(y-x) + x^0.000000000000135,
               D(y) ~ x*(ρ-z)-y,
               D(z) ~ x*y - β*z]

        @named sys = ODESystem(eqs)
        sys = structural_simplify(sys)

        u0 = [D(x) => 2.0,
              x => 1.0,
              y => 0.0,
              z => 0.0]

        p  = [σ => 28.0,
              ρ => 10.0,
              β => 8/3]

        tspan = (0.0,100.0)
        prob = ODEProblem(sys,u0,tspan,p,jac=true)

        break
    end
end


text:
using Roots, NLopt, PyPlot, LsqFit, JLD, Polynomials, JFVM, Dierckx
include("DME_solver.jl")

function rec_fact(param, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME)
  (k_o, k_w) = param
  (t, R) = DME_solver_no_gravity_no_pc_incompressible(G, K, MT, T, BC, IC, k_o, k_w, false)

  PV = t.*((u_inj./poros_val)./Lx)
  R_int = Spline1D(PV, R, k=1)
  return R_int(PV_exp_DME)
end

model = (PV_exp_DME, param) -> rec_fact(param, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME)

function error_calc(param, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME, R_exp_DME)
  return rec_fact(param, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME) - R_exp_DME
end

function error_calculation(param, param_ind, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME, R_exp_DME; grad_calc=true)
  eps1 = 1e-8
  error_vals = error_calc(param, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME, R_exp_DME)
  if grad_calc
    jac  = zeros(length(error_vals), length(param_ind))
  else
    jac = zeros(0)
  end
  if grad_calc
    for j in eachindex(param_ind)
      param2 = copy(param)
      param2[param_ind[j]]+=eps1
      error_val2 = error_calc(param2, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME, R_exp_DME)
      jac[:, j] = (error_val2 - error_vals) / eps1
    end
  end
  return error_vals, jac
end

function objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all, param_ind, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME, R_exp_DME_DME; w=1.0)
  param2 = copy(param_all)
  param2[param_ind] = param
  grad_calc = false
  if length(grad)>0
    grad_calc = true
  end

  error_val, jac = error_calculation(param2, param_ind, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME, R_exp_DME, grad_calc = grad_calc)

  obj_func_value = sum(abs2, w.*error_val)

  if length(grad)>0
      grad[:] = 2.0*sum(w.*error_val.*jac, 1)
  end
  return obj_func_value
end

obj_fun = (param, grad) -> objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all, param_ind, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp_DME, R_exp_DME, w=w)

experiment = 1

if experiment == 1
    perm_val  = 1e-15 # [m^2] permeability
elseif experiment == 2
    perm_val  = 0.8e-15 # [m^2] permeability
end

poros_val = 0.256     # [-] porosity
ρ_rock    = 2700.0     # [kg/m^3]

D_core = 0.025
A_core  = π*D_core^2/4

if experiment == 1
    injection_rate = 4.5202702702702705e-10 # 0.00011/(3600*24) # m^3/day
    u_inj   = injection_rate/A_core # Darcy velocity
elseif experiment == 2
    injection_rate = 4.5202702702702705e-10 #u_inj = 0.00015/(3600*24) m^3/day
    u_inj   = injection_rate/A_core # Darcy velocity
end

c_inj1 = 0.0 # mol/m^3 DME in injection water at initial water flooding
c_inj2 = 10.0 # mol/m^3 DME in injection water at DME flooding
c_tracer1 = 0.0
c_tracer2 = 1.0

Kval = c -> 2.0
dKval_dc = c -> zeros(size(c))
a = 2.0*ρ_rock*1000.0*(1.0-poros_val) # m^2/m^3 = 2 m^2/g * 2700 kg/m^3 * 1000 g/kg * (1-phi)
D_oil = 0.5*10^(-9.0) # m^2/s
D_water = 2*10^(-9.0) # m^2/s

## Creating DME parameters
ρ_oil = c_o -> 748.0015057043045-0.0016431963703508173.*c_o+2.9364412314092913.*10.0.^(-7.0).*c_o.^2.0-4.107225978803957.*10.0.^(-11.0).*c_o.^3.0
dρ_dc_oil = c_o -> -0.0016431963703508173+2.0.*2.9364412314092913.*10.0.^(-7.0).*c_o-3.0.*4.107225978803957.*10.0.^(-11.0).*c_o.^2.0
dρ_dp_oil = c_o -> zeros(size(c_o))
ρ_water = c_w -> 983.6319389796079.*ones(size(c_w))
dρ_dp_water = c_w -> zeros(size(c_w))
μ_oil = c_o -> 0.0010131130815849523-2.528969514668829.*10.0.^(-5.0).*c_o+2.3873417406201346.*10.0.^(-11.0).*c_o.^2.0
dμ_dc_oil = c_o -> -2.528969514668829.*10.0.^(-5.0)+2*2.3873417406201346.*10.0.^(-11.0).*c_o
μ_water = c_w -> 0.00040706344808758737+1.47269029810352587.*10.0.^(-7.0).*c_w
dμ_dc_water = c_w -> 1.47269029810352587.*10.0.^(-7.0).*ones(size(c_w))

MT = FluidParametersDME(ρ_oil, dρ_dc_oil, dρ_dp_oil, ρ_water, dρ_dp_water, μ_oil, dμ_dc_oil, μ_water, dμ_dc_water, Kval, dKval_dc, a, D_oil, D_water)

# Creating realitive permeability

if experiment == 1
    #swc  = 0.392349
    #kro0 = 0.121566
    #krw0 = 0.269128
    #n_w  = 1.88088
    #n_o  = 1.91705
    #sor_min  = 0.03
    #sor_max  = 0.250158
    (sor_max, swc, kro0min, krw0, n_o, n_w) = [0.198781  0.205555  0.249794  0.149609  3.03704  3.22222]
    sor_min  = 0.01 # equivalent of equilibrium mass frac with 0.1 mol frac DME-water
    kro0max = 0.95
elseif experiment == 2
    swc  = 0.317254
    kro0 = 0.563868
    krw0 = 0.366759
    n_o  = 2.38513
    n_w  = 2.40567
    sor_min  = 0.1
    sor_max  = 0.3184591
end

m_kro = -(kro0max-kro0min)/(sor_max-sor_min)
kro0_c = c -> m_kro*c + kro0max

c_min = 0.0 # min DME mol/m^3 in oil
c_max = Kval(0)*c_inj2 # max DME mol/m^3 fraction in oil
c_max = c_max[1]
m_sor = -(sor_max-sor_min)/(c_max-c_min)
sor_c = c -> m_sor*c+sor_max
dsor_c = c -> m_sor + 0.0*c #irrelevant

K = RelativePermeabilityDME(kro0_c, krw0, n_o, n_w, sor_c, dsor_c, swc)

# =============================================================================
# Preparing geology/mesh
# =============================================================================
dim = 1

srand(123456543)
λ_dp    = 0.8
cor_x   = 0.1
cor_y   = 0.1
cor_z   = 0.1
Lx      = 0.104 # [m] length of core
Ly      = 0.025 # [m] diameter of core
Lz      = 0.025 # [m] diameter of core
Nx      = 50 # number of grids in the x direction
Ny      = 20 # number of grids in the y direction
Nz      = 10 # number of grids in the z direction

if dim == 1
    m           = createMesh1D(Nx, Lx)
    #perm_mat    = permfieldlogrnde(Nx, perm_val, λ_dp, cor_x)
elseif dim == 2
    m           = createMesh2D(Nx, Ny, Lx, Ly)
    #perm_mat    = permfieldlogrnde(Nx, Ny, perm_val, λ_dp, cor_x, cor_y)
else
    m           = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz)
    #perm_mat    = permfieldlogrnde(Nx, Ny, Nz, perm_val, λ_dp, cor_x, cor_y, cor_z)
end

k       = createCellVariable(m, perm_val)
k_face  = harmonicMean(k)     # permeability on the cell faces
ϕ       = createCellVariable(m, poros_val)

G = Geology(m, perm_val, poros_val, k, ϕ)

# =============================================================================
# Preparing time-stepping parameters
# =============================================================================

if experiment == 1
    n_pv_water1 = 9.0 # number of injected pore volumes
elseif experiment == 2
    n_pv_water1 = 13.3
end
t_final_water1  = n_pv_water1*Lx/(u_inj/poros_val) # [s] final time for water flood
dt_water1       = t_final_water1/n_pv_water1/Nx/5 # [s] time step

if experiment == 1
    n_pv_DME = 9.0
elseif experiment == 2
    n_pv_DME = 0.3
end
t_final_DME     = t_final_water1+n_pv_DME*Lx/(u_inj/poros_val) # [s] final time for DME flood
dt_DME          = t_final_DME/(n_pv_water1 + n_pv_DME)/Nx/5 # [s] time step

if experiment == 1
    n_pv_water2 = 3.8
elseif experiment == 2
    n_pv_water2 = 3.1
end

t_final_water2  = t_final_DME+n_pv_water2*Lx/(u_inj/poros_val) # [s] final time for DME flood
dt_water2       = t_final_water2/(n_pv_water1 + n_pv_DME + n_pv_water2)/Nx/5 # [s] time step

T_water1 = SimulationTime(t_final_water1, dt_water1)
T_DME    = SimulationTime(t_final_DME, dt_DME)
T_water2 = SimulationTime(t_final_water2, dt_water2)

T = [T_water1, T_DME, T_water2]

## Boundary Conditions for water flooding
#T0 = 70 + 273.15 # [K]

c_oil0      = 0.0 # initial DME concentration in the oil phase
c_water0    = 0.0 # initial DME concentration in the water phase
p0          = 2000/14.7*1e5   # [Pa]

if experiment == 1
    sw0 = 1.0-0.667 # [vol frac] initial water saturation
elseif experiment == 2
    sw0 = 1.0-0.738
end

c_t0 = 0.0 # initial tracer concentration

BCc_oil_water     = createBC(m) # concentration (DME in oil) boundary condition
BCc_water_water   = createBC(m) # concentration (DME in water) boundary condition
BCp_water         = createBC(m) # pressure boundary condition
BCs_water         = createBC(m) # saturation boundary condition
BCt_water         = createBC(m) # tracer concentration boundary condition

BCp_water.right.a[:] = 0.0
BCp_water.right.b[:] = 1.0 #Dirichlet
BCp_water.right.c[:] = p0

BCp_water.left.a[:]  = k_face.xvalue[1,:]./μ_water(c_inj1) #Neumann
BCp_water.left.b[:]  = 0.0
BCp_water.left.c[:]  = -u_inj

BCs_water.left.a[:]  = 0.0
BCs_water.left.b[:]  = 1.0 #Dirichlet
BCs_water.left.c[:]  = 1.0 #Value S_w=1

BCc_water_water.left.a[:]  = 0.0 #
BCc_water_water.left.b[:]  = 1.0 #Dirichlet
BCc_water_water.left.c[:]  = c_inj1 #injected concentration of DME in water

BCc_oil_water.left.a[:]  = 0.0 #
BCc_oil_water.left.b[:]  = 1.0 #Dirichlet
BCc_oil_water.left.c[:]  = Kval(0)*c_inj1 #injected concentration of DME in water

BCt_water.left.a[:]  = 0.0
BCt_water.left.b[:]  = 1.0 #Dirichlet
BCt_water.left.c[:]  = c_tracer1

BC_water = BoundaryConditionDME(BCc_oil_water, BCc_water_water, BCp_water, BCs_water, BCt_water)

c_oil_init      = createCellVariable(m, c_oil0, BCc_oil_water)
c_water_init    = createCellVariable(m, c_water0, BCc_water_water)
p_init          = createCellVariable(m, p0, BCp_water)
sw_init         = createCellVariable(m, sw0, BCs_water)
c_t_init        = createCellVariable(m, c_t0, BCt_water)

IC = InitialConditionDME(c_oil_init, c_water_init, p_init, sw_init, c_t_init)

# Creating Boundary conditions for DME flooding
BCc_oil_DME     = createBC(m) # concentration (DME in oil) boundary condition
BCc_water_DME   = createBC(m) # concentration (DME in water) boundary condition
BCp_DME         = createBC(m) # pressure boundary condition
BCs_DME         = createBC(m) # saturation boundary condition
BCt_DME         = createBC(m) # tracer concentration boundary condition

BCc_water_DME.left.a[:] = 0.0
BCc_water_DME.left.b[:] = 1.0 #Dirichlet
BCc_water_DME.left.c[:] = c_inj2 #Injected concentration of DME

BCc_oil_DME.left.a[:] = 0.0
BCc_oil_DME.left.b[:] = 1.0 #Dirichlet
BCc_oil_DME.left.c[:] = Kval(0)*c_inj2 #Created from reaction at boundary

BCp_DME.right.a[:] = 0.0
BCp_DME.right.b[:] = 1.0 #Dirichlet
BCp_DME.right.c[:] = p0

BCp_DME.left.a[:]  = k_face.xvalue[1,:]/μ_water(c_inj2) #Neumann
BCp_DME.left.b[:]  = 0.0
BCp_DME.left.c[:]  = -u_inj

BCs_DME.left.a[:]  = 0.0
BCs_DME.left.b[:]  = 1.0 #Dirichlet
BCs_DME.left.c[:]  = 1.0 #Value S_w=1

BCt_DME.left.a[:] = 0.0
BCt_DME.left.b[:] = 1.0 #Dirichlet
BCt_DME.left.c[:] = c_tracer2

BC_DME = BoundaryConditionDME(BCc_oil_DME, BCc_water_DME, BCp_DME, BCs_DME, BCt_DME)

BC = [BC_water, BC_DME, BC_water]

## Experimental results
if experiment == 1
    data = [0, 0.03883495145631066,
            0.18404907975460127, 0.23300970873786409,
            0.3067484662576685, 0.4271844660194175,
            0.42944785276073616, 0.4951456310679613,
            0.8588957055214725, 0.5048543689320388,
            0.9815950920245402, 0.5242718446601942,
            1.288343558282209, 0.5339805825242718,
            1.5950920245398772, 0.5436893203883496,
            2.2085889570552153, 0.5631067961165048,
            2.883435582822086, 0.5825242718446602,
            3.619631901840491, 0.5825242718446602,
            4.294478527607362, 0.5922330097087378,
            4.969325153374233, 0.5922330097087378,
            5.644171779141105, 0.6019417475728156,
            6.380368098159509, 0.6019417475728156,
            7.0552147239263805, 0.6019417475728156,
            7.730061349693252, 0.6019417475728156,
            8.404907975460121, 0.6019417475728156,
            9.079754601226991, 0.6019417475728156,
            9.570552147239262, 0.6019417475728156,
            9.938650306748466, 0.6116504854368932,
            10.061349693251532, 0.6310679611650485,
            10.184049079754601, 0.6893203883495146,
            10.306748466257668, 0.7475728155339806,
            10.429447852760735, 0.7669902912621359,
            10.736196319018404, 0.8155339805825242,
            10.920245398773005, 0.8252427184466019,
            11.16564417177914, 0.8446601941747572,
            11.411042944785276, 0.854368932038835,
            11.717791411042944, 0.8640776699029126,
            12.085889570552146, 0.8737864077669903,
            12.760736196319018, 0.9029126213592233,
            13.43558282208589, 0.9223300970873787,
            14.478527607361961, 0.941747572815534,
            15.0920245398773, 0.9514563106796117,
            15.766871165644169, 0.9514563106796117,
            16.44171779141104, 0.9514563106796117,
            17.116564417177912, 0.9514563106796117,
            17.791411042944784, 0.9514563106796117,
            18.343558282208587, 0.9514563106796117,
            18.71165644171779, 0.9514563106796117,
            19.01840490797546, 0.9514563106796117,
            19.386503067484657, 0.9514563106796117,
            20.06134969325153, 0.9514563106796117,
            20.7361963190184, 0.9514563106796117,
            21.411042944785272, 0.9514563106796117]

elseif experiment == 2
    data = [0.046242774566473965, 0.08333333333333326,
            0.18497109826589608, 0.26041666666666663,
            0.323699421965318, 0.4375,
            0.601156069364162, 0.5,
            0.8786127167630056, 0.5,
            1.0635838150289019, 0.5,
            1.2023121387283238, 0.5312499999999999,
            1.38728323699422, 0.5312499999999999,
            1.6647398843930634, 0.5312499999999999,
            1.8497109826589597, 0.5312499999999999,
            2.0346820809248554, 0.5312499999999999,
            2.6820809248554918, 0.5520833333333333,
            3.421965317919075, 0.5520833333333333,
            4.069364161849711, 0.5520833333333333,
            4.809248554913295, 0.5520833333333333,
            5.456647398843931, 0.5520833333333333,
            6.1502890173410405, 0.5520833333333333,
            6.84393063583815, 0.5520833333333333,
            7.5375722543352595, 0.5520833333333333,
            8.23121387283237, 0.5520833333333333,
            8.92485549132948, 0.5520833333333333,
            9.61849710982659, 0.5520833333333333,
            10.312138728323701, 0.5520833333333333,
            11.00578034682081, 0.5520833333333333,
            11.653179190751446, 0.5520833333333333,
            12.346820809248555, 0.5520833333333333,
            13.086705202312139, 0.5520833333333333,
            13.271676300578035, 0.5520833333333333,
            13.45664739884393, 0.5520833333333333,
            13.641618497109828, 0.5520833333333333,
            13.826589595375722, 0.5520833333333333,
            14.011560693641618, 0.5520833333333333,
            14.150289017341041, 0.5520833333333333,
            14.335260115606937, 0.5729166666666666,
            14.47398843930636, 0.6041666666666665,
            14.612716763005782, 0.625,
            14.751445086705203, 0.6354166666666665,
            14.890173410404625, 0.6458333333333333,
            15.075144508670519, 0.6458333333333333,
            15.260115606936417, 0.65625,
            15.44508670520231, 0.65625,
            15.676300578034684, 0.65625,
            15.861271676300579, 0.65625,
            16, 0.65625,
            16.739884393063583, 0.65625]
end

N_d = length(data)
PV_exp = zeros(N_d/2)
R_exp = zeros(N_d/2)
count_PV = 1
count_R = 1
for i = 1:N_d
    if mod(i,2) == 1
        PV_exp[count_PV] = data[i]
        count_PV += 1
    elseif mod(i,2) == 0
        R_exp[count_R] = data[i]
        count_R += 1
    end
end

PV_exp_DME = zeros(0)
R_exp_DME = zeros(0)
N_obs = length(PV_exp)
for i = 1:N_obs
    if PV_exp[i] >= n_pv_water1
        PV_exp_DME = push!(PV_exp_DME, PV_exp[i])
        R_exp_DME = push!(R_exp_DME, R_exp[i])
    end
end

# =============================================================================
# Calculating heuristic for k_o and k_w
# =============================================================================
d = sqrt(perm_val) #Characteristic length, i.e. diameter of capillary in porous medium
c_o = Kval(0)*c_inj2
c_o = c_o[1]
c_w = c_inj2
v_o = u_inj/poros_val
v_w = u_inj/poros_val

k_o_ini = 0.023*D_oil/d*((ρ_oil(c_o)*v_o*d)/μ_oil(c_o))^(0.83)*(μ_oil(c_o)/(ρ_oil(c_o)*D_oil))^(1.0/3.0)
k_w_ini = 0.023*D_water/d*((ρ_water(c_w)[1]*v_w*d)/μ_water(c_w))^(0.83)*(μ_water(c_w)/(ρ_water(c_w)[1]*D_water))^(1.0/3.0)

R_heu = rec_fact([k_o_ini, k_w_ini], G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp)

param_all = [1e-3*k_o_ini, 1e-3*k_w_ini]
param_ind = [1, 2]
w = ones(length(R_exp_DME))

x_init = copy(param_all)
x_lb = [1e-13, 1e-13]
x_ub = [1e-10, 1e-10]

# =============================================================================
# Optimizing
# =============================================================================

opt_type = 2

if opt_type == 1
    fit = curve_fit(model, PV_exp_DME, R_exp_DME, w, x_init, lower = x_lb, upper = x_ub)
    print("k_o_ini=", k_o_ini, ", k_w_ini=", k_w_ini, "\n\n")
    print("k_o_opt=", fit.param[1], ", k_w_opt=", fit.param[2])

elseif opt_type == 2
    # algorithms
    # :LD_MMA
    # :LN_COBYLA
    # :LD_LBFGS
    # :GN_DIRECT
    # :GN_DIRECT_L
    opt_alg=:LD_MMA

    opt1 = Opt(opt_alg, length(x_init)) # choose the algorithm
    lower_bounds!(opt1, x_lb)
    upper_bounds!(opt1, x_ub)
    ftol_rel!(opt1, 1e-5)
    ftol_abs!(opt1, 1e-5)

    min_objective!(opt1, obj_fun)
    (fObjOpt, paramOpt, flag) = optimize(opt1, x_init)
    print("k_o_ini=", k_o_ini, ", k_w_ini=", k_w_ini, "\n\n")
    print("k_o_opt=", paramOpt[1], ", k_w_opt=", paramOpt[2])
end

R_opt = rec_fact(paramOpt, G, K, MT, T, BC, IC, Lx, u_inj, poros_val, PV_exp)

#Visualizing
figure()
plot(PV_exp, R_exp.*100)
plot(PV_exp, R_heu.*100)
plot(PV_exp, R_opt.*100)
plot([n_pv_water1, n_pv_water1], [0.0, 100.0], color=:black, linestyle="--")
plot([n_pv_water1 + n_pv_DME, n_pv_water1 + n_pv_DME], [0.0, 100.0], color=:black, linestyle="--")
axis([0, n_pv_water1 + n_pv_DME + n_pv_water2, 0, 100])
grid("on")
title("Experimental vs. Simulated Recovery Factor")
xlabel(L"$PVI$")
ylabel(L"$RF$"" [%]")
legend(["Experimental", "Heuristic", "Optimized"])
savefig(string("Optimization Results/ExpVsSim", experiment, ".png"))

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers