text:
# FT and NaN test for the structs
@testset "Photosynthesis --- structures" begin
    for FT in [Float32, Float64]
        for data_set in [ Photosynthesis.KcTDBernacchi(FT),
                          Photosynthesis.VpmaxTDBoyd(FT),
                          C3CLM(FT),
                          C4CLM(FT),
                          AirLayer{FT}(),
                          Leaf{FT}(),
                          Photosynthesis.FluorescenceVanDerTol(FT),
                          Photosynthesis.FluorescenceVanDerTolDrought(FT),
                          Photosynthesis.KoTDBernacchi(FT),
                          Photosynthesis.RespirationTDBernacchi(FT),
                          Photosynthesis.VcmaxTDBernacchi(FT),
                          Photosynthesis.VomaxTDBernacchi(FT),
                          Photosynthesis.ΓStarTDBernacchi(FT),
                          Photosynthesis.KpepTDBoyd(FT),
                          Photosynthesis.JmaxTDLeuning(FT),
                          Photosynthesis.VcmaxTDLeuning(FT),
                          Photosynthesis.JmaxTDBernacchi(FT),
                          Photosynthesis.VtoRCollatz(FT),
                          C3Bernacchi(FT),
                          Photosynthesis.Q10TDAngiosperm(FT),
                          Photosynthesis.Q10TDGymnosperm(FT) ]
            @test FT_test(data_set, FT);
            @test NaN_test(data_set);
        end
    end
end

# FT and NaN test for the structs
println();
@testset "Photosynthesis --- functions" begin
    for FT in [Float32, Float64]
        c3_set   = C3CLM(FT);
        c4_set   = C4CLM(FT);
        leaf_3   = Leaf{FT}();
        leaf_4   = Leaf{FT}();
        td_q10   = Photosynthesis.Q10TD{FT}(1, 273.15, 1.7);
        envir    = AirLayer{FT}();
        fluo_set = c3_set.Flu;
        T        = rand(FT) + 298;
        glc      = FT(0.1);
        p_i      = rand(FT) + 20;

        # temperature corrections
        Photosynthesis.photo_TD_from_set(td_q10, T);
        leaf_temperature_dependence!(c3_set, leaf_3, envir, T);
        leaf_temperature_dependence!(c4_set, leaf_4, envir, T);

        # rubisco limited rates
        Photosynthesis.rubisco_limited_rate!(c3_set, leaf_3);
        Photosynthesis.rubisco_limited_rate!(c4_set, leaf_4);
        @test NaN_test(leaf_3);
        @test NaN_test(leaf_4);
        Photosynthesis.rubisco_limited_rate!(c3_set, leaf_3, envir);
        @test NaN_test(leaf_3);

        # light limited rates
        Photosynthesis.leaf_ETR!(c3_set, leaf_3);
        Photosynthesis.leaf_ETR!(c4_set, leaf_4);
        Photosynthesis.light_limited_rate!(c3_set, leaf_3);
        Photosynthesis.light_limited_rate!(c4_set, leaf_4);
        @test NaN_test(leaf_3);
        @test NaN_test(leaf_4);
        Photosynthesis.light_limited_rate!(c3_set, leaf_3, envir);
        @test NaN_test(leaf_3);

        # product limited rates
        Photosynthesis.product_limited_rate!(c3_set, leaf_3);
        Photosynthesis.product_limited_rate!(c4_set, leaf_4);
        @test NaN_test(leaf_3);
        @test NaN_test(leaf_4);
        Photosynthesis.product_limited_rate!(c4_set, leaf_4, envir);
        @test NaN_test(leaf_4);

        # fluorescence
        leaf_photosynthesis!(c3_set, leaf_3, envir, PCO₂Mode(), FT(2));
        leaf_fluorescence!(fluo_set, leaf_3);
        leaf_photosynthesis!(c3_set, leaf_3, envir, GCO₂Mode());
        leaf_fluorescence!(fluo_set, leaf_3);
        @test NaN_test(leaf_3);

        # leaf photo from glc
        leaf_photosynthesis!(c3_set, leaf_3, envir, GCO₂Mode(), glc);
        leaf_photosynthesis!(c4_set, leaf_4, envir, GCO₂Mode(), glc);
        @test NaN_test(leaf_3);
        @test NaN_test(leaf_4);

        # leaf photo from p_i
        leaf_photosynthesis!(c3_set, leaf_3, envir, PCO₂Mode(), p_i);
        leaf_photosynthesis!(c4_set, leaf_4, envir, PCO₂Mode(), p_i);
        @test NaN_test(leaf_3);
        @test NaN_test(leaf_4);
    end
end

max_stars_repo_name: climate-machine/Land


text:
abstract type AbstractReca <: AbstractEchoStateNetwork end

struct RECA_discrete <: AbstractReca
    res_size::Int
    in_size::Int
    out_size::Int
    rule::Int
    generations::Int
    expansion_size::Int
    permutations::Int
    train_data::AbstractArray{Float64}
    nla_type::ReservoirComputing.NonLinearAlgorithm
    states::AbstractArray{Float64}
    maps::AbstractArray{Int}
end

"""
    RECA_discrete(train_data::AbstractArray{Int}, rule::Int, generations::Int,
    expansion_size::Int, permutations::Int [, nla_type])

Create a RECA struct as described in [1] and [2].

[1] Yilmaz, Ozgur. “Reservoir computing using cellular automata.” arXiv preprint arXiv:1410.0162 (2014).
[2] Nichele, Stefano, and Andreas Molund. “Deep reservoir computing using cellular automata.” arXiv preprint arXiv:1703.02806 (2017).

"""
function RECA_discrete(train_data::AbstractArray{Int},
    rule::Int,
    generations::Int,
    expansion_size::Int,
    permutations::Int;
    nla_type = NLADefault())

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)
    res_size = generations*expansion_size*permutations

    maps = init_maps(in_size, permutations, expansion_size)
    states = reca_states_discrete(train_data, expansion_size, permutations, generations, rule, maps)
    train_data = convert(AbstractArray{Float64}, train_data)

    return RECA_discrete(res_size, in_size, out_size, rule, generations, expansion_size, permutations, train_data, nla_type, states, maps)
end

"""
    RECAdirect_predict_discrete(reca::AbstractReca, W_out::AbstractArray{Float64}, test_data::AbstractArray{Int})

Given the input data return the corresponding predicted output, as described in [1].

[1] Yilmaz, Ozgur. “Reservoir computing using cellular automata.” arXiv preprint arXiv:1410.0162 (2014).
"""
function RECAdirect_predict_discrete(reca::AbstractReca,
    W_out::AbstractArray{Float64},
    test_data::AbstractArray{Int})

    predict_len = size(test_data, 2)
    output = Array{Int}(undef, size(W_out, 1), predict_len)
    init_ca = zeros(Int, reca.expansion_size*reca.permutations)

    for i=1:predict_len
        init_ca = encoding(test_data[:,i], init_ca, reca.maps)
        ca = ECA(reca.rule, init_ca, reca.generations+1)
        ca_states = ca.cells[2:end ,:]
        x = copy(reshape(transpose(ca_states), reca.generations*reca.expansion_size*reca.permutations))
        out = W_out*x
        init_ca = ca.cells[end, :]
        output[:,i] = convert(AbstractArray{Int}, out .> 0.5)
    end
    return output
end

function reca_states_discrete(input_data::AbstractArray{Int},
    expansion_size::Int,
    permutations::Int,
    generations::Int,
    rule::Int,
    maps::AbstractArray{Int})

    input_size = size(input_data, 1)
    train_time = size(input_data, 2)

    states = zeros(Int, generations*expansion_size*permutations, train_time)
    init_ca = zeros(Int, expansion_size*permutations)

    for i=1:train_time
        init_ca = encoding(input_data[:,i], init_ca, maps)
        ca = ECA(rule, init_ca, generations+1)
        ca_states = ca.cells[2:end ,:]
        states[:,i] = reshape(transpose(ca_states), generations*expansion_size*permutations)
        init_ca = ca.cells[end, :]
    end
    return convert(AbstractArray{Float64}, states)
end

function encoding(input_vector::AbstractArray{Int},
    tot_encoded_vector::AbstractArray{Int},
    maps::AbstractArray{Int})

    permutations = size(maps, 1)
    input_size = size(input_vector,1)
    single_encoded_size = Int(size(tot_encoded_vector, 1)/permutations)

    new_tot_enc_vec = copy(tot_encoded_vector)

    for i=1:permutations
       new_tot_enc_vec[(i-1)*single_encoded_size+1:i*single_encoded_size]  = single_encoding(input_vector,
            new_tot_enc_vec[(i-1)*single_encoded_size+1:i*single_encoded_size],
            maps[i,:])
    end
    return new_tot_enc_vec
end

function single_encoding(input_vector::AbstractArray{Int},
    encoded_vector::AbstractArray{Int},
    map::AbstractArray{Int})

    new_enc_vec = copy(encoded_vector)

    for i=1:size(input_vector, 1)
        new_enc_vec[map[i]] = input_vector[i]
    end
    return new_enc_vec
end

function init_maps(input_size::Int,
        permutations::Int,
        mapped_vector_size::Int)

    maps = Array{Int}(undef, permutations, input_size)
    #tot_size = input_size*permutations

    for i=1:permutations
        maps[i,:] = mapping(input_size, mapped_vector_size)
    end
    return maps
end

function mapping(input_size::Int, mapped_vector_size::Int)

    map = sample(1:mapped_vector_size, input_size, replace=false)
    return map
end

ext: jl


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

Struct holding solution information for SparseSystem. Solution
is stored in a sparse matrix structure.

This class plays well with the abstract array interface.

Fields:

$(TYPEDFIELDS)
"""
struct SparseSolutionArray{Tv,Ti} <: AbstractMatrix{Tv}

    """
    Sparse matrix holding actual data.
    """
    node_dof::SparseMatrixCSC{Tv,Ti}
end

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

Return size of sparse solution array.
"""
Base.size(a::SparseSolutionArray)=size(a.node_dof)

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

Array of values in sparse solution array.
"""
values(a::SparseSolutionArray)=a.node_dof.nzval

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

Create a copy of sparse solution array
"""
Base.copy(this::SparseSolutionArray{Tv,Ti}) where {Tv,Ti} = SparseSolutionArray{Tv,Ti}(SparseMatrixCSC(this.node_dof.m,
                                                                                this.node_dof.n,
                                                                                this.node_dof.colptr,
                                                                                this.node_dof.rowval,
                                                                                Base.copy(this.node_dof.nzval)
                                                                                )
                                                                                )

"""
$(SIGNATURES)

Create a similar unintialized sparse solution array
"""
Base.similar(this::SparseSolutionArray{Tv,Ti}) where {Tv,Ti} = SparseSolutionArray{Tv,Ti}(SparseMatrixCSC(this.node_dof.m,
                                                                                this.node_dof.n,
                                                                                this.node_dof.colptr,
                                                                                this.node_dof.rowval,
                                                                                Base.similar(this.node_dof.nzval)
                                                                                )
                                                                                )
##################################################################
"""
$(SIGNATURES)

Get number of degree of freedom. Return 0 if species is not defined in node.
"""
function dof(a::SparseSolutionArray{Tv,Ti},i::Integer, j::Integer) where {Tv,Ti}
    A=a.node_dof
    coljfirstk = Int(A.colptr[j])
    coljlastk = Int(A.colptr[j+1] - 1)
    searchk = searchsortedfirst(A.rowval, i, coljfirstk, coljlastk, Base.Order.Forward)
    if searchk <= coljlastk && A.rowval[searchk] == i
        return searchk
    end
    return 0
end

struct SparseSolutionIndices
    a::SparseSolutionArray
end

unknown_indices(a::SparseSolutionArray) = SparseSolutionIndices(a)

Base.getindex(idx::SparseSolutionIndices,i,j)=dof(idx.a,i,j)

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

Set value for degree of freedom.
"""
function setdof!(a::SparseSolutionArray,v,i::Integer)
    a.node_dof.nzval[i] = v
end

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

Return  value for degree of freedom.
"""
getdof(a::SparseSolutionArray,i::Integer) =a.node_dof.nzval[i]

Base.:-(a::SparseSolutionArray,b::SparseSolutionArray)=SparseSolutionArray(a.node_dof-b.node_dof)
Base.:+(a::SparseSolutionArray,b::SparseSolutionArray)=SparseSolutionArray(a.node_dof+b.node_dof)

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

Accessor for sparse solution array.
"""
function Base.setindex!(a::SparseSolutionArray, v, ispec::Integer, inode::Integer)
    searchk=dof(a,ispec,inode)
    if searchk>0
        setdof!(a,v,searchk)
        return a
    end
    # TODO: what is the right reacton here ?
    # Ignoring seems to be better, so we can broacast etc.
    # throw(DomainError("undefined degree of freedom"))
end

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

Accessor for sparse solution array.
"""
function Base.getindex(a::SparseSolutionArray, ispec::Integer, inode::Integer)
    searchk=dof(a,ispec,inode)
    if searchk>0
        return getdof(a,searchk)
    end
    #
    # TODO: what is the right reacton here ?
    # Actually, NaN plays well with pyplot...
    return NaN
end

#
# Accessors for node-dof based loops
#
_firstnodedof(U::SparseSolutionArray,K) =U.node_dof.colptr[K]
_lastnodedof(U::SparseSolutionArray,K) =U.node_dof.colptr[K+1]-1
_spec(U::SparseSolutionArray,idof,K) =U.node_dof.rowval[idof]
_add(U::SparseSolutionArray,idof,val)=U.node_dof.nzval[idof]+=val


text:
# deprecated functions

function getVert{T <: AbstractString}(fgl::FactorGraph, lbl::T; api::DataLayerAPI=dlapi)
  warn("IncrementalInference.getVert{T <: AbstractString}(fgl::FactorGraph, lbl::T; api::DataLayerAPI=dlapi) is deprecated, use lbl::Symbol instead")
  getVert(fgl, Symbol(lbl), api=api)
end

function getVal{T <: AbstractString}(fgl::FactorGraph, lbl::T; api::DataLayerAPI=dlapi)
  warn("IncrementalInference.getVal{T <: AbstractString}(fgl::FactorGraph, lbl::T; api::DataLayerAPI=dlapi) is deprecated, use lbl::Symbol instead")
  getVal(fgl, Symbol(lbl),api=api)
end


text:
using NQCModels: ndofs

"""
    CompositeFrictionModel{M,F} <: AdiabaticFrictionModel

Combine an `AdiabaticModel` with an `ElectronicFrictionProvider`.

This allows for arbitrary composition of potentials and friction providers,
such that any adiabatic model can be augmented with any form of electronic friction.
"""
struct CompositeFrictionModel{M<:AdiabaticModel,F} <: AdiabaticFrictionModel
    pes_model::M
    friction_model::F
end

function NQCModels.ndofs(model::CompositeFrictionModel)
    model_dofs = ndofs(model.pes_model)
    friction_dofs = ndofs(model.friction_model)

    if model_dofs == friction_dofs
        return model_dofs
    else
        throw(ArgumentError("Dimensionality of PES and friction models does not match."))
    end
end

function NQCModels.potential(model::CompositeFrictionModel, R::AbstractMatrix)
    NQCModels.potential(model.pes_model, R)
end

function NQCModels.derivative!(model::CompositeFrictionModel, D, R::AbstractMatrix)
    NQCModels.derivative!(model.pes_model, D, R)
end

function friction!(model::CompositeFrictionModel, F, R)
    friction!(model.friction_model, F, R)
end


text:
SUITE["parse"] = BenchmarkGroup()

# https://github.com/JuliaTime/TimeZones.jl/issues/25
function parse_dates(n)
    df = DateFormat("yyyymmddHH:MM:SS ZZZ")
    arr = Vector{ZonedDateTime}(undef, n)
    for (i, s) in enumerate(Iterators.repeated("2016060701:02:03 America/Toronto", n))
        arr[i] = ZonedDateTime(s,df)
    end
    return arr
end

SUITE["parse"]["issue#25"] = @benchmarkable parse_dates(1000)

SUITE["tryparsenext_fixedtz"] = BenchmarkGroup()
fixed_tz_strs = [
    "Z",
    "-06:00",
    "+06:00",
    "-0600",
    "+0600",
    "-06",
    "+06",
    "UTC",  # Invalid
]

for fixed_tz in fixed_tz_strs
    str = "2000-01-02T03:04:05" * fixed_tz
    SUITE["tryparsenext_fixedtz"][fixed_tz] = begin
        @benchmarkable TimeZones.tryparsenext_fixedtz($str, 20, $(lastindex(str)))
    end
end

SUITE["tryparsenext_tz"] = BenchmarkGroup()
tz_strs = [
    "UTC",
    "GMT",
    "EST",  # Invalid
    "America/Winnipeg",
    "America/Argentina/ComodRivadavia",
]

for tz in tz_strs
    str = "2000-01-02T03:04:05" * tz
    SUITE["tryparsenext_tz"][tz] = begin
        @benchmarkable TimeZones.tryparsenext_tz($str, 20, $(lastindex(str)))
    end
end

SUITE["parse"]["ISOZonedDateTimeFormat"] = begin
    str = "2000-01-02T03:04:05.006+0700"
    df = TimeZones.ISOZonedDateTimeFormat
    @benchmarkable parse(ZonedDateTime, $str, $df)
end

SUITE["parse"]["generic-string"] = begin
    str = GenericString("2018-01-01 00:00 UTC")
    df = dateformat"yyyy-mm-dd HH:MM ZZZ"
    @benchmarkable parse(ZonedDateTime, $str, $df)
end


text:
# Deterministic Finite Automaton
# ==============================

struct DFANode
    edges::Vector{Tuple{Edge,DFANode}}
    final::Bool
    eof_actions::ActionList
    nfanodes::Set{NFANode}
end

function DFANode(final::Bool, eof_actions::ActionList, nodes::Set{NFANode})
    return DFANode(Tuple{Edge,DFANode}[], final, eof_actions, nodes)
end

function DFANode(final::Bool, nodes::Set{NFANode})
    return DFANode(final, ActionList(), nodes)
end

function Base.show(io::IO, node::DFANode)
    print(io, summary(node), "(<", length(node.edges), " edges", '@', objectid(node), ">)")
end

struct DFA
    start::DFANode
end

# Check if the DFA is really deterministic or not.
function validate(dfa::DFA)
    is_non_deterministic(e1, e2) = !(isdisjoint(e1.labels, e2.labels) || conflicts(e1.precond, e2.precond))
    for s in traverse(dfa.start)
        for i in 1:lastindex(s.edges), j in 1:i-1
            ei = s.edges[i][1]
            ej = s.edges[j][1]
            if overlaps(ei, ej)
                error("found non-deterministic edges")
            end
        end
    end
end

function nfa2dfa(nfa::NFA)
    newnodes = Dict{Set{NFANode},DFANode}()
    new(S) = get!(newnodes, S, DFANode(nfa.final ∈ S, S))
    isvisited(S) = haskey(newnodes, S)
    S = epsilon_closure(nfa.start)
    start = new(S)
    unvisited = [S]
    while !isempty(unvisited)
        # TODO: support fail
        S = pop!(unvisited)
        S_actions = accumulate_actions(S)
        s′ = new(S)

        if s′.final
            union!(s′.eof_actions, S_actions[nfa.final])
        end

        # accumulate edges and preconditions
        labels = Vector{ByteSet}()
        preconds = Set{Symbol}()
        for s in S, (e, t) in s.edges
            if !iseps(e)
                push!(labels, e.labels)
                union!(preconds, precondition_names(e.precond))
            end
        end

        # append DFA edges and nodes
        pn = collect(preconds)
        for label in disjoint_split(labels)
            # This enumeration will not finish in reasonable time when there
            # are too many preconditions.
            edges = Dict{Tuple{DFANode,ActionList},Vector{UInt64}}()
            for pv in UInt64(0):UInt64((1 << length(pn)) - 1)
                T = Set{NFANode}()
                actions = ActionList()
                for s in S, (e, t) in s.edges
                    if !isdisjoint(e.labels, label) && satisfies(e, pn, pv)
                        push!(T, t)
                        union!(actions, e.actions)
                        union!(actions, S_actions[s])
                    end
                end
                if !isempty(T)
                    T = epsilon_closure(T)
                    if !isvisited(T)
                        push!(unvisited, T)
                    end
                    push!(get!(edges, (new(T), actions), UInt64[]), pv)
                end
            end
            for ((t′, actions), pvs) in edges
                pn′, pvs′ = remove_redundant_preconds(pn, pvs)
                for pv′ in pvs′
                    push!(s′.edges, (Edge(label, make_precond(pn′, pv′), actions), t′))
                end
            end
        end
    end
    return DFA(start)
end

function epsilon_closure(node::NFANode)
    return epsilon_closure(Set([node]))
end

function epsilon_closure(nodes::Set{NFANode})
    closure = Set{NFANode}()
    unvisited = copy(nodes)
    while !isempty(unvisited)
        s = pop!(unvisited)
        push!(closure, s)
        for (e, t) in s.edges
            if iseps(e) && t ∉ closure
                push!(unvisited, t)
            end
        end
    end
    return closure
end

function disjoint_split(sets::Vector{ByteSet})
    # TODO: maybe too slow when length(sets) is large
    cut(s1, s2) = (intersect(s1, s2), setdiff(s1, s2))
    m = typemax(UInt64)
    disjsets = [ByteSet(m, m, m, m)]
    disjsets′ = ByteSet[]
    for x in sets
        for y in disjsets
            y1, y2 = cut(y, x)
            if !all(isdisjoint(z, y1) for z in sets)
                push!(disjsets′, y1)
            end
            if !all(isdisjoint(z, y2) for z in sets)
                push!(disjsets′, y2)
            end
        end
        disjsets, disjsets′ = disjsets′, disjsets
        empty!(disjsets′)
    end
    return disjsets
end

function accumulate_actions(S::Set{NFANode})
    top = copy(S)
    for s in S
        for (e, t) in s.edges
            if iseps(e)
                delete!(top, t)
            end
        end
    end
    @assert !isempty(top)
    actions = Dict(s => ActionList() for s in S)
    visited = Set{NFANode}()
    unvisited = top
    while !isempty(unvisited)
        s = pop!(unvisited)
        push!(visited, s)
        for (e, t) in s.edges
            if iseps(e)
                @assert !isconditioned(e.precond)
                union!(actions[t], e.actions)
                union!(actions[t], actions[s])
                if t ∉ visited
                    push!(unvisited, t)
                end
            end
        end
    end
    return actions
end

function satisfies(edge::Edge, names::Vector{Symbol}, pv::UInt64)
    for (n, v) in edge.precond
        i = findfirst(isequal(n), names)
        @assert i !== nothing
        @assert 0 < i ≤ 64
        vi = bitat(pv, i)
        if !(v == BOTH || (v == TRUE && vi) || (v == FALSE && !vi))
            return false
        end
    end
    return true
end

function remove_redundant_preconds(names::Vector{Symbol}, pvs::Vector{UInt64})
    mask(n) = ((1 << n) - 1) % UInt64
    newnames = Symbol[]
    pvs = copy(pvs)
    left = length(names)
    for name in names
        sort!(pvs)
        fnd = findfirst(pv -> bitat(pv, left), pvs)
        k = ifelse(fnd == nothing, 0, fnd) # TODO: See if there is a more elegant way of doing this.
        if (k - 1) * 2 == length(pvs)
            redundant = true
            for i in 1:k-1
                m = mask(left - 1)
                if pvs[i] & m != pvs[i+k-1] & m
                    redundant = false
                    break
                end
            end
        else
            redundant = false
        end
        if redundant
            left -= 1
            for i in 1:lastindex(pvs)
                # remove the redundant bit
                pvs[i] = pvs[i] & mask(left)
            end
        else
            push!(newnames, name)
            for i in 1:lastindex(pvs)
                # circular left shift
                pvs[i] = ((pvs[i] << 1) & mask(left)) | bitat(pvs[i], left)
            end
        end
    end
    return newnames, unique(pvs)
end

function make_precond(names::Vector{Symbol}, pv::UInt64)
    precond = Precondition()
    for (i, n) in enumerate(names)
        push!(precond, n => bitat(pv, i) ? TRUE : FALSE)
    end
    return precond
end

function bitat(x::UInt64, i::Integer)
    return ((x >> (i - 1)) & 1) == 1
end

function reduce_nodes(dfa::DFA)
    Q = Set(traverse(dfa.start))
    distinct = distinct_nodes(Q)
    newnodes = Dict{Set{DFANode},DFANode}()
    new(S) = get!(newnodes, S) do
        s = first(S)
        return DFANode(s.final, s.eof_actions, foldl((x, s) -> union(x, s.nfanodes), S, init=Set{NFANode}()))
    end
    equivalent(s) = filter(t -> (s, t) ∉ distinct, Q)
    isvisited(T) = haskey(newnodes, T)
    S = equivalent(dfa.start)
    start = new(S)
    unvisited = [S]
    while !isempty(unvisited)
        S = pop!(unvisited)
        s′ = new(S)
        for (e, t) in first(S).edges
            T = equivalent(t)
            if !isvisited(T)
                push!(unvisited, T)
            end
            push!(s′.edges, (e, new(T)))
        end
    end
    return DFA(start)
end

function distinct_nodes(S::Set{DFANode})
    labels = Dict(s => foldl((x, y) -> union(x, y[1].labels), s.edges, init=ByteSet()) for s in S)
    distinct = Set{Tuple{DFANode,DFANode}}()

    for s1 in S, s2 in S
        if s1.final != s2.final || labels[s1] != labels[s2] || s1.eof_actions != s2.eof_actions
            push!(distinct, (s1, s2))
        end
    end

    converged = false
    while !converged
        converged = true
        for s1 in S, s2 in S
            if s1 == s2 || (s1, s2) ∈ distinct
                continue
            end
            @assert labels[s1] == labels[s2] && s1.eof_actions == s2.eof_actions
            for (e1, t1) in s1.edges, (e2, t2) in s2.edges
                if overlaps(e1, e2) && ((t1, t2) ∈ distinct || e1.actions != e2.actions)
                    push!(distinct, (s1, s2), (s2, s1))
                    converged = false
                    break
                end
            end
        end
    end

    return distinct
end

function overlaps(e1::Edge, e2::Edge)
    return !(isdisjoint(e1.labels, e2.labels) || conflicts(e1.precond, e2.precond))
end

function revoke_finals(p::Function, dfa::DFA)
    newnodes = Dict{DFANode,DFANode}()
    new(s) = get!(newnodes, s) do
        return DFANode(s.final && !p(s), s.eof_actions, s.nfanodes)
    end
    for s in traverse(dfa.start)
        s′ = new(s)
        for (e, t) in s.edges
            push!(s′.edges, (e, new(t)))
        end
    end
    return DFA(new(dfa.start))
end

function dfa2nfa(dfa::DFA)
    newnodes = Dict{DFANode,NFANode}()
    new(s) = get!(newnodes, s, NFANode())
    final = NFANode()
    for s in traverse(dfa.start)
        s′ = new(s)
        for (e, t) in s.edges
            push!(s′.edges, (e, new(t)))
        end
        if s.final
            push!(s′.edges, (Edge(eps, s.eof_actions), final))
        end
    end
    start = NFANode()
    push!(start.edges, (Edge(eps), new(dfa.start)))
    return NFA(start, final)
end

function remove_dead_nodes(dfa::DFA)
    backrefs = Dict(dfa.start => Set{DFANode}())
    for s in traverse(dfa.start), (_, t) in s.edges
        push!(get!(backrefs, t, Set{DFANode}()), s)
    end

    alive = Set{DFANode}()
    unvisited = [s for s in keys(backrefs) if s.final]
    while !isempty(unvisited)
        s = pop!(unvisited)
        push!(alive, s)
        for t in backrefs[s]
            if t ∉ alive
                push!(unvisited, t)
            end
        end
    end
    @assert dfa.start ∈ alive

    newnodes = Dict{DFANode,DFANode}()
    new(s) = get!(newnodes, s, DFANode(s.final, s.eof_actions, s.nfanodes))
    isvisited(s) = haskey(newnodes, s)
    unvisited = [dfa.start]
    while !isempty(unvisited)
        s = pop!(unvisited)
        s′ = new(s)
        for (e, t) in s.edges
            if t ∈ alive
                if !isvisited(t)
                    push!(unvisited, t)
                end
                push!(s′.edges, (e, new(t)))
            end
        end
    end

    return DFA(new(dfa.start))
end


text:
macro debug(expression)
    return :(@assert $expression)
    # return nothing
end


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

plain(x) = sprint(plain, x)

function plain(io::IO, content::Vector)
    isempty(content) && return
    for md in content[1:end-1]
        plain(io, md)
        println(io)
    end
    plain(io, content[end])
end

plain(io::IO, md::MD) = plain(io, md.content)

function plain(io::IO, header::Header{l}) where l
    print(io, "#"^l*" ")
    plaininline(io, header.text)
    println(io)
end

function plain(io::IO, code::Code)
    # If the code includes a fenced block this will break parsing,
    # so it must be enclosed by a longer ````-sequence.
    n = mapreduce(m -> length(m.match), max, 2, eachmatch(r"^`+"m, code.code)) + 1
    println(io, "`" ^ n, code.language)
    println(io, code.code)
    println(io, "`" ^ n)
end

function plain(io::IO, p::Paragraph)
    plaininline(io, p.content)
    println(io)
end

function plain(io::IO, list::List)
    for (i, item) in enumerate(list.items)
        print(io, isordered(list) ? "$(i + list.ordered - 1). " : "  * ")
        lines = split(rstrip(sprint(plain, item)), "\n")
        for (n, line) in enumerate(lines)
            print(io, (n == 1 || isempty(line)) ? "" : "    ", line)
            n < length(lines) && println(io)
        end
        println(io)
    end
end

function plain(io::IO, q::BlockQuote)
    s = sprint(plain, q.content)
    for line in split(rstrip(s), "\n")
        println(io, isempty(line) ? ">" : "> ", line)
    end
    println(io)
end

function plain(io::IO, f::Footnote)
    print(io, "[^", f.id, "]:")
    s = sprint(plain, f.text)
    lines = split(rstrip(s), "\n")
    # Single line footnotes are printed on the same line as their label
    # rather than taking up an additional line.
    if length(lines) == 1
        println(io, " ", lines[1])
    else
        println(io)
        for line in lines
            println(io, isempty(line) ? "" : "    ", line)
        end
        println(io)
    end
end

function plain(io::IO, md::Admonition)
    s = sprint(plain, md.content)
    title = md.title == uppercasefirst(md.category) ? "" : " \"$(md.title)\""
    println(io, "!!! ", md.category, title)
    for line in split(rstrip(s), "\n")
        println(io, isempty(line) ? "" : "    ", line)
    end
    println(io)
end

function plain(io::IO, md::HorizontalRule)
    println(io, "-" ^ 3)
end

function plain(io::IO, l::LaTeX)
    println(io, '$', '$')
    println(io, l.formula)
    println(io, '$', '$')
end

function plain(io::IO, md)
    show(io,  MIME"text/plain"(), md)
    println(io)
end

# Inline elements

plaininline(x) = sprint(plaininline, x)

function plaininline(io::IO, md...)
    for el in md
        plaininline(io, el)
    end
end

plaininline(io::IO, md::Vector) = !isempty(md) && plaininline(io, md...)

plaininline(io::IO, f::Footnote) = print(io, "[^", f.id, "]")

plaininline(io::IO, link::Link) = plaininline(io, "[", link.text, "](", link.url, ")")

plaininline(io::IO, md::Image) = plaininline(io, "![", md.alt, "](", md.url, ")")

plaininline(io::IO, s::AbstractString) = print(io, s)

plaininline(io::IO, md::Bold) = plaininline(io, "**", md.text, "**")

plaininline(io::IO, md::Italic) = plaininline(io, "*", md.text, "*")

function plaininline(io::IO, md::Code)
    if occursin("`", md.code)
        n = maximum(length(m.match) for m in eachmatch(r"(`+)", md.code))
        s = "`"^((iseven(n) ? 1 : 2) + n)
        print(io, s, Base.startswith(md.code, "`") ? " " : "")
        print(io, md.code, endswith(md.code, "`") ? " " : "", s)
    else
        print(io, "`", md.code, "`")
    end
end

plaininline(io::IO, br::LineBreak) = println(io)

plaininline(io::IO, x) = show(io, MIME"text/plain"(), x)

# show

Base.show(io::IO, md::MD) = plain(io, md)
Base.show(io::IO, ::MIME"text/markdown", md::MD) = plain(io, md)

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
## Activate project environment
import Pkg; Pkg.activate(".")

Pkg.instantiate()

## Load required packages

# Manipulating variables, data frames and files
import CSV
using DataFrames
using DelimitedFiles
using JLD

# Solving systems of equations numerically
import Ipopt
using JuMP

# Doing statistics and models
using Distributions
using LinearAlgebra
using ProgressMeter
using SparseArrays
using Statistics
using StatsBase
using Turing

# Making plots
using Plots
using StatsPlots

# Analyzing ecological networks
using EcologicalNetworks
using Mangal
using RandomBooleanMatrices

## Load custom functions
include(joinpath("code", "functions", "predict_links.jl")) # predict the number of links from the number of species
include(joinpath("code", "functions", "joint_degree_dist_maxent.jl")) # gives the joint degree distribution of maximum entropy given numbers of species and links
include(joinpath("code", "functions", "degree_dist_maxent.jl")) # gives the degree distribution of maximum entropy given numbers of species and links
include(joinpath("code", "functions", "joint_degree_seq.jl")) # gives the joint degree sequence of a unipartite network
include(joinpath("code", "functions", "joint_degree_dist.jl")) # gives the joint degree distribution of a unipartite network
include(joinpath("code", "functions", "network_maxent.jl")) # gives the joint degree distribution of a unipartite network

## Load scripts
include(joinpath("code", "01_import_mangal_metadata.jl"))
include(joinpath("code", "02_predict_networks.jl"))
include(joinpath("code", "03_make_figures.jl"))


text:
function acronym(phrase)
    phrase = replace(uppercase(phrase), "'" => "")
    phrase = replace(phrase, r"[^A-Z]+" => " ")

    out = ""
    for word in split(phrase, " ")
        out *= word[1]
    end
    return out
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
function Bounds(DF,TXTL,state_array)
FB = DF["default_flux_bounds_array"]

Vmax = 600 # 1000 s^-1 * (3600 s/h) * (0.000167 mM)

FB[1:191,1] = 0
FB[1:191,2] = Vmax

FB[192:end,1] = 0
FB[192:end,2] = 0
FB[197,2] = 100 # hydrogen transport

#==============================================TXTL=====================================================#
#RNAP_elongation_rate = TXTL["RNAP_elongation_rate"];
RIBOSOME_concentration = state_array[141]#TXTL["RIBOSOME_concentration"];
RIBOSOME_elongation_rate = TXTL["RIBOSOME_elongation_rate"];
#mRNA_degradation_rate = TXTL["mRNA_degradation_rate"]*2;
mRNA_length = TXTL["mRNA_length"];
protein_length = TXTL["protein_length"];
gene_copies = TXTL["gene_copies"];
volume = TXTL["volume"];
polysome_amplification = TXTL["polysome_gain"];
plasmid_saturation_coefficient = TXTL["plasmid_saturation_coefficient"]*1e-6;
mRNA_saturation_coefficient = TXTL["mRNA_saturation_coefficient"]
Promoter = TXTL["Promoter"]
inducer = TXTL["inducer"]
#====================================Transcription===================================================#
#Compute the promoter strength P -
P = 0.9 #0.8*(K1)/(1+K1);
RNAP_elongation_rate = TXTL["RNAP_elongation_rate"]
plasmid_saturation_coefficient = TXTL["plasmid_saturation_coefficient"]
# 1 GENE_CAT
TX_Vmax = RNAP_elongation_rate/mRNA_length*3600*P
TX = TX_Vmax*state_array[143]*state_array[1]/(plasmid_saturation_coefficient+state_array[1]);

#====================================Translation===================================================#
polysome_amplification = TXTL["polysome_amplification"]
RIBOSOME_elongation_rate = TXTL["RIBOSOME_elongation_rate"]
mRNA_saturation_coefficient = TXTL["mRNA_saturation_coefficient"]
mRNA_degradation_rate = TXTL["mRNA_degradation_rate"]
# 141 RIBOSOME
# 144 mRNA_CAT
TL_Vmax = polysome_amplification*(3*RIBOSOME_elongation_rate)*(1/mRNA_length)*3600
TL = TL_Vmax*state_array[141]*state_array[144]/(mRNA_saturation_coefficient+state_array[144])
DEG = state_array[144]*mRNA_degradation_rate
#===================================================================================================#
FB[167,1] = TX #transcriptional initiation
FB[168,1] = TX #transcriptional elongation
FB[167,2] = TX #transcriptional initiation
FB[168,2] = TX #transcriptional elongation
FB[169,1] = 0 #mRNA_degradation
FB[169,2] = DEG #mRNA_degradation
FB[170,1] = 0 #translation initiation
FB[171,1] = 0 #translation elongation
FB[170,2] = TL #translation initiation
FB[171,2] = TL #translation elongation
FB[172:191,1] = 0 #tRNA charging
FB[172:191,2] = 50*TL #tRNA charging
FB = max(0,FB)
DF["default_flux_bounds_array"] = FB
return DF
end


text:
import QuantumOptics.printing: permuted_sparsedata, permuted_densedata
export tex

"""
    tex(x, statename="ψ") # for Ket and Bra
    tex(x, statename="Operator") # for Operator

LaTeX (ASCII art) rendering for composite (SpinBasis and NLevelBasis) system's Ket, Bra, and Operator on IJulia (REPL).
If you use this on IJulia, you can change LHS state name by second argument.

# Example (REPL)
```julia
ulia> using QuantumInformation.ShortNames

julia> ψ = cnot() * (H() ⊗ id()) * qubit("00")
Ket(dim=4)
  basis: [Spin(1/2) ⊗ Spin(1/2)]
 0.707107+0.0im
      0.0+0.0im
      0.0+0.0im
 0.707107+0.0im

julia> tex(ψ)
|ψ⟩ = 0.707|00⟩ + 0.707|11⟩

julia> dm(ψ)
DenseOperator(dim=4x4)
  basis: [Spin(1/2) ⊗ Spin(1/2)]
 0.5+0.0im  0.0+0.0im  0.0+0.0im  0.5+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.5+0.0im  0.0+0.0im  0.0+0.0im  0.5+0.0im

julia> dm(ψ) |> tex
Operator = 0.5 |00⟩⟨00| +0.5 |00⟩⟨11| +0.5 |11⟩⟨00| +0.5 |11⟩⟨11|
```
"""
function tex(x::Union{Ket,Bra}, statename::String="")
    if isempty(statename) && isdefined(Main, :IJulia) && Main.IJulia.inited
        statename = "\\psi"
    elseif isempty(statename)
        statename = "ψ"
    end

    if all(x.basis.shape .== 2)
        if isdefined(Main, :IJulia) && Main.IJulia.inited
            str = _mdbit(x, statename)
            display("text/latex", "\$" * str * "\$") # for IJulia
        else
            str = _aabit(x, statename)
            println(str) # for REPL
        end
    else
        if isdefined(Main, :IJulia) && Main.IJulia.inited
            str = _mdxit(x, statename)
            display("text/latex", "\$" * str * "\$") # for IJulia
        else
            str = _aaxit(x, statename)
            println(str) # for REPL
        end
    end
end
function tex(x::Union{DenseOperator,SparseOperator}, statename::String="")
    if isempty(statename) && isdefined(Main, :IJulia) && Main.IJulia.inited
        statename = "\\mathrm{Operator}"
    elseif isempty(statename)
        statename = "Operator"
    end
    if all(x.basis_l.shape .== 2) && all(x.basis_r.shape .== 2)
        if isdefined(Main, :IJulia) && Main.IJulia.inited
            str = _mdbit(x, statename)
            display("text/latex", "\$" * str * "\$") # for IJulia
        else
            str = _aabit(x, statename)
            println(str) # for REPL
        end
    else
        if isdefined(Main, :IJulia) && Main.IJulia.inited
            str = _mdxit(x, statename)
            display("text/latex", "\$" * str * "\$") # for IJulia
        else
            str = _aaxit(x, statename)
            println(str) # for REPL
        end
    end
end

"make markdown for Bra and Ket (qubit system)"
function _mdbit(x::Union{Ket,Bra}, statename::String)
    nq = length(x.basis.shape)
    isfirstterm = true
    if nq != 1
        perm = collect(reverse(1:nq))
        x = permutesystems(x, perm)
    end
    data = x.data
    braket = ifelse(typeof(x) == Ket, ["|", "\\rangle"], ["\\langle", "|"])

    str = "$(braket[1]) $(statename) $(braket[2]) = "
    for (idx, ent) in enumerate(data)
        if ent ≈ 1
            value = 1.0
        else
            value = ent
        end
        if !(ent ≈ 0)
            if isfirstterm
                isfirstterm = false
                if value == 1.
                    str *= "$(braket[1]) $(bin(idx-1, nq)) $(braket[2])"
                else
                    str *= "$(n2s(ent)) $(braket[1]) $(bin(idx-1, nq)) $(braket[2])"
                end
            else
                if n2s(ent)[1] == '-'
                    str *= "$(n2s(ent)) $(braket[1]) $(bin(idx-1, nq)) $(braket[2])"
                else
                    str *= "+$(n2s(ent)) $(braket[1]) $(bin(idx-1, nq)) $(braket[2])"
                end
            end
        end
    end

    return str
end
function _mdxit(x::Union{Ket,Bra}, statename::String) # quxit system
    shape = x.basis.shape
    nq = length(shape)
    isfirstterm = true
    perm = collect(reverse(1:nq))
    if nq != 1
        x = permutesystems(x, perm)
    end
    data = x.data
    braket = ifelse(typeof(x) == Ket, ["|", "\\rangle"], ["\\langle", "|"])

    str = "$(braket[1]) $(statename) $(braket[2]) = "
    for (idx, ent) in enumerate(data)
        if ent ≈ 1
            value = 1.0
        else
            value = ent
        end
        if !(ent ≈ 0)
            if isfirstterm
                isfirstterm = false
                if value == 1.
                    str *= "$(braket[1]) $(xit(idx, shape)) $(braket[2])"
                else
                    str *= "$(n2s(ent)) $(braket[1]) $(xit(idx, shape)) $(braket[2])"
                end
            else
                if n2s(ent)[1] == '-'
                    str *= "$(n2s(ent)) $(braket[1]) $(xit(idx, shape)) $(braket[2])"
                else
                    str *= "+$(n2s(ent)) $(braket[1]) $(xit(idx, shape)) $(braket[2])"
                end
            end
        end
    end

    return str
end

"make markdown for Operator (qubits system)"
function _mdbit(x::Union{DenseOperator,SparseOperator}, statename::String)
    nq = x.basis_l.shape |> length
    isfirstterm = true
    if typeof(x) == DenseOperator
        data = permuted_densedata(x)
    elseif typeof(x) == SparseOperator
        data = permuted_sparsedata(x)
    end

    str = "$(statename) = "
    for i in 1:2^nq
        for j in 1:2^nq
            ent = data[i,j]
            if ent ≈ 1
                value = 1.0
            else
                value = ent
            end
            if !(ent ≈ 0)
                if isfirstterm
                    isfirstterm = false
                    if value == 1.
                        str *= "| $(bin(i-1, nq)) \\rangle \\langle $(bin(j-1, nq)) |"
                    else
                        str *= "$(n2s(ent)) | $(bin(i-1, nq)) \\rangle \\langle $(bin(j-1, nq)) |"
                    end
                else
                    if value == 1.0
                        str *= "+ | $(bin(i-1, nq)) \\rangle \\langle $(bin(j-1, nq)) |"
                    elseif n2s(ent)[1] == '-'
                        str *= "$(n2s(ent)) | $(bin(i-1, nq)) \\rangle \\langle $(bin(j-1, nq)) |"
                    else
                        str *= "+$(n2s(ent)) | $(bin(i-1, nq)) \\rangle \\langle $(bin(j-1, nq)) |"
                    end
                end
            end

        end
    end

    return str
end
function _mdxit(x::Union{DenseOperator,SparseOperator}, statename::String)
    rshape = x.basis_r.shape
    lshape = x.basis_l.shape
    ncol = prod(rshape)
    nrow = prod(lshape)

    isfirstterm = true
    if typeof(x) == DenseOperator
        data = permuted_densedata(x)
    elseif typeof(x) == SparseOperator
        data = permuted_sparsedata(x)
    end

    str = "$(statename) = "
    for i in 1:nrow
        for j in 1:ncol
            ent = data[i,j]
            if ent ≈ 1
                value = 1.0
            else
                value = ent
            end
            if !(ent ≈ 0)
                if isfirstterm
                    isfirstterm = false
                    if value == 1.
                        str *= "| $(xit(i, lshape)) \\rangle \\langle $(xit(j, rshape)) |"
                    else
                        str *= "$(n2s(ent)) | $(xit(i, lshape)) \\rangle \\langle $(xit(j, rshape)) |"
                    end
                else
                    if value == 1.0
                        str *= "+ | $(xit(i, lshape)) \\rangle \\langle $(xit(j, rshape)) |"
                    elseif n2s(ent)[1] == '-'
                        str *= "$(n2s(ent)) | $(xit(i, lshape)) \\rangle \\langle $(xit(j, rshape)) |"
                    else
                        str *= "+$(n2s(ent)) | $(xit(i, lshape)) \\rangle \\langle $(xit(j, rshape)) |"
                    end
                end
            end

        end
    end

    return str
end

"make ASCIIart for Bra and Ket (qubits system)"
function _aabit(x::Union{Ket,Bra}, statename::String="")
    nq = length(x.basis.shape)
    isfirstterm = true
    perm = collect(reverse(1:nq))
    if nq != 1
        x = permutesystems(x, perm)
    end
    data = x.data
    # braket = Dict(Ket=>["|", ">"], Bra=>["<", "|"])[typeof(x)]
    braket = ifelse(typeof(x) == Ket, ["|", "⟩"], ["⟨", "|"])

    str = "$(braket[1])$(statename)$(braket[2]) = "
    for (idx, ent) in enumerate(data)
        if ent ≈ 1
            value = 1.0
        else
            value = ent
        end
        if !(ent ≈ 0)
            if isfirstterm
                isfirstterm = false
                if value == 1.
                    str *= "$(braket[1])$(bin(idx-1, nq))$(braket[2])"
                else
                    str *= "$(n2s(ent))$(braket[1])$(bin(idx-1, nq))$(braket[2])"
                end
            else
                if n2s(ent)[1] == '-'
                    str *= " $(n2s(ent))$(braket[1])$(bin(idx-1, nq))$(braket[2])"
                else
                    str *= " + $(n2s(ent))$(braket[1])$(bin(idx-1, nq))$(braket[2])"
                end
            end
        end
    end

    return str
end
function _aaxit(x::Union{Ket,Bra}, statename::String="")
    shape = x.basis.shape
    nq = length(shape)
    isfirstterm = true
    perm = collect(reverse(1:nq))
    if nq != 1
        x = permutesystems(x, perm)
    end
    data = x.data
    braket = ifelse(typeof(x) == Ket, ["|", "⟩"], ["⟨", "|"])

    str = "$(braket[1])$(statename)$(braket[2]) = "
    for (idx, ent) in enumerate(data)
        if ent ≈ 1
            value = 1.0
        else
            value = ent
        end
        if !(ent ≈ 0)
            if isfirstterm
                isfirstterm = false
                if value == 1.
                    str *= "$(braket[1])$(xit(idx, shape))$(braket[2])"
                else
                    str *= "$(n2s(ent))$(braket[1])$(xit(idx, shape))$(braket[2])"
                end
            else
                if n2s(ent)[1] == '-'
                    str *= " $(n2s(ent))$(braket[1])$(xit(idx, shape))$(braket[2])"
                else
                    str *= " + $(n2s(ent))$(braket[1])$(xit(idx, shape))$(braket[2])"
                end
            end
        end
    end

    return str
end
"make ASCIIart for Operator"
function _aabit(x::Union{DenseOperator,SparseOperator}, statename::String="")
    nq = x.basis_l |> length |> log2 |> Int
    isfirstterm = true
    if typeof(x) == DenseOperator
        data = permuted_densedata(x)
    elseif typeof(x) == SparseOperator
        data = permuted_sparsedata(x)
    end

    str = "$(statename) = "
    for i in 1:2^nq
        for j in 1:2^nq
            ent = data[i,j]
            if ent ≈ 1
                value = 1.0
            else
                value = ent
            end
            if !(ent ≈ 0)
                if isfirstterm
                    isfirstterm = false
                    if value == 1.
                        str *= "|$(bin(i-1, nq))⟩⟨$(bin(j-1, nq))|"
                    else
                        str *= "$(n2s(ent)) |$(bin(i-1, nq))⟩⟨$(bin(j-1, nq))|"
                    end
                else
                    if value == 1.0
                        str *= " + |$(bin(i-1, nq))⟩⟨$(bin(j-1, nq))|"
                    elseif n2s(ent)[1] == '-'
                        str *= " $(n2s(ent)) |$(bin(i-1, nq))⟩⟨$(bin(j-1, nq))|"
                    else
                        str *= " +$(n2s(ent)) |$(bin(i-1, nq))⟩⟨$(bin(j-1, nq))|"
                    end
                end
            end

        end
    end

    return str
end
function _aaxit(x::Union{DenseOperator,SparseOperator}, statename::String="")
    isfirstterm = true
    rshape = x.basis_r.shape
    lshape = x.basis_l.shape
    ncol = prod(rshape)
    nrow = prod(lshape)

    isfirstterm = true
    if typeof(x) == DenseOperator
        data = permuted_densedata(x)
    elseif typeof(x) == SparseOperator
        data = permuted_sparsedata(x)
    end

    str = "$(statename) = "
    for i in 1:nrow
        for j in 1:ncol
            ent = data[i,j]
            if ent ≈ 1
                value = 1.0
            else
                value = ent
            end
            if !(ent ≈ 0)
                if isfirstterm
                    isfirstterm = false
                    if value == 1.
                        str *= "|$(xit(i, lshape))⟩⟨$(xit(j, rshape))|"
                    else
                        str *= "$(n2s(ent)) |$(xit(i, lshape))⟩⟨$(xit(j, rshape))|"
                    end
                else
                    if value == 1.0
                        str *= "+ |$(xit(i, lshape))⟩⟨$(xit(j, rshape))|"
                    elseif n2s(ent)[1] == '-'
                        str *= "$(n2s(ent)) |$(xit(i, lshape))⟩⟨$(xit(j, rshape))|"
                    else
                        str *= "+$(n2s(ent)) |$(xit(i, lshape))⟩⟨$(xit(j, rshape))|"
                    end
                end
            end

        end
    end
    return str
end

"""
    n2s(x)

number to string for tex function.
"""
function n2s(x::T) where T <: Complex
    str = ""
    str1, str2 = "", ""
    if !(x.re ≈ 0)
        str1 *= "$(round(x.re, 3))"
    end
    if !(x.im ≈ 0)
        if x.im ≈ one(x.im)
            str2 *= "i"
        elseif -x.im ≈ one(x.im)
            str2 *= "-i"
        else
            str2 *= "$(round(x.im, 3))i"
        end
    end
    if !isempty(str1) && x.im > zero(x.im)
        str2 = "+" * str2
    end
    str = str1 * str2
    if !isempty(str1) && !isempty(str2)
        str = "(" * str * ")"
    end
    return str
end

"""
    xit(n::Int, dims::Vector{Int})

# Arguments
- `n::Int`: n th row(column) of qubit/operator array.
- `dims::Vector{Int}`: N-ary array
"""
function xit(m::Int, dims::Vector{Int})
    m = m - 1
    str = ""
    nquxit = length(dims)
    product = prod(dims[2:end])
    for ith in 1:nquxit-1
        d = div(m, product)
        m = m - d * product
        product = div(product, dims[ith+1])
        str *= string(d)
    end
    str *= string(m)
    return str
end


text:
# Poly type manipulations

module Polynomial
#todo: sparse polynomials?

export Poly, polyval, polyint, polyder, poly, roots, degree
export polydir #Deprecated

import Base: length, endof, getindex, setindex!, copy, zero, one, convert
import Base: show, print, *, /, //, -, +, ==, divrem, rem, eltype

eps{T}(::Type{T}) = convert(T,0)
eps{F<:FloatingPoint}(x::Type{F}) = Base.eps(F)
eps{T}(x::Type{Complex{T}}) = eps(T)

immutable Poly{T<:Number}
    a::Vector{T}
    nzfirst::Int #for effiencicy, track the first non-zero index
    var::Symbol
    function Poly(a::Vector{T}, var)
        nzfirst = 0 #find and chop leading zeros
        for i = 1:length(a)
            if abs(a[i]) > 2*eps(T)
                break
            end
            nzfirst = i
        end
        new(a, nzfirst, symbol(var))
    end
end

Poly{T<:Number}(a::Vector{T}, var::Union(String,Symbol,Char)=:x) = Poly{T}(a, var)

convert{T}(::Type{Poly{T}}, p::Poly) = Poly(convert(Vector{T}, p.a), p.var)
promote_rule{T, S}(::Type{Poly{T}}, ::Type{Poly{S}}) = Poly{promote_type(T, S)}
eltype{T}(::Poly{T}) = T

length(p::Poly) = length(p.a)-p.nzfirst
endof(p::Poly) = length(p)
deg(p::Poly) = length(p) - 1

getindex(p::Poly, i) = p.a[i+p.nzfirst]
setindex!(p::Poly, v, i) = (p.a[i+p.nzfirst] = v)

copy(p::Poly) = Poly(copy(p.a[1+p.nzfirst:end]), p.var)

zero{T}(p::Poly{T}) = Poly([zero(T)], p.var)
zero{T}(::Type{Poly{T}}) = Poly([zero(T)])
one{T}(p::Poly{T}) = Poly([one(T)], p.var)
one{T}(::Type{Poly{T}}) = Poly([one(T)])

function show(io::IO, p::Poly)
    print(io,"Poly(")
    print(io,p)
    print(io,")")
end

function print{T<:Real}(io::IO, p::Poly{T})
    n = length(p)
    if n <= 0
        print(io,"0")
    else
        for j = 1:n
            pj = p[j]
            magpj = abs(pj)
            if magpj > 2*eps(T)
                if j == 1
                    pj < 0 && print(io, "-")    #Prepend - if first and negative
                else
                    pj < 0 ? print(io," - ") : print(" + ")
                end
                #Print pj if pj is the last coefficient, or pj is not identically 1
                if j == n || abs(magpj - 1) > 2*eps(T)
                    print(io, magpj)
                end
                exp = n-j
                if exp > 0
                    print(io, p.var)
                    if exp > 1
                        print(io, '^', exp)
                    end
                end
            end
        end
    end
end

function print{T<:ComplexPair}(io::IO, p::Poly{T})
    n = length(p)
    if n <= 0
        print(io,"0")
    else
        for j = 1:n
            pj = p[j]
            abs_repj = abs(real(pj))
            abs_impj = abs(imag(pj))
            if abs(pj) > 2*eps(T)
                if !(abs_impj > 2*eps(T))
                    if j > 1
                        real(pj) < 0 ? print(io," - ") : print(" + ")
                    else
                        real(pj) < 0 && print(io, "-")    #Prepend - if first and negative
                    end
                else
                    j > 1 && print(io, " + ")
                end
                if abs_repj > 2*eps(T)    #Real part is not 0
                    if abs_impj > 2*eps(T)    #Imag part is not 0
                        print(io,'(',pj,')')
                    elseif abs(abs_repj - 1) > 2*eps(T) || j == n
                        print(io,abs_repj)
                    end
                else
                    if abs_impj > 2*eps(T)
                        print(io,'(', imag(pj),"im)")
                    end
                end
                exp = n-j
                if exp > 0
                    print(io, p.var)
                    if exp > 1
                        print(io, '^', exp)
                    end
                end
            end
        end
    end
end

*(c::Number, p::Poly) = Poly(c * p.a[1+p.nzfirst:end], p.var)
*(p::Poly, c::Number) = Poly(c * p.a[1+p.nzfirst:end], p.var)
/(p::Poly, c::Number) = Poly(p.a[1+p.nzfirst:end] / c, p.var)
-(p::Poly) = Poly(-p.a[1+p.nzfirst:end], p.var)

-(p::Poly, c::Number) = +(p, -c)
+(c::Number, p::Poly) = +(p, c)
function +(p::Poly, c::Number)
    if length(p) < 1
        return Poly([c,], p.var)
    else
        p2 = copy(p)
        p2.a[end] += c;
        return p2;
    end
end
function -(c::Number, p::Poly)
    if length(p) < 1
        return Poly([c,], p.var)
    else
        p2 = -p;
        p2.a[end] += c;
        return p2;
    end
end

function +{T,S}(p1::Poly{T}, p2::Poly{S})
    if p1.var != p2.var
        error("Polynomials must have same variable")
    end
    R = promote_type(T,S)
    n = length(p1)
    m = length(p2)
    if n > m
        a = Array(R, n)
        for i = 1:m
            a[n-m+i] = p1[n-m+i] + p2[i]
        end
        for i = 1:n-m
            a[i] = p1[i]
        end
    else
        a = Array(R, m)
        for i = 1:n
            a[m-n+i] = p1[i] + p2[m-n+i]
        end
        for i = 1:m-n
            a[i] = p2[i]
        end
    end
    Poly(a, p1.var)
end

function -{T,S}(p1::Poly{T}, p2::Poly{S})
    if p1.var != p2.var
        error("Polynomials must have same variable")
    end
    R = promote_type(T,S)
    n = length(p1)
    m = length(p2)
    if n > m
        a = Array(R, n)
        for i = 1:m
            a[n-m+i] = p1[n-m+i] - p2[i]
        end
        for i = 1:n-m
            a[i] = p1[i]
        end
    else
        a = Array(R, m)
        for i = 1:n
            a[m-n+i] = p1[i] - p2[m-n+i]
        end
        for i = 1:m-n
            a[i] = -p2[i]
        end
    end
    Poly(a, p1.var)
end

function *{T,S}(p1::Poly{T}, p2::Poly{S})
    if p1.var != p2.var
        error("Polynomials must have same variable")
    end
    R = promote_type(T,S)
    n = length(p1)
    m = length(p2)
    if n == 0 || m == 0
        return Poly(R[], p1.var)
    end
    a = zeros(R, n+m-1)
    for i = 1:length(p1)
        for j = 1:length(p2)
            a[i+j-1] += p1[i] * p2[j]
        end
    end
    Poly(a, p1.var)
end

function divrem{T, S}(num::Poly{T}, den::Poly{S})
    if num.var != den.var
        error("Polynomials must have same variable")
    end
    m = length(den)
    if m == 0
        throw(DivideError())
    end
    R = typeof(one(T)/one(S))
    n = length(num)
    deg = n-m+1
    if deg <= 0
        return zero(Poly{R}), convert(Poly{R}, num)
    end
    d = zeros(R, n)
    q = zeros(R, deg)
    r = zeros(R, n)
    r[:] = num.a[1+num.nzfirst:end]
    for i = 1:deg
        quot = r[i] / den[1]
        q[i] = quot
        if i > 1
            d[i-1] = 0
            r[i-1] = 0
        end
        for j = 1:m
            k = i+j-1
            elem = den[j]*quot
            d[k] = elem
            r[k] -= elem
        end
    end
    return Poly(q, num.var), Poly(r, num.var)
end
/(num::Poly, den::Poly) = divrem(num, den)[1]
rem(num::Poly, den::Poly) = divrem(num, den)[2]

function ==(p1::Poly, p2::Poly)
    if length(p1) != length(p2)
        return false
    elseif p1.var != p2.var
        return false
    else
        return p1.a[1+p1.nzfirst:end] == p2.a[1+p2.nzfirst:end]
    end
end

function polyval{T}(p::Poly{T}, x::Number)
    R = promote_type(T, typeof(x))
    lenp = length(p)
    if lenp == 0
        return zero(R)
    else
        y = convert(R, p[1])
        for i = 2:lenp
            y = p[i] + x.*y
        end
        return y
    end
end

polyval(p::Poly, v::AbstractVector) = map(x->polyval(p, x), v)

function polyint{T}(p::Poly{T}, k::Number=0)
    n = length(p)
    R = typeof(one(T)/1)
    a2 = Array(R, n+1)
    for i = 1:n
        a2[i] = p[i] / (n-i+1)
    end
    a2[end] = k
    Poly(a2, p.var)
end

@deprecate polydir polyder

function polyder{T}(p::Poly{T})
    n = length(p)
    if n > 0
        a2 = Array(T, n-1)
        for i = 1:n-1
            a2[i] = p[i] * (n-i)
        end
    else
        a2 = zeros(T, 0)
    end
    Poly(a2, p.var)
end

# create a Poly object from its roots
function poly{T}(r::AbstractVector{T}, var=:x)
    n = length(r)
    c = zeros(T, n+1)
    c[1] = 1
    for j = 1:n
        c[2:j+1] = c[2:j+1]-r[j]*c[1:j]
    end
    return Poly(c, var)
end
poly(A::Matrix, var=:x) = poly(eig(A)[1], var)
poly(A::Matrix, var::String) = poly(eig(A)[1], symbol(var))
poly(A::Matrix, var::Char) = poly(eig(A)[1], symbol(var))

roots{T}(p::Poly{Rational{T}}) = roots(convert(Poly{promote_type(T, Float64)}, p))

# compute the roots of a polynomial
function roots{T}(p::Poly{T})
    R = promote_type(T, Float64)
    num_zeros = 0
    if length(p) == 0
        return zeros(R, 0)
    end
    while abs(p[end-num_zeros]) <= 2*eps(T)
        if num_zeros == length(p)-1
            return zeros(R, 0)
        end
        num_zeros += 1
    end
    n = length(p)-num_zeros-1
    if n < 1
        return zeros(R, length(p)-1)
    end
    companion = zeros(R, n, n)
    a0 = p[end-num_zeros]
    for i = 1:n-1
        companion[1,i] = -p[end-num_zeros-i] / a0
        companion[i+1,i] = 1;
    end
    companion[1,end] = -p[1] / a0
    D,V = eig(companion)
    r = zeros(eltype(D),length(p)-1)
    r[1:n] = 1./D
    return r
end

end # module Poly


text:
module Exports

###############################################################################
using ..MeshSteward: initbox, updatebox!, boundingbox, inflatebox!, inbox, boxesoverlap, intersectboxes
export initbox, updatebox!, boundingbox, inflatebox!, inbox, boxesoverlap, intersectboxes

###############################################################################
using ..MeshSteward: vselect, connectedv
export vselect, connectedv

###############################################################################
using ..MeshSteward: eselect
export eselect

###############################################################################
using ..MeshSteward: linearspace, gradedspace
export linearspace, gradedspace

###############################################################################
using ..MeshSteward: T4blockx, T4block
export T4blockx, T4block

###############################################################################
using ..MeshSteward: T3blockx, T3block, T3toT6, T6blockx, T6block, T6toT3
export T3blockx, T3block, T3toT6, T6blockx, T6block, T6toT3

###############################################################################
using ..MeshSteward: Q4blockx, Q4block, Q4quadrilateral, Q4blockwdistortion
export Q4blockx, Q4block, Q4quadrilateral, Q4blockwdistortion

###############################################################################
using ..MeshSteward: L2blockx, L2block
export L2blockx, L2block

###############################################################################
using ..MeshSteward: import_MESH, import_NASTRAN, import_ABAQUS
export import_MESH, import_NASTRAN, import_ABAQUS

###############################################################################
using ..MeshSteward: vtkwrite, export_MESH
export vtkwrite, export_MESH

###############################################################################
using ..MeshSteward: transform, vconnected, vnewnumbering, compactify, fusevertices, withvertices, renumbered, cat, mergeirs, minimize_profile
export transform, vconnected, vnewnumbering, compactify, fusevertices, withvertices, renumbered, cat, mergeirs, minimize_profile

###############################################################################
using ..MeshSteward: Mesh
export Mesh
using ..MeshSteward: load, save, increl, attach!, basecode, nspacedims, baseincrel, geometry, summary, vselect, boundary, vertices, submesh, label
export load, save, increl, attach!, basecode, nspacedims, baseincrel, geometry, summary, vselect, boundary, vertices, submesh, label

end # module


text:
include("foo.jl")

using Foo: Data

struct CData
  name::Ptr{Cwchar_t}
  valuesPtr::Ptr{Cdouble}
  valuesLen::Csize_t
end

function Base.convert(::Type{Data}, s::CData)
  values = unsafe_wrap(Array, s.valuesPtr, s.valuesLen)
  name = wchartostring(s.name)
  Data(name, values)
end

function wchartostring(p::Ptr{Cwchar_t}; maxlen = 65536)
  len = 0
  while unsafe_load(p, len + 1) != 0 && len < maxlen; len += 1; end
  transcode(String, unsafe_wrap(Array{Cwchar_t,1}, p, len))
end

# These functions just convert C objects to Julia objects and call corresponding Julia functions.

Base.@ccallable function jl_FOO_bar(dataPtr::Ptr{CData}, dataLen::Csize_t)::Cint
  data = Vector{Data}(unsafe_wrap(Array, dataPtr, dataLen))
  Foo.bar(data)
end

Base.@ccallable function jl_FOO_baz(dataPtr::Ptr{CData}, dataLen::Csize_t, logPath::Ptr{Cwchar_t})::Cint
  data = Vector{Data}(unsafe_wrap(Array, dataPtr, dataLen))
  logPathJ = wchartostring(logPath)
  Foo.baz(data, logPathJ)
end


text:
module SWIM

using DataStructures
using Distances
using Rectangle
using Distributions
using Dates
using Base

export Canvas
export initialize
export generate_nodes
export Cell
export Node
export move
export update_seen
export simulate
export close_writer

include("Canvas.jl")
include("EventWriter.jl")
include("Node.jl")
include("Simulation.jl")

end


text:
module BigNumbers

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export bignumber

register_normal_element("st__big__number", context = @__MODULE__)

"""
bignumber(label::String = "",
          number::Union{Symbol,Number,Nothing} = nothing,
          args...;
          kwargs...)

Renders a Big Number UI element.

# Arguments
* label::Union{String,Symbol}
* number::Union{String,Symbol,Nothing,String}
* icon::Union{String,Symbol}
* color::Union{String,Symbol} = "positive"|"negative"
* arrow::Union{String,Symbol} = "up"|"down"
"""
function bignumber( label::Union{String,Symbol} = "",
                    number::Union{Symbol,Number,Nothing,String} = nothing,
                    args...;
                    wrap::Function = StippleUI.DEFAULT_WRAPPER,
                    kwargs...)
  wrap() do
    st__big__number(args...; attributes([:title => label, :number => number, kwargs...])...)
  end
end

end