text:
# ---
# title: 1639. Number of Ways to Form a Target String Given a Dictionary
# id: problem1639
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/number-of-ways-to-form-a-target-string-given-a-dictionary/description/>
# hidden: true
# ---
#
# You are given a list of strings of the **same length** `words` and a string
# `target`.
#
# Your task is to form `target` using the given `words` under the following
# rules:
#
#   * `target` should be formed from left to right.
#   * To form the `ith` character ( **0-indexed** ) of `target`, you can choose the `kth` character of the `jth` string in `words` if `target[i] = words[j][k]`.
#   * Once you use the `kth` character of the `jth` string of `words`, you **can no longer** use the `xth` character of any string in `words` where `x <= k`. In other words, all characters to the left of or at index `k` become unusuable for every string.
#   * Repeat the process until you form the string `target`.
#
# **Notice**  that you can use **multiple characters** from the **same string**
# in `words` provided the conditions above are met.
#
# Return _the number of ways to form`target` from `words`_. Since the answer may
# be too large, return it **modulo** `109 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: words = ["acca","bbbb","caca"], target = "aba"
#     Output: 6
#     Explanation: There are 6 ways to form target.
#     "aba" -> index 0 (" _a_ cca"), index 1 ("b _b_ bb"), index 3 ("cac _a_ ")
#     "aba" -> index 0 (" _a_ cca"), index 2 ("bb _b_ b"), index 3 ("cac _a_ ")
#     "aba" -> index 0 (" _a_ cca"), index 1 ("b _b_ bb"), index 3 ("acc _a_ ")
#     "aba" -> index 0 (" _a_ cca"), index 2 ("bb _b_ b"), index 3 ("acc _a_ ")
#     "aba" -> index 1 ("c _a_ ca"), index 2 ("bb _b_ b"), index 3 ("acc _a_ ")
#     "aba" -> index 1 ("c _a_ ca"), index 2 ("bb _b_ b"), index 3 ("cac _a_ ")
#
#
# **Example 2:**
#
#
#
#     Input: words = ["abba","baab"], target = "bab"
#     Output: 4
#     Explanation: There are 4 ways to form target.
#     "bab" -> index 0 (" _b_ aab"), index 1 ("b _a_ ab"), index 2 ("ab _b_ a")
#     "bab" -> index 0 (" _b_ aab"), index 1 ("b _a_ ab"), index 3 ("baa _b_ ")
#     "bab" -> index 0 (" _b_ aab"), index 2 ("ba _a_ b"), index 3 ("baa _b_ ")
#     "bab" -> index 1 ("a _b_ ba"), index 2 ("ba _a_ b"), index 3 ("baa _b_ ")
#
#
# **Example 3:**
#
#
#
#     Input: words = ["abcd"], target = "abcd"
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: words = ["abab","baba","abba","baab"], target = "abba"
#     Output: 16
#
#
#
#
# **Constraints:**
#
#   * `1 <= words.length <= 1000`
#   * `1 <= words[i].length <= 1000`
#   * All strings in `words` have the same length.
#   * `1 <= target.length <= 1000`
#   * `words[i]` and `target` contain only lowercase English letters.
#
#
## @lc code=start
using LeetCode

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


text:
ipopt_optimizer =
    JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-6, "print_level" => 0)
fast_ipopt_optimizer = JuMP.optimizer_with_attributes(
    Ipopt.Optimizer,
    "print_level" => 0,
    "max_cpu_time" => 5.0,
)
# use default print_level = 5 # set to 0 to disable
GLPK_optimizer = JuMP.optimizer_with_attributes(GLPK.Optimizer)
Cbc_optimizer = JuMP.optimizer_with_attributes(Cbc.Optimizer)
OSQP_optimizer =
    JuMP.optimizer_with_attributes(OSQP.Optimizer, "verbose" => false, "max_iter" => 50000)
fast_lp_optimizer = JuMP.optimizer_with_attributes(Cbc.Optimizer, "seconds" => 3.0)
scs_solver = JuMP.optimizer_with_attributes(
    SCS.Optimizer,
    "max_iters" => 100000,
    "eps" => 1e-4,
    "verbose" => 0,
)

ext: jl


text:
function generateCompPlot( compReport::Dict{String,DataFrame},
    nodes::NTuple{N,String}, timeFactor::Float64, showPlot::Bool,
    savePlot::Bool, filename::String, plotFunction::Function ) where N

    fileroot = nothing
    extension = nothing

    # Get file name ready.
    if savePlot
        # If the user doesn't specify a filename, a random one is generated.
        if filename == ""
            dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
            fileroot = string( string( "Node composition plot (", dateStr, ")",
                labelInfix[:fluxType] ) )
            extension = "svg"
        else
            extension = split( filename, "." )
            fileroot = length( extension ) == 1 ? string( extension[1] ) :
                join( extension[1:(end - 1)], "." )
            extension = length( extension ) == 1 ? "" : extension[end]
            extension = extension ∈ extensions ? extension : "svg"
        end  # if filename == ""

        if !ispath( dirname( fileroot ) )
            mkpath( dirname( fileroot ) )
        end  # if !ispath( dirname( fileroot ) )
    end  # if savePlot

    # Generate plots.
    for node in unique( nodes )
        if !haskey( compReport, node )
            continue
        end  # if !haskey( compReport, node )

        plotfilename = string( fileroot, " ", node, ".", extension )
        plt = plotFunction( compReport[node], string( node ), timeFactor )

        # Show plot if needed.
        if showPlot
            gui( plt )
        end  # if showPlot

        # Save plot if needed.
        if savePlot
            savefig( plt, plotfilename )
        end  # if savePlot
    end  # for node in unique( nodes )

end  # generateCompPlot( compReport, nodes, timeFactor, showPlot, savePlot, fileroot, extension, plotFunction )

function normalCompPlot( compReport::DataFrame, node::String,
    timeFactor::Float64 )

    timePoints = compReport[:, :timePoint]
    labels = vcat( "Node total", string.( names( compReport )[3:end] ) )
    plotTitle = string( "Composition of node ", node )
    timePoints = compReport[:, :timePoint] / timeFactor
    dataPoints = Matrix( compReport[:, 2:end] )
    ymax = max( maximum( dataPoints[:, 1] ), 20 )

    plt = plot( timePoints, dataPoints[:, 1], show=false, size=(960, 540),
        title=plotTitle, label=labels[1], lw=3, color=:black, ylim=[0, ymax] )
    plot!( plt, timePoints, dataPoints[:, 2:end],
        label=hcat( labels[2:end]... ), lw=2 )

    return plt

end  # normalCompPlot( compReport, node, timeFactor )

function stackedCompPlot( compReport::DataFrame, node::String,
    timeFactor::Float64 )

    timePoints = compReport[:, :timePoint]
    labels = names( compReport )[3:end]
    plotTitle = string( "Composition of node ", node )
    timePoints = compReport[:, :timePoint] / timeFactor
    dataPoints = Matrix( compReport[:, 3:end] )
    nPoints, nSeries = size( dataPoints )
    dataPoints = hcat( fill( 0, nPoints ), cumsum( dataPoints, dims=2 ) )
    ymax = max( maximum( dataPoints[:, 1] ), 20 )

    plt = plot( show=false, size=(960,540), title=plotTitle, ylim=[0, ymax] )

    for ii in 1:nSeries
        plot!( plt, timePoints, dataPoints[:, ii + 1], lw=2, lalpha=1.0, fillto=dataPoints[:, ii], falpha=0.5, label=labels[ii] )
    end  # for ii in 1:nSeries

    return plt

end  # stackedCompPlot( compReport, node, timeFactor )


text:
export generate_positions, initialize_mesh, start, getstats, plotgraph

to_mW(dBm) = 10^(dBm / 10)

const ble_37_channel_wavelength = 0.12491352
const part_of_path_loss_for_37 = (ble_37_channel_wavelength / (4 * π))^2
const message_length = 312
const scanner_sensitivity = to_mW(-95)

function randexclude(range::UnitRange{T}, exclude::T) where T
    while true
        x = rand(range)
        x != exclude && return x
    end
end

function agent_step!(node::Node, model::AgentBasedModel)
    # produce a packet and store it in the queue for later
    if rand() < model.packet_emit_rate / model.n
        packet_seq = length(model.packets) + 1
        packet = Packet(packet_seq, node.id, randexclude(1:model.n, node.id), model.ttl, model.tick, 0, false)
        push!(model.packets, packet)

        # (pomdp): keeping track of nodes that will receive this packet
        model.packet_xs[packet.seq] = Int[]
        push!(node.pocket, (model.tick + rand(1:node.t_back_off_delay), packet_seq, model.ttl))
    end

    if node.state == :scanning
        node.event_start > model.tick && return
        node.channel = div(abs(model.tick - node.event_start), node.t_scan_interval) + 37

        # no packets acquired, repeat scanning
        if node.channel > 39
            node.channel = 37
            node.event_start = model.tick
        end

        # there is no more job for sinks here
        if node.role == :relay && node.packet_seq != 0 && !(node.packet_seq in node.received_packets)
            push!(node.received_packets, node.packet_seq)

            # our own packets are treated the same way as those generated by not us
            if node.packet_ttl > 1
                # put a packet in the pocket
                push!(node.pocket, (model.tick + rand(1:node.t_back_off_delay), node.packet_seq, node.packet_ttl - 1))
            end
        end
    end

    # wait until it's time to emit the first packet from the queue
    if !isempty(node.pocket) && model.tick >= node.pocket[1][1]
        node.state = :advertising
        node.event_start = model.tick
    end

    if node.state == :advertising
        node.event_start > model.tick && return
        if model.tick == node.event_start
            _, node.packet_seq, node.packet_ttl = popfirst!(node.pocket)
        end

        d, r = divrem(abs(model.tick - node.event_start), node.t_interpdu)

        # transmitting happens once for each channel
        node.transmitting = iszero(r)
        node.channel = d + 37

        # the end of one advertising cycle
        if node.channel > 39
            # bonus retransmissions for the original/retx packet
            if node.n_transmit_left > 0
                push!(node.pocket, (model.tick + node.t_back_off_delay, node.packet_seq, node.packet_ttl))
            end

            node.packet_seq = 0
            node.transmitting = false
            node.state = :scanning
            node.event_start = model.tick + 1
            node.n_transmit_left = node.n_retx_transmit_count
            node.channel = 37
        end
    end
end

# without sqrt, and taking into account overlapping devices
distance(a::Tuple{Int, Int}, b::Tuple{Int, Int}) = (a[1] - b[1]) ^ 2 + (a[2] - b[2]) ^ 2 + 1e-3

function calc_rssi!(agents, P_tx::Vector{Float64})
    n = length(agents)
    rssi_map = Matrix{Float64}(undef, n, n)

    L_p = Matrix{Float64}(undef, n, n)

    for i in 1:n
        for j in 1:(i - 1)
            L_p[i, j] = L_p[j, i]
        end

        for j in (i + 1):n
            L_p[i, j] = 1 / distance(agents[i].pos, agents[j].pos)
        end
    end

    for i in 1:n
        for j in 1:n
            rssi_map[i, j] = P_tx[j] * L_p[i, j] * part_of_path_loss_for_37 # dB
        end

        rssi_map[i, i] = 0
    end

    rssi_map
end

function calculate_par!(rssi_neighbourhood::Vector{Float64}, total_rssi::Float64)
    total_par = 0.0

    @inbounds for idx in eachindex(rssi_neighbourhood)
        rssi_neighbourhood[idx] == 0 && continue

        SINR = rssi_neighbourhood[idx] / (total_rssi - rssi_neighbourhood[idx])
        BER = 0.5 * erfc(sqrt(SINR / 2))
        PAR = (1 - BER)^message_length

        rssi_neighbourhood[idx] = PAR
        total_par += PAR
    end

    total_par
end

function select_rssi!(rssi_neighbourhood::Vector{Float64}, rssi_map::SubArray{Float64}, transmitters::CircularBuffer{Int})
    total_rssi = 0.0

    @inbounds for idx in transmitters
        rssi = rssi_map[idx]

        total_rssi += rssi
        if rssi > scanner_sensitivity
            rssi_neighbourhood[idx] = rssi
        end
    end

    total_rssi
end

function model_step!(model::AgentBasedModel)
    transmitters_count = 0

    for agent in allagents(model)
        if agent.transmitting
            push!(model.transmitters[agent.channel - 36], agent.id)
            transmitters_count += 1
        end
    end

    # (pomdp): clear previous reward nominees
    empty!(model.reward_plate)

    if iszero(transmitters_count)
        model.tick += 1
        return
    end

    rssi_neighbourhood = zeros(model.n)
    count_successes = zeros(model.n)

    for dst in allagents(model)
        dst.state == :scanning || continue

        # select RSSI of only transmitting devices reachable to here
        # begs the question whether RSSI can be calculated on the fly
        fill!(rssi_neighbourhood, 0)

        txs = model.transmitters[dst.channel - 36]
        total_rssi = select_rssi!(rssi_neighbourhood, view(model.rssi_map, :, dst.id), txs)

        iszero(total_rssi) && continue

        total_par = calculate_par!(rssi_neighbourhood, total_rssi)

        # do any of this packets pass?
        rand() <= total_par || continue

        tx_id = sample(1:length(txs), Weights(rssi_neighbourhood[1:length(txs)]))
        source_id = txs[tx_id]

        # @show model[source_id]
        packet = model.packets[model[source_id].packet_seq]
        count_successes[source_id] = 1

        if packet.dst == dst.id
            if !packet.done
                # since the packet is delivered, acknowledge everyone participating
                append!(model.reward_plate, model.packet_xs[packet.seq])

                # marking the packet as delivered
                packet.done = true
                packet.time_done = model.tick
            end
        else
            # register a touch upon this not delivered packet
            !packet.done && push!(model.packet_xs[packet.seq], dst.id)

            # copy packet's reference and source's ttl
            dst.packet_seq = packet.seq
            dst.packet_ttl = model[source_id].packet_ttl
            # dst.packet_rssi = neighbours_agents[src_i].packet_ttl
        end

        fill!(rssi_neighbourhood, 0.0)
    end

    if transmitters_count > 0
        model.packets_lost += transmitters_count - sum(count_successes)
        model.transmitters_count += transmitters_count
    end

    # reset transmitter signaling array
    for idx = 1:3
        empty!(model.transmitters[idx])
    end

    model.tick += 1
end

function plotgraph(model::AgentBasedModel)
    g = SimpleGraph(nagents(model))

    for (src_i, src) in enumerate(allagents(model)), (dst_i, dst) in enumerate(allagents(model))
        src.role == :relay || continue
        dst.role != :source || continue

        if model.rssi_map[src_i, dst_i] > model.scanner_sensitivity
            add_edge!(g, src.id, dst.id)
        end
    end

    colors = Dict(:relay => "orange", :sink => "darkblue", :source => "transparent")

    agents = sort(collect(allagents(model)), by=a -> a.id)
    filter!(a -> a.role != :source, agents)

    xs = map(a -> a.pos[1], agents)
    ys = map(a -> a.pos[2], agents)
    cs = map(a -> colors[a.role], agents)

    gplot(g, xs, ys; nodefillc = cs)
end

"""
    generate_positions(; dims = (200, 200), n = 100) → positions::Array{Tuple{Int,Int},1}

Generates random xy-coordinates for `n` nodes with `dims`
"""
generate_positions(; dims=(100, 100), n=64) = [(rand(1:dims[1]), rand(1:dims[2])) for _ in 1:n]

"""
    initialize_mesh(positions::Vector{Tuple{Int, Int}}, roles::Vector{Int}) → mesh::AgentBasedModel

Create a mesh by specifying roles of each node by `roles`=[0, 1, 0...] where 1 means that the node is a relay, sink otherwise, and including `positions`=[(Int, Int)...] as nodes' positions
"""
function initialize_mesh(positions::Vector{NTuple{2, Int}}, roles::AbstractVector)
    length(roles) == length(positions) || throw(ArgumentError("Both arguments must be of equal length"))
    n = length(roles)

    properties = Dict(
        :n => n,
        :tick => 0,
        :packet_emit_rate => 50 / 1000,
        :transmitters => [CircularBuffer{Int}(n) for _ = 1:3],
        :packets_lost => 0,
        :transmitters_count => 0,
        :ttl => 5,
        :rssi_map => Matrix{Float64}(undef, n, n),
        :tx_powers => Vector{Float64}(undef, n),
        :scanner_sensitivity => to_mW(-95),
        :shadow_d => LogNormal(0, 1),
        :multipath_d => Rayleigh(4),
        :wifi_noise_d => Normal((to_mW(-125) - to_mW(-135)) / 2 + to_mW(-135), (to_mW(-125) - to_mW(-135)) / 2),
        :gaussian_noise => Normal((to_mW(-125) - to_mW(-135)) / 2 + to_mW(-135), (to_mW(-125) - to_mW(-135)) / 2),
        :packets => Packet[],
        # packet.seq => nodes that have touched that packet
        :packet_xs => Dict{Int, Vector{Int}}(),
        # plate for the actors who contributed to packet's successful delivery
        :reward_plate => Int[]
    )

    space = GridSpace((maximum(first.(positions)), maximum(last.(positions))))
    model = ABM(Node, space; properties, warn = false)

    for i in 1:model.n
        node = Node(id = i; pos = positions[i], role = roles[i] == 1 ? :relay : :sink)

        add_agent_pos!(node, model)

        model.tx_powers[node.id] = node.tx_power
        model.tx_powers[i] = node.tx_power
    end

    model.rssi_map = calc_rssi!(collect(allagents(model)), model.tx_powers)
    model
end

function getstats(model::AgentBasedModel)
    # filter packets which were too recently produced, or produced in an unreachable position
    filter!(packet -> packet.time_start < model.tick - 500, model.packets)

    # worst case in terms of PDR to a single device
    worstnode = argmax(map(id -> count(p -> p.dst == id && p.done == false, model.packets), 1:nagents(model)))
    deprived = filter(p -> p.dst == worstnode, model.packets)
    delivered = filter(p -> p.done == true, model.packets)
    delays = map(p -> p.time_done - p.time_start, delivered)

    # probing centrality measure
    devices = collect(allagents(model))
    filter!(n -> n.role != :source, devices)

    g = SimpleDiGraph(length(devices))

    for (src_i, src) in enumerate(devices), (dst_i, dst) in enumerate(devices)
        src.role == :relay || continue

        if model.rssi_map[src_i, dst_i] > model.scanner_sensitivity
            add_edge!(g, src.id, dst.id)
        end
    end

    centrality = indegree_centrality(g, normalize=false) ./ length(devices) |> mean

    return ( pdr = length(delivered) / length(model.packets),
             worstpdr = count(p -> p.done == true, deprived) / length(deprived),
             delay = mean(delays),
             centrality = centrality,
             packetloss = model.packets_lost / model.transmitters_count )
end

"""
    start(model::AgentBasedModel, minutes = 1) → (received, produced)

Start an experiment from `model` for the given number of minutes.
"""
function start(model::AgentBasedModel; minutes = 1)
    steps = minutes * 60 * 1000

    run!(model, agent_step!, model_step!, steps)

    getstats(model)
end


text:
module InterestRates

using BusinessDays
using Dates

export
    ERF, ER, discountfactor, zero_rate, forward_rate,
    ERF_to_rate, discountfactor_to_rate,
    isnullcurve

include("splines.jl")
include("types.jl")
include("nullcurve.jl")
include("composite.jl")
include("daycount.jl")
include("discount.jl")
include("bufferedcurve.jl")

forward_rate(curve::AbstractIRCurve, forward_date::Date, maturity::Date) = ERF_to_rate(curve, ERF(curve, forward_date, maturity), yearfraction(curve_get_daycount(curve), forward_date, maturity))

# Let's use the curve's method to multiple-dispatch. Ugly methods _zero_rate are not exported.
zero_rate(curve::AbstractIRCurve, maturity::Date) = _zero_rate(curve_get_method(curve), curve, maturity)
zero_rate(curve::AbstractIRCurve, maturity_vec::Vector{Date}) = _zero_rate(curve_get_method(curve), curve, maturity_vec)

include("methods.jl") # implements various zero_rate() methods for each CurveMethod

end # module InterestRates


text:
"""
    imresize(img; ratio)

Return a copy of the image `img` scaled by `ratio`.

In case `ratio` is close to `1`, this function might also return
the original argument.

This is a replacement for a similar function from ImageTransformations.jl,
in order to avoid the loading time of that package. Currently this function might
not be as fancy as the one it replaces.
"""
function imresize(img; ratio::Float64)

    height, width = size(img)
    new_height = max(1, round(Int, ratio * height))
    new_width = max(1, round(Int, ratio * width))

    if (new_height == height && new_width == width)
        return img
    end

    new_img = similar(img, new_height, new_width)
    # TODO it would be nicer to use Lanczos interpolation
    # but that one creates negative colors at the moment
    itp = interpolate(img, BSpline(Linear()))

    # TODO StepRangeLen does not work for length == 1, maybe we can
    # solve this without type instability
    xs = (new_width <= 1) ? (1:1) : range(1, width, length=new_width)
    ys = (new_height <= 1) ? (1:1) : range(1, height, length=new_height)

    for (new_x, x) in enumerate(xs), (new_y, y) in enumerate(ys)
        new_img[new_y, new_x] = itp(y, x)
    end

    return new_img
end


text:
include("intervalbox.jl")
include("setdiff.jl")
include("arithmetic.jl")


text:
# https://software.intel.com/en-us/mkl-developer-reference-fortran-intel-mkl-pardiso-parameters-in-tabular-form

###############################################################
# MSGLVL: Pardiso verbosity
###############################################################
const MSGLVL_QUIET   = 0
const MSGLVL_VERBOSE = 1

###############################################################
# MTYPE: Pardiso matrix type
# This scalar value defines the matrix type. PARDISO supports the following matrices
###############################################################
# Real matrices
const MTYPE_REAL_STRUCTURALLY_SYMMETRIC         = 1
const MTYPE_REAL_SYMMETRIC_POSITIVE_DEFINITE    = 2
const MTYPE_REAL_SYMMETRIC_INDEFINITE           = -2
const MTYPE_REAL_NON_SYMMETRIC                  = 11
# Complex matrices
const MTYPE_COMPLEX_STRUCTURALLY_SYMMETRIC      = 3
const MTYPE_COMPLEX_HERMITIAN_POSITIVE_DEFINITE = 4
const MTYPE_COMPLEX_HERMITIAN_INDEFINITE        = -4
const MTYPE_COMPLEX_SYMMETRIC                   = 6
const MTYPE_COMPLEX_NON_SYMMETRIC               = 13

###############################################################
# PHASE: Controls the execution of the solver
#
# Usually it is a two- or three-digit integer.
# The first digit indicates the starting phase of execution and the second digit indicates the ending phase.
# Intel MKL PARDISO has the following phases of execution:
# 1. Phase 1: Fill-reduction analysis and symbolic factorization
# 2. Phase 2: Numerical factorization
# 3. Phase 3: Forward and Backward solve including iterative refinement.
#             This phase can be divided into two or three separate substitutions: forward, backward, and diagonal.
# 4. Termination and Memory Release Phase (PHASE ≤ 0)
###############################################################

const PHASE_ANALYSIS                                                    = 11
const PHASE_ANALYSIS_NUMERICAL_FACTORIZATION                            = 12
const PHASE_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT = 13
const PHASE_NUMERICAL_FACTORIZATION                                     = 22
const PHASE_SELECTED_INVERSION                                          = -22
const PHASE_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT          = 23
const PHASE_SOLVE_ITERATIVE_REFINEMENT                                  = 33
const PHASE_SOLVE_ITERATIVE_REFINEMENT_FORWARD_SUBSTITUTION             = 331
const PHASE_SOLVE_ITERATIVE_REFINEMENT_DIAGONAL_SUBSTITUTION            = 332
const PHASE_SOLVE_ITERATIVE_REFINEMENT_BACKWARD_SUBSTITUTION            = 333
const PHASE_RELEASE_INTERNAL_MEMORY                                     = 0
const PHASE_RELEASE_ALL_INTERNAL_MEMORY                                 = -1

"""
    pardiso_report_error(code::Int)
Report Pardiso error given its code.
"""
function pardiso_report_error(code::Int)

  if code < 0
    code == -1   &&
        error(string("Pardiso Error (", code, "): ", "Input inconsistent."))
    code == -2   &&
        error(string("Pardiso Error (", code, "): ", "Not enough memory."))
    code == -3   &&
        error(string("Pardiso Error (", code, "): ", "Reordering problem."))
    code == -4   &&
        error(string("Pardiso Error (", code, "): ",
            "Zero pivot, numerical fact. or iterative refinement problem."))
    code == -5   &&
        error(string("Pardiso Error (", code, "): ", "Unclassified (internal) error."))
    code == -6   &&
        error(string("Pardiso Error (", code, "): ", "Teordering failed (matrix types 11, 13 only)."))
    code == -7   &&
        error(string("Pardiso Error (", code, "): ", "Diagonal matrix is singular."))
    code == -8   &&
        error(string("Pardiso Error (", code, "): ", "32-bit integer overflow problem."))
    code == -10  &&
        error(string("Pardiso Error (", code, "): ", "Error opening OOC files."))
    code == -11  &&
        error(string("Pardiso Error (", code, "): ", "Read/write error with OOC files."))
    code == -12  &&
        error(string("Pardiso Error (", code, "): ", "pardiso_64 called from 32-bit library."))
    code == -13  &&
        error(string("Pardiso Error (", code, "): ", "Interrupted by the (user-defined) mkl_progress function."))
    code == -15  &&
        error(string("Pardiso Error (", code, "): ",
            "Internal error which can appear for iparm(24)=10 and iparm(13)=1. ",
            "Try switch matching off (set iparm(13)=0 and rerun.)."))

    error(string("Pardiso Error (", code, "): ", "Unknown error code."))
  end

end

###############################################################
# Parameter values for IPARM[12]
# IPARM_TRANSPOSED_OR_CONJUGATED_TRANSPOSED
###############################################################
const PARDISO_SOLVE_LINEAR_SYSTEM               = 0 # linear system.
const PARDISO_SOLVE_CONJUGATE_TRANSPOSED_SYSTEM = 1 # conjugate transposed system
const PARDISO_SOLVE_TRANSPOSED_SYSTEM           = 2 # transposed system

###############################################################
# Parameter values for IPARM[35]
# IPARM_ONE_OR_ZERO_BASED_INDEXING
###############################################################
const PARDISO_ONE_BASED_INDEXING                = 0 # One-based indexing
const PARDISO_ZERO_BASED_INDEXING               = 1 # Zero-based indexing

###############################################################
# Iparm parameters
#
# [I] -> input
# [O] -> output
# for iparm[i], where i is:
#
###############################################################
#  1:  [I]   Use default values
const IPARM_USE_DEFAULT_VALUES                                              = 1
#  2:  [I]   Fill-in reducing ordering for the input matrix.
const IPARM_FILL_IN_REDUCING_ORDERING                                       = 2
#  3:  [-]   Reserved. Set to zero.
#  4:  [I]   Preconditioned CGS/CG.
const IPARM_PRECONDITIONED_CGS/CG                                           = 4
#  5:  [I]   User permutation.
const IPARM_USER_PERMUTATION                                                = 5
#  6:  [I]   Write solution on x.
const IPARM_WRITE_SOLUTION_ON_X                                             = 6
#  7:  [O]   Number of iterative refinement steps performed.
const IPARM_NUMBER_ITERATIVE_REFINEMENT_STEPS                               = 7
#  8:  [I]   Iterative refinement step.
const IPARM_ITERATIVE_REFINEMENT_STEP                                       = 8
#  9:  [-]   Reserved. Set to zero.
#  10: [I]   Pivoting perturbation.
const IPARM_PIVOTING_PERTURBATION                                           = 10
#  11: [I]   Scaling vectors.
const IPARM_SCALING_VECTORS                                                 = 11
#  12: [I]   Solve with transposed or conjugate transposed matrix A.
const IPARM_TRANSPOSED_OR_CONJUGATED_TRANSPOSED                             = 12
#  13: [I]   Improved accuracy using (non-) symmetric weighted matching.
const IPARM_NON_SYMMETRIC_WEIGHTED_MATCHING                                 = 13
#  14: [O]   Number of perturbed pivots.
const IPARM_NUMBER_OF_PERTURBED_PIVOTS                                      = 14
#  15: [O]   Peak memory on symbolic factorization.
const IPARM_PEAK_MEMORY_ON_SYMBOLIC_FACTORIZATION                           = 15
#  16: [O]   Permanent memory on symbolic factorization.
const IPARM_PERMANENT_MEMORY_ON_SYMBOLIC_FACTORIZATION                      = 16
#  17: [O]   Size of factors/Peak memory on numerical factorization and solution.
const
IPARM_SIZE_OF_FACTORS_PEAK_MEMORY_ON_NUMERICAL_FACTORIZATION_AND_SOLUTION   = 17
#  18: [I/O] Report the number of non-zero elements in the factors.
const IPARM_REPORT_NUMBER_OF_NON_ZEROS_IN_FACTORS                           = 18
#  19: [I/O] Report number of floating point operations (in 106 floating point operations) that are necessary to factor the matrix A.
const IPARM_REPORT_NUMBER_OF_FLOATING_POINT_OPERATIONS                      = 19
#  20: [O]   Report CG/CGS diagnostics.
const IPARM_REPORT_CG_CGS_DIAGNOSTICS                                       = 20
#  21: [I]   Pivoting for symmetric indefinite matrices.
const IPARM_PRIVOTING_FOR_SYMMETRIC_INDEFINITE_MATRICES                     = 21
#  22: [O]   Inertia: number of positive eigenvalues.
const IPARM_NUMBER_OF_POSITIVE_EIGENVALUES                                  = 22
#  23: [O]   Inertia: number of negative eigenvalues.
const IPARM_NUMBER_OF_NEGATIVE_EIGENVALUES                                  = 23
#  24: [I]   Parallel factorization control.
const IPARM_PARALLEL_FACTORIZATION_CONTROL                                  = 24
#  25: [I]   Parallel forward/backward solve control.
const IPARM_PARALLEL_FORWARD_BACKWARD_SOLVE_CONTROL                         = 25
#  26: [-]   Reserved. Set to zero.
#  27: [I]   Matrix checker.
const IPARM_MATRIX_CHECKER                                                  = 27
#  28: [I]   Single or double precision Intel MKL PARDISO.
const IPARM_SINGLE_OR_DOUBLE_PRECISION                                      = 28
#  29: [-]   Reserved. Set to zero.
#  30: [O]   Number of zero or negative pivots.
const IPARM_NUMBER_OF_ZERO_OR_NEGATIVE_PIVOTS                               = 30
#  31: [I]   Partial solve and computing selected components of the solution vectors.
const IPARM_PARTIAL_SOLVE_AND_COMPUTING_SELECTED_COMPONENTS                 = 31
#  32: [-]   Reserved. Set to zero.
#  33: [-]   Reserved. Set to zero.
#  34: [I]   Optimal number of OpenMP threads for conditional numerical reproducibility (CNR) mode.
const IPARM_OPTIMAL_NUMBER_OF_OPENMPI_THREADS                               = 34
#  35: [I]   One- or zero-based indexing of columns and rows.
const IPARM_ONE_OR_ZERO_BASED_INDEXING                                      = 35
#  36: [I/O] Schur complement matrix computation control.
const IPARM_SCHUR_COMPLEMENT_MATRIX                                         = 36
#  37: [I]   Format for matrix storage.
const IPARM_FORMAT_FOR_MATRIX_STORAGE                                       = 37
#  38: [-]   Reserved. Set to zero.
#  39: [-]   Enable low rank update to accelerate factorization for multiple matrices with identical structure and similar values.
const IPARM_ENABLE_LOW_RANK                                                 = 39
#  40: [-]   Reserved. Set to zero.
#  41: [-]   Reserved. Set to zero.
#  42: [-]   Reserved. Set to zero.
#  43: [-]   Control parameter for the computation of the diagonal of inverse matrix.
const
IPARM_CONTROL_PARAMETER_FOR_COMPUTATION_OF_THE_DIAGONAL_OF_INVERSE_MATRIX   = 43
#  44: [-]   Reserved. Set to zero.
#  45: [-]   Reserved. Set to zero.
#  46: [-]   Reserved. Set to zero.
#  47: [-]   Reserved. Set to zero.
#  48: [-]   Reserved. Set to zero.
#  49: [-]   Reserved. Set to zero.
#  50: [-]   Reserved. Set to zero.
#  51: [-]   Reserved. Set to zero.
#  52: [-]   Reserved. Set to zero.
#  53: [-]   Reserved. Set to zero.
#  54: [-]   Reserved. Set to zero.
#  55: [-]   Reserved. Set to zero.
#  56: [-]   Diagonal and pivoting control.
#  57: [-]   Reserved. Set to zero.
#  58: [-]   Reserved. Set to zero.
#  59: [-]   Reserved. Set to zero.
#  60: [I]   Intel MKL PARDISO mode.
const IPARM_INTEL_MKL_PARDISO_MODE                                          = 60
#  61: [-]   Reserved. Set to zero.
#  62: [-]   Reserved. Set to zero.
#  63: [O]   Size of the minimum OOC memory for numerical factorization and solution.
const IPARM_SIZE_OF_THE_MINIMUM_OCC_MEMORY                                  = 63
#  64: [-]   Reserved. Set to zero.


text:
module GlnModel

using JuMP
using KNITRO

function gln_model(;
    glc_ext = 0.05,
    lac_ext = 1e-4,
    nh4_ext = 0.01,
    ac_ext = 1e-6,
    etoh_ext = 1e-6,
    formate_ext = 1e-4,
    proteome = 0.05,
    num_ms = 20,
    atp_adp_ratio = -1,
    nadh_nad_ratio = -1,
    silence=true,
)

    gln_model = Model(
        optimizer_with_attributes(
            KNITRO.Optimizer,
            "ms_enable" => 1,
            "ms_maxsolves" => num_ms,
        ),
    )
    if silence
        JuMP.set_silent(gln_model)
    end
    # Thermodynamic function
    thermo_factor(x) = tanh(-0.7 * x) # x = ΔrG/RT
    register(gln_model, :thermo_factor, 1, thermo_factor; autodiff = true)

    # Bounds
    LB_enz = 0.0 # [g enz / g DW cell]
    UB_enz = proteome # [g enz / g DW cell]
    LB_met = log(1e-9) # [log(M)]
    UB_met = log(100e-3) # [log(M)]
    LB_v = -1 # [mmol/gDW/h]
    UB_v = 100 # [mmol/gDW/h]
    LB_dg = -80 # [kJ/mol]
    UB_dg = 1 # [kJ/mol]
    mf = 4.0 # metabolite factor

    @variables gln_model begin
        0 <= mu <= UB_v # gln production rate

        # atp burn (excess)
        0 <= v_burn <= UB_v # mmol/gDW/h

        # Enzymes concentrations [g enz / g DW]
        LB_enz <= pts <= UB_enz
        LB_enz <= pgi <= UB_enz
        LB_enz <= pfk <= UB_enz
        LB_enz <= fba <= UB_enz
        LB_enz <= tpi <= UB_enz
        LB_enz <= gapd <= UB_enz
        LB_enz <= pgk <= UB_enz
        LB_enz <= pgm <= UB_enz
        LB_enz <= eno <= UB_enz
        LB_enz <= pyk <= UB_enz
        LB_enz <= ppc <= UB_enz
        LB_enz <= ldh <= UB_enz
        LB_enz <= pdh <= UB_enz
        LB_enz <= lact <= UB_enz
        LB_enz <= cs <= UB_enz
        LB_enz <= aconta <= UB_enz
        LB_enz <= acontb <= UB_enz
        LB_enz <= icdh <= UB_enz
        LB_enz <= gludy <= UB_enz
        LB_enz <= glns <= UB_enz
        LB_enz <= nh4t <= UB_enz
        LB_enz <= ptar <= UB_enz
        LB_enz <= ackr <= UB_enz
        LB_enz <= act <= UB_enz
        LB_enz <= acald <= UB_enz
        LB_enz <= alcd <= UB_enz
        LB_enz <= etoht <= UB_enz
        LB_enz <= pfl <= UB_enz
        LB_enz <= fort <= UB_enz

        # Metabolite concentrations [log(M)]
        LB_met <= g6p <= UB_met
        LB_met <= f6p <= UB_met
        LB_met <= dpg13 <= UB_met
        LB_met <= pg3 <= UB_met
        LB_met <= pg2 <= UB_met
        LB_met <= pyr <= UB_met
        LB_met <= lac <= UB_met
        LB_met <= oaa <= UB_met
        LB_met <= acon <= UB_met
        LB_met <= icit <= UB_met
        LB_met <= nh4 <= UB_met
        LB_met <= acetald <= UB_met
        LB_met <= etoh <= UB_met
        LB_met <= actp <= UB_met
        LB_met <= ac <= UB_met
        LB_met <= formate <= UB_met

        LB_met <= pep <= UB_met
        # log(1/mf * 1.8e-4) <= pep <= log(mf * 1.8e-4)

        LB_met <= dhap <= UB_met
        # log(1/mf * 3.7e-4) <= dhap <= log(mf * 3.7e-4)

        LB_met <= akg <= UB_met
        # log(1/mf * 4.4e-4) <= akg <= log(mf * 4.4e-4)

        LB_met <= accoa <= UB_met
        # log(1/mf * 6.1e-4) <= accoa <= log(mf * 6.1e-4)

        LB_met <= g3p <= UB_met
        # log(1/mf * 1.5e-3) <= g3p <= log(mf * 1.5e-3)

        LB_met <= cit <= UB_met
        # log(1/mf * 2.0e-3) <= cit <= log(mf * 2.0e-3)

        LB_met <= gln <= UB_met
        # log(1/mf * 3.8e-3) <= gln <= log(mf * 3.8e-3)

        LB_met <= fdp <= UB_met
        # log(1/mf * 1.5e-2) <= fdp <= log(mf * 1.5e-2)

        LB_met <= glu <= UB_met
        # log(1/mf * 9.6e-2) <= glu <= log(mf * 9.6e-2)

        # cofactors
        LB_met <= adp <= UB_met
        # log(mf * 5.6e-4) <= adp <= log(mf * 5.6e-4)
        LB_met <= atp <= UB_met
        # log(1/mf * 9.6e-3) <= atp <= log(mf * 9.6e-3)
        LB_met <= nad <= UB_met
        # log(1/mf * 2.6e-3) <= nad <= log(mf * 2.6e-3)
        LB_met <= nadh <= UB_met
        # log(mf * 8.3e-5) <= nadh <= log(mf * 8.3e-5)

        # Fluxes [mmol/gDW/h]
        LB_v <= v_pts <= UB_v
        LB_v <= v_nh4t <= UB_v
        LB_v <= v_lact <= UB_v
        LB_v <= v_act <= UB_v
        LB_v <= v_etoht <= UB_v

        LB_v <= v_pgi <= UB_v
        LB_v <= v_pfk <= UB_v
        LB_v <= v_fba <= UB_v
        LB_v <= v_tpi <= UB_v
        LB_v <= v_gapd <= UB_v
        LB_v <= v_pgk <= UB_v
        LB_v <= v_pgm <= UB_v
        LB_v <= v_pyk <= UB_v
        LB_v <= v_eno <= UB_v
        LB_v <= v_ppc <= UB_v
        LB_v <= v_ldh <= UB_v
        LB_v <= v_pdh <= UB_v
        LB_v <= v_cs <= UB_v
        LB_v <= v_aconta <= UB_v
        LB_v <= v_acontb <= UB_v
        LB_v <= v_icdh <= UB_v
        LB_v <= v_gludy <= UB_v
        LB_v <= v_glns <= UB_v
        LB_v <= v_ptar <= UB_v
        LB_v <= v_ackr <= UB_v
        LB_v <= v_acald <= UB_v
        LB_v <= v_alcd <= UB_v
        LB_v <= v_pfl <= UB_v
        LB_v <= v_fort <= UB_v

        # Thermodynamic variables
        LB_dg <= dg_pts <= UB_dg
        LB_dg <= dg_nh4t <= UB_dg
        LB_dg <= dg_lact <= UB_dg
        LB_dg <= dg_etoht <= UB_dg
        LB_dg <= dg_act <= UB_dg

        LB_dg <= dg_pgi <= UB_dg
        LB_dg <= dg_pfk <= UB_dg
        LB_dg <= dg_fba <= UB_dg
        LB_dg <= dg_tpi <= UB_dg
        LB_dg <= dg_gapd <= UB_dg
        LB_dg <= dg_pgk <= UB_dg
        LB_dg <= dg_pgm <= UB_dg
        LB_dg <= dg_pyk <= UB_dg
        LB_dg <= dg_eno <= UB_dg
        LB_dg <= dg_ppc <= UB_dg
        LB_dg <= dg_ldh <= UB_dg
        LB_dg <= dg_pdh <= UB_dg
        LB_dg <= dg_cs <= UB_dg
        LB_dg <= dg_aconta <= UB_dg
        LB_dg <= dg_acontb <= UB_dg
        LB_dg <= dg_icdh <= UB_dg
        LB_dg <= dg_gludy <= UB_dg
        LB_dg <= dg_glns <= UB_dg
        LB_dg <= dg_ptar <= UB_dg
        LB_dg <= dg_ackr <= UB_dg
        LB_dg <= dg_acald <= UB_dg
        LB_dg <= dg_alcd <= UB_dg
        LB_dg <= dg_pfl <= UB_dg
        LB_dg <= dg_fort <= UB_dg
    end

    @NLparameters gln_model begin
        RT == 8.3145e-3 * 298.15 # [kJ/mol]

        # Enzyme rates [mmol/g/h]
        kcat_pts == 7240.319618601774
        kcat_pgi == 7602.809491918478
        kcat_pfk == 29779.773325683316
        kcat_fba == 12070.317665674209
        kcat_tpi == 25881.04927925183
        kcat_gapd == 7096.06140643298
        kcat_pgk == 10031.027221032862
        kcat_pgm == 31309.78286506871
        kcat_pyk == 8271.195139265732
        kcat_eno == 5415.070665517126
        kcat_ppc == 1671.6634868719905
        kcat_ldh == 1604.7523387985657
        kcat_pdh == 7349.716396809779
        kcat_lact == 3963.1636080052076
        kcat_cs == 2361.1129054755233
        kcat_aconta == 1998.844393316627
        kcat_acontb == 1672.316435538672
        kcat_icdh == 865.4413532355705
        kcat_gludy == 13709.063214013708
        kcat_glns == 1708.433754802496
        kcat_nh4t == 6173.782075500597
        kcat_ptar == 23851.29662502135
        kcat_ackr == 26834.756823345528
        kcat_act == 7090.909090909091
        kcat_acald == 11454.251270331857
        kcat_alcd == 1622.4755081413768
        kcat_etoht == 7090.909090909091
        kcat_pfl == 1533.4864896342106
        kcat_fort == 7603.197337814887

        # Gibbs energy of reaction
        dg0_pts == -2.9
        dg0_nh4t == 0.0 # diffusion
        dg0_etoht == 4.0 # proton symport
        dg0_act == 4.0 # proton symport
        dg0_lact == 4.0 # proton symport
        dg0_fort == 4.0 # proton symport

        dg0_pgi == 0.9
        dg0_pfk == -14.7
        dg0_fba == 25.5
        dg0_tpi == 6.1
        dg0_gapd == 2.3
        dg0_pgk == -17.8
        dg0_pgm == -2.85
        dg0_pyk == -26.8
        dg0_eno == -4.3
        dg0_ppc == -38.7
        dg0_ldh == -23.15
        dg0_pdh == -35.9
        dg0_cs == -37.9
        dg0_aconta == -12.04
        dg0_acontb == 1.9
        dg0_icdh == 3.25
        dg0_gludy == -37.7
        dg0_glns == -16.1
        dg0_acald == 21.7
        dg0_ptar == 14.9
        dg0_ackr == -16.2
        dg0_alcd == -20.9
        dg0_pfl == -23.4

        # media conditions log[M]
        glc_e == log(glc_ext)
        lac_e == log(lac_ext)
        co2 == log(1e-4)
        nh4_e == log(nh4_ext)
        etoh_e == log(etoh_ext)
        ac_e == log(ac_ext)
        formate_e == log(formate_ext)

        # intracellular conditions log[M]
        phos == log(1e-3)
        nadph == log(1.2e-4)
        nadp == log(2.1e-6)
        coa == log(1.4e-3)
    end

    @NLconstraints gln_model begin
        # ΔG at concentrations
        dg_pts == dg0_pts + RT * (g6p + pyr - pep - glc_e)
        dg_nh4t == dg0_nh4t + RT * (nh4 - nh4_e)
        dg_lact == dg0_lact + RT * (lac_e - lac)
        dg_etoht == dg0_etoht + RT * (etoh_e - etoh)
        dg_act == dg0_act + RT * (ac_e - ac)
        dg_fort == dg0_fort + RT * (formate_e - formate)

        dg_pgi == dg0_pgi + RT * (f6p - g6p)
        dg_pfk == dg0_pfk + RT * (fdp + adp - atp - f6p)
        dg_fba == dg0_fba + RT * (g3p + dhap - fdp)
        dg_tpi == dg0_tpi + RT * (g3p - dhap)
        dg_gapd == dg0_gapd + RT * (dpg13 + nadh - nad - phos - g3p)
        dg_pgk == dg0_pgk + RT * (pg3 + atp - adp - dpg13)
        dg_pgm == dg0_pgm + RT * (pg2 - pg3)
        dg_eno == dg0_eno + RT * (pep - pg2)
        dg_pyk == dg0_pyk + RT * (pyr + atp - adp - pep)
        dg_ppc == dg0_ppc + RT * (oaa + phos - co2 - pep)
        dg_ldh == dg0_ldh + RT * (lac + nad - nadh - pyr)
        dg_pdh == dg0_pdh + RT * (accoa + nadh + co2 - nad - coa - pyr)
        dg_cs == dg0_cs + RT * (cit + coa - oaa)
        dg_aconta == dg0_aconta + RT * (acon - cit)
        dg_acontb == dg0_acontb + RT * (icit - acon)
        dg_icdh == dg0_icdh + RT * (akg + nadph + co2 - nadp - icit)
        dg_gludy == dg0_gludy + RT * (glu + nadp - nadph - nh4 - akg)
        dg_glns == dg0_glns + RT * (gln + phos + adp - atp - nh4 - glu)
        dg_acald == dg0_acald + RT * (acetald + nad - nadh - accoa)
        dg_alcd == dg0_alcd + RT * (etoh + nad - nadh - acetald)
        dg_ptar == dg0_ptar + RT * (actp + coa - phos - accoa)
        dg_ackr == dg0_ackr + RT * (ac + atp - adp - actp)
        dg_pfl == dg0_pfl + RT * (accoa + formate - coa - pyr)

        # flux to kinetics and thermodynamics
        v_pts == kcat_pts * pts * thermo_factor(dg_pts / RT)
        v_lact == kcat_lact * lact * thermo_factor(dg_lact / RT)
        v_nh4t == kcat_nh4t * nh4t * thermo_factor(dg_nh4t / RT)
        v_etoht == kcat_etoht * etoht * thermo_factor(dg_etoht / RT)
        v_act == kcat_act * act * thermo_factor(dg_act / RT)
        v_fort == kcat_fort * fort * thermo_factor(dg_fort / RT)

        v_pgi == kcat_pgi * pgi * thermo_factor(dg_pgi / RT)
        v_pfk == kcat_pfk * pfk * thermo_factor(dg_pfk / RT)
        v_fba == kcat_fba * fba * thermo_factor(dg_fba / RT)
        v_tpi == kcat_tpi * tpi * thermo_factor(dg_tpi / RT)
        v_gapd == kcat_gapd * gapd * thermo_factor(dg_gapd / RT)
        v_pgk == kcat_pgk * pgk * thermo_factor(dg_pgk / RT)
        v_pgm == kcat_pgm * pgm * thermo_factor(dg_pgm / RT)
        v_eno == kcat_eno * eno * thermo_factor(dg_eno / RT)
        v_pyk == kcat_pyk * pyk * thermo_factor(dg_pyk / RT)
        v_ppc == kcat_ppc * ppc * thermo_factor(dg_ppc / RT)
        v_ldh == kcat_ldh * ldh * thermo_factor(dg_ldh / RT)
        v_pdh == kcat_pdh * pdh * thermo_factor(dg_pdh / RT)
        v_cs == kcat_cs * cs * thermo_factor(dg_cs / RT)
        v_aconta == kcat_aconta * aconta * thermo_factor(dg_aconta / RT)
        v_acontb == kcat_acontb * acontb * thermo_factor(dg_acontb / RT)
        v_icdh == kcat_icdh * icdh * thermo_factor(dg_icdh / RT)
        v_gludy == kcat_gludy * gludy * thermo_factor(dg_gludy / RT)
        v_acald == kcat_acald * acald * thermo_factor(dg_acald / RT)
        v_alcd == kcat_alcd * alcd * thermo_factor(dg_alcd / RT)
        v_ptar == kcat_ptar * ptar * thermo_factor(dg_ptar / RT)
        v_ackr == kcat_ackr * ackr * thermo_factor(dg_ackr / RT)
        v_pfl == kcat_pfl * pfl * thermo_factor(dg_pfl / RT)

        # mass balance constraints
        v_pts - v_pgi == 0 # g6p
        v_pgi - v_pfk == 0 # f6p
        v_pfk - v_fba == 0 # fdp
        v_fba - v_tpi == 0 # dhap
        v_fba + v_tpi - v_gapd == 0 # g3p
        v_gapd - v_pgk == 0 # dpg13
        v_pgk - v_pgm == 0 # pg3
        v_pgm - v_eno == 0 # pg2
        v_eno - v_pyk - v_ppc == 0 # pep
        v_pyk - v_ldh - v_pdh - v_pfl == 0 # pyr
        v_ldh - v_lact == 0 # lac
        v_pdh - v_cs - v_acald - v_ptar + v_pfl == 0 # accoa
        v_ppc - v_cs == 0 # oaa
        v_cs - v_aconta == 0 # cit
        v_aconta - v_acontb == 0 # acon
        v_acontb - v_icdh == 0 # icit
        v_icdh - v_gludy == 0 # akg
        v_nh4t - v_gludy - v_glns == 0 # nh4
        v_gludy - v_glns == 0 # glu
        v_glns - mu == 0 # gln
        v_pgk + v_pyk - v_pfk - v_glns - v_burn + v_ackr == 0 # atp
        -v_pgk - v_pyk + v_pfk + v_glns + v_burn - v_ackr == 0 # adp
        v_ldh - v_pdh - v_gapd + v_acald + v_alcd == 0 # nad
        -v_ldh + v_pdh + v_gapd - v_acald - v_alcd == 0 # nadh
        v_ptar - v_ackr == 0 # actp
        v_ackr - v_act == 0 # ac
        v_acald - v_alcd == 0 # acetald
        v_alcd - v_etoht == 0 # etoh
        v_pfl - v_fort == 0 # formate

        # density constraint(s) (can add more, e.g. membrane)
        pts +
        pgi +
        pfk +
        fba +
        tpi +
        gapd +
        pgk +
        pgm +
        eno +
        pyk +
        ppc +
        ldh +
        pdh +
        lact +
        cs +
        aconta +
        acontb +
        icdh +
        gludy +
        glns +
        nh4t +
        ptar +
        ackr +
        acald +
        alcd +
        etoht +
        act +
        pfl +
        fort <= proteome
    end

    if atp_adp_ratio > 0
        @constraint(gln_model, atp == log(atp_adp_ratio) + adp)
    end

    if nadh_nad_ratio > 0
        @constraint(gln_model, nadh == log(nadh_nad_ratio) + nad)
    end

    return gln_model
end

function max_mu!(gln_model, prev_sol=Dict{Symbol, Float64}())

    if !isempty(prev_sol)
        for (k, v) in prev_sol
            set_start_value(gln_model[k], v)
        end
    end

    @objective(gln_model, Max, gln_model[:mu])

    optimize!(gln_model)

    if termination_status(gln_model) ∉ [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_LOCALLY_SOLVED, MOI.ALMOST_OPTIMAL]
        return false
    else
        return true
    end
end

function find_nearest!(gln_model, prev_sol=Dict{Symbol, Float64}())

    isempty(prev_sol) && return true # if no prev solution, don't do anything

    if !isempty(prev_sol)
        for (k, v) in prev_sol
            set_start_value(gln_model[k], v)
        end
    end

    mu = value(gln_model[:mu])

    JuMP.fix(gln_model[:mu], mu; force=true)

    @objective(gln_model, Min,
        sum( (gln_model[k] - v)^2 for (k, v) in prev_sol)
    )

    optimize!(gln_model)

    if termination_status(gln_model) ∉ [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_LOCALLY_SOLVED, MOI.ALMOST_OPTIMAL]
        return false
    else
        return true
    end
end

function update!(gln_model, prev_sol)
    for var in Symbol.(all_variables(gln_model))
        prev_sol[var] = value(gln_model[var])
    end
end

end #module


text:
#=
Functionalities to generate a sandbox module.
=#

"""
$SIGNATURES

Return a sandbox module name corresponding to the page at `fpath`, effectively
`FD_SANDBOX_*` where `*` is a hash of the path.
"""
modulename(fpath::AS) = "FD_SANDBOX_$(hash(fpath))"

"""
$SIGNATURES

Checks whether a name is a defined module.
"""
function ismodule(name::String)::Bool
    s = Symbol(name)
    isdefined(Main, s) || return false
    typeof(getfield(Main, s)) === Module
end

"""
$SIGNATURES

Creates a new module with a given name, if the module exists, it is wiped.
Discards the warning message that a module is replaced which may otherwise
happen. Return a handle pointing to the module.
"""
function newmodule(name::String)::Module
    mod  = nothing
    junk = tempname()

    open(junk, "w") do outf
        # discard the "WARNING: redefining module X"
        redirect_stderr(outf) do
            mod = Core.eval(Main, Meta.parse("""
                module $name
                    import Franklin
                    import Franklin: @OUTPUT, @delay, fdplotly,
                                     locvar, pagevar, globvar,
                                     fd2html, get_url
                    """ * ifelse(startswith(name, "Utils"), "",
                            ifelse(FD_ENV[:UTILS_HASH] === nothing, "", """
                    import $(Franklin.utils_name())
                    """)) * """
                    using Dates
                end
                """))
        end
    end
    return mod
end


text:
function train_don()
    # if has_cuda()
    #     @info "CUDA is on"
    #     device = gpu
    #     CUDA.allowscalar(false)
    # else
        device = cpu
    # end

    x, y = get_data_don(n=300)
    xtrain = x[1:280, :]' |> device
    xval = x[end-19:end, :]' |> device

    ytrain = y[1:280, :] |> device
    yval = y[end-19:end, :] |> device

    grid = collect(range(0, 1, length=1024))' |> device

    learning_rate = 0.001
    opt = ADAM(learning_rate)

    m = DeepONet((1024,1024,1024),(1,1024,1024),gelu,gelu)
    loss(xtrain,ytrain,sensor) = Flux.Losses.mse(m(xtrain,sensor),ytrain)
    evalcb() = @show(loss(xval,yval,grid))

    Flux.@epochs 400 Flux.train!(loss, params(m), [(xtrain,ytrain,grid)], opt, cb = evalcb)
    ỹ = m(xval, grid)

    diffvec = vec(abs.((yval .- ỹ)))
    mean_diff = sum(diffvec)/length(diffvec)
    return mean_diff
end


text:
"""
	struct FFTHolder{FP,IP}

A container for pre-allocated FFT grids and plans for a provided grid-size.
"""
struct FFTHolder{FP,IP}
    dims   :: NTuple{3,Int}
    rgrids :: Vector{FieldGrid{Float64}}
    kgrids :: Vector{FieldGrid{Complex{Float64}}}

    FT     :: FP
    iFT    :: IP
	# FT     :: FFTW.rFFTWPlan{Float64,-1,false,3}
    # iFT    :: AbstractFFTs.ScaledPlan{Complex{Float64},FFTW.rFFTWPlan{Complex{Float64},1,false,3},Float64}
end

function FFTHolder(dims::NTuple{3,<:Integer}; nthreads::Integer = -1)
	gr = zeros(Float64, dims)
	gk = zeros(Complex{Float64}, ksize(dims))

	TGR = typeof(gr)
	TGK = typeof(gk)

	# Setup multi-threads and FFT plan
	if nthreads > 0
		FFTW.set_num_threads(nthreads)
	else
		nthreads = determine_nthreads(prod(dims))
		FFTW.set_num_threads(nthreads)
	end
	FT = plan_rfft(gr; flags = FFTW.MEASURE)
	iFT = plan_irfft(gk, dims[1]; flags = FFTW.MEASURE) #inv(FT)
	plan = FFTHolder(dims, TGR[], TGK[], FT, iFT)

	return plan
end

FFTHolder(nx::Integer, ny::Integer, nz::Integer; kwargs...) = FFTHolder((nx, ny, nz); kwargs...)

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

Base.show(io::IO, plan::FFTHolder) = @printf(io, "%s(dims = %s)", typeof(plan), plan.dims)

Base.empty!(plan::FFTHolder) = (empty!(plan.rgrids); empty!(plan.kgrids);)

nrgrids(plan::FFTHolder) = length(plan.rgrids)
nkgrids(plan::FFTHolder) = length(plan.kgrids)

add_rgrid!(plan::FFTHolder) = push!(plan.rgrids, zeros(Float64, plan.dims))
add_kgrid!(plan::FFTHolder) = push!(plan.kgrids, zeros(Complex{Float64}, ksize(plan.dims)))

function rgrid(plan::FFTHolder, idx::Integer)
	while nrgrids(plan) < idx
		add_rgrid!(plan)
	end
	return plan.rgrids[idx]
end

function kgrid(plan::FFTHolder, idx::Integer)
	while nkgrids(plan) < idx
		add_kgrid!(plan)
	end
	return plan.kgrids[idx]
end


text:
"""
    reglin(X,Y)

adalah fungsi untuk mencari garis linear dari variabel dependen `Y` dan
independen `X` dengan persamaan `y=Ax+B`.

# Example
```jldoctest
julia> x = [-1, 0, 1, 2, 3, 4, 5, 6];

julia> y = [10, 9, 7, 5, 4, 3, 0, -1];

julia> A,B = reglin(x,y)
(-1.6071428571428572, 8.642857142857142)
```
return nilai koefisien `A` dan `B`.
"""
function reglin(X,Y)
    xmean = mean(X);
    ymean = mean(Y);
    # Hitung nilai jumlah dari xy dan x^2
    sumxy = (X.-xmean)'*(Y.-ymean)
    sumx2 = (X.-xmean)'*(X.-xmean)
    # Hitung nilau koefisien garis regresi linear Y=Ax+B
    A = sumxy/sumx2;
    B = ymean .- A*xmean;
    return A,B
end