text:
module DynamicNLPModels

import NLPModels
import QuadraticModels
import LinearAlgebra
import SparseArrays
import SparseArrays: SparseMatrixCSC

export LQDynamicData, SparseLQDynamicModel, DenseLQDynamicModel, get_u, get_s

abstract type AbstractLQDynData{T,V} end
"""
    LQDynamicData{T,V,M,MK} <: AbstractLQDynData{T,V}

A struct to represent the features of the optimization problem

```math
    minimize    \\frac{1}{2} \\sum_{i = 0}^{N-1}(s_i^T Q s_i + 2 u_i^T S^T x_i + u_i^T R u_i) + \\frac{1}{2} s_N^T Qf s_N
    subject to  s_{i+1} = A s_i + B u_i  for i=0, 1, ..., N-1
                u_i = Kx_i + v_i  \\forall i = 0, 1, ..., N - 1
                gl \\le E s_i + F u_i \\le gu for i = 0, 1, ..., N-1
                sl \\le s \\le su
                ul \\le u \\le uu
                s_0 = s0
```
---
Attributes include:
- `s0`: initial state of system
- `A` : constraint matrix for system states
- `B` : constraint matrix for system inputs
- `Q` : objective function matrix for system states from 1:(N-1)
- `R` : objective function matrix for system inputs from 1:(N-1)
- `N` : number of time steps
- `Qf`: objective function matrix for system state at time N
- `S` : objective function matrix for system states and inputs
- `ns`: number of state variables
- `nu`: number of input varaibles
- `E` : constraint matrix for state variables
- `F` : constraint matrix for input variables
- `K` : feedback gain matrix
- `sl`: vector of lower bounds on state variables
- `su`: vector of upper bounds on state variables
- `ul`: vector of lower bounds on input variables
- `uu`: vector of upper bounds on input variables
- `gl`: vector of lower bounds on constraints
- `gu`: vector of upper bounds on constraints

see also `LQDynamicData(s0, A, B, Q, R, N; ...)`
"""
struct LQDynamicData{T, V, M, MK} <: AbstractLQDynData{T,V}
    s0::V
    A::M
    B::M
    Q::M
    R::M
    N

    Qf::M
    S::M
    ns::Int
    nu::Int
    E::M
    F::M
    K::MK

    sl::V
    su::V
    ul::V
    uu::V
    gl::V
    gu::V
end

"""
    LQDynamicData(s0, A, B, Q, R, N; ...) -> LQDynamicData{T, V, M, MK}
A constructor for building an object of type `LQDynamicData` for the optimization problem
```math
    minimize    \\frac{1}{2} \\sum_{i = 0}^{N-1}(s_i^T Q s_i + 2 u_i^T S^T x_i + u_i^T R u_i) + \\frac{1}{2} s_N^T Qf s_N
    subject to  s_{i+1} = A s_i + B u_i  \\forall i=0, 1, ..., N - 1
                u_i = Kx_i + v_i  \\forall i = 0, 1, ..., N - 1
                gl \\le E s_i + F u_i \\le gu \\forall i = 0, 1, ..., N-1
                sl \\le s \\le su
                ul \\le u \\le uu
                s_0 = s0
```
---
- `s0`: initial state of system
- `A` : constraint matrix for system states
- `B` : constraint matrix for system inputs
- `Q` : objective function matrix for system states from 1:(N-1)
- `R` : objective function matrix for system inputs from 1:(N-1)
- `N` : number of time steps
The following attributes of the `LQDynamicData` type are detected automatically from the length of s0 and size of R
- `ns`: number of state variables
- `nu`: number of input varaibles
The following keyward arguments are also accepted
- `Qf = Q`: objective function matrix for system state at time N; dimensions must be ns x ns
- `S  = nothing`: objective function matrix for system state and inputs
- `E  = zeros(eltype(Q), 0, ns)`  : constraint matrix for state variables
- `F  = zeros(eltype(Q), 0, nu)`  : constraint matrix for input variables
- `K  = nothing`       : feedback gain matrix
- `sl = fill(-Inf, ns)`: vector of lower bounds on state variables
- `su = fill(Inf, ns)` : vector of upper bounds on state variables
- `ul = fill(-Inf, nu)`: vector of lower bounds on input variables
- `uu = fill(Inf, nu)` : vector of upper bounds on input variables
- `gl = fill(-Inf, size(E, 1))` : vector of lower bounds on constraints
- `gu = fill(Inf, size(E, 1))`  : vector of upper bounds on constraints
"""
function LQDynamicData(
    s0::V,
    A::M,
    B::M,
    Q::M,
    R::M,
    N;

    Qf::M = Q,
    S::M  = (similar(Q, size(Q, 1), size(R, 1)) .= 0),
    E::M  = (similar(Q, 0, length(s0)) .= 0),
    F::M  = (similar(Q, 0, size(R, 1)) .= 0),
    K::MK = nothing,

    sl::V = (similar(s0) .= -Inf),
    su::V = (similar(s0) .=  Inf),
    ul::V = (similar(s0, size(R, 1)) .= -Inf),
    uu::V = (similar(s0, size(R, 1)) .=  Inf),
    gl::V = (similar(s0, size(E, 1)) .= -Inf),
    gu::V = (similar(s0, size(F, 1)) .= Inf)
    ) where {T,V <: AbstractVector{T}, M <: AbstractMatrix{T}, MK <: Union{Nothing, AbstractMatrix{T}}}

    if size(Q, 1) != size(Q, 2)
        error("Q matrix is not square")
    end
    if size(R, 1) != size(R, 1)
        error("R matrix is not square")
    end
    if size(A, 2) != length(s0)
        error("Number of columns of A are not equal to the number of states")
    end
    if size(B, 2) != size(R, 1)
        error("Number of columns of B are not equal to the number of inputs")
    end
    if length(s0) != size(Q, 1)
        error("size of Q is not consistent with length of s0")
    end

    if !(sl  <= su)
        error("lower bound(s) on x is > upper bound(s)")
    end
    if !(ul <= uu)
        error("lower bound(s) on u is > upper bound(s)")
    end
    if !(s0 >= sl) || !(s0 <= su)
        error("s0 is not within the given upper and lower bounds")
    end

    if size(E, 1) != size(F, 1)
        error("E and F have different numbers of rows")
    end
    if !(gl <= gu)
        error("lower bound(s) on Es + Fu is > upper bound(s)")
    end
    if size(E, 2) != size(Q, 1)
        error("Dimensions of E are not the same as number of states")
    end
    if size(F, 2) != size(R, 1)
        error("Dimensions of F are not the same as the number of inputs")
    end
    if length(gl) != size(E, 1)
        error("Dimensions of gl do not match E and F")
    end
    if length(gu) != size(E, 1)
        error("Dimensions of gu do not match E and F")
    end
    if size(S, 1) != size(Q, 1) || size(S, 2) != size(R, 1)
        error("Dimensions of S do not match dimensions of Q and R")
    end
    if K != nothing
        if size(K, 1) != size(R, 1) || size(K, 2) != size(Q,1)
            error("Dimensions of K  do not match number of states and inputs")
        end
    end

    ns = size(Q,1)
    nu = size(R,1)

    LQDynamicData{T, V, M, MK}(
        s0, A, B, Q, R, N,
        Qf, S, ns, nu, E, F, K,
        sl, su, ul, uu, gl, gu
    )
end

abstract type AbstractDynamicModel{T,V} <: QuadraticModels.AbstractQuadraticModel{T, V} end

mutable struct SparseLQDynamicModel{T, V, M1, M2, M3, MK} <:  AbstractDynamicModel{T,V}
  meta::NLPModels.NLPModelMeta{T, V}
  counters::NLPModels.Counters
  data::QuadraticModels.QPData{T, V, M1, M2}
  dynamic_data::LQDynamicData{T, V, M3, MK}
end

struct Block_Matrices{T, M1, M2}
    A::M1
    B::M1
    Q::M2
    R::M2
    S::M2
    K::M2
    E::M1
    F::M1
    gl::M1
    gu::M1
end

function Block_Matrices(
    A::M1,
    B::M1,
    Q::M2,
    R::M2,
    S::M2,
    K::M2,
    E::M1,
    F::M1,
    gl::M1,
    gu::M1
) where {T, M1 <: AbstractMatrix{T}, M2 <: AbstractMatrix{T}}

    Block_Matrices{T, M1, M2}(
        A,
        B,
        Q,
        R,
        S,
        K,
        E,
        F,
        gl,
        gu
    )
end

mutable struct DenseLQDynamicModel{T, V, M1, M2, M3, M4, MK} <:  AbstractDynamicModel{T,V}
    meta::NLPModels.NLPModelMeta{T, V}
    counters::NLPModels.Counters
    data::QuadraticModels.QPData{T, V, M1, M2}
    dynamic_data::LQDynamicData{T, V, M3, MK}
    blocks::Block_Matrices{T, M3, M4}
end

"""
    SparseLQDynamicModel(dnlp::LQDynamicData)    -> SparseLQDynamicModel
    SparseLQDynamicModel(s0, A, B, Q, R, N; ...) -> SparseLQDynamicModel
A constructor for building a `SparseLQDynamicModel <: QuadraticModels.AbstractQuadraticModel`
Input data is for the problem of the form
```math
    minimize    \\frac{1}{2} \\sum_{i = 0}^{N-1}(s_i^T Q s_i + 2 u_i^T S^T x_i + u_i^T R u_i) + \\frac{1}{2} s_N^T Qf s_N
    subject to  s_{i+1} = A s_i + B u_i  for i=0, 1, ..., N-1
                u_i = Kx_i + v_i  \\forall i = 0, 1, ..., N - 1
                gl \\le E s_i + F u_i \\le gu for i = 0, 1, ..., N-1
                sl \\le s \\le su
                ul \\le u \\le uu
                s_0 = s0
```
---

Data is converted to the form

```math
    minimize    \\frac{1}{2} z^T H z
    subject to  lcon \\le Jz \\le ucon
                lvar \\le z \\le uvar
```
Resulting `H` and `J` matrices are stored as `QuadraticModels.QPData` within the `SparseLQDynamicModel` struct and
variable and constraint limits are stored within `NLPModels.NLPModelMeta`

If `K` is defined, then `u` variables are replaced by `v` variables, and `u` can be queried by `get_u` and `get_s` within `DynamicNLPModels.jl`
"""
function SparseLQDynamicModel(dnlp::LQDynamicData{T,V,M}) where {T, V <: AbstractVector{T}, M  <: AbstractMatrix{T}, MK <: Union{Nothing, AbstractMatrix{T}}}
        _build_sparse_lq_dynamic_model(dnlp)
end

function SparseLQDynamicModel(
    s0::V,
    A::M,
    B::M,
    Q::M,
    R::M,
    N;
    Qf::M = Q,
    S::M  = (similar(Q, size(Q, 1), size(R, 1)) .= 0),
    E::M  = (similar(Q, 0, length(s0)) .= 0),
    F::M  = (similar(Q, 0, size(R, 1)) .= 0),
    K::MK = nothing,
    sl::V = (similar(s0) .= -Inf),
    su::V = (similar(s0) .=  Inf),
    ul::V = (similar(s0, size(R, 1)) .= -Inf),
    uu::V = (similar(s0, size(R, 1)) .=  Inf),
    gl::V = (similar(s0, size(E, 1)) .= -Inf),
    gu::V = (similar(s0, size(F, 1)) .= Inf),
    dense=false
) where {T, V <: AbstractVector{T}, M <: AbstractMatrix{T}, MK <: Union{Nothing, AbstractMatrix{T}}}

    dnlp = LQDynamicData(
        s0, A, B, Q, R, N;
        Qf = Qf, S = S, E = E, F = F, K = K,
        sl = sl, su = su, ul = ul, uu = uu, gl = gl, gu = gu)

    SparseLQDynamicModel(dnlp)
end

"""
    DenseLQDynamicModel(dnlp::LQDynamicData)    -> DenseLQDynamicModel
    DenseLQDynamicModel(s0, A, B, Q, R, N; ...) -> DenseLQDynamicModel
A constructor for building a `DenseLQDynamicModel <: QuadraticModels.AbstractQuadraticModel`

Input data is for the problem of the form
```math
    minimize    \\frac{1}{2} \\sum_{i = 0}^{N-1}(s_i^T Q s_i + 2 u_i^T S^T x_i + u_i^T R u_i) + \\frac{1}{2} s_N^T Qf s_N
    subject to  s_{i+1} = A s_i + B u_i  for i=0, 1, ..., N-1
                u_i = Kx_i + v_i  \\forall i = 0, 1, ..., N - 1
                gl \\le E s_i + F u_i \\le gu for i = 0, 1, ..., N-1
                sl \\le s \\le su
                ul \\le u \\le uu
                s_0 = s0
```
---

Data is converted to the form

```math
    minimize    \\frac{1}{2} u^T H u + h^T u + h0
    subject to  Jz \\le g
                ul \\le u \\le uu
```

Resulting `H`, `J`, `h`, and `h0` matrices are stored within `QuadraticModels.QPData` as `H`, `A`, `c`, and `c0` attributes respectively

If `K` is defined, then `u` variables are replaced by `v` variables. The bounds on `u` are transformed into algebraic constraints,
and `u` can be queried by `get_u` and `get_s` within `DynamicNLPModels.jl`
"""
function DenseLQDynamicModel(dnlp::LQDynamicData{T,V,M}) where {T, V <: AbstractVector{T}, M  <: AbstractMatrix{T}, MK <: Union{Nothing, AbstractMatrix{T}}}
    _build_dense_lq_dynamic_model(dnlp)
end

function DenseLQDynamicModel(
    s0::V,
    A::M,
    B::M,
    Q::M,
    R::M,
    N;
    Qf::M = Q,
    S::M  = (similar(Q, size(Q, 1), size(R, 1)) .= 0),
    E::M  = (similar(Q, 0, length(s0)) .= 0),
    F::M  = (similar(Q, 0, size(R, 1)) .= 0),
    K::MK = nothing,
    sl::V = (similar(s0) .= -Inf),
    su::V = (similar(s0) .=  Inf),
    ul::V = (similar(s0, size(R, 1)) .= -Inf),
    uu::V = (similar(s0, size(R, 1)) .=  Inf),
    gl::V = (similar(s0, size(E, 1)) .= -Inf),
    gu::V = (similar(s0, size(F, 1)) .= Inf),
    dense=false
) where {T, V <: AbstractVector{T}, M <: AbstractMatrix{T}, MK <: Union{Nothing, AbstractMatrix{T}}}

    dnlp = LQDynamicData(
        s0, A, B, Q, R, N;
        Qf = Qf, S = S, E = E, F = F, K = K,
        sl = sl, su = su, ul = ul, uu = uu, gl = gl, gu = gu)

    DenseLQDynamicModel(dnlp)
end

function _build_sparse_lq_dynamic_model(dnlp::LQDynamicData{T, V, M, MK}) where {T, V <: AbstractVector{T}, M  <: AbstractMatrix{T}, MK <: Nothing}
    s0 = dnlp.s0
    A  = dnlp.A
    B  = dnlp.B
    Q  = dnlp.Q
    R  = dnlp.R
    N  = dnlp.N

    Qf = dnlp.Qf
    S  = dnlp.S
    ns = dnlp.ns
    nu = dnlp.nu
    E  = dnlp.E
    F  = dnlp.F

    sl = dnlp.sl
    su = dnlp.su
    ul = dnlp.ul
    uu = dnlp.uu
    gl = dnlp.gl
    gu = dnlp.gu

    H   = _build_H(Q, R, N; Qf = Qf, S = S)
    J1  = _build_sparse_J1(A, B, N)
    J2  = _build_sparse_J2(E, F, N)

    J = vcat(J1, J2)

    c0  = zero(eltype(s0))

    nvar = ns * (N + 1) + nu * N
    c  = similar(s0, nvar); fill!(c, 0)

    lvar  = similar(s0, nvar); fill!(lvar, 0)
    uvar  = similar(s0, nvar); fill!(uvar, 0)

    lvar[1:ns] .= s0
    uvar[1:ns] .= s0

    lcon  = similar(s0, ns * N + N * length(gl)); fill!(lcon, 0)
    ucon  = similar(s0, ns * N + N * length(gl)); fill!(ucon, 0)

    ncon  = size(J, 1)
    nnzj = length(J.rowval)
    nnzh = length(H.rowval)

    for i in 1:N
        lvar[(i * ns + 1):((i + 1) * ns)] .= sl
        uvar[(i * ns + 1):((i + 1) * ns)] .= su

        lcon[(ns * N + 1 + (i -1) * length(gl)):(ns * N + i * length(gl))] .= gl
        ucon[(ns * N + 1 + (i -1) * length(gl)):(ns * N + i * length(gl))] .= gu
    end

    for j in 1:N
        lvar[((N + 1) * ns + (j - 1) * nu + 1):((N + 1) * ns + j * nu)] .= ul
        uvar[((N + 1) * ns + (j - 1) * nu + 1):((N + 1) * ns + j * nu)] .= uu
    end

    SparseLQDynamicModel(
        NLPModels.NLPModelMeta(
        nvar,
        x0   = (similar(s0, nvar) .= 0),
        lvar = lvar,
        uvar = uvar,
        ncon = ncon,
        lcon = lcon,
        ucon = ucon,
        nnzj = nnzj,
        nnzh = nnzh,
        lin = 1:ncon,
        islp = (ncon == 0);
        ),
        NLPModels.Counters(),
        QuadraticModels.QPData(
        c0,
        c,
        H,
        J
        ),
        dnlp
    )

end

function _build_sparse_lq_dynamic_model(dnlp::LQDynamicData{T, V, M, MK}) where {T, V <: AbstractVector{T}, M  <: AbstractMatrix{T}, MK <: AbstractMatrix}
    s0 = dnlp.s0
    A  = dnlp.A
    B  = dnlp.B
    Q  = dnlp.Q
    R  = dnlp.R
    N  = dnlp.N

    Qf = dnlp.Qf
    S  = dnlp.S
    ns = dnlp.ns
    nu = dnlp.nu
    E  = dnlp.E
    F  = dnlp.F
    K  = dnlp.K

    sl = dnlp.sl
    su = dnlp.su
    ul = dnlp.ul
    uu = dnlp.uu
    gl = dnlp.gl
    gu = dnlp.gu

    # Transform u variables to v variables
    new_Q = similar(Q)
    new_S = similar(S)
    new_A = similar(A)
    new_E = similar(E)

    LinearAlgebra.copyto!(new_Q, Q)
    LinearAlgebra.copyto!(new_S, S)
    LinearAlgebra.copyto!(new_A, A)
    LinearAlgebra.copyto!(new_E, E)

    KTR  = similar(Q, size(K, 2), size(R, 2))
    LinearAlgebra.mul!(KTR, K', R)
    LinearAlgebra.axpy!(1, KTR, new_S)

    SK   = similar(Q, size(S, 1), size(K, 2))
    KTRK = similar(Q, size(K, 2), size(K, 2))
    LinearAlgebra.mul!(SK, S, K)
    LinearAlgebra.mul!(KTRK, KTR, K)
    LinearAlgebra.axpy!(1, SK, new_Q)
    LinearAlgebra.axpy!(1, SK', new_Q)
    LinearAlgebra.axpy!(1, KTRK, new_Q)

    BK    = similar(Q, size(B, 1), size(K, 2))
    LinearAlgebra.mul!(BK, B, K)
    LinearAlgebra.axpy!(1, BK, new_A)

    FK    = similar(Q, size(F, 1), size(K, 2))
    LinearAlgebra.mul!(FK, F, K)
    LinearAlgebra.axpy!(1, FK, new_E)

    # Get H and J matrices from new matrices
    H   = _build_H(new_Q, R, N; Qf = Qf, S = new_S)
    J1  = _build_sparse_J1(new_A, B, N)
    J2  = _build_sparse_J2(new_E, F, N)
    J3, lcon3, ucon3  = _build_sparse_J3(K, N, uu, ul)

    J = vcat(J1, J2)
    J = vcat(J, J3)

    nvar = ns * (N + 1) + nu * N

    lvar  = similar(s0, nvar); fill!(lvar, -Inf)
    uvar  = similar(s0, nvar); fill!(uvar, Inf)

    lvar[1:ns] .= s0
    uvar[1:ns] .= s0

    lcon  = similar(s0, ns * N + N * length(gl) + length(lcon3)); fill!(lcon, 0)
    ucon  = similar(s0, ns * N + N * length(gl) + length(lcon3)); fill!(ucon, 0)

    ncon  = size(J, 1)
    nnzj = length(J.rowval)
    nnzh = length(H.rowval)

    for i in 1:N
        lvar[(i * ns + 1):((i + 1) * ns)] .= sl
        uvar[(i * ns + 1):((i + 1) * ns)] .= su

        lcon[(ns * N + 1 + (i -1) * length(gl)):(ns * N + i * length(gl))] .= gl
        ucon[(ns * N + 1 + (i -1) * length(gl)):(ns * N + i * length(gl))] .= gu
    end

    if length(lcon3) > 0
        lcon[(1 + ns * N + N * length(gl)):end] .= lcon3
        ucon[(1 + ns * N + N * length(gl)):end] .= ucon3
    end

    c0 = zero(eltype(s0))
    c  = similar(s0, nvar); fill!(c, 0)

    SparseLQDynamicModel(
        NLPModels.NLPModelMeta(
        nvar,
        x0   = (similar(s0, nvar) .= 0),
        lvar = lvar,
        uvar = uvar,
        ncon = ncon,
        lcon = lcon,
        ucon = ucon,
        nnzj = nnzj,
        nnzh = nnzh,
        lin = 1:ncon,
        islp = (ncon == 0);
        ),
        NLPModels.Counters(),
        QuadraticModels.QPData(
        c0,
        c,
        H,
        J
        ),
        dnlp
    )
end

function _build_dense_lq_dynamic_model(dnlp::LQDynamicData{T,V,M,MK}) where {T, V <: AbstractVector{T}, M <: AbstractMatrix{T}, MK <: Nothing}
    s0 = dnlp.s0
    A  = dnlp.A
    B  = dnlp.B
    Q  = dnlp.Q
    R  = dnlp.R
    N  = dnlp.N

    Qf = dnlp.Qf
    S  = dnlp.S
    ns = dnlp.ns
    nu = dnlp.nu
    E  = dnlp.E
    F  = dnlp.F
    K  = dnlp.K

    sl = dnlp.sl
    su = dnlp.su
    ul = dnlp.ul
    uu = dnlp.uu
    gl = dnlp.gl
    gu = dnlp.gu

    dense_blocks = _build_block_matrices(s0, Q, R, A, B, E, F, N, gu, gl, K; Qf = Qf, S = S)
    block_A  = dense_blocks.A
    block_B  = dense_blocks.B
    block_Q  = dense_blocks.Q
    block_R  = dense_blocks.R
    block_E  = dense_blocks.E
    block_F  = dense_blocks.F
    block_S  = dense_blocks.S
    block_K  = dense_blocks.K
    block_gl = dense_blocks.gl
    block_gu = dense_blocks.gu

    H_blocks = _build_H_blocks(block_Q, block_R, block_A, block_B, block_S, block_K, s0, N, K)

    H  = H_blocks.H
    c0 = H_blocks.c0

    G_blocks = _build_G_blocks(block_A, block_B, block_E, block_F, block_K, block_gl, block_gu, s0, N)

    J1    = G_blocks.J
    lcon1 = G_blocks.lcon
    ucon1 = G_blocks.ucon
    As0   = G_blocks.As0

    lvar = similar(s0, nu * N); fill!(lvar, -Inf)
    uvar = similar(s0, nu * N); fill!(uvar, Inf)

    for i in 1:(N)
        lvar[((i - 1) * nu + 1):(i * nu)] = ul
        uvar[((i - 1) * nu + 1):(i * nu)] = uu
    end

    # Convert state variable constraints to algebraic constraints
    bool_vec        = (su .!= Inf .|| sl .!= -Inf)
    num_real_bounds = sum(bool_vec)

    J2         = similar(Q, num_real_bounds * N, nu * N); fill!(J2, 0)
    As0_bounds = similar(s0, num_real_bounds * N); fill!(As0_bounds, 0)

    if num_real_bounds == length(sl)
        J2        .= block_B[(1 + ns):end,:]
        As0_bounds .= As0[(1 + ns):end,1]
    else
        for i in 1:N
            row_range = (1 + (i - 1) * num_real_bounds):(i * num_real_bounds)
            J2[row_range, :] .= block_B[(1 + ns * i): (ns * (i + 1)), :][bool_vec, :]
            As0_bounds[row_range] .= As0[(1 + ns * i):(ns * (i + 1))][bool_vec]
        end

        sl = sl[bool_vec]
        su = su[bool_vec]
    end

    lcon2 = similar(s0, size(J2, 1)); fill!(lcon2, 0)
    ucon2 = similar(s0, size(J2, 1)); fill!(ucon2, 0)

    for i in 1:N
        lcon2[((i - 1) * length(su) + 1):(i * length(su)),1] = sl
        ucon2[((i - 1) * length(su) + 1):(i * length(su)),1] = su
    end

    LinearAlgebra.axpy!(-1, As0_bounds, ucon2)
    LinearAlgebra.axpy!(-1, As0_bounds, lcon2)

    lcon = similar(s0, length(lcon1) + length(lcon2)); fill!(lcon, 0)
    ucon = similar(s0, length(ucon1) + length(ucon2)); fill!(ucon, 0)

    lcon[1:length(lcon1)] .= lcon1
    ucon[1:length(ucon1)] .= ucon1

    if length(lcon2) > 0
        lcon[(1 + length(lcon1)):end] .= lcon2
        ucon[(1 + length(ucon1)):end] .= ucon2
    end

    J   = similar(Q, size(J1, 1) + size(J2, 1), size(J1,2)); fill!(J, 0)
    J[1:(size(J1, 1)), :]   .= J1
    if size(J2, 1) > 0
        J[(1 + size(J1, 1)):end, :] .= J2
    end

    nvar = nu * N
    nnzj = size(J, 1) * size(J, 2)
    nnzh = sum(LinearAlgebra.LowerTriangular(H) .!= 0)
    ncon = size(J, 1)

    c = similar(s0, nvar)
    c .= H_blocks.c

    DenseLQDynamicModel(
        NLPModels.NLPModelMeta(
        nvar,
        x0   = (similar(s0, nvar) .= 0),
        lvar = lvar,
        uvar = uvar,
        ncon = ncon,
        lcon = lcon,
        ucon = ucon,
        nnzj = nnzj,
        nnzh = nnzh,
        lin = 1:ncon,
        islp = (ncon == 0);
        ),
        NLPModels.Counters(),
        QuadraticModels.QPData(
        c0[1,1],
        c,
        H,
        J
        ),
        dnlp,
        dense_blocks
    )
end

function _build_dense_lq_dynamic_model(dnlp::LQDynamicData{T,V,M,MK}) where {T, V <: AbstractVector{T}, M <: AbstractMatrix{T}, MK <: AbstractMatrix{T}}
    s0 = dnlp.s0
    A  = dnlp.A
    B  = dnlp.B
    Q  = dnlp.Q
    R  = dnlp.R
    N  = dnlp.N

    Qf = dnlp.Qf
    S  = dnlp.S
    ns = dnlp.ns
    nu = dnlp.nu
    E  = dnlp.E
    F  = dnlp.F
    K  = dnlp.K

    sl = dnlp.sl
    su = dnlp.su
    ul = dnlp.ul
    uu = dnlp.uu
    gl = dnlp.gl
    gu = dnlp.gu

    dense_blocks = _build_block_matrices(s0, Q, R, A, B, E, F, N, gu, gl, K; Qf = Qf, S = S)
    block_A  = dense_blocks.A
    block_B  = dense_blocks.B
    block_Q  = dense_blocks.Q
    block_R  = dense_blocks.R
    block_E  = dense_blocks.E
    block_F  = dense_blocks.F
    block_S  = dense_blocks.S
    block_K  = dense_blocks.K
    block_gl = dense_blocks.gl
    block_gu = dense_blocks.gu

    H_blocks = _build_H_blocks(block_Q, block_R, block_A, block_B, block_S, block_K, s0, N, K)

    H  = H_blocks.H
    c0 = H_blocks.c0

    G_blocks = _build_G_blocks(block_A, block_B, block_E, block_F, block_K, block_gl, block_gu, s0, N)

    J1    = G_blocks.J
    lcon1 = G_blocks.lcon
    ucon1 = G_blocks.ucon
    As0   = G_blocks.As0

    # Convert state variable constraints to algebraic constraints
    bool_vec_s        = (su .!= Inf .|| sl .!= -Inf)
    num_real_bounds   = sum(bool_vec_s)

    J2         = similar(Q, num_real_bounds * N, nu * N); fill!(J2, 0)
    As0_bounds = similar(s0, num_real_bounds * N); fill!(As0_bounds, 0)

    if num_real_bounds == length(sl)
        J2         .= block_B[(1 + ns):end,:]
        As0_bounds .= As0[(1 + ns):end,1]
    else
        for i in 1:N
            row_range = (1 + (i - 1) * num_real_bounds):(i * num_real_bounds)
            J2[row_range, :] .= block_B[(1 + ns * i): (ns * (i + 1)), :][bool_vec_s, :]
            As0_bounds[row_range] .= As0[(1 + ns * i):(ns * (i + 1))][bool_vec_s]
        end

        sl = sl[bool_vec_s]
        su = su[bool_vec_s]
    end

    # Convert bounds on u to algebraic constraints
    bool_vec_u       = (ul .!= -Inf .|| uu .!= Inf)
    num_real_bounds = sum(bool_vec_u)

    KBI  = similar(Q, nu * N, nu * N)
    KAs0 = similar(s0, nu * N)

    I_mat = similar(Q, nu * N, nu * N); fill!(I_mat, 0)

    for i in 1:(nu * N)
        I_mat[i, i] = 1.0
    end

    LinearAlgebra.mul!(KBI, block_K, block_B)
    LinearAlgebra.axpy!(1, I_mat, KBI)
    LinearAlgebra.mul!(KAs0, block_K, As0)

    J3          = similar(Q, num_real_bounds * N, nu * N); fill!(J3, 0)
    KAs0_bounds = similar(s0, num_real_bounds * N); fill!(KAs0_bounds, 0)

    if num_real_bounds == length(ul)
        J3 .= KBI
        KAs0_bounds .= KAs0
    else
        for i in 1:N
            row_range   = (1 + (i - 1) * num_real_bounds):(i * num_real_bounds)
            J3[row_range, :] .= KBI[(1 + nu * (i - 1)):(nu * i), :][bool_vec_u, :]
            KAs0_bounds[row_range]      .= KAs0[(1 + nu * (i - 1)):(nu * i)][bool_vec_u]
        end

        ul = ul[bool_vec_u]
        uu = uu[bool_vec_u]
    end

    lcon2 = similar(s0, size(J2, 1)); fill!(lcon2, 0)
    ucon2 = similar(s0, size(J2, 1)); fill!(ucon2, 0)

    lcon3 = similar(s0, size(J3, 1)); fill!(lcon3, 0)
    ucon3 = similar(s0, size(J3, 1)); fill!(ucon3, 0)

    for i in 1:N
        lcon2[((i - 1) * length(su) + 1):(i * length(su)),1] = sl
        ucon2[((i - 1) * length(su) + 1):(i * length(su)),1] = su

        lcon3[((i - 1) * length(uu) + 1):(i * length(uu)),1] = ul
        ucon3[((i - 1) * length(uu) + 1):(i * length(uu)),1] = uu
    end

    LinearAlgebra.axpy!(-1, As0_bounds, lcon2)
    LinearAlgebra.axpy!(-1, As0_bounds, ucon2)

    LinearAlgebra.axpy!(-1, KAs0_bounds, lcon3)
    LinearAlgebra.axpy!(-1, KAs0_bounds, ucon3)

    J = similar(Q, (size(J1, 1) + size(J2, 1) + size(J3, 1)), size(J1, 2)); fill!(J, 0)
    J[1:size(J1, 1), :] .= J1

    if size(J2, 1) > 0
        J[(size(J1, 1) + 1):(size(J1, 1) + size(J2, 1)), :] .= J2
    end

    if size(J3, 1) >0
        J[(size(J1, 1) + size(J2, 1) + 1):(size(J1, 1) + size(J2, 1) + size(J3, 1)), :] .= J3
    end

    lcon = similar(s0, size(J, 1)); fill!(lcon, 0)
    ucon = similar(s0, size(J, 1)); fill!(ucon, 0)

    lcon[1:length(lcon1)] .= lcon1
    ucon[1:length(ucon1)] .= ucon1

    if length(lcon2) > 0
        lcon[(length(lcon1) + 1):(length(lcon1) + length(lcon2))] .= lcon2
        ucon[(length(ucon1) + 1):(length(ucon1) + length(ucon2))] .= ucon2
    end

    if length(lcon3) > 0
        lcon[(length(lcon1) + length(lcon2) + 1):(length(lcon1) + length(lcon2) + length(lcon3))] .= lcon3
        ucon[(length(ucon1) + length(ucon2) + 1):(length(ucon1) + length(ucon2) + length(ucon3))] .= ucon3
    end

    nvar = nu * N
    nnzj = size(J, 1) * size(J, 2)
    nnzh = sum(LinearAlgebra.LowerTriangular(H) .!= 0)
    ncon = size(J, 1)

    c = similar(s0, nvar)
    c .= H_blocks.c

    DenseLQDynamicModel(
        NLPModels.NLPModelMeta(
        nvar,
        x0   = (similar(s0, nvar) .= 0),
        ncon = ncon,
        lcon = lcon,
        ucon = ucon,
        nnzj = nnzj,
        nnzh = nnzh,
        lin = 1:ncon,
        islp = (ncon == 0);
        ),
        NLPModels.Counters(),
        QuadraticModels.QPData(
        c0[1,1],
        c,
        H,
        J
        ),
        dnlp,
        dense_blocks
    )
end

function _build_block_matrices(
    s0, Q, R, A, B, E, F, N, gu, gl, K;
    Qf = Q,
    S = (similar(Q, size(Q, 1), size(R, 1)) .= 0)
    )

    ns = size(Q, 1)
    nu = size(R, 1)

    if K == nothing
        K = similar(Q, nu, ns); fill!(K, 0)
    end

    # Define block matrices
    block_B = similar(Q, ns * (N + 1), nu * N); fill!(block_B, 0)
    block_A = similar(Q, ns * (N + 1), ns); fill!(block_A, 0)
    block_Q = SparseArrays.sparse([],[], eltype(Q)[], ns * (N + 1), ns * (N + 1))
    block_R = SparseArrays.sparse([],[], eltype(R)[], nu * N, nu * N)
    block_S = SparseArrays.sparse([],[], eltype(S)[], ns * (N + 1), nu * N)
    block_K = SparseArrays.sparse([],[], eltype(R)[], nu * N, ns * (N + 1))

    nE1 = size(E, 1)
    nE2 = size(E, 2)
    nF1 = size(F, 1)
    nF2 = size(F, 2)

    block_E  = similar(Q, nE1 * N, nE2 * (N + 1)); fill!(block_E, 0)
    block_F  = similar(Q, nF1 * N, nF2 * N); fill!(block_F, 0)
    block_gl = similar(Q, nE1 * N, 1); fill!(block_gl, 0)
    block_gu = similar(Q, nE1 * N, 1); fill!(block_gu, 0)

    # Build E, F, and d (gl and gu) blocks
    for i in 1:N
        block_E[((i - 1) * nE1 + 1):(i * nE1), ((i - 1) * nE2 + 1):(i * nE2)] = E
        block_F[((i - 1) * nF1 + 1):(i * nF1), ((i - 1) * nF2 + 1):(i * nF2)] = F
        block_gl[((i - 1) * nE1  + 1):(i * nE1)]  = gl
        block_gu[((i - 1) * nE1  + 1):(i * nE1)]  = gu
    end

    # Add diagonal of Bs and fill Q, R, S, and K block matrices
    for j in 1:N
        B_row_range = (j * ns + 1):((j + 1) * ns)
        B_col_range = ((j - 1) * nu + 1):(j * nu)
        block_B[B_row_range, B_col_range] = B

        block_Q[((j - 1) * ns + 1):(j * ns), ((j - 1) * ns + 1):(j * ns)] = Q
        block_R[((j - 1) * nu + 1):(j * nu), ((j - 1) * nu + 1):(j * nu)] = R

        S_row_range = (1 + ns * (j - 1)):(ns * j)
        S_col_range = (1 + nu * (j - 1)):(nu * j)

        block_S[S_row_range, S_col_range] = S

        K_row_range = ((j - 1) * nu + 1):(j * nu)
        block_K[K_row_range, ((j - 1) * ns + 1):(j * ns)] = K
    end

    for i in 1:ns
        block_A[i, i] = 1.0
    end

    A_k = copy(A)
    BK  = similar(Q, size(B, 1), size(K, 2))
    LinearAlgebra.mul!(BK, B, K)
    LinearAlgebra.axpy!(1, BK, A_k)

    # Define matrices for mul!
    A_klast  = copy(A_k)
    A_knext  = copy(A_k)
    AB_klast = similar(Q, size(B, 1), size(B, 2)); fill!(AB_klast, 0)
    AB_k     = similar(Q, size(B, 1), size(B, 2)); fill!(AB_k, 0)

    # Fill the A and B matrices
    for i in 1:(N - 1)
        if i == 1
            block_A[(ns + 1):ns*2, :] = A_k
            LinearAlgebra.mul!(AB_k, A_k, B)
            for k in 1:(N-i)
                row_range = (1 + (k + 1) * ns):((k + 2) * ns)
                col_range = (1 + (k - 1) * nu):(k * nu)
                block_B[row_range, col_range] = AB_k
            end
            AB_klast = copy(AB_k)
        else
            LinearAlgebra.mul!(AB_k, A_k, AB_klast)
            LinearAlgebra.mul!(A_knext, A_k, A_klast)
            block_A[(ns * i + 1):ns * (i + 1),:] = A_knext

            for k in 1:(N-i)
                row_range = (1 + (k + i) * ns):((k + i + 1) * ns)
                col_range = (1 + (k - 1) * nu):(k * nu)
                block_B[row_range, col_range] = AB_k
            end

            AB_klast = copy(AB_k)
            A_klast  = copy(A_knext)
        end
    end

    LinearAlgebra.mul!(A_knext, A_k, A_klast)

    block_A[(ns * N + 1):ns * (N + 1), :] = A_knext
    block_Q[(ns * N + 1):((N + 1) * ns), (N * ns + 1):((N + 1) * ns)] = Qf

    Block_Matrices(
        block_A,
        block_B,
        block_Q,
        block_R,
        block_S,
        block_K,
        block_E,
        block_F,
        block_gl,
        block_gu
    )
end

function _build_H_blocks(block_Q, block_R, block_A::M, block_B::M, block_S, block_K, s0, N, K::MK) where {T, M <: AbstractMatrix{T}, MK <: Nothing}
    As0      = similar(s0, size(block_A, 1))
    QB       = similar(block_A, size(block_Q, 1), size(block_B, 2))
    STB      = similar(block_A, size(block_S, 2), size(block_B, 2))
    B_Q_B    = similar(block_A, size(block_B, 2), size(block_B, 2))

    LinearAlgebra.mul!(As0, block_A, s0)
    LinearAlgebra.mul!(QB, block_Q, block_B)
    LinearAlgebra.mul!(STB, block_S', block_B)
    LinearAlgebra.mul!(B_Q_B, block_B', QB)

    # Define Hessian term so that H = B_Q_B
    LinearAlgebra.axpy!(1, block_R, B_Q_B)
    LinearAlgebra.axpy!(1, STB, B_Q_B)
    LinearAlgebra.axpy!(1, STB', B_Q_B)

    # Define linear term so that c = h
    h = similar(s0, size(block_B, 2))
    LinearAlgebra.axpy!(1, block_S, QB)
    LinearAlgebra.mul!(h, QB', As0)

    # Define linear term so that c0 = h0
    h0   = similar(block_A, 1,1)
    QAs0 = similar(block_A, size(block_Q, 1), 1)
    LinearAlgebra.mul!(QAs0, block_Q, As0)
    LinearAlgebra.mul!(h0, As0', QAs0)

    return (H = B_Q_B, c = h, c0 = h0 ./ T(2))
end

function _build_H_blocks(block_Q, block_R, block_A, block_B, block_S, block_K, s0, N, K::MK) where {T, MK <: AbstractMatrix{T}}

    As0      = similar(s0, size(block_A, 1))
    RK       = similar(block_A, size(block_R, 1), size(block_K, 2))
    RKB      = similar(block_A, size(block_R, 1), size(block_B, 2))
    SK       = similar(block_A, size(block_S, 1), size(block_K, 2))
    SKB      = similar(block_A, size(block_S, 1), size(block_B, 2))
    STB      = similar(block_A, size(block_S, 2), size(block_B, 2))
    KTSTB    = similar(block_A, size(block_K, 2), size(block_B, 2))
    KTRK     = similar(block_A, size(block_K, 2), size(block_K, 2))
    QB       = similar(block_A, size(block_Q, 1), size(block_B, 2))
    KTRK_B   = similar(block_A, size(block_K, 2), size(block_B, 2))
    B_Q_B    = similar(block_A, size(block_B, 2), size(block_B, 2))

    LinearAlgebra.mul!(As0, block_A, s0)
    LinearAlgebra.mul!(RK, block_R, block_K)
    LinearAlgebra.mul!(RKB, RK, block_B)
    LinearAlgebra.mul!(SK, block_S, block_K)
    LinearAlgebra.mul!(SKB, SK, block_B)
    LinearAlgebra.mul!(STB, block_S', block_B)
    LinearAlgebra.mul!(KTSTB, block_K', STB)
    LinearAlgebra.mul!(QB, block_Q, block_B)
    LinearAlgebra.mul!(KTRK, block_K', RK)
    LinearAlgebra.mul!(KTRK_B, KTRK, block_B)

    LinearAlgebra.axpy!(1, KTRK_B, QB)
    LinearAlgebra.axpy!(1, SKB, QB)
    LinearAlgebra.axpy!(1, KTSTB, QB)  # QB now equals QB + KTRKB + SKB + KTSTB

    # Define Hessian term so that H = B_Q_B
    LinearAlgebra.mul!(B_Q_B, block_B', QB)
    LinearAlgebra.axpy!(1, block_R, B_Q_B)
    LinearAlgebra.axpy!(1, RKB', B_Q_B)
    LinearAlgebra.axpy!(1, RKB, B_Q_B)
    LinearAlgebra.axpy!(1, STB, B_Q_B)
    LinearAlgebra.axpy!(1, STB', B_Q_B)

    # Define linear term so that c = h
    h = similar(s0, size(block_B, 2))
    LinearAlgebra.axpy!(1, block_S, QB)
    LinearAlgebra.axpy!(1, RK', QB)
    LinearAlgebra.mul!(h, QB', As0)

    # Define constant term sot hat c0 = h0
    hR_term = similar(block_A, 1, 1) # = s0^T A^T K^T R K A s0
    hS_term = similar(block_A, 1, 1) # = s0^T A^T K^T S^T A s0 = s0^T A^T S K A s0
    hQ_term = similar(block_A, 1, 1) # = s0^T A^T Q A s0

    KTRKAs0 = similar(block_A, size(block_K, 2), 1)
    SKAs0   = similar(block_A, size(block_S, 1), 1)
    QAs0    = similar(block_A, size(block_Q, 1), 1)

    LinearAlgebra.mul!(KTRKAs0, KTRK, As0)
    LinearAlgebra.mul!(SKAs0, SK, As0)
    LinearAlgebra.mul!(QAs0, block_Q, As0)

    LinearAlgebra.mul!(hR_term, As0', KTRKAs0)
    LinearAlgebra.mul!(hS_term, As0', SKAs0)
    LinearAlgebra.mul!(hQ_term, As0', QAs0)

    h0 = hR_term ./ T(2) .+ hQ_term ./ T(2) .+ hS_term

    return (H = B_Q_B, c = h, c0 = h0)
end

function _build_G_blocks(block_A, block_B, block_E, block_F, block_K, block_gl, block_gu, s0, N)

    G = similar(block_A, size(block_F, 1), size(block_F, 2))

    As0  = similar(s0, size(block_A, 1))
    EAs0 = similar(s0, size(block_E, 1))
    FK   = similar(block_A, size(block_F, 1), size(block_K, 2))

    block_E_mul = similar(block_E)
    LinearAlgebra.copyto!(block_E_mul, block_E)

    LinearAlgebra.mul!(FK, block_F, block_K)
    LinearAlgebra.axpy!(1, FK, block_E_mul)
    LinearAlgebra.mul!(G, block_E_mul, block_B)
    LinearAlgebra.axpy!(1, block_F, G)

    LinearAlgebra.mul!(As0, block_A, s0)
    LinearAlgebra.mul!(EAs0, block_E_mul, As0)
    LinearAlgebra.axpy!(-1, EAs0, block_gl)
    LinearAlgebra.axpy!(-1, EAs0, block_gu)

    return (J = G, lcon = vec(block_gl), ucon = vec(block_gu), As0 = As0)
end

"""
    get_u(solution_ref, lqdm::SparseLQDynamicModel) -> u <: vector
    get_u(solution_ref, lqdm::DenseLQDynamicModel) -> u <: vector

Query the solution `u` from the solver. If `K = nothing`, the solution for `u` is queried from `solution_ref.solution`

If `K <: AbstractMatrix`, `solution_ref.solution` returns `v`, and `get_u` solves for `u` using the `K` matrix (and the `A` and `B` matrices if `lqdm <: DenseLQDynamicModel`)
"""
function get_u(
    solver_status,
    lqdm::SparseLQDynamicModel{T, V, M1, M2, M3, MK}
    ) where {T, V <: AbstractVector{T}, M1 <: AbstractMatrix{T}, M2 <: AbstractMatrix{T}, M3 <: AbstractMatrix{T}, MK <: AbstractMatrix{T}}

    solution = solver_status.solution
    ns       = lqdm.dynamic_data.ns
    nu       = lqdm.dynamic_data.nu
    N        = lqdm.dynamic_data.N
    K        = lqdm.dynamic_data.K

    u = zeros(T, nu * N)

    for i in 1:N
        start_v = (i - 1) * nu + 1
        end_v   = i * nu
        start_s = (i - 1) * ns + 1
        end_s   = i * ns

        Ks = zeros(T, size(K, 1), 1)

        s = solution[start_s:end_s]
        v = solution[(ns * (N + 1) + start_v):(ns * (N + 1) + end_v)]

        LinearAlgebra.mul!(Ks, K, s)
        LinearAlgebra.axpy!(1, v, Ks)

        u[start_v:end_v] = Ks
    end
    return u
end

function get_u(
    solver_status,
    lqdm::DenseLQDynamicModel{T, V, M1, M2, M3, M4, MK}
    ) where {T, V <: AbstractVector{T}, M1 <: AbstractMatrix{T}, M2 <: AbstractMatrix{T}, M3 <: AbstractMatrix{T}, M4 <: AbstractMatrix{T}, MK <: AbstractMatrix{T}}

    dnlp = lqdm.dynamic_data
    block_A = lqdm.blocks.A
    block_B = lqdm.blocks.B
    block_K = lqdm.blocks.K

    v = solver_status.solution

    As0 = zeros(T, size(block_A, 1))
    Bv  = zeros(T, size(block_B, 1))

    LinearAlgebra.mul!(As0, block_A, dnlp.s0)
    LinearAlgebra.mul!(Bv, block_B, v)
    LinearAlgebra.axpy!(1, As0, Bv)

    Ks = zeros(T, size(block_K, 1))

    LinearAlgebra.mul!(Ks, block_K, Bv)

    u = copy(v)
    LinearAlgebra.axpy!(1, Ks, u)

    return vec(u)
end

function get_u(
    solver_status,
    lqdm::SparseLQDynamicModel{T, V, M1, M2, M3, MK}
    ) where {T, V <: AbstractVector{T}, M1 <: AbstractMatrix{T}, M2 <: AbstractMatrix{T}, M3 <: AbstractMatrix{T}, MK <: Nothing}

    solution = solver_status.solution
    ns       = lqdm.dynamic_data.ns
    nu       = lqdm.dynamic_data.nu
    N        = lqdm.dynamic_data.N

    u = solution[(ns * (N + 1) + 1):end]
    return u
end

function get_u(
    solver_status,
    lqdm::DenseLQDynamicModel{T, V, M1, M2, M3, M4, MK}
    ) where {T, V <: AbstractVector{T}, M1 <: AbstractMatrix{T}, M2 <: AbstractMatrix{T}, M3 <: AbstractMatrix{T}, M4 <: AbstractMatrix{T}, MK <: Nothing}

    return copy(solver_status.solution)
end

"""
    get_s(solution_ref, lqdm::SparseLQDynamicModel) -> s <: vector
    get_s(solution_ref, lqdm::DenseLQDynamicModel) -> s <: vector

Query the solution `s` from the solver. If `lqdm <: SparseLQDynamicModel`, the solution is queried directly from `solution_ref.solution`
If `lqdm <: DenseLQDynamicModel`, then `solution_ref.solution` returns `u` (if `K = nothing`) or `v` (if `K <: AbstactMatrix`), and `s` is found form
transforming `u` or `v` into `s` using `A`, `B`, and `K` matrices.
"""
function get_s(
    solver_status,
    lqdm::SparseLQDynamicModel{T, V, M1, M2, M3, MK}
    ) where {T, V <: AbstractVector{T}, M1 <: AbstractMatrix{T}, M2 <: AbstractMatrix{T}, M3 <: AbstractMatrix{T}, MK <: Union{Nothing, AbstractMatrix}}

    solution = solver_status.solution
    ns       = lqdm.dynamic_data.ns
    N        = lqdm.dynamic_data.N

    s = solution[1:(ns * (N + 1))]
    return s
end
function get_s(
    solver_status,
    lqdm::DenseLQDynamicModel{T,V, M1, M2, M3, M4, MK}
    ) where {T, V <: AbstractVector{T}, M1 <: AbstractMatrix{T}, M2 <: AbstractMatrix{T}, M3 <: AbstractMatrix{T}, M4 <: AbstractMatrix{T}, MK <: Union{Nothing, AbstractMatrix}}

    dnlp = lqdm.dynamic_data
    block_A = lqdm.blocks.A
    block_B = lqdm.blocks.B

    v = solver_status.solution

    As0 = zeros(T, size(block_A, 1))
    Bv  = zeros(T, size(block_B, 1))

    LinearAlgebra.mul!(As0, block_A, dnlp.s0)
    LinearAlgebra.mul!(Bv, block_B, v)
    LinearAlgebra.axpy!(1, As0, Bv)

    return vec(Bv)
end

for field in fieldnames(LQDynamicData)
    method = Symbol("get_", field)
    @eval begin
        @doc """
            $($method)(LQDynamicData)
            $($method)(SparseLQDynamicModel)
            $($method)(DenseLQDynamicModel)
        Return the value $($(QuoteNode(field))) from LQDynamicData or SparseLQDynamicModel.dynamic_data or DenseLQDynamicModel.dynamic_data
        """
        $method(dyn_data::LQDynamicData) = getproperty(dyn_data, $(QuoteNode(field)))
    end
    @eval $method(dyn_model::SparseLQDynamicModel) = $method(dyn_model.dynamic_data)
    @eval $method(dyn_model::DenseLQDynamicModel)  = $method(dyn_model.dynamic_data)
    @eval export $method
end

for field in [:A, :B, :Q, :R, :Qf, :E, :F, :S, :K]
    method = Symbol("set_", field, "!")
    @eval begin
        @doc """
            $($method)(LQDynamicData, row, col, val)
            $($method)(SparseLQDynamicModel, row, col, val)
            $($method)(DenseLQDynamicModel, row, col, val)
        Set the value of entry $($(QuoteNode(field)))[row, col] to val for LQDynamicData, SparseLQDynamicModel.dynamic_data, or DenseLQDynamicModel.dynamic_data
        """
        $method(dyn_data::LQDynamicData, row, col, val) = (dyn_data.$field[row, col] = val)
    end
    @eval $method(dyn_model::SparseLQDynamicModel, row, col, val) = (dyn_model.dynamic_data.$field[row, col] = val)
    @eval $method(dyn_model::DenseLQDynamicModel, row, col, val)  = (dyn_model.dynamic_data.$field[row, col] = val)
    @eval export $method
end

for field in [:s0, :sl, :su, :ul, :uu, :gl, :gu]
    method = Symbol("set_", field, "!")
    @eval begin
        @doc """
            $($method)(LQDynamicData, index, val)
            $($method)(SparseLQDynamicModel, index, val)
            $($method)(DenseLQDynamicModel, index, val)
        Set the value of entry $($(QuoteNode(field)))[index] to val for LQDynamicData, SparseLQDynamicModel.dynamic_data, or DenseLQDynamicModel.dynamic_data
        """
        $method(dyn_data::LQDynamicData, index, val) = (dyn_data.$field[index] = val)
    end
    @eval $method(dyn_model::SparseLQDynamicModel, index, val) = (dyn_model.dynamic_data.$field[index] = val)
    @eval $method(dyn_model::DenseLQDynamicModel, index, val)  = (dyn_model.dynamic_data.$field[index] = val)
    @eval export $method
end

function fill_structure!(S::SparseMatrixCSC, rows, cols)
    count = 1
    @inbounds for col = 1:size(S, 2), k = S.colptr[col]:(S.colptr[col + 1] - 1)
        rows[count] = S.rowval[k]
        cols[count] = col
        count += 1
    end
end

function fill_coord!(S::SparseMatrixCSC, vals, obj_weight)
    count = 1
    @inbounds for col = 1:size(S, 2), k = S.colptr[col]:(S.colptr[col + 1] - 1)
        vals[count] = obj_weight * S.nzval[k]
        count += 1
    end
end

function NLPModels.hess_structure!(
    qp::SparseLQDynamicModel{T, V, M1, M2, M3},
    rows::AbstractVector{<:Integer},
    cols::AbstractVector{<:Integer},
) where {T, V, M1 <: SparseMatrixCSC, M2 <: SparseMatrixCSC, M3<: AbstractMatrix}
    fill_structure!(qp.data.H, rows, cols)
    return rows, cols
end

function NLPModels.hess_structure!(
    qp::DenseLQDynamicModel{T, V, M1, M2, M3},
    rows::AbstractVector{<:Integer},
    cols::AbstractVector{<:Integer},
) where {T, V, M1 <: Matrix, M2<: Matrix, M3<: Matrix}
    count = 1
    for j = 1:(qp.meta.nvar)
        for i = j:(qp.meta.nvar)
            rows[count] = i
            cols[count] = j
            count += 1
        end
    end
    return rows, cols
end

function NLPModels.hess_coord!(
    qp::SparseLQDynamicModel{T, V, M1, M2, M3},
    x::AbstractVector{T},
    vals::AbstractVector{T};
    obj_weight::Real = one(eltype(x)),
) where {T, V, M1 <: SparseMatrixCSC, M2 <: SparseMatrixCSC, M3 <: Matrix}
    NLPModels.increment!(qp, :neval_hess)
    fill_coord!(qp.data.H, vals, obj_weight)
    return vals
end

function NLPModels.hess_coord!(
    qp::DenseLQDynamicModel{T, V, M1, M2, M3},
    x::AbstractVector{T},
    vals::AbstractVector{T};
    obj_weight::Real = one(eltype(x)),
  ) where {T, V, M1 <: Matrix, M2 <: Matrix, M3 <: Matrix}
    NLPModels.increment!(qp, :neval_hess)
    count = 1
    for j = 1:(qp.meta.nvar)
      for i = j:(qp.meta.nvar)
        vals[count] = obj_weight * qp.data.H[i, j]
        count += 1
      end
    end
    return vals
  end

NLPModels.hess_coord!(
    qp::SparseLQDynamicModel,
    x::AbstractVector,
    y::AbstractVector,
    vals::AbstractVector;
    obj_weight::Real = one(eltype(x)),
) = NLPModels.hess_coord!(qp, x, vals, obj_weight = obj_weight)

NLPModels.hess_coord!(
    qp::DenseLQDynamicModel,
    x::AbstractVector,
    y::AbstractVector,
    vals::AbstractVector;
    obj_weight::Real = one(eltype(x)),
) = NLPModels.hess_coord!(qp, x, vals, obj_weight = obj_weight)

function NLPModels.jac_structure!(
    qp::SparseLQDynamicModel{T, V, M1, M2, M3},
    rows::AbstractVector{<:Integer},
    cols::AbstractVector{<:Integer},
) where {T, V, M1 <: SparseMatrixCSC, M2 <: SparseMatrixCSC, M3<: AbstractMatrix}
    fill_structure!(qp.data.A, rows, cols)
    return rows, cols
end

function NLPModels.jac_structure!(
    qp::DenseLQDynamicModel{T, V, M1, M2, M3},
    rows::AbstractVector{<:Integer},
    cols::AbstractVector{<:Integer},
  ) where {T, V, M1<: Matrix, M2 <: Matrix, M3 <: Matrix}
    count = 1
    for j = 1:(qp.meta.nvar)
      for i = 1:(qp.meta.ncon)
        rows[count] = i
        cols[count] = j
        count += 1
      end
    end
    return rows, cols
end

function NLPModels.jac_coord!(
    qp::SparseLQDynamicModel{T, V, M1, M2, M3},
    x::AbstractVector,
    vals::AbstractVector,
) where {T, V, M1 <: SparseMatrixCSC, M2 <: SparseMatrixCSC, M3 <: AbstractMatrix}
    NLPModels.increment!(qp, :neval_jac)
    fill_coord!(qp.data.A, vals, one(T))
    return vals
end

function NLPModels.jac_coord!(
    qp::DenseLQDynamicModel{T, V, M1, M2, M3},
    x::AbstractVector,
    vals::AbstractVector,
    ) where {T, V, M1 <: Matrix, M2 <: Matrix, M3 <: Matrix}
    NLPModels.increment!(qp, :neval_jac)
    count = 1
    for j = 1:(qp.meta.nvar)
      for i = 1:(qp.meta.ncon)
        vals[count] = qp.data.A[i, j]
        count += 1
      end
    end
    return vals
end

"""
    _build_H(Q, R, N; Qf = []) -> H

Build the (sparse) `H` matrix from square `Q` and `R` matrices such that
 z^T H z = sum_{i=1}^{N-1} s_i^T Q s + sum_{i=1}^{N-1} u^T R u + s_N^T Qf s_n .

# Examples
```julia-repl
julia> Q = [1 2; 2 1]; R = ones(1,1); _build_H(Q, R, 2)
6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 9 stored entries:
 1.0  2.0   ⋅    ⋅    ⋅    ⋅
 2.0  1.0   ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅   1.0  2.0   ⋅    ⋅
  ⋅    ⋅   2.0  1.0   ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅     ⋅
```

If `Qf` is not given, then `Qf` defaults to `Q`
"""
function _build_H(
    Q::M, R::M, N;
    Qf::M = Q,
    S::Union{M, Nothing} = nothing) where M <: AbstractMatrix
    ns = size(Q, 1)
    nu = size(R, 1)

    H = SparseArrays.sparse([],[],eltype(Q)[],(ns * (N + 1) + nu * N), (ns * (N+1) + nu * N))

    for i in 1:N
        range_Q = (1 + (i - 1) * ns): (i * ns)
        range_R = (ns * (N + 1) + 1 + (i - 1) * nu):(ns * (N + 1) + i * nu)
        H[range_Q, range_Q] = Q
        H[range_R, range_R] = R
        if S != nothing
            H[range_Q, range_R] = S
            H[range_R, range_Q] = S'
        end
    end

    H[(N * ns + 1):( N * ns + ns), (N * ns + 1):(N * ns + ns)] = Qf

    return H
end

"""
    _build_sparse_J1(A, B, N) -> J

Build the (sparse) `J` matrix or a linear model from `A` and `B` matrices such that
0 <= Jz <= 0 is equivalent to s_{i+1} = As_i + Bs_i for i = 1,..., N-1

# Examples
```julia-repl
julia> A = [1 2 ; 3 4]; B = [5 6; 7 8]; _build_J(A,B,3)
4×12 SparseArrays.SparseMatrixCSC{Float64, Int64} with 20 stored entries:
 1.0  2.0  -1.0    ⋅     ⋅     ⋅   5.0  6.0   ⋅    ⋅    ⋅    ⋅
 3.0  4.0    ⋅   -1.0    ⋅     ⋅   7.0  8.0   ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    1.0   2.0  -1.0    ⋅    ⋅    ⋅   5.0  6.0   ⋅    ⋅
  ⋅    ⋅    3.0   4.0    ⋅   -1.0   ⋅    ⋅   7.0  8.0   ⋅    ⋅
```
"""
function _build_sparse_J1(A::M, B::M, N) where M <: AbstractMatrix

    ns = size(A, 2)
    nu = size(B, 2)

    J1 = SparseArrays.sparse([], [], eltype(A)[], ns * N, (ns* (N + 1) + nu * N))

    neg_ones = .-Matrix(LinearAlgebra.I, ns, ns)

    for i in 1:N
        row_range  = (ns * (i - 1) + 1):(i * ns)
        Acol_range = (ns * (i - 1) + 1):(i * ns)
        Bcol_range = (ns * (N + 1) + 1 + (i - 1) * nu):(ns * (N + 1) + i * nu)
        J1[row_range, Acol_range] = A
        J1[row_range, Bcol_range] = B

        Icol_range = (ns * i + 1):(ns * (i + 1))

        J1[row_range, Icol_range] = neg_ones
    end

    return J1
end

function _build_sparse_J2(E, F, N)
    ns = size(E, 2)
    nu = size(F, 2)
    nc = size(E, 1)

    J2 = SparseArrays.sparse([],[], eltype(E)[], N * nc, ns * (N + 1) + nu * N)

    if nc != 0
        for i in 1:N
            row_range   = (1 + nc * (i - 1)):(nc * i)
            col_range_E = (1 + ns * (i - 1)):(ns * i)
            col_range_F = (ns * (N + 1) + 1 + nu * (i - 1)):(ns * (N + 1) + nu * i)

            J2[row_range, col_range_E] = E
            J2[row_range, col_range_F] = F
        end
    end

    return J2

end

function _build_sparse_J3(K, N, uu, ul)

    # Remove algebraic constraints if u variable is unbounded on both upper and lower ends
    nu = length(ul)
    ns = size(K, 2)

    bool_vec        = (ul .!= -Inf .|| uu .!= Inf)
    num_real_bounds = sum(bool_vec)

    J3 = SparseArrays.sparse([],[],eltype(K)[], nu * N, ns * (N + 1) + nu * N)
    I_mat = Matrix(LinearAlgebra.I, nu, nu)

    full_bool_vec = fill(true, nu * N)

    lcon3 = similar(ul, nu * N); fill!(lcon3, 0)
    ucon3 = similar(ul, nu * N); fill!(ucon3, 0)

    for i in 1:N
        row_range   = (nu * (i - 1) + 1):(nu * i)
        K_col_range = (ns * (i - 1) + 1):(ns * i)
        I_col_range = (ns * (N + 1) + 1 + nu * (i - 1)):(ns * (N + 1) + nu * i)
        J3[row_range, K_col_range] = K
        J3[row_range, I_col_range] = I_mat

        lcon3[row_range] .= ul
        ucon3[row_range] .= uu
        full_bool_vec[row_range] = bool_vec
    end

    return J3[full_bool_vec, :], lcon3[full_bool_vec], ucon3[full_bool_vec]
end

end # module


text:
using PyPlot

function test_sampling(graph = Graph(1000, 1000))
    const stepsize = 0.05
    const range = 0.0:stepsize:1.0

    dist1 = makedistributions(-2 ./ degree(graph) + 2 * degree(graph))
    dist2 = φ -> fill(Bernoulli(φ), nv(graph))

    # println(map(φ -> params(dist1(φ)[1]), range))
    # println(10 ./ degree(graph))

    results1 = samplepercolations(graph, dist1, 10, stepsize) do g
        reduce(max, 0, length.(connected_components(g)))
    end

    println(percolation_area(results1))
    # println(results1[1, :])
    plot(0.0:stepsize:1.0, smoothcurves(normalizecurves(results1)), color = "r")

    results2 = samplepercolations(graph, dist2, 10, stepsize) do g
        reduce(max, 0, length.(connected_components(g)))
    end

    println(percolation_area(results2))
    # println(results2[1, :])
    plot(0.0:stepsize:1.0, smoothcurves(normalizecurves(results2)), color = "g")
end

function test_annealing(graph = Graph(1000, 1000))
    error_trace = Float64[]
    infos = [:degree, :local_clustering]

    function log_annealing(;kwargs...)
        args = Dict(kwargs)
        push!(error_trace, args[:f_new])

        if args[:steps] % 10 == 0
            println(args[:steps], " steps, T = ", args[:T])
            println("\tParameters: ", args[:x_new])
            println("\tError: ", args[:f_new])

            if args[:updated]
                dist = informed_probabilities(graph, args[:x_new], infos)
                curves = samplepercolations(graph, dist, 1, 0.05) do g
                    reduce(max, 0, length.(connected_components(g)))
                end
                plot(0.0:0.05:1.0, smoothcurves(normalizecurves(curves)))
            end
        end
    end

    params = annealing(percolation_loss(graph, infos, 10), fill(0.0, length(infos)),
                       10.0, 0.9, 0.01, (s, v) -> s > 500;
                       debug_callback = log_annealing)

    figure()
    plot(error_trace)
end

ext: jl


text:
include("acoustics.jl")
include("circle.jl")
include("concentric_capsule.jl")
include("sphere.jl")
include("source.jl")
include("boundary_data.jl")


text:
@testset "307.range-sum-query-mutable.jl" begin
    @testset "307 case 1" begin
        ST = SegmentTree([1, 3, 5])
        @test sum_range(ST, 1, 3) == 9
        update!(ST, 2, 2)
        @test sum_range(ST, 1, 3) == 8
    end
end


text:
function saveSlam(slamwrapper::SLAMWrapper; filename::AbstractString="tempSlam.jld")
  saveslam = deepcopy(slamwrapper);
  @warn "Saving factor graph and landmark index -- unresolved issue with saving current tree, but can be reccomputed from factor graph with IncrementalInference.wipeBuildNewTree!(...)."
  saveslam.tree = Union{}
  jldopen(filename,"w") do file
    write(file, "slam", saveslam)
  end
  nothing
end

function loadSlam(; filename::AbstractString="tempSlam.jld")
  d = jldopen(filename,"r") do file
    read(file, "slam", slam)
  end
  return d
end

function haselement(arr::Vector{T}, val::T) where {T}
  for a in arr
    if a == val
      return true
    end
  end
  return false
end


text:
using Manopt, Manifolds, Test
@testset "Alternating Gradient Descent" begin
    # Note that this is mereely an alternating gradient descent toy example
    M = Sphere(2)
    N = M × M
    data = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
    function F(N, x)
        return 1 / 2 *
               (distance(N[1], x[N, 1], data[1])^2 + distance(N[2], x[N, 2], data[2])^2)
    end
    gradF1(N, x) = -log(N[1], x[N, 1], data[1])
    gradF1!(N, Y, x) = (Y .= -log(N[1], x[N, 1], data[1]))
    gradF2(N, x) = -log(N[2], x[N, 2], data[2])
    gradF2!(N, Y, x) = (Y .= -log(N[2], x[N, 2], data[2]))
    gradF(N, x) = ProductRepr([-log(N[i], x[N, i], data[i]) for i in [1, 2]]...)
    function gradF!(N, Y, x)
        log!(N[1], Y[N, 1], x[N, 1], data[1])
        log!(N[2], Y[N, 2], x[N, 2], data[2])
        return Y .*= -1
    end
    x = ProductRepr([0.0, 0.0, 1.0], [0.0, 0.0, 1.0])

    @testset "Test gradient access" begin
        Pf = AlternatingGradientProblem(N, F, gradF)
        Pv = AlternatingGradientProblem(N, F, [gradF1, gradF2])
        Pf! = AlternatingGradientProblem(N, F, gradF!; evaluation=MutatingEvaluation())
        Pv! = AlternatingGradientProblem(
            N, F, [gradF1!, gradF2!]; evaluation=MutatingEvaluation()
        )
        for P in [Pf, Pv, Pf!, Pv!]
            Y = zero_vector(N, x)
            @test get_gradient(P, x)[N, 1] == gradF(N, x)[N, 1]
            @test get_gradient(P, x)[N, 2] == gradF(N, x)[N, 2]
            get_gradient!(P, Y, x)
            @test Y[N, 1] == gradF(N, x)[N, 1]
            @test Y[N, 2] == gradF(N, x)[N, 2]
            @test get_gradient(P, 1, x) == gradF(N, x)[N, 1]
            @test get_gradient(P, 2, x) == gradF(N, x)[N, 2]
            Y = zero_vector(N, x)
            get_gradient!(P, Y[N, 1], 1, x)
            @test Y[N, 1] == gradF(N, x)[N, 1]
            get_gradient!(P, Y[N, 2], 2, x)
            @test Y[N, 2] == gradF(N, x)[N, 2]
        end
    end
    @testset "Test high level interface" begin
        y2 = allocate(x)
        copyto!(N, y2, x)
        y3 = allocate(x)
        copyto!(N, y3, x)
        y = alternating_gradient_descent(
            N, F, [gradF1!, gradF2!], x; order_type=:Linear, evaluation=MutatingEvaluation()
        )
        alternating_gradient_descent!(
            N,
            F,
            [gradF1!, gradF2!],
            y2;
            order_type=:Linear,
            evaluation=MutatingEvaluation(),
        )
        @test isapprox(N, y, y2)
        o = alternating_gradient_descent!(
            N,
            F,
            [gradF1!, gradF2!],
            y3;
            evaluation=MutatingEvaluation(),
            order_type=:Linear,
            return_options=true,
        )
        @test isapprox(N, y, o.x; atol=10^-3)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c930a0fa-4ef8-4963-b4b2-0d576ecdde7e
using Pkg, DrWatson

# ╔═╡ 919b41e5-7e43-45e0-8676-373b77ac6208
begin
	@quickactivate "NumericalMethodsforEngineers"
	using Roots, Plots
	using Symbolics, SymbolicUtils, Latexify, LaTeXStrings, ModelingToolkit
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ c34bdc36-5013-4096-9cb2-581ad9e02239
md" ## f-03-00.jl"

# ╔═╡ 63d01a6a-e8ae-4a94-9ecf-ad5768ff9ab0
begin
	@variables t z
	z = t^3 - t - 1
	L(z)
end

# ╔═╡ 4640e464-f75c-40fd-a26c-1f90363a6424
L(substitute(z, Dict([t => 2])))

# ╔═╡ 0cdd79ec-8e5c-4861-8b4d-88e1a8a426c8
f_j = build_function(z, [t], expression=Val{false});

# ╔═╡ c2d89d0e-95c2-4b25-8fff-68dca427c87d
f_j(2)

# ╔═╡ e626a4f9-97d5-4213-9f66-a17c24e37148
begin
	x = -2.0:0.01:2.0
	plot(x, f_j.(x))
end

# ╔═╡ d0d36a61-76f8-4017-9919-74bacb62e882
begin
	D = Differential(t) # symbolic representation of derivative
	L(D(z))
end

# ╔═╡ 08be6434-a03d-4167-abac-a54dec3ae240
begin
	fp_s = expand_derivatives(D(z))
	L(fp_s)
end

# ╔═╡ 4caaff50-acc2-438c-929b-68445d1b3d67
fp_j = build_function(fp_s, [t], expression=Val{false});

# ╔═╡ 5595773c-5510-48bf-94fd-216822419b2d
fp_j(2)

# ╔═╡ 7e58479c-48b5-4330-b312-1eb5f7c8796b
begin
	fj_s = Symbolics.jacobian([z], [t]; simplify=true)
	L(fj_s)
end

# ╔═╡ 57abf1cc-b6b5-4bb3-b80b-a3b8e41a8609
fj_j = build_function(fj_s[1], [t], expression=Val{false});

# ╔═╡ fe9b5719-017f-4f3b-b78f-16680a155e44
fj_j(2)

# ╔═╡ 1dd78a44-b910-4bb2-8af0-9a1ab3ee4a51
Roots.newton(f_j, fp_j, 8)        # 8.613169456441398

# ╔═╡ 26cd9127-ae9e-4f09-97e4-d10f514e2164
md" ## End of f-03-00.jl"

# ╔═╡ Cell order:
# ╟─c34bdc36-5013-4096-9cb2-581ad9e02239
# ╠═c930a0fa-4ef8-4963-b4b2-0d576ecdde7e
# ╠═919b41e5-7e43-45e0-8676-373b77ac6208
# ╠═63d01a6a-e8ae-4a94-9ecf-ad5768ff9ab0
# ╠═4640e464-f75c-40fd-a26c-1f90363a6424
# ╠═0cdd79ec-8e5c-4861-8b4d-88e1a8a426c8
# ╠═c2d89d0e-95c2-4b25-8fff-68dca427c87d
# ╠═e626a4f9-97d5-4213-9f66-a17c24e37148
# ╠═d0d36a61-76f8-4017-9919-74bacb62e882
# ╠═08be6434-a03d-4167-abac-a54dec3ae240
# ╠═4caaff50-acc2-438c-929b-68445d1b3d67
# ╠═5595773c-5510-48bf-94fd-216822419b2d
# ╠═7e58479c-48b5-4330-b312-1eb5f7c8796b
# ╠═57abf1cc-b6b5-4bb3-b80b-a3b8e41a8609
# ╠═fe9b5719-017f-4f3b-b78f-16680a155e44
# ╠═1dd78a44-b910-4bb2-8af0-9a1ab3ee4a51
# ╟─26cd9127-ae9e-4f09-97e4-d10f514e2164


text:
using PooksoftOptionsKit
using Dates

# what is the stock price?
stock_price = 54.17
implied_volatility = 0.67
DTE = (31.0/365.0)
number_of_levels = 12
risk_free_rate = 0.15 # this is in percent
dividend_rate = 0.0 # this is in percent
strike_price = 60.0
premium_value = 2.06

# setup options calculation -
amd_option_parameters = PSOptionKitPricingParameters(implied_volatility, DTE, number_of_levels, risk_free_rate, dividend_rate)

# setup asset contract -
callOptionContract = PSCallOptionContract("AMD", Date(2020,6,25), strike_price, premium_value, 1)

# estimate the implied volatility -
iv = estimate_implied_volatility(callOptionContract,amd_option_parameters,stock_price)


text:
export intersects2D
function intersects2D(lattice2D::Lattice; eps = 1E-10)
    lines = [edgeline(lattice2D.locations, i) for i in edges(lattice2D.graph)]
    nedges = ne(lattice2D.graph)
    buffer = []
    for i = 1:nedges-1
        for j = i+1:nedges
            isinter, point = intersects(lines[i], lines[j])
            if isinter
                tp = (lines[i][1], lines[i][2], lines[j][1], lines[j][2])
                if norm(tp[1] - point) > eps &&
                   norm(tp[2] - point) > eps &&
                   norm(tp[3] - point) > eps &&
                   norm(tp[4] - point) > eps
                    push!(buffer, point)
                end
            end
        end
    end
    return buffer
end

max_stars_repo_name: qiyang-ustc/FTensorKit

max_issues_repo_name: qiyang-ustc/FTensorKit

max_forks_repo_name: qiyang-ustc/FTensorKit


text:
export Dykstra_prox_parallel
function Dykstra_prox_parallel(
                                        x::Vector{TF},
                                        P::Vector{Any},
                                        P_sub,
                                        TD_OP,
                                        closed_form,
                                        maxit_dyk,
                                        dyk_feas_tol,
                                        obj_dyk_tol
                                        ) where {TF<:Real}

#this script computes the projection of x0 onto the intersection of and
#arbitrary number of closed and convex sets. May still work in case
#nonconvex sets are used.
# Solves: min_x 0.5||x-x0||_2^2 s.t. x in intersection of m convex sets.
#
# This algorithm allows for parallel evaluation of the prox part
# (projections). This version does it in serial anyway, so no Matlab
# Parallel computing toolbox is required.
#
# input:
#       x                   -   vector to be projection onto the intersection
#       P                   -   contains projectors. P{1},P{2}...P{m} as generated by setup_constraint.m in folder :
#       options_PARSDMM.
#           options_PARSDMM.tol     -   tolerance based on residual as stopping condition
#           options_PARSDMM.maxIt   -   maximum number of Parallel proximal Dykstra-like iterations
#           options_PARSDMM.minIt   -   minimum number of Parallel proximal Dykstra-like iterations
#           options_PARSDMM.log_vec -   if (=1), saves all vectors at every iteration for analysis purposes.
#           options_PARSDMM.feas_tol-   feasibility tolerance for warning message, not for stopping condition
#           options_PARSDMM.evol_rel_tol- tolerance on relative evolution between iterations: exit if norm(x-x_old)/norm(x_old) becomes too small
#
# output:
#       x           -   result
#       res         -   relative residual
#       prog_rel    -   relative evolution (movement) per iteration

# Author: Bas Peters

# Initialize
N       = length(x)
m = length(P) #number of constraints
ptp   = Vector{TF}(undef,N)
x_old = Vector{TF}(undef,N)
copyto!(ptp,x)
copyto!(x_old,x)
z=Vector{Vector{TF}}(undef,m)
p=Vector{Vector{TF}}(undef,m)

obj=zeros(TF,maxit_dyk)
rel_feasibility_err=zeros(TF,maxit_dyk+1,m)
evol_x = zeros(TF,maxit_dyk+1,m);

l1_P=zeros(Int64,maxit_dyk)
bounds_P=zeros(Int64,maxit_dyk)
ARADMM_it=zeros(Int64,maxit_dyk)
cg_it=zeros(Int64,maxit_dyk)

for i=1:m
    z[i]=zeros(TF,N)
    copyto!(z[i],x)
    p[i]=zeros(TF,N)
end
omega   = 1.0 ./ m; #weights are hardcored here

#initial feasibility (transform-domain, relative)
#check distance to feasibility
for i=1:m
  rel_feasibility_err[1,i]=norm(P_sub[i](TD_OP[i]*x)-TD_OP[i]*x)/norm(TD_OP[i]*x)
end

## Main loop
for n=1:maxit_dyk

    #evaluate prox
    temp1=0
    temp2=0
    temp3=0
    temp4=0
    for i=1:m #make this loop parallel for actual parallel proximal evaluation
        if closed_form[i]==false
          (p[i],log_PARSDMM) = P[i](z[i])
          if m==2; temp1= length(log_PARSDMM.obj); end
          if m==3; temp2= length(log_PARSDMM.obj); end
          if m==2; temp3= sum(log_PARSDMM.cg_it); end
          if m==3; temp4= sum(log_PARSDMM.cg_it); end
          if m==2; bounds_P[n]= length(log_PARSDMM.obj); end
          if m==3; l1_P[n] = length(log_PARSDMM.obj); end
        else
          copyto!(p[i],P[i](z[i]))
        end
  end
#   [ (p[i],log_PARSDMM) = P[i](z[i])
#         if m==2; temp1= length(log_PARSDMM.obj); end
#         if m==3; temp2= length(log_PARSDMM.obj); end
#         if m==2; temp3= sum(log_PARSDMM.cg_it); end
#         if m==3; temp4= sum(log_PARSDMM.cg_it); end
#         if m==2; bounds_P[n]= length(log_PARSDMM.obj); end
#         if m==3; l1_P[n] = length(log_PARSDMM.obj); end
#       else
#         copy!(p[i],P[i](z[i]))
#       end
# end
  ARADMM_it[n] = max(temp1,temp2)
  cg_it[n] = max(temp3,temp4)

    #averaging step
    if n>1; copyto!(x_old,x); end;
    fill!(x,TF(0.0))
    for i=1:m
        x .= x .+ omega.*p[i];
    end

    #log objective
    obj[n] = 0.5*norm(ptp-x,2).^2;
    evol_x[n]=norm(x_old-x)/norm(x)

    #updating
    for i=1:m
        z[i] .= x .+ z[i] .- p[i];
    end

    #check distance to feasibility
    for i=1:m
      rel_feasibility_err[n+1,i]=norm(P_sub[i](TD_OP[i]*x)-TD_OP[i]*x)/norm(TD_OP[i]*x)
    end

    #stop if objective value does not change and x is sufficiently feasible for all sets
    if n>1 && (maximum(rel_feasibility_err[n+1,:])<dyk_feas_tol) && ( ( (obj[n]-obj[n-1])./obj[n] ) < obj_dyk_tol )
        println("stationary objective and reached feasibility, exiting Dykstra (iteration ",n,")")
        obj=obj[1:n]
        cg_it=cg_it[1:n]
        l1_P=l1_P[1:n]
        ARADMM_it=ARADMM_it[1:n]
        evol_x = evol_x[1:n]
        rel_feasibility_err=rel_feasibility_err[1:n+1,:]
        return x,obj,rel_feasibility_err,cg_it,ARADMM_it,l1_P,bounds_P,evol_x
    end
end

#obj=obj[1:n]
#rel_feasibility_err=rel_feasibility_err[1:n,:]

return x,obj,rel_feasibility_err,cg_it,ARADMM_it,l1_P,bounds_P,evol_x
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module AlphaVantage

const HTTP_SCHEME = "https"
const HTTP_HOST = "alphavantage.co/query"
const ENV_VARIABLE_API_KEY = "ALPHA_VANTAGE_API_KEY"
const DEFAULT_API_KEY = "demo"

using HTTP
using CSV, JSON
using DataFrames
using ArgCheck

include("av_types.jl")
include("av_http.jl")
include("av_handlers.jl")

export
  AlphaVantageClient

export
  AVFunction,
    AVTimeSeries,
      TIME_SERIES_INTRADAY,
      TIME_SERIES_INTRADAY_EXTENDED,
      TIME_SERIES_DAILY,
      TIME_SERIES_DAILY_ADJUSTED,
      TIME_SERIES_WEEKLY,
      TIME_SERIES_WEEKLY_ADJUSTED,
      TIME_SERIES_MONTHLY,
      TIME_SERIES_MONTHLY_ADJUSTED,
      GLOBAL_QUOTE,
      SYMBOL_SEARCH,
    AVFundamental,
      INCOME_STATEMENT,
      BALANCE_SHEET,
      CASH_FLOW,
      EARNINGS,
      OVERVIEW,
      LISTING_STATUS,
      EARNINGS_CALENDAR,
      IPO_CALENDAR,
    AVCurrency,
      CURRENCY_EXCHANGE_RATE,
      FX_INTRADAY,
      FX_DAILY,
      FX_WEEKLY,
      FX_MONTHLY,
      AVC_PRICE,
      CRPYTO_INTRADAY,
      DIGITAL_CURRENCY_DAILY,
      DIGITAL_CURRENCY_WEEKLY,
      DIGITAL_CURRENCY_MONTHLY,
    AVEconomicIndicator,
      REAL_GDP,
      REAL_GDP_PER_CAPITA,
      TREASURY_YIELD,
      FEDERAL_FUNDS_RATE,
      INFLATION,
      INFLATION_EXPECTATION,
      CONSUMER_SENTIMENT,
      RETAIL_SALES,
      DURABLES,
      UNEMPLOYMENT,
      NONFARM_PAYROLL,
    AVTechnicalIndicator,
      SMA,
      EMA,
      WMA,
      DEMA,
      TEMA,
      TRIMA,
      KAMA,
      MAMA,
      VWAP,
      T3,
      MACD,
      MACDEXT,
      STOCH,
      STOCHF,
      RSI,
      STOCHRSI,
      WILLR,
      ADX,
      ADXR,
      APO,
      PPO,
      MOM,
      BOP,
      CCI,
      CMO,
      ROC,
      ROCR,
      AROON,
      AROONOSC,
      MFI,
      TRIX,
      ULTOSC,
      DX,
      MINUS_DI,
      PLUS_DI,
      MINUS_DM,
      PLUS_DM,
      BBANDS,
      MIDPOINT,
      MIDPRICE,
      SAR,
      TRANGE,
      ATR,
      NATR,
      AD,
      ADOSC,
      OBV,
      HT_TRENDLINE,
      HT_SINE,
      HT_TRENDMODE,
      HT_DCPERIOD,
      HT_DCPHASE,
      HT_PHASOR

end # module


text:
# Usage:
# julia --project=benchmark
# ```julia
# using PkgBenchmark, Convolution2D
# rst = PkgBenchmark.benchmarkpkg(Convolution2D)
# open("benchmark_report.md", "w") do io
#     PkgBenchmark.export_markdown(io, rst)
# end
# ```
using BenchmarkTools
using ImageFiltering
using ImageTransformations
using TestImages
using Convolution2D

const SUITE = BenchmarkGroup()

X = Float32.(imresize(testimage("cameraman"), (255, 255)))
kern = [0 1 0
    1 -4 1
    0 1 0]
SUITE["Convolution2D_255x255_gray"] = @benchmarkable conv2D($X, $kern)
SUITE["ImageFiltering_255x255_gray"] = @benchmarkable imfilter($X, $kern)


text:
using Documenter, Persa

makedocs(
    modules = [Persa],
    sitename = "Persa.jl",
    authors = "Filipe Braida and contributors.",
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        canonical = "https://juliadocs.github.io/Documenter.jl/stable/",
        assets = ["assets/favicon.ico"],
        analytics = "UA-128580038-1",
        highlights = ["yaml"],
    ),
    pages = [
        "Home" => "index.md",
        "Getting Started" => "man/getting_started.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaRecsys/Persa.jl.git",
    target = "build",
    push_preview = true
)


text:
function hscale{T<:Number}(alpha, A::HMat{T})
    C = HMat{T}();
    C.height = A.height;
    C.width = A.width;
    C.trg = A.trg;
    C.src = A.src;
    C.level = A.level;
    C.blockType = A.blockType;

    C.EPS = A.EPS;
    C.MAXRANK = A.MAXRANK;
    C.MINN = A.MINN;

    if A.blockType == LOWRANK
        C.UMat = sign(alpha)*sqrt(abs(alpha))*A.UMat;
        C.VMat = sqrt(abs(alpha))*A.VMat;
    elseif A.blockType == DENSE
        C.DMat = alpha*A.DMat;
    elseif A.blockType == HMAT
        C.blockType = HMAT;
        C.childHMat = Array(HMat{T},4,4);
        for i = 1:4, j = 1:4
            C.childHMat[i,j] = hscale(alpha,A.childHMat[i,j]);
        end
    end
    return C;
end

function hscale!(alpha, A::HMat)
    if A.blockType == LOWRANK
        A.UMat *= sign(alpha)*sqrt(abs(alpha));
        A.VMat *= sqrt(abs(alpha));
    elseif A.blockType == DENSE
        A.DMat *= alpha;
    elseif A.blockType == HMAT
        for C in A.childHMat
            hscale!(alpha,C);
        end
    end
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
using Turing
using Turing.RandomMeasures
using DynamicPPL
using Random
using AutoGibbs
using StatsPlots

include("models.jl")

const DATA_RNG = MersenneTwister(424242)
const MODEL_SETUPS = Dict([
    "GMM" => (gmm_generate, gmm_example, gmm_tarray_example, :z, (:w, :μ)),
    "HMM" => (hmm_generate, hmm_example, hmm_tarray_example, :s, (:t, :m)),
    "IMM" => (imm_stick_generate, imm_stick_example, imm_stick_tarray_example, :z, (:μ, :v))
])
const DATA_SIZE = 10
const CHAIN_LENGTH = 500

function evaluate_graphic(
    model::String,
    samplers::String,
    stepsize::String,
    steps::String,
    particles::String
)
    data, chain = evaluate(
        model,
        samplers,
        parse(Float64, stepsize),
        parse(Int, steps),
        parse(Int, particles)
    )
    @show data
    @show chain
    plot(chain)
end

function evaluate_grid(
    model,
    samplers;
    stepsizes=(0.01, 0.025, 0.05, 0.1),
    steps=(5, 10),
    particles=(20, 40, 60),
    N=CHAIN_LENGTH)

    for ss in stepsizes, s in steps, p in particles
        println("Evaluating $model with $samplers")
        println("using stepsize $ss, $s steps, and $p particles")
        println("for chain of length $N")
        data, chain = evaluate(model, samplers, ss, s, p, N)
        show(ess(chain))
    end
end

function evaluate_marginalized(
    stepsizes=(0.01, 0.025, 0.05, 0.1),
    steps=(5, 10),
    N=CHAIN_LENGTH
)
    dataname, data = gmm_marginalized_generate(DATA_RNG, 10)
    model = gmm_marginalized_example(; dataname => data)
    @show model
    for ss in stepsizes, s in steps
        println("Using stepsize $ss, $s steps, for chain of length $N")
        chain = sample(model, HMC(ss, s), N)
        show(ess(chain))
    end
end

"""
    evaluate(model, samplers, stepsize, steps, particles)

Try out the model with specified samplers and parameters.  `stepsize` and `steps` are for HMC,
`particles` is for PG.  `samplers` is a string containing either `nuts` or `hmc`, and either
`ag` or `pg`.  Unneeded parameters are ignored, but must be present for the signature.

Example:
```julia
evaluate_grid("GMM", "ag,nuts", stepsizes=(nothing,), steps=(nothing,))
```
"""
function evaluate(model, samplers, stepsize, steps, particles, N)
    generate, example, tarray_example, p_discrete, p_continuous = MODEL_SETUPS[model]

    dataname, data = generate(DATA_RNG, DATA_SIZE)
    model_pg = tarray_example(; dataname => data)

    if occursin(r"nuts"i, samplers)
        println("NUTS for continuous sampler.")
        dummy_nuts = NUTS()
        ad, _, metric = typeof(dummy_nuts).parameters
        sampler_continuous = NUTS{ad, (p_continuous...,), metric}(
            map(f -> getfield(dummy_nuts, f), fieldnames(typeof(dummy_nuts)))...
        )
    elseif occursin(r"hmc"i, samplers)
        println("HMC for continuous sampler.")
        sampler_continuous = HMC(stepsize, steps, p_continuous...)
    else
        error("Unknown continuous sampler")
    end

    if occursin(r"ag"i, samplers)
        println("AG for discrete sampler. Compiling static conditional...")
        sampler_discrete = StaticConditional(model_pg, p_discrete)
        println("Done.")
    elseif occursin(r"pg"i, samplers)
        println("PG for discrete sampler.")
        sampler_discrete = PG(particles, p_discrete)
    else
        error("Unknown discrete sampler")
    end

    chain = sample(model_pg, Gibbs(sampler_continuous, sampler_discrete), N)
    return data, chain
end


text:
"""
# module NumerovMethod

# Examples

```jldoctest
julia>
```
"""
module NumerovMethod

using NumericalMethodsInQuantumMechanics.EigenvalueProblems: InitialCondition

export integrate

"""
    numerov_iter(y_prev, y, dx, gvec, svec)

Return the ``y[n + 1]`` step for the Numerov's method, given the current and previous step.

# Arguments
- `y_prev::Real`: the `y[n - 1]` element.
- `y::Real`: the `y[n]` element.
- `dx::Real`: the step length ``h``, need to be small.
- `gvec::AbstractArray{<:Real}`: stores `g[n - 1]`, `g[n]` and `g[n + 1]`.
- `svec::AbstractArray{<:Real}`: stores `s[n - 1]`, `s[n]` and `s[n + 1]`.
"""
function iter(yᵢ₋₁, yᵢ, dx, gvec, svec)
    gᵢ₋₁, gᵢ, gᵢ₊₁ = gvec
    sᵢ₋₁, sᵢ, sᵢ₊₁ = svec
    coeffᵢ₋₁ = -(1 + dx^2 / 12 * gᵢ₋₁)
    coeffᵢ = 2(1 - 5dx^2 / 12 * gᵢ)
    coeffᵢ₊₁ = 1 + dx^2 / 12 * gᵢ₊₁
    coeff_sᵢ = dx^2 / 12 * (sᵢ₋₁ + 10sᵢ + sᵢ₊₁)
    return (coeffᵢ * yᵢ + coeffᵢ₋₁ * yᵢ₋₁) / coeffᵢ₊₁
end
"""
    numerov_iter(y_prev, y, dx, gvec)

Same as `numerov_iter(y_prev, y, dx, gvec, svec)`, if ``s(x) ≡ 0`` on the domain.

# Arguments
- `y_prev::Real`: the `y[n - 1]` element.
- `y::Real`: the `y[n]` element.
- `dx::Real`: the step length, need to be small.
- `gvec::AbstractArray{<:Real}`: stores `g[n - 1]`, `g[n]` and `g[n + 1]`.
"""
iter(yᵢ₋₁, yᵢ, dx, gvec) = 2(12 - 5dx^2 * gvec[2]) / (12 + dx^2 * gvec[3]) * yᵢ - yᵢ₋₁

"""
    integrate(ic, r, gvec, svec)

Do the Numerov's method integration, return the solution ``y(x)``, given that ``g(x)`` and ``s(x)`` already
as vectors (already applied on ``x``).

# Arguments
- `ic::InitialCondition`: the initial condition ``y(x_0) = y_0`` and ``y'(x_0) = y_0'``, could be a guess.
- `r::AbstractRange{<:Real}`: a range, the domain ``x``.
- `gvec::AbstractArray{<:Real}`: the result of function ``g`` applied on ``x`` (range `r`).
- `svec::AbstractArray{<:Real}`: the result of function ``s`` applied on ``x`` (range `r`).
"""
function integrate(gvec, svec, ic::InitialCondition)
    if length(gvec) != length(svec)
        throw(DimensionMismatch("Integ"))
    end
    N = length(gvec)
    dx = inv(N)
    ϕ₀, ϕ′₀ = ic
    ϕ = [ϕ₀, ϕ′₀ * dx]  # ϕ₀, ϕ₁
    for i in 1:(N-2)
        ϕᵢ₊₂ = iter(ϕ[i], ϕ[i+1], dx, gvec[i:(i+2)], svec[i:(i+2)])
        push!(ϕ, ϕᵢ₊₂)
    end
    return ϕ
end
"""
    integrate(ic, r, g, s)

Same as `integrate(ic, r, gvec, svec)`, but `g` and `s` are two functions.

# Arguments
- `ic::InitialCondition`: the initial condition ``y(x_0) = y_0`` and ``y'(x_0) = y_0'``, could be a guess.
- `r::AbstractRange{<:Real}`: a range, the domain ``x``.
- `g::Function`: the function ``g``.
- `s::Function`: the function ``s``.
"""
function integrate(g::Function, s::Function, ic::InitialCondition, dx)
    vec = 0:dx:1
    gvec, svec = map(g, vec), map(s, vec)
    return integrate(gvec, svec, ic)
end

end


text:
module SliderStrings

Base.Experimental.@optlevel 1

using Redux
using CImGui
import ..ReduxImGui: get_label, get_string

# actions
abstract type AbstractSliderStringAction <: AbstractSyncAction end

"""
    Rename(label, new_label)
Change label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractSliderStringAction
    label::String
    new_label::String
end

"""
    ChangeStrings(label, strs)
Change widget strings to `strs`.
"""
struct ChangeStrings <: AbstractSliderStringAction
    label::String
    strs::Vector{String}
end

"""
    ChangeIndex(label, idx)
Change current strings to the one indexed by `idx`.
"""
struct ChangeIndex <: AbstractSliderStringAction
    label::String
    idx::Cint
end

# state
struct State <: AbstractImmutableState
    label::String
    strs::Vector{String}
    idx::Ref{Cint}
end
State(label::AbstractString, strs::Vector{<:AbstractString}) = State(label, strs, Ref{Cint}(1))

# reducers
slider_string(state::AbstractState, action::AbstractAction) = state
slider_string(state::Vector{<:AbstractState}, action::AbstractAction) = state
slider_string(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

function slider_string(state::Dict{String,State}, action::AbstractSliderStringAction)
    s = Dict{String,State}()
    for (k, v) in state
        s[k] = get_label(v) == action.label ? slider_string(v, action) : v
    end
    return s
end

slider_string(s::State, a::Rename) = State(a.new_label, s.strs, s.idx)
slider_string(s::State, a::ChangeStrings) = State(s.label, a.strs, s.idx)
slider_string(s::State, a::ChangeIndex) = State(s.label, s.strs, Ref{Cint}(a.idx))

# helper
"""
    SliderString(store::AbstractStore, get_state=Redux.get_state) -> Bool
Return `true` when triggered.
"""
function SliderString(store::AbstractStore, get_state = Redux.get_state)
    s = get_state(store)
    min, max = 1, length(s.strs)
    # str = min ≤ s.idx[] ≤ max ? s.strs[s.idx[]] : "Unknown"
    is_triggered = CImGui.SliderInt(get_label(s), s.idx, min, max, get_string(s))
    return is_triggered
end

get_label(s::State) = s.label
get_string(s::State) = s.strs[s.idx[]]

end # module


text:
push!(LOAD_PATH,"C:\\julia_code")

using GPbox
using Distributions
using Gadfly

temp = linspace(-5,5,50);
x = Array{Float64}(1,length(temp));
for i in eachindex(temp)
    x[1,i] = temp[i];
end
len = .3;
trainx = [-4, -1.6, -.8, 2.5, 3.9]
trainf = [-1, -.1, .4, .2, -1.2]

trainx = reshape(trainx,1,length(trainx));
trainf = reshape(trainf,1,length(trainf));

gpk = GPkernel("squared exp",Dict("len"=>len))
pred_mn, pred_cv = estGP(x,trainx,trainf,gpk);
pred_std = sqrt.(diag(pred_cv));

l1 = layer(x=x,
            y=pred_mn,
            Geom.line,
            Theme(default_color=colorant"violet"));

l2 = layer(
            x=trainx,
            y=trainf,
            Geom.point,
            style(default_color=colorant"green"));

l3 = layer(x=x,
            ymin=pred_mn-pred_std,
            ymax=pred_mn+pred_std,
            Geom.ribbon);

lyrs = PlotGPsamples1D(x,5,pred_mn,pred_cv);
plot(l1,l2,l3,lyrs...)


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

JLLWrappers.@generate_main_file_header("CImPlot")
JLLWrappers.@generate_main_file("CImPlot", UUID("278f1f02-c5a0-563f-9174-e8b244dda450"))
end  # module CImPlot_jll


text:
@inline function getVelocityDelta(joint::Translational3, body1::AbstractBody, body2::Body{T}, v::Union{Nothing,SVector{0,T}}) where T
    Δv = @SVector zeros(T,3)
    return Δv
end

@inline function getPositionDelta(joint::Translational3, body1::AbstractBody, body2::Body{T}, x::Union{Nothing,SVector{0,T}}) where T
    Δx = @SVector zeros(T,3)
    return Δx
end

@inline function minimalCoordinates(joint::Translational3, body1::AbstractBody{T}, body2::Body, No) where T
    SVector{0,T}()
end

@inline function g(joint::Translational3, body1::Body, body2::Body, Δt, No)
    vertices = joint.vertices
    getx3(body2, Δt) + vrotate(vertices[2], getq3(body2, Δt)) - (getx3(body1, Δt) + vrotate(vertices[1], getq3(body1, Δt)))
end

@inline function g(joint::Translational3, body1::Origin, body2::Body, Δt, No)
    vertices = joint.vertices
    getx3(body2, Δt) + vrotate(vertices[2], getq3(body2, Δt)) - vertices[1]
end

@inline function ∂g∂posa(joint::Translational3{T}, body1::Body, body2::Body, No) where T
    if body2.id == joint.cid
        q1 = body1.q[No]

        X = SMatrix{3,3,T,9}(-I)
        R = -2 * VRᵀmat(q1) * Rmat(Quaternion(joint.vertices[1])) * LVᵀmat(q1)

        return [X R]
    else
        return ∂g∂posa(joint)
    end
end

@inline function ∂g∂posb(joint::Translational3{T}, body1::AbstractBody, body2::Body, No) where T
    if body2.id == joint.cid
        q2 = body2.q[No]

        X = SMatrix{3,3,T,9}(I)
        R = 2 * VRᵀmat(q2) * Rmat(Quaternion(joint.vertices[2])) * LVᵀmat(q2)

        return [X R]
    else
        return ∂g∂posb(joint)
    end
end

@inline function ∂g∂vela(joint::Translational3{T}, body1::Body, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        q1 = body1.q[No]

        V = SMatrix{3,3,T,9}(-Δt * I)
        Ω = -2 * VRᵀmat(q1) * Lmat(q1) * Rᵀmat(ωbar(body1, Δt)) * Rmat(Quaternion(joint.vertices[1])) * derivωbar(body1, Δt)

        return [V Ω]
    else
        return ∂g∂vela(joint)
    end
end

@inline function ∂g∂velb(joint::Translational3{T}, body1::AbstractBody, body2::Body, Δt, No) where T
    if body2.id == joint.cid
        q2 = body2.q[No]

        V = SMatrix{3,3,T,9}(Δt * I)
        Ω = 2 * VRᵀmat(q2) * Lmat(q2) * Rᵀmat(ωbar(body2, Δt)) * Rmat(Quaternion(joint.vertices[2])) * derivωbar(body2, Δt)

        return [V Ω]
    else
        return ∂g∂velb(joint)
    end
end


text:
# find zero for probability density function f(x)
# where x > 0
# start at x0
function my_bisection(f, a, b)
    tol = 1e-4
    eps = 1e-5
    kmax = 100
    k = 1
    c = 0.
    while k < kmax
        # println("ITERATION $k")
        # flush(stdout)
        c = (a+b)/2
        if abs(f(c)) < eps || abs(a-b) < tol
            return c
        elseif sign(f(c)) == sign(f(a))
            a = c
            # println("Current status: [$a, $b], \n $(f(a)), $(f(b))")
            # flush(stdout)
        elseif sign(f(c)) == sign(f(b))
            b = c
            # println("Current status: [$a, $b], \n $(f(a)), $(f(b))")
            # flush(stdout)
        end
        k += 1
        # println("========= FINISHING =========")
        # flush(stdout)

    end
    return c
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# Type and constructors
"""
	(type) Circle{T<:AnyComplex} in the complex plane

Each `Circle` type is parameterized according to the common type of its complex input arguments.
"""
struct Circle{T<:AnyComplex} <: AbstractClosedCurve
	center::T
	radius::Float64
	ccw::Bool
end
"""
	Circle(zc,r,ccw=true)

Construct the circle with given center `zc`, radius `r`, and orientation (defaults to counterclockwise).

	Circle(a,b,c)

Construct the circle passing through the given three numbers. Orientation is determined so that the values are visited in the given order. If the three points are collinear (including when one of the given values is infinite), a `Line` is returned instead.
"""
Circle(z::AnyComplex,r::Real,ccw::Bool=true) = Circle{typeof(z)}(z,r,ccw)
Circle(z::Number,r::Real,ccw::Bool=true) = Circle(complex(float(z)),r,ccw)

# Construction by three points
function Circle(a::Number,b::Number,c::Number)
	a,b,c = promote( complex.(float.([a,b,c]))... )
	Circle(a,b,c)
end
function Circle(a::T,b::T,c::T) where {T<:AnyComplex}
	isinf(a) && return Line(b,c)
	isinf(b) && return Line(c,a)
	isinf(c) && return Line(a,b)
	# Use intersection of chord bisectors to find the center of the circle.
	w = (a-c)/2
	d1,d2 = a-b,c-b
	M = SMatrix{2,2}(real(d1),imag(d1),real(d2),imag(d2))
	if cond(M) > 0.1/eps(typeof(float(real(a))))
		# Collinear points
		return Line(a,b)
	else
		p =  M \ SVector(imag(w),-real(w))
		cen = (a+b)/2 - 1im*p[1]*d1
		ccw = isccw(a-cen,b-cen,c-cen)
		return Circle{T}(cen,abs(a-cen),ccw)
	end
end

# Complex type converters
for ctype in [:Spherical,:Polar,:Complex]
	@eval begin
		function $ctype(C::Circle{T}) where T<:AnyComplex
			Circle($ctype(C.center),C.radius,C.ccw)
		end
	end
end

# Required methods
function point(C::Circle,t::Real)
	z = C.ccw ? exp(2im*pi*t) : exp(-2im*pi*t)
	C.center + typeof(C.center)(z*C.radius)
end
(C::Circle)(t::Real) = point(C,t)

ispositive(C::Circle) = C.ccw

arclength(C::Circle) = 2π*C.radius

"""
	arg(C::Circle,z)

Find the parameter argument `t` such that `C(t)==z` is true.

This gives undefined results if `z` is not actually on the circle.
"""
function arg(C::Circle,z::Number)
	α = angle(z-C.center)/(2π)
	C.ccw ? mod(α,1) : mod(-α,1)
end

function unittangent(C::Circle{T},t::Real) where T <: AnyComplex
	C.ccw ? T( 1im*exp(2im*pi*t) ) : T( -1im*exp(-2im*pi*t) )
end
function tangent(C::Circle{T},t::Real) where T <: AnyComplex
	T( 2π*C.radius*unittangent(C,t) )
end

+(C::Circle,z::Number) = Circle(C.center+z,C.radius,C.ccw)
-(C::Circle) = Circle(-C.center,C.radius,C.ccw)
*(C::Circle,z::Number) = Circle(C.center*z,C.radius*abs(z),C.ccw)

"""
	inv(C)
Invert the circle `C` through the origin. In general the inverse is a `Circle`, though the result is a `Line` if `C` passes through the origin.
"""
function inv(C::Circle)
	w = 1 ./ point(C,[0,0.25,0.5])
	Circle(w...)
end

# Other methods
isfinite(::Circle) = true
conj(C::Circle) = Circle(conj(C.center),C.radius,!C.ccw)
reverse(C::Circle) = Circle(C.center,C.radius,!C.ccw)

"""
	isapprox(C1::Circle,C2::Circle; tol=<default>)
	C1 ≈ C2
Determine if `C1` and `C2` represent the same circle, irrespective of the type or values of its parameters. Identity is determined by agreement within `tol`, which is interpreted as the weaker of absolute and relative differences.
"""
function isapprox(C1::Circle,C2::Circle;tol=DEFAULT[:tol])
	return isapprox(C1.center,C2.center,rtol=tol,atol=tol) &&
	isapprox(C1.radius,C2.radius,rtol=tol,atol=tol)
end

function winding(C::Circle,z::Number)
	w =  abs(z-C.center) < C.radius ? 1 : 0
	C.ccw ? w : -w
end

"""
	dist(z,C::Circle)
Compute the distance from number `z` to the circle `C`.
"""
dist(z::Number,C::Circle) = abs( abs(z-C.center) - C.radius )

"""
	closest(z,C::Circle)
Find the point on circle `C` that lies closest to `z`.
"""
closest(z::Number,C::Circle) =	C.center + C.radius*sign(z - C.center)

"""
	reflect(z,C::Circle)
Reflect the value `z` across the circle `C`. (For reflection of a circle through a point, use translation and negation.)
"""
function reflect(z::Number,C::Circle)
	ζ = z-C.center
	ζ==0 ? convert(typeof(float(z)),Inf) : C.center + C.radius^2/conj(ζ)
end

function show(io::IO,C::Circle)
	orient = C.ccw ? "ccw" : "cw"
	print(IOContext(io,:compact=>true),"Circle(",C.center,",",C.radius,",",orient,")")
end
function show(io::IO,::MIME"text/plain",C::Circle{T}) where {T}
	orient = C.ccw ? "positively" : "negatively"
	print(io,"Circle{$T} in the complex plane:\n   centered at (",C.center,") with radius $(C.radius), $(orient) oriented")
end


text:
using PyPlot
using JLD
X = load("../data/erg_avg/Jar1_res.jld")
z_rho = X["z_rho"].s
rho = X["rho"]
X = load("../data/erg_avg/Jar2_res.jld")
z_rho = [z_rho, X["z_rho"].s]
rho = [rho, X["rho"]]
z_rho = vcat(z_rho...)
rho = vcat(rho...)

X = load("../data/erg_avg/Jar3_res.jld")
z_rho = [z_rho, X["z_rho"].s]
rho = [rho, X["rho"]]

z_rho = vcat(z_rho...)
rho = vcat(rho...)

X = load("../data/erg_avg/Jar4_res.jld")
z_rho = [z_rho, X["z_rho"].s]
rho = [rho, X["rho"]]

z_rho = vcat(z_rho...)
rho = vcat(rho...)

max_stars_repo_name: nishaChandramoorthy/Lipschitz

max_issues_repo_name: nishaChandramoorthy/Lipschitz

max_forks_repo_name: nishaChandramoorthy/Lipschitz


text:
@testset "B" begin
    @testset "b1900" begin
        @test b1900() == 2415020.31352
    end
    @testset "b1900" begin
        @test b1950() == 2433282.42345905
    end
    @testset "badkpv" begin
        try
            furnsh(path(CORE, :pck))
            @test !badkpv("SPICE.jl", "BODY399_RADII", "=", 3, 1, "N")
            @test_throws SpiceError badkpv("SPICE.jl", "BODY399_RADII", "=", 3, 1, "C")
        finally
            kclear()
        end
    end
    @testset "bltfrm" begin
        idset = bltfrm(-1)
        @test length(idset) >= 126
    end
    @testset "bodc2n" begin
        @test bodc2n(399) == "EARTH"
        @test bodc2n(typemax(Cint)) === nothing
    end
    @testset "bodc2s" begin
        @test bodc2s(399) == "EARTH"
        @test bodc2s(typemax(Cint)) == string(typemax(Cint))
    end
    @testset "boddef" begin
        boddef("Spaceball One", 999999)
        @test bodc2n(999999) == "Spaceball One"
    end
    @testset "bodfnd" begin
        try
            furnsh(path(CORE, :pck))
            @test bodfnd(399, "RADII") == true
            @test bodfnd(999999, "RADII") == false
        finally
            kclear()
        end
    end
    @testset "bodn2c" begin
        @test bodn2c("EARTH") == 399
        @test bodn2c("Norbert") === nothing
    end
    @testset "bods2c" begin
        @test bods2c("EARTH") == 399
        @test bods2c("888888") == 888888
        @test bodn2c("Norbert") === nothing
    end
    @testset "bodvcd" begin
        try
            furnsh(path(CORE, :pck))
            @test bodvcd(399, "RADII") == [6378.1366, 6378.1366, 6356.7519]
            @test bodvcd(399, "RADII") == [6378.1366, 6378.1366, 6356.7519]
            @test_throws SpiceError bodvcd(399, "Norbert")
            @test_throws SpiceError bodvcd(999999, "RADII")
        finally
            kclear()
        end
    end
    @testset "bodvrd" begin
        try
            furnsh(path(CORE, :pck))
            @test bodvrd("EARTH", "RADII") == [6378.1366, 6378.1366, 6356.7519]
            @test bodvrd("EARTH", "RADII") == [6378.1366, 6378.1366, 6356.7519]
            @test_throws SpiceError bodvrd("EARTH", "Norbert")
            @test_throws SpiceError bodvrd("Norbert", "RADII")
        finally
            kclear()
        end
    end
    @testset "brcktd/brckti" begin
        @test SPICE._brcktd(4.0, 1.0, 3.0) == clamp(4.0, 1.0, 3.0)
        @test SPICE._brckti(4, 1, 3) == clamp(4, 1, 3)
        arr = ["FEYNMAN", "BOHR", "EINSTEIN", "NEWTON", "GALILEO"]
        ord = sortperm(arr)
        for name in arr
            @test SPICE._bschoc(name, arr, ord) == findfirst(arr .== name)
            @test SPICE._bsrchc(name, sort(arr)) == findfirst(sort(arr) .== name)
        end
    end
    @testset "bschoi/bsrchd/bsrchi" begin
        @test SPICE._bschoi(1, [3, 1, 2], [2, 3, 1]) == findfirst([3, 1, 2] .== 1)
        @test SPICE._bsrchd(3.0, [1.0, 2.0, 3.0]) == findfirst([1.0, 2.0, 3.0] .== 3.0)
        @test SPICE._bsrchi(3, [1, 2, 3]) == findfirst([1, 2, 3] .== 3)
    end
end


text:
function sptPlot(path::ASCIIString, it::Int64; name="NULL", attribute="p", xlabel="Length", xunits="(m)", ylabel="Depth", yunits="(m)", ox=0, dx=5, oy=0, dy=5, cmap="seismic", clip=1.0, dpi=100, aspect="auto", wbox=6, hbox=6, cbar=false, labelsize=15, ticksize=15)
    spt = readSnapShot(path, it)
    nz = spt.nz ; nx    = spt.nx   ;
    ext= spt.ext; iflag = spt.iflag;
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    fig = plt.figure(figsize=(wbox, hbox), dpi=dpi, facecolor="w", edgecolor="k")
    if attribute == "p"
       d = reshape(spt.pz+spt.px, Nz, Nx)
    elseif attribute == "vz"
       d = reshape(spt.vz, Nz, Nx)
    elseif attribute == "vx"
       d = reshape(spt.vx, Nz, Nx)
    else
       error("wrong attribute")
    end
    vmax = quantile(vec(abs(d)), clip); vmin = -vmax
    im = plt.imshow(d, cmap=cmap, vmin=vmin, vmax=vmax, extent=[ox,ox+(size(d,2)-1)*dx,oy+(size(d,1)-1)*dy,oy], aspect=aspect)
    plt.xlabel(join([xlabel " " xunits]), fontsize=labelsize)
    plt.ylabel(join([ylabel " " yunits]), fontsize=labelsize)
    if cbar == true
      plt.colorbar()
    end
    ax = plt.gca()
    plt.setp(ax[:get_xticklabels](), fontsize=ticksize)
    plt.setp(ax[:get_yticklabels](), fontsize=ticksize)
    if (name == "NULL")
  		plt.show()
  	else
  		plt.savefig(name,dpi=dpi)
  		plt.close()
  	end
    return nothing
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
#: by David Thierry
using JuMP
using Clp
using DataFrames
using StatsPlots
tHorz = 30
m = Model()

#vCapCombTurb = 3.
vCapSteamTurb = 2.
vCapTransInter = 5.
vCapPcc = 20.
vCapComp = 1000.

rnum = rand(Float64, tHorz)

#@variable(m, 0 <= vCapCombTurb)
#@variable(m, 0 <= vPowCombTurb[0:tHorz - 1] <= vCapCombTurb)
@variable(m, 0 <= vPowCombTurb[0:tHorz - 1])

#@variable(m, 0 <= ePowSteamTurb[0:tHorz - 1] <= vCapSteamTurb)
@variable(m, 0 <= ePowSteamTurb[0:tHorz - 1])
#@variable(m, 0 <= vCapSteamTurb)

#@variable(m, rnum[i+1] <= ePowOut[i = 0:tHorz - 1] <= vCapTransInter)
@variable(m, ePowOut[i = 0:tHorz - 1] >= rnum[i + 1])
#@variable(m, 0 <= ePowOut[i = 0:tHorz - 1])
@variable(m, 0 <= ePowHp[0:tHorz - 1])
@variable(m, 0 <= ePowIp[0:tHorz - 1])
#@variable(m, 0 <= vCapTransInter)

@variable(m, 0 <= eFuel[0:tHorz - 1])

@variable(m, 0 <= eSteamIpLp[0:tHorz - 1])
@variable(m, 0 <= eSteamLp1[0:tHorz - 1])

@variable(m, 0 <= ePowLp1[0:tHorz - 1])
@variable(m, 0 <= ePowLp2[0:tHorz - 1])

@variable(m, 0 <= eCo2Fuel[0:tHorz - 1])

# Pcc
#@variable(m, 0 <= eCo2CapPcc[0:tHorz - 1] <= vCapPcc)
@variable(m, 0 <= eCo2CapPcc[0:tHorz - 1])
@variable(m, 0 <= eSteamUsePcc[0:tHorz - 1])
@variable(m, 0 <= ePowUsePcc[0:tHorz - 1])
@variable(m, 0 <= eCo2PccOut[0:tHorz - 1])
#@variable(m, 0 <= vCo2PccVent[0:tHorz - 1] <= 0.1)
@variable(m, 0 <= vCo2PccVent[0:tHorz - 1])
#vCo2PccVent = 0.0
@variable(m, 0 <= eCo2DacFlueIn[0:tHorz - 1])

#@variable(m, 0 <= a0[0:tHorz - 1])  # We could generalize this for any number of hours.
#@variable(m, 0 <= a1[0:tHorz])
#@variable(m, 0 <= a2[0:tHorz - 1])

#@variable(m, 0 <= aR0[0:tHorz - 1])
#@variable(m, 0 <= aR1[0:tHorz - 1])

#@variable(m, 0 <= sF[0:tHorz])
#@variable(m, 0 <= sS[0:tHorz])

@variable(m, 0 <= eCo2DacVent[0:tHorz - 1])

@variable(m, 0 <= eCo2Comp[0:tHorz - 1])
#@variable(m, 0 <= ePowUseComp[0:tHorz - 1] <= vCapComp)
@variable(m, 0 <= ePowUseComp[0:tHorz - 1])
#@variable(m, 0 <= vCapComp)
@variable(m, 0 <= eCo2Vent[0:tHorz - 1])

# Parameters
pHeatRateCombTur = 1.
# pSteamRate = 1.
pEmissFactor = 200.1

# pCo2CapRatePcc = 1.

#: Turbine Parameters
pFuelPowHp = 1.
pFuelPowIp = 1.
pFuelSteamIpLp = 1.
pHeatRateLp1 = 1.
pFuelPowLp = 1.
#
#
pSteamUseRatePcc = 0.001
pPowUseRatePcc = 0.022

pSteamUseRateDacAir = 0.00001
pSteamUseRateDacFlue = 0.00001

pPowUseRateDacAir = 0.0001
pPowUseRateDacFlue = 0.02

pHeatRateSteamTur = 1.

pCo2PccCapRate = .1
pSorbCo2Cap = 0.1

pPowUseRateDac = 0.001

pPowUseRateComp = 0.001
#

pCostInvCombTurb = 1e+02
pCostInvSteaTurb = 1e+02
pCostInvTransInter = 1e+02
pCostInvPcc = 1e+02
pCostInvDac = 1e+03
pCostInvComp = 1e+01

# Cost parameters.
pCostFuel = 100.
pEmissionPrice = 1.
pCo2TranspPrice = 1e+03

# Constraints
# 1
#@constraint(m, capCombTurIn[i = 0:tHorz - 1], vPowCombTurb[i] <= vCapCombTurb)
# 2
# @constraint(m, capSteamTurbIn[i = 0:tHorz - 1], ePowSteamTurb[i] <= vCapSteamTurb)
# 3
# @constraint(m, capTransmIntercIn[i = 0:tHorz - 1], ePowOut[i] <= vCapTransInter)
# 4
@constraint(m, fuelEq[i = 0:tHorz - 1], eFuel[i] == pHeatRateCombTur * vPowCombTurb[i])
# 5
@constraint(m, co2FuelEq[i = 0:tHorz - 1], eCo2Fuel[i] == pEmissFactor * vPowCombTurb[i])
@variable(m, 0 <= ePowGross[0:tHorz - 1])
@constraint(m, powGrossEq[i = 0:tHorz - 1], ePowGross[i] == vPowCombTurb[i] + ePowSteamTurb[i])
# 6
@constraint(m, powOutEq[i = 0:tHorz - 1], ePowOut[i] ==  ePowGross[i] - ePowUsePcc[i] - ePowUseComp[i])

# Steam
# 7a
@constraint(m, powHpEq[i = 0:tHorz - 1], ePowHp[i] == pFuelPowHp * eFuel[i])
# 7b
@constraint(m, powIpEq[i = 0:tHorz - 1], ePowIp[i] == pFuelPowIp * eFuel[i])

# 8
@constraint(m, stamIpLpEq[i = 0:tHorz - 1], eSteamIpLp[i] == pFuelSteamIpLp * eFuel[i])
# 9
@constraint(m, steamLp1[i = 0:tHorz - 1], eSteamLp1[i] == eSteamIpLp[i] - eSteamUsePcc[i])
# 10
@constraint(m, powLp1Eq[i = 0:tHorz - 1], ePowLp1[i] == eSteamLp1[i] / pHeatRateLp1)
# 11
@constraint(m, powLp2Eq[i = 0:tHorz - 1], ePowLp2[i] == pFuelPowLp * eFuel[i])
# 12
@constraint(m, powerSteamEq[i = 0:tHorz - 1], ePowSteamTurb[i] == ePowHp[i] + ePowIp[i] + ePowLp1[i] + ePowLp2[i])

# PCC
# 13
#@constraint(m, co2CapPccEq[i = 0:tHorz - 1], eCo2CapPcc[i] == pCo2PccCapRate * eCo2Fuel[i])
@constraint(m, co2CapPccEq[i = 0:tHorz - 1], eCo2CapPcc[i] == 0.95 * eCo2Fuel[i])
# 14
@constraint(m, co2PccOutEq[i = 0:tHorz - 1], eCo2PccOut[i] == eCo2Fuel[i] - eCo2CapPcc[i])
# 15
@constraint(m, co2DacFlueInEq[i = 0:tHorz - 1], eCo2DacFlueIn[i] == eCo2PccOut[i] - vCo2PccVent[i])
# 16
# @constraint(m, co2CapPccIn[i = 0:tHorz - 1], eCo2CapPcc[i] <= vCapPcc)
# 17
@constraint(m, steamUsePccEq[i = 0:tHorz - 1], eSteamUsePcc[i] == pSteamUseRatePcc * eCo2CapPcc[i])
# 18
@constraint(m, powerUsePccEq[i = 0:tHorz - 1], ePowUsePcc[i] == pPowUseRatePcc * eCo2CapPcc[i])

@constraint(m, co2DacFlueVentEq[i = 0:tHorz - 1], eCo2DacVent[i] == eCo2DacFlueIn[i])
# Co2 Compression
# 27
@constraint(m, co2CompEq[i = 0:tHorz - 1], eCo2Comp[i] == eCo2CapPcc[i])
# 28
@constraint(m, powUseCompEq[i = 0:tHorz - 1], ePowUseComp[i] == pPowUseRateComp * eCo2Comp[i])
# 29
# @constraint(m, powUseCompIn[i = 0:tHorz - 1], ePowUseComp[i] <= vCapComp)

@constraint(m, co2VentEq[i = 0:tHorz - 1], eCo2Vent[i] == vCo2PccVent[i] + eCo2DacVent[i])

@expression(m, eObjfExpr, sum(pCostFuel * eFuel[i] +
                              pEmissionPrice * eCo2Vent[i] +
                              pCo2TranspPrice * eCo2Comp[i] -
                              ePowOut[i] * 10 for i in 0:tHorz - 1))
#@expression(m, eObjfExpr, sum(-ePowOut[i] for i in 0:tHorz - 1))

@objective(m, Min, eObjfExpr)

set_optimizer(m, Clp.Optimizer)
set_optimizer_attribute(m, "LogLevel", 3)
set_optimizer_attribute(m, "PresolveType", 1)

# set_start_value(vCapCombTurb, 3.)
# set_start_value(vCapSteamTurb, 2.)
# set_start_value(vCapTransInter, 5.)
# set_start_value(vCapPcc, 5.)
# set_start_value(vCapComp, 8.)

# fix(vCapCombTurb, 3., force = true)
# fix(vCapSteamTurb, 2., force = true)
# fix(vCapTransInter, 5., force = true)
# fix(vCapPcc, 5., force = true)
# fix(vCapComp, 8., force = true)

# fix(a0[0], 0., force = true)
# fix(aR0[0], 0., force = true)
# fix(vAbs[0], 0., force = true)
# fix(vReg[0], 0., force = true)

#fix.(a0, 0., force = true)
#fix.(aR0, 0., force = true)
#fix.(vAbs, 0., force = true)
#fix.(vReg, 0., force = true)

optimize!(m)
println(termination_status(m))

#f = open("model.lp", "w")
#print(f, m)
#close(f)

write_to_file(m, "lp_mk0.mps")
#format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC

# Design decisions.
# vCapCombTurb
# vCapSteamTurb
# vCapTransInter
# vCapPcc
# vCapDac
# vCapComp
#
# Raw materials.
# eFuel
# sF0
# sS0
#
# penalties
# vCo2Vent
#
#
## vCo2PccVent * 2N

# Actual variables
# vPowCombTurb * 2

df = DataFrame(Symbol("eFuel[i]") => Float64[],
               Symbol("eCo2Fuel[i]") => Float64[],
               Symbol("eCo2CapPcc[i]") => Float64[],
               Symbol("eCo2Vent") => Float64[],
               Symbol("vCo2PccVent") => Float64[],
               Symbol("ePowUsePcc[i]") => Float64[],
               Symbol("ePowUseComp[i]") => Float64[],
               Symbol("demand") => Float64[],
               Symbol("ePowGross") => Float64[],
               Symbol("ePowOut") => Float64[])
for i in 0:tHorz-1
    push!(df, (value(eFuel[i]),
               value(eCo2Fuel[i]),
               value(eCo2CapPcc[i]),
               value(eCo2Vent[i]),
               value(vCo2PccVent[i]),
               value(ePowUsePcc[i]),
               value(ePowUseComp[i]),
               rnum[i+1],
               value(ePowGross[i]),
               value(ePowOut[i])))
end

print(df)

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
include("../../../good_plots/iron_line_heatmap.jl")

vars = LinRange(5.0, 85.0, 500)

hmp = iron_line_heatmap(tolerance=1e-12,
                        dtmax=5,
                        size_multiplier=5,
                        nbins=300,
                        vars=vars)
png(hmp, "iron_heatmap.png")