text:
module CbcCInterface
using Compat: Cvoid, @compat
using Compat.SparseArrays

import ..Cbc

export CbcModel,
    loadProblem,
    readMps,
    writeMps,
    setInitialSolution,
    problemName,
    setProblemName,
    getNumElements,
    getVectorStarts,
    getIndices,
    getElements,
    getNumRows,
    getNumCols,
    setObjSense,
    getObjSense,
    getRowLower,
    getRowUpper,
    getRowActivity,
    getColLower,
    getColUpper,
    getObjCoefficients,
    getColSolution,
    setRowUpper,
    setRowLower,
    setObjCoeff,
    setColLower,
    setColUpper,
    isInteger,
    setContinuous,
    setInteger,
    setParameter,
    solve,
    sumPrimalInfeasibilities,
    numberPrimalInfeasibilities,
    getIterationCount,
    optimizeNotCalled,
    isAbandoned,
    isProvenOptimal,
    isProvenInfeasible,
    isContinuousUnbounded,
    isNodeLimitReached,
    isSecondsLimitReached,
    isSolutionLimitReached,
    isInitialSolveAbandoned,
    isInitialSolveProvenOptimal,
    isInitialSolveProvenPrimalInfeasible,
    getObjValue,
    getBestPossibleObjValue,
    getNodeCount,
    status,
    secondaryStatus,
    addSOS

# helper macros/functions

macro cbc_ccall(func, args...)
    args = map(esc,args)
    f = "Cbc_$(func)"
    quote
        ccall(($f,Cbc.libcbcsolver), $(args...))
    end
end

mutable struct CbcModel
    p::Ptr{Cvoid}
    function CbcModel()
        p = @cbc_ccall newModel Ptr{Cvoid} ()
        prob = new(p)
        @compat finalizer(deleteModel, prob)
        return prob
    end
end

# This is defined so that we can pass a CbcModel as an argument to ccall without worrying about it being garbage collected during the call.
Base.unsafe_convert(::Type{Ptr{Cvoid}}, prob::CbcModel) = prob.p

function deleteModel(prob::CbcModel)
    if prob.p == C_NULL
        return
    end
    @cbc_ccall deleteModel Cvoid (Ptr{Cvoid},) prob
    prob.p = C_NULL
    return
end

function check_problem(prob::CbcModel)
   if prob.p == C_NULL
       error("Invalid CbcModel")
   end
   return true
end

macro getproperty(T, name)
    @eval function ($name)(prob::CbcModel)
        check_problem(prob)
        @cbc_ccall $name $T (Ptr{Cvoid},) prob
    end
end

# Note: we assume COIN_BIG_INDEX and COIN_BIG_DOUBLE
# were not defined when compiling Cbc (which is the
# default)
const CoinBigIndex = Int32
const CoinBigDouble = Float64

# copied from Clp.jl
const VecOrNothing = Union{Vector,Cvoid}
function vec_or_null(::Type{T}, a::VecOrNothing, len::Integer) where T
    if isequal(a, nothing) || length(a) == 0
        return C_NULL
    else # todo: helpful message if convert fails
        if length(a) != len
            error("Expected vector to have length $len")
        end
        return convert(Vector{T},a)
    end
end

function getVersion()
    s = @cbc_ccall getVersion Ptr{UInt8} ()
    return string(s)
end

function loadProblem(prob::CbcModel,
    constraint_matrix::AbstractMatrix,
    col_lb::VecOrNothing, col_ub::VecOrNothing,
    obj::VecOrNothing,
    row_lb::VecOrNothing, row_ub::VecOrNothing)
    check_problem(prob)

    mat = convert(SparseMatrixCSC{Float64,Int32},constraint_matrix)
    nrow,ncol = size(mat)

    @cbc_ccall loadProblem Cvoid (Ptr{Cvoid}, Int32, Int32, Ptr{CoinBigIndex},
        Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
        Ptr{Float64}, Ptr{Float64}, Ptr{Float64}) prob ncol nrow mat.colptr.-Int32(1) mat.rowval.-Int32(1) mat.nzval vec_or_null(Float64, col_lb, ncol) vec_or_null(Float64, col_ub, ncol) vec_or_null(Float64, obj, ncol) vec_or_null(Float64, row_lb, nrow) vec_or_null(Float64, row_ub, nrow)
end

function readMps(prob::CbcModel, filename::String)
    check_problem(prob)
    @assert isascii(filename)
    @cbc_ccall readMps Cint (Ptr{Cvoid}, Ptr{UInt8}) prob filename
end

function writeMps(prob::CbcModel, filename::String)
    check_problem(prob)
    @assert isascii(filename)
    @cbc_ccall writeMps Cint (Ptr{Cvoid}, Ptr{UInt8}) prob filename
end

function setInitialSolution(prob::CbcModel, array::Vector{Float64})
    check_problem(prob)
    @cbc_ccall setInitialSolution Cvoid (Ptr{Cvoid}, Ptr{Float64}) prob array
end

function problemName(prob::CbcModel)
    check_problem(prob)
    a = Array(UInt8, 100)
    @cbc_ccall problemName Cvoid (Ptr{Cvoid},Cint,Ptr{UInt8}) prob 100 a
    return string(a)
end

#function setProblemName(prob::CbcModel)
#    check_problem(prob)
#    @cbc_ccall setProblemName Cint (Ptr{Cvoid},Ptr{UInt8}) prob bytestring(a)
#end

function getNumElements(prob::CbcModel)
    check_problem(prob)
    @cbc_ccall getNumElements Cint (Ptr{Cvoid},) prob
end

function getVectorStarts(prob::CbcModel)
    check_problem(prob)
    p = @cbc_ccall getVectorStarts Ptr{CoinBigIndex} (Ptr{Cvoid},) prob
    num_cols = Int(getNumCols(prob))
    return copy(unsafe_wrap(Array,p,(num_cols+1,)))
end

function getIndices(prob::CbcModel)
    check_problem(prob)
    p = @cbc_ccall getIndices Ptr{Cint} (Ptr{Cvoid},) prob
    nnz = Int(getNumElements(prob))
    return copy(unsafe_wrap(Array,p,(nnz,)))
end

function getElements(prob::CbcModel)
    check_problem(prob)
    p = @cbc_ccall getElements Ptr{Float64} (Ptr{Cvoid},) prob
    nnz = Int(getNumElements(prob))
    return copy(unsafe_wrap(Array,p,(nnz,)))
end

# TODO: maxNameLength getRowName getColName setColName setRowName

@getproperty Cint getNumRows
@getproperty Cint getNumCols

# 1 : minimize, -1 : maximize
function setObjSense(prob::CbcModel, sense)
    check_problem(prob)
    @cbc_ccall setObjSense Cvoid (Ptr{Cvoid}, Float64) prob sense
end

@getproperty Float64 getObjSense

for s in (:getRowLower, :getRowUpper, :getRowActivity)
    @eval function ($s)(prob::CbcModel)
        check_problem(prob)
        nrow = Int(getNumRows(prob))
        p = @cbc_ccall $s Ptr{Float64} (Ptr{Cvoid},) prob
        return copy(unsafe_wrap(Array,p,(nrow,)))
    end
end

for s in (:getColLower, :getColUpper, :getObjCoefficients, :getColSolution)
    @eval function ($s)(prob::CbcModel)
        check_problem(prob)
        ncol = Int(getNumCols(prob))
        p = @cbc_ccall $s Ptr{Float64} (Ptr{Cvoid},) prob
        return copy(unsafe_wrap(Array,p,(ncol,)))
    end
end

for s in (:setRowUpper, :setRowLower, :setObjCoeff, :setColLower, :setColUpper)
    @eval function($s)(prob::CbcModel, index::Integer, value::Float64)
        check_problem(prob)
        @cbc_ccall $s Cvoid (Ptr{Cvoid}, Cint, Float64) prob index value
    end
end

function isInteger(prob::CbcModel, index::Integer)
    check_problem(prob)
    v = @cbc_ccall isInteger Cint (Ptr{Cvoid},Cint) prob index
    return v == 1
end

function setContinuous(prob::CbcModel, index::Integer)
    check_problem(prob)
    @cbc_ccall setContinuous Cvoid (Ptr{Cvoid},Cint) prob index
end

function setInteger(prob::CbcModel, index::Integer)
    check_problem(prob)
    @cbc_ccall setInteger Cvoid (Ptr{Cvoid},Cint) prob index
end

function Base.copy(prob::CbcModel)
    p = @cbc_ccall clone Ptr{Cvoid} (Ptr{Cvoid},) prob
    prob = CbcModel(p)
    @compat finalizer(deleteModel, prob)
    return prob
end

function setParameter(prob::CbcModel, name::String, value::String)
    @assert isascii(name)
    @assert isascii(value)
    @cbc_ccall setParameter Cvoid (Ptr{Cvoid},Ptr{UInt8},Ptr{UInt8}) prob name value
end

# TODO: registerCallBack clearCallBack

function solve(prob::CbcModel)
    @cbc_ccall solve Cint (Ptr{Cvoid},) prob
end

@getproperty Float64 sumPrimalInfeasibilities
@getproperty Cint numberPrimalInfeasibilities

# TODO: checkSolution

@getproperty Cint getIterationCount

for s in (:isAbandoned, :isProvenOptimal, :isProvenInfeasible, :isContinuousUnbounded,
    :isNodeLimitReached, :isSecondsLimitReached, :isSolutionLimitReached, :isInitialSolveAbandoned, :isInitialSolveProvenOptimal, :isInitialSolveProvenPrimalInfeasible)
    @eval function ($s)(prob::CbcModel)
        check_problem(prob)
        v = @cbc_ccall $s Cint (Ptr{Cvoid},) prob
        return v != 0
    end
end

@getproperty Float64 getObjValue
@getproperty Float64 getBestPossibleObjValue
@getproperty Cint getNodeCount

# TODO: printSolution

function addSOS(prob::CbcModel, numRows::Integer, rowStarts::Vector{Cint},
    colIndices::Vector{Cint}, weights::Vector{Float64}, typ::Integer)

    @cbc_ccall(addSOS,Cvoid,(Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint},
                            Ptr{Float64}, Cint), prob, numRows,
                            rowStarts .- convert(Cint,1),
                            colIndices .- convert(Cint,1),
                            weights, typ)
end

# see Cbc_C_Interface.h documentation
@getproperty Cint status
@getproperty Cint secondaryStatus

# Both values are set to `-1` in the constructor and `resetModel`
function optimizeNotCalled(prob::CbcModel)
    return status(prob) == -1 && secondaryStatus(prob) == -1
end

end


text:
using Reduce

Reduce.@subtype FakeReal <: Real
load_package(:scope)

R(x) = FakeReal(RExpr(x))

macro R(x)
    :($(esc(x)) = FakeReal(RExpr($(QuoteNode(x)))))
end

export reduce

function reduce(m::JointDistribution,x)
    return _reduce(getmoduletypencoding(m.model), m.model, m.args, x)
end

@gg M function _reduce(_::Type{M}, _m::Model, _args, _data) where M <: TypeLevel{Module}
    Expr(:let,
        Expr(:(=), :M, from_type(M)),
        type2model(_m) |> sourceReduce() |> loadvals(_args, _data))
end

function sourceReduce()
    function(_m::Model)
        function proc(_m, st :: Assign)
            @q begin
                $(st.x) = RExpr($(st.rhs))
            end
        end

        function proc(_m, st :: Sample)
            @q begin
                $(st.x) = FakeReal(RExpr($(QuoteNode(st.x))))
                _ℓ += $(reduce())
            end
                 = :(_ℓ += reduce($(st.rhs), $(st.x)))
        proc(_m, st :: Return)     = nothing
        proc(_m, st :: LineNumber) = nothing

        wrap(kernel) = @q begin
            _ℓ = 0.0
            $kernel
            return _ℓ
        end

        buildSource(_m, proc, wrap) |> flatten
    end
end

# julia> Algebra.optimize(:(z = a^2*b^2+10*a^2*m^6+a^2*m^2+2*a*b*m^4+2*b^2*m^6+b^2*m^2))
# quote
#     g23 = b * a
#     g27 = m * m
#     g24 = g27 * b * b
#     g25 = g27 * a * a
#     g26 = g27 * g27
#     z = g24 + g25 + g23 * (2g26 + g23) + g26 * (2g24 + 10g25)
# end

# Algebra.sum(:(1/((p+(k-1)*q)*(p+k*q))),:k,1,:(n+1))

ext: jl


text:
__precompile__()

module Gems
using DataFrames
using CSV

gems = Gems
export gems

function load_gems9_10()
        file_path = join([@__DIR__, "gems9_10.csv"], "/")
        df = CSV.File(file_path) |> DataFrame
        df = rename(df, Dict(:source => :icd9, :target => :icd10))
    return df
end

function load_gems10_9()
        file_path = join([@__DIR__, "gems10_9.csv"], "/")
        df = CSV.File(file_path) |> DataFrame
        df = rename(df, Dict(:source => :icd10, :target => :icd9))
    return df
end

"""
    filter_flags(df, flag_type)

Limit GEM mapping table to a specific mapping relationship.
"""
function filter_flags(df::DataFrames.DataFrame,
        flag_type:: String;)
    if length(flag_type) > 0
        return df[df[Symbol(flag_type)] .== 1, :]
    end
end

function forward_mapping(icd_code:: String;
                flag_type:: String = "",
                hide_flags:: Bool = true)
        df = load_gems9_10()
        df = df[df[:icd9] .== icd_code,:]
        df = filter_flags(df, flag_type)
        if hide_flags
                df = df[:,[:icd9,:icd10,:target_descriptions,:source_descriptions]]
        end
    return df
end

function backward_mapping(icd_code:: String;
                flag_type:: String = "",
                hide_flags:: Bool = true)
        df = load_gems10_9()
        df = df[df[:icd10] .== icd_code,:]
        if hide_flags
                df = df[:,[:icd10,:icd9,:target_descriptions,:source_descriptions]]
        end
    return df
end

#function gems(icd_code:: String;
#        map_to:: String = "icd10",
#        flag_type:: String = "",
#        show_flags:: Bool = false)
#
#    df = load_lookups(map_to)
#
#    if length(flag_type) > 0
#        df = set_map_type(df,flag_type)
#    end
#
#    if show_flags
#        df = df[df[:source] .== icd_code, names(df)]
#        df = map_direction(df, map_to)
#    else
#        df = df[df[:source] .== icd_code,
#            [:source,:target,:descriptions ]]
#        df = map_direction(df, map_to)
#    end
#    return df
#end

end

max_stars_repo_name: pkmklong/Gems

max_issues_repo_name: pkmklong/Gems

max_forks_repo_name: pkmklong/Gems


text:
function HexaGrid(nrows, ncols)::Grid{Hexa,NonPeriodic}
    return Grid(Hexa, NonPeriodic, nrows, ncols)
end
function PeriodicHexaGrid(nrows, ncols)::Grid{Hexa,Periodic}
    @assert iseven(ncols)
    return Grid(Hexa, Periodic, nrows, ncols)
end


text:
using ITensorNetworkAD
using AutoHOOT, ITensors, Zygote
using ITensorNetworkAD.ITensorNetworks: GeneralMPSTensor

const itensorah = ITensorNetworkAD.ITensorAutoHOOT

@testset "test MPSTensor" begin
  i = Index(2, "i")
  j = Index(3, "j")
  k = Index(2, "k")
  l = Index(4, "l")
  m = Index(5, "m")

  A = randomITensor(i, j, k)
  B = randomITensor(k, l, m)
  C = randomITensor(i, j, l, m)
  mps_A = GeneralMPSTensor(MPS(A, inds(A)))
  mps_B = GeneralMPSTensor(MPS(B, inds(B)))
  mps_C = GeneralMPSTensor(MPS(C, inds(C)))

  out = A * B
  network = [mps_A, mps_B]
  nodes, dict = itensorah.generate_einsum_expr([network])
  out_list = itensorah.compute_graph(nodes, dict; cutoff=1e-15, maxdim=1000)
  @test isapprox(out, ITensor(out_list[1]))

  out = A * B * C
  out2 = contract(mps_A, mps_B, mps_C; cutoff=1e-15, maxdim=1000)
  @test isapprox(out, ITensor(out2))
end

@testset "test batch_tensor_contraction" begin
  i = Index(2, "i")
  j = Index(3, "j")
  k = Index(2, "k")
  A = randomITensor(i, j)
  B = randomITensor(j, k)
  C = randomITensor(k, i)

  function network(A)
    mps_A = GeneralMPSTensor(A; cutoff=1e-15, maxdim=1000)
    mps_B = GeneralMPSTensor(B; cutoff=1e-15, maxdim=1000)
    mps_C = GeneralMPSTensor(C; cutoff=1e-15, maxdim=1000)
    tensor_network = [mps_A, mps_B, mps_C]
    out = itensorah.batch_tensor_contraction(
      [tensor_network], mps_A; cutoff=1e-15, maxdim=1000
    )
    return sum(out)[]
  end
  grad_A = gradient(network, A)
  @test isapprox(grad_A[1], B * C)
end


text:
using SplitApplyCombine

@time groupreduce(x->x[1],x->x[2],(x,y) -> x+y, zip(id6,v1)); #51 seconds
@time sumby(id6,v1); # 48 seconds


text:
module TimeZoneLookup

using Shapefile
using Setfield
using AbstractTrees
import AbstractTrees: children, printnode
import Base: isless, zero, similar, show, insert!, delete!, pop!
using DataStructures
import DataStructures: insert_node!, RBTreeNode, fix_insert!, search_node, rb_transplant, minimum_node, delete_fix

export buildedges, mergeedges, buildedges!, buildmega, boundingbox, TrapezoidalSearchNode, TrapezoidData, query, Trapezoid, E, V, HE, Face
export faceof, EData, lhs, rhs, lp, rp, update!, sideof, buildsearch

export HE, Face, DCEL, twin, next, prev, loaddcel

include("structs.jl")
include("utils.jl")
include("pointlocator.jl")
include("dcel.jl")
include("sweep.jl")

end # module


text:
function _display_model(
    val::Dict{Tuple{String, Symbol}, T},
    field::Symbol,
    io::IO,
) where {T <: ServiceModel}
    field = titlecase(string(field))
    println(io, "$(field) Models:\n")
    if isempty(val)
        println("\t No Models Specified\n")
        return
    end
    for (i, ix) in val
        println(io, "\tType: $(ix.component_type)\n \tFormulation: $(ix.formulation)\n")
        if ix.use_service_name
            println(io, "\tName specific Model\n")
        end
    end
end

function _display_model(
    val::Dict{Symbol, T},
    field::Symbol,
    io::IO,
) where {T <: DeviceModel}
    field = titlecase(string(field))
    println(io, "$(field) Models: \n")
    if isempty(val)
        println("\t No Models Specified\n")
        return
    end
    for (i, ix) in val
        println(io, "\tType: $(ix.component_type)\n \tFormulation: $(ix.formulation)\n")
    end
end

"""
    Base.show(io::IO, ::MIME"text/plain", op_problem::OperationsProblem)

This function goes through the fields in OperationsProblem and then in OperationsProblemTemplate,
if the field contains a Device model dictionary, it calls organize_device_model() &
prints the data by field, key, value. If the field is not a Device model dictionary,
and a value exists for that field it prints the value.

"""
function Base.show(io::IO, m::MIME"text/plain", op_problem::OperationsProblem)
    show(io, m, op_problem.template)
end

function Base.show(io::IO, ::MIME"text/plain", template::OperationsProblemTemplate)
    println(io, "\nOperations Problem Specification")
    println(io, "============================================")

    for field in fieldnames(OperationsProblemTemplate)
        val = getfield(template, Symbol(field))
        if field == :transmission
            println(io, "Transmission: $val")
        elseif typeof(val) <: Dict{Symbol, <:DeviceModel}
            println(io, "============================================")
            _display_model(val, field, io)
        elseif typeof(val) <: Dict{Tuple{String, Symbol}, <:ServiceModel}
            println(io, "============================================")
            _display_model(val, field, io)
        else
            println(io, "")
        end
    end
    println(io, "============================================")
end

function Base.show(io::IO, optimization_container::OptimizationContainer)
    show(io, get_jump_model(optimization_container))
end

function Base.show(io::IO, sim::Simulation)
    println(io, "Simulation()")
end

function Base.show(io::IO, ::MIME"text/plain", results::SimulationResults)
    for res in values(results.problem_results)
        show(io, MIME"text/plain"(), res)
    end
end

function Base.show(io::IO, ::MIME"text/plain", results::SimulationProblemResults)
    title = results.problem * " Results"
    println(io, "\n$title")
    bars = join(("=" for _ in 1:length(title)))
    println(io, "$bars\n")
    timestamps = get_existing_timestamps(results)
    println(io, "Start: $(timestamps.start)")
    println(io, "End: $(timestamps.stop)")
    println(io, "Resolution: $(timestamps.step)")
    println(io, "\n")
    println(io, "Variables")
    println(io, "=========\n")
    for v in get_existing_variables(results)
        println(io, "$(v)")
    end
    println(io, "\n")
    parameters = get_existing_parameters(results)
    duals = get_existing_duals(results)
    for val in [("Parameters", parameters), ("Duals", duals)]
        if !isempty(val[2])
            println(io, "$(val[1])")
            println(io, "==========\n")
            for v in val[2]
                println(io, "$(v)")
            end
            println(io, "\n")
        end
    end
end
function Base.show(io::IO, ::MIME"text/html", results::PSIResults)
    println(io, "
Results
")
    timestamps = get_existing_timestamps(results)
    println(io, "
 Start: $(timestamps.start)
")
    println(io, "
 End: $(timestamps.stop)
")
    println(io, "
 Resolution: $(timestamps.step)
")
    println(io, "
Variables
")
    times = IS.get_timestamp(results)
    variables = IS.get_variables(results)
    if (length(keys(variables)) > 5)
        for (k, v) in variables
            println(io, "
$k: $(size(v))
")
        end
    else
        for (k, v) in IS.get_variables(results)
            if size(times, 1) == size(v, 1)
                var = hcat(times, v)
            else
                var = v
            end
            (l, w) = size(var)
            if w < 6
                println(io, "$(k)")
                println(io, "
$("-" ^ length("$k"))
")
                show(io, MIME"text/html"(), var)
            else
                println(io, "
$(k)  size ($l, $w)
")
            end
        end
    end
    parameters = IS.get_parameters(results)
    if !isempty(parameters)
        println(io, "
Parameters
")
        for (k, v) in parameters
            if size(times, 1) == size(v, 1)
                var = hcat(times, v)
            else
                var = v
            end
            (l, w) = size(var)
            if w < 6
                println(io, "$(k)")
                println(io, "
$("-" ^ length("$k"))
")
                show(io, MIME"text/html"(), var)
            else
                println(io, "
$(k)  size ($l, $w)
")
            end
        end
    end
    println(io, "
Optimizer Log
")
    for (k, v) in results.optimizer_stats
        if !(v === nothing)
            println(io, "
        $(k) = $(v)
")
        end
    end
    println(io, "\n")
    for (k, v) in results.total_cost
        println(io, "
Total Cost: $(v)
")
    end
end

function Base.show(io::IO, stage::OperationsProblem)
    println(io, "OperationsProblem()")
end

function Base.show(io::IO, ::MIME"text/plain", results::ProblemResults)
    println(io, "ProblemResults:")
    println(io, "  Base power: $(results.base_power)")
    vars = join(keys(results.variable_values), " ")
    println(io, "  Variables: $vars")
    duals = join(keys(results.dual_values), " ")
    println(io, "  Duals: $duals")
    params = join(keys(results.parameter_values), " ")
    println(io, "  Parameters: $params")
end

function Base.show(io::IO, ::MIME"text/html", services::Dict{Symbol, PSI.ServiceModel})
    println(io, "
Services
")
    for (k, v) in services
        println(io, "
$(k)
")
        println(io, "
$(v)
")
    end
end

function _count_stages(sequence::Array)
    stages = Dict{Int, Int}()
    stage = 1
    count = 0
    for i in 1:length(sequence)
        if sequence[i] == stage
            count += 1
            if i == length(sequence)
                stages[stage] = count
            end
        else
            stages[stage] = count
            stage += 1
            count = 1
        end
    end
    return stages
end

function _print_feedforward(io::IO, feed_forward::Dict, to::Array, from::Any)
    for (keys, sync) in feed_forward
        period = sync.periods
        stage1 = keys[1]
        stage2 = keys[2]
        spaces = " "^(length(stage2) + 2)
        dashes = "-"^(length(stage2) + 2)
        if period <= 12
            times = period
            line5 = string("└─$stage2 "^times, "($period) to : $to")
        else
            times = 12
            line5 = string("└─$stage2 "^times, "... (x$period) to : $to")
        end
        if times == 1
            line1 = "$stage1--┐ from : $from"
            println(io, "$line1\n$spaces|\n$spaces$line5\n")
        else
            if times == 2
                line3 = string("┌", string(dashes, "┤"))
                spacing = 0
            elseif times == 3
                line3 = string("┌", string(dashes, "┼"), string(dashes, "┐"))
                spacing = 0
            elseif iseven(times)
                spacing = (Int(times / 2) - 2)
                line3 = string(
                    "┌",
                    string(dashes, "┬")^spacing,
                    "----",
                    "┼",
                    string(dashes, "┬")^(spacing + 1),
                    "----┐",
                )
            else
                spacing = Int((times / 2) - 1.5)
                line3 = string(
                    "┌",
                    string(dashes, "┬")^spacing,
                    "----",
                    "┼",
                    string(dashes, "┬")^(spacing),
                    "----┐",
                )
            end
            line1 = string("     "^(spacing), " $stage1--┐ from : $from")
            line2 = string("     "^(spacing), " "^length(stage1), "   |")
            line4 = string("|", string(spaces, "|")^(times - 2), "    |")
            println(io, "$line1\n$line2\n$line3\n$line4\n$line4\n$line5\n")
        end
    end
end
function _print_inter_stages(io::IO, stages::Dict{Int, Int})
    list = sort!(collect(keys(stages)))
    for i in list
        num = stages[i]
        total = length(list)
        if total > 5
            println(io, "Too many stages to print.")
            break
        else
            if length("$num") == 1
                num = " (x0$num)"
            elseif length("$num") == 3
                num = " ($num)"
            elseif length("$num") == 4
                num = "($num)"
            else
                num = " (x$num)"
            end
            if i == 1
                if total == 2
                    println(io, "$i")
                else
                    println(io, "$i\n|")
                end
            else
                N = 2^(i - 2)
                space_count = 10 * (2^(total - i))
                if i == total
                    print1 = "|             ┌----/"^(N - 1)
                    print2 = "|             |     "^(N - 1)
                    print3 = "$i --> $i ...$num   "^N
                    println(io, "$print1|\n$print2|\n$print3")
                else
                    spaces = " "^(space_count - 11)
                    up = Int(space_count / 2)
                    print = string("$i ", " "^(space_count - 4), "  $i ...$num", spaces)^N
                    println(io, "$print")
                    if i !== total - 1
                        indent1 = string(
                            string("|", " "^(up + 7), "┌", "-"^(up - 10), "/")^(2 * N - 1),
                            "|",
                        )
                        indent2 = string("|", " "^(up + 7), "|", " "^(up - 9))^(2 * N - 1)
                        println(io, "$indent1\n$indent2|")
                    end
                end
            end
        end
    end
end

function _print_intra_stages(io::IO, stages::Dict{Int, Int})
    list = sort!(collect(keys(stages)))
    for i in list
        num = stages[i]
        total = length(list)
        if total > 5
            println(io, "Too many stages to print.")
            break
        else
            if length("$num") == 1
                num = " (x0$num)"
            elseif length("$num") == 3
                num = " ($num)"
            elseif length("$num") == 4
                num = "($num)"
            else
                num = " (x$num)"
            end
            if i == 1
                println(io, "$i\n")
            else
                N = 2^(i - 2)
                space_count = 10 * (2^(total - i))
                if i == total
                    print = "$i --> $i ...$num   "^N
                    println(io, "$print\n\n")
                else
                    spaces = " "^(space_count - 11)
                    print = string("$i ", "-"^(space_count - 4), "> $i ...$num", spaces)^N
                    indent = string(" "^(space_count))^(N * 2)
                    println(io, "$print\n$indent")
                end
            end
        end
    end
end

function Base.show(io::IO, sequence::SimulationSequence)
    stages = _count_stages(sequence.execution_order)
    println(io, "Feed Forward Chronology")
    println(io, "-----------------------\n")
    to = []
    from = String("")
    for (k, v) in sequence.feedforward
        println(io, "$(k[1]): $(typeof(v)) -> $(k[3])\n")
        to = String.(v.affected_variables)
        if isa(v, SemiContinuousFF)
            from = String.(v.binary_source_problem)
        elseif isa(v, RangeFF)
            from = String.([v.variable_source_problem_ub, v.variable_source_problem_lb])
        else
            from = String.(v.variable_source_problem)
        end
        _print_feedforward(io, sequence.feedforward_chronologies, to, from)
    end
    println(io, "Initial Condition Chronology")
    println(io, "----------------------------\n")
    if sequence.ini_cond_chronology == IntraProblemChronology()
        _print_intra_stages(io, stages)
    elseif sequence.ini_cond_chronology == InterProblemChronology()
        _print_inter_stages(io, stages)
    end
end


text:
using Parameters
using DifferentialEquations
using ForwardDiff

abstract type AbstractReactor end
export AbstractReactor

struct Reactor{D<:AbstractDomain} <: AbstractReactor
    domain::D
    ode::ODEProblem
end

function Reactor(domain::T,y0::Array{W,1},tspan::Tuple) where {T<:AbstractDomain,W<:Real}
    dydt(y::Array{T,1},p::V,t::Q) where {T<:Real,Q<:Real,V} = dydtreactor!(y,t,domain)
    ode = ODEProblem(dydt,y0,tspan)
    return Reactor(domain,ode)
end
export Reactor

@inline function getrate(rxn::T,cs::Array{W,1},kfs::Array{Q,1},krevs::Array{Q,1}) where {T<:AbstractReaction,Q,W<:Real}
    Nreact = length(rxn.reactantinds)
    Nprod = length(rxn.productinds)
    R = 0.0
    if Nreact == 1
        @fastmath @inbounds R += kfs[rxn.index]*cs[rxn.reactantinds[1]]
    elseif Nreact == 2
        @fastmath @inbounds R += kfs[rxn.index]*cs[rxn.reactantinds[1]]*cs[rxn.reactantinds[2]]
    elseif Nreact == 3
        @fastmath @inbounds R += kfs[rxn.index]*cs[rxn.reactantinds[1]]*cs[rxn.reactantinds[2]]*cs[rxn.reactantinds[3]]
    end

    if Nprod == 1
        @fastmath @inbounds R -= krevs[rxn.index]*cs[rxn.productinds[1]]
    elseif Nprod == 2
        @fastmath @inbounds R -= krevs[rxn.index]*cs[rxn.productinds[1]]*cs[rxn.productinds[2]]
    elseif Nprod == 3
        @fastmath @inbounds R -= krevs[rxn.index]*cs[rxn.productinds[1]]*cs[rxn.productinds[2]]*cs[rxn.productinds[3]]
    end

    return R
end
export getrate

@inline function addreactionratecontributions!(dydt::Array{Q,1},rarray::Array{UInt16,2},cs::Array{W,1},kfs::Array{Z,1},krevs::Array{Z,1}) where {Q<:Real,Z<:Real,T<:Integer,W<:Real}
    @inbounds @simd for i = 1:size(rarray)[2]
        if @inbounds rarray[2,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
        elseif @inbounds rarray[3,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
        else
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
        end
        if @inbounds rarray[5,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
        elseif @inbounds rarray[6,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
        else
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
        end
        @fastmath R = fR - rR
        @inbounds @fastmath dydt[rarray[1,i]] -= R
        if @inbounds rarray[2,i] != 0
            @inbounds @fastmath dydt[rarray[2,i]] -= R
            if @inbounds rarray[3,i] != 0
                @inbounds @fastmath dydt[rarray[3,i]] -= R
            end
        end
        @inbounds @fastmath dydt[rarray[4,i]] += R
        if @inbounds rarray[5,i] != 0
            @inbounds @fastmath dydt[rarray[5,i]] += R
            if @inbounds rarray[6,i] != 0
                @inbounds @fastmath dydt[rarray[6,i]] += R
            end
        end
    end
end
export addreactionratecontributions!

@inline function dydtreactor!(y::Array{U,1},t::Z,domain::Q;sensitivity::Bool=true) where {Z<:Real,U<:Real,J<:Integer,Q<:AbstractDomain}
    dydt = zeros(U,length(y))
    if sensitivity #if sensitivity isn't explicitly set to false set it to domain.sensitivity
        sensitivity = domain.sensitivity
    end
    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave = calcthermo(domain,y,t)
    addreactionratecontributions!(dydt,domain.rxnarray,cs,kfs,krevs)
    dydt *= V
    if sensitivity && isa(domain,ConstantVDomain)
        wV = copy(dydt[domain.indexes[1]:domain.indexes[2]])
    else
        wV = Array{U,1}()
    end
    calcdomainderivatives!(domain,dydt;t=t,T=T,P=P,Us=Us,V=V,C=C,ns=ns,N=N,Cvave=Cvave)
    if sensitivity
        Nspcs = length(cs)
        Nrxns = length(domain.phase.reactions)
        jacobian!(domain.jacobian,y,nothing,t,domain)
        dgdk = ratederivative(domain;cs=cs,V=V,T=T,kfs=kfs,Us=Us,Cvave=Cvave,N=N,krevs=krevs,wV=wV,sparse=false)
        for j  = 1:Nspcs+Nrxns #kfs and Gfs
            for i = 1:(domain.indexes[end]-domain.indexes[1]+1) #species
                for z in 1:(domain.indexes[end]-domain.indexes[1]+1)
                    dydt[j*Nspcs+i] += domain.jacobian[i,z]*y[j*Nspcs+z]
                end
                dydt[j*Nspcs+i] += dgdk[i,j]
            end
        end
    end
    return dydt
end
export dydtreactor!

function jacobian!(J::Q,y::U,p::W,t::Z,domain::V) where {Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
    if domain.t[1] == t && domain.jacuptodate == true
        return
    else
        f(y::Array{T,1}) where {T<:Real} = dydtreactor!(y,domain.t[1],domain;sensitivity=false)
        ForwardDiff.jacobian!(domain.jacobian,f,y[1:domain.indexes[end]-domain.indexes[1]+1])
        domain.jacuptodate[1] = true
        return
    end
end
export jacobian!

function ratederivative(d::W; cs::Q,V::Y,T::Y2,Us::Z3,Cvave::Z4,N::Z5,kfs::Z,krevs::X,wV::Q2,sparse::Bool=false) where {W<:Union{ConstantTPDomain,ConstantTVDomain},Z4<:Real,Z5<:Real,Z3<:AbstractArray,Q2<:AbstractArray,Q<:AbstractArray,Y2<:Real,Y<:Real,Z<:AbstractArray,X<:AbstractArray}
    Nspcs = length(cs)
    rxns = d.phase.reactions
    Nrxns = length(rxns)

    if sparse
        ratederiv = spzeros(Nspcs,Nspcs+Nrxns)
    else
        ratederiv = zeros(Nspcs,Nspcs+Nrxns)
    end
    RTinv = 1.0/(R*T)

    for (j,rxn) in enumerate(rxns)
        Nreact = length(rxn.reactantinds)
        Nprod = length(rxn.productinds)

        if Nreact == 1
            rind1 = rxn.reactantinds[1]
            fderiv = cs[rind1]
        elseif Nreact == 2
            rind1,rind2 = rxn.reactantinds
            fderiv = cs[rind1]*cs[rind2]
        else
            rind1,rind2,rind3 = rxn.reactantinds
            fderiv = cs[rind1]*cs[rind2]*cs[rind3]
        end

        if Nprod == 1
            pind1 = rxn.productinds[1]
            rderiv = krevs[j]/kfs[j]*cs[pind1]
        elseif Nprod == 2
            pind1,pind2 = rxn.productinds
            rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]
        else
            pind1,pind2,pind3 = rxn.productinds
            rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]*cs[pind3]
        end

        flux = fderiv-rderiv
        gderiv = rderiv*kfs[j]*RTinv

        deriv = zeros(Nspcs)

        deriv[rind1] += gderiv
        if Nreact > 1
            deriv[rind2] += gderiv
            if Nreact > 2
                deriv[rind3] == gderiv
            end
        end

        deriv[pind1] -= gderiv
        if Nprod > 1
            deriv[pind2] -= gderiv
            if Nprod > 2
                deriv[pind3] -= gderiv
            end
        end

        ratederiv[rind1,j] -= flux
        ratederiv[rind1,Nrxns+1:Nrxns+Nspcs] .-= deriv
        if Nreact > 1
            ratederiv[rind2,j] -= flux
            ratederiv[rind2,Nrxns+1:Nrxns+Nspcs] .-= deriv
            if Nreact > 2
                ratederiv[rind3,j] -= flux
                ratederiv[rind3,Nrxns+1:Nrxns+Nspcs] .-= deriv
            end
        end

        ratederiv[pind1,j] += flux
        ratederiv[pind1,Nrxns+1:Nrxns+Nspcs] .+= deriv
        if Nprod > 1
            ratederiv[pind2,j] += flux
            ratederiv[pind2,Nrxns+1:Nrxns+Nspcs] .+= deriv
            if Nprod > 2
                ratederiv[pind3,j] += flux
                ratederiv[pind3,Nrxns+1:Nrxns+Nspcs] .+= deriv
            end
        end
    end
    return V*ratederiv
end

function ratederivative(d::ConstantVDomain; cs::Q,V::Y,T::Y2,Us::Z3,Cvave::Y3,N::Y2,kfs::Z,krevs::X,wV::Q2,sparse::Bool=false) where {Z3<:AbstractArray,Q<:AbstractArray,Q2<:AbstractArray,Y3<:Real,Y2<:Real,Y<:Real,Z<:AbstractArray,X<:AbstractArray}
    Nspcs = length(cs)
    rxns = d.phase.reactions
    Nrxns = length(rxns)

    if sparse
        ratederiv = spzeros(Nspcs+1,Nspcs+Nrxns)
    else
        ratederiv = zeros(Nspcs+1,Nspcs+Nrxns)
    end
    RTinv = 1.0/(R*T)

    for (j,rxn) in enumerate(rxns)
        Nreact = length(rxn.reactantinds)
        Nprod = length(rxn.productinds)

        if Nreact == 1
            rind1 = rxn.reactantinds[1]
            fderiv = cs[rind1]
        elseif Nreact == 2
            rind1,rind2 = rxn.reactantinds
            fderiv = cs[rind1]*cs[rind2]
        else
            rind1,rind2,rind3 = rxn.reactantinds
            fderiv = cs[rind1]*cs[rind2]*cs[rind3]
        end

        if Nprod == 1
            pind1 = rxn.productinds[1]
            rderiv = krevs[j]/kfs[j]*cs[pind1]
        elseif Nprod == 2
            pind1,pind2 = rxn.productinds
            rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]
        else
            pind1,pind2,pind3 = rxn.productinds
            rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]*cs[pind3]
        end

        flux = fderiv-rderiv
        gderiv = rderiv*kfs[j]*RTinv

        deriv = zeros(Nspcs)

        deriv[rind1] += gderiv
        if Nreact > 1
            deriv[rind2] += gderiv
            if Nreact > 2
                deriv[rind3] == gderiv
            end
        end

        deriv[pind1] -= gderiv
        if Nprod > 1
            deriv[pind2] -= gderiv
            if Nprod > 2
                deriv[pind3] -= gderiv
            end
        end

        ratederiv[rind1,j] -= flux
        ratederiv[rind1,Nrxns+1:Nrxns+Nspcs] .-= deriv
        if Nreact > 1
            ratederiv[rind2,j] -= flux
            ratederiv[rind2,Nrxns+1:Nrxns+Nspcs] .-= deriv
            if Nreact > 2
                ratederiv[rind3,j] -= flux
                ratederiv[rind3,Nrxns+1:Nrxns+Nspcs] .-= deriv
            end
        end

        ratederiv[pind1,j] += flux
        ratederiv[pind1,Nrxns+1:Nrxns+Nspcs] .+= deriv
        if Nprod > 1
            ratederiv[pind2,j] += flux
            ratederiv[pind2,Nrxns+1:Nrxns+Nspcs] .+= deriv
            if Nprod > 2
                ratederiv[pind3,j] += flux
                ratederiv[pind3,Nrxns+1:Nrxns+Nspcs] .+= deriv
            end
        end
    end
    ratederiv *= V
    #Temperature stuff
    @views ratederiv[end,:] += (Us'*ratederiv[1:end-1,:])[1,:]
    ratederiv[end,1:Nspcs] += wV
    ratederiv[end,:] /= (N*Cvave)
    return ratederiv
end
export ratederivative


text:
using Distributions
using StatsBase
include("utils.jl")

# We assume here that proj is implemented symmetrically for primal and dual variables
function extraGrad(A::Array{Float64, 2}, proj::Function, z0::Array{Float64, 1},  τ::Float64, iter::Int64; tol=1e-6)
    m, n = size(A)
    @assert  size(z0)[1] == m + n
    energy = Float64[]
    x, y = z0[1:n], z0[(n+1):end]
    iter = ceil(Int, iter / 2) # since every iteration is 2 epochs.

    for i in 1:iter
        x, y, energy = extraGradUpdate(A, proj, x, y, energy, τ)
        if energy[end] < tol
            println("The extragradient algorithm achieves $tol accuracy in $(i-1) iterations")
            break
        end
    end
    return energy, x, y
end

"""
    extraGradUpdate: main update for extragradient method for matrix games on a simplex.
"""
function extraGradUpdate(A, proj, x, y, energy, τ)

    Ax, ATy =  A * x, A'* y
    gap = maximum(Ax) - minimum(ATy)
    append!(energy, gap)
    x_ = proj(x .- τ .* ATy)
    y_ = proj(y .+ τ .* Ax)
    # AA: added below to avoid multiplying by 2 when plotting
    Ax, ATy =  A * x_, A'* y_
    gap = maximum(Ax) - minimum(ATy)
    append!(energy, gap)
    x = proj(x .- τ .* ATy)
    y = proj(y .+ τ .* Ax)
    return x, y, energy
end

####################### Extragrad with Bregman update #######
function extraGrad_bregman(A::Array{T, 2}, z0::Array{T, 1},  τ::T, iter::Int64; tol=1e-6) where {T<:Real}
    m, n = size(A)
    iter = ceil(Int, iter / 2) # since every iteration is 2 epochs.
    if size(z0)[1] != m + n
        println("The dimension of a starting point doesn't match the ones of the matrix A")
    else
        energy = Float64[]
        x, y = z0[1:n], z0[(n+1):end]
        for i in 1:iter
            Ax, ATy =  A * x, A'* y
            gap = maximum(Ax) - minimum(ATy)
            append!(energy, gap)
            if gap > tol
                x_ = x .* exp.(-τ * ATy)
                x_ ./= sum(x_)
                y_ = y .* exp.(τ * Ax)
                y_ ./= sum(y_)
                Ax, ATy =  A * x_, A'* y_
                gap = maximum(Ax) - minimum(ATy)
                append!(energy, gap)
                x = x .* exp.(-τ .* ATy)
                x ./= sum(x)
                y  = y .* exp.(τ .* Ax)
                y ./= sum(y)
            else
                println("The extragradient-bregman algorithm achieves $tol accuracy in $(i-1) iterations")
                break
            end
        end
        return energy, x, y
    end
end

function extraGradBregman(A::Array{T, 2}, z0::Array{T, 1},
                          τ::T, iter::Int64; tol=1e-6) where {T<:Real}
    m, n = size(A)
    @assert m + n == size(z0)[1]
    energy = Float64[]
    x, y = z0[1:n], z0[(n+1):end]
    X, Y = zeros(Float64, n), zeros(Float64, m)
    for i in 1:iter
        x, y, X, Y, energy = extraGradBregmanUpdate(A, x, y, X, Y, energy, τ)
        if energy[end] < tol
            println("The extragradient algorithm achieves $tol accuracy in $(i-1) iterations")
            break
        end
    end
    return energy, x, y
end

function extraGradBregmanUpdate(A, x, y, X, Y, energy, τ)
    Ax, ATy =  A * x, A'* y
    gap = maximum(Ax) - minimum(ATy)
    append!(energy, gap)
    X_ = X - τ * ATy
    Y_ = Y + τ * Ax
    X = X - τ * (A' * softmax(Y_))
    Y = Y + τ * (A * softmax(X_))
    x = softmax(X)
    y = softmax(Y)
    return x, y, X, Y, energy
end

######################### ExtraGrad ######################################################
"""
    stochastic ExtraGrad with variance reduction, loopless variant.

bregman = true: uses mirror projection automatically. It ignores input function `proj`
bregman = false: uses Euclidean projection, given by `proj`
distr = false: uses uniform sampling for rows and columns of A
distr = true: uses weighted sampling. If `bregma=true` it uses l_1 sampling which is
    computed in every iteration. If `bregman=false` it uses a sampling |A[i,:]|^2/|A|_F and similarly for columns.
"""

function stochExtraGradLoopless(A::Array{Float64, 2}, proj::Function,
                           z0::Array{Float64, 1},  τ::Float64,
                           α::Float64, p::Float64,
                           max_epoch::Int64; bregman=false, distr=false, tol=1e-6)
    m, n = size(A)
    @assert m + n == size(z0)[1]
    iter = ceil(Int, max_epoch * 2 * m * n / (p * 2 * m * n + m + n ))
    cheap_update = (m + n) / (2 * m * n)
    update_wArray = rand(Bernoulli(p), iter)

    # check if sampling is uniform or with weights |A_i|^2/|A|_F
    if distr && !bregman
        arrayI, arrayJ, rows_weights, columns_weights = sample_with_Frobenius(A, iter)
    elseif !bregman
        arrayI = rand(1:m, iter)
        arrayJ = rand(1:n, iter)
    end
    x, y = z0[1:n], z0[(n+1):end]
    wx, wy = x, y
    wx_old, wy_old = x, y
    Awx, Awy = A * wx, A' * wy
    energy = [maximum(Awx) - minimum(Awy)]
    epoch = [0.0]
    epoch_count = 0.0

    for k in 1:iter
        update_w = update_wArray[k]
        # choose a distance for update: KL or euclidean
        if bregman
            x, y, wx, wy,  Awx, Awy, energy, epoch, epoch_count =
                stochExtraGradLooplessBregman_update(A, x, y, wx, wy, Awx, Awy, energy,
                                                     epoch, epoch_count, α, τ, cheap_update, update_w, distr)
        else
            i, j = arrayI[k], arrayJ[k]
            x, y, wx, wy, Awx, Awy, energy, epoch, epoch_count =
                stochExtraGradLoopless_update(A, proj, x, y, wx, wy, Awx, Awy,
                                         energy, epoch, epoch_count, α, τ, i, j, cheap_update,
                                         update_w, rows_weights, columns_weights)
        end
        if energy[end] < tol
            total_epoch = ceil(sum(epoch))
            println("StocExtraGrad-VR achieves $tol accuracy in $(total_epoch) epochs")
            break
        end
    end
    gap = maximum(A * x) - minimum(A' * y)
    append!(energy, gap)
    append!(epoch, 1.0)
    return energy, x, y, cumsum(epoch)
end

function stochExtraGradLoopless_update(A, proj, x, y, wx, wy, Awx, Awy, energy, epoch,
                                       epoch_count, α, τ, i, j, cheap_update, update_w,
                                       rows_weights, columns_weights)

    Ai = A[i, :]
    ATj = A[:, j]
    x_ = α .* x .+ (1-α) .* wx
    y_ = α .* y .+ (1-α) .* wy
    xx = proj(x_ .- τ .* Awy)
    yy = proj(y_ .+ τ .* Awx)
    x = proj(x_ .- τ .* (Awy  .+ (1/rows_weights[i] * (yy[i] - wy[i])) .* Ai))
    y = proj(y_ .+ τ .* (Awx  .+ (1/columns_weights[i] * (xx[j] - wx[j])) .* ATj))
    epoch_count += cheap_update
    if update_w
        wx, wy = x, y
        Awx, Awy = A * wx, A' * wy
        gap = maximum(Awx) - minimum(Awy)
        append!(energy, gap)
        epoch_count += 1.0
        append!(epoch, epoch_count)
        epoch_count = 0.0
    end

    return x, y, wx, wy,  Awx, Awy, energy, epoch, epoch_count
end

function stochExtraGradLooplessBregman_update(A, x, y, wx, wy, Awx, Awy, energy, epoch,
                                              epoch_count, α, τ, cheap_update, update_w, distr)
    # Main update
    x_ = x.^α .* wx.^(1-α)
    y_ = y.^α .* wy.^(1-α)
    xx = x_ .* exp.(-τ .*  Awy)
    yy = y_ .* exp.(τ .* Awx)
    xx ./= sum(xx)
    yy ./= sum(yy)

    i, weights_yi = sampling_Bregman(yy, wy, distr)
    j, weights_xj = sampling_Bregman(xx, wx, distr)
    Ai = A[i, :]
    ATj = A[:, j]

    x = x_ .* exp.(-τ .* (Awy  .+ ((yy[i] - wy[i]) / weights_yi) .* Ai))
    y = y_ .* exp.(τ .* (Awx  .+ ((xx[j] - wx[j]) / weights_xj) .* ATj))
    x ./= sum(x)
    y ./= sum(y)
    epoch_count += cheap_update
    if update_w
        wx, wy = x, y
        Awx, Awy = A * wx, A' * wy
        gap = maximum(Awx) - minimum(Awy)
        append!(energy, gap)
        epoch_count += 1
        append!(epoch, epoch_count)
        epoch_count = 0.0
    end
    return x, y, wx, wy, Awx, Awy, energy, epoch, epoch_count
end

##########################################################################################

"""
    stochastic FoRB with variance reduction, loopless variant.

bregman = true: uses mirror projection automatically. It ignores input function `proj`
bregman = false: uses Euclidean projection, given by `proj`
distr = false: uses uniform sampling for rows and columns of A
distr = true: uses weighted sampling. If `bregma=true` it uses l_1 sampling which is
    computed in every iteration. If `bregman=false` it uses a sampling |A[i,:]|^2/|A|_F and similarly for columns.
"""

function stochForbLoopless(A::Array{Float64, 2}, proj::Function,
                           z0::Array{Float64, 1},  τ::Float64,
                           α::Float64, p::Float64,
                           max_epoch::Int64; bregman=false, distr=false, tol=1e-6)
    m, n = size(A)
    @assert m + n == size(z0)[1]
    x_weights_ar, y_weights_ar = Float64[], Float64[]
    # Input is max_epoch now, let us compute approximate # of iterations
    # every iteration costs  p * 2mn + m+n and our overall budget is iter * 2mn
    # This will not be exact, but it seems close enough
    iter = ceil(Int, max_epoch * 2 * m * n / (p * 2 * m * n + m + n ))
    cheap_update = (m + n) / (2 * m * n)

    update_wArray = rand(Bernoulli(p), iter)
    # check if sampling is uniform or with weights |A_i|^2/|A|_F
    if distr && !bregman
        arrayI, arrayJ, rows_weights, columns_weights = sample_with_Frobenius(A, iter)
    elseif !bregman
        arrayI = rand(1:m, iter)
        arrayJ = rand(1:n, iter)
    end
    x, y = z0[1:n], z0[(n+1):end]
    wx, wy = x, y
    wx_old, wy_old = x, y
    Awx, Awy = A * wx, A' * wy
    energy = [maximum(Awx) - minimum(Awy)]
    epoch = [0.0]
    epoch_count = 0.0

    for k in 1:iter
        update_w = update_wArray[k]
        # choose a distance for update: KL or euclidean
        if bregman
            x, y, wx, wy, wx_old, wy_old, Awx, Awy, energy, epoch, epoch_count, x_weights_ar, y_weights_ar =
                stochForbLooplessBregman_update(A, x, y, wx, wy, wx_old, wy_old,
                                                Awx, Awy, energy, epoch, epoch_count,
                                                α, τ, cheap_update, update_w, distr,x_weights_ar, y_weights_ar)
        else
            i, j = arrayI[k], arrayJ[k]
            x, y, wx, wy, wx_old, wy_old, Awx, Awy, energy, epoch, epoch_count =
                stochForbLoopless_update(A, proj, x, y, wx, wy, wx_old, wy_old, Awx, Awy,
                                         energy, epoch, epoch_count, α, τ, i, j, cheap_update,
                                         update_w, rows_weights, columns_weights)
        end
        if energy[end] < tol
            total_epoch = ceil(sum(epoch))
            println("StochFoRB-VR algorithm achieves $tol accuracy in $(total_epoch) epochs")
            break
        end
    end
    gap = maximum(A * x) - minimum(A' * y)
    append!(energy, gap)
    append!(epoch, 1.0) ## not sure about this
    return energy, x, y, cumsum(epoch), x_weights_ar, y_weights_ar
end

function sample_with_Frobenius(A, iter)
    m, n = size(A)
    frobenius_norm = norm(A)
    rows_norm = [norm(A[i, :]) for i in 1:m]
    columns_norm = [norm(A[:, j]) for j in 1:n]
    rows_weights = rows_norm.^2 / frobenius_norm^2
    columns_weights = columns_norm.^2 / frobenius_norm^2
    arrayI = sample((1:m), Weights(rows_weights), iter)
    arrayJ = sample((1:n), Weights(columns_weights), iter)
    return arrayI, arrayJ, rows_weights, columns_weights
end

function stochForbLoopless_update(A, proj, x, y, wx, wy, wx_old, wy_old, Awx, Awy,
                                  energy, epoch, epoch_count, α, τ, i, j, cheap_update,
                                  update_w, rows_weights, columns_weights)
    Ai = A[i, :]
    ATj = A[:, j]
    x = proj(α .* x .+ (1-α) .* wx .- τ .* (Awy  .+ (1/rows_weights[i] * (y[i] - wy_old[i])) .* Ai))
    y = proj(α .* y .+ (1-α) .* wy .+ τ .* (Awx  .+ (1/columns_weights[i] * (x[j] - wx_old[j])) .* ATj))

    # Cost is (m+n) / (2mn). ONLY for dense case, I will do sparse case later.
    #cheap_update = (m + n) / (2 * m * n) I defined it in the main program
    epoch_count += cheap_update
    wx_old, wy_old = wx, wy
    if update_w
        wx, wy = x, y
        Awx, Awy = A * wx, A' * wy
        gap = maximum(Awx) - minimum(Awy)
        append!(energy, gap)
        epoch_count += 1.0
        append!(epoch, epoch_count)
        epoch_count = 0.0
    end
    return x, y, wx, wy, wx_old, wy_old, Awx, Awy, energy, epoch, epoch_count
end

function sampling_Bregman(x, wx, distr)
    # Sampling procedure
    n = length(x)
    # what is cheaper: to pass n as an argument or every time take length?
    if distr == true
        diff_x = abs.(x .- wx)
        sum_diff_x = sum(diff_x)
        if sum_diff_x < 1e-12
            j = 1 # anyway the vector x - wx is zero
            weights_xj = 1.0 / n
        else
            weights_x = diff_x ./ sum_diff_x
            j = sample((1:n), Weights(weights_x))
            weights_xj = weights_x[j]
        end
        @assert weights_xj != 0
    else # Uniform sampling
        j = sample((1:n))
        weights_xj = 1.0 / n
    end
    return j, weights_xj
end

function stochForbLooplessBregman_update(A, x, y, wx, wy, wx_old, wy_old,
                                         Awx, Awy, energy, epoch, epoch_count,
                                         α, τ, cheap_update, update_w, distr, x_weights_ar, y_weights_ar)

    i, weights_yi = sampling_Bregman(y, wy_old, distr)
    j, weights_xj = sampling_Bregman(x, wx_old, distr)
    append!(x_weights_ar, weights_xj)
    append!(y_weights_ar, weights_yi)

    # Main update
    Ai = A[i, :]
    ATj = A[:, j]
    xx = x.^α .* wx.^(1-α) .* exp.(-τ .* (Awy  .+ ((y[i] - wy_old[i]) / weights_yi) .* Ai))
    yy = y.^α .* wy.^(1-α) .* exp.(τ .* (Awx  .+ ((x[j] - wx_old[j]) / weights_xj) .* ATj))
    x = xx ./ sum(xx)
    y = yy ./ sum(yy)
    epoch_count +=  cheap_update
    wx_old, wy_old = wx, wy
    if update_w
        wx, wy = x, y
        Awx, Awy = A * wx, A' * wy
        gap = maximum(Awx) - minimum(Awy)
        append!(energy, gap)
        epoch_count += 1
        append!(epoch, epoch_count)
        epoch_count = 0.0
    end
    return x, y, wx, wy, wx_old, wy_old, Awx, Awy, energy, epoch, epoch_count, x_weights_ar, y_weights_ar
end

"""
Carmon et al. variant
"""
function stochMPCarmon(A::Array{Float64, 2}, proj::Function,
                       z0::Array{Float64, 1},  α::Float64,
                       η::Float64, max_epoch::Int64;
                       bregman=false, distr=false, tol=1e-6)

    m, n = size(A)
    @assert m + n == size(z0)[1]
    T = ceil(Int, 4 / (η * α))
    outer_iter = ceil(Int, max_epoch * 2 * m * n / (4 * m * n + T * (m + n)))
    cheap_update = (m + n) / (2 * m * n)
    # averaged cost for each time we save data. In [Carmon et al], every outer iteration
    # does two full updates, thus 2 times save data. Cost of one outer iteration
    # is (2 + K*cheap_update) so cost per save is (2 + K*cheap_update)/2
    cost_per_save = (2 + T * cheap_update)/2

    # check if sampling is uniform or with weights |A_i|^2/|A|_F
    if distr && !bregman
        arrayI, arrayJ, rows_weights, columns_weights = sample_with_Frobenius(A, T * outer_iter)
    elseif !bregman
        rows_weights = 1 / m * ones(Float64, n)
        columns_weights = 1 / n * ones(Float64, m)
        arrayI = rand(1:m, T * outer_iter) # AA: fixing error for iter
        arrayJ = rand(1:n, T * outer_iter)
    end
    x, y = z0[1:n], z0[(n+1):end]
    if bregman
        Awx, Awy = A * x, A' * y
        x0, y0 = ones(n) + log.(x), ones(m) + log.(y)
        x, y = x0, y0
    else
        x0, y0 = x, y
        Awx, Awy = A * x, A' * y
    end
    energy = [maximum(Awx) - minimum(Awy)]

    x_avg, y_avg = zeros(Float64, n), zeros(Float64, m)
    for k in 1:outer_iter
        if bregman
            Awx, Awy, energy = compute_full_operator(A, softmax(x0), softmax(y0), energy)
        else
            Awx, Awy, energy = compute_full_operator(A, x0, y0, energy)
        end
        if energy[end] < tol
            break
        end
        if bregman
            x, y, x_avg, y_avg =
                stochMPCarmon_update_bregman(A, x, y, x_avg, y_avg, x0, y0, Awx, Awy,
                                     α, η, m, n, T, k)
        else
            x, y, x_avg, y_avg =
                stochMPCarmon_update_euclidean(A, proj, x, y, x_avg, y_avg, x0, y0, Awx, Awy,
                                     α, η, arrayI, arrayJ, m, n, rows_weights, columns_weights, T, k)
        end
        # x_avg is in the primal space, so no need to softmax here.
        Awx, Awy, energy = compute_full_operator(A, x_avg, y_avg, energy)

        if energy[end] < tol
            #total_epoch = ceil(sum(epoch))
            println("Stochastic Carmon algorithm achieved $tol accuracy")
            break
        end

        if bregman
            # x0 in the dual space
            x0, y0 = stochMPCarmon_fullupdate_bregman(x0, y0, Awx, Awy, α)
        else
            x0, y0 = stochMPCarmon_fullupdate_euclidean(proj, x0, y0, Awx, Awy, α)
        end
    end
    running_cost = Array(1:length(energy)) * cost_per_save
    return energy, x, y, running_cost
end

function stochMPCarmon_update_euclidean(A, proj, x, y, x_avg, y_avg, x0, y0, Awx, Awy,
                              α, η, arrayI, arrayJ, m, n, rows_weights, columns_weights, T, k)
    x, y = x0, y0
    for t in 1:T
        i = arrayI[(k - 1) * T + t]
        j = arrayJ[(k - 1) * T + t]
        Ai = A[i, :]
        ATj = A[:, j]

        x = proj((x .+ (η * α / 2) .* x0 .- η .* (Awy .+ (1/rows_weights[i] * (y[i] - y0[i])) .* Ai)) ./ (1 + η*α/2) )
        y = proj((y .+ (η * α / 2) .* y0 .+ η .* (Awx .+ (1/columns_weights[j] * (x[j] - x0[j])) .* ATj)) ./ (1 + η*α/2) )
        x_avg = (1 / t) .* x .+ (1 - 1 / t) .* x_avg
        y_avg = (1 / t) .* y .+ (1 - 1 / t) .* y_avg
    end
    return x, y, x_avg, y_avg
end

function stochMPCarmon_update_bregman(A, x, y, x_avg, y_avg, x0, y0, Awx, Awy,
                                    α, η, m, n, T, k)
    x, y = x0, y0
    for t in 1:T
        oracle_x = MirProxVR_stoch_oracle(y, softmax(y0), A, 1)
        oracle_y = MirProxVR_stoch_oracle(x, softmax(x0), A, 2)

        x = (x .+ (η * α / 2) .* x0 .- η .* (Awy .+ oracle_x)) ./ (1 + η * α / 2)
        y = (y .+ (η * α / 2) .* y0 .+ η .* (Awx .+ oracle_y)) ./ (1 + η * α / 2)
        x_avg = (1 / t) .* softmax(x) .+ (1 - 1 / t) .* x_avg
        y_avg = (1 / t) .* softmax(y) .+ (1 - 1 / t) .* y_avg
    end
    return x, y, x_avg, y_avg
end

function stochMPCarmon_fullupdate_euclidean(proj, x0, y0, Awx, Awy, α)
    x = proj(x0 .- (Awy ./ α))
    y = proj(y0 .+ (Awx ./ α))
return x, y
end

function stochMPCarmon_fullupdate_bregman(x0, y0, Awx, Awy, α)
    x = x0 .- (Awy ./ α)
    y = y0 .+ (Awx ./ α)
return x, y
end

function compute_full_operator(A, x_avg, y_avg, energy)
    Awx, Awy = A * x_avg, A' * y_avg
    append!(energy, maximum(Awx) - minimum(Awy))
    return Awx, Awy, energy
end

#### Stochastic Bregman Extragradient-VR (mirror prox) for Matrix Games.

"""
Stochastic Extragradient with variance reduction, two-loop variant.
X = nabla h_1(x)
u, v = w_k
u_new, v_new = w_{k+1}

U, V = nabla h(w_bar_k)
U_new, V_new = nabla h(w_bar_{k+1})
"""
function MirProxVR_stoch_oracle(X, u, A, flag)
    N = length(X)
    x = softmax(X)
    diff = x .- u
    abs_diff = abs.(diff)
    norm_diff = sum(abs_diff)
    i = sample((1:N), Weights(abs_diff/norm_diff))
    if flag == 1
        row_or_column = A[i, :]
    else
        row_or_column = A[:, i]
    end
    return (norm_diff * sign(diff[i])) .* row_or_column
end

function MirProxVR_update(Z, w, W_, w_next, W_next, Fw, A, α, τ, k, n)
    Z_ = α .* Z .+ (1 - α) .* W_ .- τ .* Fw
    oracle_x = MirProxVR_stoch_oracle(Z_[(n+1):end], w[(n+1):end], A, 1)
    oracle_y = MirProxVR_stoch_oracle(Z_[1:n], w[1:n], A, 2)
    Z .= Z_ .- τ .* [oracle_x; -oracle_y]
    z = [softmax(Z[1:n]); softmax(Z[(n+1):end])]

    w_next = (z .+ (k-1) .* w_next) ./ k
    W_next = (Z .+ (k-1) .* W_next) ./ k
    return Z, w_next, W_next
end

function MirProxVR(A::Array{Float64, 2}, z0::Array{Float64, 1},
                        τ::Float64, α::Float64, K::Int64, max_epoch::Int64;  tol=1e-6)
    m, n = size(A)
    @assert m + n == size(z0)[1]
    cheap_update = (m + n) / (2 * m * n)
    cost_per_epoch = 1 + K * cheap_update
    S = ceil(Int, max_epoch * 2 * m * n / (2 * m * n + K * (m + n)))
    z = copy(z0)
    Z = log.(z) + ones(m+n)
    w = z
    W_ = Z
    w_next = zeros(Float64, m+n)
    W_next = zeros(Float64, m+n)
    Fw = [A' * w[(n+1):end];  -A * w[1:n]]
    energy = [maximum(-Fw[(n+1):end]) - minimum(Fw[1:n])]
    for s in 1:S
        for k in 1:K
            Z, w_next, W_next =
                MirProxVR_update(Z, w, W_, w_next, W_next, Fw, A, α, τ, k, n)
        end
        W_ = W_next
        w = w_next
        # no need to update w_next, W_next, since for k=1, they are just zero anyways
        Fw = [A' * w[(n+1):end];  -A * w[1:n]]
        gap = maximum(-Fw[(n+1):end]) - minimum(Fw[1:n])
        append!(energy, gap)
        if gap < tol
            println("StochMirProx-VR algorithm achieved $tol accuracy")
            break
        end
    end
    # compute how much we spent
    running_cost = Array(1:length(energy)) * cost_per_epoch
    z = softmax(Z)
    return energy, z[1:n], z[n:end], running_cost
end

######################################################
######################################################
######################################################
################# Extragrad - Looped #################
######################################################
######################################################
######################################################
function stochExtraGradLooped(A::Array{Float64, 2}, proj::Function,
                           z0::Array{Float64, 1},  τ::Float64,
                           α::Float64, T::Int64,
                           max_epoch::Int64; bregman=false, distr=false, tol=1e-6)
    m, n = size(A)
    @assert m + n == size(z0)[1]
    outer_iter = ceil(Int, max_epoch * 2 * m * n / (2 * m * n + T * (m + n)))
    cheap_update = (m + n) / (2 * m * n)
    cost_per_epoch = 1 + T * cheap_update

    if distr && !bregman
        arrayI, arrayJ, rows_weights, columns_weights = sample_with_Frobenius(A, T * outer_iter)
    elseif !bregman
        rows_weights = 1 / m * ones(Float64, n)
        columns_weights = 1 / n * ones(Float64, m)
        arrayI = rand(1:m, T * outer_iter)
        arrayJ = rand(1:n, T * outer_iter)
    end

    x, y = z0[1:n], z0[(n+1):end]
    x0, y0 = x, y
    wx, wy = x, y
    Awx, Awy = A * x, A' * y
    epoch = [0.0]
    energy = [maximum(Awx) - minimum(Awy)]
    epoch_count = 0.0

    x_avg, y_avg = zeros(Float64, n), zeros(Float64, m)

    for k in 1:outer_iter
        x, y, x_avg, y_avg =
            stochExtraGradLooped_update(A, proj, x, y, x_avg, y_avg, wx, wy, Awx, Awy,
                                 α, τ, arrayI, arrayJ, m, n, rows_weights, columns_weights, T, k)

        Awx, Awy, energy = compute_full_operator(A, x_avg, y_avg, energy)
        wx, wy = x_avg, y_avg
        if energy[end] < tol
            println("StocExtraGrad-VR-Looped achieved $tol accuracy")
            break
        end
    end
    # x_avg is on the primal space, so no need to softmax here.
    running_cost = Array(1:length(energy)) * cost_per_epoch
    return energy, x, y, running_cost
end

function stochExtraGradLooped_update(A, proj, x, y, x_avg, y_avg, wx, wy, Awx, Awy, α, τ,
                                        arrayI, arrayJ, m, n, rows_weights, columns_weights, T, k)
    for t in 1:T
        i = arrayI[(k - 1) * T + t]
        j = arrayJ[(k - 1) * T + t]
        Ai = A[i, :]
        ATj = A[:, j]
        x_ = α .* x .+ (1-α) .* wx
        y_ = α .* y .+ (1-α) .* wy
        xx = proj(x_ .- τ .* Awy)
        yy = proj(y_ .+ τ .* Awx)
        x = proj(x_ .- τ .* (Awy  .+ (1/rows_weights[i] * (yy[i] - wy[i])) .* Ai))
        y = proj(y_ .+ τ .* (Awx  .+ (1/columns_weights[i] * (xx[j] - wx[j])) .* ATj))

        x_avg = 1 / t .* x .+ (1 - 1 / t) .* x_avg
        y_avg = 1 / t .* y .+ (1 - 1 / t) .* y_avg
    end

    return x, y, x_avg, y_avg
end


text:
using DataFrames
using CSV

SAVE_PATH = "./data/raw_data/gems/"
URL_CMS = "https://www.cms.gov/Medicare/Coding/"

URL_GEMS = string(URL_CMS, "ICD10/Downloads/2018-ICD-10-CM-General-Equivalence-Mappings.zip")
URL_ICD10  = string(URL_CMS, "ICD10/Downloads/2018-ICD-10-Code-Descriptions.zip")
URL_ICD9 =  string(URL_CMS, "ICD9ProviderDiagnosticCodes/Downloads/ICD-9-CM-v32-master-descriptions.zip")

URL_GEMS_PCS = string(URL_CMS, "ICD10/Downloads/ProcedureGEMs-2014.zip")
URL_GEMS_PCS_DESC = string(URL_CMS, "ICD10/Downloads/2014-PCS-long-and-abbreviated-titles.zip")

cd("../data/")

gems_files = download(URL_GEMS,"gems.zip")
icd10_files = download(URL_ICD10,"icd10.zip")
icd9_files = download(URL_ICD9, "icd9.zip")
gems_pcs = download(URL_GEMS_PCS, "gems_pcs.zip")
pcs_desc = download(URL_GEMS_PCS_DESC, "pcs_desc.zip")

run(`unzip gems.zip`);
run(`unzip icd10.zip`);
run(`unzip icd9.zip`);
run(`unzip gems_pcs.zip`);
run(`unzip pcs_desc`);


text:
"""
    authorize()

Get and store client credentials. Any other credentials will be dropped.
"""
function authorize()
    @info "Retrieving client credentials, which typically lasts 1 hour."
    SPOTCRED[] = get_spotify_credentials()
    if string(SPOTCRED[]) != string(SpotifyCredentials())
        @info "Expires at $(SPOTCRED[].expires_at). Access `spotcred()`, refresh with `refresh_spotify_credentials()`."
    else
        @info "When configured, `refresh_spotify_credentials()`."
    end
end
function refresh_spotify_credentials()
    SPOTCRED[] = get_spotify_credentials()
    if string(SPOTCRED[]) != string(SpotifyCredentials())
        @info "Expires at $(SPOTCRED[].expires_at)."
    else
        @info "When configured, `refresh_spotify_credentials()`."
    end
end

function get_spotify_credentials()
    c = get_init_file_spotify_credentials()
    if string(c) == string(SpotifyCredentials())
        return c
    end
    j = get_authorization_token(c)
    if isnothing(j)
        return SpotifyCredentials()
    end
    c.access_token = j.access_token
    c.token_type = j.token_type
    c.expires_at = string(Dates.now() + Dates.Second(j.expires_in))
    c
end

function get_authorization_token(sc_tokenless::SpotifyCredentials)
    refreshtoken = sc_tokenless.encoded_credentials
    headers = ["Authorization" => "Basic $refreshtoken",
              "Accept" => "*/*",
              "Content-Type" => "application/x-www-form-urlencoded"]
    body = "grant_type=client_credentials"
    resp = HTTP.Messages.Response()
    try
        resp = HTTP.post(AUTH_URL, headers, body)
    catch e
        request = "HTTP.post call: AUTH_URL = $AUTH_URL\n  headers = $headers \n  body = $body"
        @error request
        @error e
        return nothing
    end
    response_body = resp.body |> String
    response_body |> JSON3.read
end

function get_init_file_spotify_credentials()
    id, secret, redirect = get_id_secret_redirect()
    if id == NOT_ACTUAL || secret == NOT_ACTUAL
        @warn "User needs to configure $(_get_ini_fnam())"
        SpotifyCredentials()
    else
        enc_cred = base64encode(id * ":" * secret)
        c = SpotifyCredentials(client_id = id, client_secret =secret, encoded_credentials = enc_cred)
        c.redirect = redirect
        c
    end
end

"Get id and secret as 32-byte string, no encryption"
function get_id_secret_redirect()
    container = read(Inifile(), _get_ini_fnam())
    id = get(container, "Spotify developer's credentials", "CLIENT_ID",  NOT_ACTUAL)
    secret = get(container, "Spotify developer's credentials", "CLIENT_SECRET",  NOT_ACTUAL)
    redirect = get(container, "Spotify developer's credentials", "REDIRECT_URI",  DEFAULT_REDIRECT_URI)
    id, secret, redirect
end
#=
function get_working_browser_cmd()
    container = read(Inifile(), _get_ini_fnam())
    cmds = get(container, "Working browser command (auto filled in)", "CMD",  "")
    @cmd cmds
end
=#
#=
function set_working_browser_cmd(cmd)
    container = read(Inifile(), _get_ini_fnam())
    set(container, "Working browser command (auto filled in)", "CMD",  string(cmd))
end
=#
function get_user_name()
    container = read(Inifile(), _get_ini_fnam())
    get(container, "Spotify user id", "user_name",  "")
end

"Get an existing, readable ini file name, create it if necessary"
function _get_ini_fnam()
    fna = joinpath(homedir(), "spotify_credentials.ini")
    if !isfile(fna)
        open(fna, "w+") do io
            _prepare_init_file_with_instructions(io)
        end
        if Sys.iswindows()
            run(`cmd /c $fna`; wait = false)
        end
        println("Instructions in $fna")
    end
    fna
end

function _prepare_init_file_with_instructions(io)
    conta = Inifile()
    set(conta, "User procedure:", "1: How to get the CLIENT_ID", "Developer.spotify.com -> Dashboard -> log in -> Create an app -> Fill out name and purpose -> Copy to the field below.")
    set(conta, "User procedure:", "2: How to get CLIENT_SECRET", "When you have client id, press 'Show client secret' -> Copy to below. No need for quotation marks.")
    set(conta, "User procedure:", "3: How to give REDIRECT_URL", "Still in the app dashboard:
    'Edit settings' -> Redirect uris -> http://127.0.0.1:8080 -> Save changes")
    set(conta, "Spotify developer's credentials", "CLIENT_ID", NOT_ACTUAL)
    set(conta, "Spotify developer's credentials", "CLIENT_SECRET", NOT_ACTUAL)
    set(conta, "Spotify developer's credentials", "REDIRECT_URI", "http://127.0.0.1:8080")
    set(conta, "Spotify user id", "user_name", "Slartibartfast")
    write(io, conta)
end

client_credentials_still_valid() = now() < parse(DateTime, spotcred().expires_at)


text:
module RaviartThomasRefFEsTest

using Test
using Gridap.Polynomials
using Gridap.Fields
using Gridap.TensorValues
using Gridap.Fields: MockField
using Gridap.ReferenceFEs

p = QUAD
D = num_dims(QUAD)
et = Float64
order = 0

reffe = RaviartThomasRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 4
@test get_order(get_prebasis(reffe)) == 0
@test num_dofs(reffe) == 4
@test get_default_conformity(reffe) == DivConformity()

p = QUAD
D = num_dims(QUAD)
et = Float64
order = 1

reffe = RaviartThomasRefFE(et,p,order)
test_reference_fe(reffe)
@test num_terms(get_prebasis(reffe)) == 12
@test num_dofs(reffe) == 12
@test get_order(get_prebasis(reffe)) == 1

prebasis = get_prebasis(reffe)
dof_basis = get_dof_basis(reffe)

v = VectorValue(3.0,0.0)
field = MockField{D}(v)

cache = dof_cache(dof_basis,field)
r = evaluate_dof!(cache, dof_basis, field)
test_dof(dof_basis,field,r)

cache = dof_cache(dof_basis,prebasis)
r = evaluate_dof!(cache, dof_basis, prebasis)
test_dof(dof_basis,prebasis,r)

end # module


text:
using DataFrames

"""
This function takes in a file path/name and exports a dataframe consisting of the circuit parameters
Specifically the file should be from LTSPICE (or it wont work)

If no input is given, the user will be prompted for a FileName
"""
function LTSpiceLoad(FileName=nothing)
        # if FileName===nothing
        #         FileName = open_dialog("Pick a file")
        # end
        f =open(FileName)
         Lines =  readlines(f)

         close(f)

         SPICE_DF = DataFrame((Type = Any[], Value=Any[], Node1=Any[], Node2=Any[], Name=Any[],ESR = Any[]))
        for i in 1:length(Lines)

                Spaces = vcat(findall(" ",Lines[i])...)

                if ~(Lines[i][1]=='*')
                        # println(Lines[i])
                end
                if (Lines[i][1]=='R')|(Lines[i][1]=='L')|(Lines[i][1]=='C')
                        if length(Spaces)==3
                                Value = Lines[i][Spaces[3]+1:end]
                        else
                                Value = Lines[i][Spaces[3]+1:Spaces[4]-1]
                        end
                        Node1 = Lines[i][Spaces[1]+1:Spaces[2]-1]
                        Node2 = Lines[i][Spaces[2]+1:Spaces[3]-1]
                        Name = Lines[i][1:Spaces[1]-1]
                        RSerInds = findall("Rser=", Lines[i])
                        if length(RSerInds)>0
                                RSerIndEnd =RSerInds[1][end]+1
                                if RSerInds[1][end]>maximum(Spaces)
                                        RSerVal = MakeNumericalVals(Lines[i][RSerIndEnd:end])
                                else
                                        Tmp = findfirst(x-> x>RSerIndEnd,Spaces)
                                        RSerVal = MakeNumericalVals(Lines[i][RSerIndEnd:Spaces[Tmp]])
                                end
                        else
                                RSerVal=0
                        end
                        push!(SPICE_DF,[Lines[i][1] MakeNumericalVals(Value) Node1 Node2 Name RSerVal])

                elseif (Lines[i][1]=='V')
                        Node1 = Lines[i][Spaces[1]+1:Spaces[2]-1]
                        Node2 = Lines[i][Spaces[2]+1:Spaces[3]-1]
                        Name = Lines[i][1:Spaces[1]-1]
                        push!(SPICE_DF,[Lines[i][1] 1 Node1 Node2 Name 0])

                end

        end

        return SPICE_DF
end

"""
A simple string to numeric conversion for metric prefixes
Input: String
Output Float64

e.g.
        MakeNumericalVals("1n")
         > 1.0e-9
        MakeNumericalVals.(["1n" , "2m"])
                2-element Array{Float64,1}:
                1.0e-9
                0.002
"""
function MakeNumericalVals(ValString::String)
        NewString = replace(ValString,"Meg"=>"e6")
        NewString = replace(NewString,"k"=>"e3")
        NewString = replace(NewString,"m"=>"e-3")
        NewString = replace(NewString,"u"=>"e-6")
        NewString = replace(NewString, "n"=>"e-9")
        NewString = replace(NewString, "p"=>"e-12")
        return parse(Float64,NewString)
end


text:
# form intervals
function endPoints( A, V, M )
    nvec = size(V,2);
    p = kpm(A, V, M);

    ∂p = Polynomials.derivative(p);
    ∂2p = Polynomials.derivative(∂p);

    r_∂p = Polynomials.roots(∂p) |>
        x -> filter( isreal, x ) |>
        real                     |>
        x-> filter(
            item -> -1 < item < 1,
            x )

    r_∂2p = Polynomials.roots(∂2p) |>
        x -> filter( isreal, x )   |>
        real                       |>
        x-> filter(
            item -> -1 < item < 1,
            x )

    roots =  [ r_∂p; r_∂2p]
    roots,p
end


text:
@testset "Elastic net ($T)" for T in [Float32, Float64, ComplexF32, ComplexF64]

    using ProximalOperators
    using ProximalAlgorithms
    using LinearAlgebra
    using Random

    Random.seed!(0)

    # TODO: generate problem with known solution instead

    m, n = 50, 200

    A = randn(T, m, n)
    b = randn(T, m)

    R = real(T)

    reg = ElasticNet(R(1), R(1))
    reg1 = NormL1(R(1))
    reg2 = SqrNormL2(R(1))
    loss = Translate(SqrNormL2(R(1)), -b)

    L = opnorm(A)^2

    x0 = zeros(T, n)
    solver = ProximalAlgorithms.PANOC{R}(tol=eps(R))
    x_panoc, it = solver(x0, f=loss, A=A, g=reg)

    @testset "DavisYin" begin

        x0 = zeros(T, n)
        solver = ProximalAlgorithms.DavisYin{R}(tol=R(1e-6))
        xf_dys, xg_dys, it_dys = solver(
            x0, f=reg1, g=reg2, h=loss, A=A, L=L
        )
        @test eltype(xf_dys) == T
        @test eltype(xg_dys) == T
        @test norm(xf_dys - x_panoc, Inf) <= 1e-3
        @test norm(xg_dys - x_panoc, Inf) <= 1e-3
        @test it_dys <= 1900

        x0 = randn(T, n)
        solver = ProximalAlgorithms.DavisYin{R}(tol=R(1e-6))
        xf_dys, xg_dys, it_dys = solver(
            x0, f=reg1, g=reg2, h=loss, A=A, L=L
        )
        @test eltype(xf_dys) == T
        @test eltype(xg_dys) == T
        @test norm(xf_dys - x_panoc, Inf) <= 1e-3
        @test norm(xg_dys - x_panoc, Inf) <= 1e-3
        @test it_dys <= 2600

    end

    afba_test_params = [
        (R(2), R(0), 230),
        (R(1), R(1), 230),
        (R(0), R(1), 470),
        (R(0), R(0), 500),
        (R(1), R(0), 230)
    ]

    @testset "AFBA" for (theta, mu, maxit) in afba_test_params

        x0 = randn(T, n)
        y0 = randn(T, m)

        solver = ProximalAlgorithms.AFBA{R}(theta=theta, mu=mu, tol=R(1e-6))
        x_afba, y_afba, it_afba = solver(
            x0, y0, f=reg2, g=reg1, h=loss, L=A, betaQ=R(1),
        )
        @test eltype(x_afba) == T
        @test eltype(y_afba) == T
        @test norm(x_afba - x_panoc, Inf) <= 1e-4
        @test it_afba <= maxit

    end
end


text:
## Exercise 2-3
## Practice using the Julia REPL as a calculator:

## 1. The volume of a sphere with radius r is (4/3)πr^3. What is the volume of a sphere with radius 5?
println("Ans 1:")
println("the volume of sphere with radius 5 is: " * string((4/3) * π * 5 ^ 3))

## 2. Suppose the cover price of a book is $ 24.95, but bookstores get a 40 % discount. Shipping costs $ 3 for the first copy and 75 cents for each additional copy. What is the total wholesale cost for 60 copies?
println("Ans 2:")
cover_price = 24.95
n_copies = 60
shipping_cost = 3 + (n_copies - 1) * 0.75
books_price = n_copies * cover_price * 0.6
total = shipping_cost + books_price
println(total)

## 3. If I leave my house at 6:52 am and run 1 mile at an easy pace (8:15 per mile), then 3 miles at tempo (7:12 per mile) and 1 mile at easy pace again, what time do I get home for breakfast?
println("Ans 3:")
easy_pace_time = 8*60 + 15
tempo_pace_time = 7*60 + 12

total_time = 2 * easy_pace_time + 3 * tempo_pace_time

seconds = total_time % 60
minutes = total_time ÷ 60

println("you will reach home at 6:52 + ", minutes, " minutes and ", seconds, " seconds")
println("End.")


text:
# Multivariate LogNormal distribution

###########################################################
#
#   Abstract base class for multivariate lognormal
#
#   Each subtype should provide the following methods:
#
#   - length(d)         vector dimension
#   - params(d)         Get the parameters from the underlying Normal distribution
#   - location(d)       Location parameter
#   - scale(d)          Scale parameter
#   - _rand!(d, x)      Sample random vector(s)
#   - _logpdf(d,x)      Evaluate logarithm of pdf
#   - _pdf(d,x)         Evaluate the pdf
#   - mean(d)           Mean of the distribution
#   - median(d)         Median of the distribution
#   - mode(d)           Mode of the distribution
#   - var(d)            Vector of element-wise variance
#   - cov(d)            Covariance matrix
#   - entropy(d)        Compute the entropy
#
#
###########################################################

@compat abstract type AbstractMvLogNormal <: ContinuousMultivariateDistribution end

function insupport{T<:Real,D<:AbstractMvLogNormal}(::Type{D},x::AbstractVector{T})
    for i=1:length(x)
      @inbounds 0.0<x[i]<Inf?continue:(return false)
    end
    true
end
insupport{T<:Real}(l::AbstractMvLogNormal,x::AbstractVector{T}) = insupport(typeof(l),x)
assertinsupport{D<:AbstractMvLogNormal}(::Type{D},m::AbstractVector) = @assert insupport(D,m) "Mean of LogNormal distribution should be strictly positive"

###Internal functions to calculate scale and location for a desired average and covariance
function _location!{D<:AbstractMvLogNormal}(::Type{D},::Type{Val{:meancov}},mn::AbstractVector,S::AbstractMatrix,μ::AbstractVector)
    @simd for i=1:length(mn)
      @inbounds μ[i] = log(mn[i]/sqrt(1+S[i,i]/mn[i]/mn[i]))
    end
    μ
end

function _scale!{D<:AbstractMvLogNormal}(::Type{D},::Type{Val{:meancov}},mn::AbstractVector,S::AbstractMatrix,Σ::AbstractMatrix)
    for j=1:length(mn)
      @simd for i=j:length(mn)
        @inbounds Σ[i,j] = Σ[j,i] = log(1 + S[j,i]/mn[i]/mn[j])
      end
    end
    Σ
end

function _location!{D<:AbstractMvLogNormal}(::Type{D},::Type{Val{:mean}},mn::AbstractVector,S::AbstractMatrix,μ::AbstractVector)
    @simd for i=1:length(mn)
      @inbounds μ[i] = log(mn[i]) - S[i,i]/2
    end
    μ
end

function _location!{D<:AbstractMvLogNormal}(::Type{D},::Type{Val{:median}},md::AbstractVector,S::AbstractMatrix,μ::AbstractVector)
    @simd for i=1:length(md)
      @inbounds μ[i] = log(md[i])
    end
    μ
end

function _location!{D<:AbstractMvLogNormal}(::Type{D},::Type{Val{:mode}},mo::AbstractVector,S::AbstractMatrix,μ::AbstractVector)
    @simd for i=1:length(mo)
      @inbounds μ[i] = log(mo[i]) + S[i,i]
    end
    μ
end

###Functions to calculate location and scale for a distribution with desired :mean, :median or :mode and covariance
function location!{D<:AbstractMvLogNormal}(::Type{D},s::Symbol,m::AbstractVector,S::AbstractMatrix,μ::AbstractVector)
  @assert size(S) == (length(m),length(m)) && length(m) == length(μ)
  assertinsupport(D,m)
  _location!(D,Val{s},m,S,μ)
end

function location{D<:AbstractMvLogNormal}(::Type{D},s::Symbol,m::AbstractVector,S::AbstractMatrix)
    @assert size(S) == (length(m),length(m))
    assertinsupport(D,m)
    _location!(D,Val{s},m,S,similar(m))
end

function scale!{D<:AbstractMvLogNormal}(::Type{D},s::Symbol,m::AbstractVector,S::AbstractMatrix,Σ::AbstractMatrix)
    @assert size(S) == size(Σ) == (length(m),length(m))
    assertinsupport(D,m)
    _scale!(D,Val{s},m,S,Σ)
end

function scale{D<:AbstractMvLogNormal}(::Type{D},s::Symbol,m::AbstractVector,S::AbstractMatrix)
    @assert size(S) == (length(m),length(m))
    assertinsupport(D,m)
    _scale!(D,Val{s},m,S,similar(S))
end

params!{D<:AbstractMvLogNormal}(::Type{D},m::AbstractVector,S::AbstractMatrix,μ::AbstractVector,Σ::AbstractMatrix) = location!(D,:meancov,m,S,μ),scale!(D,:meancov,m,S,Σ)
params{D<:AbstractMvLogNormal}(::Type{D},m::AbstractVector,S::AbstractMatrix) = params!(D,m,S,similar(m),similar(S))

#########################################################
#
#   MvLogNormal
#
#   Multivariate lognormal distribution based on MvNormal
#
#########################################################
immutable MvLogNormal{T<:Real,Cov<:AbstractPDMat,Mean<:Union{Vector, ZeroVector}} <: AbstractMvLogNormal
    normal::MvNormal{T,Cov,Mean}
end

#Constructors mirror the ones for MvNormmal
MvLogNormal(μ::Union{Vector,ZeroVector},Σ::AbstractPDMat) = MvLogNormal(MvNormal(μ,Σ))
MvLogNormal(Σ::AbstractPDMat) = MvLogNormal(MvNormal(ZeroVector(eltype(Σ),dim(Σ)),Σ))
MvLogNormal(μ::Vector,Σ::Matrix) = MvLogNormal(MvNormal(μ,Σ))
MvLogNormal(μ::Vector,σ::Vector) = MvLogNormal(MvNormal(μ,σ))
MvLogNormal(μ::Vector,s::Real) = MvLogNormal(MvNormal(μ,s))
MvLogNormal(Σ::Matrix) = MvLogNormal(MvNormal(Σ))
MvLogNormal(σ::Vector) = MvLogNormal(MvNormal(σ))
MvLogNormal(d::Int,s::Real) = MvLogNormal(MvNormal(d,s))

### Conversion
function convert{T<:Real}(::Type{MvLogNormal{T}}, d::MvLogNormal)
    MvLogNormal(convert(MvNormal{T}, d.normal))
end
function convert{T<:Real}(::Type{MvLogNormal{T}}, pars...)
    MvLogNormal(convert(MvNormal{T}, MvNormal(pars...)))
end

length(d::MvLogNormal) = length(d.normal)
params(d::MvLogNormal) = params(d.normal)
@inline partype{T<:Real}(d::MvLogNormal{T}) = T
location(d::MvLogNormal) = mean(d.normal)
scale(d::MvLogNormal) = cov(d.normal)

#See https://en.wikipedia.org/wiki/Log-normal_distribution
mean(d::MvLogNormal) = @compat(exp.(mean(d.normal) + var(d.normal)/2))
median(d::MvLogNormal) = @compat(exp.(mean(d.normal)))
mode(d::MvLogNormal) = @compat(exp.(mean(d.normal) - var(d.normal)))
function cov(d::MvLogNormal)
    m = mean(d)
    return m*m'.*(@compat(exp.(cov(d.normal))) - 1)
  end
var(d::MvLogNormal) = diag(cov(d))

#see Zografos & Nadarajah (2005) Stat. Prob. Let 71(1) pp71-84 DOI: 10.1016/j.spl.2004.10.023
entropy(d::MvLogNormal) = length(d)*(1+log2π)/2 + logdetcov(d.normal)/2 + sum(mean(d.normal))

#See https://en.wikipedia.org/wiki/Log-normal_distribution
_rand!{T<:Real}(d::MvLogNormal, x::AbstractVecOrMat{T}) = @compat(exp!(_rand!(d.normal, x)))
@compat _logpdf{T<:Real}(d::MvLogNormal, x::AbstractVecOrMat{T}) = insupport(d, x) ? (_logpdf(d.normal, log.(x)) - sum(log.(x))) : -Inf
_pdf{T<:Real}(d::MvLogNormal, x::AbstractVecOrMat{T}) = insupport(d,x) ? _pdf(d.normal, @compat(log.(x)))/prod(x) : 0.0

Base.show(io::IO,d::MvLogNormal) = show_multline(io, d, [(:dim, length(d)), (:μ, mean(d.normal)), (:Σ, cov(d.normal))])


text:
module ReinforcementLearningCore

using ReinforcementLearningBase
using Reexport

const RLCore = ReinforcementLearningCore

export RLCore

@reexport using ReinforcementLearningTrajectories

include("core/core.jl")
include("policies/policies.jl")
include("utils/utils.jl")

end # module


text:
export Grid
export dx, dy, dh, xc, yc, xe, ye
export coord2ind, x2ind, y2ind

const DEFAULT_L₀ = 1e-6;

struct Grid{D}
    L::SVector{D,Float} # Length of dimensions
    L₀::Float # Length unit
    N::SVector{D,Int64} # Number of cells
    Npml::SVector{D,Int64} # Number of PML cells on each end for each dimension
    bounds::Tuple2{SVector{D,Float}} # Boundary coordinates, accessed as bounds[edge][dim]
end

Base.ndims(::Grid{D}) where {D} = D
Base.length(g::Grid) = prod(g.N)
Base.size(g::Grid{2}) = (g.N[1], g.N[2])
Base.size(g::Grid{1}) = (g.N[1])
Base.size(g::Grid{D}, i::Integer) where {D} = i > D ? 1 : g.N[i]

"""
    Grid(dh::Number, Npml::AbstractArray{<:Integer}, xrange::AbstractArray{<:Real}, yrange::AbstractArray{<:Real}; L₀=DEFAULT_L₀)

Constructs a 2D Grid.
"""
function Grid(dh::Number, Npml::AbstractArray{<:Integer}, xrange::AbstractArray{<:Real}, yrange::AbstractArray{<:Real}; L₀=DEFAULT_L₀)
    L = SVector{2}([ Float(xrange[2]-xrange[1]), Float(yrange[2]-yrange[1])]);
    N = SVector{2}(Int.(round.(L/dh)));
    bounds1 = SVector{2}([Float(xrange[1]), Float(yrange[1])]);
    bounds2 = SVector{2}([Float(xrange[2]), Float(yrange[2])]);
    @debug "Grid size: $N"
    return Grid(L, L₀, N, SVector{2}(Npml), (bounds1, bounds2));
end

"""
    Grid(dh::Number, Npml::Integer, xrange::AbstractArray{<:Real}; L₀=DEFAULT_L₀)

Constructs a 1D Grid.
"""
function Grid(dh::Number, Npml::Integer, xrange::AbstractArray{<:Real}; L₀=DEFAULT_L₀)
    L = SVector{1}(Float(xrange[2]-xrange[1]));
    N = SVector{1}(Int.(round.(L/dh)));
    bounds1 = SVector{1}(Float(xrange[1]));
    bounds2 = SVector{1}(Float(xrange[2]));
    return Grid(L, L₀, N, SVector{1}(Npml), (bounds1, bounds2));
end

"""
    Grid(N::Integer, xrange::AbstractArray{<:Real}; L₀=DEFAULT_L₀)

Constructs a 1D Grid.
"""
function Grid(N::Integer, xrange::AbstractArray{<:Real}; L₀=DEFAULT_L₀)
    L = SVector{1}(Float(xrange[2]-xrange[1]));
    N = SVector{1}(N);
    bounds1 = SVector{1}(Float(xrange[1]));
    bounds2 = SVector{1}(Float(xrange[2]));
    return Grid(L, L₀, SVector{1}(N), SVector{1}(0), (bounds1, bounds2));
end

"Returns the cell size for the grid `g` in the direction specified by `w`"
function dh(g::Grid, w::Direction)
	w == x̂ && return dx(g)
    w == ŷ && return dy(g)
end

function dx(g::Grid)
	return (g.bounds[2][1]-g.bounds[1][1])/g.N[1]
end

function dy(g::Grid{2})
	return (g.bounds[2][2]-g.bounds[1][2])/g.N[2]
end

function xc(g::Grid)
	return g.bounds[1][1] .+ dx(g)*(0.5:1:g.N[1])
end

function yc(g::Grid{2})
	return g.bounds[1][2] .+ dy(g)*(0.5:1:g.N[2])
end

xc(g::Grid, i::CartesianIndex{2}) = xc(g)[i.I[1]];
yc(g::Grid, i::CartesianIndex{2}) = yc(g)[i.I[2]];
xc(g::Grid, i::Int64) = xc(g)[ CartesianIndices(size(g))[i][1] ];
yc(g::Grid, i::Int64) = yc(g)[ CartesianIndices(size(g))[i][2] ];

function xe(g::Grid)
	return g.bounds[1][1] .+ dx(g)*(0:1:g.N[1])
end

function ye(g::Grid{2})
	return g.bounds[1][2] .+ dy(g)*(0:1:g.N[2])
end

xe(g::Grid, i::CartesianIndex{2}) = xe(g)[i.I[1]];
ye(g::Grid, i::CartesianIndex{2}) = ye(g)[i.I[2]];
xe(g::Grid, i::Int64) = xe(g)[ CartesianIndices(size(g))[i][1] ];
ye(g::Grid, i::Int64) = ye(g)[ CartesianIndices(size(g))[i][2] ];

"    coord2ind(g::Grid{D}, pt::Point)"
function coord2ind(g::Grid{D}, pt::Point) where {D}
    indx = x2ind(g, pt.x);
    if D == 1
    	return indx
    end
    return (indx, y2ind(g, pt.y))
end

"    x2ind(g::Grid, x::Real)"
function x2ind(g::Grid, x::Real)
    ind = Int(round((x-g.bounds[1][1])/g.L[1]*size(g,1))+1);
    ind < 1 && return 1
    ind > g.N[1] && return g.N[1]
    return ind
end

"    y2ind(g::Grid, y::Real)"
function y2ind(g::Grid, y::Real)
    ind = Int(round((y-g.bounds[1][2])/g.L[2]*size(g,2))+1);
    ind < 1 && return 1
    ind > g.N[2] && return g.N[2]
    return ind
end

"    δ(w::Direction, s::DerivativeDirection, g::Grid)"
function δ(w::Direction, s::DerivativeDirection, g::Grid{K}) where {K}
    if K == 2
        (Nx, Ny) = size(g);
    else
        Nx = size(g, 1);
        Ny = 1;
    end

    if w == x̂
        if s == Forward
            δxf = 1/dx(g)*spdiagm(1 => ones(Nx-1), 0 => -ones(Nx), -Nx+1 => [1]);
            return kron(sparse(1.0I, Ny, Ny), δxf)
        else
            δxb = 1/dx(g)*spdiagm(-1 => -ones(Nx-1), 0 => ones(Nx), Nx-1 => [-1]);
            return kron(sparse(1.0I, Ny, Ny), δxb)
        end
    end
    if w == ŷ
        if s == Forward
            δyf = 1/dy(g)*spdiagm(1 => ones(Ny-1), 0 => -ones(Ny), -Ny+1 => [1]);
            return kron(δyf, sparse(1.0I, Nx, Nx))
        else
            δyb = 1/dy(g)*spdiagm(-1 => -ones(Ny-1), 0 => ones(Ny), Ny-1 => [-1]);
            return kron(δyb, sparse(1.0I, Nx, Nx))
        end
    end
end

"    grid_average(centerarray::AbstractArray, w::Direction)"
function grid_average(centerarray::AbstractArray, w::Direction)
    ndims(centerarray) == 1 && return (centerarray+circshift(centerarray, (1)))/2
    w == x̂ && return (centerarray+circshift(centerarray, (1, 0)))/2;
    w == ŷ && return (centerarray+circshift(centerarray, (0, 1)))/2;
    return centerarray
end


text:
using FiniteVolumes, OrdinaryDiffEq, Plots

mesh = PeriodicCartesianMesh(100)
u0 = [sin(2π*x) for x in cell_centers(mesh)]
F(u, n) = u^2 / 2 * n

dudt(u, t, p) = -FiniteVolumes.div(F, mesh, u, Upwind())
problem = ODEProblem(dudt, u0, (0.0, 0.5))
sol = solve(problem, Euler(), dt=1e-3, saveat=0.05)

plot(mesh, [u0 sol(0.25) sol(0.5)],
     label=["initial condition" "middle state" "final state"])


text:
# Define block-arrays
module BlockArrays

using LinearAlgebra

export RealOrComplex,
       BlockArray,
       blocksize,
       blockeltype,
       blocklength,
       blockvecnorm,
       blockmaxabs,
       blocksimilar,
       blockcopy,
       blockcopy!,
       blockset!,
       blockvecdot,
       blockzeros,
       blockones,
       blockaxpy!,
       blockscale!,
       blockcumscale!,
       blockiszero

const RealOrComplex{R} = Union{R, Complex{R}}
const BlockArray{R} = Union{
	AbstractArray{C, N} where {C <: RealOrComplex{R}, N},
	Tuple{Vararg{AbstractArray{C, N} where {C <: RealOrComplex{R}, N}}}
}

# Operations on block-arrays

blocksize(x::Tuple) = blocksize.(x)
blocksize(x::AbstractArray) = size(x)

blockeltype(x::Tuple) = blockeltype.(x)
blockeltype(x::AbstractArray) = eltype(x)

blocklength(x::Tuple) = sum(blocklength.(x))
blocklength(x::AbstractArray) = length(x)

blockvecnorm(x::Tuple) = sqrt(real(blockvecdot(x, x)))
blockvecnorm(x::AbstractArray{R}) where {R <: Number} = norm(x)

blockmaxabs(x::Tuple) = maximum(blockmaxabs.(x))
blockmaxabs(x::AbstractArray{R}) where {R <: Number}= maximum(abs, x)

blocksimilar(x::Tuple) = blocksimilar.(x)
blocksimilar(x::AbstractArray) = similar(x)

blockcopy(x::Tuple) = blockcopy.(x)
blockcopy(x::AbstractArray) = copy(x)

blockcopy!(y::Tuple, x::Tuple) = blockcopy!.(y, x)
blockcopy!(y::AbstractArray, x::AbstractArray) = copyto!(y, x)

blockset!(y::Tuple, x) = blockset!.(y, x)
blockset!(y::AbstractArray, x) = (y .= x)

blockvecdot(x::T1, y::T2) where {T1 <: Tuple, T2 <: Tuple} = sum(blockvecdot.(x,y))
blockvecdot(x::AbstractArray{R1}, y::AbstractArray{R2}) where {R1 <: Number, R2 <: Number} = dot(x, y)

blockzeros(t::Tuple, s::Tuple) = blockzeros.(t, s)
blockzeros(t::Type, n::NTuple{N, Integer} where {N}) = zeros(t, n)
blockzeros(t::Tuple) = blockzeros.(t)
blockzeros(n::NTuple{N, Integer} where {N}) = zeros(n)
blockzeros(n::Integer) = zeros(n)
blockzeros(a::AbstractArray) = fill!(similar(a),zero(eltype(a)))

blockones(t::Tuple, s::Tuple) = blockones.(t, s)
blockones(t::Type, n::NTuple{N, Integer} where {N}) = ones(t, n)
blockones(t::Tuple) = blockones.(t)
blockones(n::NTuple{N, Integer} where {N}) = ones(n)
blockones(n::Integer) = ones(n)
blockones(a::AbstractArray) = fill!(a,one(eltype(a)))

blockscale!(z::Tuple, alpha::Real, y::Tuple) = blockscale!.(z, alpha, y)
blockscale!(z::AbstractArray, alpha::Real, y::AbstractArray) = (z .= alpha.*y)

blockcumscale!(z::Tuple, alpha::Real, y::Tuple) = blockcumscale!.(z, alpha, y)
blockcumscale!(z::AbstractArray, alpha::Real, y::AbstractArray) = (z .+= alpha.*y)

blockaxpy!(z::Tuple, x::Tuple, alpha::Real, y::Tuple) = blockaxpy!.(z, x, alpha, y)
blockaxpy!(z::AbstractArray, x::AbstractArray, alpha::Real, y::AbstractArray) = (z .= x .+ alpha.*y)

blockiszero(x::AbstractArray) = iszero(x)
blockiszero(x::Tuple) = all(iszero.(x))

# Define broadcast

import Base: broadcast!

function broadcast!(f::Any, dest::Tuple, op1::Tuple)
   for k = eachindex(dest)
       broadcast!(f, dest[k], op1[k])
   end
   return dest
end

function broadcast!(f::Any, dest::Tuple, op1::Tuple, op2::Tuple)
   for k = eachindex(dest)
       broadcast!(f, dest[k], op1[k], op2[k])
   end
   return dest
end

function broadcast!(f::Any, dest::Tuple, coef::Number, op2::Tuple)
   for k = eachindex(dest)
       broadcast!(f, dest[k], coef, op2[k])
   end
   return dest
end

function broadcast!(f::Any, dest::Tuple, op1::Tuple, coef::Number, op2::Tuple)
   for k = eachindex(dest)
       broadcast!(f, dest[k], op1[k], coef, op2[k])
   end
   return dest
end

function broadcast!(f::Any, dest::Tuple, op1::Tuple, coef::Number, op2::Tuple, op3::Tuple)
   for k = eachindex(dest)
       broadcast!(f, dest[k], op1[k], coef, op2[k], op3[k])
   end
   return dest
end

end


text:
#
# Factors.jl Tests
#

using Factors
using Base.Test
using DataFrames

include("test_dimensions.jl")

include("test_factors.jl")
include("test_factors_broadcast.jl")
include("test_factors_index.jl")
include("test_factors_iter.jl")
include("test_factors_join.jl")
include("test_factors_map.jl")
include("test_factors_reduce.jl")

# include("test_dataframes.jl") # tested elsewhere ...
include("test_negatives.jl")


text:
"""
    CircShiftedArray(parent::AbstractArray, shifts)

Custom `AbstractArray` object to store an `AbstractArray` `parent` circularly shifted by `shifts` steps (where `shifts` is
a `Tuple` with one `shift` value per dimension of `parent`). Note that `shift` is modified with a modulo operation and does
not store the passed value but instead a positive number which leads to an equivalent shift.
Use `copy` to collect the values of a `CircShiftedArray` into a normal `Array`.

# Examples

```jldoctest circshiftedarray
julia> v = [1, 3, 5, 4];

julia> s = CircShiftedArray(v, (1,))
4-element CircShiftedVector{Int64, Vector{Int64}}:
 4
 1
 3
 5

julia> copy(s)
4-element Vector{Int64}:
 4
 1
 3
 5
```
"""
struct CircShiftedArray{T, N, S<:AbstractArray} <: AbstractArray{T, N}
    parent::S
    # the field `shifts` stores the circular shifts modulo the size of the parent array
    shifts::NTuple{N, Int}
    function CircShiftedArray(p::AbstractArray{T, N}, n = ()) where {T, N}
        shifts = map(mod, padded_tuple(p, n), size(p))
        return new{T, N, typeof(p)}(p, shifts)
    end
end

"""
    CircShiftedVector{T, S<:AbstractArray}

Shorthand for `CircShiftedArray{T, 1, S}`.
"""
const CircShiftedVector{T, S<:AbstractArray} = CircShiftedArray{T, 1, S}

CircShiftedVector(v::AbstractVector, n = ()) = CircShiftedArray(v, n)

size(s::CircShiftedArray) = size(parent(s))
axes(s::CircShiftedArray) = axes(parent(s))

@inline function bringwithin(ind_with_offset::Int, ranges::AbstractUnitRange)
    return ifelse(ind_with_offset < first(ranges), ind_with_offset + length(ranges), ind_with_offset)
end

@inline function getindex(s::CircShiftedArray{T, N}, x::Vararg{Int, N}) where {T, N}
    @boundscheck checkbounds(s, x...)
    v, ind = parent(s), offset(shifts(s), x)
    i = map(bringwithin, ind, axes(s))
    return @inbounds v[i...]
end

@inline function setindex!(s::CircShiftedArray{T, N}, el, x::Vararg{Int, N}) where {T, N}
    @boundscheck checkbounds(s, x...)
    v, ind = parent(s), offset(shifts(s), x)
    i = map(bringwithin, ind, axes(s))
    @inbounds v[i...] = el
    return s
end

parent(s::CircShiftedArray) = s.parent

"""
    shifts(s::CircShiftedArray)

Return amount by which `s` is shifted compared to `parent(s)`.
"""
shifts(s::CircShiftedArray) = s.shifts