text:
mutable struct SimpleRingBuffer
    buffer::Vector{UInt8}
    offset::Int
end

SimpleRingBuffer(buf_size::Integer) = SimpleRingBuffer(Vector{UInt8}(undef, buf_size), 0)
Base.pointer(buf::SimpleRingBuffer) = pointer(buf.buffer) + buf.offset
reset!(buf::SimpleRingBuffer) = buf.offset = 0

function copy_data!(dest::SimpleRingBuffer, src::Memory, data_size::Integer)
    checkbounds(src, data_size)
    if length(dest.buffer) < data_size
        throw(ArgumentError(
            "Cannot store $data_size bytes in buffer of size $(length(dest.buffer))"
        ))
    end

    if dest.offset + data_size > length(dest.buffer)
        dest.offset = 0
    end
    data_start = pointer(dest)
    unsafe_copyto!(data_start, src.ptr, data_size)

    dest.offset += data_size
    return data_start
end

mutable struct LZ4HCCompressor <: TranscodingStreams.Codec
    streamptr::Ptr{LZ4_streamHC_t}
    compressionlevel::Cint
    block_size::Int
    buffer::SimpleRingBuffer
end

"""
    LZ4HCCompressor(; kwargs...)

Creates an LZ4 HC compression codec.

# Keywords
- `compressionlevel::Integer=$LZ4HC_CLEVEL_DEFAULT`: compression level
- `block_size::Integer=1024`: The size in bytes to encrypt into each block. (Max 4MB)
"""
function LZ4HCCompressor(;
    compressionlevel::Integer=LZ4HC_CLEVEL_DEFAULT,
    block_size::Integer=1024
)
    if block_size > LZ4_MAX_INPUT_SIZE
        throw(ArgumentError("`block_size` larger than $LZ4_MAX_INPUT_SIZE."))
    end

    return LZ4HCCompressor(
        C_NULL,
        compressionlevel,
        block_size,
        SimpleRingBuffer(4 * block_size),
    )
end

const LZ4HCCompressorStream{S} = TranscodingStream{LZ4HCCompressor,S} where S<:IO

"""
    LZ4HCCompressorStream(stream::IO; kwargs...)

Creates an LZ4 HC compression stream.
See `LZ4HCCompressorStream()` and `TranscodingStream()` for arguments.
"""
function LZ4HCCompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:compressionlevel, :block_size))
    return TranscodingStream(LZ4HCCompressor(; x...), stream; y...)
end

"""
   TranscodingStreams.initialize(codec::LZ4HCCompressor)

Initializes the LZ4 HC compression Codec.
"""
function TranscodingStreams.initialize(codec::LZ4HCCompressor)::Nothing
    codec.streamptr = LZ4_createStreamHC()
    nothing
end

"""
    TranscodingStreams.expectedsize(codec::Union{LZ4FastCompressor, LZ4HCCompressor}, input::Memory)

Returns the expected size of the transcoded data.
"""
function TranscodingStreams.expectedsize(codec::LZ4HCCompressor, input::Memory)::Int
    max_compressed_size(length(input), codec.block_size)
end

"""
   TranscodingStreams.minoutsize(codec::Union{LZ4FastCompressor, LZ4HCCompressor}, input::Memory)

Returns the minimum output size of `process`.
"""
function TranscodingStreams.minoutsize(codec::LZ4HCCompressor, input::Memory)::Int
    LZ4_compressBound(input.size) + CINT_SIZE
end

"""
    TranscodingStreams.finalize(codec::LZ4HCCompressor)

Finalizes the LZHC Compression Codec.
"""
function TranscodingStreams.finalize(codec::LZ4HCCompressor)::Nothing
    LZ4_freeStreamHC(codec.streamptr)
    codec.streamptr = C_NULL
    nothing
end

"""
    TranscodingStreams.startproc(codec::LZ4HCCompressor, mode::Symbol, error::Error)

Starts processing with the codec
"""
function TranscodingStreams.startproc(
    codec::LZ4HCCompressor,
    mode::Symbol,
    error::Error
)::Symbol

    try
        LZ4_resetStreamHC(codec.streamptr, codec.compressionlevel)
        reset!(codec.buffer)
        :ok
    catch err
        error[] = err
        :error
    end
end

"""
    TranscodingStreams.process(codec::LZ4HCCompressor, input::Memory, output::Memory, error::Error)

Compresses the data from `input` and writes to `output`.
The process follows the compression example from
https://github.com/lz4/lz4/blob/dev/examples/HCStreaming_ringBuffer.c
wherein each block of encoded data is prefixed by the number of bytes contained in that block.
Decoding can be done through the LZ4SafeDecompressor.
"""
function TranscodingStreams.process(
    codec::LZ4HCCompressor,
    input::Memory,
    output::Memory,
    error::Error
)::Tuple{Int,Int,Symbol}

    input.size == 0 && return (0, 0, :end)
    try
        data_size = min(input.size, codec.block_size)

        in_buffer = copy_data!(codec.buffer, input, data_size)
        out_buffer = Vector{UInt8}(undef, LZ4_compressBound(data_size))

        compressed_size = LZ4_compress_HC_continue(
            codec.streamptr,
            in_buffer,
            pointer(out_buffer),
            data_size,
            length(out_buffer)
        )

        checkbounds(output, compressed_size + CINT_SIZE)
        writeint(output, compressed_size)
        unsafe_copyto!(output.ptr+CINT_SIZE, pointer(out_buffer), compressed_size)

        return (data_size, compressed_size + CINT_SIZE, :ok)

    catch err
        error[] = err
        return (0, 0, :error)
    end
end


text:
using RNLA, Test, LinearAlgebra

@testset "rcur_test1" begin
    m = 100
    n = 50
    k = 4

    C = rand(m, k)
    U = rand(k, k)
    R = rand(k, n)
    A = C * U * R

    C1, U1, R1 = rcur(A)

    @test norm(C - C1) <= sqrt(eps())
    @test norm(U - U1) <= sqrt(eps())
    @test norm(R - R1) <= sqrt(eps())
end

@testset "rcur_test2" begin
end

ext: jl


text:
export cnAdjLoadChkpt, cnAdjCalcdRdu, cnAdjDirect

"""
function cnAdjLoadChkpt:
  loads a checkpoint from disk based on desired forward sweep index

  returns eqn object, which has only the necessary correct fields for the unsteady adjoint
"""
# function cnAdjLoadChkpt(mesh, sbp, opts, adj, physics_func, i_fwd, t)
function cnAdjLoadChkpt(mesh, sbp, opts, adj, physics_func, i_fwd)

  # initialize dummy eqn object for jacobian calculation use
  # Note: can't just call this:
  #   eqn_dummy = AdvectionData_{Tsol, Tres, Tdim, Tmsh}(mesh, sbp, opts)
  # because it needs AdvectionEquationMod loaded, which needs PDESolver loaded, which can't be loaded more than once.
  eqn_dummy = eqn_deepcopy(mesh, sbp, adj, opts)                     # allocate a dummy eqn object

  # println("in cnAdjLoadChkpt: pointer(adj.q_vec): ", pointer(adj.q_vec))
  # println("in cnAdjLoadChkpt: pointer(adj.res_vec): ", pointer(adj.res_vec))
  # println("in cnAdjLoadChkpt: pointer(eqn_dummy.q_vec):   ", pointer(eqn_dummy.q_vec))
  # println("in cnAdjLoadChkpt: pointer(eqn_dummy.res_vec): ", pointer(eqn_dummy.res_vec))

  qvec_filename = string("qvec_for_adj-", i_fwd, ".dat")
  println("  Calculating Jac using forward sweep data from: ", qvec_filename)
  q_vec_with_complex = readdlm(qvec_filename)
  eqn_dummy.q_vec = q_vec_with_complex[:,1]     # because readdlm gives a field to the zero-valued complex part
  array1DTo3D(mesh, sbp, eqn_dummy, opts, eqn_dummy.q_vec, eqn_dummy.q)

  check_q_qvec_consistency(mesh, sbp, eqn_dummy, opts)
  println("  eqn_dummy.q_vec loaded. i_fwd: ", i_fwd)
  # print_qvec_coords(mesh, sbp, eqn_dummy, opts)

  vis_filename = string("solution_loadedfromdisk_ifwd-", i_fwd)
  saveSolutionToMesh(mesh, real(eqn_dummy.q_vec))
  writeVisFiles(mesh, vis_filename)

  # sync up eqn_dummy.q and eqn_dummy.q_vec

  # TODO: is this needed here? YES.
  #     explain why here
  # TODO new: I don't think so now that eqn_deepcopy is properly implemented
  # eqn_dummy.q = reshape(eqn_dummy.q_vec, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
  # eqn_dummy.res = reshape(eqn_dummy.res_vec, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)

  #=
  # make eqn object consistent
  evalResidual(mesh, sbp, eqn_dummy, opts)
  array3DTo1D(mesh, sbp, eqn_dummy, opts, eqn_dummy.res, eqn_dummy.res_vec)

  # TODO: fix ordering & repetition of these
  eqn_dummy.q = reshape(eqn_dummy.q_vec, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
  eqn_dummy.res = reshape(eqn_dummy.res_vec, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
  =#
  # looks like the above commented section doesn't make a difference

  # use startDataExchange to sync up q/q_vec and q_face_send/recv
  if opts["parallel_type"] == 2 && mesh.npeers > 0
    startDataExchange(mesh, opts, eqn_dummy.q, eqn_dummy.q_face_send, eqn_dummy.q_face_recv, eqn_dummy.params.f)
  end

  return eqn_dummy

end   # end function cnAdjLoadChkpt

"""
function cnAdjCalcdRdu:
  calculates the Jacobian dRdu at an arbitrary time t given an associated eqn_dummy object

  returns a jac
"""
function cnAdjCalcdRdu(mesh, sbp, opts, eqn_dummy, physics_func, i_fwd, t)
  # Note: probably don't need to allocate another jac, but this should cause no problems aside from allocation cost
  newton_data_discard, jac, rhs_vec_discard = setupNewton(mesh, mesh, sbp, eqn_dummy, opts, physics_func, () )

  # make sure we're doing complex step! since res_copy is zeros, it would mess up the FD calc
  assert(opts["jac_method"] == 2)
  epsilon = opts["epsilon"]
  pert = complex(0, epsilon)

  # checked: eqn_dummy is loaded properly
  # checked: cjc args are correct
  # checked: norm of eqn_dummy.q_vec in cJC matches it here

  println("  calculating jac for loaded eqn_dummy. i_fwd = ", i_fwd,", t = ", t)
  calcJacobianComplex(newton_data_discard, mesh, sbp, eqn_dummy, opts, physics_func, pert, jac, t)

  println(" ::::::::::::: in cnAdjCalcdRdu, after calcJacobianComplex, norm(jac): ", norm(jac, 1))

  # filename = string("jac_from_cnAdjCalcdRdu-i_fwd-",i_fwd,".dat")
  # writedlm(filename, jac)

  return jac
end   # end function cnAdjCalcdRdu

"""
function cnAdjDirect:
  performs a direct solve of the linear adjoint PDE

  (*)_(i+1) corresponds to t
  (*)_i corresponds to t_nextstep

  given psi_(i+1)

  returns psi_i
"""
function cnAdjDirect(mesh, sbp, opts, adj, physics_func, jac, i_fwd, h, t_steps, t, dRdu_global_rev)
  # Note: t is passed in (instead of t_nextstep), but t_nextstep is calc'd & used below
  #
  # Note: adj is not altered in this function. It is only used for passing into cnAdjLoadChkpt for an eqn_deepcopy

# adj_nextstep.q_vec = cnAdjDirect(mesh, sbp, opts, adj_nextstep, jac, i_fwd, t_nextstep)

  t_nextstep = t - h
  t_nextstep = negativeZeroCheck(t_nextstep)   # ensure negative zero is changed to zero

  # Load checkpoint at time t, aka i+1
  eqn_dummy = cnAdjLoadChkpt(mesh, sbp, opts, adj, physics_func, i_fwd)   # adj is only used for a eqn_deepcopy within
                                                                          #   cnAdjLoadCheckpt
  # checked: eqn_dummy is loaded properly
  # Use time t's eqn.q_vec and t_nextstep, aka i, to compute the jac (dRdu) at t_nextstep aka i
  jac = cnAdjCalcdRdu(mesh, sbp, opts, eqn_dummy, physics_func, i_fwd, t_nextstep)
  dRdu_i = transpose(jac)       # see derivation for transpose justification

  # TODO: double check that there is not an off by one error on dRdu:
  #       even though calculated at t_nextstep, am I loading q_vec at i+1 instead of i?

  I = eye(length(adj.q_vec))

  if opts["uadj_global"]
    ### GLOBAL dRdu section
    B1_pos_for_dRdu_global =  1.0*I - 0.5*h*dRdu_i
    B1_neg_for_dRdu_global = -1.0*I - 0.5*h*dRdu_i

    blksz = mesh.numDof
    # blksz = 3   # testing 44
    # B1_pos_for_dRdu_global = ones(blksz, blksz)*44
    # B1_neg_for_dRdu_global = ones(blksz, blksz)*55

    println(" GLOBAL: forming global dRdu (rev), i_fwd = $i_fwd")
    println(" GLOBAL: size of B1_pos_for_dRdu_global: ", size(B1_pos_for_dRdu_global))
    row_ix_start = (i_fwd-2)*blksz + 1
    row_ix_end = (i_fwd-2)*blksz + blksz
    col_ix_start = (i_fwd-2)*blksz + 1
    col_ix_end = (i_fwd-2)*blksz + blksz
    println(" GLOBAL: diagonal block ix's: [",row_ix_start,":",row_ix_end,", ",col_ix_start,":", col_ix_end,"]")
    dRdu_global_rev[row_ix_start:row_ix_end, col_ix_start:col_ix_end] = B1_pos_for_dRdu_global

    if i_fwd != t_steps + 1
      col_ix_start = (i_fwd-2)*blksz + blksz + 1
      col_ix_end = (i_fwd-2)*blksz + 2*blksz
      row_ix_start = (i_fwd-2)*blksz + 1
      row_ix_end = (i_fwd-2)*blksz + blksz
      println(" GLOBAL: off-diagonal block ix's: [",row_ix_start,":",row_ix_end,", ",col_ix_start,":", col_ix_end,"]")

      dRdu_global_rev[row_ix_start:row_ix_end, col_ix_start:col_ix_end] = B1_neg_for_dRdu_global
    end   # end of "if i_fwd != t_steps + 1"
    ### end GLOBAL
  end # end of "if opts["uadj_global"]

  dJdu_i = calcdJdu_CS(mesh, sbp, eqn_dummy, opts, h, t)  # obtain dJdu at time step i

  # 20170831. Looks like there's a transpose missing here on B1.
  # See notes, LHS, 08/30/2017, for rigorous explanation.
  # B1 = I - 0.5*h*dRdu_i
  # B2 = adj.q_vec + 0.5*h*dRdu_i*adj.q_vec - dJdu_i
  B1 = transpose(I - 0.5*h*dRdu_i)
  B2 = adj.q_vec + 0.5*h*dRdu_i*adj.q_vec - dJdu_i

  nextstep_q_vec = zeros(adj.q_vec)
  nextstep_q_vec = B1\B2

  # filename = string("eqn_dummyqvec_from_cnAdjDirect-i_fwd-",i_fwd,".dat")
  # writedlm(filename, eqn_dummy.q_vec)

  return nextstep_q_vec, jac

end     # end function cnAdjDirect


text:
using AbstractAlgebra
using Base.Threads

# import Singular

include("../src/Groebner.jl")

function generatesystems(ground)
    [
        # ("noon 3", Groebner.noon3(ground=ground)),
        ("cyclic 4", Groebner.cyclicn(4, ground=ground)),
        ("cyclic 5", Groebner.cyclicn(5, ground=ground)),
        ("cyclic 6", Groebner.cyclicn(6, ground=ground)),
        ("cyclic 7", Groebner.cyclicn(7, ground=ground)),
        ("cyclic 8", Groebner.cyclicn(8, ground=ground)),
        ("cyclic 9", Groebner.cyclicn(9, ground=ground)),
        # ("eco 11", Groebner.eco11(ground=ground)),
        # ("eco 12", Groebner.eco12(ground=ground)),
        # ("eco 13", Groebner.eco13(ground=ground))
        #("noon 6", Groebner.noonn(6, ground=ground)),
        #("noon 7", Groebner.noonn(7, ground=ground)),
        #("noon 8", Groebner.noonn(8, ground=ground)),
        # ("noon 9", Groebner.noonn(9, ground=ground)),
        #("cyclic 10", Groebner.rootn(10, ground=ground)),
        #("cyclic 11", Groebner.rootn(11, ground=ground)),
        #("cyclic 12", Groebner.rootn(12, ground=ground)),
        #("henrion 5", Groebner.henrion5(ground=ground)),
        #("henrion 6", Groebner.henrion6(ground=ground)),
        # ("henrion 7", Groebner.henrion7(ground=ground))
    ]
end

function iszerodim(gb)
    R = parent(first(gb))
    vs = gens(R)
    for v in vs
        flag = false
        for lt in map(leading_monomial, gb)
            if exponent_vector(v, 1) .* total_degree(lt) == exponent_vector(lt, 1)
                flag = true
            end
        end
        if !flag
            return false
        end
    end
    return true
end

function nsols(gb)
    R = parent(first(gb))
    vs = gens(R)
    sols = 1
    for v in vs
        for lt in map(leading_monomial, gb)
            if exponent_vector(v, 1) .* total_degree(lt) == exponent_vector(lt, 1)
                sols *= total_degree(lt)
            end
        end
    end
    return sols
end

#=
function systeminfo(system)
    R = AbstractAlgebra.parent(system[1])
    modulo = AbstractAlgebra.characteristic(R)
    n = AbstractAlgebra.nvars(R)
    ground_s = Singular.N_ZpField(modulo)
    R_s, _ = Singular.PolynomialRing(ground_s, ["x$i" for i in 1:n], ordering=:degrevlex)

    system_s = map(
        f -> AbstractAlgebra.change_base_ring(
                    ground_s,
                    AbstractAlgebra.map_coefficients(c -> ground_s(c.d), f),
                    parent=R_s),
        system)

    ideal_s = Singular.Ideal(R_s, system_s)
    gb = Singular.std(ideal_s)

    # println("Is zerodim ", Singular.is_zerodim(ideal_s))
    println("dim R / I ", length(Singular.gens(Singular.kbase(gb))))
end
=#

function systeminfo(system)
    println(iszerodim(Groebner.groebner(system, ordering=:degrevlex)))
end

function runall(ground)

    println("-"^20)
    for (name, system) in generatesystems(ground)
        println(name)
        systeminfo(system)
        println("-"^20)
    end
end

ground = GF(2^31-1)
runall(ground)


text:
#=
1D Euler equations
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVeuler1d");

useLog("FVeuler1dlog", level=3)

# Configuration setup
domain(1) # 1D
solverType(FV)
timeStepper(RK4)
timeInterval(0.2)

# Mesh
ord = 1 # order used for flux
n = Int(ceil(150/ord)) # number of elements scaled for subdividing for higher orders
mesh(LINEMESH, elsperdim=n) # A uniform 1D mesh of the unit interval.

if ord > 1
    finiteVolumeOrder(ord)
end

# Primitive variables
r = variable("r", location=CELL)
u = variable("u", location=CELL)
p = variable("p", location=CELL)
U = [r, u, p];

# Conserved variables
q1 = variable("q1", location=CELL)
q2 = variable("q2", location=CELL)
q3 = variable("q3", location=CELL)
Q = [q1, q2, q3];

# Transformations between the variable sets
U2Q = variableTransform(U, Q,
(U) -> (
    gamma = 1.4;
    Q = [U[1],
         U[1] * U[2],
         U[3] / (gamma-1) + 0.5 * U[1] * U[2] * U[2]];
    return Q;
))
Q2U = variableTransform(Q, U,
(Q) -> (
    gamma = 1.4;
    U = [Q[1],
         Q[2] / Q[1],
         (gamma-1) * (Q[3] - 0.5 * Q[2] * Q[2] / Q[1])];
    return U;
))

# Shock tube initial conditions
initial(r, "x<0.3 ? 1 : 0.125")
initial(u, "x<0.3 ? 0.75 : 0")
initial(p, "x<0.3 ? 1 : 0.1")
Finch.eval_initial_conditions(); # set U initial conditions
transformVariable(U2Q); # set Q initial conditions

# Boundary conditions
boundary(q1, 1, NO_BC) # equivalent to transmissive bdry
boundary(q2, 1, NO_BC)
boundary(q3, 1, NO_BC)

# The commented function is a compact version for inspection.
# The longer function below that is similar, but arranged to work
# as a callback in the flux function.
#=
function flux_u(U)
    local r = U[1];
    local u = U[2];
    local p = U[3];
    local gamma = 1.4;

    FU = [ r * u,
           p + r * u*u,
           u * (p + p/(gamma-1) + r*u*u*0.5)];

    # Steger-Warming splitting
    c = sqrt(gamma*p/r)
    M = u/c
    if M < -1
        Fm = FU;
        Fp = zeros(3);
    elseif M > 1
        Fm = zeros(3);
        Fp = FU;
    elseif M >= -1 && M < 0
        Fp = r*(u+c)/(2*gamma) .* [1, (u+c), (u+c)^2/2+(3-gamma)/(gamma-1)*c^2/2];
        Fm = FU - Fp;
    else # 0 < M < 1
        Fm = r*(u-c)/(2*gamma) .* [1, (u-c), (u-c)^2/2+(3-gamma)/(gamma-1)*c^2/2];
        Fp = FU - Fm;
    end

    return (Fm, Fp);
end
=#
function flux_u2(r1,u1,p1, r2,u2,p2, comp, normal)
    gamma = 1.4;
    # The normal points from cell 1 to cell 2
    # So "left" is currently cell 1
    if normal < 0
        # swap left and right
        tmpr = r1; tmpu = u1; tmpp = p1;
        r1 = r2;   u1 = u2;   p1 = p2;
        r2 = tmpr; u2 = tmpu; p2 = tmpp;
    end
    c1 = sqrt(gamma*p1/r1);
    M1 = u1/c1;
    c2 = sqrt(gamma*p2/r2);
    M2 = u2/c2;

    if comp == 1 # rho
        if M1 < -1
            Fp = 0;
        elseif M1 > 1
            Fp = r1 * u1;
        elseif M1 >= -1 && M1 < 0
            Fp = r1*(u1+c1)/(2*gamma);
        else # 0 < M < 1
            Fp = r1*u1 - r1*(u1-c1)/(2*gamma);
        end

        if M2 < -1
            Fm = r2*u2;
        elseif M2 > 1
            Fm = 0;
        elseif M2 >= -1 && M2 < 0
            Fm = r2*u2 - r2*(u2+c2)/(2*gamma);
        else # 0 < M < 1
            Fm = r2*(u2-c2)/(2*gamma);
        end

    elseif comp == 2 # u
        if M1 < -1
            Fp = 0;
        elseif M1 > 1
            Fp = p1 + r1 * u1*u1;
        elseif M1 >= -1 && M1 < 0
            Fp = r1*(u1+c1)/(2*gamma) * (u1+c1);
        else # 0 < M < 1
            Fp = p1 + r1 * u1*u1 - r1*(u1-c1)/(2*gamma) * (u1-c1);
        end

        if M2 < -1
            Fm = p2 + r2 * u2*u2;
        elseif M2 > 1
            Fm = 0;
        elseif M2 >= -1 && M2 < 0
            Fm = p2 + r2 * u2*u2 - r2*(u2+c2)/(2*gamma) * (u2+c2);
        else # 0 < M < 1
            Fm = r2*(u2-c2)/(2*gamma) * (u2-c2);
        end

    elseif comp == 3 # p
        if M1 < -1
            Fp = 0;
        elseif M1 > 1
            Fp = u1 * (p1 + p1/(gamma-1) + r1*u1*u1*0.5);
        elseif M1 >= -1 && M1 < 0
            Fp = r1*(u1+c1)/(2*gamma) * ((u1+c1)^2/2+(3-gamma)/(gamma-1)*c1^2/2);
        else # 0 < M < 1
            Fp = u1 * (p1 + p1/(gamma-1) + r1*u1*u1*0.5) - (r1*(u1-c1)/(2*gamma) * ((u1-c1)^2/2+(3-gamma)/(gamma-1)*c1^2/2));
        end

        if M2 < -1
            Fm = u2 * (p2 + p2/(gamma-1) + r2*u2*u2*0.5);
        elseif M2 > 1
            Fm = 0;
        elseif M2 >= -1 && M2 < 0
            Fm = u2 * (p2 + p2/(gamma-1) + r2*u2*u2*0.5) - (r2*(u2+c2)/(2*gamma) * ((u2+c2)^2/2+(3-gamma)/(gamma-1)*c2^2/2));
        else # 0 < M < 1
            Fm = r2*(u2-c2)/(2*gamma) * ((u2-c2)^2/2+(3-gamma)/(gamma-1)*c2^2/2);
        end
    end
    return (Fp + Fm) * normal;
end

callbackFunction(flux_u2);

flux(Q, ["flux_u2(left(r),left(u),left(p), right(r),right(u),right(p), $i, normal())" for i in 1:3])
#=
The line above is in a compact form. A verbose form would be:
flux([q1,q2,q3], ["flux_u2(left(r),left(u),left(p), right(r),right(u),right(p), 1, normal())"
                  "flux_u2(left(r),left(u),left(p), right(r),right(u),right(p), 2, normal())"
                  "flux_u2(left(r),left(u),left(p), right(r),right(u),right(p), 3, normal())"])

Note:
To get a reconstructed value for a variable r,
    "r"         -> A centered interpolation
    "left(r)"   -> extrapolated from cells on side 1 of the face (the normal vector points from side 1 to side 2)
    "right(r)"  -> similar, but for side 2
    "central(r)"-> similar to "r"
In all cases, neighborhoods that are truncated by a boundary will simply use fewer cells.
Boundary faces will use one ghost cell identical to the interior cell.
The side1/2 orientation is only tied to the normal vector, so if the physical orientation is needed,
pass the normal vector as well with "normal()".
=#

# For inspecting the generated code, look in "fveuler1dcode.jl"
# If you wish to modify it and reimport, use importCode("fveuler1dcode") instead.
# Importing will override anything previously generated.
exportCode("fveuler1dcode")

# Need to update primitive variables each time conserved variables are changed.
@postStepFunction(
    transformVariable(Q2U)
);

# Solve for the conserved variables because the equations were defined for them.
# Since primitive variables are also updated each step, they are also solved for.
solve(Q)

finalize_finch()

##### Uncomment below to plot

x = Finch.fv_info.cellCenters[:];

using Plots
pyplot();
display(plot([x x x], [r.values[:] u.values[:] p.values[:]], markershape=:circle, label=["r" "u" "p"]))

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
function isisogram(s)

end


text:
# ## Contour plot
# ### Contour plot

#md # ![FILE_NAME.png](images/FILE_NAME.png)

f(x,y) = x^2 + y^2
x = repeat(range(0, stop=2, length=50), 1, 50)
y = repeat(range(0, stop=2, length=50), 1, 50)'

contour( f.(x, y) )


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DelayEmbeddings
using DynamicalSystemsBase
using Test
using DelimitedFiles
using Random
using PyPlot
pygui(true)

## We here focus on the dependence of the L-statistic on its parameter Tw and
# exemplary look at the Lorenz63 and the Lorenz96 system, as well as on random
# uniformly distributed noise.

## Generate Data:
# Lorenz63

# The time series have been computed in the script `/test/timeseries/produce_timeseries.jl`
s = readdlm("./test/timeseries/lorenz_pecora_uni_x.csv")
s1 = vec(s[1:5000]) # input timeseries = x component of lorenz
w1 = estimate_delay(vec(s), "mi_min")    # 1st minimum of mutual inf

# Lorenz96
idx = 487
N = 8 # number of oscillators
Fs = 3.5:0.002:5 # parameter spectrum
F = Fs[idx]
dt = 0.1 # sampling time
total = 5000  # time series length

# Parameters analysis:
ε = 0.05  # recurrence threshold
dmax = 10   # maximum dimension for traditional tde
lmin = 2   # minimum line length for RQA
trials = 80 # trials for MCDTS
taus = 0:100 # possible delays

#pick one time series
t_idx = 2

# init Lorenz96
u0 = [0.590; 0.766; 0.566; 0.460; 0.794; 0.854; 0.200; 0.298]
lo96 = Systems.lorenz96(N, u0; F = F)

#set_parameter!(lo96, 1, F)
data = trajectory(lo96, total*dt; dt = dt, Ttr = 2500 * dt)
s2 = data[1:5000,t_idx]
w2 = estimate_delay(vec(s2), "mi_min")    # 1st minimum of mutual inf

# Random numbers
Random.seed!(1234)
s3 = rand(1000)
w3 = estimate_delay(vec(s3), "mi_min")    # 1st minimum of mutual inf

# normalize data
D1 = regularize(Dataset(s1))
D2 = regularize(Dataset(s2))
D3 = regularize(Dataset(s3))

# generate embedding candidate
Y_t1 = DelayEmbeddings.hcat_lagged_values(D1, vec(Matrix(D1)), w1)
Y_t2 = DelayEmbeddings.hcat_lagged_values(D2, vec(Matrix(D2)), w2)
Y_t3 = DelayEmbeddings.hcat_lagged_values(D3, vec(Matrix(D3)), w3)

Tw = 100    # maximum Tw
L11 = zeros(Tw)
L12 = zeros(Tw)
L21 = zeros(Tw)
L22 = zeros(Tw)
L31 = zeros(Tw)
L32 = zeros(Tw)
for T = 1:Tw
    L11[T] = uzal_cost(D1; w = w1, samplesize = 1, Tw = T)
    L12[T] = uzal_cost(Y_t1; w = w1, samplesize = 1, Tw = T)
    L21[T] = uzal_cost(D2; w = w2, samplesize = 1, Tw = T)
    L22[T] = uzal_cost(Y_t2; w = w2, samplesize = 1, Tw = T)
    L31[T] = uzal_cost(D3; w = w3, samplesize = 1, Tw = T)
    L32[T] = uzal_cost(Y_t3; w = w3, samplesize = 1, Tw = T)
end

# Plotting

panelnames = ["A" "B" "C"]

fig = figure(figsize=[15,5])

lwg = 1.5         # linewidth of the graph
lwa = 2         # linewidth of the axis
fsp = 18        # Fontsize of panelnames
fsa = 12        # Fontsize of the axis
fsl = 10        # Fontsize of the legendentries
fst = 20        # Fontsize of title
axislabelsize = 12 # axislabelsize
ticklabelsize = 10  # labelsize of ticks
ms = 20         # markersize of maxima
ms_style = "*"  # markerstyle

subplots_adjust(
#left = 0.07,
#right = 0.96,
#bottom = 0.05,
#top = 0.96,
wspace = 0.35,
#hspace = 0.55)
)

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
font0 = Dict(
        "font.size" => fsa,
        "font.weight" => "bold",
        "axes.labelweight" => "bold",
        "axes.titleweight" => "bold",
        "axes.labelsize" => axislabelsize,
        "axes.linewidth" => lwa,
        "xtick.labelsize" => ticklabelsize,
        "ytick.labelsize" => ticklabelsize,
        "legend.fontsize" => fsl)
        #"text.usetex" => true)
merge!(rcParams, font0)

ax = fig.add_subplot(1, 3, 1)
p1 = plot(1:Tw, L11, label="time series", linewidth=lwg)
color1 = p1[1].get_color()
scatter(1:Tw, L11, ms, label="")
p2 = plot(1:Tw, L12, label="2-dim embedding", linewidth=lwg)
color2 = p2[1].get_color()
scatter(1:Tw, L12, ms, label="")
ylims = ax.get_ylim()
ax.fill_between(1:Tw, L11, ylims[1], where = L11 .>= L12, facecolor = "gray", alpha = 0.5)
legend()
ax.text(-0.2, 1.1, panelnames[1], transform=ax.transAxes,
     fontsize=fsp, fontweight="bold", va="top")
title("Lorenz 63")
xlabel(L"T_M")
ylabel("L")
grid()

ax = fig.add_subplot(1, 3, 2)
p1 = plot(1:Tw, L21, label="time series", linewidth=lwg)
color1 = p1[1].get_color()
scatter(1:Tw, L21, ms, label="")
p2 = plot(1:Tw, L22, label="2-dim embedding", linewidth=lwg)
color2 = p2[1].get_color()
scatter(1:Tw, L22, ms, label="")
ylims = ax.get_ylim()
ax.fill_between(1:Tw, L21, ylims[1], where = L21 .>= L22, facecolor = "gray", alpha = 0.5)
#legend()
ax.text(-0.2, 1.1, panelnames[2], transform=ax.transAxes,
     fontsize=fsp, fontweight="bold", va="top")
title("Lorenz 96")
xlabel(L"T_M")
grid()

ax = fig.add_subplot(1, 3, 3)
p1 = plot(1:Tw, L31, label="time series", linewidth=lwg)
color1 = p1[1].get_color()
scatter(1:Tw, L31, ms, label="")
p2 = plot(1:Tw, L32, label="2-dim embedding", linewidth=lwg)
color2 = p2[1].get_color()
scatter(1:Tw, L32, ms, label="")
ylims = ax.get_ylim()
ax.fill_between(1:Tw, L31, ylims[1], where=L31 .>= L32, facecolor = "gray", alpha = 0.5)
#legend()
ax.text(-0.2, 1.1, panelnames[3], transform=ax.transAxes,
     fontsize=fsp, fontweight="bold", va="top")
title("Random numbers")
xlabel(L"T_M")
grid()


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -0.27492067550218746 <= q <= 2.084201684402027)

                     add_NL_constraint(m, :(softplus(-0.0416569430659357 + 0.13675627334278007*softplus(0.1752604319892206 + -0.8557289632032381*softplus(-0.028885055816546323 + 0.15091441453709065*softplus(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*softplus(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.21149394336428795*softplus(-0.29077418744555983 + -0.11119024401630329*softplus(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*softplus(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3])))) + 0.703085236566185*softplus(-0.4914467550189556 + 0.5078124686195609*softplus(-0.028885055816546323 + 0.15091441453709065*softplus(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*softplus(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.7523099505151869*softplus(-0.29077418744555983 + -0.11119024401630329*softplus(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*softplus(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) + softplus(0.6991595157542867 + 0.5132010992296534*softplus(0.1752604319892206 + -0.8557289632032381*softplus(-0.028885055816546323 + 0.15091441453709065*softplus(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*softplus(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.21149394336428795*softplus(-0.29077418744555983 + -0.11119024401630329*softplus(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*softplus(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3])))) + 0.6682334830395109*softplus(-0.4914467550189556 + 0.5078124686195609*softplus(-0.028885055816546323 + 0.15091441453709065*softplus(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*softplus(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.7523099505151869*softplus(-0.29077418744555983 + -0.11119024401630329*softplus(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*softplus(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Autogenerated wrapper script for Birch_Standard_jll for x86_64-apple-darwin
export libbirch_standard, libbirch_standard_debug

using boost_jll
using CompilerSupportLibraries_jll
using Eigen_jll
using LibBirch_jll
using LibYAML_jll
JLLWrappers.@generate_wrapper_header("Birch_Standard")
JLLWrappers.@declare_library_product(libbirch_standard, "@rpath/libbirch-standard.0.dylib")
JLLWrappers.@declare_library_product(libbirch_standard_debug, "@rpath/libbirch-standard-debug.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, CompilerSupportLibraries_jll, Eigen_jll, LibBirch_jll, LibYAML_jll)
    JLLWrappers.@init_library_product(
        libbirch_standard,
        "lib/libbirch-standard.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbirch_standard_debug,
        "lib/libbirch-standard-debug.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
Source:
  S. Ahmed and R. Garcia. "Dynamic Capacity Acquisition and Assignment under Uncertainty," Annals of Operations Research, vol.124, pp. 267-283, 2003

Input:
  nR: number of resources
  nN: number of tasks
  nT: number of time periods
  nS: number of scenarios

Sets:
  sR: resources
  sN: tasks
  sT: time periods
  sS: scenarios

Variables (1st Stage):
  x[i,t]: capacity acquired for resource i at period t
  u[i,t]: 1 if x[i,t] > 0, 0 otherwise

Variables (2nd Stage):
  y[i,j,t]: 1 if resource i is assigned to task j in period t, 0 otherwise

Parameters (general):
  a[i,t]: linear component of expansion cost for resource i at period t
  b[i,t]: fixed component of expansion cost for resource i at period t
  c[i,j,t,s]: cost of assigning resource i to task j in period t
  c0[j,t,s]: penalty incurred if task j in period t is not served

Parameters (scenario):
  d[j,t,s]: capacity required for to perform task j in period t in scenario s
=#
using DualDecomposition
using JuMP, Ipopt, GLPK
using Random

const DD = DualDecomposition

function main_dcap(nR::Int, nN::Int, nT::Int, nS::Int, seed::Int=1)

  Random.seed!(seed)

  sR = 1:nR
  sN = 1:nN
  sT = 1:nT
  sS = 1:nS

  ## parameters
  a = rand(nR, nT) * 5 .+ 5
  b = rand(nR, nT) * 40 .+ 10
  c = rand(nR, nN, nT, nS) * 5 .+ 5
  c0 = rand(nN, nT, nS) * 500 .+ 500
  d = rand(nN, nT, nS) .+ 0.5
  Pr = ones(nS)/nS

  # This creates a Lagrange dual problem for each scenario s.
  function create_scenario_model(s::Int64)

      # construct JuMP.Model
      model = Model(GLPK.Optimizer)

      ## 1st stage
      @variable(model, x[i=sR,t=sT] >= 0)
      @variable(model, u[i=sR,t=sT], Bin)
      @variable(model, y[i=sR,j=sN,t=sT], Bin)
      @variable(model, z[j=sN,t=sT], Bin)
      @objective(model, Min,
            sum(a[i,t]*x[i,t] + b[i,t]*u[i,t] for i in sR for t in sT)
          + sum(c[i,j,t,s]*y[i,j,t] for i in sR for j in sN for t in sT)
          + sum(c0[j,t,s]*z[j,t] for j in sN for t in sT))
      @constraint(model, [i=sR,t=sT], x[i,t] - u[i,t] <= 0)
      @constraint(model, [i=sR,t=sT], -sum(x[i,tau] for tau in 1:t) + sum(d[j,t,s]*y[i,j,t] for j in sN) <= 0)
      @constraint(model, [j=sN,t=sT], sum(y[i,j,t] for i in sR) + z[j,t] == 1)

      return model
  end

  # Create DualDecomposition instance.
  algo = DD.LagrangeDual()

  # Add Lagrange dual problem for each scenario s.
  models = Dict{Int,JuMP.Model}(s => create_scenario_model(s) for s in sS)
  for s in sS
      DD.add_block_model!(algo, s, models[s])
  end

  coupling_variables = Vector{DD.CouplingVariableRef}()
  for s in sS
      model = models[s]
      xref = model[:x]
      for i in sR, t in sT
          push!(coupling_variables, DD.CouplingVariableRef(s, (i,t), xref[i,t]))
      end
      uref = model[:u]
      for i in sR, t in sT
          push!(coupling_variables, DD.CouplingVariableRef(s, (i,t), uref[i,t]))
      end
  end

  # Set nonanticipativity variables as an array of symbols.
  DD.set_coupling_variables!(algo, coupling_variables)

  # Solve the problem with the solver; this solver is for the underlying bundle method.
  DD.run!(algo, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
end

main_dcap(2,3,3,20)


text:
module HMatrices

const PROJECT_ROOT =  pkgdir(HMatrices)

using StaticArrays
using LinearAlgebra
using Statistics:median
using TimerOutputs
using Printf
using RecipesBase
using Distributed
using Base.Threads

using WavePropBase
using WavePropBase.Trees
using WavePropBase.Utils
using WavePropBase.Geometry

import AbstractTrees
import LinearAlgebra: mul!, lu!, lu, LU, ldiv!, rdiv!, axpy!, rank, rmul!, lmul!
import Base: Matrix, adjoint

"""
    const ALLOW_GETINDEX

If set to false, the `getindex(H,i,j)` method will throw an error on
`AbstractHMatrix`.
"""
const ALLOW_GETINDEX = Ref(true)

include("utils.jl")
include("kernelmatrix.jl")
include("rkmatrix.jl")
include("compressor.jl")
include("hmatrix.jl")
include("dhmatrix.jl")
include("addition.jl")
include("multiplication.jl")
include("triangular.jl")
include("lu.jl")

export Utils,
    Geometry,
    Trees,
    # types (re-exported)
    ClusterTree,
    CardinalitySplitter,
    DyadicSplitter,
    GeometricSplitter,
    GeometricMinimalSplitter,
    HyperRectangle,
    # abstract types
    AbstractKernelMatrix,
    # types
    HMatrix,
    KernelMatrix,
    StrongAdmissibilityStd,
    WeakAdmissibilityStd,
    PartialACA,
    ACA,
    TSVD,
    # functions
    compression_ratio,
    print_tree,
    assemble_hmat,
    # macros
    @hprofile

end

max_stars_repo_name: WaveProp/HMatrices

max_issues_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

__precompile__(true)

module SuiteSparse

import Base: \
import LinearAlgebra: ldiv!, rdiv!

## Functions to switch to 0-based indexing to call external sparse solvers

# Convert from 1-based to 0-based indices
function decrement!(A::AbstractArray{T}) where T<:Integer
    for i in 1:length(A); A[i] -= oneunit(T) end
    A
end
decrement(A::AbstractArray{<:Integer}) = decrement!(copy(A))

# Convert from 0-based to 1-based indices
function increment!(A::AbstractArray{T}) where T<:Integer
    for i in 1:length(A); A[i] += oneunit(T) end
    A
end
increment(A::AbstractArray{<:Integer}) = increment!(copy(A))

if Base.USE_GPL_LIBS
    include("umfpack.jl")
    include("cholmod.jl")
    include("spqr.jl")
    include("deprecated.jl")
end

end # module SuiteSparse


text:
function nlsolve!(f, ∇f!, scheme, x0; inverse=true, c=0.001, g_tol=1e-8, maxiter=10^4, show_trace=false)
    n = length(x0)
    x_curr = copy(x0)
    ∇f_next, ∇f_curr = similar(x0), similar(x0)
    ∇f!(∇f_curr, x_curr)

    B = initial_approximation(inverse, scheme, n)

    d = .-solve(B, ∇f_curr)
    α, f_α, ls_success = backtracking(f, x_curr, d, ∇f_curr; α_0=1.0, c=c)
    Δx = α .* d
    x_next = x_curr .+ Δx

    ∇f!(∇f_next, x_next)
    y = ∇f_next .- ∇f_curr
    update!(scheme, B, Δx, y)

    for i = 1:maxiter
        copyto!(x_curr, x_next)
        ∇f!(∇f_curr, x_curr)

        d = solve!(d, B, ∇f_curr)

        α, f_α, ls_success = backtracking(f, x_curr, d, ∇f_curr; α_0=1.0, c=c, verbose=show_trace)
        Δx .= α .* d
        x_next .= x_curr .+ Δx
        ∇f!(∇f_next, x_next)
        y .= ∇f_next .- ∇f_curr

        if norm(∇f_next) < g_tol
            return x_next, ∇f_next
        end

        update!(scheme, B, Δx, y)
    end
    return x_next, ∇f_next
end

function nlsolve(f, ∇f, scheme, x0; inverse=true, c=0.001, g_tol=1e-8, maxiter=10^4, show_trace=false)
    n = length(x0)
    x_curr = copy(x0)
    ∇f_curr = ∇f(x_curr)

    B = initial_approximation(inverse, scheme, n)

    d = -solve(B, ∇f_curr)
    α, f_α, ls_success = backtracking(f, x_curr, d, ∇f_curr; α_0=1.0, c=c)
    Δx = α .* d
    x_next = x_curr + Δx
    ∇f_next = ∇f(x_next)
    y = ∇f_next - ∇f_curr
    update!(scheme, B, Δx, y)

    for i = 1:maxiter
        x_curr = copy(x_next)
        ∇f_curr = ∇f(x_curr)
        d = .- solve(B, ∇f_curr)
        α, f_α, ls_success = backtracking(f, x_curr, d, ∇f_curr; α_0=1.0, c=c)
        Δx = α .* d
        x_next = x_curr + Δx
        ∇f_next = ∇f(x_next)
        y = ∇f_next - ∇f_curr

        if norm(∇f_next) < g_tol
            return x_next, ∇f_next
        end
        if show_trace
            println("Objective value (curr): ", f(x_next))
            println("Objective value (prev): ", f(x_curr))
            println("Step size: ", α)
        end
        update!(scheme, B, Δx, y)
    end
    return x_next, ∇f_next
end


text:
"Autum Expressions"
module AExpressions

using MLStyle
export AExpr

export istypesymbol,
       istypevarsymbol,
       args,
       arg,
       wrap,
       showstring

const autumngrammar = """
x           := a | b | ... | aa ...
program     := statement*
statement   := externaldecl | assignexpr | typedecl | typedef

typedef     := type fields  #FIXME
typealias   := "type alias" type fields
fields      := field | fields field
field       := constructor | constructor typesymbol*
cosntructor := typesymbol

typedecl    := x : typeexpr
externaldecl:= external typedecl

assignexpr  := x = valueexpr

typeexpr    := typesymbol | paramtype | typevar | functiontype
funtype     := typeexpr -> typeexpr
producttype := typeexpr × typexexpr × ...
typesymbol  := primtype | customtype
primtype    := Int | Bool | Float
customtype  := A | B | ... | Aa | ...

valueexpr   := fappexpr | lambdaexpr | iteexpr | initnextexpr | letexpr |
               this | lambdaexpr
iteexpr     := if expr then expr else expr
intextexpr  := init expr next expr
fappexpr    := valueexpr valueexpr*
letexpr     := let x = valueexpr in valueexpr
lambdaexpr  := x --> expr
"""

"Autumn Expression"
struct AExpr
  head::Symbol
  args::Vector{Any}
  AExpr(head::Symbol, @nospecialize args...) = new(head, [args...])
end
"Arguements of expression"
function args end

args(aex::AExpr) = aex.args
head(aex::AExpr) = aex.head
args(ex::Expr) = ex.args

"Expr in ith location in arg"
arg(aex, i) = args(aex)[i]

Base.Expr(aex::AExpr) = Expr(aex.head, aex.args...)

# wrap(expr::Expr) = AExpr(expr)
# wrap(x) = x

# AExpr(xs...) = AExpr(Expr(xs...))

# function Base.getproperty(aexpr::AExpr, name::Symbol)
#   expr = getfield(aexpr, :expr)
#   if name == :expr
#     expr
#   elseif name == :head
#     expr.head
#   elseif name == :args
#     expr.args
#   else
#     error("no property $name of AExpr")
#   end
# end

# Expression types
"Is `sym` a type symbol"
istypesymbol(sym) = (q = string(sym); length(q) > 0 && isuppercase(q[1]))
istypevarsymbol(sym) = (q = string(sym); length(q) > 0 && islowercase(q[1]))

# ## Printing
isinfix(f::Symbol) = f ∈ [:+, :-, :/, :*, :&&, :||, :>=, :<=, :>, :<, :(==)]
isinfix(f) = false

"Pretty print"
function showstring(expr::Expr)
  @match expr begin
    Expr(:program, statements...) => join(map(showstring, expr.args), "\n")
    Expr(:producttype, ts) => join(map(showstring, ts), "×")
    Expr(:functiontype, int, outt) => "($(showstring(int)) -> $(showstring(outt)))"
    Expr(:functiontype, vars...) => reduce(((x, y) -> string(showstring(x), " -> ", showstring(y))), vars)
    Expr(:typedecl, x, val) => "$x : $(showstring(val))"
    Expr(:externaldecl, x, val) => "external $x : $(showstring(val))"
    Expr(:external, val) => "external $(showstring(val))"
    Expr(:assign, x, val) => "$x $(needequals(val)) $(showstring(val))"
    Expr(:if, i, t, e) => "if ($(showstring(i))) then ($(showstring(t))) else ($(showstring(e)))"
    Expr(:initnext, i, n) => "init $(showstring(i)) next $(showstring(n))"
    Expr(:args, args...) => join(map(showstring, args), " ")
    Expr(:call, f, arg1, arg2) && if isinfix(f) end => "$(showstring(arg1)) $f $(showstring(arg2))"
    Expr(:call, f, args...) => "($(join(map(showstring, [f ; args]), " ")))"
    Expr(:let, vars...) => "let \n\t$(join(map(showstring, vars[1:end-1]), "\n\t"))\nin\n\t$(showstring(vars[end]))"
    Expr(:paramtype, type, param) => string(type, " ", param)
    Expr(:paramtype, type) => string(type)
    Expr(:case, type, vars...) => string("\n\tcase $(showstring(type)) of \n\t\t", join(map(showstring, vars), "\n\t\t"))
    Expr(:casevalue, type, value) => "$(showstring(type)) => $value"
    Expr(:casetype, type, vars...) => "$type $(join(vars, " "))"
    Expr(:type, vars...) => "type $(vars[1]) $(vars[2]) = $(join(map(showstring, vars[3:length(vars)]), " | "))"
    Expr(:typealias, var, val) => "type alias $(showstring(var)) = $(showstring(val))"
    Expr(:fn, params, body) => "fn $(showstring(params)) ($(showstring(body)))"
    Expr(:list, vals...) => "($(join(vals, ", ")))"
    Expr(:field, var, field) => "$(showstring(var)).$(showstring(field))"
    Expr(:typealiasargs, vals...) => "$(string("{ ", join(map(showstring, vals), ", ")," }"))"
    Expr(:lambda, var, val) => "($(showstring(var)) -> ($(showstring(val))))"
    x                       => "Fail $x"

  end
end

showstring(lst::Array{}) = join(map(string, lst), " ")
showstring(str::String) = str

function needequals(val)
  if typeof(val) == Expr && val.head == :fn
    ""
   else
    "="
  end
end

showstring(aexpr::AExpr) = showstring(Expr(aexpr))

showstring(s::Union{Symbol, Integer}) = s
showstring(s::Type{T}) where {T <: Number} = s
Base.show(io::IO, aexpr::AExpr) = print(io, showstring(aexpr))

end


text:
## In this notebook, we'll look at stochastic gradient descent on a least squares
## system

##
using Plots
pyplot()
theme(:dark)
##
# https://vincentarelbundock.github.io/Rdatasets/csv/carData/Davis.csv
using CSV
df = CSV.read(download("https://vincentarelbundock.github.io/Rdatasets/csv/carData/Davis.csv"))
X = Float64.([df[:weight] df[:height]])
# Filter the data
goodpts = X[:,1] .<= 115
Xf = X[goodpts,:]
##
scatter(Xf[:,1],Xf[:,2],label="",xlabel="weight (kg)", ylabel="height (cm)")
##
x = Xf[:,1]
A = [x ones(length(x))]
b = Xf[:,2]
ab = A\b # solve the least squares problem

##
using Printf
using LinearAlgebra
""" Run stochastic gradient descent to solve a least squares problem. """
function stochastic_gradient_descent_leastsquares(A,b;
            step=0.0001, maxiter=10*size(A,1),tol=1e-6)
  x = zeros(size(A,2)) # start at 0
  for i=1:1000000
    si = rand(1:size(A,1))
    ai = A[si,:]
    bi = b[si]
    g = 2.0.*(ai'*x - bi)*ai
    x -= step*g
  end
  return x
end
stochastic_gradient_descent_leastsquares(A,b)
## Animate the algorithm
quadratic_leastsq(A,b) = x-> (0.5*norm(A*x - b)^2)/size(A,1)
ezsurf(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot(x, y, Z, st=:surface)
end
ezsurf(-3:0.1:3,-50:0.5:250,
  quadratic_leastsq(A,b))
plot!(colorbar=false,size=(450,450),dpi=300)
plot!(camera=(170,30))
x = zeros(size(A,2)) # initialize
p = plot3d!(1, label="", marker = 2)
step = 0.0001
anim = Animation()
for i=0:1000000
  global x
  g = (A'*A*x - A'*b)/size(A,1)
  x -= step*g
  if ((i <= 10) || ((i <= 30) & (rem(i,2) == 0)) || ((i <= 100) && rem(i,5) == 0) ||
      ((i <= 1000) && rem(i,10) == 0) || ((i <=10000) && rem(i,100)==0) ||
      ((i <= 100000) && rem(i,1000)==0) || rem(i,10000)==0)
    push!(p[1][2], x[1], x[2], quadratic_leastsq(A,b)(x))
    title!(@sprintf("Iteration %3i  Residual = %.3e", i, norm(g)))
    @printf("Iteration %3i  Residual = %.3e\n", i, norm(g))
    frame(anim)
  end
end
gif(anim, "stochatic-grad-descent-1.gif",fps=30)

##
## Animate the stochastic algorithm
using Random
Random.seed!(0)
quadratic_leastsq(A,b) = x-> (0.5*norm(A*x - b)^2)/size(A,1)
ezsurf(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot(x, y, Z, st=:surface)
end
ezsurf(-3:0.1:3,-50:0.5:250,
  quadratic_leastsq(A,b))
plot!(colorbar=false,size=(450,450),dpi=300)
plot!(camera=(170,30))
x = zeros(size(A,2)) # initialize
xgd = zeros(size(A,2)) # initialize
p = plot3d!(1, label="", color=3)
p2 = plot3d!(1, label="", color=2)
step = 0.0001
anim = Animation()
for i=0:1000000
  global x, xgd
  si = rand(1:size(A,1))
  ai = A[si,:]
  bi = b[si]
  g = (ai'*x - bi)*ai
  x -= step*g

  xgd -= step*(A'*A*xgd - A'*b)/size(A,1)
  if ((i <= 10) || ((i <= 30) & (rem(i,2) == 0)) || ((i <= 100) && rem(i,5) == 0) ||
      ((i <= 1000) && rem(i,10) == 0) || ((i <=10000) && rem(i,100)==0) ||
      ((i <= 100000) && rem(i,1000)==0) || rem(i,10000)==0)
    push!(p[1][2], x[1], x[2], quadratic_leastsq(A,b)(x))
    push!(p[1][3], xgd[1], xgd[2], quadratic_leastsq(A,b)(xgd))
    title!(@sprintf("Iteration %3i  Residual = %.3e", i, norm(g)))
    @printf("Iteration %3i  Residual = %.3e\n", i, norm(g))
    frame(anim)
  end
end
gif(anim, "stochatic-grad-descent-2.gif",fps=30)

## Show the quadratic surface
quadratic_leastsq(A,b) = x-> (0.5*norm(A*x - b)^2)
ezsurf(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot(x, y, Z, st=:surface)
end
ezsurf(-3:0.1:3,-50:0.5:250,
  quadratic_leastsq(A,b))
##
plot!(colorbar=false,size=(300,300),dpi=300)
anim = @animate for i=1:360
  plot!(camera=(i,30))
end
gif(anim, "leastsq-quadratic.gif")

##
## Show the quadratic surface with alpha-blending
quadratic_leastsq(A,b) = x-> (0.5*norm(A*x - b)^2)
ezsurf(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot()
    plot!(x, y, Z, st=:surface)
end
ezsurf(-3:0.1:3,-50:0.5:250,
  quadratic_leastsq(A,b))

##
## Show k random samples of the quadratic surface
plotlyjs()
quadratic_leastsq(A,b) = x-> (0.5*norm(A*x - b)^2)/size(A,1)
ezsurf!(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot!(x, y, Z, st=:surface, alpha=0.1)
end
orig = ezsurf(-3:0.25:3,-50:1:250,
  quadratic_leastsq(A,b))
nsamp = 25
k = 25
for i=1:k
  subset = randperm(size(A,1))[1:25] # choose 25 of 200 samples.
  ezsurf!(-3:0.1:3,-50:0.5:250,
    quadratic_leastsq(A[subset,:],b[subset,:]))
end
plot!() # get the plot to show
## Show k random samples of the quadratic surface
using Random
using LinearAlgebra
plotlyjs()
quadratic_leastsq(A,b) = x-> (0.5*norm(A*x - b)^2)/size(A,1)
ezsurf!(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot!(x, y, Z, st=:surface, alpha=0.4)
end
plot(framestyle=:none,colorbar=false,size=(300,300),dpi=300)
nsamp = 10
k = 25
for i=1:k
  subset = randperm(size(A,1))[1:nsamp] # choose 25 of 200 samples.
  ezsurf!(-3:0.25:3,-50:2:250,
    quadratic_leastsq(A[subset,:],b[subset,:]))
end
plot!() # get the plot to show
##
#plot!(colorbar=false,size=(300,300),dpi=300)
plot!(colorbar=false,size=(900,900))
anim = @animate for i=2:2:360
  plot!(camera=(i,30))
end
gif(anim, "leastsq-quadratic-samples.gif")

## Show the same thing with a log10
using Random
using LinearAlgebra
plotlyjs()
quadratic_leastsq(A,b) = x-> log10(0.5*norm(A*x - b)^2)/size(A,1)
plot(framestyle=:none,colorbar=false,size=(300,300),dpi=300)
nsamp = 10
k = 25
for i=1:k
  subset = randperm(size(A,1))[1:nsamp] # choose 25 of 200 samples.
  ezsurf!(-3:0.1:3,-50:0.5:250,
    quadratic_leastsq(A[subset,:],b[subset,:]))
end
plot!() # get the plot to show
##
#plot!(colorbar=false,size=(300,300),dpi=300)
plot!(colorbar=false,size=(900,900))
anim = @animate for i=2:2:360
  plot!(camera=(i,30))
end
gif(anim, "leastsq-quadratic-samples-log.gif")
##


text:
## source from: http://www.pjm.com/markets-and-operations/ops-analysis/historical-load-data.aspx
function load24_data(filename_Load24)
    temp_data = readcsv(filename_Load24,header=true)[1][:,3:end]
    raw_data = Array{Float64}(size(temp_data))
    normalized_data = Array{Float64}(size(temp_data))
    raw_data[:] = temp_data[:]
    N_days = size(raw_data,1)
    for k = 1:N_days
        normalized_data[k,:] = raw_data[k,:] / maximum(raw_data[k,:])
    end
    load24 = normalized_data
    return load24
end


text:
using Test
using KernelDensity

X = randn(100)
Y = randn(100)

k = kde(X)
@test pdf(k, k.x) ≈ k.density

k = kde((X,Y))
@test pdf(k, k.x, k.y) ≈ k.density

# Try to evaluate the KDE outside the interpolation domain
# The KDE is allowed to be zero, but it should not be greater than the exact solution
k = kde([0.0], bandwidth=1.0)
@test pdf(k, k.x) ≈ k.density
@test pdf(k, -10.0) ≤ pdf(Normal(), -10.0)
@test pdf(k, +10.0) ≤ pdf(Normal(), +10.0)

k = kde(([0.0],[0.0]), bandwidth=(1.0, 1.0))
@test pdf(k, k.x, k.y) ≈ k.density
@test pdf(k, -10.0, 0.0) ≤ pdf(MvNormal(2, 1.0), [-10.0, 0.0])
@test pdf(k, +10.0, 0.0) ≤ pdf(MvNormal(2, 1.0), [+10.0, 0.0])
@test pdf(k, 0.0, -10.0) ≤ pdf(MvNormal(2, 1.0), [0.0, -10.0])
@test pdf(k, 0.0, +10.0) ≤ pdf(MvNormal(2, 1.0), [0.0, +10.0])


text:
"""
    add_reactions!(model::CoreModel, rxns::Vector{Reaction})

Add `rxns` to `model` efficiently. The model must already contain the metabolites used by
`rxns` in the model.
"""
function add_reactions!(model::CoreModel, rxns::Vector{Reaction})
    I = Int64[] # rows
    J = Int64[] # cols
    V = Float64[] # values
    cs = zeros(length(rxns))
    lbs = zeros(length(rxns))
    ubs = zeros(length(rxns))
    for (j, rxn) in enumerate(rxns)
        req = rxn.metabolites
        for (i, v) in zip(indexin(keys(req), metabolites(model)), values(req))
            push!(J, j)
            push!(I, i)
            push!(V, v)
        end
        push!(model.rxns, rxn.id)
        lbs[j] = rxn.lb
        ubs[j] = rxn.ub
        cs[j] = rxn.objective_coefficient
    end
    Sadd = sparse(I, J, V, n_metabolites(model), length(rxns))
    model.S = [model.S Sadd]
    model.c = dropzeros([model.c; cs])
    model.xu = ubs
    model.xl = lbs
    return nothing
end

"""
    add_reaction!(model::CoreModel, rxn::Reaction)

Add `rxn` to `model`. The model must already contain the metabolites used by
`rxn` in the model.
"""
add_reaction!(model::CoreModel, rxn::Reaction) = add_reactions!(model, [rxn])

"""
    add_reactions(
        m::CoreModel,
        s::VecType,
        b::VecType,
        c::AbstractFloat,
        xl::AbstractFloat,
        xu::AbstractFloat;
        check_consistency = false,
    )

Add reaction(s) to a `CoreModel` model `m`.
"""
function add_reactions(
    m::CoreModel,
    s::VecType,
    b::VecType,
    c::AbstractFloat,
    xl::AbstractFloat,
    xu::AbstractFloat;
    check_consistency = false,
)
    return add_reactions(
        m,
        sparse(reshape(s, (length(s), 1))),
        sparse(b),
        sparse([c]),
        sparse([xl]),
        sparse([xu]),
        check_consistency = check_consistency,
    )
end

"""
    add_reactions(
        m::CoreModel,
        s::VecType,
        b::VecType,
        c::AbstractFloat,
        xl::AbstractFloat,
        xu::AbstractFloat,
        rxn::String,
        mets::K;
        check_consistency = false,
    )
"""
function add_reactions(
    m::CoreModel,
    s::VecType,
    b::VecType,
    c::AbstractFloat,
    xl::AbstractFloat,
    xu::AbstractFloat,
    rxn::String,
    mets::StringVecType;
    check_consistency = false,
)
    return add_reactions(
        m,
        sparse(reshape(s, (length(s), 1))),
        sparse(b),
        sparse([c]),
        sparse([xl]),
        sparse([xu]),
        [rxn],
        mets,
        check_consistency = check_consistency,
    )
end

"""
    add_reactions(
        m::CoreModel,
        Sp::MatType,
        b::VecType,
        c::VecType,
        xl::VecType,
        xu::VecType;
        check_consistency = false,
    )
"""
function add_reactions(
    m::CoreModel,
    Sp::MatType,
    b::VecType,
    c::VecType,
    xl::VecType,
    xu::VecType;
    check_consistency = false,
)
    rxns = ["r$x" for x = length(m.rxns)+1:length(m.rxns)+length(xu)]
    mets = ["m$x" for x = length(m.mets)+1:length(m.mets)+size(Sp)[1]]
    return add_reactions(
        m,
        Sp,
        b,
        c,
        xl,
        xu,
        rxns,
        mets,
        check_consistency = check_consistency,
    )
end

"""
    add_reactions(m1::CoreModel, m2::CoreModel; check_consistency = false)

Add all reactions from `m2` to `m1`.
"""
function add_reactions(m1::CoreModel, m2::CoreModel; check_consistency = false)
    return add_reactions(
        m1,
        m2.S,
        m2.b,
        m2.c,
        m2.xl,
        m2.xu,
        m2.rxns,
        m2.mets,
        check_consistency = check_consistency,
    )
end

"""
    add_reactions(
        m::CoreModel,
        Sp::MatType,
        b::VecType,
        c::VecType,
        xl::VecType,
        xu::VecType,
        rxns::StringVecType,
        mets::StringVecType;
        check_consistency = false,
    )
"""
function add_reactions(
    m::CoreModel,
    Sp::MatType,
    b::VecType,
    c::VecType,
    xl::VecType,
    xu::VecType,
    rxns::StringVecType,
    mets::StringVecType;
    check_consistency = false,
)
    Sp = sparse(Sp)
    b = sparse(b)
    c = sparse(c)
    xl = collect(xl)
    xu = collect(xu)

    all([length(b), length(mets)] .== size(Sp, 1)) ||
        throw(DimensionMismatch("inconsistent number of metabolites"))
    all(length.([c, xl, xu, rxns]) .== size(Sp, 2)) ||
        throw(DimensionMismatch("inconsistent number of reactions"))

    new_reactions = findall(Bool[!(rxn in m.rxns) for rxn in rxns])
    new_metabolites = findall(Bool[!(met in m.mets) for met in mets])

    if check_consistency
        (newReactions1, newMetabolites1) = verify_consistency(
            m,
            Sp,
            b,
            c,
            xl,
            xu,
            rxns,
            mets,
            new_reactions,
            new_metabolites,
        )
    end

    new_mets = vcat(m.mets, mets[new_metabolites])

    zero_block = spzeros(length(new_metabolites), n_reactions(m))
    ext_s = vcat(sparse(m.S), zero_block)

    mapping = [findfirst(isequal(met), new_mets) for met in mets]
    (I, J, elements) = findnz(sparse(Sp[:, new_reactions]))
    ext_sp = spzeros(length(new_mets), length(new_reactions))
    for (k, i) in enumerate(I)
        new_i = mapping[i]
        ext_sp[new_i, J[k]] = elements[k]
    end

    new_s = hcat(ext_s, ext_sp)
    newb = vcat(m.b, b[new_metabolites])
    newc = vcat(m.c, c[new_reactions])
    newxl = vcat(m.xl, xl[new_reactions])
    newxu = vcat(m.xu, xu[new_reactions])
    new_rxns = vcat(m.rxns, rxns[new_reactions])
    new_lp = CoreModel(new_s, newb, newc, newxl, newxu, new_rxns, new_mets)

    if check_consistency
        return (new_lp, new_reactions, new_metabolites)
    else
        return new_lp
    end
end

"""
    verify_consistency(
        m::CoreModel,
        Sp::M,
        b::V,
        c::V,
        xl::B,
        xu::B,
        names::K,
        mets::K,
        new_reactions,
        new_metabolites,
    ) where {M<:MatType,V<:VecType,B<:VecTypeK<:StringVecType}

Check the consistency of given reactions with existing reactions in `m`.

TODO: work in progress, doesn't return consistency status.
"""
function verify_consistency(
    m::CoreModel,
    Sp::M,
    b::V,
    c::V,
    xl::B,
    xu::B,
    names::K,
    mets::K,
    new_reactions,
    new_metabolites,
) where {M<:MatType,V<:VecType,B<:VecType,K<:StringVecType}

    if !isempty(new_reactions)
        statuses = Vector{ReactionStatus}(undef, length(names))
        for (i, name) in enumerate(names)
            rxn_index = findfirst(isequal(name), m.rxns)
            reaction = Sp[:, i]
            stoich_index = findfirst(Bool[reaction == m.S[:, j] for j = 1:size(m.S, 2)])
            if isnothing(rxn_index) & isnothing(stoich_index)
                statuses[i] = ReactionStatus(false, 0, "new")
            end

            if !isnothing(rxn_index) & isnothing(stoich_index)
                statuses[i] = ReactionStatus(true, 0, "same name")
            end

            if isnothing(rxn_index) & !isnothing(stoich_index)
                statuses[i] = ReactionStatus(true, 0, "same stoichiometry")
            end

            if !isnothing(rxn_index) & !isnothing(stoich_index)
                statuses[i] = ReactionStatus(true, 0, "same name, same stoichiometry")
            end
        end
    end

    return (new_reactions, new_metabolites)
end

@_change_bounds_fn CoreModel Int inplace begin
    isnothing(lower) || (model.xl[rxn_idx] = lower)
    isnothing(upper) || (model.xu[rxn_idx] = upper)
    nothing
end

@_change_bounds_fn CoreModel Int inplace plural begin
    for (i, l, u) in zip(rxn_idxs, lower, upper)
        change_bound!(model, i, lower = l, upper = u)
    end
end

@_change_bounds_fn CoreModel Int begin
    change_bounds(model, [rxn_idx], lower = [lower], upper = [upper])
end

@_change_bounds_fn CoreModel Int plural begin
    n = copy(model)
    n.xl = copy(n.xl)
    n.xu = copy(n.xu)
    change_bounds!(n, rxn_idxs, lower = lower, upper = upper)
    n
end

@_change_bounds_fn CoreModel String inplace begin
    change_bounds!(model, [rxn_id], lower = [lower], upper = [upper])
end

@_change_bounds_fn CoreModel String inplace plural begin
    change_bounds!(
        model,
        Vector{Int}(indexin(rxn_ids, reactions(model))),
        lower = lower,
        upper = upper,
    )
end

@_change_bounds_fn CoreModel String begin
    change_bounds(model, [rxn_id], lower = [lower], upper = [upper])
end

@_change_bounds_fn CoreModel String plural begin
    change_bounds(
        model,
        Int.(indexin(rxn_ids, reactions(model))),
        lower = lower,
        upper = upper,
    )
end

@_remove_fn reaction CoreModel Int inplace begin
    remove_reactions!(model, [reaction_idx])
end

@_remove_fn reaction CoreModel Int inplace plural begin
    mask = .!in.(1:n_reactions(model), Ref(reaction_idxs))
    model.S = model.S[:, mask]
    model.c = model.c[mask]
    model.xl = model.xl[mask]
    model.xu = model.xu[mask]
    model.rxns = model.rxns[mask]
    nothing
end

@_remove_fn reaction CoreModel Int begin
    remove_reactions(model, [reaction_idx])
end

@_remove_fn reaction CoreModel Int plural begin
    n = copy(model)
    remove_reactions!(n, reaction_idxs)
    return n
end

@_remove_fn reaction CoreModel String inplace begin
    remove_reactions!(model, [reaction_id])
end

@_remove_fn reaction CoreModel String inplace plural begin
    remove_reactions!(model, Int.(indexin(reaction_ids, reactions(model))))
end

@_remove_fn reaction CoreModel String begin
    remove_reactions(model, [reaction_id])
end

@_remove_fn reaction CoreModel String plural begin
    remove_reactions(model, Int.(indexin(reaction_ids, reactions(model))))
end

@_remove_fn metabolite CoreModel Int inplace begin
    remove_metabolites!(model, [metabolite_idx])
end

@_remove_fn metabolite CoreModel Int plural inplace begin
    remove_reactions!(
        model,
        [
            ridx for ridx in 1:n_reactions(model) if
            any(in.(findnz(model.S[:, ridx])[1], Ref(metabolite_idxs)))
        ],
    )
    mask = .!in.(1:n_metabolites(model), Ref(metabolite_idxs))
    model.S = model.S[mask, :]
    model.b = model.b[mask]
    model.mets = model.mets[mask]
    nothing
end

@_remove_fn metabolite CoreModel Int begin
    remove_metabolites(model, [metabolite_idx])
end

@_remove_fn metabolite CoreModel Int plural begin
    n = deepcopy(model) #everything gets changed anyway
    remove_metabolites!(n, metabolite_idxs)
    return n
end

@_remove_fn metabolite CoreModel String inplace begin
    remove_metabolites!(model, [metabolite_id])
end

@_remove_fn metabolite CoreModel String inplace plural begin
    remove_metabolites!(model, Int.(indexin(metabolite_ids, metabolites(model))))
end

@_remove_fn metabolite CoreModel String begin
    remove_metabolites(model, [metabolite_id])
end

@_remove_fn metabolite CoreModel String plural begin
    remove_metabolites(model, Int.(indexin(metabolite_ids, metabolites(model))))
end

"""
    change_objective!(
        model::CoreModel,
        rxn_idxs::Vector{Int};
        weights = ones(length(rxn_idxs)),
    )

Change the objective to reactions at given indexes, optionally specifying their
`weights` in the same order. By default, all set weights are 1.
"""
function change_objective!(
    model::CoreModel,
    rxn_idxs::Vector{Int};
    weights = ones(length(rxn_idxs)),
)
    model.c = spzeros(length(model.c))
    model.c[rxn_idxs] .= weights
    nothing
end

"""
    change_objective!(model::CoreModel, rxn_idx::Int)

Change objective function of a CoreModel to a single `1` at reaction index
`rxn_idx`.
"""
change_objective!(model::CoreModel, rxn_idx::Int) = change_objective!(model, [rxn_idx])

"""
    change_objective!(
        model::CoreModel,
        rxn_ids::Vector{String};
        weights = ones(length(rxn_ids)),
    )

Change objective of given reaction IDs, optionally specifying objective
`weights` in the same order as `rxn_ids`. By default, all set weights are 1.
"""
function change_objective!(
    model::CoreModel,
    rxn_ids::Vector{String};
    weights = ones(length(rxn_ids)),
)
    idxs = indexin(rxn_ids, reactions(model))
    any(isnothing(idx) for idx in idxs) &&
        throw(DomainError(rxn_ids, "Some reaction ids not found in the model"))
    change_objective!(model, Int.(idxs); weights)
end

"""
    change_objective!(model::CoreModel, rxn_id::String)

Change objective function of a CoreModel to a single `1` at the given reaction
ID.
"""
change_objective!(model::CoreModel, rxn_id::String) = change_objective!(model, [rxn_id])


text:
using Jecco
using LinearAlgebra
using SparseArrays
using Plots

#= solving equation

 f_xx + f_yy + x y f_xy + x f_x + y f_y + (x^2 + y^2) f
  + (f_x)^2 + (f_y)^2 + x y f_x f_y = Source(x,y)

for the source below, the solution should be given by

 f = exp(-x^2 - y^2)
=#

source(x,y) = exp(-x^2 - y^2) * (-4 + 3 * (x^2 + y^2) + 4 * x^2 * y^2) +
    4 * exp(-2*x^2 - 2*y^2) * (x^2 + y^2 + x^2 * y^2)

function compute_residual!(res::Array, f::Array, x::Vector, y::Vector,
                           Dxx_op::Jecco.FiniteDiffDeriv, Dyy_op::Jecco.FiniteDiffDeriv,
                           Dx_op::Jecco.FiniteDiffDeriv, Dy_op::Jecco.FiniteDiffDeriv)
    Nx = length(x)
    Ny = length(y)

    f_xx = Dxx_op * f
    f_yy = Dyy_op * f
    f_x  = Dx_op  * f
    f_y  = Dy_op  * f
    f_xy = Dx_op  * f_y

    for j in 1:Ny, i in 1:Nx
        x1 = x[i]
        y1 = y[j]

        res[i,j] = f_xx[i,j] + f_yy[i,j] + x1 * y1 * f_xy[i,j] +
            x1 * f_x[i,j] + y1 * f_y[i,j] +
            (x1*x1 + y1*y1) * f[i,j] + f_x[i,j]^2 + f_y[i,j]^2 +
            x1 * y1 * f_x[i,j] * f_y[i,j] - source(x1,y1)
    end

    res
end

function deriv_operators(hx, hy, Nx::Int, Ny::Int, ord::Int)
    Dx_op  = CenteredDiff{1}(1, ord, hx, Nx)
    Dxx_op = CenteredDiff{1}(2, ord, hx, Nx)

    Dy_op  = CenteredDiff{2}(1, ord, hy, Ny)
    Dyy_op = CenteredDiff{2}(2, ord, hy, Ny)

    Dxx_op, Dyy_op, Dx_op, Dy_op
end

#=
use the Kronecker product (kron) to build the 2-dimensional derivation matrices
from the 1-dimensional ones. see for instance:

  https://en.wikipedia.org/wiki/Kronecker_product

  https://arxiv.org/pdf/1801.01483.pdf (section 5)
=#
function deriv_matrices(Dxx_op::Jecco.FiniteDiffDeriv, Dyy_op::Jecco.FiniteDiffDeriv,
                        Dx_op::Jecco.FiniteDiffDeriv, Dy_op::Jecco.FiniteDiffDeriv)
    Dx  = kron(I(Ny), SparseMatrixCSC(Dx_op))
    Dxx = kron(I(Ny), SparseMatrixCSC(Dxx_op))
    Dy  = kron(SparseMatrixCSC(Dy_op), I(Nx))
    Dyy = kron(SparseMatrixCSC(Dyy_op), I(Nx))

    Dxx, Dyy, Dx * Dy, Dx, Dy
end

# returns axx Dxx + ayy Dyy + axy Dxy + bx Dx + by Dy + cc. note that this
# function overwrites the input matrices to save memory
function build_jacobian(Dxx::SparseMatrixCSC, Dyy::SparseMatrixCSC, Dxy::SparseMatrixCSC,
                        Dx::SparseMatrixCSC, Dy::SparseMatrixCSC,
                        axx::Diagonal, ayy::Diagonal, axy::Diagonal,
                        bx::Diagonal, by::Diagonal, cc::Diagonal)
    Dxx_ = similar(Dxx)
    Dyy_ = similar(Dyy)
    Dxy_ = similar(Dxy)
    Dx_  = similar(Dx)
    Dy_  = similar(Dy)

    mul!(Dxx_, axx, Dxx)
    mul!(Dyy_, ayy, Dyy)
    mul!(Dxy_, axy, Dxy)
    mul!(Dx_, bx, Dx)
    mul!(Dy_, by, Dy)

    Dxx_ + Dyy_ + Dxy_ + Dx_ + Dy_ + cc
end

x_min    = -5.0
x_max    =  5.0
x_nodes  =  128
# x_nodes  =  256
y_min    = -5.0
y_max    =  5.0
y_nodes  =  64
# y_nodes  =  256

ord = 4
# ord = 2

xcoord  = Cartesian{1}("x", x_min, x_max, x_nodes, endpoint=false)
ycoord  = Cartesian{2}("y", y_min, y_max, y_nodes, endpoint=false)

hx = Jecco.delta(xcoord)
hy = Jecco.delta(ycoord)

Nx = xcoord.nodes
Ny = ycoord.nodes

f_exact = [exp(-xcoord[i]^2 - ycoord[j]^2) for i in 1:Nx, j in 1:Ny]

Dxx_op, Dyy_op, Dx_op, Dy_op = deriv_operators(hx, hy, Nx, Ny, ord)
Dxx, Dyy, Dxy, Dx, Dy = deriv_matrices(Dxx_op, Dyy_op, Dx_op, Dy_op)

# initial guess
# fsol   = zeros(Nx,Ny)
fsol   = copy(f_exact)

res    = zeros(Nx,Ny)
ind2D  = LinearIndices(fsol)

M = Nx * Ny

b_vec   = zeros(M)
axx     = ones(M)
ayy     = ones(M)
axy     = zeros(M)
bx      = zeros(M)
by      = zeros(M)
cc      = zeros(M)
f0      = zeros(M)

itmax   = 8
epsilon = 1e-12

for it in 1:itmax

    compute_residual!(res, fsol, xcoord[:], ycoord[:], Dxx_op, Dyy_op, Dx_op, Dy_op)
    max_res = maximum(abs.(res))

    @show it, max_res
    if max_res < epsilon
        break
    end

    # build (linearized) operator
    # A = Dxx + Dyy + x y Dxy + (x + 2 f0_x + x y f0_y) Dx + (y + 2 f0_y + x y f0_x) Dy + (x^2 + y^2)

    @inbounds for idx in eachindex(f0)
        f0[idx] = fsol[idx]
    end
    f0_x = Dx * f0
    f0_y = Dy * f0

    for j in 1:Ny, i in 1:Nx
        idx = ind2D[i,j]

        x1  = xcoord[i]
        y1  = ycoord[j]

        axy[idx] = x1 * y1
        bx[idx]  = x1 + 2 * f0_x[idx] + x1 * y1 * f0_y[idx]
        by[idx]  = y1 + 2 * f0_y[idx] + x1 * y1 * f0_x[idx]
        cc[idx]  = x1^2 + y1^2

        b_vec[idx] = -res[idx]
    end

    A_mat = build_jacobian(Dxx, Dyy, Dxy, Dx, Dy,
                           Diagonal(axx), Diagonal(ayy), Diagonal(axy),
                           Diagonal(bx), Diagonal(by), Diagonal(cc))

    A_fact = factorize(A_mat)
    ldiv!(f0, A_fact, b_vec)

    # update solution
    @inbounds for idx in eachindex(fsol)
        fsol[idx] += f0[idx]
    end

end

j_slice = div(Ny,2) + 1
x = xcoord[:]

plot(x, f_exact[:,j_slice])
scatter!(x, fsol[:,j_slice])


text:
@testset "Interval matrix construction" begin
    m1 = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.8..4.2 0.0..0.9])
    m2 = IntervalMatrix{Float64}(undef, 2, 2)
    @test m2 isa IntervalMatrix{Float64} && size(m2) == (2, 2)
    m3 = similar(m1)
    @test m3 isa IntervalMatrix{Float64} && size(m3) == size(m1)
    m = [1.0 2.0; 3.0 4.0]
    mint = IntervalMatrix([Interval(1) Interval(2); Interval(3) Interval(4)])
    @test IntervalMatrix(m) == mint

    A = [1 2; 3 4]
    B = [1 2; 4 5]

    @test IntervalMatrix(A, B) == IntervalMatrix([1..1 2..2; 3..4 4..5])

    @test A ± B == IntervalMatrix([0..2 0..4;-1..7 -1..9])
end

@testset "Interval matrix midpoint_radius" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    c, s = midpoint_radius(m)
    @test c ≈ [-0.1 -4.; 4. -0.1]
    @test s ≈ [1. 0.1; 0.1 1.]
end

@testset "Complex interval matrices" begin
    m1 = IntervalMatrix([(1..2)+im*(3..4) 1;2 3])
    @test m1 isa IntervalMatrix{Float64}
    @test eltype(m1) == Complex{Interval{Float64}}

    rp = real(m1)
    ip = imag(m1)
    @test rp isa IntervalMatrix{Float64}
    @test eltype(rp) == Interval{Float64}

    @test ip isa IntervalMatrix{Float64}
    @test eltype(ip) == Interval{Float64}

    m2 = IntervalMatrix([1+im 2+im;3+im 4+im])
    @test m2 isa IntervalMatrix{Float64}
    @test eltype(m2) == Complex{Interval{Float64}}
end


text:
mutable struct Model{T, TV<:AbstractVector{T}, TC<:AbstractVector{<:Function}}
    dim::Int
    objective::Function
    ineq_constraints::TC
    box_max::TV
    box_min::TV
end
GPUUtils.whichdevice(m::Model) = whichdevice(m.box_max)

dim(m::Model) = m.dim
min(m::Model, i::Integer) = m.box_min[i]
max(m::Model, i::Integer) = m.box_max[i]
min(m::Model)= m.box_max
max(m::Model) = m.box_min
objective(m::Model) = m.objective
constraints(m::Model) = m.ineq_constraints
constraint(m::Model, i::Integer) = m.ineq_constraints[i]

eval_objective(m, x::AbstractVector{T}) where {T} = eval_objective(m, x, T[])
eval_objective(m, x, ∇g) = eval_objective(whichdevice(objective(m)), m, x, ∇g)
function eval_objective(::CPU, m, x::AbstractVector{T}, ∇g) where {T}
    return T(m.objective(x, ∇g))
end
eval_objective(::GPU, m, x::GPUVector{T}, ∇g) where {T} = T(m.objective(x, ∇g))
function eval_objective(::GPU, m, x::AbstractVector{T}, ∇g) where {T}
    x_gpu = CuArray(x)
    ∇g_gpu = CuArray(∇g)
    obj = T(m.objective(x_gpu, ∇g_gpu))
    copyto!(∇g, ∇g_gpu)
    return obj
end
function eval_objective(::CPU, m, x::CuVector{T}, ∇g) where {T}
    error("Optimization on the GPU with the objective evaluation on the CPU is weird!")
end

eval_constraint(m, i, x::AbstractVector{T}) where {T} = eval_constraint(m, i, x, T[])
eval_constraint(m, i, x, ∇g) = eval_constraint(whichdevice(constraint(m, i)), m, i, x, ∇g)
eval_constraint(::CPU, m, i, x::AbstractVector{T}, ∇g) where T = T(constraint(m, i)(x, ∇g))
eval_constraint(::GPU, m, i, x::GPUVector{T}, ∇g) where T = T(constraint(m, i)(x, ∇g))
function eval_constraint(::GPU, m, i, x::AbstractVector{T}, ∇g) where T
    x_gpu = CuArray(x)
    ∇g_gpu = CuArray(∇g)
    constr = T(constraint(m, i)(x_gpu, ∇g_gpu))
    copyto!(∇g, ∇g_gpu)
    return constr
end
function eval_constraint(::CPU, m, i, x::GPUVector, ∇g)
    error("Optimization on the GPU with the constraint evaluation on the CPU is weird!")
end

ftol(m) = m.ftol
xtol(m) = m.xtol
grtol(m) = m.grtol
ftol!(m, v) = m.ftol = v
xtol!(m, v) = m.xtol = v
grtol!(m, v) = m.grtol = v

Model(args...; kwargs...) = Model{CPU}(args...; kwargs...)
Model{T}(args...; kwargs...) where T = Model(T(), args...; kwargs...)
Model(::CPU, args...; kwargs...) = Model{Float64, Vector{Float64}, Vector{Function}}(args...; kwargs...)
Model(::GPU, args...; kwargs...) = Model{Float64, CuVector{Float64}, Vector{Function}}(args...; kwargs...)

function Model{T, TV, TC}(dim, objective::Function) where {T, TV, TC}
    mins = ninfsof(TV, dim)
    maxs = infsof(TV, dim)
    Model{T, TV, TC}(dim, objective, Function[],
             mins, maxs)
end

# Box constraints
function box!(m::Model, i::Integer, minb::T, maxb::T) where {T}
    if !(1 <= i <= dim(m))
        throw(ArgumentError("box constraint need to applied to an existing variable"))
    end
    m.box_min[i] = minb
    m.box_max[i] = maxb
end

function box!(m::Model, minb::T, maxb::T) where {T}
    nv = dim(m)
    m.box_min[1:nv] .= minb
    m.box_max[1:nv] .= maxb
end

function box!(m::Model, minbs::AbstractVector{T}, maxbs::AbstractVector{T}) where {T}
    if (length(minbs) != dim(m)) || (length(minbs) != dim(m))
        throw(ArgumentError("box constraint vector must have same size as problem dimension"))
    end
    nv = dim(m)
    map!(identity, m.box_min, minbs)
    map!(identity, m.box_max, maxbs)
end

function ineq_constraint!(m::Model, f::Function)
    push!(m.ineq_constraints, f)
end

function ineq_constraint!(m::Model, fs::Vector{Function})
    for f in fs
        push!(m.ineq_constraints, f)
    end
end


text:
export iauRy
"""
Rotate an r-matrix about the y-axis.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   theta  double          angle (radians)

Given and returned:
   r      double[3][3]    r-matrix, rotated

Notes:

   1. Calling this function with positive theta incorporates in the
   supplied r-matrix r an additional rotation, about the y-axis,
   anticlockwise as seen looking towards the origin from positive y.

   2. The additional rotation can be represented by this matrix:

      (  + cos(theta)     0      - sin(theta)  )
      (                                        )
      (       0           1           0        )
      (                                        )
      (  + sin(theta)     0      + cos(theta)  )

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauRy(double theta, double r[3][3])
function iauRy(theta::Real, r::AbstractMatrix{<:Real})

   # Transpose matrix upfront
   r = Array{Float64, 2}(r') # Transpose input up front

   ccall((:iauRy, libsofa_c), Cvoid,
           (Cdouble, Ptr{Cdouble}),
           theta, pointer(r))

   return convert(Array{Float64, 2}, r')
end


text:
##################################################################
"""
$(TYPEDEF)

Structure holding data for finite volume system solution.
Information on species distribution is kept in dense
matrices, and the solution array is of type Array{2}.

Unlike in the SparseSystem, the system matrix handles exactly those
degrees of freedom which correspond to unknowns, and dummy
degrees of freedom where unknowns are not defined. Handling
of the sparse matrix structures for the bookeeping of the unknowns
has less overhead, but additional dummy equations are added
to the system matrix.

$(TYPEDFIELDS)
"""
mutable struct DenseSystem{Tv,Ti, Tm} <: AbstractSystem{Tv,Ti, Tm}
    """
    Grid
    """
    grid

    """
    Physics data
    """
    physics::Physics

    """
    Array of boundary values
    """
    boundary_values::Array{Tv,2}

    """
    Array of boundary factors
    """
    boundary_factors::Array{Tv,2}

    """
    Full matrix containing species numbers for inner regions
    """
    region_species::Array{Int8,2}

    """
    Full matrix containing species numbers for boundary regions
    """
    bregion_species::Array{Int8,2}

    """
    Full matrix containing degree of freedom numbers for each node
    """
    node_dof::Array{Int8,2}

    """
    Jacobi matrix for nonlinear problem
    """
    matrix::ExtendableSparseMatrix{Tv,Tm}

    """
    Matrix factorization
    """
    factorization::Union{Nothing,ExtendableSparse.AbstractFactorization{Tv,Tm}}

    """
    Flag which says if the number of unknowns per node is constant
    """
    species_homogeneous::Bool

    """
    Solution vector holding Newton update
    """
    update::DenseSolutionArray{Tv}

    """
    Solution vector holding Newton residual
    """
    residual::DenseSolutionArray{Tv}

    """
    Precomputed geometry factors for cell nodes
    """
    cellnodefactors::Array{Tv,2}

    """
    Precomputed geometry factors for cell edges
    """
    celledgefactors::Array{Tv,2}

    """
    Precomputed geometry factors for boundary nodes
    """
    bfacenodefactors::Array{Tv,2}

    """
    Precomputed geometry factors for boundary edges
    """
    bfaceedgefactors::Array{Tv,2}

    """
    Sparse matrix for generic operator handling
    """
    generic_matrix::SparseMatrixCSC

    """
    Sparse matrix colors for generic operator handling
    """
    generic_matrix_colors::Vector

    """
    Hash value of latest unknowns vector the assembly was called with
    """
    uhash::UInt64

    """
    Data for allocation check
    """
    allocs::Int

    DenseSystem{Tv,Ti, Tm}() where {Tv,Ti, Tm} = new()
end
##################################################################
"""
$(SIGNATURES)

Constructor for DenseSystem.
"""
function  DenseSystem(grid,physics::Physics;matrixindextype=Int64)
    Tv      = coord_type(grid)
    Ti      = index_type(grid)
    Tm      = matrixindextype
    system  = DenseSystem{Tv,Ti,Tm}()
    maxspec = physics.num_species

    system.grid                = grid
    system.physics             = physics
    system.region_species      = spzeros(Int8,Int16,maxspec,num_cellregions(grid))
    system.bregion_species     = spzeros(Int8,Int16,maxspec,num_bfaceregions(grid))
    system.node_dof            = spzeros(Int8,Int32,maxspec,num_nodes(grid))
    system.boundary_values     = zeros(Tv,maxspec,num_bfaceregions(grid))
    system.boundary_factors    = zeros(Tv,maxspec,num_bfaceregions(grid))
    system.species_homogeneous = false
    system.uhash               = 0x0
    system.allocs              = -1000
    system.factorization       = nothing
    return system
end

##################################################################
"""
$(SIGNATURES)

Number of degrees of freedom for system.
"""
num_dof(system::DenseSystem)= length(system.node_dof)

##################################################################
"""
$(SIGNATURES)

Create a solution vector for dense system.
If inival is not specified, the entries of the returned vector are undefined.
"""
unknowns(system::DenseSystem{Tv};inival=undef) where Tv = unknowns(Tv,system,inival=inival)

##################################################################
"""
$(SIGNATURES)

Create a solution vector for dense system with elements of type `Tu`.
If inival is not specified, the entries of the returned vector are undefined.
"""
function unknowns(Tu::Type, sys::DenseSystem{Tv}; inival=undef) where Tv
    a=Array{Tu}(undef,size(sys.node_dof)...)
    if inival!=undef
        fill!(a,inival)
    end
    return a
end

Base.reshape(v::DenseSolutionArray{Tu},system::DenseSystem{Tv}) where {Tu,Tv}=v

##################################################################
"""
$(SIGNATURES)

Reshape vector to fit as solution to system.
"""
function Base.reshape(v::AbstractVector{Tu}, sys::DenseSystem{Tv}) where {Tu,Tv}
    @assert  length(v)==num_dof(sys)
    nspec=num_species(sys)
    reshape(v,Int64(nspec),Int64(length(v)/nspec))
end

#
# Assemble dummy equations for inactive species
#
function _eval_and_assemble_inactive_species(system::DenseSystem,U,Uold,F)
    if system.species_homogeneous
        return
    end
    for inode=1:size(system.node_dof,2)
        for ispec=1:size(system.node_dof,1)
            if !isdof(system,ispec,inode)
                F[ispec,inode]+= U[ispec,inode]-Uold[ispec,inode];
                idof=dof(F,ispec,inode)
                system.matrix[idof,idof]+=1.0
            end
        end
    end
end

#
# Initialize values in inactive dof for dense system
#
function _initialize_inactive_dof!(U::DenseSolutionArray{Tv},system::DenseSystem{Tv}) where {Tv}
    if system.species_homogeneous
        return
    end
    for inode=1:size(system.node_dof,2)
        for ispec=1:size(system.node_dof,1)
            if !isdof(system,ispec,inode)
                U[ispec,inode]=0
            end
        end
    end
end

"""
$(SIGNATURES)

Calculate norm, paying attention to species distribution over regions
"""
function LinearAlgebra.norm(system::DenseSystem,u,p)
    _initialize_inactive_dof!(u,system)
    norm(u,p)
end


text:
"""
module FymEnvs

`FymEnvs` is developed by research project `NRF`
of `FDCL` (fligth dynamics and control laboratory)
in `SNU` (seoul national university).

The origin version of this package, `Fym`, has been developed in `Python`.
This is a Julia version of the `Fym`.

In this module,
some parts follows the convention of `Fym`,
e.g., `BaseEnv` and `BaseSystem`.
"""
module FymEnvs

using Reexport

# to avoid conflict; same name functions in different modules
function close! end
function record end

# Note: .FymLogging ->(reexport) .FymCore ->(reexport) .FymModels
# include("FymCore.jl")
# include("FymLogging.jl")
include("FymModels.jl")

# @reexport using .FymCore
# @reexport using .FymLogging
@reexport using .FymModels

end


text:
# BED Writer
# ==========

"""
    BEDWriter(output::IO)

Create a data writer of the BED file format.

# Arguments
* `output`: data sink
"""
type BEDWriter{T<:IO} <: Bio.IO.AbstractWriter
    output::T

    # number of written additional fields (-1: inferred from the first interval)
    n_fields::Int
end

function Bio.IO.stream(writer::BEDWriter)
    return writer.output
end

function BEDWriter(output::IO, n_fields::Integer=-1)
    return BEDWriter{typeof(output)}(output, n_fields)
end

function Base.write(writer::BEDWriter, interval::BEDInterval)
    if writer.n_fields == -1
        writer.n_fields = interval.metadata.used_fields
    end

    if interval.metadata.used_fields < writer.n_fields
        error("interval doesn't have enough additional fields")
    end

    output = writer.output
    n_fields = writer.n_fields
    n = 0
    n += write(
        output,
        interval.seqname, '\t',
        dec(interval.first - 1), '\t',
        dec(interval.last))

    if n_fields ≥ 1
        n += write(output, '\t', interval.metadata.name)
    end

    if n_fields ≥ 2
        n += write(output, '\t', dec(interval.metadata.score))
    end

    if n_fields ≥ 3
        n += write(output, '\t', Char(interval.strand))
    end

    if n_fields ≥ 4
        n += write(output, '\t', dec(interval.metadata.thick_first - 1))
    end

    if n_fields ≥ 5
        n += write(output, '\t', dec(interval.metadata.thick_last))
    end

    if n_fields ≥ 6
        rgb = interval.metadata.item_rgb
        n += write(
            output, '\t',
            dec(round(Int, 255 * rgb.r)), ',',
            dec(round(Int, 255 * rgb.g)), ',',
            dec(round(Int, 255 * rgb.b)))
    end

    if n_fields ≥ 7
        n += write(output, '\t', dec(interval.metadata.block_count))
    end

    if n_fields ≥ 8
        block_sizes = interval.metadata.block_sizes
        if !isempty(block_sizes)
            n += write(output, '\t', dec(block_sizes[1]))
            for i in 2:length(block_sizes)
                n += write(output, ',', dec(block_sizes[i]))
            end
        end
    end

    if n_fields ≥ 9
        block_firsts = interval.metadata.block_firsts
        if !isempty(block_firsts)
            n += write(output, '\t', dec(block_firsts[1] - 1))
            for i in 2:length(block_firsts)
                n += write(output, ',', dec(block_firsts[i] - 1))
            end
        end
    end

    n += write(output, '\n')
    return n
end

function Base.write(out::IO, interval::BEDInterval)
    print(out, interval.seqname, '\t', interval.first - 1, '\t', interval.last)
    write_optional_fields(out, interval)
    write(out, '\n')
end

function write_optional_fields(out::IO, interval::BEDInterval, leadingtab::Bool=true)
    if interval.metadata.used_fields >= 1
        if leadingtab
            write(out, '\t')
        end
        print(out, interval.metadata.name)
    else return end

    if interval.metadata.used_fields >= 2
        print(out, '\t', interval.metadata.score)
    else return end

    if interval.metadata.used_fields >= 3
        print(out, '\t', interval.strand)
    else return end

    if interval.metadata.used_fields >= 4
        print(out, '\t', interval.metadata.thick_first - 1)
    else return end

    if interval.metadata.used_fields >= 5
        print(out, '\t', interval.metadata.thick_last)
    else return end

    if interval.metadata.used_fields >= 6
        item_rgb = interval.metadata.item_rgb
        print(out, '\t',
              round(Int, 255 * item_rgb.r), ',',
              round(Int, 255 * item_rgb.g), ',',
              round(Int, 255 * item_rgb.b))
    else return end

    if interval.metadata.used_fields >= 7
        print(out, '\t', interval.metadata.block_count)
    else return end

    if interval.metadata.used_fields >= 8
        block_sizes = interval.metadata.block_sizes
        if !isempty(block_sizes)
            print(out, '\t', block_sizes[1])
            for i in 2:length(block_sizes)
                print(out, ',', block_sizes[i])
            end
        end
    else return end

    if interval.metadata.used_fields >= 9
        block_firsts = interval.metadata.block_firsts
        if !isempty(block_firsts)
            print(out, '\t', block_firsts[1] - 1)
            for i in 2:length(block_firsts)
                print(out, ',', block_firsts[i] - 1)
            end
        end
    end
end


text:
using Statistics

(X, Y) = bloobs(randomize = false)

μ0 = vec(mean(X, dims = 1))
κ0 = 5.0
ν0 = 9.0
Σ0 = cov(X)

H = WishartGaussian(μ0, κ0, ν0, Σ0)

model = DPM(H)

initialisation = RandomInitialisation(k = 10)
modelBuffer = init(X, model, initialisation)

model0 = BayesianNonparametrics.extractpointestimate(modelBuffer)
model1 = train(modelBuffer, DPMHyperparam(), Gibbs(maxiter = 1))[end]

initialisation = KMeansInitialisation(k = 10)
modelBuffer = init(X, model, initialisation)

model0 = BayesianNonparametrics.extractpointestimate(modelBuffer)
model1 = train(modelBuffer, DPMHyperparam(), Gibbs(maxiter = 1))[end]


text:
# Creates the expression that should replace the original call
function jitexpr(funexpr, callname, argname, fixtypes)
    retval = Expr(:block)
    push!(retval.args, profileexpr(callname, argname))
    next_target = retval.args
    for fixtype in fixtypes
        @assert fixtype isa Type
        ex = Expr(:if, :($argname isa $fixtype), :($funexpr))
        push!(next_target, ex)
        next_target = ex.args
    end
    push!(next_target, funexpr)
    return retval
end

function rreplace(expr::Expr, query, transform; symbols=false)
    if query(expr)
        return transform(expr)
    end
    args = [rreplace(subexpr, query, transform; symbols=symbols) for subexpr in expr.args]
    return Expr(expr.head, args...)
end
rreplace(expr::QuoteNode, query, transform; symbols=false) = QuoteNode(rreplace(expr.value, query, transform; symbols=symbols))
rreplace(expr::Symbol, query, transform; symbols=false) = symbols && query(expr) ? transform(expr) : expr
rreplace(expr, query, transform; symbols=false) = expr

function inject_jit(expr, jittedcallname, jittedarg, fixtypes)
    return rreplace(
        expr,
        e -> e isa Expr && e.head == :call && e.args[1] == jittedcallname,
        e -> jitexpr(e, jittedcallname, jittedarg, fixtypes)
    )
end

macro jit(jittedcallname::Symbol, jittedarg::Symbol, expr::Expr)
    @assert expr.head == :function "@jit only supports long-form function declarations"
    signature = expr.args[1]
    body = Expr(:block, expr.args[2].args...)
    explore = exploreexpr(jittedcallname, jittedarg)
    newbody = quote
        Base.@_inline_meta
        if ($explore) # A known, already jitted call
            callctx = Catwalk.callctx(jitctx, $jittedcallname)
            decoded_fixtypes = Catwalk.decode(Catwalk.fixtypes(callctx))
            return Catwalk.inject_jit(
                $(Expr(:quote, body)),
                $(Expr(:quote, jittedcallname)),
                $(Expr(:quote, jittedarg)),
                decoded_fixtypes
            )
        else
            return $(Expr(:quote, body))
        end
    end
    newfn = Expr(:macrocall, Symbol("@generated"), LineNumberNode(0), Expr(:function, signature, newbody))
    return newfn |> esc
end


text:
using SeisNetwork
using Test

@testset "SeisNetwork.jl" begin
    # Write your own tests here.
end


text:
using Pkg.Artifacts

using ClimateMachine.ArtifactWrappers

# Get dycoms dataset folder:
dycoms_dataset = ArtifactWrapper(
    joinpath(@__DIR__, "Artifacts.toml"),
    "dycoms",
    ArtifactFile[ArtifactFile(
        url = "https://caltech.box.com/shared/static/bxau6i46y6ikxn2sy9krgz0sw5vuptfo.nc",
        filename = "test_data_PhaseEquil.nc",
    ),],
)
dycoms_dataset_path = get_data_folder(dycoms_dataset)

@testset "Data tests" begin
    FT = Float64
    data = joinpath(dycoms_dataset_path, "test_data_PhaseEquil.nc")
    ds_PhaseEquil = Dataset(data, "r")
    e_int = Array{FT}(ds_PhaseEquil["e_int"][:])
    ρ = Array{FT}(ds_PhaseEquil["ρ"][:])
    q_tot = Array{FT}(ds_PhaseEquil["q_tot"][:])

    ts = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot, 4)
    # ts = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot, 3) # Fails
end


text:
export NN

"""
Neural Network structure
"""
mutable struct NN
    layers
    tmp
end

NN(layers=[SingleLayer();SingleLayer()]) = NN(layers,Array{AbstractArray}(undef,length(layers),2))

nLayers(N::NN) = length(N.layers)

"""
evaluate layer for current weights Θ=(K,b)
"""
function (N::NN)(S::AbstractArray{R},Θ) where R <: Real
    for k=1:nLayers(N)
        N.tmp[k,1] = S
        S = N.layers[k](S,Θ[k])
    end
    return S
end

"""
compute matvec J_S N(S,Θ)'*Z
"""
function getJSTmv(N::NN,Z::AbstractArray{R},S::AbstractArray{R},Θ) where R <: Real
    for k=nLayers(N):-1:1
        N.tmp[k,2] = Z
        Z = getJSTmv(N.layers[k],Z,N.tmp[k,1],Θ[k])
    end
    return Z
end

function getGradAndHessian(N::NN,dZ::AbstractArray{R},S::AbstractArray{R},Θ) where R <: Real
    N.tmp[end,2] = dZ
    dZ, d2Z = getGradAndHessian(N.layers[end],dZ,N.tmp[end,1],Θ[end])
    # dZ  = getJSTmv(N.layers[end],dZ,N.tmp[end],Θ[end])

    for k=nLayers(N)-1:-1:1
        N.tmp[k,2] = dZ
        dZ,d2Z = getGradAndHessian(N.layers[k],dZ,d2Z,N.tmp[k,1],Θ[k])
        # dZ  = getJSTmv(N.layers[k],dZ,N.tmp[k],Θ[k])
    end
    return dZ, d2Z
end

"""
compute matvec J_S(J_S N(S,Θ)'*Z(S))

here we use product rule

J_S N(S,Θ)'*dZ + J_S(N(S,Θ)'*Zfix)
"""
function getJSJSTmv(N::NN,dZ::AbstractArray{R},S::AbstractArray{R},Θ) where R <: Real
    d2Z = getJSJSTmv(N.layers[end],dZ,N.tmp[end,1],Θ[end])
    dZ  = getJSTmv(N.layers[end],dZ,N.tmp[end,1],Θ[end])

    for k=nLayers(N)-1:-1:1
        d2Z = getJSJSTmv(N.layers[k],dZ,d2Z,N.tmp[k,1],Θ[k])
        dZ  = getJSTmv(N.layers[k],dZ,N.tmp[k,1],Θ[k])
    end
    return d2Z,dZ
end

function getJSJSTmv(N::NN,dZ::AbstractArray{R},d2Z::AbstractArray{R},S::AbstractArray{R},Θ) where R <: Real
    for k=nLayers(N):-1:1
        d2Z = getJSJSTmv(N.layers[k],dZ,d2Z,N.tmp[k,1],Θ[k])
        # if k>1
            dZ  = getJSTmv(N.layers[k],dZ,N.tmp[k,1],Θ[k])
        # end
    end
    return d2Z,dZ
end

function getHessmv(N::NN,dZ::AbstractArray{R},S::AbstractArray{R},Θ,k=1) where R <: Real
    H1 = getHessmv(N.layers[k],N.tmp[k,2],dZ,N.tmp[k,1],Θ[k])
    if k < nLayers(N)
        dZ = getJSmv(N.layers[k],dZ,N.tmp[k,1],Θ[k])
        dZ = getHessmv(N,dZ,N.tmp[k+1,1],Θ,k+1)
        dZ = getJSTmv(N.layers[k],dZ,N.tmp[k,1],Θ[k])
        return H1 + dZ
    else
        return H1
    end
end

function getDiagHess(N::NN,dZ::AbstractArray{R},S::AbstractArray{R},Θ,k=1) where R <: Real
    H1 = getDiagHess(N.layers[k],N.tmp[k,2],dZ,N.tmp[k,1],Θ[k])
    if k < nLayers(N)
        dZ = getJSmv(N.layers[k],dZ,N.tmp[k,1],Θ[k])
        H2 = getDiagHess(N,dZ,N.tmp[k+1,1],Θ,k+1)
        # sum(dZ.*getHessmv(N,dZ,N.tmp[k+1,1],Θ,k+1),dims=1)
        # dZ = getJSTmv(N.layers[k],dZ,N.tmp[k,1],Θ[k])
        return H1 + H2
    else
        return H1
    end
end

function getTraceHess(N::NN,S::AbstractArray{R},Θ) where R <: Real
    trH1,Jac = getTraceHessAndGrad(N.layers[1],N.tmp[1,2],N.tmp[1,1],Θ[1])
    trH2 = getTraceHess(N,[],Jac,N.tmp[2,1],Θ,2)
    return trH1 + trH2
end

function getTraceHess(N::NN,w,Jac::AbstractArray{R},S::AbstractArray{R},Θ,k::Int=1) where R <: Real
    # FEATURE: second input, w, is not used here.
    if k < nLayers(N)
        trH1,Jac = getTraceHessAndGrad(N.layers[k],N.tmp[k,2],Jac,N.tmp[k,1],Θ[k])
        trH2,Jac = getTraceHessAndGrad(N,[],Jac,N.tmp[k+1,1],Θ,k+1)
        return trH1 + trH2
    else
        trH1 = getTraceHess(N.layers[k],N.tmp[k,2],Jac,N.tmp[k,1],Θ[k])
        return trH1
    end
end

function getTraceHessAndGrad(N::NN,S::AbstractArray{R},Θ) where R <: Real
    trH1,Jac = getTraceHessAndGrad(N.layers[1],N.tmp[1,2],N.tmp[1,1],Θ[1])
    trH2, Jac = getTraceHessAndGrad(N,[],Jac,N.tmp[2,1],Θ,2)
    return trH1 + trH2, Jac
end

function getTraceHessAndGrad(N::NN,w,Jac::AbstractArray{R},S::AbstractArray{R},Θ,k::Int=1) where R <: Real
    # FEATURE: second input, w, is not used here.
    trH1,Jac = getTraceHessAndGrad(N.layers[k],N.tmp[k,2],Jac,N.tmp[k,1],Θ[k])
    if k < nLayers(N)
        trH2,Jac = getTraceHessAndGrad(N,[],Jac,N.tmp[k+1,1],Θ,k+1)
        return trH1 + trH2, Jac
    else
        return trH1, Jac
    end
end


text:
using EconModels, Base.Test

###############################################
## 1. Define a model
###############################################

type TestModel <: ModelInterface
    spec::AbstractString
    subspec::AbstractString
    settings::Dict{Symbol,Setting}
    test_settings::Dict{Symbol,Setting}
    testing::Bool
    another_field::Vector{Symbol}
end

###############################################
## 2. Construct model
###############################################

# Define settings field
settings = Dict{Symbol, Setting}()
settings[:data_vintage] = Setting(:data_vintage, "170217", true, "vint", "Data vintage") # full constructor
settings[:use_parallel_workers] = Setting(:use_parallel_workers, false) # short constructor
tdir = tempdir()
settings[:saveroot] = Setting(:saveroot, tdir)
settings[:dataroot] = Setting(:dataroot, tdir)

# Construct test_settings field
test_settings = Dict{Symbol, Setting}()
test_settings[:data_vintage] = Setting(:data_vintage, "160217", true, "vint", "Test data vintage")

# Random other field
other_field = [:foo, :bar]

# Construct model
m = TestModel("m01", "ss1", settings, test_settings, false, other_field)

###############################################
## 3. Check that everything works
###############################################

# Getting non-test settings
@test get_setting(m, :data_vintage) == "170217"
@test get_setting(m, :use_parallel_workers) == false

# Getting test settings
m.testing = true
@test get_setting(m, :data_vintage) == "160217"
@test get_setting(m, :use_parallel_workers) == false
@test filestring(m) == "_test"

# Overwriting test settings
a = gensym() # unlikely to clash
b = gensym()
m <= Setting(a, 0, true, "abcd", "a")
m <= Setting(a, 1)
@test m.test_settings[a].value == 1
@test m.test_settings[a].print == true
@test m.test_settings[a].code == "abcd"
@test m.test_settings[a].description == "a"
m <= Setting(b, 2, false, "", "b")
m <= Setting(b, 3, true, "abcd", "b1")
@test m.test_settings[b].value == 3
@test m.test_settings[b].print == true
@test m.test_settings[b].code == "abcd"
@test m.test_settings[b].description == "b1"

# Builtin functions
@test spec(m) == "m01"
@test subspec(m) == "ss1"
@test m.another_field == other_field

# Test filestrings
m.testing = false
m <= Setting(:test, 5, true, "test", "a test setting")
@test m.settings[:test].value == 5
@test ismatch(r"^\s*_test=5_vint=(\d{6})", filestring(m))
@test ismatch(r"^\s*_key=val_test=5_vint=(\d{6})", filestring(m, "key=val"))
@test ismatch(r"^\s*_foo=bar_key=val_test=5_vint=(\d{6})", filestring(m, ["key=val", "foo=bar"]))
m.testing = true

# Model paths. all this should work without errors
m.testing = true
addl_strings = ["foo=bar", "hat=head", "city=newyork"]
for fn in [:rawpath, :workpath, :tablespath, :figurespath]
    @eval $(fn)(m, "test")
    @eval $(fn)(m, "test", "temp")
    @eval $(fn)(m, "test", "temp", addl_strings)
end