text:
include("common_micro_bm.jl")

# JET does not like KernelAbstractions.@print calls,
# so we disable them here.
TD.print_warning() = false

function jet_thermo_states()
    ArrayType = Array{Float64}
    profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)

    @testset "JET tests" begin
        for C in (
            TD.PhaseEquil_ρeq,
            TD.PhaseEquil_ρTq,
            TD.PhaseEquil_pTq,
            TD.PhaseEquil_peq,
            TD.PhaseEquil_phq,
            TD.PhaseEquil_ρθq,
            TD.PhaseEquil_pθq,
            TD.PhaseEquil_ρpq,
        )
            for cond in conditions(C)
                args = sample_args(profiles, cond, C)
                JET.@test_opt C(param_set, args...)
            end
        end
    end
end

jet_thermo_states()


text:
function chwrite(sock::ClickHouseSock, x::VarUInt)
    mx::UInt64 = x
    while mx >= 0x80
        write(sock.io, UInt8(mx & 0xFF) | 0x80)
        mx >>= 7
    end
    write(sock.io, UInt8(mx & 0xFF))
end

function chread(sock::ClickHouseSock, ::Type{VarUInt})::VarUInt
    x::UInt64 = 0
    s::UInt32 = 0
    i::UInt64 = 0
    while true
        b = read(sock.io, UInt8)
        if b < 0x80
            if i > 9 || (i == 9 && b > 1)
                throw(OverflowError("varint would overflow"))
            end
            return x | UInt64(b) << s
        end

        x |= UInt64(b & 0x7F) << s
        s += 7
        i += 1
    end
end

chread(sock::ClickHouseSock, ::Type{T}) where T <: Number =
    read(sock.io, T)

chread(sock::ClickHouseSock, x::UInt64)::Vector{UInt8} =
    read(sock.io, x)

function chread(sock::ClickHouseSock, ::Type{String})::String
    len = chread(sock, VarUInt) |> UInt64
    chread(sock, len) |> String
end

# Vector reads
function chread(
    sock::ClickHouseSock,
    ::Type{Vector{T}},
    count::VarUInt,
)::Vector{T} where T <: Number
    data = Vector{T}(undef, UInt64(count))
    read!(sock.io, data)
    data
end

chread(
    sock::ClickHouseSock,
    ::Type{Vector{String}},
    count::VarUInt,
)::Vector{String} = [chread(sock, String) for _ ∈ 1:UInt64(count)]

# Scalar writes
chwrite(sock::ClickHouseSock, x::Number) =
    write(sock.io, x)

function chwrite(sock::ClickHouseSock, x::String)
    chwrite(sock, x |> sizeof |> VarUInt)
    chwrite(sock, x |> Array{UInt8})
end

# Vector writes
chwrite(sock::ClickHouseSock, x::AbstractVector{T}) where T <: Number =
    write(sock.io, x)

chwrite(sock::ClickHouseSock, x::AbstractVector{String}) =
    foreach(x -> chwrite(sock, x), x)

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct SecurityGroupViewResult <: SwaggerModel
    networkInterfaces::Any # spec type: Union{ Nothing, Vector{SecurityGroupNetworkInterface} } # spec name: networkInterfaces

    function SecurityGroupViewResult(;networkInterfaces=nothing)
        o = new()
        validate_property(SecurityGroupViewResult, Symbol("networkInterfaces"), networkInterfaces)
        setfield!(o, Symbol("networkInterfaces"), networkInterfaces)
        o
    end
end # type SecurityGroupViewResult

const _property_map_SecurityGroupViewResult = Dict{Symbol,Symbol}(Symbol("networkInterfaces")=>Symbol("networkInterfaces"))
const _property_types_SecurityGroupViewResult = Dict{Symbol,String}(Symbol("networkInterfaces")=>"Vector{SecurityGroupNetworkInterface}")
Base.propertynames(::Type{ SecurityGroupViewResult }) = collect(keys(_property_map_SecurityGroupViewResult))
Swagger.property_type(::Type{ SecurityGroupViewResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_SecurityGroupViewResult[name]))}
Swagger.field_name(::Type{ SecurityGroupViewResult }, property_name::Symbol) =  _property_map_SecurityGroupViewResult[property_name]

function check_required(o::SecurityGroupViewResult)
    true
end

function validate_property(::Type{ SecurityGroupViewResult }, name::Symbol, val)
end


text:
@doc raw"""
    BayesianSVM() -> BernoulliLikelihood

The [Bayesian SVM](https://arxiv.org/abs/1707.05532) is a Bayesian interpretation of the classical SVM.
```math
p(y|f) \propto \exp(2 \max(1-yf, 0))
```

---

For the analytic version of the likelihood, it is augmented via:

```math
p(y|f, ω) = \frac{1}{\sqrt(2\pi\omega)} \exp\left(-\frac{(1+\omega-yf)^2}{2\omega})\right)
```

where ``ω \sim 1[0,\infty)`` has an improper prior (his posterior is however has a valid distribution, a Generalized Inverse Gaussian). For reference [see this paper](http://ecmlpkdd2017.ijs.si/papers/paperID502.pdf).
"""
BayesianSVM() = BernoulliLikelihood(SVMLink())

struct SVMLink <: AbstractLink end

implemented(::BernoulliLikelihood{<:SVMLink}, ::AnalyticVI) = true

function (::SVMLink)(f::Real)
    return svmlikelihood(f)
end

# Return likelihood equivalent to SVM hinge loss
function svmlikelihood(f::Real)
    pos = svmpseudolikelihood(f)
    return pos ./ (pos .+ svmpseudolikelihood(-f))
end

# Return the pseudo likelihood of the SVM hinge loss
function svmpseudolikelihood(f::Real)
    return exp(-2 * max(1 - f, zero(f)))
end

function local_updates!(
    local_vars,
    ::BernoulliLikelihood{<:SVMLink},
    y::AbstractVector,
    μ::AbstractVector{T},
    diagΣ::AbstractVector,
) where {T}
    map!(local_vars.c, μ, diagΣ, y) do μ, σ², y
        abs2(one(T) - y * μ) + σ²
    end
    map!(inv ∘ sqrt, local_vars.θ, local_vars.c)
    return local_vars
end

@inline function ∇E_μ(
    ::BernoulliLikelihood{<:SVMLink}, ::AOptimizer, y::AbstractVector, state
)
    return (y .* (state.θ .+ one(eltype(state.θ))),)
end

@inline function ∇E_Σ(
    ::BernoulliLikelihood{<:SVMLink}, ::AOptimizer, ::AbstractVector, state
)
    return (state.θ / 2,)
end

## ELBO

function expec_loglikelihood(
    ::BernoulliLikelihood{<:SVMLink},
    ::AnalyticVI,
    y,
    μ::AbstractVector,
    diag_cov::AbstractVector,
    state,
)
    tot = -length(y) * logtwo / 2
    tot += dot(μ, y)
    tot += -dot(state.θ, diag_cov) / 2 + dot(state.θ, abs2.(one(eltype(μ)) .- y .* μ))
    return tot
end

function AugmentedKL(l::BernoulliLikelihood{<:SVMLink}, state, ::Any)
    Zygote.@ignore(GIGEntropy(l, state))
end

function GIGEntropy(::BernoulliLikelihood{<:SVMLink}, state)
    return sum(log.(state.c)) / 2 + sum(log.(2.0 * besselk.(0.5, sqrt.(state.c)))) -
           sum(sqrt.(state.c)) / 2
end


text:
using PooksoftOptionsKit
using Documenter

# call to make docs -
makedocs(sitename="PooksoftOptionsKit.jl",

    pages = [
        "index.md",
        "Option pricing models" => [
            "binomial.md"
        ],
        "Profit and loss simulations" => [
            "profit_and_loss.md"
        ],
        "The Greeks" => [
            "greeks.md"
        ]
    ]
)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""Contains example usage of the features in the OpticSim.jl package."""
module Examples
using ..OpticSim
using ..OpticSim.Vis
using ..OpticSim.Geometry
using ..OpticSim.Emitters
using ..OpticSim.GlassCat

using StaticArrays
using DataFrames
using Images
using Unitful
using Plots
using LinearAlgebra

include("docs_examples.jl")
include("other_examples.jl")

end #module Examples
export Examples


text:
using QuOptimalControl
using Test
using LinearAlgebra
using StaticArrays
using Optim

include("setup_tests.jl")
include("state_transfer_tests.jl")
include("unitary_gate_tests.jl")


text:
#=
Type definition and methods containing the overarching backtesting object fueling the engine
=#

mutable struct Strategy
    universe::Universe
    indicator::Indicator
    rules::Tuple{Vararg{Rule}}
    portfolio::Portfolio
    results::Results
    function Strategy(universe::Universe,
                      indicator::Indicator,
                      rules::Tuple{Vararg{Rule}},
                      portfolio::Portfolio=Portfolio(universe))
        return new(universe, indicator, rules, portfolio, Results())
    end
end

function generate_trades(strat::Strategy; verbose::Bool=true)::Dict{String,TS}
    all_trades = Dict{String,TS}()
    for asset in strat.universe.assets
        verbose ? print("Generating trades for asset $asset...") : nothing
        trades = TS(falses(size(strat.universe.data[asset],1), length(strat.rules)),
                    strat.universe.data[asset].index)
        local indicator_data = calculate(strat.indicator, strat.universe.data[asset])
        for (i,rule) in enumerate(strat.rules);
            trades[:,i] = rule.trigger.fun(indicator_data)
        end
        all_trades[asset] = trades
        verbose ? print("Done.\n") : nothing
    end
    return all_trades
end

function generate_trades!(strat::Strategy; args...)::Void
    strat.results.trades = generate_trades(strat; args...)
    return nothing
end

function backtest(strat::Strategy; px_trade::Symbol=:Open, px_close::Symbol=:Settle, verbose::Bool=true)::Dict{String,TS{Float64}}
    if isempty(strat.results.trades)
        generate_trades!(strat, verbose=verbose)
    end
    result = Dict{String,TS}()
    for asset in strat.universe.assets
        verbose ? print("Running backtest for asset $asset...") : nothing
        trades = strat.results.trades[asset].values
        N = size(trades, 1)
        summary_ts = strat.universe.data[asset]
        #TODO: add setindex! method for TS objects using Symbol and Vector to assign inplace
        #TODO: generalize this logic to incorporate order types
        #FIXME: generalize this logic to use the actual rules (this is a temporary quickfix)
        trade_price = summary_ts[px_trade].values
        close_price = summary_ts[px_close].values
        pos = zeros(Float64, N)
        pnl = zeros(Float64, N)
        do_trade = false
        for t in 2:N
            for (i,rule) in enumerate(strat.rules)
                if trades[t-1,i] != 0
                    do_trade = true
                    #TODO: fill out this logic with the various order types
                    order_side = rule.action in (long,buy) ? 1 : rule.action in (short,sell) ? -1 : 0
                    #TODO: add logic here for the int vs. float argument type to order function
                    (order_qty,) = rule.args
                    #if isa(order_qty, Int); else FIXME: portfolio adjustment logic; end
                    pos[t] = order_qty * order_side
                    pnl[t] = pos[t] * (close_price[t] - trade_price[t])
                end
            end
            if !do_trade
                pos[t] = pos[t-1]
                pnl[t] = pos[t] * (close_price[t]-close_price[t-1])
            end
            do_trade = false
        end
        summary_ts = [summary_ts TS([pos pnl cumsum(pnl)], summary_ts.index, [:Pos,:PNL,:CumPNL])]
        result[asset] = summary_ts
        verbose ? print("Done.\n") : nothing
    end
    return result
end

function backtest!(strat::Strategy; args...)::Void
    strat.results.backtest = backtest(strat; args...)
    return nothing
end

Base.copy(strat::Strategy) = Strategy(strat.universe, strat.indicator, strat.rules)

# define matrix row iterator protocol
# this allows us to `enumerate(EachRow(M))`
# thereby getting the count of the iteration as well as the row
immutable EachRow{T<:AbstractMatrix}
    A::T
end
Base.start(::EachRow) = 1
Base.next(itr::EachRow, s) = (itr.A[s,:], s+1)
Base.done(itr::EachRow, s) = s > size(itr.A,1)

#TODO: more meaningful progres information
#TODO: parallel processing
#TODO: streamline this so that it doesnt run so slow (seems to be recompiling at each run)
function optimize(strat::Strategy; samples::Int=0, seed::Int=0, verbose::Bool=true, summary_fun::Function=cum_pnl, args...)::Matrix
    strat_save = copy(strat)
    n_runs = get_n_runs(strat.indicator.paramset)
    idx_samples::Vector{Int} = collect(1:n_runs)
    if samples > 0
        srand(seed)
        idx_samples = rand(idx_samples, samples)
    else
        samples = n_runs
    end
    combos = get_param_combos(strat.indicator.paramset, n_runs=n_runs)[idx_samples,:]
    result = zeros(samples)
    for (run, combo) in enumerate(idx_samples)
        verbose ? println("Run $run/$samples") : nothing
        strat.indicator.paramset.arg_defaults = combo
        generate_trades!(strat, verbose=false)
        backtest!(strat, verbose=false; args...)
        result[run] = summary_fun(strat.results.backtest)
    end
    # prevent out-of-scope alteration of strat object
    strat = strat_save
    return result
end

# TODO: implement function to edit results member of strat in place
function optimize!(strat::Strategy; samples::Int=0, seed::Int=0, verbose::Bool=true, summary_fun::Function=cum_pnl, args...)::Void
    n_runs = get_n_runs(strat.indicator.paramset)
    idx_samples::Vector{Int} = collect(1:n_runs)
    if samples > 0
        if seed >= 0
            srand(seed)
        end
        idx_samples = rand(idx_samples, samples)
    else
        samples = n_runs
    end
    combos = get_param_combos(strat.indicator.paramset, n_runs=n_runs)[idx_samples,:]
    strat.results.optimization = zeros(samples,1)
    for (run, combo) in enumerate(EachRow(combos))
        verbose ? println("Run $run/$samples") : nothing
        strat.indicator.paramset.arg_defaults = combo
        generate_trades!(strat, verbose=false)
        backtest!(strat, verbose=false; args...)
        strat.results.optimization[run] = summary_fun(strat.results)
    end
    strat.results.optimization = [combos strat.results.optimization]
    return nothing
end


text:
module MathieuFunctions
using LinearAlgebra

  include("char-values.jl")

end


text:
################################################################################
#                             BASIC REFERENCE METHODS
################################################################################
"""
    JuMP.owner_model(cref::InfOptConstraintRef)::InfiniteModel

Extend `JuMP.owner_model` to return the infinite model associated with `cref`.

**Example**
```julia-repl
julia> model = owner_model(cref)
An InfiniteOpt Model
Minimization problem with:
Finite Parameters: 0
Infinite Parameters: 3
Variables: 3
Derivatives: 0
Measures: 0
Objective function type: GeneralVariableRef
`GenericAffExpr{Float64,GeneralVariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 1 constraint
Names registered in the model: g, t, h, x
Optimizer model backend information:
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
```
"""
JuMP.owner_model(cref::InfOptConstraintRef)::InfiniteModel = cref.model

"""
    JuMP.index(cref::InfOptConstraintRef)::InfOptConstraintIndex

Extend `JuMP.index` to return the index of an `InfiniteOpt` constraint `cref`.

**Example**
```julia-repl
julia> index(cref)
InfOptConstraintIndex(2)
```
"""
JuMP.index(cref::InfOptConstraintRef)::InfOptConstraintIndex = cref.index

# Extend Base and JuMP functions
function Base.:(==)(v::InfOptConstraintRef, w::InfOptConstraintRef)::Bool
    return v.model === w.model && v.index == w.index
end
Base.broadcastable(cref::InfOptConstraintRef) = Ref(cref)

################################################################################
#                             CORE OBJECT METHODS
################################################################################
# Extend _add_data_object
function _add_data_object(
    model::InfiniteModel,
    object::ConstraintData
    )::InfOptConstraintIndex
    return MOIUC.add_item(model.constraints, object)
end

# Extend _data_dictionary
function _data_dictionary(
    cref::InfOptConstraintRef
    )::MOIUC.CleverDict{InfOptConstraintIndex, <:ConstraintData}
    return JuMP.owner_model(cref).constraints
end

# Extend _data_object
function _data_object(cref::InfOptConstraintRef)::ConstraintData
    object = get(_data_dictionary(cref), JuMP.index(cref), nothing)
    if object === nothing
        error("Invalid constraint reference, cannot find corresponding ",
              "constraint in the model. This is likely caused by using the ",
              "reference of a deleted constraint.")
    end
    return object
end

# Return the core constraint object
function _core_constraint_object(
    cref::InfOptConstraintRef
    )::JuMP.AbstractConstraint
    return _data_object(cref).constraint
end

## Set helper methods for adapting data_objects with parametric changes
# No change needed
function _adaptive_data_update(
    cref::InfOptConstraintRef,
    c::C,
    data::ConstraintData{C}
    )::Nothing where {C <: JuMP.AbstractConstraint}
    data.constraint = c
    return
end

# Reconstruction is necessary
function _adaptive_data_update(
    cref::InfOptConstraintRef,
    c::C1,
    data::ConstraintData{C2}
    )::Nothing  where {C1, C2}
    new_data = ConstraintData(c, data.object_nums, data.name,
                              data.measure_indices, data.is_info_constraint)
    _data_dictionary(cref)[JuMP.index(cref)] = new_data
    return
end

# Update the core constraint object
function _set_core_constraint_object(
    cref::InfOptConstraintRef,
    constr::JuMP.AbstractConstraint
    )::Nothing
    _adaptive_data_update(cref, constr, _data_object(cref))
    set_optimizer_model_ready(JuMP.owner_model(cref), false)
    return
end

# Extend _object_numbers
function _object_numbers(cref::InfOptConstraintRef)::Vector{Int}
    return _data_object(cref).object_nums
end

# Extend _measure_dependencies
function _measure_dependencies(
    cref::InfOptConstraintRef
    )::Vector{MeasureIndex}
    return _data_object(cref).measure_indices
end

# Return if this constraint is an info constraint
function _is_info_constraint(cref::InfOptConstraintRef)::Bool
    return _data_object(cref).is_info_constraint
end

# Extend _delete_data_object
function _delete_data_object(cref::InfOptConstraintRef)::Nothing
    delete!(_data_dictionary(cref), JuMP.index(cref))
    return
end

################################################################################
#                             DEFINITION METHODS
################################################################################
# Check that domain_restrictions argument is valid
function _check_restrictions(
    restrictions::DomainRestrictions{GeneralVariableRef};
    _error = error
    )::Nothing
    depend_counter = Dict{DependentParameterRef, Int}()
    for (pref, domain) in restrictions
        # check that pref is an infinite parameter
        if !(_index_type(pref) <: InfiniteParameterIndex)
            _error("Can only specify infinite parameters for domain restrictions.")
        end
        # check that respects lower bound
        if JuMP.has_lower_bound(pref) && (JuMP.lower_bound(domain) < JuMP.lower_bound(pref))
                _error("Specified parameter lower bound exceeds that defined " *
                       "for $pref.")
        end
        # check that respects upper bound
        if JuMP.has_upper_bound(pref) && (JuMP.upper_bound(domain) > JuMP.upper_bound(pref))
                _error("Specified parameter upper bound exceeds that defined " *
                       "for $pref.")
        end
        # keep track of dependent parameters with equality restrictions to ensure completeness
        if (_index_type(pref) == DependentParameterIndex) &&
           (JuMP.lower_bound(domain) == JuMP.upper_bound(domain))
            idx = DependentParameterIndex(DependentParametersIndex(_raw_index(pref)), 1)
            dumby_pref = dispatch_variable_ref(pref.model, idx)
            if haskey(depend_counter, dumby_pref)
                depend_counter[dumby_pref] += 1
            else
                depend_counter[dumby_pref] = 1
            end
        end
    end
    # check dimensions of dependent parameter equalities (if any are provided)
    for (pref, amt) in depend_counter
        if _num_parameters(pref) != amt
            _error("Cannot specify equality domain restrictions for a subset of ",
                   "dependent infinite parameters.")
        end
    end
    return
end

"""
    JuMP.build_constraint(_error::Function, func, set,
                          restrictions::DomainRestrictions{GeneralVariableRef}
                          )::DomainRestrictedConstraint

Extend `JuMP.buid_constraint` to handle including `restrictions` to its inherit
infinite parameter domains in addition to the traditional `func` in `set` setup.
This returns a `DomainRestrictedConstraint` that can then
be added via `JuMP.add_constraint`. Errors if the restrictions are incompadible
with infinite parameter domains.

**Example**
```julia-repl
julia> restrictions = DomainRestrictions(t => 0)
Subdomain restrictions (1): t = 0

julia> con = build_constraint(error, y + 2, MOI.LessThan(0.0), restrictions);
```
"""
function JuMP.build_constraint(
    _error::Function,
    func,
    set,
    restrictions::DomainRestrictions
    )::DomainRestrictedConstraint
    # make the constraint and check the domain restrictions
    constr = JuMP.build_constraint(_error, func, set)
    _check_restrictions(restrictions, _error = _error)
    return DomainRestrictedConstraint(constr, restrictions)
end

# Validate domain restrictions and add support(s) if needed
function _validate_restrictions(
    model::InfiniteModel,
    restrictions::DomainRestrictions
    )::Nothing
    depend_supps = Dict{DependentParametersIndex, Matrix{Float64}}()
    for (pref, domain) in restrictions
        # check validity
        JuMP.check_belongs_to_model(pref, model)
        # ensure has a support if a point constraint was given
        if (JuMP.lower_bound(domain) == JuMP.upper_bound(domain))
            if _index_type(pref) == IndependentParameterIndex
                # label will be UserDefined
                add_supports(pref, JuMP.lower_bound(domain), check = false)
            else
                idx = DependentParametersIndex(_raw_index(pref))
                if !haskey(depend_supps, idx)
                    dumby_pref = dispatch_variable_ref(model, DependentParameterIndex(idx, 1))
                    depend_supps[idx] = Matrix{Float64}(undef, _num_parameters(dumby_pref), 1)
                end
                depend_supps[idx][_param_index(pref)] = JuMP.lower_bound(domain)
            end
        end
    end
    # add dependent supports if any are given
    for (idx, supp) in depend_supps
        prefs = [dispatch_variable_ref(model, DependentParameterIndex(idx, i))
                 for i in 1:length(supp)]
        add_supports(prefs, supp, check = false) # label will be UserDefined
    end
    return
end

# Used to update the model.var_to_constrs field
function _update_var_constr_mapping(
    vrefs::Vector{GeneralVariableRef},
    cref::InfOptConstraintRef
    )::Nothing
    for vref in vrefs
        dvref = dispatch_variable_ref(vref)
        push!(_constraint_dependencies(dvref), JuMP.index(cref))
        if dvref isa MeasureRef
            push!(_measure_dependencies(cref), JuMP.index(dvref))
        end
        # TODO maybe make mapping for derivatives...
    end
    return
end

"""
    JuMP.add_constraint(model::InfiniteModel, c::JuMP.AbstractConstraint,
                        [name::String = ""])::InfOptConstraintRef

Extend `JuMP.add_constraint` to add a constraint `c` to an infinite model
`model` with name `name`. Returns an appropriate constraint reference whose type
depends on what variables are used to define the constraint. Errors if any
variables do not belong to `model`. This is primarily used as an internal
method for the constraint macros.

**Example**
```julia-repl
julia> @infinite_parameter(model, t in [0, 10]);

julia> @variable(model, g, Infinite(t));

julia> @variable(model, x);

julia> constr = build_constraint(error, g + x, MOI.EqualTo(42));

julia> cref = add_constraint(model, constr, "name")
name : g(t) + x = 42.0, ∀ t ∈ [0, 10]
```
"""
function JuMP.add_constraint(
    model::InfiniteModel,
    c::JuMP.AbstractConstraint,
    name::String = "";
    is_info_constr::Bool = false
    )::InfOptConstraintRef
    # gather the unique list of variable references for testing and mapping
    vrefs = _all_function_variables(JuMP.jump_function(c))
    # test in the model
    for vref in vrefs
        JuMP.check_belongs_to_model(vref, model)
    end
    # get the parameter object numbers
    object_nums = sort!(_object_numbers(vrefs))
    # add the constaint to the model
    constr_object = ConstraintData(c, object_nums, name, MeasureIndex[],
                                   is_info_constr)
    cindex = _add_data_object(model, constr_object)
    cref = InfOptConstraintRef(model, cindex)
    # update the variable mappings and model status
    _update_var_constr_mapping(vrefs, cref)
    set_optimizer_model_ready(model, false)
    # clear out the name dictionary
    model.name_to_constr = nothing
    # return the constraint reference
    return cref
end

# Domain restricted constraints
function JuMP.add_constraint(
    model::InfiniteModel,
    c::DomainRestrictedConstraint,
    name::String = ""
    )::InfOptConstraintRef
    # test domain restrictions and add needed supports
    _validate_restrictions(model, c.restrictions)
    # add the underlying constraint
    cref = JuMP.add_constraint(model, c.constraint, name)
    # add the domain restrictions
    model.constraint_restrictions[JuMP.index(cref)] = c.restrictions
    # return the constraint reference
    return cref
end

################################################################################
#                                 JuMP METHODS
################################################################################
"""
    JuMP.is_valid(model::InfiniteModel, cref::InfOptConstraintRef)::Bool

Extend `JuMP.is_valid` to return `Bool` whether an `InfiniteOpt` constraint
reference is valid.

**Example**
```julia-repl
julia> is_valid(model, cref)
true
```
"""
function JuMP.is_valid(
    model::InfiniteModel,
    cref::InfOptConstraintRef
    )::Bool
    return (model === JuMP.owner_model(cref) &&
            JuMP.index(cref) in keys(_data_dictionary(cref)))
end

"""
    JuMP.constraint_object(cref::InfOptConstraintRef)::JuMP.AbstractConstraint

Extend `JuMP.constraint_object` to return the constraint object associated with
`cref`.

**Example**
```julia-repl
julia> @infinite_parameter(model, t in [0, 10]);

julia> @variable(model, x <= 1);

julia> cref = UpperBoundRef(x);

julia> obj = constraint_object(cref)
ScalarConstraint{GeneralVariableRef,MathOptInterface.LessThan{Float64}}(x,
MathOptInterface.LessThan{Float64}(1.0))
```
"""
function JuMP.constraint_object(
    cref::InfOptConstraintRef
    )::JuMP.AbstractConstraint
    return _core_constraint_object(cref)
end

"""
    JuMP.name(cref::InfOptConstraintRef)::String

Extend `JuMP.name` to return the name of an `InfiniteOpt` constraint.

**Example**
```julia-repl
julia> name(cref)
"constr_name"
```
"""
function JuMP.name(cref::InfOptConstraintRef)::String
    object = get(_data_dictionary(cref), JuMP.index(cref), nothing)
    return object === nothing ? "" : object.name
end

"""
    JuMP.set_name(cref::InfOptConstraintRef, name::String)::Nothing

Extend `JuMP.set_name` to specify the name of a constraint `cref`.

**Example**
```julia-repl
julia> set_name(cref, "new_name")

julia> name(cref)
"new_name"
```
"""
function JuMP.set_name(
    cref::InfOptConstraintRef,
    name::String
    )::Nothing
    _data_object(cref).name = name
    JuMP.owner_model(cref).name_to_constr = nothing
    return
end

# Return a constraint set with an updated value
function _set_set_value(
    set::S,
    value::Real
    ) where {T, S <: Union{MOI.LessThan{T}, MOI.GreaterThan{T}, MOI.EqualTo{T}}}
    return S(convert(T, value))
end

# Enforce that the MOI set is a traditional scalar one
function _enforce_rhs_set(
    set::Union{MOI.LessThan{T}, MOI.GreaterThan{T}, MOI.EqualTo{T}}
    )::Nothing where {T}
    return
end
function _enforce_rhs_set(set)
    error("The right hand side value is not well defined for constraints ",
          "with sets of type `$(typeof(set))`.")
end

"""
    JuMP.set_normalized_rhs(cref::InfOptConstraintRef, value::Real)::Nothing

Set the right-hand side term of `constraint` to `value`.
Note that prior to this step, JuMP will aggregate all constant terms onto the
right-hand side of the constraint. For example, given a constraint `2x + 1 <=
2`, `set_normalized_rhs(con, 4)` will create the constraint `2x <= 4`, not `2x +
1 <= 4`.

```julia-repl
julia> @constraint(model, con, 2x + 1 <= 2)
con : 2 x ≤ 1.0

julia> set_normalized_rhs(con, 4)

julia> con
con : 2 x ≤ 4.0
```
"""
function JuMP.set_normalized_rhs(
    cref::InfOptConstraintRef,
    value::Real
    )::Nothing
    old_constr = JuMP.constraint_object(cref)
    _enforce_rhs_set(JuMP.moi_set(old_constr))
    new_set = _set_set_value(JuMP.moi_set(old_constr), value)
    new_constr = JuMP.ScalarConstraint(JuMP.jump_function(old_constr), new_set)
    _set_core_constraint_object(cref, new_constr)
    return
end

"""
    JuMP.normalized_rhs(cref::InfOptConstraintRef)::Float64

Return the right-hand side term of `cref` after JuMP has converted the
constraint into its normalized form.
"""
function JuMP.normalized_rhs(cref::InfOptConstraintRef)::Float64
    constr = JuMP.constraint_object(cref)
    set = JuMP.moi_set(constr)
    _enforce_rhs_set(set)
    return MOI.constant(set)
end

"""
    JuMP.add_to_function_constant(cref::InfOptConstraintRef, value::Real)::Nothing

Add `value` to the function constant term.
Note that for scalar constraints, JuMP will aggregate all constant terms onto the
right-hand side of the constraint so instead of modifying the function, the set
will be translated by `-value`. For example, given a constraint `2x <=
3`, `add_to_function_constant(c, 4)` will modify it to `2x <= -1`.
```
"""
function JuMP.add_to_function_constant(
    cref::InfOptConstraintRef,
    value::Real
    )::Nothing
    current_value = JuMP.normalized_rhs(cref)
    JuMP.set_normalized_rhs(cref, current_value - value)
    return
end

"""
    JuMP.set_normalized_coefficient(cref::InfOptConstraintRef,
                                    variable::GeneralVariableRef,
                                    value::Real)::Nothing

Set the coefficient of `variable` in the constraint `constraint` to `value`.
Note that prior to this step, JuMP will aggregate multiple terms containing the
same variable. For example, given a constraint `2x + 3x <= 2`,
`set_normalized_coefficient(con, x, 4)` will create the constraint `4x <= 2`.

```julia-repl
julia> con
con : 5 x ≤ 2.0

julia> set_normalized_coefficient(con, x, 4)

julia> con
con : 4 x ≤ 2.0
```
"""
function JuMP.set_normalized_coefficient(
    cref::InfOptConstraintRef,
    variable::GeneralVariableRef,
    value::Real
    )::Nothing
    # update the constraint expression and update the constraint
    old_constr = JuMP.constraint_object(cref)
    new_func = _set_variable_coefficient!(JuMP.jump_function(old_constr),
                                          variable, value) # checks valid
    new_constr = JuMP.ScalarConstraint(new_func, JuMP.moi_set(old_constr))
    _set_core_constraint_object(cref, new_constr)
    return
end

"""
    JuMP.normalized_coefficient(cref::InfOptConstraintRef,
                                variable::GeneralVariableRef)::Float64

Return the coefficient associated with `variable` in `constraint` after JuMP has
normalized the constraint into its standard form.
"""
function JuMP.normalized_coefficient(
    cref::InfOptConstraintRef,
    variable::GeneralVariableRef
    )::Float64
    constr = JuMP.constraint_object(cref)
    func = JuMP.jump_function(constr)
    return _affine_coefficient(func, variable) # checks valid
end

# Return the appropriate constraint reference given the index and model
function _make_constraint_ref(model::InfiniteModel,
    index::InfOptConstraintIndex
    )::InfOptConstraintRef
    return InfOptConstraintRef(model, index)
end

"""
    JuMP.constraint_by_name(model::InfiniteModel,
                            name::String)::Union{InfOptConstraintRef, Nothing}

Extend `JuMP.constraint_by_name` to return the constraint reference
associated with `name` if one exists or returns nothing. Errors if more than
one constraint uses the same name.

**Example**
```julia-repl
julia> constraint_by_name(model, "constr_name")
constr_name : x + pt = 3.0
```
"""
function JuMP.constraint_by_name(
    model::InfiniteModel,
    name::String
    )::Union{InfOptConstraintRef, Nothing}
    if model.name_to_constr === nothing
        # Inspired from MOI/src/Utilities/model.jl
        model.name_to_constr = Dict{String, Int}()
        for (index, data_object) in model.constraints
            constr_name = data_object.name
            if haskey(model.name_to_constr, constr_name)
                model.name_to_constr[constr_name] = InfOptConstraintIndex(-1)
            else
                model.name_to_constr[constr_name] = index
            end
        end
    end
    index = get(model.name_to_constr, name, nothing)
    if index === nothing
        return nothing
    elseif index == InfOptConstraintIndex(-1)
        error("Multiple constraints have the name $name.")
    else
        return _make_constraint_ref(model, index)
    end
end

"""
    JuMP.num_constraints(model::InfiniteModel, [function_type], [set_type])::Int

Extend `JuMP.num_constraints` to return the number of constraints with a
partiuclar function type and set type.

**Example**
```julia-repl
julia> num_constraints(model, FiniteVariableRef, MOI.LessThan)
1

julia> num_constraints(model, FiniteVariableRef)
3

julia> num_constraints(model, MOI.LessThan)
2

julia> num_constraints(model)
4
```
"""
function JuMP.num_constraints(
    model::InfiniteModel,
    function_type,
    set_type
    )::Int
    counter = 0
    for (index, data_object) in model.constraints
        if isa(JuMP.jump_function(data_object.constraint), function_type) &&
           isa(JuMP.moi_set(data_object.constraint), set_type)
            counter += 1
        end
    end
    return counter
end

# Function type only
function JuMP.num_constraints(
    model::InfiniteModel,
    function_type
    )::Int
    return JuMP.num_constraints(model, function_type, MOI.AbstractSet)
end

# Set type only
function JuMP.num_constraints(
    model::InfiniteModel,
    set_type::Type{<:MOI.AbstractSet}
    )::Int
    return JuMP.num_constraints(model, Any, set_type)
end

# All the constraints
function JuMP.num_constraints(model::InfiniteModel)::Int
    return length(model.constraints)
end

"""
    JuMP.all_constraints(model::InfiniteModel, [function_type], [set_type])::Vector{InfOptConstraintRef}

Extend `JuMP.all_constraints` to return a list of all the constraints with a
particular function type and set type.

**Example**
```julia-repl
julia> all_constraints(model, GeneralVariableRef, MOI.LessThan)
1-element Array{InfOptConstraintRef,1}:
 x ≤ 1.0

julia> all_constraints(model, GeneralVariableRef)
3-element Array{InfOptConstraintRef,1}:
 x ≥ 0.0
 x ≤ 3.0
 x integer

julia> all_constraints(model, MOI.GreaterThan)
3-element Array{InfOptConstraintRef,1}:
 x ≥ 0.0
 g(t) ≥ 0.0, ∀ t ∈ [0, 6]
 g(0.5) ≥ 0.0

julia> all_constraints(model)
5-element Array{InfOptConstraintRef,1}:
 x ≥ 0.0
 x ≤ 3.0
 x integer
 g(t) ≥ 0.0, ∀ t ∈ [0, 6]
 g(0.5) ≥ 0.0
```
"""
function JuMP.all_constraints(
    model::InfiniteModel,
    function_type,
    set_type
    )::Vector{InfOptConstraintRef}
    constr_list = Vector{InfOptConstraintRef}(undef,
                           JuMP.num_constraints(model, function_type, set_type))
    counter = 1
    for (index, object) in model.constraints
        if isa(JuMP.jump_function(object.constraint), function_type) &&
           isa(JuMP.moi_set(object.constraint), set_type)
            constr_list[counter] = _make_constraint_ref(model, index)
            counter += 1
        end
    end
    return constr_list
end

# Function type only
function JuMP.all_constraints(
    model::InfiniteModel,
    function_type
    )::Vector{InfOptConstraintRef}
    return JuMP.all_constraints(model, function_type, MOI.AbstractSet)
end

# Set type only
function JuMP.all_constraints(
    model::InfiniteModel,
    set_type::Type{<:MOI.AbstractSet}
    )::Vector{InfOptConstraintRef}
    return JuMP.all_constraints(model, JuMP.AbstractJuMPScalar, set_type)
end

# All the constraints
function JuMP.all_constraints(
    model::InfiniteModel
    )::Vector{InfOptConstraintRef}
    return [_make_constraint_ref(model, idx) for (idx, _) in model.constraints]
end

"""
    JuMP.list_of_constraint_types(model::InfiniteModel)::Vector{Tuple{DataType, DataType}}

Extend `JuMP.list_of_constraint_types` to return a list of tuples that contain
all the used combinations of function types and set types in the model.

**Example**
```julia-repl
julia> all_constraints(model)
3-element Array{Tuple{DataType,DataType},1}:
 (GeneralVariableRef, MathOptInterface.LessThan{Float64})
 (GeneralVariableRef, MathOptInterface.GreaterThan{Float64})
 (GeneralVariableRef, MathOptInterface.Integer)
```
"""
function JuMP.list_of_constraint_types(
    model::InfiniteModel
    )::Vector{Tuple{DataType, DataType}}
    type_set = Set{Tuple{DataType, DataType}}()
    for (index, object) in model.constraints
        push!(type_set, (typeof(JuMP.jump_function(object.constraint)),
                         typeof(JuMP.moi_set(object.constraint))))
    end
    return collect(type_set)
end

################################################################################
#                         PARAMETER REFERENCE METHODS
################################################################################
"""
    parameter_refs(cref::InfOptConstraintRef)::Tuple

Return the tuple of infinite parameter references that determine the infinite
dependencies of `cref`.

**Example**
```julia-repl
julia> parameter_refs(cref)
(t,)
```
"""
function parameter_refs(cref::InfOptConstraintRef)::Tuple
    model = JuMP.owner_model(cref)
    obj_indices = _param_object_indices(model)[_object_numbers(cref)]
    return Tuple(_make_param_tuple_element(model, idx) for idx in obj_indices)
end

################################################################################
#                         DOMAIN RESTRICTION METHODS
################################################################################
"""
    has_domain_restrictions(cref::InfOptConstraintRef)::Bool

Return a `Bool` indicating if `cref` is limited to a sub-domain as defined
by a [`DomainRestrictions`](@ref) object.

**Example**
```julia-repl
julia> has_domain_restrictions(cref)
true
```
"""
function has_domain_restrictions(cref::InfOptConstraintRef)::Bool
    return !isempty(domain_restrictions(cref))
end

"""
    domain_restrictions(cref::InfOptConstraintRef)::DomainRestrictions{GeneralVariableRef}

Return the [`DomainRestrictions`](@ref) object associated with the constraint
`cref`.

**Example**
```julia-repl
julia> domain_restrictions(cref)
Subdomain restrictions (1): t ∈ [0, 2]
```
"""
function domain_restrictions(
    cref::InfOptConstraintRef
    )::DomainRestrictions{GeneralVariableRef}
    return get(JuMP.owner_model(cref).constraint_restrictions, JuMP.index(cref),
               DomainRestrictions())
end

"""
    set_domain_restrictions(cref::InfOptConstraintRef,
                         restrictions:DomainRestrictions{GeneralVariableRef};
                         [force::Bool = false])::Nothing

Specify a new [`DomainRestrictions`](@ref) object `restrictions` for the
constraint `cref`. Errors if `cref` already has restrictions and `force = false`.
Where possible it is recommended to use [`add_domain_restrictions`](@ref) instead.

**Example**
```julia-repl
julia> set_domain_restrictions(cref, DomainRestrictions(t => [0, 2]))

julia> domain_restrictions(cref)
Subdomain restrictions (1): t ∈ [0, 2]
```
"""
function set_domain_restrictions(
    cref::InfOptConstraintRef,
    restrictions::DomainRestrictions{GeneralVariableRef};
    force::Bool = false
    )::Nothing
    if has_domain_restrictions(cref) && !force
        error("$cref already has domain restrictions. Consider adding more using " *
               "`add_domain_restrictions` or overwriting them by setting " *
               "the keyword argument `force = true`")
    end
    # check that restrictions are valid and add support(s) if necessary
    _check_restrictions(restrictions)
    model = JuMP.owner_model(cref)
    _validate_restrictions(model, restrictions)
    # set the new restrictions
    model.constraint_restrictions[JuMP.index(cref)] = restrictions
    # update status
    set_optimizer_model_ready(JuMP.owner_model(cref), false)
    return
end

# Update the current restrictions to overlap with the new ones if possible
function _update_restrictions(
    old::DomainRestrictions{GeneralVariableRef},
    new::DomainRestrictions{GeneralVariableRef}
    )::Nothing
    # check each new restriction
    for (pref, domain) in new
        # we have a new restriction
        if !haskey(old, pref)
            old[pref] = domain
        # the previous domain and the new one do not overlap
        elseif (JuMP.lower_bound(domain) > JuMP.upper_bound(old[pref])) ||
               (JuMP.upper_bound(domain) < JuMP.lower_bound(old[pref]))
            error("The new domain restrictions are incompatible with the ",
                  "existing ones.")
        # we have an existing restriction
        else
            # we have a new stricter lower bound to update with
            if JuMP.lower_bound(domain) > JuMP.lower_bound(old[pref])
                old[pref] = IntervalDomain(JuMP.lower_bound(domain),
                                            JuMP.upper_bound(old[pref]))
            end
            # we have a new stricter upper bound to update with
            if JuMP.upper_bound(domain) < JuMP.upper_bound(old[pref])
                old[pref] = IntervalDomain(JuMP.lower_bound(old[pref]),
                                            JuMP.upper_bound(domain))
            end
        end
    end
    return
end

"""
    add_domain_restrictions(cref::InfOptConstraintRef,
                         new_restrictions::DomainRestrictions{GeneralVariableRef}
                         )::Nothing

Add additional domain restrictions to `cref` such that it is defined over the
sub-domain based on `pref` from `lower` to `upper`.

```julia-repl
julia> add_domain_restrictions(cref, DomainRestrictions(t => [0, 2]))

julia> domain_restrictions(cref)
Subdomain restrictions (1): t ∈ [0, 2]
```
"""
function add_domain_restrictions(
    cref::InfOptConstraintRef,
    new_restrictions::DomainRestrictions{GeneralVariableRef}
    )::Nothing
    # check the new restrictions
    _check_restrictions(new_restrictions)
    model = JuMP.owner_model(cref)
    _validate_restrictions(model, new_restrictions)
    # add the restrictions
    if has_domain_restrictions(cref)
        _update_restrictions(domain_restrictions(cref), new_restrictions)
    else
        model.constraint_restrictions[JuMP.index(cref)] = new_restrictions
    end
    # update the optimizer model status
    set_optimizer_model_ready(JuMP.owner_model(cref), false)
    return
end

"""
    delete_domain_restrictions(cref::InfOptConstraintRef)::Nothing

Delete all the domain restrictions of the constraint `cref`. Note any restrictions that
are needed for finite variables inside in `cref` will be unaffected.

**Example**
```julia-repl
julia> @constraint(model, c1, y <= 42, DomainRestrictions(x => 0))
c1 : y(x) ≤ 42, ∀ x[1] = 0, x[2] = 0

julia> delete_domain_restrictions(c1)

julia> c1
c1 : y(x) ≤ 42, ∀ x[1] ∈ [-1, 1], x[2] ∈ [-1, 1]
```
"""
function delete_domain_restrictions(
    cref::InfOptConstraintRef
    )::Nothing
    # delete the restrictions if there are any
    delete!(JuMP.owner_model(cref).constraint_restrictions, JuMP.index(cref))
    # update status
    set_optimizer_model_ready(JuMP.owner_model(cref), false)
    return
end

################################################################################
#                                 DELETION
################################################################################
"""
    JuMP.delete(model::InfiniteModel, cref::InfOptConstraintRef)::Nothing

Extend `JuMP.delete` to delete an `InfiniteOpt` constraint and all associated
information. Errors if `cref` is invalid.

**Example**
```julia-repl
julia> print(model)
Min measure(g(t)*t) + z
Subject to
 z ≥ 0.0
 g(t) + z ≥ 42.0, ∀ t ∈ [0, 6]

julia> delete(model, cref)

julia> print(model)
Min measure(g(t)*t) + z
Subject to
 z ≥ 0.0
```
"""
function JuMP.delete(
    model::InfiniteModel,
    cref::InfOptConstraintRef
    )::Nothing
    # check valid reference
    @assert JuMP.is_valid(model, cref) "Invalid constraint reference."
    # update variable dependencies
    constr = _core_constraint_object(cref)
    all_vrefs = _all_function_variables(JuMP.jump_function(constr))
    for vref in all_vrefs
        filter!(e -> e != JuMP.index(cref), _constraint_dependencies(vref))
    end
    # delete any restrictions there are
    delete_domain_restrictions(cref)
    # delete constraint information
    _delete_data_object(cref)
    # reset optimizer model status
    set_optimizer_model_ready(model, false)
    return
end


text:
####################################################################################################
#
# gas_water_trans.jl - a julia-language script to model the (1) coupled advection of gas and
# water, and (2) diffusive-dispersive component transport in partially water-saturated porous media
#
# model is not currently set up to handle fully water-saturated conditions
#
# model definition can be multi-dimensional, but cross-term on dispersivity tensor is not
# considered (i.e., only longitudinal dispersion is simulated for aqueous-phase transport)
#
# density-driven gas-phase (and water-phase) advection is not considered
#
####################################################################################################

### constants
const g = 9.807             # gravitational acceleration (m/sec^2)
const R = 8.314             # universal gas constant (J K-1 mol-1)
const T = 298.15 			# ambient temperature (25 C)
const rho_w = 1000.         # density of water (kg/m^3)
const u_w = 8.9e-4          # viscosity of water (Pa*sec)
const u_g = 1.84e-5         # gas viscosity

########################
# data types
########################

type Node
	x::Float64 										# element's node point location
	y::Float64
	z::Float64
	vol::Float64 									# volume of element
	k::Float64 										# absolute permeability
	kr_w::Float64 									# relative permeability, water abd gas phases
	kr_g::Float64
	phi::Float64 									# porosity
	n::Float64 										# Van Genuchten relative permeability curve factors
	m::Float64
	alpha::Float64 									# Van Genuchten capillary pressure curve factor
	Sr::Float64 									# residual water saturation
	Q_index::Int64 									# water source index number
	Q::Float64 										# external water volumetric flux
	F_index::Int64 									# gas source index number
	F::Float64 										# external total (molar) gas mass flux
	alpha_L::Float64 								# longitudinal dispersivity
	P::Float64 										# absolute pressure
	Pc::Float64 									# capillary pressure
	S::Float64 										# water saturation (calculated)
	Xg::Array{Float64, 1} 							# mole fractions of volatile components in gas phase
	C_g::Array{Float64, 1} 							# concentration (moles/vol) of gas components in gas phase
	C_aq::Array{Float64, 1} 						# concentration (moles/vol) of components in aqueous phase
	connect_list::Array{Tuple{Int64, Int64}} 		# list of connected elements/nodes (connection number, connecting node index number)
	sigma_w::Float64 								# water and gas total conductance (updated each time step)
	sigma_g::Float64
	sigma_c::Float64 								# aqueous component transport total conductance (updated each time step)
	sigma_x::Array{Float64, 1}	 					# gas component(s) transport total conductance (updated each time step)
	tot_Qw_in::Float64 								# total flow balances for water from volume element by advection, per time step
	tot_Qw_out::Float64
	tot_Qg_in::Float64 								# total flow balances for gas from volume element by advection, per time step
	tot_Qg_out::Float64
end

type Connect
	node_1::Int64 									# index numbers for connected nodes
	node_2::Int64
	dx::Float64										# inter-node distance (for computing gradients)
	area::Float64 									# volume element connection interfacial area
	conduct_w::Float64 								# connection conductances for water and gas phases
	conduct_g::Float64
	conduct_c::Float64 								# connection conductances for aqueous component transport
	conduct_x::Array{Float64, 1}	 				# connection conductances for gas component(s) transport
	Qw::Float64 									# advective water flux
	Qg::Float64 									# advective gas flux
end

type Component
	name::AbstractString 							# name of component
	volatile::Bool 									# if true, component is present in gas phase
	D::Float64 										# gaseous diffusion coefficient
	MW::Float64 									# molecular weight (for quantifying gas-phase source contribution)
end

type Knobs
	gamma::Float64 									# time-stepping weighting factor for implicit solution scheme
	dt_init::Float64 								# initial, minimum, and maximum time step size
	dt_min::Float64
	dt_max::Float64
	dP_max::Float64 								# maximum change in absolute (gas) pressure, per time step
	dPc_max::Float64 								# maximum change in capillary pressure, per time step
	dt_decrease::Float64							# time step reduction and increase factors
	dt_increase::Float64
	upw_w::Float64									# upstream weighting factors for water and gas component transport
	upw_g::Float64
end

#############################
# input functions
#############################

function ReadSource(source_flag::Int64)::Array{Float64, 2}
	# read source file (for gases or aqueous components) and return distribution array
	if source_flag == 1
		infile = "source_gas.txt"
	else
		infile = "source_aqueous.txt"
	end
	data = readdlm(infile, '\t', header=true)
	nrows, ncols = size(data[1])
	source = zeros(Float64, nrows, ncols-1)
	for i = 1:nrows
		for j = 2:ncols
			source[i, j-1] = Float64(data[1][i, j])
		end
	end
	return source
end

function ReadComps()
	# read in component properties, return separate items for gases and aqueous components
	comp = Component[]
	gas = Component[]
	data = readdlm("components.txt", '\t', header=true)
	for i = 1:size(data[1], 1)
		name = data[1][i, 1]
		volatile = Bool(data[1][i, 2])
		D = Float64(data[1][i, 3])
		MW = Float64(data[1][i, 4])
		push!(comp, Component(name, volatile, D, MW))
		if volatile == true
			push!(gas, Component(name, volatile, D, MW))
		end
	end
	num_comps = length(comp)
	num_gas = length(gas)
	println("Read chemical component properties.")
	return comp, gas, num_comps, num_gas
end

function GetKnobs()::Knobs
	# read numerical model "knobs" from file
	data = readdlm("knobs.txt", '\t', header=false)
	gamma = Float64(data[1, 2])
	dt_init = Float64(data[2, 2])
	dt_min = Float64(data[3, 2])
	dt_max = Float64(data[4, 2])
	dP_max = Float64(data[5, 2])
	dPc_max = Float64(data[6, 2])
	dt_decrease = Float64(data[7, 2])
	dt_increase = Float64(data[8, 2])
	upw_w = Float64(data[9, 2])
	upw_g = Float64(data[10, 2])
	knobs = Knobs(gamma, dt_init, dt_min, dt_max, dP_max, dPc_max, dt_decrease, dt_increase, upw_w, upw_g)
	println("Read in computational knobs.")
	return knobs
end

function ReadNodes(num_comps::Int64, num_gas::Int64, S_gMW::Array{Float64,1})

	# read in nodes file and populate node type array
	node = Node[]
	num_nodes = 0
	sigma_x = zeros(Float64, num_gas) 					# initialize total gas component transport conductance array (with respect to each gas phase)
	data = readdlm("nodes.txt", '\t', header=true)
	for i = 1:size(data[1], 1)
		x = Float64(data[1][i, 1])
		y = Float64(data[1][i, 2])
		z = Float64(data[1][i, 3])
		num = Int64(data[1][i, 4])
		x_step = Float64(data[1][i, 5])
		y_step = Float64(data[1][i, 6])
		z_step = Float64(data[1][i, 7])
		vol = Float64(data[1][i, 8])
		k = Float64(data[1][i, 9])
		phi = Float64(data[1][i, 10])
		n = Float64(data[1][i, 11])
		m = 1. - 1/n
		alpha = Float64(data[1][i, 12])
		Sr = Float64(data[1][i, 13])
		Q_index = Int64(data[1][i, 14])
		Q = Float64(data[1][i, 15])
		F_index = Int64(data[1][i, 16])
		F_blended_mass = Float64(data[1][i, 17])
		F = F_blended_mass/S_gMW[F_index]
		alpha_L = Float64(data[1][i, 18])
		P = Float64(data[1][i, 19]) 					# initial condition: pressure
		S = Float64(data[1][i, 20]) 					# initial condition: water saturation
		Xg = zeros(Float64, num_gas)
		C_g = zeros(Float64, num_comps)
		C_aq = zeros(Float64, num_comps)
		for j = 1:num_gas  								# initial gas composition (mole fraction and molar concentration)
			Xg[j] = Float64(data[1][i, 20+j])
			C_g[j] = Xg[j]*P/(R*T)
		end
		for j = 1:num_comps 							# initial aqueous composition
			C_aq[j] = Float64(data[1][i, 20+num_gas+j])
		end
		for j = 1:num
			# for each like node; connecting node list and sigma terms will be created later from connections
			push!(node, Node(x+(j-1)*x_step, y+(j-1)*y_step, z+(j-1)*z_step, vol, k, 0., 0., phi, n, m, alpha, Sr, Q_index,
				Q, F_index, F, alpha_L, P, 0., S, copy(Xg), copy(C_g), copy(C_aq), Tuple{Int64, Int64}[], 0., 0., 0., copy(sigma_x), 0., 0., 0., 0.))
			node[end].Pc = PCap(node[end])
			num_nodes += 1
		end
	end

	# write to node geometry summary file (k and S also listed as proxies for properties and initial conditions)
	fname = "nodes_geometry.csv"
	csvfile = open(fname,"w")
	line_out = "node" * "," * "x" * "," * "y" * "," * "z" * "," * "volume" * "," * "permeability" *
		"," * "saturation"
	println(csvfile, line_out)
	for i = 1:num_nodes
		line_out = string(i) * "," * string(node[i].x) * "," * string(node[i].y) * "," * string(node[i].z) *
			"," * string(node[i].vol) *	"," * string(node[i].k) * "," * string(node[i].S)
		println(csvfile,line_out)
	end
	close(csvfile)
	println("Processed nodes.")

	return node, num_nodes

end

function ReadConnects(node::Array{Node,1}, num_nodes::Int64, gas::Array{Component,1}, num_gas::Int64, knobs::Knobs)

	# read in connections file and populate node type array
	connect = Connect[]
	num_connects = 0
	data = readdlm("connects.txt", '\t', header=true)
	for i = 1:size(data[1], 1)
		node_1 = Int64(data[1][i, 1])
		node_2 = Int64(data[1][i, 2])
		num = Int64(data[1][i, 3])
		step_1 = Int64(data[1][i, 4])
		step_2 = Int64(data[1][i, 5])
		dx = Float64(data[1][i, 6])
		area = Float64(data[1][i, 7])
		for j = 1:num
			conduct_w = Conduct_w(node, node_1, node_2, dx, area)									# water-phase conductance
			conduct_g = Conduct_g(node, node_1, node_2, dx, area) 									# gas-phase conductance
			conduct_c = Conduct_c(conduct_w, node_1, node_2, node, dx, area, knobs, num_nodes) 		# aqueous chemical transport conductance
			conduct_x = Conduct_x(node, i, j, dx, area, gas, num_gas) 								# gaseous chemical transport conductance (array of size num_gas)
			# for each like connection; note that advective fluxes are assigned later (set to zero initially, here)
			push!(connect, Connect(node_1+(j-1)*step_1, node_2+(j-1)*step_2, dx, area, conduct_w, conduct_g, conduct_c, conduct_x, 0., 0.))
			num_connects += 1
		end
	end

	# write to connection summary file ...
	fname = "connects_geometry.csv"
	csvfile = open(fname,"w")
	line_out = "connection" * "," * "node_1" * "," * "node_2" * "," * "dx" * "," * "area"
	println(csvfile, line_out)
	for i = 1:num_connects
		line_out = string(i) * "," * string(connect[i].node_1) * "," * string(connect[i].node_2) *
			"," * string(connect[i].dx) * "," * string(connect[i].area)
		println(csvfile,line_out)
	end
	close(csvfile)
	println("Processed connections.")

	return connect, num_connects

end

#############################
# unsaturated flow functions
#############################

function SatEff(node_i::Node)::Float64
	# effective water saturation
	return (node_i.S - node_i.Sr)/(1. - node_i.Sr)
end

function PCap(node_i::Node)::Float64
	# capillary pressure as a function of water saturation (used to calculate initial capillary pressures from saturations)
	Se = SatEff(node_i)
	return -exp(log(exp(-log(Se)/node_i.m) - 1.0)/node_i.n - log(node_i.alpha))
end

function Sat(node_i::Node)::Float64
	# water saturation as a function of capillary pressure (update after solving pressure equations)
	Se = (1 + abs(node_i.alpha * node_i.Pc)^node_i.n)^(-node_i.m)
	return Se * (1.0 - node_i.Sr) + node_i.Sr
end

function dSdP(node_i::Node)::Float64
	# change in water saturation per change in capillary pressure
	S_prime = -(node_i.alpha*abs(node_i.Pc))^node_i.n * ((node_i.alpha*abs(node_i.Pc))^node_i.n + 1.0)^(-node_i.m-1.0) * node_i.m * node_i.n/node_i.Pc
	return S_prime * (1.0 - node_i.Sr)
end

function PermRel_w(node_i::Node)::Float64
	# water relative permeability
	Se = SatEff(node_i)
	return Se^0.5 * (1. - (1. - Se^(1.0/node_i.m))^node_i.m)^2
end

function PermRel_g(node_i::Node)::Float64
	# gas relative permeability
	Se = SatEff(node_i)
	return (1. - Se)^0.5 * (1. - Se^(1.0/node_i.m))^(2. * node_i.m)
end

######################################
# utility functions
######################################

function GetMeanGasMW(gas::Array{Component, 1}, num_gas::Int64, S_gas::Array{Float64, 2})::Array{Float64, 1}
	# compute the weighted mean molecular weight for the gas associated with each source term
	MW_gas = zeros(Float64, num_gas)				# gas molecular weight array
	for (i, gs) in enumerate(gas)
		MW_gas[i] = gs.MW
	end
	num_sources = size(S_gas, 1)
	S_gMW = zeros(Float64, num_sources)
	for i = 1:num_sources
		S_gMW[i] = dot(MW_gas, S_gas[i, :])
	end
	return S_gMW
end

function MapNodeConnects(node::Array{Node, 1}, connect::Array{Connect, 1})::Array{Node, 1}
	# create list of node connections per each volume element; will be used to construct matrix of flow equations
	for (i, cn) in enumerate(connect)
		push!(node[cn.node_1].connect_list, (i, cn.node_2))
		push!(node[cn.node_2].connect_list, (i, cn.node_1))
	end
	return node
end

function UpdateRelPerm(node::Array{Node,1})::Array{Node, 1}
	# update relative permeabilities of both water and gases phases across all nodes
	for nd in node
		nd.kr_w = PermRel_w(nd)
		nd.kr_g = PermRel_g(nd)
	end
	return node
end

function Conduct_w(node::Array{Node, 1}, i::Int64, j::Int64, dx::Float64, area::Float64)::Float64
	# compute/update conductance (with respect to unsaturated water flow) between two nodes
	return mean([node[i].kr_w * node[i].k, node[j].kr_w * node[j].k]) / u_w * area / dx
end

function Conduct_g(node::Array{Node, 1}, i::Int64, j::Int64, dx::Float64, area::Float64)::Float64
	# compute/update conductance (with respect to flow of the gas phase) between two nodes
	return mean([node[i].kr_g * node[i].k * node[i].P, node[j].kr_g * node[j].k * node[j].P]) / (u_g * R * T) * area / dx
end

function Conduct_c(conduct_w::Float64, i::Int64, j::Int64, node::Array{Node, 1}, dx::Float64, area::Float64, knobs::Knobs, num_nodes::Int64)::Float64
	# compute/update chemical transport conductance (with respect to the water phase) between two nodes
	dP = zeros(Float64, 2*num_nodes) 																				# empty placeholder for the WaterDarcy function
	v = WaterDarcy(conduct_w, i, j, node, dP, knobs, num_nodes) / (area * mean([node[i].phi, node[j].phi])) 		# pore velocity across interface
	D = abs(v) * mean([node[i].alpha_L, node[j].alpha_L]) 															# (longitudinal) dispersion coefficient
	return mean([node[i].phi * node[i].S, node[j].phi * node[j].S]) * D * area / dx
end

function Conduct_x(node::Array{Node, 1}, i::Int64, j::Int64, dx::Float64, area::Float64, gas::Array{Component, 1}, num_gas::Int64)::Array{Float64, 1}
	# compute/update gas-phase diffusive conductance between two nodes (returns an array of size num_gas)
	cx = zeros(Float64, num_gas)
	for k = 1:num_gas
		cx[k] = mean([node[i].phi * (1. - node[i].S), node[j].phi * (1. - node[j].S)]) * gas[k].D * area / dx
	end
	return cx
end

function UpdateTotConduct(node::Array{Node, 1}, connect::Array{Connect, 1}, num_gas::Int64)::Array{Node, 1}
	# update node total conductances (i.e., across all connections) for water and gas phases
	for nd in node
		nd.sigma_w = 0. 									# clear old values
		nd.sigma_g = 0.
		nd.sigma_c = 0.
		for i = 1:num_gas
			nd.sigma_x[i] = 0.
		end
	end
	for cn in connect
		node[cn.node_1].sigma_w += cn.conduct_w
		node[cn.node_2].sigma_w += cn.conduct_w
		node[cn.node_1].sigma_g += cn.conduct_g
		node[cn.node_2].sigma_g += cn.conduct_g
		node[cn.node_1].sigma_c += cn.conduct_c
		node[cn.node_2].sigma_c += cn.conduct_c
		for i = 1:num_gas
			node[cn.node_1].sigma_x[i] += cn.conduct_x[i]
			node[cn.node_2].sigma_x[i] += cn.conduct_x[i]
		end
	end
	return node
end

function UpdateTotAdvect(node::Array{Node, 1}, connect::Array{Connect, 1})::Array{Node, 1}
	# update node total conductances (i.e., across all connections) for water and gas phases
	for nd in node
		nd.tot_Qw_in = 0. 									# clear old values
		nd.tot_Qw_out = 0.
		nd.tot_Qg_in = 0.
		nd.tot_Qg_out = 0.
	end
	for cn in connect 										# flow is always defined relative to node_1 to node_2

		# water
		if cn.Qw > 0.
			node[cn.node_2].tot_Qw_in += cn.Qw 				# an inflow (positive value)
			node[cn.node_1].tot_Qw_out -= cn.Qw 			# an outflow (negative value)
		else
			node[cn.node_2].tot_Qw_out += cn.Qw 			# an outflow ...
			node[cn.node_1].tot_Qw_in -= cn.Qw 				# an inflow
		end

		# gas
		if cn.Qg > 0.
			node[cn.node_2].tot_Qg_in += cn.Qg
			node[cn.node_1].tot_Qg_out -= cn.Qg
		else
			node[cn.node_2].tot_Qg_out += cn.Qg
			node[cn.node_1].tot_Qg_in -= cn.Qg
		end

	end
	return node
end

function UpdateAdvect(node::Array{Node, 1}, connect::Array{Connect, 1}, dP::Array{Float64, 1}, knobs::Knobs, num_nodes::Int64)::Array{Connect, 1}
	# calculate advective fluxes (mean values per time step, as constrained by knobs.gamma)
	for cn in connect
		cn.Qw = WaterDarcy(cn.conduct_w, cn.node_1, cn.node_2, node, dP, knobs, num_nodes)
		cn.Qg = GasDarcy(cn, node, dP, knobs)
	end
	return connect
end

function GasDarcy(cn::Connect, node::Array{Node, 1}, dP::Array{Float64, 1}, knobs::Knobs)::Float64
	# mean gas volumetric flux across volume element interface (from node_1 into node_2) per current time step
	i = cn.node_1
	j = cn.node_2
	tot_P_i = node[i].P + knobs.gamma*dP[i]
	tot_P_j = node[j].P + knobs.gamma*dP[j]
	return R*T*cn.conduct_g/mean([node[i].P, node[j].P]) * (tot_P_i - tot_P_j)
end

function WaterDarcy(conduct_w::Float64, i::Int64, j::Int64, node::Array{Node, 1}, dP::Array{Float64, 1}, knobs::Knobs, num_nodes::Int64)::Float64
	# mean Darcy flux across volume element interface (from node i into node j) per current time step
	# this function is structured differently than gas-phase equivalent because pore velocity is needed elsewhere in the code
	tot_h_i = node[i].P + knobs.gamma*dP[i] + node[i].Pc + knobs.gamma*dP[i+num_nodes] + rho_w*g*node[i].z
	tot_h_j = node[j].P + knobs.gamma*dP[j] + node[j].Pc + knobs.gamma*dP[j+num_nodes] + rho_w*g*node[j].z
	return conduct_w * (tot_h_i - tot_h_j)
end

######################################
# functions to write to output files
######################################

function WriteFinalState(node::Array{Node, 1}, num_nodes::Int64, gas::Array{Component, 1}, num_gas::Int64, comp::Array{Component, 1}, num_comps::Int64)

	# model-wide output file: header
	fname = "final_state.csv"
	csvfile = open(fname,"w")
	line_out = "node" * "," * "x" * "," * "y" * "," * "z" * "," * "P" * "," * "Pc" * "," * "S"
	for i = 1:num_gas
		line_out = line_out * "," * "X_" * gas[i].name
	end
	for i = 1:num_comps
		line_out = line_out * "," * "C_" * comp[i].name
	end
	println(csvfile, line_out)

	# node-by-node summary
	for i = 1:num_nodes
		line_out = string(i) * "," * string(node[i].x) * "," * string(node[i].y) * "," * string(node[i].z) * "," * string(node[i].P) * "," * string(node[i].Pc) *
			"," * string(node[i].S)
		for j = 1:num_gas
			line_out = line_out * "," * string(node[i].Xg[j])
		end
		for j = 1:num_comps
			line_out = line_out * "," * string(node[i].C_aq[j])
		end
		println(csvfile, line_out)
	end

	close(csvfile)
	println("Wrote final model state to output file.")
	return
end

function MonHeader(csvfile::IOStream, gas::Array{Component, 1}, num_gas::Int64, comp::Array{Component, 1}, num_comps::Int64)
	line_out = "time" * "," * "P" * "," * "Pc" * "," * "S"
	for i = 1:num_gas
		line_out = line_out * "," * "X_" * gas[i].name
	end
	for i = 1:num_comps
		line_out = line_out * "," * "C_" * comp[i].name
	end
	println(csvfile, line_out)
end

function Monitor(t::Float64, node_i::Node, num_comps::Int64, num_gas::Int64, csvfile::IOStream)
	# append current state of monitor node to monitoring well file
	line_out = string(t) * "," * string(node_i.P) * "," * string(node_i.Pc) * "," * string(node_i.S)
	for j = 1:num_gas
		line_out = line_out * "," * string(node_i.Xg[j])
	end
	for j = 1:num_comps
		line_out = line_out * "," * string(node_i.C_aq[j])
	end
	println(csvfile, line_out)
	return
end

############################################
# linear algebra functions: flow/pressure
############################################

function LHSmatrixFlow(connect::Array{Connect,1}, node::Array{Node,1}, knobs::Knobs, dt::Float64, num_nodes::Int64)

	# fill out the LHS of the flow equation matrix and record row-column index positions

	row_index = Int64[] 					# indexing system for sparse matrix
	col_index = Int64[]
	data = Float64[]

	# diagonal elements, or those elements involving the same node but phase interactions
	for (i, nd) in enumerate(node)

		push!(row_index, i)
		push!(col_index, i)
		push!(data, knobs.gamma*nd.sigma_g + nd.phi*(1.-nd.S)*nd.vol/(R*T*dt))			# gas-phase

		push!(row_index, i + num_nodes)
		push!(col_index, i + num_nodes)
		push!(data, knobs.gamma*nd.sigma_w + nd.phi*nd.vol*dSdP(nd)/dt) 				# water-phase

		push!(row_index, i + num_nodes)
		push!(col_index, i)
		push!(data, knobs.gamma * nd.sigma_w)	 										# influence of gas-phase on water

	end

	# non-diagonal elements
	for (i, nd) in enumerate(node)
		for cn in nd.connect_list

			push!(row_index, i)
			push!(col_index, cn[2])
			push!(data, -knobs.gamma * connect[cn[1]].conduct_g) 			# gas-gas (upper left)

			push!(row_index, i + num_nodes)
			push!(col_index, cn[2] + num_nodes)
			push!(data, -knobs.gamma * connect[cn[1]].conduct_w) 			# water-water (lower right)

			push!(row_index, i + num_nodes)
			push!(col_index, cn[2])
			push!(data, -knobs.gamma * connect[cn[1]].conduct_w) 			# gas-water	(lower left)

			# no equations for upper right, since gas mass balance does not (directly) depend on capillary pressure

		end
	end

	return data, row_index, col_index

end

function UpdateFlowLHS(connect::Array{Connect,1}, node::Array{Node,1}, knobs::Knobs, dt::Float64, num_nodes::Int64)

	# fill out the LHS of the flow equation matrix (row-column index positions already have been recorded)

	data = Float64[]

	# diagonal elements, or those elements involving the same-node phase interactions
	for (i, nd) in enumerate(node)
		push!(data, knobs.gamma*nd.sigma_g + nd.phi*(1.-nd.S)*nd.vol/(R*T*dt))			# gas-phase
		push!(data, knobs.gamma*nd.sigma_w + nd.phi*nd.vol*dSdP(nd)/dt) 				# water-phase
		push!(data, knobs.gamma * nd.sigma_w)	 										# influence of gas-phase on water
	end

	# non-diagonal elements
	for (i, nd) in enumerate(node)
		for cn in nd.connect_list
			push!(data, -knobs.gamma * connect[cn[1]].conduct_g) 			# gas-gas (upper left)
			push!(data, -knobs.gamma * connect[cn[1]].conduct_w) 			# water-water (lower right)
			push!(data, -knobs.gamma * connect[cn[1]].conduct_w) 			# gas-water	(lower left)
		end
	end

	return data

end

function RHSvectorFlow(connect::Array{Connect,1}, node::Array{Node,1}, num_nodes::Int64)::Array{Float64, 1}

	# construct explicit flow matrix (run for each time step)
	b = zeros(Float64, 2*num_nodes)

	for i = 1:num_nodes
		b[i] = node[i].F - node[i].sigma_g * node[i].P
		b[i+num_nodes] = node[i].Q - node[i].sigma_w * (rho_w*g*node[i].z + node[i].P + node[i].Pc)
	end

	for cn in connect

		# gas balance
        b[cn.node_1] += cn.conduct_g * node[cn.node_2].P
        b[cn.node_2] += cn.conduct_g * node[cn.node_1].P

		# water balance
        b[cn.node_1 + num_nodes] += cn.conduct_w *
			(rho_w*g*node[cn.node_2].z + node[cn.node_2].P + node[cn.node_2].Pc)
        b[cn.node_2 + num_nodes] += cn.conduct_w *
			(rho_w*g*node[cn.node_1].z + node[cn.node_1].P + node[cn.node_1].Pc)

	end

	return b

end

###############################################
# linear algebra functions: aqueous transport
###############################################

function LHSmatrixTrans(connect::Array{Connect,1}, node::Array{Node,1}, knobs::Knobs, dt::Float64, num_nodes::Int64)

	# fill out the LHS of the advection-dispersion equation matrix and record row-column index positions

	row_index = Int64[] 					# indexing system for sparse matrix
	col_index = Int64[]
	data = Float64[]

	# diagonal elements
	for (i, nd) in enumerate(node)
		push!(row_index, i)
		push!(col_index, i)
		push!(data, nd.phi*nd.S*nd.vol/dt + knobs.gamma*(nd.sigma_c - nd.tot_Qw_in + knobs.upw_w*nd.tot_Qw_in - knobs.upw_w*nd.tot_Qw_out
		+ (node[i].tot_Qw_in + node[i].tot_Qw_out)))
	end

	# non-diagonal elements
	for (i, nd) in enumerate(node)
		for cn in nd.connect_list
			push!(row_index, i)
			push!(col_index, cn[2])
			term = -knobs.gamma * connect[cn[1]].conduct_c
			if (((connect[cn[1]].Qw > 0.) & (connect[cn[1]].node_2 == i)) | ((connect[cn[1]].Qw < 0.) & (connect[cn[1]].node_1 == i)))
				# flow is from node j into node i: this is an inflow-type connection for i
				term += -knobs.gamma * knobs.upw_w * abs(connect[cn[1]].Qw)
			else
				# flow is from node i into node j: this is an outflow-type connection for i
				term += -knobs.gamma * -abs(connect[cn[1]].Qw)  * (1 - knobs.upw_w)
			end
			push!(data, term)
		end
	end

	return data, row_index, col_index

end

function UpdateTransAqLHS(connect::Array{Connect,1}, node::Array{Node,1}, knobs::Knobs, dt::Float64, num_nodes::Int64)

	# fill out the LHS of the advection-dispersion equation matrix and record row-column index positions

	data = Float64[]

	# diagonal elements
	for (i, nd) in enumerate(node)
		push!(data, nd.phi*nd.S*nd.vol/dt + knobs.gamma*(nd.sigma_c - nd.tot_Qw_in + knobs.upw_w*nd.tot_Qw_in - knobs.upw_w*nd.tot_Qw_out
		+ (node[i].tot_Qw_in + node[i].tot_Qw_out)))
	end

	# non-diagonal elements
	for (i, nd) in enumerate(node)
		for cn in nd.connect_list
			term = -knobs.gamma * connect[cn[1]].conduct_c
			if (((connect[cn[1]].Qw > 0.) & (connect[cn[1]].node_2 == i)) | ((connect[cn[1]].Qw < 0.) & (connect[cn[1]].node_1 == i)))
				# flow is from node j into node i: this is an inflow-type connection for i
				term += -knobs.gamma * knobs.upw_w * abs(connect[cn[1]].Qw)
			else
				# flow is from node i into node j: this is an outflow-type connection for i
				term += -knobs.gamma * -abs(connect[cn[1]].Qw)  * (1 - knobs.upw_w)
			end
			push!(data, term)
		end
	end

	return data
end

function RHSvectorTransAq(connect::Array{Connect,1}, node::Array{Node,1}, S_aq::Array{Float64, 2}, num_nodes::Int64, icomp::Int64, knobs::Knobs)::Array{Float64, 1}

	# construct explicit advective-dispersive matrix for aqueous species chemical transport (run for each time step)
	b = zeros(Float64, num_nodes)

	for i = 1:num_nodes 			# terms written with respect to concentration in node i
		b[i] = node[i].Q*S_aq[node[i].Q_index, icomp] + node[i].C_aq[icomp]*(-node[i].sigma_c + (1-knobs.upw_w)*node[i].tot_Qw_in + knobs.upw_w*node[i].tot_Qw_out
			- (node[i].tot_Qw_in + node[i].tot_Qw_out)) 								# special provision for change in water saturation
	end

	for cn in connect 				# terms written with respect to concentration in node j
		if cn.Qw > 0.
			# flow is from node_1 into node_2
			b[cn.node_1] += (cn.conduct_c + (1-knobs.upw_w)*(-cn.Qw)) * node[cn.node_2].C_aq[icomp] 			# only outflows (negative) from node i (node_1)
			b[cn.node_2] += (cn.conduct_c + knobs.upw_w*cn.Qw) * node[cn.node_1].C_aq[icomp] 					# only inflows (positive) into node i (node_2)
		else
			# flow is from node_2 into node_1
			b[cn.node_1] += (cn.conduct_c + knobs.upw_w*(-cn.Qw)) * node[cn.node_2].C_aq[icomp] 				# only inflows (positive) into node i (node_2)
			b[cn.node_2] += (cn.conduct_c + (1-knobs.upw_w)*cn.Qw) * node[cn.node_1].C_aq[icomp] 				# only outflows (negative) from node i (node_1)
		end
	end

	return b

end

###############################################
# linear algebra functions: gas transport
###############################################

function UpdateTransGasLHS(connect::Array{Connect,1}, node::Array{Node,1}, knobs::Knobs, dt::Float64, num_nodes::Int64, igas::Int64)

	# fill out the LHS of the advection-dispersion equation matrix for a gas component

	data = Float64[]

	# diagonal elements
	for (i, nd) in enumerate(node)
		push!(data, nd.phi*(1.-nd.S)*nd.vol/dt + knobs.gamma*(nd.sigma_x[igas] - nd.tot_Qg_in + knobs.upw_g*nd.tot_Qg_in - knobs.upw_g*nd.tot_Qg_out
		- (node[i].tot_Qw_in + node[i].tot_Qw_out)))
	end

	# non-diagonal elements
	for (i, nd) in enumerate(node)
		for cn in nd.connect_list
			term = -knobs.gamma * connect[cn[1]].conduct_x[igas]
			if (((connect[cn[1]].Qg > 0.) & (connect[cn[1]].node_2 == i)) | ((connect[cn[1]].Qg < 0.) & (connect[cn[1]].node_1 == i)))
				# flow is from node j into node i: this is an inflow-type connection w.r.t. node i
				term += -knobs.gamma * knobs.upw_g * abs(connect[cn[1]].Qg)
			else
				# flow is from node i into node j: this is an outflow-type connection w.r.t. node i
				term += -knobs.gamma * -abs(connect[cn[1]].Qg)  * (1 - knobs.upw_g)
			end
			push!(data, term)
		end
	end

	return data

end

function RHSvectorTransGas(connect::Array{Connect,1}, node::Array{Node,1}, S_gas::Array{Float64, 2}, num_nodes::Int64, gas::Array{Component, 1}, igas::Int64, knobs::Knobs)::Array{Float64, 1}

	# construct explicit advective-dispersive matrix for gas species chemical transport (run for each time step)
	b = zeros(Float64, num_nodes)

	for i = 1:num_nodes
		b[i] = node[i].F*S_gas[node[i].F_index, igas] + node[i].C_g[igas]*(-node[i].sigma_x[igas] + (1-knobs.upw_g)*node[i].tot_Qg_in + knobs.upw_g*node[i].tot_Qg_out
		+ (node[i].tot_Qw_in + node[i].tot_Qw_out))
	end

	for cn in connect
		if cn.Qg > 0.
			# flow is from node_1 into node_2
			b[cn.node_1] += (cn.conduct_x[igas] + (1-knobs.upw_g)*(-cn.Qg)) * node[cn.node_2].C_g[igas] 		# only outflows (negative) from node i (node_1)
			b[cn.node_2] += (cn.conduct_x[igas] + knobs.upw_g*cn.Qg) * node[cn.node_1].C_g[igas] 				# only inflows (positive) into node i (node_2)
		else
			# flow is from node_2 into node_1
			b[cn.node_1] += (cn.conduct_x[igas] + knobs.upw_g*(-cn.Qg)) * node[cn.node_2].C_g[igas] 				# only inflows (positive) into node i (node_2)
			b[cn.node_2] += (cn.conduct_x[igas] + (1-knobs.upw_g)*cn.Qg) * node[cn.node_1].C_g[igas] 				# only outflows (negative) from node i (node_1)
		end
	end

	return b

end

######################################
# main script
######################################

function gas_water_trans(t_end::Float64, mon_well::Int64)

	knobs = GetKnobs() 																# read in computational parameters
	comp, gas, num_comps, num_gas = ReadComps() 									# read chemical component properties
	S_gas = ReadSource(1)															# read sources
	S_aq = ReadSource(2)
	S_gMW = GetMeanGasMW(gas, num_gas, S_gas)										# weighted average gas molecular weight (array, one per source)
	node, num_nodes = ReadNodes(num_comps, num_gas, S_gMW)							# read nodes file and organize
	node = UpdateRelPerm(node) 														# populate relative permeabilities for both water and gas phases
	connect, num_connect = ReadConnects(node, num_nodes, gas, num_gas, knobs)		# read connections file and organize
	node = MapNodeConnects(node, connect) 											# create lists of connected nodes, per volume element
	node = UpdateTotConduct(node, connect, num_gas)									# populate total conductances, per node

	PP = zeros(Float64, num_gas) 													# temporary array for gas partial pressures, per node
	t = 0.
	dt = knobs.dt_init

	# open monitor file and write header
	fname = "monitor.csv"
	csvfile = open(fname,"w")
	MonHeader(csvfile, gas, num_gas, comp, num_comps)

	# set up matrix indexing system for left-hand-side of flow/pressure balance equations, "p_data" is really a placeholder
	p_data, p_row_index, p_col_index = LHSmatrixFlow(connect, node, knobs, dt, num_nodes)

	# set up matrix indexing system for left-hand-side of transport equations, "c_data" is really a placeholder
	# same indexing system is used for aqueous component as well as gaseous component transport equations
	c_data, c_row_index, c_col_index = LHSmatrixTrans(connect, node, knobs, dt, num_nodes)

    while (t < t_end)

		### solve the flow/pressure portion ###

        advect_complete = false

        while advect_complete == false

            p_data = UpdateFlowLHS(connect, node, knobs, dt, num_nodes) 				# update LHS elements (new dt, new conductances)
 			A = sparse(p_row_index, p_col_index, p_data, 2*num_nodes, 2*num_nodes)		# update sparse flow equation matrix
            b = RHSvectorFlow(connect, node, num_nodes)									# construct explicit vector
			global dP = \(A, b)	 														# solve equation set (1st half of dP is gas pressure, 2nd is capillary pressure)

            # check maximum pressure criterion (applied to both water and gas phases) at this time step size
			sum_complete = 0
			for i = 1:num_nodes
				sum_complete += (abs(dP[i]) > knobs.dP_max) + (abs(dP[i+num_nodes]) > knobs.dPc_max)
			end
			advect_complete = 1 - sign(sum_complete)
            if advect_complete == false 				# reduce time step size if pressure change criterion not satisfied
				dt *= knobs.dt_decrease
				assert(dt > knobs.dt_min)
			end

		end

		# update advective flows; implemented here to use most current P and dP estimates
		connect = UpdateAdvect(node, connect, dP, knobs, num_nodes)
		node = UpdateTotAdvect(node, connect)

		### solve the advection-dispersion equation: aqueous components ###
		c_data = UpdateTransAqLHS(connect, node, knobs, dt, num_nodes) 						# update LHS elements (new dt, new conductances)
		A = sparse(c_row_index, c_col_index, c_data, num_nodes, num_nodes)					# update sparse transport equation matrix
		for icomp = 1:num_comps
			b = RHSvectorTransAq(connect, node, S_aq, num_nodes, icomp, knobs)				# construct explicit vector
			dC = \(A, b)
			# update concentrations
			for (i, nd) in enumerate(node)
				nd.C_aq[icomp] += dC[i]
			end
		end

		### solve the advection-diffusion equation: gaseous components ###
		for igas = 1:num_gas 																# note: LHS for gas equations is component-specific because of differing diffusion coefficients
			c_data = UpdateTransGasLHS(connect, node, knobs, dt, num_nodes, igas) 			# update LHS elements (new dt, new conductances)
			A = sparse(c_row_index, c_col_index, c_data, num_nodes, num_nodes)				# update sparse transport equation matrix
			b = RHSvectorTransGas(connect, node, S_gas, num_nodes, gas, igas, knobs)		# construct explicit vector
			dC = \(A, b)
			# update concentrations
			for (i, nd) in enumerate(node)
				nd.C_g[igas] += dC[i]
			end
		end

        # apply updates after concluding current time step
        t += dt 																# increment simulation time
		for (i, nd) in enumerate(node)
			nd.Pc += dP[i+num_nodes] 											# update capillary pressure
			nd.S = Sat(nd)		 												# update water saturation implied by capillary pressure
			nd.P = 0.
			for igas = 1:num_gas 												# calculate gas partial pressures from ideal gas law
				PP[igas] = nd.C_g[igas]*R*T 									# from P = (n/V)*R*T
				nd.P += PP[igas]
			end
			for igas = 1:num_gas 												# update gas mole fractions (for model output)
				nd.Xg[igas] = PP[igas]/nd.P
			end
		end
		node = UpdateRelPerm(node) 												# relative permeabilities for water and gas phases
		Monitor(t, node[mon_well], num_comps, num_gas, csvfile) 				# append to monitoring well report

		# update conductances (prior to next time step)
		for cn in connect
			cn.conduct_w = Conduct_w(node, cn.node_1, cn.node_2, cn.dx, cn.area)
			cn.conduct_g = Conduct_g(node, cn.node_1, cn.node_2, cn.dx, cn.area)
			cn.conduct_c = Conduct_c(cn.conduct_w, cn.node_1, cn.node_2, node, cn.dx, cn.area, knobs, num_nodes)
			cn.conduct_x = Conduct_x(node, cn.node_1, cn.node_2, cn.dx, cn.area, gas, num_gas)
		end

		# update total conductances, per node
		node = UpdateTotConduct(node, connect, num_gas)

        # update time step
        dt *= knobs.dt_increase
        dt = min(dt, knobs.dt_max, t_end - t)

	end

	close(csvfile) 																# close monitor file
	WriteFinalState(node, num_nodes, gas, num_gas, comp, num_comps) 			# write final state of system to output file

	println("Done.")

end

### run script as gas_water_trans(t_end, mon_well)

# t_end = model end-time
# mon_well = node index number corresponding to monitor location (node no.)

gas_water_trans(8640., 2)


text:
using IntcodeVM, Sockets, OffsetArrays

const code = IntcodeVM.load_intcode("input.txt")
const queues = OffsetArray([[(i, -1)] for i=0:49], 0:49)

# Spin up 50 instances of the VM.
foreach(i -> IntcodeVM.run_async(code, 60000 + i), 0:49)

# Connect to all 50 VMs on TCP.
const sockets = OffsetArray([connect(60000 + i) for i=0:49], 0:49)

# Part 1 asks for the value when we first encounter destination address 255.
first_255_seen = false

# The current (x,y) values sent to 255. This value can be overwritten.
nat = (missing, missing)

# The most recent NAT value sent to address 0 when the network was idle.
last_nat = undef

# Switch values read from a VM to the appropriate destination.
# These writes are not atomic. I could envision a race condition where
# the main loop sends a packet to a VM, then moves to the next queue
# and sends another packet to a different VM, and the second response
# makes it to the switching thread first. In practice, however,
# this function is good enough to reliably get the correct answer.
function switch(socket::IO)
    @async begin
        while isopen(socket)
            (dst, x, y) = [parse(Int, readline(socket)) for i=1:3]
            if dst != 255
                push!(queues[dst], (x,y))
            elseif !first_255_seen && dst == 255
                println("Day 23 Part 1: $(y)")
                global first_255_seen = true
                global nat = (x,y)
            else
                global nat = (x,y)
            end
        end
    end
end

# Start the switching threads.
foreach(switch, sockets)

# Assign addresses to each VM.
foreach(i -> println(sockets[i], i), 0:49)

# Iterate across all queues repeatedly.
# Terminate the outer loop when we find the part 2 answer.
for i=0:1000
    idle = true

    # Iterate over each queue, sending (x,y) "packets" if present and -1 otherwise.
    for j=0:49
        if !isempty(queues[j])
            (x,y) = popfirst!(queues[j])
            println(sockets[j], x)
            println(sockets[j], y)
            idle = false
        else
            println(sockets[j], -1)
        end
    end

    # The network is considered "idle" when all queues are empty.
    if idle
        if nat == last_nat
            println("Day 23 Part 2: $(nat[2])")
            break
        else
            push!(queues[0], nat)
            global last_nat = nat
        end
    end
end

foreach(close, sockets)


text:
import Base: ==, >, <, >=, <=, !=

for (op, pyop) in [(:(==), :__eq__), (:>, :__gt__), (:<, :__lt__), (:>=, :__ge__), (:<=, :__le__), (:!=, :__ne__)]
    @eval function Base.broadcast(::typeof($op), s::PandasWrapped, x)
        method = s.pyo.$(QuoteNode(pyop))
        pandas_wrap(pycall(method, PyObject, x))
    end
end


text:
include("diffusion.jl")

function assembleDiffusionMatrix(gamma, phi, gradPhi)
    println("Assembling the diffusion matrix")
    fluxCf, fluxFf, fluxVf = assembleDiffusionTerm(gamma, phi, gradPhi)

    A = zeros(mesh.nCells, mesh.nCells)
    B = zeros(mesh.nCells)

    for iFace in 1:mesh.nFaces
        iOwner = mesh.faces[iFace].iOwner
        iNeighbour = mesh.faces[iFace].iNeighbour
        A[iOwner, iOwner] += fluxCf[iFace]
        B[iOwner] += -fluxVf[iFace]
        if iNeighbour != -1
            A[iNeighbour, iNeighbour] += fluxCf[iFace]
            A[iOwner, iNeighbour] += fluxFf[iFace]
            A[iNeighbour, iOwner] += fluxFf[iFace]
        end
    end
    return A, B
end


text:
include("../common.jl")

# Dependencies that must be installed before this package can be built

dependencies = [
    BuildDependency(PackageSpec(name="LLVM_full_jll", version=v"9.0.1")),
#    Dependency(PackageSpec(name="libLLVM_jll", version=v"9.0.1"))
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"8", julia_compat="1.5")

max_stars_repo_name: Affie/Yggdrasil

max_issues_repo_name: Affie/Yggdrasil

max_forks_repo_name: Affie/Yggdrasil


text:
# ## Building

# ## Model

module building #jl

using ReachabilityAnalysis, MAT
using ReachabilityModels: @relpath

file = matopen(@relpath "building.mat")

## system matrix
A = read(file, "A")

## input matrix
B = read(file, "B")

## state domain
X = Universe(48)

## input domain
U = BallInf([0.5], 0.3)

function model(X0)
    S = @system(x' = Ax + Bu, x ∈ X, u ∈ U)
    return IVP(S, X0)
end

end  #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
"""
	Settings(;arg=val)

Creates a COSMO settings object that is used to pass user settings to the solver.

Argument | Description | Values (default)
--- | --- | ---
rho | ADMM rho step | 0.1
sigma | ADMM sigma step | 1e-6
alpha | Relaxation parameter | 1.6
eps_abs | Absolute residual tolerance | 1e-4
eps_rel | Relative residual tolerance | 1e-4
eps\\_prim\\_inf | Primal infeasibility tolerance | 1e-4
eps\\_dual\\_inf | Dual infeasibility tolerance | 1e-4
max_iter | Maximum number of iterations | 2500
verbose | Verbose printing | false
verbose_timing | Verbose timing | false
check_termination | Check termination interval | 40
check_infeasibility | Check infeasibility interval | 40
scaling | Number of scaling iterations | 10
adaptive_rho | Automatic adaptation of step size parameter | true
time_limit | set solver time limit in s | 0
"""
mutable struct Settings
	rho::Float64
	sigma::Float64
	alpha::Float64
	eps_abs::Float64
	eps_rel::Float64
	eps_prim_inf::Float64
	eps_dual_inf::Float64
	max_iter::Int64
	verbose::Bool
	check_termination::Int64
	check_infeasibility::Int64
	scaling::Int64
	MIN_SCALING::Float64
	MAX_SCALING::Float64
	adaptive_rho::Bool
	adaptive_rho_interval::Int64
	adaptive_rho_tolerance::Float64
	verbose_timing::Bool
	RHO_MIN::Float64
	RHO_MAX::Float64
	RHO_TOL::Float64
	time_limit::Float64
	obj_true::Float64
	obj_true_tol::Float64
	#constructor
	function Settings(;
		rho=0.1,
		sigma=1e-6,
		alpha=1.6,
		eps_abs=1e-4,
		eps_rel=1e-4,
		eps_prim_inf=1e-6,
		eps_dual_inf=1e-4,
		max_iter=2500,
		verbose=false,
		check_termination=40,
		check_infeasibility=40,
		scaling=10,
		MIN_SCALING = 1e-4,
		MAX_SCALING = 1e4,
		adaptive_rho = true,
		adaptive_rho_interval = 40,
		adaptive_rho_tolerance = 5,
		verbose_timing = false,
		RHO_MIN = 1e-6,
		RHO_MAX = 1e6,
		RHO_TOL = 1e-4,
		time_limit = 0.0,
		obj_true = NaN,
		obj_true_tol = 1e-3
		)
	new(rho, sigma, alpha, eps_abs, eps_rel, eps_prim_inf, eps_dual_inf, max_iter, verbose,  check_termination, check_infeasibility, scaling, MIN_SCALING, MAX_SCALING, adaptive_rho, adaptive_rho_interval, adaptive_rho_tolerance, verbose_timing, RHO_MIN, RHO_MAX, RHO_TOL, time_limit, obj_true, obj_true_tol)
end
end


text:
using CSV, JSON, DataFrames

# --- Utility ---

flatten(x::Array{<:Array, 1}) = Iterators.flatten(x)|> collect|> flatten
flatten(x::Array{<:Number, 1}) = x

shape(x::Array{<:Array, 1}) = vcat(shape(first(x)), [length(x)])
shape(x::Array{<:Number, 1}) = [length(x)]

mdim(x::Array{<:Array, 1}) = reshape(flatten(x), shape(x)...)
mdim(x::Array{<:Number, 1}) = x

transform(x::Array, ::Type{Array{T, N}}) where T <: Number where N = mdim(x)
transform(x::Array, ::Type{Array{T, 1}}) where T <: Array = x
transform(x::Number, ::Type{T}) where T <: Number = x

function convert_type(::Type{Array{T, N}}) where T <: Number where N
    reduce((t,_) -> Array{t, 1}, 1:N; init=T)
end
convert_type(t::Type{Array{T, 1}}) where T <: Array = t
convert_type(t::Type{T}) where T <: Number = t

# -- IO ---

"""Saves any JSON serializable `object` into JSON file in `filepath`."""
function save_json(object::Any, filepath::AbstractString)
    open(filepath, "w") do io
        JSON.print(io, object)
    end
end

"""Loads values from JSON file in `filepath` to DataType supplied by `dtype`."""
function load_json(dtype::DataType, filepath::AbstractString)
    objects = JSON.parsefile(filepath)
    fields = (
        objects[string(s)] |>
        v -> convert(convert_type(t), v) |>
        v -> transform(v, t)
        for (s, t) in zip(fieldnames(dtype), fieldtypes(dtype)))
    dtype(fields...)
end

"""Load sets, subsets and parameters from CSV files.

Each CSV file in `product_path` contains following attributes:

- `product_id`
- `category_id`
- `brand_id`
- `width`
- `height`
- `depth`
- `weight`
- `monthly_demand`
- `replenishment_interval`
- `price`
- `unit_margin`
- `blocking_field`
- `min_facing`
- `max_facing`
- `max_stack`
- `up_down_order_criteria`

Each CSV file in `shelf_path` contains following attributes:

- `module`
- `id`
- `level`
- `total_width`
- `total_height`
- `total_length`
- `product_min_unit_weight`
- `product_max_unit_weight`

Examples of input parameters can found inside `examples/instances` directory.

# Arguments
- `product_path::AbstractString`
- `shelf_path::AbstractString`
"""
function Params(product_path::AbstractString, shelf_path::AbstractString)
    # Load data from CSV files. Data is read into a DataFrame
    product_data = CSV.read(product_path, DataFrame)
    shelf_data = CSV.read(shelf_path, DataFrame)

    # Sets and Subsets
    products = 1:size(product_data, 1)
    shelves = 1:size(shelf_data, 1)

    # Blocks
    bfs = product_data.blocking_field
    P_b = [collect(products)[bfs .== bf] for bf in unique(bfs)]
    blocks = 1:size(P_b, 1)

    # Modules
    mds = shelf_data.module
    S_m = [collect(shelves)[mds .== md] for md in unique(mds)]
    modules = 1:size(S_m, 1)

    # Return Params struct
    Params(
        products = products,
        shelves = shelves,
        blocks = blocks,
        modules = modules,
        P_b = P_b,
        S_m = S_m,
        G_p = product_data.unit_margin,
        H_s = shelf_data.total_height,
        L_p = product_data.up_down_order_criteria,
        P_ps = transpose(shelf_data.total_length) ./ product_data.depth,
        D_p = product_data.monthly_demand,
        N_p_min = product_data.min_facing,
        N_p_max = product_data.max_facing,
        W_p = product_data.width,
        W_s = shelf_data.total_width,
        M_p = product_data.weight,
        M_s_min = shelf_data.product_min_unit_weight,
        M_s_max = shelf_data.product_max_unit_weight,
        R_p = product_data.replenishment_interval,
        L_s = shelf_data.level,
        H_p = product_data.height,
        SK_p = product_data.max_stack
    )
end


text:
"""
Allow to work as a server
"""

using HTTP
using ZipFile
using Sockets
# using Base64

using MoodleQuestions
using SimpleTranslations

"""
get_params(req::HTTP.Request)

Returns the POST parameters in a Dictionary.
"""
function get_params(req::HTTP.Request)
    data = String(req.body)
    params = Dict{String,String}()

    headers = filter(x->x[1] == "Content-Type", req.headers)
    # Check there is content-type
    if (isempty(headers))
        content_type = ""
    else
        content_type = only(headers)[2]
    end

    if (content_type == "application/x-www-form-urlencoded")
        pairs = split.(split(data, "&"), Ref("="))
        params = Dict(name => HTTP.URIs.unescapeuri(value) for (name, value) in pairs)
    elseif (startswith(content_type, "multipart/form-data"))
        body = IOBuffer(data)

        # First line is the stopping line
        stopping_line = replace(readline(body), "-" => "")

        while(!eof(body))
            line = readline(body)
            text_head = r"Content-Disposition: form-data; name=\"(.*)\""
            m = match(text_head, line)

            if !isnothing(m)
                name = m.captures[1]
                _ = readline(body)
                line = readline(body)
                value = ""

                if (!occursin(stopping_line, line))
                    value = "$(line)"
                end

                while (!eof(body) && !occursin(stopping_line, line))
                    line = readline(body)

                    if (!occursin(stopping_line, line))
                        value *= "\n$(line)"
                    end
                end

                params[name] = value
            end
        end
    end
    return params
end

"""
replace_accent(string)

Return the string without accents (for XML filename saving).
"""
function replace_accent(string)
    string = replace(string, "á" => "a")
    string = replace(string, "é" => "e")
    string = replace(string, "í" => "i")
    string = replace(string, "ó" => "o")
    string = replace(string, "ú" => "u")
    return string
end

function response_error(msg)
    return HTTP.Response(500, [],  body=msg)
end

function handle(req::HTTP.Request)
    params = get_params(req)
    penalty_boolean = tryparse(Float32, get(params, "penalty_boolean", "0"))
    penalty_options = tryparse(Float32, get(params, "penalty_options", "0"))

    language = get(params, "language", "es")
    set_language!(language)

    if (isnothing(penalty_boolean) || isnothing(penalty_options))
        return response_error
    end

    content = get(params, "text", "")

    try
        quiz = read_txt(IOBuffer(content))
       dir = mktempdir()
    # Save quiz to temp dir
    save_to_moodle(quiz, joinpath(dir, "quiz.xml"), penalty_options=penalty_options, penalty_boolean=penalty_boolean)
    # Get list
    files = readdir(dir, join=true)

    if length(files) == 1
        fname = files[1]
        content_type = "application/xml"
    else
        fname = joinpath(dir, "quiz.zip")
        zipfile = ZipFile.Writer(fname)

        for fname in files
            nameinzip = replace_accent(basename(fname))
            fileinzip = ZipFile.addfile(zipfile, nameinzip)
            open(fname) do file
                data = read(file, String)
                write(fileinzip, data)
            end
        end

        close(zipfile)
        content_type = "application/zip"
    end

    filename = basename(fname)

    open(fname) do file
        content = read(file, String)
    end

    headers = ["Pragma" => "public", "Expires" => "0", "Content-Description" => "File Transfer",
                    "Content-type" => content_type, "Content-Transfer-Encoding" => "binary",
               "Content-Length" => "$(sizeof(content))", "Content-Disposition" => "attachment; filename=\"$(filename)\""]
        return HTTP.Response(200, headers, body=content)

    catch e
        if (:msg in propertynames(e))
            return response_error(e.msg)
        else
            @show e
            return response_error("Unknown error")
        end
    end
end

"""
serve_quiz(port)

Run the text file in the txt format, and return the XML file.
If there is only one category only one parameter is defined.

serve_quiz(port)
"""
function serve_quiz(port = 8100)
    router = HTTP.Router()
    fmessages = joinpath(dirname(pathof(MoodleQuestions)), "messages.ini")
    loadmsgs!(fmessages, strict_mode=true)
    HTTP.@register(router, "POST", "/*", handle)
    HTTP.serve(router, Sockets.getipaddr(), port)
end

function test_serve(port = 8080)
    open("/tmp/recibido.zip", "w") do io
        HTTP.request("POST", "http://127.0.0.1:$port/send", Dict("ver"=>"texto"), response_stream=io)
    end
end

max_stars_repo_name: dmolina/MoodleQuestions

max_issues_repo_name: dmolina/MoodleQuestions

max_forks_repo_name: dmolina/MoodleQuestions