text:
# Copyright 2015, 2016 Martin Holters
# See accompanying license file.

VERSION >= v"0.4.0-dev+6521" && __precompile__()

module ACME

export Circuit, add!, connect!, DiscreteModel, run!, steadystate, steadystate!

using ProgressMeter
using Compat

import Base.getindex

include("kdtree.jl")
include("solvers.jl")

type Element
  mv :: SparseMatrixCSC{Real,Int}
  mi :: SparseMatrixCSC{Real,Int}
  mx :: SparseMatrixCSC{Real,Int}
  mxd :: SparseMatrixCSC{Real,Int}
  mq :: SparseMatrixCSC{Real,Int}
  mu :: SparseMatrixCSC{Real,Int}
  u0 :: SparseMatrixCSC{Real,Int}
  pv :: SparseMatrixCSC{Real,Int}
  pi :: SparseMatrixCSC{Real,Int}
  px :: SparseMatrixCSC{Real,Int}
  pxd :: SparseMatrixCSC{Real,Int}
  pq :: SparseMatrixCSC{Real,Int}
  nonlinear_eq :: Expr
  pins :: Dict{Symbol, Vector{@compat Tuple{Int, Int}}}

  function Element(;args...)
    sizes = @compat Dict{Symbol,Int}(:n0 => 1)

    function update_sizes(mat, syms)
      for (sym, s) in zip(syms, size(mat))
        if !haskey(sizes, sym)
          sizes[sym] = s
        elseif sizes[sym] ≠ s
          error("Inconsistent sizes for ", sym)
        end
      end
    end

    function make_pin_dict(syms)
      dict = @compat Dict{Symbol,Vector{@compat Tuple{Int, Int}}}()
      for i in 1:length(syms)
        branch = div(i+1, 2)
        polarity = 2mod(i, 2) - 1
        push!(get!(dict, @compat(Symbol(syms[i])), []), (branch, polarity))
      end
      dict
    end
    make_pin_dict(dict::Dict) = dict

    const mat_dims =
        @compat  Dict( :mv => (:nl,:nb), :mi => (:nl,:nb), :mx => (:nl,:nx),
                       :mxd => (:nl,:nx), :mq => (:nl,:nq), :mu => (:nl,:nu),
                       :u0 => (:nl, :n0),
                       :pv => (:ny,:nb), :pi => (:ny,:nb), :px => (:ny,:nx),
                       :pxd => (:ny,:nx), :pq => (:ny,:nq) )

    elem = new()
    for (key, val) in args
      if haskey(mat_dims, key)
        val = convert(SparseMatrixCSC{Real}, sparse(hcat(val))) # turn val into a sparse matrix whatever it is
        update_sizes(val, mat_dims[key])
      elseif key == :pins
        val = make_pin_dict(val)
      end
      setfield!(elem, key, val)
    end
    for (m, ns) in mat_dims
      if !isdefined(elem, m)
        setfield!(elem, m, spzeros(Real, get(sizes, ns[1], 0), get(sizes, ns[2], 0)))
      end
    end
    if !isdefined(elem, :nonlinear_eq)
      elem.nonlinear_eq = Expr(:block)
    end
    if !isdefined(elem, :pins)
      elem.pins = make_pin_dict(1:2nb(elem))
    end
    elem
  end
end

for (n,m) in @compat Dict(:nb => :mv, :nx => :mx, :nq => :mq, :nu => :mu)
  @eval ($n)(e::Element) = size(e.$m, 2)
end
nl(e::Element) = size(e.mv, 1)
ny(e::Element) = size(e.pv, 1)
nn(e::Element) = nb(e) + nx(e) + nq(e) - nl(e)

# a Pin combines an element with a branch/polarity list
typealias Pin @compat Tuple{Element, Vector{Tuple{Int,Int}}}

# allow elem[:pin] notation to get an elements pin
getindex(e::Element, p) = (e, e.pins[@compat Symbol(p)])

include("elements.jl")

typealias Net Vector{@compat Tuple{Int,Int}} # each net is a list of branch/polarity pairs

type Circuit
    elements :: Vector{Element}
    nets :: Vector{Net}
    net_names :: Dict{Symbol, Net}
    Circuit() = new([], [], Dict{Symbol, Net}())
end

for n in [:nb; :nx; :nq; :nu; :nl; :ny; :nn]
    @eval ($n)(c::Circuit) = sum([$n(elem) for elem in c.elements])
end

for mat in [:mv; :mi; :mx; :mxd; :mq; :mu; :pv; :pi; :px; :pxd; :pq]
    # blkdiag() does not work, so include an empty matrix of desired type in
    # case c.elements is empty
    @eval ($mat)(c::Circuit) = blkdiag(spzeros(Real, 0, 0),
                                       [elem.$mat for elem in c.elements]...)
end

u0(c::Circuit) = vcat([elem.u0 for elem in c.elements]...)

function incidence(c::Circuit)
    i = sizehint!(Int[], 2nb(c))
    j = sizehint!(Int[], 2nb(c))
    v = sizehint!(Int[], 2nb(c))
    for (row, pins) in enumerate(c.nets), (branch, polarity) in pins
        push!(i, row)
        push!(j, branch)
        push!(v, polarity)
    end
    # ensure zeros due to short-circuited branches are removed, hence the
    # additional sparse(findnz(...))
    sparse(findnz(sparse(i,j,v))..., length(c.nets), nb(c))
end

function nonlinear_eq(c::Circuit)
    # construct a block expression containing all element's expressions after
    # offsetting their indexes into q, J and res

    row_offset = 0
    col_offset = 0
    nl_expr = Expr(:block)
    for elem in c.elements
        index_offsets = @compat Dict( :q => (col_offset,),
                                      :J => (row_offset, col_offset),
                                      :res => (row_offset,) )

        function offset_indexes(expr::Expr)
            ret = Expr(expr.head)
            ret.typ = expr.typ
            if expr.head == :ref && haskey(index_offsets, expr.args[1])
                push!(ret.args, expr.args[1])
                offsets = index_offsets[expr.args[1]]
                length(expr.args) == length(offsets) + 1 ||
                    throw(ArgumentError(string(expr.args[1], " must be indexed with exactly ", length(offsets), " index(es)")))
                for i in 1:length(offsets)
                    push!(ret.args,
                          :($(offsets[i]) + $(offset_indexes(expr.args[i+1]))))
                end
            else
                append!(ret.args, map(offset_indexes, expr.args))
            end
            ret
        end

        function offset_indexes(s::Symbol)
            haskey(index_offsets, s) && throw(ArgumentError(string(s, " used without indexing expression")))
            s
        end

        offset_indexes(x::Any) = x

        # wrap into a let to keep variables local
        push!(nl_expr.args, :( let; $(offset_indexes(elem.nonlinear_eq)) end))

        row_offset += nn(elem)
        col_offset += nq(elem)
    end
    nl_expr
end

function add!(c::Circuit, elem::Element)
    elem ∈ c.elements && return
    b_offset = nb(c)
    push!(c.elements, elem)
    for branch_pols in values(elem.pins)
        push!(c.nets, [(b_offset + b, pol) for (b, pol) in branch_pols])
    end
    nothing
end

add!(c::Circuit, elems::Element...) = for elem in elems add!(c, elem) end

function branch_offset(c::Circuit, elem::Element)
    offset = 0
    for el in c.elements
        el == elem && return offset
        offset += nb(el)
    end
    throw(ArgumentError("Element not found in circuit"))
end

function netfor!(c::Circuit, p::Pin)
    element = p[1]
    add!(c, element)
    b_offset = branch_offset(c, element)
    local net
    for (branch, pol) in p[2], net in c.nets
        (branch + b_offset, pol) ∈ net && break
    end
    @assert isdefined(net)
    net
end

function netfor!(c::Circuit, name::Symbol)
    haskey(c.net_names, name) || push!(c.nets, get!(c.net_names, name, []))
    c.net_names[name]
end

function connect!(c::Circuit, pins::(@compat Union{Pin,Symbol})...)
    nets = unique([netfor!(c, pin) for pin in pins])
    for net in nets[2:end]
        append!(nets[1], net)
        deleteat!(c.nets, findfirst(c.nets, net))
        for (name, named_net) in c.net_names
            if named_net == net
                c.net_names[name] = nets[1]
            end
        end
    end
end

# lines marked with !SV avoid creation of SparseVector by indexing with Ranges
# instead of Ints; a better way for cross-julia-version compatibilty would be
# nice; maybe Compat helps in the future...
function topomat!{T<:Integer}(incidence::SparseMatrixCSC{T})
    @assert all(abs(nonzeros(incidence)) .== 1)
    @assert all(sum(incidence, 1) .== 0)

    t = falses(size(incidence)[2]);

    row = 1;
    for col = 1:size(incidence)[2]
        rows = filter(r -> r ≥ row, find(incidence[:,col:col])) # !SV
        @assert length(rows) ≤ 2

        isempty(rows) && continue
        t[col] = true;

        if rows[1] ≠ row
            swaprows!(incidence, rows[1], row)
        end
        if length(rows) == 2
            @assert incidence[row, col] + incidence[rows[2], col] == 0
            incidence[rows[2],:] = incidence[rows[2],:] + incidence[row,:]
        end
        if incidence[row, col] < 0
            cols = find(incidence[row:row,:]) # !SV
            incidence[row,cols] = -incidence[row,cols]
        end
        rows = find(incidence[1:row-1,col:col] .== 1) # !SV
        incidence[rows,:] = broadcast(-, incidence[rows, :], incidence[row:row,:]) # !SV
        rows = find(incidence[1:row-1,col:col] .== -1) # !SV
        incidence[rows,:] = broadcast(+, incidence[rows, :], incidence[row:row,:]) # !SV
        row += 1
    end

    if row > 1
        ti = incidence[1:row-1, :]
    else
        ti = spzeros(T, 0, size(incidence)[2])
    end

    if all(t)
        dl = spzeros(T, row-1, 0)
        tv = spzeros(T, 0, size(incidence)[2])
    else
        dl = ti[:, ~t]
        tv = spzeros(T, size(dl)[2], size(incidence)[2])
        if ~all(~t)
            tv[:,find(t)] = -dl' # with julia 0.3.2, sparse([1 -1]).' -> [1, 0], hence use of ' (conjugate transpose)
        end
        tv[:,find(~t)] = speye(T,size(dl)[2])
    end

    tv, ti
end

topomat{T<:Integer}(incidence::SparseMatrixCSC{T}) = topomat!(copy(incidence))
topomat(c::Circuit) = topomat!(incidence(c))

type DiscreteModel{Solver}
    a::Matrix{Float64}
    b::Matrix{Float64}
    c::Matrix{Float64}
    x0::Vector{Float64}
    pexp::Matrix{Float64}
    dq::Matrix{Float64}
    eq::Matrix{Float64}
    fq::Matrix{Float64}
    q0::Vector{Float64}
    dy::Matrix{Float64}
    ey::Matrix{Float64}
    fy::Matrix{Float64}
    y0::Vector{Float64}

    nonlinear_eq :: Expr

    solver::Solver
    x::Vector{Float64}

    function DiscreteModel(circ::Circuit, t::Float64)
        model = new()

        mats = model_matrices(circ, t)
        for mat in [:a, :b, :c, :pexp, :dq, :eq, :fq, :dy, :ey, :fy]
            setfield!(model, mat, convert(Matrix{Float64}, full(mats[mat])))
        end
        for vect in [:x0, :q0, :y0]
            setfield!(model, vect, convert(Vector{Float64}, vec(full(mats[vect]))))
        end

        @assert nn(circ) == nn(model)

        model.x = zeros(nx(model))

        model.nonlinear_eq = quote
            #copy!(q, q0 + pexp * p + fq * z)
            copy!(q, q0)
            BLAS.gemv!('N',1.,pexp,p,1.,q)
            BLAS.gemv!('N',1.,fq,z,1.,q)
            let J=Jq
                $(nonlinear_eq(circ))
            end
            #copy!(J, Jq*model.fq)
            BLAS.gemm!('N', 'N', 1., Jq, fq, 0., J)
            #copy!(Jp, Jq*pexp)
            BLAS.gemm!('N', 'N', 1., Jq, pexp, 0., Jp)
        end

        init_z = initial_solution(model.nonlinear_eq, model.q0, model.fq)
        nonlinear_eq_func = eval(quote
            if VERSION < v"0.5.0-dev+2396"
                # wrap up in named function because anonymous functions are slow
                # in old Julia versions
                function $(gensym())(res, J, Jp, p, z)
                    q0=$(model.q0)
                    pexp=$(model.pexp)
                    q=$(zeros(nq(model)))
                    Jq=$(zeros(nn(model), nq(model)))
                    fq=$(model.fq)
                    $(model.nonlinear_eq)
                    return nothing
                end
            else
                (res, J, Jp, p, z) ->
                    let q0=$(model.q0), pexp=$(model.pexp),
                        q=$(zeros(nq(model))),
                        Jq=$(zeros(nn(model), nq(model))), fq=$(model.fq)
                        $(model.nonlinear_eq)
                        return nothing
                    end
            end
        end)
        model.solver =
            Solver(ParametricNonLinEq(nonlinear_eq_func, nn(model), np(model)),
                   zeros(np(model)), init_z)
        return model
    end
end

DiscreteModel(circ::Circuit, t::Float64) =
    DiscreteModel{HomotopySolver{CachingSolver{SimpleSolver}}}(circ, t)

function model_matrices(circ::Circuit, t)
    x, f = gensolve([mv(circ) mi(circ) 1/t*mxd(circ)+0.5*mx(circ) mq(circ);
                     blkdiag(topomat(circ)...) spzeros(nb(circ), nx(circ) + nq(circ))],
                    [u0(circ) mu(circ) 1/t*mxd(circ)-0.5*mx(circ);
                     spzeros(nb(circ), 1+nu(circ)+nx(circ))])

    rowsizes = [nb(circ); nb(circ); nx(circ); nq(circ)]
    res = Dict{Symbol,AbstractArray}(zip([:fv; :fi; :c; :fq], matsplit(f, rowsizes)))

    indeterminates = zeros(size(f,1), 0)
    if size(res[:fq], 2) > nn(circ)
        fq = full(res[:fq])
        q, r, piv = qr(fq.', Val{true}, thin=false)
        # fq[piv,:]*q ~ r'
        indeterminates = f * q[:,nn(circ)+1:end]
        if sumabs2(res[:c] * q[:,nn(circ)+1:end]) > 1e-20
            warn("State update depends on indeterminate quantity")
        end
        f *= q[:,1:nn(circ)]
        for k in [:fv; :fi; :c; :fq]
            res[k] *= q[:,1:nn(circ)]
        end
    end

    # choose particular solution such that the rows corresponding to q are
    # column-wise orthogonal to the column space of fq (and hence have a column
    # space of minimal dimension)
    fq = res[:fq]
    x = x - full(f)/(full(fq)'*fq)*fq'*x[end-nq(circ)+1:end,:]

    merge!(res, Dict(zip([:v0 :ev :dv; :i0 :ei :di; :x0 :b :a; :q0 :eq_full :dq_full],
                         matsplit(x, rowsizes, [1; nu(circ); nx(circ)]))))

    dq_full = res[:dq_full]
    eq_full = res[:eq_full]
    if size(dq_full)[1] > 0
        # decompose [dq_full eq_full] into pexp*[dq eq] with [dq eq] having minimum
        # number of rows based on the QR factorization
        pexp, r, piv = qr([dq_full eq_full], Val{true})
        ref_err = vecnorm([dq_full eq_full][:,piv] - pexp*r)
        local rowcount
        for rowcount=size(r)[1]:-1:0
            err = vecnorm([dq_full eq_full][:,piv] - pexp[:,1:rowcount]*r[1:rowcount,:])
            if err > ref_err
                rowcount += 1
                break
            end
        end
        res[:dq], res[:eq] =
            matsplit(r[1:rowcount,sortperm(piv)], [rowcount], [nx(circ); nu(circ)])
        res[:pexp] = pexp[:,1:rowcount]
    else
        res[:dq] = zeros(0, nx(circ))
        res[:eq] = zeros(0, nu(circ))
        res[:pexp] = zeros(0, 0)
    end

    p = [pv(circ) pi(circ) 0.5*px(circ)+1/t*pxd(circ) pq(circ)]
    if sumabs2(p * indeterminates) > 1e-20
        warn("Model output depends on indeterminate quantity")
    end
    res[:dy] = p * x[:,2+nu(circ):end] + 0.5*px(circ)-1/t*pxd(circ)
    #          p * [dv; di; a;  dq_full] + 0.5*px(circ)-1/t*pxd(circ)
    res[:ey] = p * x[:,2:1+nu(circ)] # p * [ev; ei; b;  eq_full]
    res[:fy] = p * f                 # p * [fv; fi; c;  fq]
    res[:y0] = p * x[:,1]            # p * [v0; i0; x0; q0]

    return res
end

function initial_solution(nleq, q0, fq)
    # determine an initial solution with a homotopy solver that may vary q0
    # between 0 and the true q0 -> q0 takes the role of p
    nq, nn = size(fq)
    init_nl_eq_func = eval(quote
        (res, J, Jp, p, z) ->
            let q0=$(zeros(nq)), pexp=$(eye(nq)), q=$(zeros(nq)),
                Jq=$(zeros(nn, nq)), fq=$(fq)
                $(nleq)
                return nothing
            end
    end)
    init_nleq = ParametricNonLinEq(init_nl_eq_func, nn, nq)
    init_solver = HomotopySolver{SimpleSolver}(init_nleq, zeros(nq), zeros(nn))
    init_z = solve(init_solver, q0)
    if !hasconverged(init_solver)
        error("Failed to find initial solution")
    end
    return init_z
end

nx(model::DiscreteModel) = length(model.x0)
nq(model::DiscreteModel) = length(model.q0)
np(model::DiscreteModel) = size(model.dq, 1)
nu(model::DiscreteModel) = size(model.eq, 2)
ny(model::DiscreteModel) = length(model.y0)
nn(model::DiscreteModel) = size(model.fq, 2)

function steadystate(model::DiscreteModel, u=zeros(nu(model)))
    IA_LU = lufact(eye(nx(model))-model.a)
    steady_nl_eq_func = eval(quote
        (res, J, Jp, p, z) ->
            let q0=$(zeros(nq(model))), pexp=$(eye(nq(model))),
                q=$(zeros(nq(model))), Jq=$(zeros(nn(model), nq(model))),
                fq=$(model.pexp*model.dq/IA_LU*model.c + model.fq)
                $(model.nonlinear_eq)
                return nothing
            end
    end)
    steady_nleq = ParametricNonLinEq(steady_nl_eq_func, nn(model), nq(model))
    steady_solver = HomotopySolver{SimpleSolver}(steady_nleq, zeros(nq(model)),
                                                 zeros(nn(model)))
    set_resabs2tol!(steady_solver, 1e-30)
    steady_q0 = model.q0 + model.pexp*(model.dq/IA_LU*model.b + model.eq)*u +
                model.pexp*model.dq/IA_LU*model.x0
    steady_z = solve(steady_solver, steady_q0)
    if !hasconverged(steady_solver)
        error("Failed to find steady state solution")
    end
    return IA_LU\(model.b*u + model.c*steady_z + model.x0)
end

function steadystate!(model::DiscreteModel, u=zeros(nu(model)))
    x_steady = steadystate(model, u)
    copy!(model.x, x_steady)
    return x_steady
end

function run!(model::DiscreteModel, u::AbstractMatrix{Float64})
    y = Array(Float64, ny(model), size(u)[2])
    ucur = Array(Float64, nu(model))
    p = Array(Float64, np(model))
    ycur = Array(Float64, ny(model))
    xnew = Array(Float64, nx(model))
    @showprogress 1 "Running model: " for n = 1:size(u)[2]
        copy!(ucur, u[:,n])
        # copy!(p, model.dq * model.x + model.eq * u[:,n])
        BLAS.gemv!('N', 1., model.dq, model.x, 0., p)
        BLAS.gemv!('N', 1., model.eq, ucur, 1., p)
        z = solve(model.solver, p)
        if ~hasconverged(model.solver)
            if all(isfinite(z))
                warn("Failed to converge while solving non-linear equation.")
            else
                error("Failed to converge while solving non-linear equation, got non-finite result.")
            end
        end
        #y[:,n] = model.dy * model.x + model.ey * u[:,n] + model.fy * z + model.y0
        copy!(ycur, model.y0)
        BLAS.gemv!('N', 1., model.dy, model.x, 1., ycur)
        BLAS.gemv!('N', 1., model.ey, ucur, 1., ycur)
        BLAS.gemv!('N', 1., model.fy, z, 1., ycur)
        y[:,n] = ycur
        #model.x = model.a * model.x + model.b * u[:,n] + model.c * z + model.x0
        copy!(xnew, model.x0)
        BLAS.gemv!('N', 1., model.a, model.x, 1., xnew)
        BLAS.gemv!('N', 1., model.b, ucur, 1.,xnew)
        BLAS.gemv!('N', 1., model.c, z, 1., xnew)
        copy!(model.x, xnew)
    end
    return y
end

function swaprows!(a::SparseMatrixCSC, row1, row2)
    # This sometimes gives a wrong result with julia 0.3.2:
    #    a[[row1, row2],:] = a[[row2, row1],:]
    i, j, v = findnz(a)
    row1idx = i.==row1
    row2idx = i.==row2
    i[row1idx] = row2
    i[row2idx] = row1
    a[:,:] = sparse(i, j, v, size(a)...)
end

# lines marked with !SV avoid creation of SparseVector by indexing with Ranges
# instead of Ints; a better way for cross-julia-version compatibilty would be
# nice; maybe Compat helps in the future...
function gensolve(a::SparseMatrixCSC, b, x, h, thresh=0.1)
    m = size(a)[1]
    t = sortperm(vec(sum(spones(a),2))) # row indexes in ascending order of nnz
    for i in 1:m
        ait = a[t[i:i],:] # ait is a row of the a matrix # !SV
        s = ait * h;
        if nnz(s) == 0
            continue
        end
        inz, jnz, nz_vals = findnz(s)
        nz_abs_vals = abs(nz_vals)
        jat = jnz[nz_abs_vals .≥ thresh*maximum(nz_abs_vals)] # cols above threshold
        j = jat[indmin(sum(spones(h[:,jat])))]
        q = h[:,j:j] # !SV
        # ait*q only has a single element!
        x = x + convert(typeof(x), q * ((b[t[i:i],:] - ait*x) * (1 / (ait*q)[1]))); # !SV
        if size(h)[2] > 1
            h = h[:,[1:j-1;j+1:end]] - convert(typeof(h), q * s[1:1,[1:j-1;j+1:end]]*(1/s[1,j])) # !SV
        else
            h = similar(h, eltype(h), (size(h)[1], 0))
        end
    end
    return x, h
end

gensolve(a, b, thresh=0.1) =
    gensolve(a, b, spzeros(size(a)[2], size(b)[2]), speye(size(a)[2]), thresh)

consecranges(lengths) = map(range, cumsum([1; lengths[1:end-1]]), lengths)

if VERSION < v"0.4.0-dev+2833"
    # work around a bug when indexing with empty ranges
    function matsplit(m, rowsizes, colsizes=[size(m)[2]])
        fixed_getindex(A, rs, cs) = isempty(rs) || isempty(cs) ?
            similar(A, length(rs), length(cs)) : return A[rs, cs]
        return [fixed_getindex(m, rs, cs) for rs in consecranges(rowsizes),
                                              cs in consecranges(colsizes)]
    end
else
    matsplit(m, rowsizes, colsizes=[size(m)[2]]) =
        [m[rs, cs] for rs in consecranges(rowsizes), cs in consecranges(colsizes)]
end

end # module


text:
using Distributions: Normal, MixtureModel, pdf
using Makie

function main()
    d = Normal(0, 1)
    a=0.5
    b=3.0
    N(x) = pdf(d, x)
    p(x) = (1 - a) * N(x) + a *N(x - b)
    mixmodel = MixtureModel(
        Normal[
            Normal(0., 1.0),
            Normal(b, 1.0)
        ],
        [1-a, a]
    )
    xs = -5:0.01:5
    s=lines(xs, p.(xs))
    lines!(s, xs, pdf.(mixmodel, xs))
end

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using SpecialGraphs
using Test

import LightGraphs
const LG = LightGraphs
using LightGraphs: Edge, edges

@testset "WheelGraph" begin
    @test !LG.is_directed(WheelGraph)
    @test !LG.is_directed(WheelGraph{Int})
    @test !LG.is_directed(WheelGraph{UInt})
    wg = WheelGraph(10)
    wgref = LG.WheelGraph(10)
    @test LG.edgetype(wg) <: Edge
    @test LG.has_vertex(wg, 10)
    @test !LG.has_vertex(wg, 11)
    @test !LG.has_edge(wg, 10, 11)
    @test LG.has_edge(wg, 10, 2) # boundary condition
    @test all(LG.pagerank(wg) ≈ LG.pagerank(wgref))
    @test LG.ne(wg) == LG.ne(wgref)
    for ninit in (1, 2, 5)
        @test all(LG.bfs_parents(wg, ninit) .== LG.bfs_parents(wgref, ninit))
        @test all(LG.dfs_parents(wg, ninit) .== LG.dfs_parents(wgref, ninit))
    end
    for i in 2:LG.nv(wg)-1
        e = LG.edges(wg)
        @test Edge(i, i+1) in e
        @test LG.has_edge(wg, i, i+1)
        @test Edge(1, i) in e
        @test LG.has_edge(wg, 1, i)
    end
    @test !LG.has_edge(wg, 2, 4)
end

@testset "PathGraph" begin
    pg = PathGraph(10)
    @test !LG.is_directed(PathGraph)
    @test !LG.is_directed(PathGraph{Int})
    @test !LG.is_directed(PathGraph{UInt})
    @test LG.edgetype(pg) <: Edge
    pgref = LG.PathGraph(10)
    @test LG.has_vertex(pg, 10)
    @test !LG.has_vertex(pg, 11)
    @test !LG.has_edge(pg, 10, 11)
    @test all(LG.pagerank(pg) ≈ LG.pagerank(pgref))
    @test LG.ne(pg) == LG.ne(pgref)
    for ninit in (1, 2, 5)
        @test all(LG.dfs_parents(pg, ninit) .== LG.dfs_parents(pgref, ninit))
        @test all(LG.dfs_parents(pg, ninit) .== LG.dfs_parents(pgref, ninit))
    end
    e = edges(pg)
    for v in 2:9
        @test length(LG.outneighbors(pg, v)) == 2
        @test LG.has_edge(pg, v-1, v)
        @test Edge(v-1, v) in e
    end
    @test LG.has_edge(pg, 9, 10)
end

@testset "CompleteGraph" begin
    cg = CompleteGraph(10)
    @test !LG.is_directed(CompleteGraph)
    @test !LG.is_directed(CompleteGraph{Int})
    @test !LG.is_directed(CompleteGraph{UInt})
    @test LG.edgetype(cg) <: Edge
    cgref = LG.CompleteGraph(10)
    @test !LG.has_vertex(cg, 0)
    @test !LG.has_vertex(cg, 11)
    @test !LG.has_edge(cg, 10, 11)
    @test LG.has_vertex(cg, 10)
    @test LG.ne(cg) == LG.ne(cgref)
    @test all(LG.pagerank(cg) ≈ LG.pagerank(cgref))
    for ninit in (1, 2, 5)
        @test all(LG.dfs_parents(cg, ninit) .== LG.dfs_parents(cgref, ninit))
        @test all(LG.dfs_parents(cg, ninit) .== LG.dfs_parents(cgref, ninit))
    end
    for v in LG.vertices(cg)
        @test length(LG.outneighbors(cg, v)) == LG.nv(cg) - 1
        e = LG.edges(cg)
        for v2 in LG.vertices(cg)
            if v == v2
                @test !LG.has_edge(cg, v, v2)
                @test !in(Edge(v, v2), e)
            else
                @test LG.has_edge(cg, v, v2)
                if v < v2
                    @test Edge(v, v2) in e
                end
            end
        end
    end
end


text:
struct MaintenanceSchedulingData{T <: AbstractFloat} <: AbstractModelData
    hydrodata::HydroPlantCollection{T,2}
    bidlevels::Vector{Vector{T}}
    reschedule::Bool

    function MaintenanceSchedulingData(plantdata::HydroPlantCollection{T,2},
                                       bidlevels::Vector{Vector{T}},
                                       reschedule::Bool) where T <: AbstractFloat
        length(bidlevels) == 24 || error("Supply exactly 24 bidlevel sets")
        all(length.(bidlevels) .== length(bidlevels[1])) || error("All bidlevel sets must be of the same length.")
        return new{T}(plantdata, bidlevels, reschedule)
    end
end

function MaintenanceSchedulingData(plantfilename::String, bidlevelsets::AbstractMatrix; reschedule::Bool = false)
    size(bidlevelsets, 2) == 24 || error(" ")
    plantdata = HydroPlantCollection(plantfilename)
    bidlevels = [bidlevelsets[:,i] for i = 1:size(bidlevelsets, 2)]
    for i in eachindex(bidlevels)
        prepend!(bidlevels[i], -500.0)
        push!(bidlevels[i], 3000.0)
    end
    MaintenanceSchedulingData(plantdata, bidlevels, reschedule)
end

function MaintenanceSchedulingData(plantfilename::String, bidlevels::AbstractVector; reschedule::Bool = false)
    plantdata = HydroPlantCollection(plantfilename)
    prepend!(bidlevels, -500.0)
    push!(bidlevels, 3000.0)
    MaintenanceSchedulingData(plantdata, fill(bidlevels, 24), reschedule)
end


text:
module CosmoParticles

using FastUnitfulOperations
using FillArrays
using LazyArrays
using LinearAlgebra
using Statistics
using StatsBase
using Unitful

export
    # types
    AbstractParticles,
    Particles,
    AllParticles,
    AbstractParticleCollection,
    ParticleCollection,
    RedshiftParticleCollection,

    # geometry
    AbstractCosmoGeometry,
    CosmoHyperrectangle,
    CosmoCuboid,
    CosmoRectangle,
    CosmoHypercube,
    CosmoCube,
    CosmoSquare,
    CosmoHypersphere,
    CosmoSphere,
    CosmoCircle,
    CosmoCylinder,
    CosmoStandingCylinder,

    # getters
    redshift,

    # operations
    delete,

    # transformations
    rotate,
    rotate!,
    translate,
    translate!,
    translate_periodic,
    translate_periodic!,
    translate_periodic_to_center,
    translate_periodic_to_center!,
    to_comoving,
    to_comoving!,
    to_physical,
    to_physical!,

    # properties
    meanprop,
    sumprop,
    meanpos,
    meanvel,
    colnorm,
    colnorm2,
    colcross,
    coldot,
    angmom,
    angmomtot,
    angmomtot_stable

include("abstractparticles.jl")
include("particles.jl")
include("abstractparticlecollection.jl")
include("particlecollection.jl")
include("allparticles.jl")
include("geometry.jl")

include("utils.jl")
include("transformations.jl")
include("miscoperations.jl")
include("properties.jl")

end


text:
using PyCall
using Plots

np = pyimport("numpy")

function main()

    # x = 5
    # x = [π 3.0 0 12312;
    #      5 12. 1 3412.]
    # display(np.sin(x))

    x = [1, 2, 3, 4, 5]
    y = [5, 3, 1, 3, 5]

    p = np.polyfit(x, y, 2)
    # display(p)

    x1 = 1:0.1:5
    p1 = np.polyval(p, x1)

    scatter(x, y, label="data")
    plot!(x1, p1, label="fit")

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
#=doc
Utilities
=#

export crop_array_border

"""
throwaway the border region of an array (ndims > 3), currently only works for 3D cropsize.
"""
function crop_array_border(arr::Array, cropsize::Union{Vector,Tuple})
    @assert ndims(arr) >= 3
    sz = size(arr)
    @assert sz[1]>cropsize[1]*2 &&
            sz[2]>cropsize[2]*2 &&
            sz[3]>cropsize[3]*2
    return arr[ cropsize[1]+1:sz[1]-cropsize[1],
                cropsize[2]+1:sz[2]-cropsize[2],
                cropsize[3]+1:sz[3]-cropsize[3]]
end


text:
module ArrayStreams

using ImageMetadata, ImageCore, StaticArrays, Distributions, ColorTypes
import AxisArrays
using AxisArrays: AxisArray, axisnames, permutation
using GeometryTypes: Triangle, Point
using Rotations: Quat

export AbstractArrayStream, ArrayStream, stream

function fixendian(A::AbstractArray, needswap::Bool)
    if needswap
        return mappedarray((ntoh, hton), A)
    else
        return A
    end
end

expand2indices(A::AbstractArray{T,NA}, idx::Tuple{Vararg{<:Any,NI}}) where {T,NA,NI} =
    to_indices(A,(idx..., ntuple(i->Colon(), NA-NI)...))
expand2indices(A::AbstractArray{T,N}, idx::Tuple{Vararg{<:Any,N}}) where {T,N} =
    to_indices(A,idx)
BitTypes = Union{Integer,AbstractFloat}

include("iotypes/iotypes.jl")
include("abstractarraystream.jl")
include("arraystream.jl")
include("imagemeta.jl")
include("axisarray.jl")

# TODO
# - write IOBuffer (i.e. write(s) where there's not Array parameters but buffer contains
# - documentation
#   - write should be explained in the context of sinks because ArrayStream acts as the sink in this case
# - check extractquat for proper output shape
#   - Distributions, Quat, SArray
# - iterators for chunkedstream
# - constructors for chunkedstream
# - Impliment Colors
# - append! operations (maybe, as ArrayStreams are more like StaticArrays)

end


text:
import Base: undef_ref_str

const dump_indent = "    "

richdump(io::IO, dom::Vector{Node}, @nospecialize(x)) = richdump(IOContext(io), dom, x)

function richdump(io::IOContext, dom::Vector{Node}, @nospecialize(x))
    let id = startNode(dom, io, x, "")
        T = typeof(x)
        isa(x, Function) && print(io, " (function of type ")
        richprint(io, dom, T)
        isa(x, Function) && print(io, ")")
        nf = nfields(x)
        if nf > 0
            if !Base.show_circular(io, x)
                recur_io = IOContext(io, Pair{Symbol,Any}(:SHOWN_SET, x))
                for field in 1:nf
                    println(io)
                    nextChild(dom, io, id)
                    fname = string(fieldname(T, field))
                    print(io, dump_indent, fname, ": ")
                    if isdefined(x, field)
                        richprint(recur_io, dom, getfield(x, field))
                    else
                        print(io, undef_ref_str)
                    end
                end
            end
        else
            if !isa(x, Function)
                print(io, dump_indent)
                richprint(io, dom, x)
            end
        end
        endNode(dom, io, id)
    end
end

function richdump(io::IOContext, dom::Vector{Node}, x::Core.SimpleVector)
    let id = startNode(dom, io, x, "SimpleVector")
        if isempty(x)
            print(io, "empty SimpleVector")
        else
            print(io, "SimpleVector")
            for i = 1:length(x)
                println(io)
                nextChild(dom, io, id)
                print(io, dump_indent, i, ": ")
                if isassigned(x, i)
                    richprint(io, dom, x[i])
                else
                    print(io, undef_ref_str)
                end
            end
        end
        endNode(dom, io, id)
    end
end

richdump(io::IOContext, dom::Vector{Node}, x::Module) =
    let id = startNode(dom, io, x, "Module")
        print(io, "Module ", x)
        endNode(dom, io, id)
    end
richdump(io::IOContext, dom::Vector{Node}, x::String) =
    let id = startNode(dom, io, x, "String")
        print(io, "String ")
        show(io, x)
        endNode(dom, io, id)
    end
richdump(io::IOContext, dom::Vector{Node}, x::Symbol) =
    let id = startNode(dom, io, x, "Symbol")
        print(io, "Symbol ", x)
        endNode(dom, io, id)
    end
richdump(io::IOContext, dom::Vector{Node}, x::Ptr) =
    let id = startNode(dom, io, x, "Ptr")
        print(io, x) # this avoids print the type twice
        endNode(dom, io, id)
    end

function dump_elts(io::IOContext, dom::Vector{Node}, id::Int, x::Array, i0, i1)
    for i in i0:i1
        println(io)
        nextChild(dom, io, id)
        print(io, dump_indent, i, ": ")
        if !isassigned(x, i)
            print(io, undef_ref_str)
        else
            richprint(io, dom, x[i])
        end
    end
end

function richdump(io::IOContext, dom::Vector{Node}, x::Array{T}) where {T}
    let id = startNode(dom, io, x, "Array")
        print(io, "Array{")
        richprint(io, dom, T)
        print(io, "}(", string(size(x)), ")")
        if isprimitivetype(T) # was T <: Number
            print(io, " ")
            richprint(io, dom, x)
        else
            if !isempty(x) && !Base.show_circular(io, x)
                recur_io = IOContext(io, :SHOWN_SET => x)
                lx = length(x)
                if get(io, :limit, false) && lx > 10
                    dump_elts(recur_io, dom, id, x, 1, 5)
                    println(io)
                    print(io, "  ...")
                    dump_elts(recur_io, dom, id, x, lx - 4, lx)
                else
                    dump_elts(recur_io, dom, id, x, 1, lx)
                end
            end
        end
        endNode(dom, io, id)
    end
end

# Types
function richdump(io::IOContext, dom::Vector{Node}, x::DataType)
    let id = startNode(dom, io, x, "DataType")
        richprint(io, dom, x)
        if x !== Any
            print(io, " <: ")
            richprint(io, dom, supertype(x))
        end
        if !(x <: Tuple) && !x.abstract
            tvar_io::IOContext = io
            for tparam in x.parameters
                # approximately recapture the list of tvar parameterization
                # that may be used by the internal fields
                if isa(tparam, TypeVar)
                    tvar_io = IOContext(tvar_io, :unionall_env => tparam)
                end
            end
            fieldtypes = x.types
            for idx in 1:length(fieldtypes)
                println(io)
                nextChild(dom, io, id)
                print(io, dump_indent, fieldname(x, idx), "::")
                richprint(tvar_io, dom, fieldtypes[idx])
            end
        end
        endNode(dom, io, id)
    end
end

function richdump(io::IOContext, dom::Vector{Node}, x::UnionAll)
    let id = startNode(dom, io, x, "UnionAll")
        print(io, "UnionAll where var: {")
        tvar_io::IOContext = io
        while true
            println(io)
            nextChild(dom, io, id)
            print(io, dump_indent)
            richprint(tvar_io, dom, x.var)
            tvar_io = IOContext(tvar_io, :unionall_env => x.var)
            x = x.body
            if x isa UnionAll
                print(io, ",")
            else
                break
            end
        end
        println(io)
        println(io, "}")
        print(io, "body: ")
        nextChild(dom, io, id)
        richdump(tvar_io, dom, x)
        endNode(dom, io, id)
    end
end

function richdump(io::IOContext, dom::Vector{Node}, x::Union)
    let id = startNode(dom, io, x, "Union")
        richprint(io, dom, x)
        endNode(dom, io, id)
    end
end

max_stars_repo_name: vtjnash/Firelight


text:
using Downloads

const current_dir = @__DIR__
const number_trees = [1, 1, 2, 3, 6, 11, 23, 47, 106, 235, 551, 1301, 3159]

list_g = open(dirname(current_dir)*"/src/list_trees.txt")
for n in 1:4 #d=4

    Dim_dir::String = dirname(current_dir)*"/Data/Dimension$n" #path of the folder containing the colorations
    mkpath(Dim_dir) #create the folder

    seekstart(list_g)
    #open the file containing the list of Prufer sequences of graphs

    for v in 2:(4+1) #run the computation among all graphs with fixed number of vertices

        n_trees_nv = number_trees[v - 1]  #we known how many graphs there are with fixed number of vertices

        for _ in 1:n_trees_nv  #run the computation for a fixed graph

            str = readline(list_g) #read a new line, we expect a Prufer seq plus the number of automorphisms

            name_file = string(split(str, ',')[1],"0.gz")
            if !(name_file in readdir(Dim_dir))
                url = "https://raw.githubusercontent.com/mgemath/Colorations/main/Dimension$n/$name_file"
                dest = Dim_dir*"/$name_file"
                try
                    Downloads.download(url, dest)
                catch e
                    printstyled(stderr,"ERROR: ", bold=true, color=:red)
                    printstyled(stderr,sprint(showerror,e), color=:light_red)
                    println(stderr)
                    final_state = false
                    break  #end for n_g
                end
            end

        end

    end
end
close(list_g)


text:
tic()
include("Include.jl")

# load the data dictionary -
data_dictionary = DataDictionary(0,0,0)
TXTL_dictionary = TXTLDictionary(0,0,0)
number_of_fluxes = length(data_dictionary["default_flux_bounds_array"][:,1])
number_of_species = length(data_dictionary["species_bounds_array"][:,1])
#Set objective reaction
data_dictionary["objective_coefficient_array"][171] = -1;
#-------------------------------------------------------------------------------------------------#Define case number
# 1 = Amino Acid Uptake & Synthesis
# 2 = Amino Acid Uptake w/o Synthesis
# 3 = Amino Acid Synthesis w/o Uptake
case = 1
if case == 1
  data_dictionary["AASyn"] = 100;
  data_dictionary["AAUptake"] = 30
  data_dictionary["AASecretion"] = 0;
end
if case == 2
  data_dictionary["AASyn"] = 0;
  data_dictionary["AAUptake"] = 30
  data_dictionary["AASecretion"] = 0;
end
if case == 3
  data_dictionary["AASyn"] = 100;
  data_dictionary["AAUptake"] = 0
  data_dictionary["AASecretion"] = 100;
end
#-------------------------------------------------------------------------------------------------
volume = TXTL_dictionary["volume"]
plasmid_concentration = 5;
gene_copy_number = (volume/1e9)*(6.02e23)*plasmid_concentration;
TXTL_dictionary["gene_copies"] = gene_copy_number
data_dictionary = DataDictionary(0,0,0)
#Set objective reaction
data_dictionary["objective_coefficient_array"][194] = -1;
data_dictionary["AASyn"] = 100;
data_dictionary["AAUptake"] = 30
data_dictionary["AASecretion"] = 0;
data_dictionary["GlcUptake"] = 30#*rand(1)[1];
data_dictionary["Oxygen"] = 30#+2*rand(1)[1];
TXTL_dictionary["RNAP_concentration_nM"] = 75 #+ (80-60)*rand(1)[1];
TXTL_dictionary["RNAP_elongation_rate"] = 25 #+ (30-20)*rand(1)[1];
TXTL_dictionary["RIBOSOME_concentration"] = 0.0016# + (0.0018-0.0012)*rand(1)[1];
TXTL_dictionary["RIBOSOME_elongation_rate"] = 2# + (3-1.5)*rand(1)[1];
data_dictionary["Oxygen"] = 100;
data_dictionary["AAUptake"] = 1 #+ (1.5-0.5)*rand(1)[1];
data_dictionary["GlcUptake"] = 30# + (38-28)*rand(1)[1];
data_dictionary["AC"] = 8# + (12-6)*rand(1)[1];
data_dictionary["SUCC"] = 2# + (3-1)*rand(1)[1];
data_dictionary["PYR"] = 7# + (10-5)*rand(1)[1];
data_dictionary["MAL"] = 3 #+ (5-2)*rand(1)[1];
data_dictionary["LAC"] = 8 #+ (12-8)*rand(1)[1];
data_dictionary["ENERGY"] = 1#1*rand(1)[1];
data_dictionary["ALA"] = 2 #+ (4-3)*rand(1)[1];
data_dictionary["ASP"] = 1#+rand(1)[1];
data_dictionary["GLN"] = 0.25#*rand(1)[1];
data_dictionary["LysUptake"] = 1.5#*rand(1)[1];
data_dictionary["LysSecretion"] = 0.5#1*rand(1)[1];
#TXTL_dictionary["RNAP_concentration_nM"] = 60 + (80-60)*rand(1)[1];
#TXTL_dictionary["RNAP_elongation_rate"] = 20 + (30-20)*rand(1)[1];
#TXTL_dictionary["RIBOSOME_concentration"] = 0.0012 + (0.0018-0.0012)*rand(1)[1];
#TXTL_dictionary["RIBOSOME_elongation_rate"] = 1.5 + (3-1.5)*rand(1)[1];
#-------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------
t0_exp = 0
tf_exp = 3
tstep_exp = 0.01
experimental_time_exp = [t0_exp:tstep_exp:tf_exp;]
#Mean = Unpack("./data/Mean")
#Std = Unpack("./data/Std")
#Upper = Unpack("./data/Upper")
#Lower = Unpack("./data/Lower")

species_list = readdlm("./data/species_list")
Mean_raw = readdlm("./data/Mean")
Std_raw = readdlm("./data/Std")
Upper_raw = readdlm("./data/Upper")
Lower_raw = readdlm("./data/Lower")
Mean = Dict()
Std = Dict()
Upper = Dict()
Lower = Dict()
for i in 1:length(species_list)
	rxn_idx = species_list[i]
	Mean[rxn_idx] = Mean_raw[:,i]
	Std[rxn_idx] = Std_raw[:,i]
	Upper[rxn_idx] = Upper_raw[:,i]
	Lower[rxn_idx] = Lower_raw[:,i]
end

#amplitude = 2, decay_rate = 1
data_dictionary["species_drift_amplitude"] = 2 # percent allowed to drift for unconstraint species for each iteration (0.01 hr)
data_dictionary["species_drift_decay_rate"] = 1
metabolite_list = data_dictionary["list_of_metabolite_symbols"]
rxn_list = data_dictionary["list_of_reaction_strings"]

(n_species,n_rxn) = size(data_dictionary["stoichiometric_matrix"])
t0 = 0
tf = 3
tstep = 0.01
experimental_time = t0:tstep:tf
length_time = length(experimental_time)
data_dictionary["tstep"] = tstep

#notes: two possible implementations ---
#  1. select subset of species to be constraint to synthetic data
#    have everything else that are not constrained to data be open or a function of concentration
#  2. select subset of species to be constraint to synthetic data
#    have everything else that are not constrained to data constrained at steady state (sv=0)

tRNA_index = Int64[]
for i in 1:length(metabolite_list)
	rxn_idx = metabolite_list[i]
	if length(rxn_idx) > 5
		if rxn_idx[end-4:end] == "_tRNA"
			push!(tRNA_index,i)
		end
	end
end

initial_conditions = zeros(n_species,1)
for species_index in 1:n_species
	if in(species_index,tRNA_index)
		initial_conditions[species_index] = Mean[metabolite_list[species_index]][2]
	else
		initial_conditions[species_index] = Mean[metabolite_list[species_index]][1]
	end
end
initial_conditions[142] = .001
initial_conditions[144] = .0001

# this determines which parameters are fitted to data
dataset_index = Array{Int64}(readdlm("FVA_SA/species_boolean_array"))
(n_species,n_ens) = size(dataset_index)
non_dataset_index = collect(1:number_of_species)
species_list = data_dictionary["list_of_metabolite_symbols"]
species_constraint_index = copy(dataset_index)

plot_color = "orangered"

data_dictionary["initial_conditions"] = initial_conditions
# Original values:
GENE_IC = 5.19e-6
RNAP_IC = 6.75e-5
RIBO_IC = 0.002
mRNA_IC = 0.0001
tRNA_IC = 1.6
TX_Vmax = 122.7
plasmid_saturation_coefficient = 3.5e-6
TL_Vmax = 163.6
mRNA_saturation_coefficient = 0.045
mRNA_degradation_rate = 20
params = vcat(GENE_IC, RNAP_IC, RIBO_IC, mRNA_IC, tRNA_IC, TX_Vmax, plasmid_saturation_coefficient, TL_Vmax, mRNA_saturation_coefficient, mRNA_degradation_rate)
params[8] = 40
RNAP_IC = 1e-3
RIBO_IC = 0.0021
RNAP_elongation_rate = 25
plasmid_saturation_coefficient = 116e-6
polysome_amplification = 5
RIBOSOME_elongation_rate = 2
mRNA_saturation_coefficient = 0.045
mRNA_degradation_rate = 10
params = vcat(RNAP_IC, RIBO_IC, RNAP_elongation_rate, plasmid_saturation_coefficient, polysome_amplification, RIBOSOME_elongation_rate, mRNA_saturation_coefficient, mRNA_degradation_rate)
data_dictionary["initial_conditions"][143] = params[1]
data_dictionary["initial_conditions"][141] = params[2]
TXTL_dictionary["RNAP_elongation_rate"] = params[3]
TXTL_dictionary["plasmid_saturation_coefficient"] = params[4]
TXTL_dictionary["polysome_amplification"] = params[5]
TXTL_dictionary["RIBOSOME_elongation_rate"] = params[6]
TXTL_dictionary["mRNA_saturation_coefficient"] = params[7]
TXTL_dictionary["mRNA_degradation_rate"] = params[8]

single_species_additions = [3 "M_13dpg_c"; 5 "M_2pg_c"; 6 "M_3pg_c"; 44 "M_dhap_c"; 48 "M_f6p_c"; 50 "M_fdp_c"; 55 "M_g3p_c"; 56 "M_g6p_c"; 105 "M_pep_c"; 4 "M_2ddg6p_c"; 11 "M_6pgc_c"; 12 "M_6pgl_c"; 46 "M_e4p_c"; 117 "M_r5p_c"; 118 "M_ru5p_D_c"; 119 "M_s7p_c"; 138 "M_xu5p_D_c"; 16 "M_accoa_c"; 21 "M_akg_c"; 36 "M_cit_c"; 54 "M_fum_c"; 65 "M_glx_c"; 78 "M_icit_c"; 102 "M_oaa_c"; 124 "M_succoa_c"]

species_constraint_array = dataset_index
Labels = String[]
for i in 1:n_ens
	#species_constraint_array[i,1:end] = dataset_index

	push!(Labels,"SA_$i")
end

#num_constraint_sets = size(species_constraint_array,1)
num_constraint_sets = n_ens
# not currently used
parameter = ones(202,2)
parameter[:,1] = parameter[:,1]*110
parameter[:,2] = parameter[:,2]*0.03

time_state_array = zeros(n_species,length_time)
time_state_array[:,1] = initial_conditions
time_flux_array = zeros(n_rxn,length_time)
FBA = data_dictionary["default_flux_bounds_array"]
SBA = data_dictionary["species_bounds_array"]

EXIT_FLAG = zeros(length_time-1,num_constraint_sets)

overall_state_array = zeros()
# loop for all constraint subsets
FVA_run = true

constraint_index_array = collect(1:num_constraint_sets)

#Accuracy = zeros(num_constraint_sets,1)
#Precision = zeros(num_constraint_sets,1)
for constraint_index in constraint_index_array[1:end] ;println(Labels[constraint_index])
	constraint_subset = Labels[constraint_index] #"base_case"
	if !isdir("FVA_SA/$constraint_subset")
		mkdir("FVA_SA/$constraint_subset")
	end

	# discretized dfba
	state_array = deepcopy(initial_conditions)
	index = 1

	species_constraint_index = find(Array{Int64}(species_constraint_array[:,constraint_index]))
  println(species_constraint_index)
#	species_constraint_index = copy(dataset_index)

	syn_data_upper = zeros(n_species,length_time)
	syn_data_lower = zeros(n_species,length_time)
	syn_mean = zeros(n_species,length_time)
	for syn_index in species_constraint_index
		itp = interpolate((experimental_time_exp,),Mean[metabolite_list[syn_index]],Gridded(Linear()))
		Mean_itp = itp[experimental_time]
		itp = interpolate((experimental_time_exp,),Std[metabolite_list[syn_index]],Gridded(Linear()))
		Std_itp = itp[experimental_time]
		itp = interpolate((experimental_time_exp,),Lower[metabolite_list[syn_index]],Gridded(Linear()))
		Lower_itp = itp[experimental_time]
		itp = interpolate((experimental_time_exp,),Upper[metabolite_list[syn_index]],Gridded(Linear()))
		Upper_itp = itp[experimental_time]
		syn_data_upper[syn_index,:] = Upper_itp
		syn_data_lower[syn_index,:] = Lower_itp
		syn_mean[syn_index,:] = Mean_itp
	end
	syn_data = Dict{AbstractString,Any}()
	syn_data["mean"] = syn_mean
	syn_data["upper"] = syn_data_upper
	syn_data["lower"] = syn_data_lower
    syn_data["lower"][3,:] = syn_data_lower[3,:]/5.3

    Exit_flag = Int64[]
    uptake_array = 0
    flux_array = 0

	FVA_min = Dict()
	FVA_max = Dict()
	Percentage_failed_tps_min = Dict()
	Percentage_failed_tps_max = Dict()
	for rxn_idx in 1:194
		FVA_min[rxn_idx] = zeros(length(t0:tstep:tf-tstep),194) # Only record 194 fluxes
		FVA_max[rxn_idx] = zeros(length(t0:tstep:tf-tstep),194) # Only record 194 fluxes
	end

	#dfba
	fva_time = experimental_time[1:end-1]
	for t in 1:length(fva_time)
		time_index = fva_time[t] ;println(time_index)
		# constraint
		data_dictionary["objective_coefficient_array"] = zeros(number_of_fluxes)
		data_dictionary["objective_coefficient_array"][171] = -1

		data_dictionary = Bounds(data_dictionary,TXTL_dictionary,state_array)
		species_constraints = calculate_constraints(state_array,parameter,species_constraint_index,syn_data,index,tstep)
		data_dictionary["species_bounds_array"] = species_constraints
		#  data_dictionary["default_flux_bounds_array"] = flux_constraints
		data_dictionary["state_array"] = state_array

		# solve the lp problem -
		(objective_value, flux_array, dual_array, uptake_array, exit_flag) = FluxDriver(data_dictionary)
		push!(Exit_flag,exit_flag)

		# mass balance for concentration
		state_array = state_array + (uptake_array * tstep) # cell free mass balanace
		state_array[101] = initial_conditions[101]
		Z0 = objective_value
		if FVA_run && in(time_index,collect(0.0:0.1:3.0))
			for flux_index in 1:194
				#constraint objective flux to value determined
				data_dictionary["default_flux_bounds_array"][171] = 1*Z0
				data_dictionary["default_flux_bounds_array"][171] = Z0

				# minimize flux
				data_dictionary["objective_coefficient_array"] = zeros(number_of_fluxes)
				data_dictionary["objective_coefficient_array"][flux_index] = -1
				(objective_value, flux_array, dual_array, uptake_array, exit_flag) = FluxDriver(data_dictionary)

				rxn_idx = copy(flux_index)

				FVA_min[rxn_idx][t,:] = flux_array[1:194]
				Percentage_failed_tps_min[rxn_idx] = sum(5-exit_flag)/length(experimental_time[1:end-1])

				# maximize flux
				data_dictionary["objective_coefficient_array"] = zeros(number_of_fluxes)
				data_dictionary["objective_coefficient_array"][flux_index] = 1
				(objective_value, flux_array, dual_array, uptake_array, exit_flag) = FluxDriver(data_dictionary)

				FVA_max[rxn_idx][t,:] = flux_array[1:194]
				Percentage_failed_tps_max[rxn_idx] = sum(5-exit_flag)/length(experimental_time[1:end-1])
			end
		end # if FVA_run

		index = index + 1
		time_state_array[:,index] = state_array
		time_flux_array[:,index] = flux_array

	end # for time_index in [t0:tstep:tf-tstep;]

	writedlm("FVA_SA/$constraint_subset/time_state_array.txt",time_state_array)
    writedlm("FVA_SA/$constraint_subset/constraint_subset",species_constraint_index)

	for rxn_idx in 1:194
		if !isdir("FVA_SA/$constraint_subset/$rxn_idx")
			mkdir("FVA_SA/$constraint_subset/$rxn_idx")
		end
		writedlm("FVA_SA/$constraint_subset/$rxn_idx/FVA_min",FVA_min[rxn_idx])
		writedlm("FVA_SA/$constraint_subset/$rxn_idx/FVA_max",FVA_max[rxn_idx])
		writedlm("FVA_SA/$constraint_subset/$rxn_idx/Percentage_failed_tps_min",Percentage_failed_tps_min[rxn_idx])
		writedlm("FVA_SA/$constraint_subset/$rxn_idx/Percentage_failed_tps_max",Percentage_failed_tps_max[rxn_idx])

	end

end # for constraint_index in constraint_index_array

plot_flag = false
if plot_flag
	Plot(plot_color)
	PlotSeparate(plot_color)
end

time_elapsed = toc()


text:
module DataFileNames

    import Printf: @sprintf
    import RegexTools: hex_escape

    export dfname
    export parse_dfname, tryparse_dfname
    export dfheads, dfparams, dfparam, dfext
    export isvalid_dfname

    # include("utils.jl")
    include("globals.jl")
    include("dfname.jl")
    include("parse_regex.jl")
    include("parse_dfname.jl")

    VERSION < v"1.6" && include("compat.jl")

    function __init__()
        _set_default_SEPS!()
        _set_regexs!()
    end
end


text:
var_a = -100
var_b = -2

if var_a > 1 || !(var_b < 2)
    println("It's complicated")
end

(var_a > 1 || !(var_b < 2)) && println("It's complicated")

!isdir("some/folder") && mkpath("some/folder")


text:
# Autogenerated wrapper script for MathGL_jll for x86_64-w64-mingw32-cxx03
export mgl, mgl_qt, mgl_qt5

using Qt_jll
using libpng_jll
using Zlib_jll
using GSL_jll
using Libglvnd_jll
JLLWrappers.@generate_wrapper_header("MathGL")
JLLWrappers.@declare_library_product(mgl, "libmgl.dll")
JLLWrappers.@declare_library_product(mgl_qt, "libmgl-qt.dll")
JLLWrappers.@declare_library_product(mgl_qt5, "libmgl-qt5.dll")
function __init__()
    JLLWrappers.@generate_init_header(Qt_jll, libpng_jll, Zlib_jll, GSL_jll, Libglvnd_jll)
    JLLWrappers.@init_library_product(
        mgl,
        "bin\\libmgl.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        mgl_qt,
        "bin\\libmgl-qt.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        mgl_qt5,
        "bin\\libmgl-qt5.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#<fermions>:
before_combiner_signs(T, labelsT, indsT, C, labelsC, indsC, labelsR, indsR) = T
after_combiner_signs(R, labelsR, indsR, C, labelsC, indsC) = R

function contract(T::BlockSparseTensor, labelsT, C::CombinerTensor, labelsC)
  #@timeit_debug timer "Block sparse (un)combiner" begin
  # Get the label marking the combined index
  # By convention the combined index is the first one
  # TODO: consider storing the location of the combined
  # index in preperation for multiplce combined indices
  cpos_in_labelsC = 1
  clabel = labelsC[cpos_in_labelsC]
  c = combinedindex(C)
  labels_uc = deleteat(labelsC, cpos_in_labelsC)
  if labelsC[1] ∉ labelsT
    # Combine
    labelsRc = contract_labels(labelsC, labelsT)
    cpos_in_labelsRc = findfirst(==(clabel), labelsRc)
    labelsRuc = insertat(labelsRc, labels_uc, cpos_in_labelsRc)
    indsRc = contract_inds(inds(C), labelsC, inds(T), labelsT, labelsRc)

    #<fermions>:
    T = before_combiner_signs(T, labelsT, inds(T), C, labelsC, inds(C), labelsRc, indsRc)

    perm = getperm(labelsRuc, labelsT)
    ucpos_in_labelsT = Tuple(findall(x -> x in labels_uc, labelsT))
    Rc = permutedims_combine(T, indsRc, perm, ucpos_in_labelsT, blockperm(C), blockcomb(C))
    return Rc
  else
    # Uncombine
    labelsRc = labelsT
    cpos_in_labelsRc = findfirst(==(clabel), labelsRc)
    # Move combined index to first position
    if cpos_in_labelsRc != 1
      labelsRc_orig = labelsRc
      labelsRc = deleteat(labelsRc, cpos_in_labelsRc)
      labelsRc = insertafter(labelsRc, clabel, 0)
      cpos_in_labelsRc = 1
      perm = getperm(labelsRc, labelsRc_orig)
      T = permutedims(T, perm)
      labelsT = permute(labelsT, perm)
    end
    labelsRuc = insertat(labelsRc, labels_uc, cpos_in_labelsRc)
    indsRuc = contract_inds(inds(C), labelsC, inds(T), labelsT, labelsRuc)

    # <fermions>:
    T = before_combiner_signs(T, labelsT, inds(T), C, labelsC, inds(C), labelsRuc, indsRuc)

    Ruc = uncombine(T, indsRuc, cpos_in_labelsRc, blockperm(C), blockcomb(C))

    # <fermions>:
    Ruc = after_combiner_signs(Ruc, labelsRuc, indsRuc, C, labelsC, inds(C))

    return Ruc
  end
  #end # @timeit
end

function contract(C::CombinerTensor, labelsC, T::BlockSparseTensor, labelsT)
  return contract(T, labelsT, C, labelsC)
end

# Special case when no indices are combined
# XXX: no copy
contract(T::BlockSparseTensor, labelsT, C::CombinerTensor{<:Any,0}, labelsC) = copy(T)


text:
@testset "chainrules.jl" begin
    @testset "BlockDiagonal" begin
        x = [randn(1, 2), randn(2, 2)]
        test_rrule(BlockDiagonal, x)
        test_rrule(BlockDiagonal, x; output_tangent=Tangent{BlockDiagonal}(;blocks=x))

        b = BlockDiagonal(x)
        m = Matrix(b)
        y, pb = rrule(BlockDiagonal, x)
        # want to test `output_tangent=m`, but can't do it directly because of
        # https://github.com/JuliaDiff/ChainRulesTestUtils.jl/issues/199
        # so just compare the tangent to the tangent we get from the same BlockDiagonal
        @test pb(b) == pb(m)
    end

    @testset "Matrix" begin
        D = BlockDiagonal([randn(1, 2), randn(2, 2)])
        test_rrule(Matrix, D)
    end

    @testset "BlockDiagonal * Vector" begin
        D = BlockDiagonal([rand(2, 3), rand(3, 3)])
        v = rand(6)
        test_rrule(*, D, v)
    end

    @testset "ProjectTo" begin
        bd = BlockDiagonal([ones(2, 2), ones(3, 3)])
        project = ProjectTo(bd)
        @test project(ones(5, 5)) == bd
        @test project(adjoint(ones(5, 5))) == bd
        @test project(Diagonal(ones(5))) isa BlockDiagonal
        @test project(Diagonal(ones(5))) == Diagonal(ones(5))
    end
end


text:
using FileUtils

template_package_name = "PACKAGE"

package_names = [
  "ITensors",
  "NDTensors",
  "ITensorGPU",
  "ITensorGaussianMPS",
  "ITensorVisualizationBase",
  "ITensorUnicodePlots",
  "ITensorMakie",
  "ITensorGLMakie",
]

package_ordering = Dict([
  "ITensors" => 1,
  "NDTensors" => 2,
  "ITensorGPU" => 3,
  "ITensorGaussianMPS" => 4,
  "ITensorVisualizationBase" => 5,
  "ITensorUnicodePlots" => 6,
  "ITensorMakie" => 7,
  "ITensorGLMakie" => 8,
])

function bug_report_file(package_name::String)
  return "$(package_name)_bug_report.md"
end
function feature_request_file(package_name::String)
  return "$(package_name)_feature_request.md"
end

for package_name in package_names
  @show package_name

  order = lpad(package_ordering[package_name], 2, "0")

  template_bug_report = bug_report_file(template_package_name)
  new_bug_report = order * "_" * bug_report_file(package_name)

  if isfile(new_bug_report)
    println("File $new_bug_report already exists, skipping")
  else
    println("Copying $template_bug_report to $new_bug_report")
    cp(template_bug_report, new_bug_report)

    println("Replace $template_package_name with $package_name in $new_bug_report")
    replace_in_file(new_bug_report, template_package_name => package_name)

    mv(new_bug_report, joinpath("..", new_bug_report); force=true)
  end

  template_feature_request = feature_request_file(template_package_name)
  new_feature_request = order * "_" * feature_request_file(package_name)

  if isfile(new_feature_request)
    println("File $new_feature_request already exists, skipping")
  else
    println("Copying $template_feature_request to $new_feature_request")
    cp(template_feature_request, new_feature_request)

    println("Replace $template_package_name with $package_name in $new_feature_request")
    replace_in_file(new_feature_request, template_package_name => package_name)

    mv(new_feature_request, joinpath("..", new_feature_request); force=true)
  end
end


text:
@testset "Constant variables: $solver" for solver in solvers

    @testset "Issue #166" begin
        # Issue #166
        α = Variable(5)
        fix!(α, ones(5,1))

        # has const vexity, but not at the head
        c = (rand(5,5) * α) * ones(1,5)

        β = Variable(5)
        β.value = ones(5)

        problem = minimize(sum(c * β), [β >= 0])
        solve!(problem, solver)
        @test problem.optval ≈ evaluate(sum(c * β)) atol=TOL
        @test problem.optval ≈ 0.0 atol=TOL
        @test β.value ≈ zeros(5) atol=TOL
    end

    @testset "Issue #228" begin
        x = Variable(2)
        y = Variable(2)
        fix!(x, [1 1]')
        prob = minimize(y'*(x+[2 2]'), [y>=0])
        solve!(prob, solver)
        @test prob.optval ≈ 0.0 atol = TOL

        prob = minimize(x'*y, [y>=0])
        solve!(prob, solver)
        @test prob.optval ≈ 0.0 atol = TOL
    end

    @testset "Test double `fix!`" begin
        x = Variable()
        y = Variable()
        fix!(x, 1.0)
        prob = minimize(y*x, [y >= x, x >= 0.5])
        solve!(prob, solver)
        @test prob.optval ≈ 1.0 atol = TOL

        fix!(x, 2.0)
        solve!(prob, solver)
        @test prob.optval ≈ 4.0 atol = TOL

        free!(x)
        fix!(y, 1.0)
        solve!(prob, solver)
        @test prob.optval ≈ 0.5 atol = TOL
    end

    @testset "fix! and multiply" begin
        p = Variable()
        fix!(p, 1.0)
        x = Variable(2,2)
        prob = minimize( tr(p*x), [ x >= 1 ])
        solve!(prob, solver)
        @test prob.optval ≈ 2.0 atol = TOL
        @test evaluate( tr(p*x) ) ≈ 2.0 atol = TOL
    end

    @testset "fix! with complex numbers" begin
        x = ComplexVariable()
        fix!(x, 1.0 + im*1.0)
        y = Variable()
        prob = minimize( real(x*y), [ y >= .5, real(x) >= .5, imag(x) >= 0])
        solve!(prob, solver)
        @test prob.optval ≈ .5 atol=TOL
        @test evaluate(real(x*y)) ≈ .5 atol=TOL
        @test evaluate(y) ≈ 0.5 atol=TOL

        free!(x)
        fix!(y)
        solve!(prob, solver)
        @test prob.optval ≈ 0.25 atol=TOL
        @test evaluate(real(x*y)) ≈ 0.25 atol=TOL
        @test real(evaluate(x)) ≈ 0.5 atol=TOL
        @test evaluate(y) ≈ 0.5 atol=TOL

        @test_throws DimensionMismatch fix!(x, rand(2))
        @test_throws DimensionMismatch fix!(x, rand(2,2))
    end

    @testset "fix! with vectors" begin
        x = ComplexVariable(5)
        fix!(x, ones(5) + im*ones(5))
        y = Variable()
        prob = minimize( real(y*sum(x)), [ y >= .5, real(x) >= .5, imag(x) >= 0])
        solve!(prob, solver)
        @test prob.optval ≈ 2.5 atol=TOL
        @test evaluate(real(y*sum(x))) ≈ 2.5 atol=TOL
        @test evaluate(y) ≈ 0.5 atol=TOL

        free!(x)
        fix!(y)
        solve!(prob, solver)
        @test prob.optval ≈ 1.25 atol=TOL
        @test evaluate(real(y*sum(x))) ≈ 1.25 atol=TOL
        @test real(evaluate(x)) ≈ 0.5*ones(5) atol=TOL
        @test evaluate(y) ≈ 0.5 atol=TOL

        @test_throws DimensionMismatch fix!(x, rand(5,5))
        @test_throws DimensionMismatch fix!(x, rand(4))

    end

end


text:
using MachineLearning
using Test

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


text:
#=
fld-write.jl
Jeff Fessler, University of Michigan
=#

export fld_write

import FileIO # temporary

@deprecate fld_write(args... ; kwargs...) FileIO.save(args... ; kwargs...)


text:
# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. The
# `stepper` keyword defines the time-stepper to be used.
prob = TwoDNavierStokes.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, ν=ν, dt=dt, calcF=calcF!, stepper="ETDRK4")

# Next we define some shortcuts for convenience.
sol, cl, vs, gr = prob.sol, prob.clock, prob.vars, prob.grid
x, y = gr.x, gr.y

TwoDNavierStokes.set_ζ!(prob, ζ₀)

energy_initial = energy(prob)

U = sqrt(2 * energy_initial)

Re = U * Ly / ν

# ## Diagnostics

"""
    vorticityL4(prob)
Returns the domain-averaged enstrophy, ∫ ζ⁴ dxdy / (Lx Ly), for the solution in `sol`.
"""
@inline function vorticityL4(prob)
  sol, vars, grid = prob.sol, prob.vars, prob.grid
  @. vars.ζh = sol
  ldiv!(vars.ζ, grid.rfftplan, vars.ζh)
  return sum(@. vars.ζ^4) * grid.dx * grid.dy / (grid.Lx * grid.Ly)
end

"""
    palinstrophy(prob)
Returns the domain-averaged palinstrophy, ∫ |∇ζ|² dxdy / (Lx Ly), for the solution in `sol`.
"""
@inline function palinstrophy(prob)
  sol, vars, grid = prob.sol, prob.vars, prob.grid
  palinstrophyh = vars.uh # use vars.uh as scratch variable

  @. palinstrophyh = grid.Krsq * abs2(sol)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(palinstrophyh, grid)
end

# Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top.
E = Diagnostic(energy, prob; nsteps=nsteps)
Z2 = Diagnostic(enstrophy, prob; nsteps=nsteps)
Z4 = Diagnostic(vorticityL4, prob; nsteps=nsteps)
P = Diagnostic(palinstrophy, prob; nsteps=nsteps)
diags = [E, Z2, Z4, P] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "./data/"
filename = joinpath(filepath, "euler_longtime.jld2")
filename_diags = joinpath(filepath, "euler_longtime_diags.jld2")

filename = FourierFlows.uniquepath(filename)
@info "Output will be saved at $filename."

filename_diags = FourierFlows.uniquepath(filename_diags)
@info "Diagnostics will be saved at $filename_diags."

# Do some basic file management
if isfile(filename); rm(filename); end
if isfile(filename_diags); rm(filename_diags); end

# And then create Output
get_sol(prob) = sol # extracts the Fourier-transformed solution
out = Output(prob, filename, (:zetah, get_sol))
saveproblem(out)

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

saveoutput(out) # save initial condition

file = jldopen(filename, "a+")
FourierFlows.savefield(file, "params/nsubs", nsubs)
FourierFlows.savefield(file, "params/k₀", k₀)
close(file)

@info "Starting simulation..."

startwalltime = time()

for j=0:Int(round(nsteps/nsubs))-1

  if j%(1000/nsubs)==0
    TwoDNavierStokes.updatevars!(prob)

    cfl = cl.dt * maximum([maximum(vs.u) / gr.dx, maximum(vs.v) / gr.dy])

    estimated_remaining_walltime = (time()-startwalltime)/60 / cl.step * (nsteps-cl.step)
    log = @sprintf("step: %04d, t: %d, tν : %.4f, E: %.3e, Zₗ₂: %.3e, Zₗ₄: %.3e, P: %.3e, cfl: %.4f, walltime: %.2f min, estimated remaining walltime: %.2f min",
      cl.step, cl.t, ν*k₀^2*cl.t, E.data[E.i], Z2.data[Z2.i], Z4.data[Z4.i], P.data[P.i], cfl, (time()-startwalltime)/60, estimated_remaining_walltime)
    println(log)
  end

  stepforward!(prob, diags, nsubs)
  dealias!(sol, gr)

  saveoutput(out)
end

@info @sprintf("Simulation finished after %.2f minutes.", (time()-startwalltime)/60)

savediagnostic(E, "energy", filename_diags)
savediagnostic(Z2, "enstrophyL2", filename_diags)
savediagnostic(Z4, "enstrophyL4", filename_diags)
savediagnostic(P, "palinstrophy", filename_diags)

@info "Run visualize_simulation.jl after prescribing the two filenames used for saving output at the top the visualize_simulation.jl script."

@info "Output for this simulation was saved at $filename and $filename_diags."

max_stars_repo_name: navidcy/TwoDTurbulence-EulerLimit

max_issues_repo_name: navidcy/TwoDTurbulence-EulerLimit

max_forks_repo_name: navidcy/TwoDTurbulence-EulerLimit


text:
#####
##### routines to compute warp-and-blend interpolation nodes for triangles, tets,
##### and pyramid elements. Wedges/hexes use tensor product nodes
#####

"""
    get_edge_list(elem::AbstractElemShape)

Returns list of edges for a specific element (elem = Tri(), Pyr(), or Tet()).
"""
get_edge_list(elem::AbstractElemShape)

get_edge_list(elem::Tri) = [1,2],[2,3],[3,1]
get_edge_list(elem::Tet) = [1,4],[4,3],[3,1],[1,2],[3,2],[4,2]
get_edge_list(elem::Pyr) = [1,2],[2,4],[3,4],[3,1],[1,5],[2,5],[3,5],[4,5]

get_vertices(elem::AbstractElemShape) = equi_nodes(elem,1)
get_vertex_fxns(elem::AbstractElemShape) =
    (rst...)->vandermonde(elem,1,rst...)/vandermonde(elem,1,equi_nodes(elem,1)...)

# assumes r1D in [-1,1], v1,v2 are vertices
function interp_1D_to_line(r1D,v1,v2)
    runit = @. (1+r1D)/2
    return map((a,b)->(b-a)*runit .+ a,v1,v2)
end

"""
    interp_1D_to_edges(elem::AbstractElemShape, r1D)

Interpolates points r1D to the edges of an element (elem = :Tri, :Pyr, or :Tet)
"""
function interp_1D_to_edges(elem::AbstractElemShape, r1D)
    v = get_vertices(elem)
    edges = get_edge_list(elem)
    edge_pts = ntuple(x->zeros(length(r1D),length(edges)),length(v))
    for (i,e) in enumerate(edges)
        edge_pts_i = interp_1D_to_line(r1D,getindex.(v,e[1]),getindex.(v,e[2]))
        setindex!.(edge_pts,edge_pts_i,:,i)
    end
    return vec.(edge_pts)
end

"""
    edge_basis(elem::AbstractElemShape, N, rst...)

returns generalized Vandermonde matrix evaluated using an edge basis.
"""
function edge_basis(elem::AbstractElemShape, N, rst...)
    vertices = get_vertices(elem)
    edges = get_edge_list(elem)
    vertex_functions = get_vertex_fxns(elem)
    return edge_basis(N, vertices, edges, (N,r)->basis(Line(),N,r), vertex_functions, rst...)
end

"""
    edge_basis(N, vertices, edges, basis, vertex_functions, rst...)

Computes edge basis given vertex functions and 1D basis.
"""
function edge_basis(N, vertices, edges, basis1D, vertex_functions, rst...)
    V1 = vertex_functions(rst...)
    if N < 2
        return V1
    end

    V = zeros(length(first(rst)),length(edges)*(N-1) + length(first(vertices)))
    V[:,1:size(V1,2)] .= V1 # initialize vertex functions

    id = size(V1,2)+1
    for e in edges
        # edge parametrization from paper:
        # "A Comparison of High-Order Interpolation Nodes for the Pyramid"
        r1D_edge = V1[:,e[1]]-V1[:,e[2]]
        V1D,_ = basis1D(N-2,r1D_edge)
        for i = 1:size(V1D,2)
            V[:,id] = V1D[:,i].*V1[:,e[1]].*V1[:,e[2]]
            id += 1
        end
    end
    return V
end

"""
    build_warped_nodes(elem::AbstractElemShape,N,r1D)

Computes degree N warp-and-blend interpolation nodes for elem = Tri(), Pyr(), or
Tet() based on the 1D node set "r1D". Returns a tuple "rst" containing arrays of
interpolation points.
"""
function build_warped_nodes(elem::AbstractElemShape,N,r1D)
    r1D_equi = equi_nodes(Line(),N)
    rst_edge_equi = interp_1D_to_edges(elem,r1D_equi)
    V_edge = edge_basis(elem,N,rst_edge_equi...)
    rst_edge = interp_1D_to_edges(elem,r1D)

    c = (x->V_edge\x).(rst_edge) # should be lsq solve

    rst_equi = equi_nodes(elem,N)
    return (x->edge_basis(elem,N,rst_equi...)*x).(c)
end

# function compute_lobatto_interp(elem::AbstractElemShape,N,r1D)
#     r1D_equi = equi_nodes(Line(),N)
#     rst_edge_equi = interp_1D_to_edges(elem,r1D_equi)
#     V_edge = edge_basis(elem,N,rst_edge_equi...)
#     rst_edge = interp_1D_to_edges(elem,r1D)
#
#     c = (x->V_edge\x).(rst_edge)
#     # @show ((x->V_edge*x).(c) .- rst_edge)
#     return c
# end


text:
function build_mn_data(file_name)
    mp_data = PowerModels.parse_file(file_name)
    mp_data1 = InfrastructureModels.replicate(mp_data, 2, Set{String}(["source_type", "name", "source_version", "per_unit"]))
    PowerModelsACDC.process_additional_data!(mp_data1)
    return mp_data1
end


text:
abstract type Language; end

# Languages currently supported by Snowball
#Danish, Dutch, English, Finnish, French, German, Hungarian, Italian, Norwegian,
# Portuguese, Romanian, Russian, Spanish, Swedish, Turkish

# These are ISO 639-2T alpha-3 and ISO 639-3 codes
isocode(lang::T) where {T<:Language} = isocode(T)
name(lang::T) where {T<:Language} = name(T)
english_name(lang::T) where {T<:Language} = english_name(T)

struct Esperanto <: Language; end;	english_name(::Type{Esperanto}) = "Esperanto";	name(::Type{Esperanto}) = "Esperanto";	isocode(::Type{Esperanto}) = "epo";
struct English <: Language; end;	english_name(::Type{English}) = "English";	name(::Type{English}) = "English";	isocode(::Type{English}) = "eng"
struct Russian <: Language; end;	english_name(::Type{Russian}) = "Russian";	name(::Type{Russian}) = "Русский";	isocode(::Type{Russian}) = "rus"
struct Mandarin <: Language; end;	english_name(::Type{Mandarin}) = "Mandarin";	name(::Type{Mandarin}) = "官话";	isocode(::Type{Mandarin}) = "cmn"
struct Spanish <: Language; end;	english_name(::Type{Spanish}) = "Spanish";	name(::Type{Spanish}) = "Español";	isocode(::Type{Spanish}) = "spa"
struct Portuguese <: Language; end;	english_name(::Type{Portuguese}) = "Portuguese";	name(::Type{Portuguese}) = "Português";	isocode(::Type{Portuguese}) = "por"
struct Italian <: Language; end;	english_name(::Type{Italian}) = "Italian";	name(::Type{Italian}) = "Italiano";	isocode(::Type{Italian}) = "ita"
struct Bengali <: Language; end;	english_name(::Type{Bengali}) = "Bengali";	name(::Type{Bengali}) = "বাংলা";	isocode(::Type{Bengali}) = "ben"
struct French <: Language; end;	english_name(::Type{French}) = "French";	name(::Type{French}) = "Français";	isocode(::Type{French}) = "fra"
struct German <: Language; end;	english_name(::Type{German}) = "German";	name(::Type{German}) = "Deutsch";	isocode(::Type{German}) = "deu"
struct Ukrainian <: Language; end;	english_name(::Type{Ukrainian}) = "Ukrainian";	name(::Type{Ukrainian}) = "Українська";	isocode(::Type{Ukrainian}) = "ukr"
struct Georgian <: Language; end;	english_name(::Type{Georgian}) = "Georgian";	name(::Type{Georgian}) = "ქართული";	isocode(::Type{Georgian}) = "kat"
struct Arabic <: Language; end;	english_name(::Type{Arabic}) = "Arabic";	name(::Type{Arabic}) = "العربية";	isocode(::Type{Arabic}) = "arb"
struct Hindi <: Language; end;	english_name(::Type{Hindi}) = "Hindi";	name(::Type{Hindi}) = "हिन्दी";	isocode(::Type{Hindi}) = "hin"
struct Japanese <: Language; end;	english_name(::Type{Japanese}) = "Japanese";	name(::Type{Japanese}) = "日本語";	isocode(::Type{Japanese}) = "jpn"
struct Hebrew <: Language; end;	english_name(::Type{Hebrew}) = "Hebrew";	name(::Type{Hebrew}) = "עברית";	isocode(::Type{Hebrew}) = "heb"
struct Yiddish <: Language; end;	english_name(::Type{Yiddish}) = "Yiddish";	name(::Type{Yiddish}) = "ייִדיש";	isocode(::Type{Yiddish}) = "ydd"
struct Polish <: Language; end;	english_name(::Type{Polish}) = "Polish";	name(::Type{Polish}) = "Polski";	isocode(::Type{Polish}) = "pol"
struct Amharic <: Language; end;	english_name(::Type{Amharic}) = "Amharic";	name(::Type{Amharic}) = "አማርኛ";	isocode(::Type{Amharic}) = "amh"
struct Tigrinya <: Language; end;	english_name(::Type{Tigrinya}) = "Tigrinya";	name(::Type{Tigrinya}) = "ትግርኛ";	isocode(::Type{Tigrinya}) = "tir"
struct Javanese <: Language; end;	english_name(::Type{Javanese}) = "Javanese";	name(::Type{Javanese}) = "Basa Jawa";	isocode(::Type{Javanese}) = "jav"
struct Korean <: Language; end;	english_name(::Type{Korean}) = "Korean";	name(::Type{Korean}) = "한국어";	isocode(::Type{Korean}) = "kor"
struct Bokmal <: Language; end;	english_name(::Type{Bokmal}) = "Bokmal";	name(::Type{Bokmal}) = "Bokmål";	isocode(::Type{Bokmal}) = "nob"
struct Nynorsk <: Language; end;	english_name(::Type{Nynorsk}) = "Nynorsk";	name(::Type{Nynorsk}) = "Nynorsk";	isocode(::Type{Nynorsk}) = "nno"
struct Danish <: Language; end;	english_name(::Type{Danish}) = "Danish";	name(::Type{Danish}) = "Dansk";	isocode(::Type{Danish}) = "dan"
struct Swedish <: Language; end;	english_name(::Type{Swedish}) = "Swedish";	name(::Type{Swedish}) = "Svenska";	isocode(::Type{Swedish}) = "swe"
struct Finnish <: Language; end;	english_name(::Type{Finnish}) = "Finnish";	name(::Type{Finnish}) = "Suomi";	isocode(::Type{Finnish}) = "fin"
struct Turkish <: Language; end;	english_name(::Type{Turkish}) = "Turkish";	name(::Type{Turkish}) = "Türkçe";	isocode(::Type{Turkish}) = "tur"
struct Dutch <: Language; end;	english_name(::Type{Dutch}) = "Dutch";	name(::Type{Dutch}) = "Nederlands";	isocode(::Type{Dutch}) = "nld"
struct Hungarian <: Language; end;	english_name(::Type{Hungarian}) = "Hungarian";	name(::Type{Hungarian}) = "Magyar";	isocode(::Type{Hungarian}) = "hun"
struct Czech <: Language; end;	english_name(::Type{Czech}) = "Czech";	name(::Type{Czech}) = "Čeština";	isocode(::Type{Czech}) = "ces"
struct Greek <: Language; end;	english_name(::Type{Greek}) = "Greek";	name(::Type{Greek}) = "Ελληνικά";	isocode(::Type{Greek}) = "ell"
struct Bulgarian <: Language; end;	english_name(::Type{Bulgarian}) = "Bulgarian";	name(::Type{Bulgarian}) = "Български";	isocode(::Type{Bulgarian}) = "bul"
struct Belarusian <: Language; end;	english_name(::Type{Belarusian}) = "Belarusian";	name(::Type{Belarusian}) = "Беларуская";	isocode(::Type{Belarusian}) = "bel"
struct Marathi <: Language; end;	english_name(::Type{Marathi}) = "Marathi";	name(::Type{Marathi}) = "मराठी";	isocode(::Type{Marathi}) = "mar"
struct Kannada <: Language; end;	english_name(::Type{Kannada}) = "Kannada";	name(::Type{Kannada}) = "ಕನ್ನಡ";	isocode(::Type{Kannada}) = "kan"
struct Romanian <: Language; end;	english_name(::Type{Romanian}) = "Romanian";	name(::Type{Romanian}) = "Română";	isocode(::Type{Romanian}) = "ron"
struct Slovene <: Language; end;	english_name(::Type{Slovene}) = "Slovene";	name(::Type{Slovene}) = "Slovenščina";	isocode(::Type{Slovene}) = "slv"
struct Croatian <: Language; end;	english_name(::Type{Croatian}) = "Croatian";	name(::Type{Croatian}) = "Hrvatski";	isocode(::Type{Croatian}) = "hrv"
struct Serbian <: Language; end;	english_name(::Type{Serbian}) = "Serbian";	name(::Type{Serbian}) = "Српски";	isocode(::Type{Serbian}) = "srp"
struct Macedonian <: Language; end;	english_name(::Type{Macedonian}) = "Macedonian";	name(::Type{Macedonian}) = "Македонски";	isocode(::Type{Macedonian}) = "mkd"
struct Lithuanian <: Language; end;	english_name(::Type{Lithuanian}) = "Lithuanian";	name(::Type{Lithuanian}) = "Lietuvių";	isocode(::Type{Lithuanian}) = "lit"
struct Latvian <: Language; end;	english_name(::Type{Latvian}) = "Latvian";	name(::Type{Latvian}) = "Latviešu";	isocode(::Type{Latvian}) = "lav"
struct Estonian <: Language; end;	english_name(::Type{Estonian}) = "Estonian";	name(::Type{Estonian}) = "Eesti";	isocode(::Type{Estonian}) = "est"
struct Tamil <: Language; end;	english_name(::Type{Tamil}) = "Tamil";	name(::Type{Tamil}) = "தமிழ்";	isocode(::Type{Tamil}) = "tam"
struct Vietnamese <: Language; end;	english_name(::Type{Vietnamese}) = "Vietnamese";	name(::Type{Vietnamese}) = "Tiếng Việt";	isocode(::Type{Vietnamese}) = "vie"
struct Urdu <: Language; end;	english_name(::Type{Urdu}) = "Urdu";	name(::Type{Urdu}) = "اُردُو";	isocode(::Type{Urdu}) = "urd"
struct Thai <: Language; end;	english_name(::Type{Thai}) = "Thai";	name(::Type{Thai}) = "ภาษาไทย";	isocode(::Type{Thai}) = "tha"
struct Gujarati <: Language; end;	english_name(::Type{Gujarati}) = "Gujarati";	name(::Type{Gujarati}) = "ગુજરાતી";	isocode(::Type{Gujarati}) = "guj"
struct Uzbek <: Language; end;	english_name(::Type{Uzbek}) = "Uzbek";	name(::Type{Uzbek}) = "Oʻzbekcha";	isocode(::Type{Uzbek}) = "uzb"
struct Punjabi <: Language; end;	english_name(::Type{Punjabi}) = "Punjabi";	name(::Type{Punjabi}) = "ਪੰਜਾਬੀ";	isocode(::Type{Punjabi}) = "pan"
struct Azerbaijani <: Language; end;	english_name(::Type{Azerbaijani}) = "Azerbaijani";	name(::Type{Azerbaijani}) = "Azərbaycanca";	isocode(::Type{Azerbaijani}) = "azj"
struct Indonesian <: Language; end;	english_name(::Type{Indonesian}) = "Indonesian";	name(::Type{Indonesian}) = "Bahasa Indonesia";	isocode(::Type{Indonesian}) = "ind"
struct Telugu <: Language; end;	english_name(::Type{Telugu}) = "Telugu";	name(::Type{Telugu}) = "తెలుగు";	isocode(::Type{Telugu}) = "tel"
struct Persian <: Language; end;	english_name(::Type{Persian}) = "Persian";	name(::Type{Persian}) = "فارسی";	isocode(::Type{Persian}) = "pes"
struct Malayalam <: Language; end;	english_name(::Type{Malayalam}) = "Malayalam";	name(::Type{Malayalam}) = "മലയാളം";	isocode(::Type{Malayalam}) = "mal"
struct Hausa <: Language; end;	english_name(::Type{Hausa}) = "Hausa";	name(::Type{Hausa}) = "Hausa";	isocode(::Type{Hausa}) = "hau"
struct Oriya <: Language; end;	english_name(::Type{Oriya}) = "Oriya";	name(::Type{Oriya}) = "ଓଡ଼ିଆ";	isocode(::Type{Oriya}) = "ori"
struct Burmese <: Language; end;	english_name(::Type{Burmese}) = "Burmese";	name(::Type{Burmese}) = "မြန်မာစာ";	isocode(::Type{Burmese}) = "mya"
struct Bhojpuri <: Language; end;	english_name(::Type{Bhojpuri}) = "Bhojpuri";	name(::Type{Bhojpuri}) = "भोजपुरी";	isocode(::Type{Bhojpuri}) = "bho"
struct Tagalog <: Language; end;	english_name(::Type{Tagalog}) = "Tagalog";	name(::Type{Tagalog}) = "Tagalog";	isocode(::Type{Tagalog}) = "tgl"
struct Yoruba <: Language; end;	english_name(::Type{Yoruba}) = "Yoruba";	name(::Type{Yoruba}) = "Yorùbá";	isocode(::Type{Yoruba}) = "yor"
struct Maithili <: Language; end;	english_name(::Type{Maithili}) = "Maithili";	name(::Type{Maithili}) = "मैथिली";	isocode(::Type{Maithili}) = "mai"
struct Oromo <: Language; end;	english_name(::Type{Oromo}) = "Oromo";	name(::Type{Oromo}) = "Oromoo";	isocode(::Type{Oromo}) = "orm"
struct Igbo <: Language; end;	english_name(::Type{Igbo}) = "Igbo";	name(::Type{Igbo}) = "Igbo";	isocode(::Type{Igbo}) = "ibo"
struct Cebuano <: Language; end;	english_name(::Type{Cebuano}) = "Cebuano";	name(::Type{Cebuano}) = "Cebuano";	isocode(::Type{Cebuano}) = "ceb"
struct Kurdish <: Language; end;	english_name(::Type{Kurdish}) = "Kurdish";	name(::Type{Kurdish}) = "Kurdî";	isocode(::Type{Kurdish}) = "kur"
struct Malagasy <: Language; end;	english_name(::Type{Malagasy}) = "Malagasy";	name(::Type{Malagasy}) = "Malagasy";	isocode(::Type{Malagasy}) = "mlg"
struct Saraiki <: Language; end;	english_name(::Type{Saraiki}) = "Saraiki";	name(::Type{Saraiki}) = "سرائیکی";	isocode(::Type{Saraiki}) = "skr"
struct Nepali <: Language; end;	english_name(::Type{Nepali}) = "Nepali";	name(::Type{Nepali}) = "नेपाली";	isocode(::Type{Nepali}) = "nep"
struct Sinhalese <: Language; end;	english_name(::Type{Sinhalese}) = "Sinhalese";	name(::Type{Sinhalese}) = "සිංහල";	isocode(::Type{Sinhalese}) = "sin"
struct Khmer <: Language; end;	english_name(::Type{Khmer}) = "Khmer";	name(::Type{Khmer}) = "ភាសាខ្មែរ";	isocode(::Type{Khmer}) = "khm"
struct Turkmen <: Language; end;	english_name(::Type{Turkmen}) = "Turkmen";	name(::Type{Turkmen}) = "Türkmençe";	isocode(::Type{Turkmen}) = "tuk"
struct Somali <: Language; end;	english_name(::Type{Somali}) = "Somali";	name(::Type{Somali}) = "Soomaaliga";	isocode(::Type{Somali}) = "som"
struct Chewa <: Language; end;	english_name(::Type{Chewa}) = "Chewa";	name(::Type{Chewa}) = "Chichewa";	isocode(::Type{Chewa}) = "nya"
struct Akan <: Language; end;	english_name(::Type{Akan}) = "Akan";	name(::Type{Akan}) = "Akan";	isocode(::Type{Akan}) = "aka"
struct Zulu <: Language; end;	english_name(::Type{Zulu}) = "Zulu";	name(::Type{Zulu}) = "IsiZulu";	isocode(::Type{Zulu}) = "zul"
struct Kinyarwanda <: Language; end;	english_name(::Type{Kinyarwanda}) = "Kinyarwanda";	name(::Type{Kinyarwanda}) = "Kinyarwanda";	isocode(::Type{Kinyarwanda}) = "kin"
struct HaitianCreole <: Language; end;	english_name(::Type{HaitianCreole}) = "HaitianCreole";	name(::Type{HaitianCreole}) = "Kreyòl ayisyen";	isocode(::Type{HaitianCreole}) = "hat"
struct Ilocano <: Language; end;	english_name(::Type{Ilocano}) = "Ilocano";	name(::Type{Ilocano}) = "Ilokano";	isocode(::Type{Ilocano}) = "ilo"
struct Rundi <: Language; end;	english_name(::Type{Rundi}) = "Rundi";	name(::Type{Rundi}) = "Ikirundi";	isocode(::Type{Rundi}) = "run"
struct Shona <: Language; end;	english_name(::Type{Shona}) = "Shona";	name(::Type{Shona}) = "ChiShona";	isocode(::Type{Shona}) = "sna"
struct Uyghur <: Language; end;	english_name(::Type{Uyghur}) = "Uyghur";	name(::Type{Uyghur}) = "ئۇيغۇرچە";	isocode(::Type{Uyghur}) = "uig"
struct Swahili <: Language; end;	english_name(::Type{Swahili}) = "Swahili";	name(::Type{Swahili}) = "Kiswahili";	isocode(::Type{Swahili}) = "swh"

global const code_to_lang = Dict{String, Language}(
    "epo"  => Esperanto(),
    "eng"  => English(),
    "rus"  => Russian(),
    "cmn"  => Mandarin(),
    "spa"  => Spanish(),
    "por"  => Portuguese(),
    "ita"  => Italian(),
    "ben"  => Bengali(),
    "fra"  => French(),
    "deu"  => German(),
    "ukr"  => Ukrainian(),
    "kat"  => Georgian(),
    "arb"  => Arabic(),
    "hin"  => Hindi(),
    "jpn"  => Japanese(),
    "heb"  => Hebrew(),
    "ydd"  => Yiddish(),
    "pol"  => Polish(),
    "amh"  => Amharic(),
    "tir"  => Tigrinya(),
    "jav"  => Javanese(),
    "kor"  => Korean(),
    "nob"  => Bokmal(),
    "nno"  => Nynorsk(),
    "dan"  => Danish(),
    "swe"  => Swedish(),
    "fin"  => Finnish(),
    "tur"  => Turkish(),
    "nld"  => Dutch(),
    "hun"  => Hungarian(),
    "ces"  => Czech(),
    "ell"  => Greek(),
    "bul"  => Bulgarian(),
    "bel"  => Belarusian(),
    "mar"  => Marathi(),
    "kan"  => Kannada(),
    "ron"  => Romanian(),
    "slv"  => Slovene(),
    "hrv"  => Croatian(),
    "srp"  => Serbian(),
    "mkd"  => Macedonian(),
    "lit"  => Lithuanian(),
    "lav"  => Latvian(),
    "est"  => Estonian(),
    "tam"  => Tamil(),
    "vie"  => Vietnamese(),
    "urd"  => Urdu(),
    "tha"  => Thai(),
    "guj"  => Gujarati(),
    "uzb"  => Uzbek(),
    "pan"  => Punjabi(),
    "azj"  => Azerbaijani(),
    "ind"  => Indonesian(),
    "tel"  => Telugu(),
    "pes"  => Persian(),
    "mal"  => Malayalam(),
    "hau"  => Hausa(),
    "ori"  => Oriya(),
    "mya"  => Burmese(),
    "bho"  => Bhojpuri(),
    "tgl"  => Tagalog(),
    "yor"  => Yoruba(),
    "mai"  => Maithili(),
    "orm"  => Oromo(),
    "ibo"  => Igbo(),
    "ceb"  => Cebuano(),
    "kur"  => Kurdish(),
    "mlg"  => Malagasy(),
    "skr"  => Saraiki(),
    "nep"  => Nepali(),
    "sin"  => Sinhalese(),
    "khm"  => Khmer(),
    "tuk"  => Turkmen(),
    "som"  => Somali(),
    "nya"  => Chewa(),
    "aka"  => Akan(),
    "zul"  => Zulu(),
    "kin"  => Kinyarwanda(),
    "hat"  => HaitianCreole(),
    "ilo"  => Ilocano(),
    "run"  => Rundi(),
    "sna"  => Shona(),
    "uig"  => Uyghur(),
    )

function from_code(code::String)
    return get(code_to_lang, lowercase(code), nothing)
end

abstract type Script; end

struct ArabicScript <: Script; end
struct BengaliScript <: Script; end
struct CyrillicScript <: Script; end
struct DevanagariScript <: Script; end
struct EthiopicScript <: Script; end
struct GeorgianScript <: Script; end
struct GreekScript <: Script; end
struct GujaratiScript <: Script; end
struct GurmukhiScript <: Script; end
struct HangulScript <: Script; end
struct HebrewScript <: Script; end
struct HiraganaScript <: Script; end
struct KannadaScript <: Script; end
struct KatakanaScript <: Script; end
struct KhmerScript <: Script; end
struct LatinScript <: Script; end
struct MalayalamScript <: Script; end
struct MandarinScript <: Script; end
struct MyanmarScript <: Script; end
struct OriyaScript <: Script; end
struct SinhalaScript <: Script; end
struct TamilScript <: Script; end
struct TeluguScript <: Script; end
struct ThaiScript <: Script; end

name(::ArabicScript) = "Arabic"
name(::BengaliScript) = "Bengali"
name(::CyrillicScript) = "Cyrillic"
name(::DevanagariScript) = "Devanagari"
name(::EthiopicScript) = "Ethiopic"
name(::GeorgianScript) = "Georgian"
name(::GreekScript) = "Greek"
name(::GujaratiScript) = "Gujarati"
name(::GurmukhiScript) = "Gurmukhi"
name(::HangulScript) = "Hangul"
name(::HebrewScript) = "Hebrew"
name(::HiraganaScript) = "Hiragana"
name(::KannadaScript) = "Kannada"
name(::KatakanaScript) = "Katakana"
name(::KhmerScript) = "Khmer"
name(::LatinScript) = "Latin"
name(::MalayalamScript) = "Malayalam"
name(::MandarinScript) = "Mandarin"
name(::MyanmarScript) = "Myanmar"
name(::OriyaScript) = "Oriya"
name(::SinhalaScript) = "Sinhala"
name(::TamilScript) = "Tamil"
name(::TeluguScript) = "Telugu"
name(::ThaiScript) = "Thai"


text:
"""
    WaldA(;ν, k, A, θ)

Constructor for Wald distribution

# Fields

- `ν`: drift rate
- `k`: k = b - A where b is the decision threshold, and A is the maximum starting point
- `A`: the maximum starting point diffusion process, sampled from Uniform distribution
- `θ`: a encoding-motor time offset

## Usage

```julia
using SequentialSamplingModels
dist = WaldA(ν=.5, σ=1.0, ϕ=.3)
data = rand(dist, 10)
like = pdf.(dist, data)
loglike = logpdf.(dist, data)
```
# References

Tillman, G., Van Zandt, T., & Logan, G. D. (2020). Sequential sampling models without random between-trial variability:
The racing diffusion model of speeded decision making. Psychonomic Bulletin & Review, 27, 911-936.
"""
struct WaldA{T1,T2,T3,T4} <: ContinuousUnivariateDistribution
    ν::T1
    k::T2
    A::T3
    θ::T4
end

Broadcast.broadcastable(x::WaldA) = Ref(x)

WaldA(;ν, k, A, θ) = WaldA(ν, k, A, θ)

Φ(x) = cdf(Normal(0, 1), x)
ϕ(x) = pdf(Normal(0, 1), x)

function pdf(d::WaldA, rt::Float64)
    @unpack ν, k, A, θ = d
    t = rt - θ
    α = (k - t * ν) / √(t)
    β = (A + k - t * ν) / √(t)
    return (-ν * Φ(α) + 1 / √(t) * ϕ(α) +
        ν * Φ(β) - 1 / √(t) * ϕ(β)) / A
end

logpdf(d::WaldA, rt::Float64) = log(pdf(d, rt))

function cdf(d::WaldA, rt::Float64)
    @unpack ν, k, A, θ = d
    t = rt - θ
    b = A + k
    α1 = √(2) * ((ν * t - b)/√(2 * t))
    α2 = √(2) * ((ν * t - k)/√(2 * t))
    β1 = √(2) * ((-ν * t - b)/√(2 * t))
    β2 = √(2) * ((-ν * t - k)/√(2 * t))
    v1 = (1 / (2 * ν * A)) * (Φ(α2) - Φ(α1))
    v2 = (√(t) / A) * (α2 * Φ(α2) - α1 * Φ(α1))
    v3 = -(1 / (2 * ν * A)) * (exp(2 * ν * k) * Φ(β2) - exp(2 * ν * b) * Φ(β1))
    v4 = (√(t) / A) * (ϕ(α2) - ϕ(α1))
    return v1 + v2 + v3 + v4
end

logccdf(d::WaldA, rt::Float64) = log(1 - cdf(d, rt))

# function rand(d::WaldA, α)
#     @unpack ν, θ = d
#     return rand(InverseGaussian(α / ν, α^2)) + θ
# end

function rand(d::WaldA)
    @unpack ν, k, A, θ = d
    z = rand(Uniform(0, A))
    α = k + A - z
    return rand(InverseGaussian(α / ν, α^2)) + θ
end

"""
    DiffusionRace(;ν, k, A, θ)

An object for the racing diffusion model.

# Fields

- `ν`: a vector of drift rates
- `k`: k = b - A where b is the decision threshold, and A is the maximum starting point
- `A`: the maximum starting point diffusion process, sampled from Uniform distribution
- `θ`: a encoding-motor time offset

# Example

```julia
using SequentialSamplingModels
dist = DiffusionRace(;ν=[1,2], k=.3, A=.7, θ=.2)
data = rand(dist, 10)
like = pdf.(dist, data)
loglike = logpdf.(dist, data)
```
# References

Tillman, G., Van Zandt, T., & Logan, G. D. (2020). Sequential sampling models without random between-trial variability:
The racing diffusion model of speeded decision making. Psychonomic Bulletin & Review, 27, 911-936.
"""
struct DiffusionRace{T1,T2,T3,T4} <: SequentialSamplingModel
    ν::T1
    k::T2
    A::T3
    θ::T4
end

Broadcast.broadcastable(x::DiffusionRace) = Ref(x)

DiffusionRace(;ν, k, A, θ) = DiffusionRace(ν, k, A, θ)

function rand(dist::DiffusionRace)
    @unpack ν, A, k, θ = dist
    z = rand(Uniform(0, A))
    α = k + A - z
    x = @. rand(WaldA(ν, k, A, θ))
    rt,resp = findmin(x)
    return resp,rt
end

rand(dist::DiffusionRace, N::Int) = [rand(dist) for i in 1:N]

function logpdf(d::DiffusionRace, r::Int, rt::Float64)
    @unpack ν, k, A, θ = d
    LL = 0.0
    for (i,m) in enumerate(ν)
        if i == r
            LL += logpdf(WaldA(m, k, A, θ), rt)
        else
            LL += logccdf(WaldA(m, k, A, θ), rt)
        end
    end
    return LL
end

logpdf(d::DiffusionRace, data::Tuple) = logpdf(d, data...)

pdf(d::DiffusionRace, data::Tuple) = pdf(d, data...)

function pdf(d::DiffusionRace, r::Int, rt::Float64)
    @unpack ν, k, A, θ = d
    like = 1.0
    for (i,m) in enumerate(ν)
        if i == r
            like *= pdf(WaldA(m, k, A, θ), rt)
        else
            like *= (1 - cdf(WaldA(m, k, A, θ), rt))
        end
    end
    return like
end

# using Distributions, Parameters
# import Base.rand
# import Distributions: pdf, logpdf, cdf
# struct Race{T1,T2,T3} <: ContinuousUnivariateDistribution
#     ν::T1
#     α::T2
#     θ::T3
# end

# Broadcast.broadcastable(x::Race) = Ref(x)

# Race(;ν, α, θ) = Race(ν, α, θ)

# function rand(dist::Race)
#     @unpack ν, α, θ = dist
#     x = @. rand(Wald(ν, α, θ))
#     rt,resp = findmin(x)
#     return resp,rt
# end

# rand(dist::Race, N::Int) = [rand(dist) for i in 1:N]

# function logpdf(d::Race, r::Int, rt::Float64)
#     @unpack ν, α, θ = d
#     LL = 0.0
#     for (i,m) in enumerate(ν)
#         if i == r
#             LL += logpdf(Wald(;ν, α, θ), rt)
#         else
#             LL += logccdf(Wald(;ν, α, θ), rt)
#         end
#     end
#     return LL
# end

# logpdf(d::Race, data::Tuple) = logpdf(d, data...)

# pdf(d::Race, data::Tuple) = pdf(d, data...)

# function pdf(d::Race, r::Int, rt::Float64)
#     @unpack ν, α, θ = d
#     like = 1.0
#     for (i,m) in enumerate(ν)
#         if i == r
#             like *= pdf(Wald(;ν=m, α, θ), rt)
#         else
#             like *= (1 - cdf(Wald(;ν=m, α, θ), rt))
#         end
#     end
#     return like
# end


text:
font_awesome = joinpath(InteractBase.external_assets_path, "all.js")
prism_js = joinpath(InteractBase.external_assets_path, "prism.js")
prism_css = joinpath(InteractBase.external_assets_path, "prism.css")
highlight_css = joinpath(@__DIR__, "..", "assets", "highlight.css")
nouislider_min_js = joinpath(InteractBase.external_assets_path, "nouislider.min.js")
nouislider_min_css = joinpath(InteractBase.external_assets_path, "nouislider.min.css")
style_css = joinpath(@__DIR__, "..", "assets", "style.css")

@testset "deps" begin
    @test isfile(font_awesome)
    @test isfile(prism_js)
    @test isfile(prism_css)
    @test isfile(highlight_css)
    @test isfile(nouislider_min_js)
    @test isfile(nouislider_min_css)
    @test isfile(style_css)
end


text:
# ================================================
# Lorenz
#
# system type: continuous blackbox system
# state dimension: 3
#
# See https://en.wikipedia.org/wiki/Lorenz_system.
# ================================================
using MathematicalSystems, TaylorIntegration

@taylorize function lorenz!(dx, x, params, t)
    local σ = 10.0
    local β = 8.0 / 3.0
    local ρ = 28.0
    dx[1] = σ * (x[2] - x[1])
    dx[2] = x[1] * (ρ - x[3]) - x[2]
    dx[3] = x[1] * x[2] - β * x[3]
    return dx
end

function lorenz_model()
    # continuous blackbox system
    S = @system(x' = lorenz!(x), dim:3)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/UMAT.jl/blob/master/LICENSE

using UMAT, Test, DelimitedFiles, Tensors

dtime = 0.25
mat = GursonMaterial()
times = [mat.drivers.time]
stresses = [copy(tovoigt(mat.variables.stress))]
dstrain11 = 1e-3*dtime
dtimes = [dtime, dtime, dtime, dtime, 1.0]
dstrains11 = [dstrain11, dstrain11, dstrain11, -dstrain11, -4*dstrain11]
# stresses_expected = [[38.4091, 0.0, 0.0, 0.0, 0.0, 0.0],
#                      [100.0, 0.0, 0.0, 0.0, 0.0, 0.0],
#                      [100.0, 0.0, 0.0, 0.0, 0.0, 0.0],
#                      [50.0, 0.0, 0.0, 0.0, 0.0, 0.0],
#                      [-100.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
# strains_expected = [[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0.0, 0.0, 0.0],
#                     [2*dstrain11, -0.3*dstrain11*2, -0.3*dstrain11*2, 0.0, 0.0, 0.0],
#                     [3*dstrain11, -0.3*dstrain11*2 - 0.5*dstrain11, -0.3*dstrain11*2 - 0.5*dstrain11, 0.0, 0.0, 0.0],
#                     [2*dstrain11, -0.3*dstrain11 - 0.5*dstrain11, -0.3*dstrain11 - 0.5*dstrain11, 0.0, 0.0, 0.0],
#                     [-2*dstrain11, 0.3*dstrain11*2, 0.3*dstrain11*2, 0.0, 0.0, 0.0]]
#for i in 1:length(dtimes)
    #dstrain11 = dstrains11[i]
    #dtime = dtimes[i]
    uniaxial_increment!(mat, dstrain11, dtime)
    update_material!(mat)
    # push!(times, mat.drivers.time)
    # push!(stresses, copy(tovoigt(mat.variables.stress)))
    #@info(tovoigt(mat.variables.stress), stresses_expected[i])
    #@test isapprox(tovoigt(mat.variables.stress), stresses_expected[i])
    #@test isapprox(tovoigt(mat.drivers.strain; offdiagscale=2.0), strains_expected[i])
#end

# Define strain history
#e11 = 0.5.*vcat(Array(range(0, stop=1.5e-2, length=30)), Array(range(1.5e-2, stop=-1.5e-2, length=30)), Array(range(-1.5e-2, stop=3e-2, length=40)))
#e22 = -e11/2
#e33 = e22
#strains = [[e11[i], e22[i], e33[i], 0.0, 0.0, 0.0] for i in 1:100]

#s11 = [material.variables.STRESS[1]]
#s22 = [material.variables.STRESS[2]]
#s33 = [material.variables.STRESS[3]]

# dtime = 0.01
#
# for i=2:100
#     dstrain = strains[i]-strains[i-1]
#     material.ddrivers = UmatDriverState(NTENS=6, STRAN=dstrain)
#     UMAT.integrate_material!(material)
#     push!(s11, material.variables.STRESS[1])
#     push!(s22, material.variables.STRESS[2])
#     push!(s33, material.variables.STRESS[3])
#     Materials.update_material!(material)
# end
#
# for sig in ["s11", "s22", "s33"]
#     ref_stress = readdlm("test_gurson/ref_stresses_" * sig *".txt")
#     for (comp,ref) in zip(eval(Symbol(sig)), ref_stress)
#         @test isapprox(comp,ref,atol=sqrt(eps()))
#     end
# end


text:
import Base.copy
using Parameters

export
    UnconstrainedSolver,
    ConstrainedSolver,
    DirectSolver

export
    solver_name,
    cost,
    max_violation,
    options,
    stats,
    iterations,
    states,
    controls,
    initial_controls!,
    initial_states!,
    initial_trajectory!,
    rollout!,
    set_initial_state!,
    get_trajectory,
    get_model,
    get_times,
    initialize!

""" $(TYPEDEF)
Abstract solver for trajectory optimization problems

Any type that inherits from `AbstractSolver` must define the following methods:
```julia
model = get_model(::AbstractSolver)::AbstractModel
obj = get_objective(::AbstractSolver)::AbstractObjective
Z = get_trajectory(::AbstractSolver)::Traj
n,m,N = Base.size(::AbstractSolver)
x0 = get_initial_state(::AbstractSolver)::SVector
solve!(::AbstractSolver)
```
"""
abstract type AbstractSolver{T} <: MOI.AbstractNLPEvaluator end

"$(TYPEDEF) Unconstrained optimization solver. Will ignore
any constraints in the problem"
abstract type UnconstrainedSolver{T} <: AbstractSolver{T} end

"""$(TYPEDEF)
Abstract solver for constrained trajectory optimization problems

In addition to the methods required for `AbstractSolver`, all `ConstrainedSolver`s
    must define the following method
```julia
get_constraints(::ConstrainedSolver)::ConstrainSet
```
"""
abstract type ConstrainedSolver{T} <: AbstractSolver{T} end

""" $(TYPEDEF)
Solve the trajectory optimization problem by computing search directions using the joint
state vector, often solving the KKT system directly.
"""
abstract type DirectSolver{T} <: ConstrainedSolver{T} end

include("solvers/direct/dircol_ipopt.jl")
include("solvers/direct/dircol_snopt.jl")

@inline options(solver::AbstractSolver) = solver.opts
@inline stats(solver::AbstractSolver) = solver.stats
iterations(solver::AbstractSolver) = stats(solver).iterations

function cost(solver::AbstractSolver)
    obj = get_objective(solver)
    Z = get_trajectory(solver)
    cost(obj, Z)
end

function rollout!(solver::AbstractSolver)
    Z = get_trajectory(solver)
    model = get_model(solver)
    x0 = get_initial_state(solver)
    rollout!(model, Z, x0)
end

set_initial_state!(solver, x0) = copyto!(get_initial_state(solver), x0)

states(solver::AbstractSolver) = [state(z) for z in get_trajectory(solver)]

function controls(solver::AbstractSolver)
    N = size(solver)[3]
    Z = get_trajectory(solver)
    [control(Z[k]) for k = 1:N-1]
end

@inline initial_states!(solver::AbstractSolver, X0) = set_states!(get_trajectory(solver), X0)
@inline initial_controls!(solver::AbstractSolver, U0) = set_controls!(get_trajectory(solver), U0)
function initial_trajectory!(solver::AbstractSolver, Z0::Traj)
    Z = get_trajectory(solver)
    for k in eachindex(Z)
        Z[k].z = copy(Z0[k].z)
    end
end

@inline get_trajectory(solver::AbstractSolver) = solver.Z
@inline get_times(solver::AbstractSolver) = get_times(get_trajectory(solver))

# ConstrainedSolver methods
num_constraints(solver::AbstractSolver) = num_constraints(get_constraints(solver))

function max_violation(solver::ConstrainedSolver, Z::Traj)
    update_constraints!(solver, Z)
    max_violation(solver)
end

function max_violation(solver::ConstrainedSolver)
    conSet = get_constraints(solver)
    max_violation!(conSet)
    return maximum(conSet.c_max)
end

@inline findmax_violation(solver::ConstrainedSolver) =
    findmax_violation(get_constraints(solver))

""" $(SIGNATURES)
Calculate all the constraint values given the trajectory `Z`
"""
function update_constraints!(solver::ConstrainedSolver, Z::Traj=get_trajectory(solver))
    conSet = get_constraints(solver)
    evaluate!(conSet, Z)
end

function update_active_set!(solver::ConstrainedSolver, Z=get_trajectory(solver))
    conSet = get_constraints(solver)
    update_active_set!(conSet, Z, Val(solver.opts.active_set_tolerance))
end

""" $(SIGNATURES)
Calculate all the constraint Jacobians given the trajectory `Z`
"""
function constraint_jacobian!(solver::ConstrainedSolver, Z=get_trajectory(solver))
    conSet = get_constraints(solver)
    jacobian!(conSet, Z)
    return nothing
end

# Logging
function set_verbosity!(opts)
    log_level = opts.log_level
    if opts.verbose
        set_logger()
        Logging.disable_logging(LogLevel(log_level.level-1))
    else
        Logging.disable_logging(log_level)
    end
end

function clear_cache!(opts)
    if opts.verbose
        log_level = opts.log_level
        clear_cache!(global_logger().leveldata[log_level])
    end
end