text:
"""
    StatisticalSubstitutionModel

Subtype of `StatsBase.StatisticalModel`, to fit discrete data to a model
of trait substitution along a network.
See [`fitDiscrete`](@ref) to fit a trait substitution model to discrete data.
It returns an object of type `StatisticalSubstitutionModel`, to which standard
functions can be applied, like `loglikelihood(object)`, `aic(object)` etc.
"""
mutable struct StatisticalSubstitutionModel{T} <: StatsBase.StatisticalModel
    model::TraitSubstitutionModel{T} # T: type of trait labels
    net::HybridNetwork
    # data: trait[i] for leaf with n.number = i
    #       type Int: for indices of trait labels in model.label
    #       allows missing, but not half-ambiguous and not multi-state
    trait::Vector{Vector{Union{Missings.Missing,Int}}}
    ntraits::Int
    loglik::Union{Missings.Missing, Float64}
    # log of transition probabilities: where e goes from X to Y
    # logtrans[i,j,e] = P{Y=j|X=i} where i=start_state, j=end_state, e=edge.number
    # all displayed trees use same edge numbers and same logtrans as in full network
    logtrans::Array{Float64,3}   # size: k,k, net.numEdges where k=nStates(model)
    activetrait::Int
    # type based on extracting displayed trees
    displayedtree::Vector{HybridNetwork}
    priorltw::Vector{Float64} # prior log tree weight: log product of γ's
    postltw::Vector{Float64}  # posterior log tree weight: P{tree and data}
    # partial likelihoods for active trait, at indices [i, n.number or e.number ,t] :
    # - forward likelihood: log P{data below node n in tree t given state i at n}
    # - direct likelihood:  log P{data below edge e in tree t given i at parent of e}
    # - backward likelihood:log P{data at all non-descendants of n and state i at n}
    forwardlik::Array{Float64,3} # size: k, net.numNodes, number of displayed trees
    directlik::Array{Float64,3}  # size: k, net.numEdges, number of displayed trees
    backwardlik::Array{Float64,3}# size: k, net.numNodes, number of displayed trees
end
const SSM = StatisticalSubstitutionModel{T} where T

StatsBase.loglikelihood(obj::SSM) = obj.loglik
StatsBase.islinear(SSM) = false
StatsBase.dof(obj::SSM) = nparams(obj.model)
function Base.show(io::IO, obj::SSM)
    disp =  "$(typeof(obj)):\n"
    disp *= string(obj.model)
    disp *= "$(obj.ntraits) traits, $(length(obj.trait)) species, "
    disp *= "on a network with $(obj.net.numHybrids) reticulations"
    if !ismissing(obj.loglik)
        disp *= "\nlog-likelihood: $(obj.loglik)"
    end
    print(io, disp)
end
# nobs: ntraits * nspecies, minus any missing, but ignores correlation between species
# fixit: extend the StatsBase methods
# coef, coefnames, coeftable, confint,
# deviance (not from loglik in StatsBase), nulldeviance (from intercept only. here?)
# nullloglikelihood (from intercept only. here?)
# score (gradient of loglik)
# informationmatrix: Fisher by default, observed info matrix if expected=false
# stderror, vcov, params
# weights

"""
    fitDiscrete(net, model, tipdata)
    fitDiscrete(net, model, species, traits)

Calculate the maximum likelihood (ML) score of a network given
one or more discrete characters at the tips. Along each edge, transitions
are modelled with a continous time Markov `model`, whose parameters are
estimated (by maximizing the likelihood). At each hybrid node,
the trait is assumed to be inherited from either of the two immediate
parents according to the parents' average genetic contributions
(inheritance γ). The model ignores incomplete lineage sorting.
The algorithm extracts all trees displayed in the network.

Data can given in one of the following:
- `tipdata`: dictionary taxon => state label, for a single trait.
- `tipdata`: data frame for a single trait, in which case the taxon names
  are to appear in column 1 or in a column named "taxon" or "species", and
  trait *labels* are to appear in column 2 or in a column named "trait".
  Here, trait labels should be as they appear in `model.label`.
- `species`: vector of strings, and `traits`: DataFrame of traits,
  with rows in the order corresponding to the order of species names.
  Again, trait labels should be as they appear in `model.label`.
  All traits are assumed to follow the same model, with same parameters.

Optional arguments (default):
- `fixedparam` (false): should model rate parameters be fixed, or should they be optimized?
- `NLoptMethod` (`:LN_COBYLA`, derivative-free) for the optimization algorithm.
  For other options, see the
  [NLopt](https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/).
- tolerance values to control when the optimization is stopped:
  `ftolRel` (1e-12), `ftolAbs` (1e-10) on the likelihood, and
  `xtolRel` (1e-10), `xtolAbs` (1e-10) on the model parameters.
- `verbose` (false): if true, more information is output.

# examples:

```julia-repl
julia> net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);");

julia> m1 = BinaryTraitSubstitutionModel([0.1, 0.1], ["lo", "hi"]);

julia> dat = DataFrame(species=["C","A","B","D"], trait=["hi","lo","lo","hi"]);

julia> fit1 = fitDiscrete(net, m1, dat; fixedparam=true)
PhyloNetworks.StatisticalSubstitutionModel{String}:
Binary Trait Substitution Model:
rate lo→hi α=0.1
rate hi→lo β=0.1
1 traits, 4 species, on a network with 1 reticulations
log-likelihood: -3.107539646785388

julia> PhyloNetworks.fit!(fit1; fixedparam=false)
PhyloNetworks.StatisticalSubstitutionModel{String}:
Binary Trait Substitution Model:
rate lo→hi α=0.2722215661432007
rate hi→lo β=0.3498103666174014
1 traits, 4 species, on a network with 1 reticulations
log-likelihood: -2.727701700712135

julia> tips = Dict("A" => "lo", "B" => "lo", "C" => "hi", "D" => "hi");

julia> fit2 = fitDiscrete(net, m1, tips; xtolRel=1e-16, xtolAbs=1e-16, ftolRel=1e-16)
PhyloNetworks.StatisticalSubstitutionModel{String}:
Binary Trait Substitution Model:
rate lo→hi α=0.2722215661432007
rate hi→lo β=0.3498103666174014
1 traits, 4 species, on a network with 1 reticulations
log-likelihood: -2.727701700712135
```

Note that a copy of the network is stored in the fitted object,
but the internal representation of the network may be different in
`fit1.net` and in the original network `net`:

```julia-repl
julia> [n.number for n in fit2.net.node]
9-element Array{Int64,1}:
 1
 2
 9
 8
 3
 7
 6
 4
 5

julia> [n.number for n in net.node]
9-element Array{Int64,1}:
  1
  2
  3
 -4
  4
 -6
 -3
  5
 -2
```
"""
function fitDiscrete(net, model, tips::Dict; kwargs...)
    species = Array{String}(0)
    dat = Vector{Vector{Int}}(0) # indices of trait labels
    for (k,v) in tips
        !ismissing(v) || continue
        push!(species, k)
        vi = findfirst(model.label, v) # value index in model labels
        vi > 0 || error("trait $v not found in model")
        push!(dat, [vi])
    end
    o, net = check_matchtaxonnames!(species, dat, net)
    # dat[o] would make a shallow copy only
    StatsBase.fit(StatisticalSubstitutionModel, net, model, view(dat, o); kwargs...)
end

function fitDiscrete(net, model, dat::DataFrame; kwargs...)
    i = findfirst(DataFrames.names(dat), :taxon)
    if i==0 i = findfirst(DataFrames.names(dat), :species); end
    if i==0 i=1; end # first column if not column named "taxon" or "species"
    j = findfirst(DataFrames.names(dat), :trait)
    if j==0 j=2; end
    if i==j
        error("""expecting taxon names in column 'taxon', or 'species' or column 1,
              and trait values in column 'trait' or column 2.""")
    end
    species = copy(dat[i])    # modified in place later
    dat = traitlabels2indices(dat[j], model)   # vec of vec, indices
    o, net = check_matchtaxonnames!(species, dat, net)
    StatsBase.fit(StatisticalSubstitutionModel, net, model, view(dat, o); kwargs...)
end

function fitDiscrete(net, model, species::Array{String}, dat::DataFrame; kwargs...)
    dat2 = traitlabels2indices(dat, model) # vec of vec, indices
    o, net = check_matchtaxonnames!(copy(species), dat2, net)
    StatsBase.fit(StatisticalSubstitutionModel, net, model, view(dat2, o); kwargs...)
end

"""
    fit(StatisticalSubstitutionModel, net, model, traits; kwargs...)
    fit!(StatisticalSubstitutionModel; kwargs...)

Internal function called by [`fitDiscrete`](@ref): with same key word arguments `kwargs`.
But dangerous: `traits` should be a vector of vectors as for [`fitDiscrete`](@ref)
**but** here `traits` need to contain the *indices* of trait values corresponding
to the indices in `model.label`, and species should appear in `traits` in the
order corresponding to the node numbers in `net`.
See [`traitlabels2indices`](@ref) to convert trait labels to trait indices.

**Warning**: does *not* perform checks. [`fitDiscrete`](@ref) calls this function
after doing checks, preordering nodes in the network, making sure nodes have
consecutive numbers, species are matched between data and network etc.
"""
function StatsBase.fit(::Type{SSM}, net::HybridNetwork, model::TraitSubstitutionModel,
    trait::AbstractVector; kwargs...)
    T = eltype(model.label)
    # extract displayed trees
    trees = displayedTrees(net, 0.0; keepNodes=true)
    nnodes = length(net.node)
    for tree in trees
        preorder!(tree) # no need to call directEdges! before: already done on net
        length(tree.nodes_changed) == nnodes ||
            error("displayed tree with too few nodes: $(writeTopology(tree))")
        length(tree.edge) == length(net.edge)-net.numHybrids ||
            error("displayed tree with too few edges: $(writeTopology(tree))")
    end
    ntrees = length(trees)
    # log tree weights: sum log(γ) over edges, for each displayed tree
    priorltw = inheritanceWeight.(trees)
    k = nStates(model)
    # fixit: use SharedArray's below to parallelize things
    logtrans   = zeros(Float64, k,k,length(net.edge))
    forwardlik = zeros(Float64, k, nnodes,           ntrees)
    directlik  = zeros(Float64, k, length(net.edge), ntrees)
    backwardlik= zeros(Float64, k, nnodes,           ntrees)
    postltw    = Vector{Float64}(ntrees)
    # create new model object then fit:
    fit!(StatisticalSubstitutionModel{T}(deepcopy(model),
            net, trait, length(trait[1]), missing,
            logtrans, 1, trees, priorltw, postltw, forwardlik, directlik, backwardlik);
        kwargs...)
end

function fit!(obj::SSM; fixedparam=false::Bool, verbose=false::Bool,
      NLoptMethod=:LD_MMA::Symbol, ftolRel=fRelBL::Float64, ftolAbs=fAbsBL::Float64,
      xtolRel=xRelBL::Float64, xtolAbs=xAbsBL::Float64)

    if fixedparam
        discrete_corelikelihood!(obj)
        verbose && println("loglik = $(loglikelihood(obj)) under fixed parameters, no optimization")
    else
        # set-up optimization object
        NLoptMethod=:LN_COBYLA # no gradient
        # :LN_COBYLA for (non)linear constraits, :LN_BOBYQA for bound constraints
        npar = nparams(obj.model)
        opt = NLopt.Opt(NLoptMethod, npar)
        NLopt.ftol_rel!(opt,ftolRel) # relative criterion
        NLopt.ftol_abs!(opt,ftolAbs) # absolute criterion
        NLopt.xtol_rel!(opt,xtolRel)
        NLopt.xtol_abs!(opt,xtolAbs)
        NLopt.maxeval!(opt,1000) # max number of iterations
        # NLopt.maxtime!(opt, t::Real)
        NLopt.lower_bounds!(opt, zeros(Float64, npar))
        # fixit: set upper bound depending on branch lengths in network?
        counter = [0]
        function loglikfun(x::Vector{Float64}, grad::Vector{Float64}) # modifies obj
            counter[1] += 1
            obj.model.rate[:] = x
            res = discrete_corelikelihood!(obj)
            verbose && println("loglik: $res, model rates: $x")
            length(grad) == 0 || error("gradient not implemented")
            return res
        end
        NLopt.max_objective!(opt, loglikfun)
        fmax, xmax, ret = NLopt.optimize(opt, obj.model.rate) # optimization here!
        verbose && println("got $(round(fmax,5)) at $(round.(xmax,5)) after $(counter[1]) iterations (return code $(ret))")
    end
    # return fmax,xmax,ret
    return obj
end

"""
    discrete_corelikelihood!(obj::StatisticalSubstitutionModel; whichtrait=:all)
    discrete_corelikelihood_tree!(obj, t::Integer, traitrange::AbstractArray)

Calculate the likelihood and update `obj.loglik` for discrete characters on a network
(or on a single tree: `t`th tree displayed in the network, for the second form).
Update forward and direct partial likelihoods while doing so.
The algorithm extracts all displayed trees and weights the likelihood under all these trees.
"""

function discrete_corelikelihood!(obj::SSM; whichtrait=:all::Union{Symbol,Integer})
    if whichtrait == :all
        traitrange = 1:obj.ntraits
    elseif isinteger(whichtrait) && whichtrait > 0 && whichtrait <= obj.ntraits
        obj.activetrait = whichtrait
        traitrange = range(whichtrait, 1)
    elseif whichtrait == :active
        traitrange = range(obj.activetrait, 1)
    else
        error("'whichtrait' should be :all or :active or an integer in the correct range")
    end
    for edge in obj.net.edge # update logtrans: same for all displayed trees, all traits
        obj.logtrans[:,:,edge.number] = log.(P(obj.model, edge.length)) # element-wise
    end
    for t in 1:length(obj.displayedtree) # calculate P{data | tree t} & store in obj.postltw[t]
        discrete_corelikelihood_tree!(obj, t, traitrange)
    end
    # fixit: paralellize with
    # ll = pmap(t -> discrete_corelikelihood_tree!(obj,t), 1:ntrees)
    obj.postltw .+= obj.priorltw # P{tree t and data} .+= not += to re-use memory
    res = StatsFuns.logsumexp(obj.postltw)
    obj.loglik = res
    obj.postltw .-= res # now P{tree t | data}
    # fixit: function to get these posterior probabilities (just take exp.)
    return res
end

@doc (@doc discrete_corelikelihood!) discrete_corelikelihood_tree!
function discrete_corelikelihood_tree!(obj::SSM, t::Integer, traitrange::AbstractArray)
    tree = obj.displayedtree[t]
    # info("tree: $(writeTopology(tree))")
    forwardlik = view(obj.forwardlik, :,:,t)
    directlik  = view(obj.directlik,  :,:,t)
    k = nStates(obj.model)   # also = size(logtrans,1)
    fullloglik = 0.0
    for ci in traitrange     # ci = character index
      fill!(forwardlik, 0.0) # re-initialize for each trait, each iteration
      fill!(directlik,  0.0)
      for ni in reverse(1:length(tree.nodes_changed)) # post-order
        n = tree.nodes_changed[ni]
        nnum = n.number # same n.number across trees for a given node
        if n.leaf # need forwardlik initialized at 0: keep at 0 = log(1) if no data
            state = obj.trait[nnum][ci] # here: data assumed in a row n.number
            if !ismissing(state)
                for i in 1:k
                    forwardlik[i,nnum] = -Inf64 # log(0) = -Inf if i != observed state
                end
                forwardlik[state, nnum] = 0.
            end
        else # forward likelihood = product of direct likelihood over all children edges
            for e in n.edge
                n == getParent(e) || continue # to next edge if n is not parent of e
                forwardlik[:,nnum] .+= view(directlik, :,e.number)
            end
        end
        if ni==1 # root is first index in nodes changed
            logprior = [-log(k) for i in 1:k] # uniform prior at root
            loglik = logsumexp(logprior + view(forwardlik, :,nnum)) # log P{data for ci | tree t}
            fullloglik += loglik # add loglik of character ci
            break # out of loop over nodes
        end
        # if we keep going: n is not the root
        # calculate direct likelihood on the parent edge of n
        for e in n.edge
            if n == getChild(e)
                lt = view(obj.logtrans, :,:,e.number)
                for i in 1:k # state at parent node
                    directlik[i,e.number] = logsumexp(view(lt,i,:) + view(forwardlik,:,nnum))
                end
                break # we visited the parent edge: break out of for loop
            end
        end
      end # of loop over nodes
    end # of loop over traits
    obj.postltw[t] = fullloglik
    return fullloglik
end

"""
    traitlabels2indices(data, model::TraitSubstitutionModel)

Check that the character states in `data` are compatible with (i.e. subset of)
the trait labels in `model`. All columns are used.
`data` can be a DataFrame or a Matrix (multiple traits), or a Vector (one trait).

Return a vector of vectors (one per species) with integer entries,
where each state (label) is replaced by its index in `model`.
"""
function traitlabels2indices(data::AbstractVector, model::TraitSubstitutionModel)
    A = Vector{Vector{Union{Missings.Missing,Int}}}(0) # indices of trait labels
    labs = model.label
    for l in data
        vi = missing
        if !ismissing(l)
            vi = findfirst(model.label, l) # value index in model labels
            vi > 0 || error("trait $l not found in model")
        end
        push!(A, [vi])
    end
    return A
end
function traitlabels2indices(data::Union{AbstractMatrix,DataFrame},
                             model::TraitSubstitutionModel)
    A = Vector{Vector{Union{Missings.Missing,Int}}}(0) # indices of trait labels
    labs = model.label
    ntraits = size(data,2)
    for i in 1:size(data,1) # go row by row
        V = Vector{Union{Missings.Missing,Int}}(ntraits)
        for j in 1:ntraits
            vi = missing # value index
            @inbounds l = data[i,j] # value label
            if !ismissing(l)
                vi = findfirst(labs, l)
                vi > 0 || error("trait $l not found in model")
            end
            V[j] = vi
        end
        push!(A, V)
      end
    return A
end

"""
    check_matchtaxonnames!(species, data, net)

Modify `species` and `dat` by removing the species (rows) absent from the network.
Return a new network (`net` is *not* modified) with tips matching those in species:
if some species in `net` have no data, these species are pruned from the network.
The network also has its node names reset, such that leaves have nodes have
consecutive numbers starting at 1, with leaves first.
Used by [`fitDiscrete`](@ref) to build a new [`StatisticalSubstitutionModel`](@ref).
"""
function check_matchtaxonnames!(species::AbstractVector, dat::AbstractVector, net::HybridNetwork)
    # 1. basic checks for dimensions and types
    eltt = eltype(dat)
    @assert eltt <: AbstractVector "traits should be a vector of vectors"
    @assert Missings.T(eltype(eltt)) <: Integer "traits should be integers (label indices)"
    @assert !isempty(dat) "empty data vector!"
    ntraits = length(dat[1])
    for d in dat
        @assert length(d)==ntraits "all species should have the same number of traits"
    end
    @assert length(dat) == length(species) "need as many species as rows in trait data"
    # 2. match taxon labels between data and network
    netlab = tipLabels(net)
    ind2notinnet = find(x -> x ∉ netlab, species) # species not in network
    deleteat!(species, ind2notinnet)
    deleteat!(dat,     ind2notinnet)
    nvalues = [sum(.!ismissing.(d)) for d in dat] # species with completely missing data
    indmissing = find(x -> x==0, nvalues)
    deleteat!(species, indmissing)
    deleteat!(dat,     indmissing)
    indnotindat = find(x -> x ∉ species, netlab) # species not in data
    net = deepcopy(net)
    if !isempty(indnotindat)
        warn("the network contains taxa with no data: those will be pruned")
        for i in indnotindat
            deleteleaf!(net, netlab[i])
        end
    end
    # 3. calculate order of rows to have species with node.number i on ith row
    resetNodeNumbers!(net; checkPreorder=true) # tip species: now with numbers 1:n
    resetEdgeNumbers!(net) # to use edge as indices: 1:numEdges
    netlab = [n.name for n in sort(net.leaf, by = x -> x.number)]
    nspecies = length(netlab)
    o = Vector{Int}(nspecies)
    for i in 1:nspecies
        @inbounds o[i] = findfirst(species, netlab[i])
    end
    countnz(o) == nspecies || # number of non-zeros should be total size
        error("weird: even after pruning, species in network have no data")
    return (o,net)
end

"""
    ancestralStateReconstruction(obj::SSM, trait::Integer)
    ancestralStateReconstruction(obj::SSM)

Estimate the marginal probability of ancestral states for discrete character
number `trait`, or for the active trait if `trait` is unspecified: `obj.activetrait`.
The parameters of the [`StatisticalSubstitutionModel`](@ref) object `obj`
must first be fitted using [`fitDiscrete`](@ref), and ancestral state reconstruction
is conditional on the estimated parameters. If these parameters were estimated
using all traits, they are used as is to do ancestral state reconstruction of the
particular `trait` of interest.

**output**: data frame with a first column for the node numbers, a second column for
the node labels, and a column for each possible state: the entries in these columns
give the marginal probability that a given node has a given state.

warnings

- node numbers and node labels refer to those in `obj.net`, which might
  have a different internal representation of nodes than the original network
  used to build `obj`.
- `obj` is modified: its likelihood fields (forward, directional & backward)
  are updated to make sure that they correspond to the current parameter values
  in `obj.model`, and to the `trait` of interest.

See also [`discrete_backwardlikelihood_tree!`](@ref) to update `obj.backwardlik`.

# examples

```julia-repl
julia> net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);");

julia> m1 = BinaryTraitSubstitutionModel([0.1, 0.1], ["lo", "hi"]);

julia> dat = DataFrame(species=["C","A","B","D"], trait=["hi","lo","lo","hi"]);

julia> fit1 = fitDiscrete(net, m1, dat);

julia> asr = ancestralStateReconstruction(fit1)
9×4 DataFrames.DataFrame
│ Row │ nodenumber │ nodelabel │ lo       │ hi       │
├─────┼────────────┼───────────┼──────────┼──────────┤
│ 1   │ 1          │ A         │ 1.0      │ 0.0      │
│ 2   │ 2          │ B         │ 1.0      │ 0.0      │
│ 3   │ 3          │ C         │ 0.0      │ 1.0      │
│ 4   │ 4          │ D         │ 0.0      │ 1.0      │
│ 5   │ 5          │ 5         │ 0.286019 │ 0.713981 │
│ 6   │ 6          │ 6         │ 0.319454 │ 0.680546 │
│ 7   │ 7          │ 7         │ 0.168549 │ 0.831451 │
│ 8   │ 8          │ 8         │ 0.76736  │ 0.23264  │
│ 9   │ 9          │ #H1       │ 0.782777 │ 0.217223 │

julia> exp.(fit1.postltw) # marginal (posterior) probability that the trait evolved on each displayed tree
2-element Array{Float64,1}:
 0.919831
 0.0801689

julia> using PhyloPlots

julia> plot(fit1.net, :R, nodeLabel = asr[[:nodenumber, :lo]], tipOffset=0.2); # pp for "lo" state
```
"""
ancestralStateReconstruction(obj::SSM) = ancestralStateReconstruction(obj, obj.activetrait)
function ancestralStateReconstruction(obj::SSM, trait::Integer)
    # posterior probability of state i at node n: proportional to
    # sum_{trees t} exp( ltw[t] + backwardlik[i,n,t] + forwardlik[i,n,t] )
    bkd = obj.backwardlik
    fill!(bkd, 0.0) # initialize
    discrete_corelikelihood!(obj; whichtrait=trait) # update forward, direct, logtrans, postltw, loglik
    for t in 1:length(obj.displayedtree)
        discrete_backwardlikelihood_tree!(obj, t, trait)
    end
    # fixit: paralellize with
    # ll = pmap(t -> discrete_backwardlikelihood_tree!(obj,t, trait), 1:ntrees)
    k = nStates(obj.model)
    nnodes = length(obj.net.node)
    res = Array{Float64}((k,nnodes))
    frd = obj.forwardlik
    ltw = obj.priorltw
    for i in 1:k
        for n in 1:nnodes
            res[i,n] = logsumexp(ltw + view(bkd, i,n,:) + view(frd, i,n,:))
        end
    end
    ll = obj.loglik
    map!(x -> exp(x - ll), res, res)  # to normalize: p_i / sum(p_j over all states j)
    # alternative syntax: res .= exp.(res .- obj.loglik)
    nodestringlabels = Vector{String}(nnodes)
    for n in obj.net.node
        nodestringlabels[n.number] = (n.name == "" ? string(n.number) : n.name)
    end
    dat = DataFrame(transpose(res), Symbol.(obj.model.label))
    insert!(dat, 1, collect(1:nnodes), :nodenumber, makeunique=true)
    insert!(dat, 2, nodestringlabels,  :nodelabel,  makeunique=true)
    return dat
end

"""
    discrete_backwardlikelihood_tree!(obj::SSM, tree::Integer, trait::Integer)

Update `obj.backwardlik`; assume correct forward likelihood, directional likelihood
and transition probabilities.
"""
function discrete_backwardlikelihood_tree!(obj::SSM, t::Integer, trait::Integer)
    tree = obj.displayedtree[t]
    frdlik = view(obj.forwardlik, :,:,t)
    dirlik = view(obj.directlik , :,:,t)
    bkdlik = view(obj.backwardlik,:,:,t)
    k = nStates(obj.model)
    bkwtmp = Vector{Float64}(k) # to hold bkw lik without parent edge transition
    logprior = [-log(k) for i in 1:k]
    for ni in 1:length(tree.nodes_changed) # pre-order traversal to calculate backwardlik
        n = tree.nodes_changed[ni]
        nnum = n.number
        if ni == 1 # n is the root
            bkdlik[:,nnum] = logprior
        else
            pe = getMajorParentEdge(n)
            pn = getParent(pe)
            bkwtmp[:] = bkdlik[:,pn.number] # use bktmp's original memory
            for se in pn.edge
                if se != pe && pn == getParent(se) # then se is sister edge to pe
                    bkwtmp .+= view(dirlik, :,se.number)
                end
            end
            lt = view(obj.logtrans, :,:,pe.number)
            for j in 1:k # state at node n
                bkdlik[j,nnum] = logsumexp(bkwtmp + view(lt,:,j))
            end
        end
    end
    return nothing
end

# fixit: new type for two (dependent) binary traits
# need new algorithm for model at hybrid nodes: displayed trees aren't enough


text:
# wrap in a module for easy testing at the REPL
module NG

abstract type AbstractInventory end
abstract type AbstractPriceProcess end
abstract type AbstractMarket end

struct Inventory{T} <: AbstractInventory
    terminal_inventory::T
    min_inventory::T
    max_inventory::T
end

struct PriceProcess{T} <: AbstractPriceProcess
    initial_price::T
    volatility::T
    mean_level::T
    mean_rev_rate::T
end

struct Market{T, U} <: AbstractMarket
    horizon::U
    periods::U
    runs::U
    risk_free_rate::T
end

#function main()
    Δt = 3 * 1 / 1000
    market = Market(3,1000,10000,10/100)
    inventory = vcat(0:25:100,150,200:100:900,950:50:1100,1200:100:1800,1850,1900:25:2000)
    max_inventory = maximum(inventory)
    min_inventory = minimum(inventory)

    withdrawn_gas = withdraw.(inventory, Δt, min_inventory)
    purchased_gas = inject(inventory, time_step)

#end

function simulate(p::PriceProcess, m::Market)
    Δt = time_step(m)
    # initialize return value. Julia is column major and
    # access to period values per run will be need all at once,
    # hence they are the columns.
    retval = randn(Float64, m.runs, m.periods)
    for period in 1:m.periods
        prev_price = period == 1 ? fill(p.initial_price, m.runs) : retval[:,period-1]
        retval[:,period] .= prev_price .+ p.mean_rev_rate .* (p.mean_level .- prev_price) .*
            Δt .+ p.volatility .* prev_price .* sqrt(Δt) .* retval[:,period]
    end
    return retval
end

function time_step(m::Market)
    return m.horizon * 1 / m.periods
end

function terminal_value(prices,inventory, term_inv)
    return -2 * prices[:,end] .* max.(term_inv .- inventory, 0)
end

"""
    inject
    The rate of gas injected into the reservoir as a function of the current inventory.
"""
function inject(inventory, time_step)
    # hydro gas constants
    k2 = 730000
    k3 = 500
    k4 = 2500

    inven_factor = 1 / (inventory + k3)
    deduct_factor = 1 / k4
    return k2 * sqrt(inven_factor - deduct_factor) * time_step
end

"""
    withdraw releases gas from the reservoir based on the given inventory levels and
    the gas constants.
"""
function withdraw(inventory, time_step, min_inventory)
    k1 = 2040.41
    return max(k1 * sqrt(inventory) * time_step, min_inventory - inventory)
end

end

end #module

ext: jl

max_stars_repo_name: danzoppo/DataScienceExamples

max_issues_repo_name: danzoppo/DataScienceExamples

max_forks_repo_name: danzoppo/DataScienceExamples


text:
"""
    Meshes

- domain
- topology
- coordinates
- metric terms (inverse partial derivatives)
- quadrature rules and weights

## References / notes
 - [ceed](https://ceed.exascaleproject.org/ceed-code/)
 - [QA](https://github.com/CliMA/ClimateMachine.jl/blob/ans/sphere/test/Numerics/DGMethods/compressible_navier_stokes_equations/sphere/sphere_helper_functions.jl)

"""
module Spaces

import ..slab, ..column, ..level
import ..Utilities: PlusHalf
import ..DataLayouts, ..Geometry, ..Domains, ..Meshes, ..Topologies
using StaticArrays, ForwardDiff, LinearAlgebra, UnPack

abstract type AbstractSpace end

undertype(space::AbstractSpace) =
    Geometry.undertype(eltype(local_geometry_data(space)))

coordinates_data(space::AbstractSpace) = local_geometry_data(space).coordinates

include("quadrature.jl")
import .Quadratures

include("spectralelement.jl")
include("finitedifference.jl")
include("hybrid.jl")
include("triangulation.jl")
include("dss.jl")

end # module


text:
using GPlot, Colors, Random

continuous_preview(false)

begin # GENERATION OF FIGURES
    start = time()
    print("Generating example figures...")
    function gen(s::String; format="svg")
        include(joinpath(@__DIR__, "src", "exgen", s*".jl"))
        savefig(gcf(), s; format=format, path=joinpath(@__DIR__, "src", "exgen", "out"), res=200)
    end

    f = Figure(size=(10,8))

    # -------------------------
    # quickstart
    gen("qs_ex1")
    gen("qs_ex2")

    # -------------------------
    # line-scatter
    gen("ls_demo")
    gen("ls_ex1")
    gen("ls_ex2")
    gen("ls_ex3")
    gen("ls_ex4")
    gen("ls_ex5")
    gen("ls_ex5b")
    gen("ls_ex6")
    gen("ls_ex7")
    cla(); gen("ls_ex8")
    gen("ls_ex9")
    gen("ls_ex10")
    gen("ls_ex11")

    # -------------------------
    # hist-bar
    gen("hb_demo")
    gen("hb_ex1")
    gen("hb_ex2")
    gen("hb_ex3")
    gen("hb_ex4")
    gen("hb_ex5")
    gen("hb_ex6")
    gen("hb_ex7")
    gen("hb_ex8")

    # -------------------------
    gen("l_demo")
    gen("l_ex1")
    gen("l_ex1b")
    gen("l_ex1c")
    gen("l_ex2")
    gen("l_ex3")

    # -------------------------
    # fill-between
    gen("fb_demo")
    gen("fb_ex1")
    gen("fb_ex2")
    gen("fb_ex3")

    # -------------------------
    # global-parameters
    cla(); gen("glob_1")

    println(" done ✅  [$(round(time()-start,digits=1))s]")
end

begin
    start = time()
    print("Generating doc pages...")
    cdir = @__DIR__

    for fname = ["quickstart.md",
                 "line-scatter.md",
                 "hist-bar.md",
                 "legend.md",
                 "styling.md",
                 "figure-axes-axis.md",
                 "fill-between.md",
                 "boxplot.md",
                 "heatmap.md",
                 "annotations.md",
                 "global-settings.md"]
        open(joinpath(cdir, "src/man/$fname"), "w") do outf
            inf = read(joinpath(cdir, "src/man/_$fname"), String)
            matches = eachmatch(r"@@[A-Z]+:(.*\b)", inf)

            if isempty(matches)
                write(outf, inf)
            else
                head = 1
                for m ∈ matches
                    # PRE = write everything up to HEAD *not included*
                    write(outf, inf[head:prevind(inf, m.offset)])
                    # MOVE THE HEAD TO AFTER THE EXPRESSION
                    head = m.offset + length(m.match)

                    # WHAT KIND OF BLOCK IS IT?
                    if startswith(m.match, "@@CODEIMG")
                        # FIND THE FILE
                        name = m.captures[1]
                        incf = read(joinpath(cdir, "src/exgen/$name.jl"), String)
                        # WRITE IT APPROPRIATELY GUARDED
                        write(outf, "\n```julia\n$incf```\n")
                        # WRITE THE IMG APPROPRIATELY GUARDED
                        write(outf, "\n![](../exgen/out/$name.svg)\n")
                    elseif startswith(m.match, "@@CODE")
                        # FIND THE FILE
                        name = m.captures[1]
                        incf = read(joinpath(cdir, "src/exgen/$name.jl"), String)
                        # WRITE IT APPROPRIATELY GUARDED
                        write(outf, "\n```julia\n$incf```\n")
                    elseif startswith(m.match, "@@IMG")
                        # FIND THE IMAGE (FOR NOW ASSUMED SVG FORMAT)
                        name = m.captures[1]
                        # WRITE IT APPROPRIATELY GUARDED
                        write(outf, "\n![](../exgen/out/$name.svg)\n")
                    end
                end
                write(outf, inf[head:end])
            end
        end
    end
    println("       done ✅  [$(round(time()-start,digits=1))s]")
    println("Generating docs...")
    include("make.jl")
end # RUN PREMAKE
# reset the preview in case we keep using GPlot
continuous_preview(true)


text:
"""
Handling convergence ciretron for iterative algorithms

@author : Spencer Lyon <spencer.lyon@stern.nyu.edu>
@date : 2015-04-13 11:29:12

"""
module IterationManagers
import LinearAlgebra
using LinearAlgebra: norm
using Printf

abstract type IterationManager end
abstract type IterationState{T} end

export
# types
    IterationManager, IterationState, IterTolManager,  DefaultManager,
    IterManager, DefaultState, ExtraState,

# functions
    finished, update!, managed_iteration,
    verbose, print_now, display_iter, num_iter

include("states/default.jl")
include("states/extra.jl")

include("managers/itertol.jl")
include("managers/iter.jl")
include("managers/tol.jl")

include("api.jl")

end # module


text:
"""
Simulator for FMUs of
FMI 2.0 for Model Exchange Standard
"""

using Libdl         # For using dlopen, dlclose and so on
using LightXML      # For parsing XML files

export main

include("FMIWrapper.jl")

# Macro to identify logger library
macro libLogger()
    if Sys.iswindows()
        return joinpath(dirname(dirname(Base.source_path())),"bin", "win64", "logger.dll")
    elseif Sys.islinux()
        return joinpath(dirname(dirname(Base.source_path())),"bin", "unix64", "logger.so")
    else
        error("OS not supportet")
    end
end

"""
Parse modelDescription.xml
"""
function readModelDescription(pathToModelDescription::String)

    md = ModelDescription()

    if !isfile(pathToModelDescription)
        error("File $pathToModelDescription does not exist.")
    elseif basename(pathToModelDescription) != "modelDescription.xml"
        error("File name is not equal to \"modelDescription.xml\" but $(basename(pathToModelDescription))" )
    end

    # Parse modelDescription
    xdoc = parse_file(pathToModelDescription)

    # Check root tag
    xroot = root(xdoc)
    if name(xroot)!="fmiModelDescription"
        error("modelDescription.xml root is not \"fmiModelDescription\".")
    end

    # Get attributes of tag <fmiModelDescription>
    attributes = attributes_dict(xroot)
    try
        md.fmiVersion = attributes["fmiVersion"]
        md.modelName = attributes["modelName"]
        md.guid = attributes["guid"]

        md.description = get(attributes, "description", "")
        md.author = get(attributes, "author", "")
        md.version = get(attributes, "version", "")
        md.copyright = get(attributes, "copyright", "")
        md.license = get(attributes, "license", "")
        md.generationTool = get(attributes, "generationTool", "")
        md.generationDateAndTime = get(attributes, "generationDateAndTime", "")
        md.variableNamingConvention = NamingConvention(get(attributes, "variableNamingConvention", "flat"))
        md.numberOfEventIndicators = parse(Int, attributes["numberOfEventIndicators"])

        # Get attributes of tag <ModelExchange>
        elementModelExchange = find_element(xroot, "ModelExchange")

        # Get attributes of tag <CoSimulation>
        elementCoSimulation = find_element(xroot, "CoSimulation")

        try
            if elementModelExchange != nothing
                md.isModelExchange = true
                md.modelIdentifier = attribute(elementModelExchange, "modelIdentifier"; required=true)
            else
                md.isModelExchange = false
                error("While parsing modelDescription: Only ModelExchange is supported yet but
                    FMU does not support ModelExchange.")
            end

            if elementCoSimulation != nothing
                md.isCoSimulation = true
                md.modelIdentifier = attribute(elementModelExchange, "modelIdentifier"; required=true)
            else
                md.isCoSimulation = false
            end
        catch err
            if isa(err, LightXML.XMLAttributeNotFound)
                error("While parsing modelDescription: Non-optinal element \"modelIdentifier\" not found")
            else
                rethrow(err)
            end
        end

        if elementCoSimulation == nothing && elementModelExchange == nothing
            error("modelDescription.xml is missing ModelExchange and CoSimulation tags.")
        end

        # ToDo: Add stuff for tag UnitDefinitions

        # ToDo: Add stuff for tag TypeDefinitions

        # Get attributes of tag LogCategories
        elementLogCategories = find_element(xroot, "LogCategories")
        if elementLogCategories != nothing
            numCategories = 0
            for element in child_elements(elementLogCategories)
                numCategories += 1
            end
            md.logCategories = Array{LogCategory}(undef, numCategories)

            for (index, element) in enumerate(child_elements(elementLogCategories))
                tmp_name = attribute(element, "name"; required=false)
                tmp_description = attribute(element, "description"; required=false)
                if tmp_description == nothing
                    md.logCategories[index] = LogCategory(tmp_name)
                else
                    md.logCategories[index] = LogCategory(tmp_name, tmp_description)
                end
            end
        end

        # Get attributes of tag DefaultExperiment
        elementDefaultExperiment = find_element(xroot, "DefaultExperiment")
        if elementDefaultExperiment != nothing
            startTime = attribute(elementDefaultExperiment, "startTime"; required=false)
            if startTime != nothing
                startTime = parse(Float64, startTime)
            else
                startTime = 0
            end

            stopTime = attribute(elementDefaultExperiment, "stopTime"; required=false)
            if startTime != nothing
                stopTime = parse(Float64, stopTime)
            else
                stopTime = startTime+1
            end

            tolerance = attribute(elementDefaultExperiment, "tolerance"; required=false)
            if tolerance != nothing
                tolerance = parse(Float64, tolerance)
            else
                tolerance = 1e-6
            end

            stepSize = attribute(elementDefaultExperiment, "stepSize"; required=false)
            if stepSize != nothing
                stepSize = parse(Float64, stepSize)
            else
                stepSize = 2e-3
            end

            md.defaultExperiment = ExperimentData(startTime, stopTime, tolerance, stepSize)
        else
            md.defaultExperiment = ExperimentData()
        end

        # ToDo: Add stuff for tag VendorAnnotations

        # Get attributes of tag ModelVariables
        elementModelVariables = find_element(xroot, "ModelVariables")
        if elementModelVariables == nothing
            error("modelDescription.xml is missing ModelVariables tag")
        end

        numberOfVariables = 0
        for element in get_elements_by_tagname(elementModelVariables, "ScalarVariable")
            numberOfVariables += 1
        end
        scalarVariables = Array{ScalarVariable}(undef, numberOfVariables)

        for (index, element) in enumerate(get_elements_by_tagname(elementModelVariables, "ScalarVariable"))
            # Get attributes name, valueReference, variability, causality, initial
            tmp_name = attribute(element, "name"; required=true)
            tmp_valueReference = parse(Int, attribute(element, "valueReference"; required=true))

            tmp_description = attribute(element, "description"; required=false)
            if tmp_description == nothing
                tmp_description=""
            end
            tmp_variability = attribute(element, "variability"; required=false)
            if tmp_variability == nothing
                tmp_variability=""
            end
            tmp_causality = attribute(element, "causality"; required=false)
            if tmp_causality == nothing
                tmp_causality=""
            end
            tmp_initial = attribute(element, "initial"; required=false)
            if tmp_initial == nothing
                tmp_initial=""
            end
            tmp_canHandleMultipleSetPerTimelnstant = attribute(element, "canHandleMultipleSetPerTimelnstant"; required=false)
            if tmp_canHandleMultipleSetPerTimelnstant == nothing
                tmp_canHandleMultipleSetPerTimelnstant = false
            else
                tmp_canHandleMultipleSetPerTimelnstant = parse(Bool,
                    tmp_canHandleMultipleSetPerTimelnstant)
            end

            # Get child node for typeSpecificProperties
            tmp_typeSpecificProperties = nothing
            for child in child_elements(element)
                if is_elementnode(child)
                    if name(child)=="Real"
                        tmp_declaredType = "Real"
                        tmp_variableAttributes = RealAttributes()   # TODO implement
                        tmp_start = attribute(child, "start"; required=false)
                        if tmp_start == nothing
                            tmp_start = Float64(0)
                        else
                            tmp_start = parse(Float64, tmp_start)
                        end
                        tmp_derivative = attribute(child, "derivative"; required=false)
                        if tmp_derivative == nothing
                            tmp_derivative = UInt(0)
                        else
                            tmp_derivative = parse(UInt, tmp_derivative)
                        end
                        tmp_reinit = attribute(child, "reinit"; required=false)
                        if tmp_reinit == nothing
                            tmp_reinit = false
                        end
                        tmp_typeSpecificProperties = RealProperties(tmp_declaredType, tmp_variableAttributes, tmp_start, tmp_derivative, tmp_reinit)
                    elseif name(child)=="Integer"
                        tmp_declaredType = "Integer"
                        tmp_variableAttributes = IntegerAttributes()   # TODO implement
                        tmp_start = attribute(child, "start"; required=false)
                        if tmp_start == nothing
                            tmp_start = Int(0)
                        else
                            tmp_start = parse(Int, tmp_start)
                        end
                        tmp_typeSpecificProperties = IntegerProperties(tmp_declaredType, tmp_variableAttributes, tmp_start)
                    elseif name(child)=="Boolean"
                        tmp_declaredType = "Boolean"
                        tmp_start = attribute(child, "start"; required=false)
                        if tmp_start == nothing
                            tmp_start = false
                        else
                            tmp_start = parse(Bool, tmp_start)
                        end
                        tmp_typeSpecificProperties = BooleanProperties(tmp_declaredType, tmp_start)
                    elseif name(child)=="String"
                        error("Type \"String\" not implemented.")
                    else
                        error("Unknown type \"$(name(child))\" of ScalarVariable")
                    end
                end
            end

            scalarVariables[index] = ScalarVariable(tmp_name,
                tmp_valueReference, tmp_description, tmp_causality,
                tmp_variability, tmp_initial,
                tmp_canHandleMultipleSetPerTimelnstant,
                tmp_typeSpecificProperties)
        end
        md.modelVariables = scalarVariables

        # Get attributes of tag ModelStructure

    catch err
        if isa(err, KeyError)
            error("While parsing modelDescription: Non-optional element \"$(err.key)\" not found")
        else
            rethrow(err)
        end
    end

    # Free memory
    LightXML.free(xdoc)

    return md
end

function modelDescriptionToModelData(modelDescription::ModelDescription)

    modelData = ModelData()

    for var in modelDescription.modelVariables
        if typeof(var.typeSpecificProperties)==RealProperties
            if var.typeSpecificProperties.derivative > 0
                modelData.numberOfStates += 1
                modelData.numberOfDerivatives += 1
            end
            modelData.numberOfReals += 1
        elseif typeof(var.typeSpecificProperties)==IntegerProperties
            modelData.numberOfInts += 1
        elseif typeof(var.typeSpecificProperties)==BooleanProperties
            modelData.numberOfBools += 1
        elseif typeof(var.typeSpecificProperties)==StringProperties
            modelData.numberOfStrings += 1
        else
            error()
        end

    modelData.numberOfEventIndicators = modelDescription.numberOfEventIndicators
    end

    return modelData
end

"""
Helper function to initialize modelVars in simulationData with start values from
modelDescription and set valueReference and name.
"""
function initializeSimulationData(modelDescription::ModelDescription,
    modelData::ModelData)

    prevVars=0

    simulationData = SimulationData(modelData.numberOfReals,
        modelData.numberOfInts, modelData.numberOfBools,
        modelData.numberOfStrings, 0, modelData.numberOfEventIndicators)    # TODO Add number of enumerations here!

    i_real = i_int = i_bool = i_string = i_enumertion = 0

    # Fill simulation data with start value, value reference and name for all
    # scalar variables
    for scalarVar in modelDescription.modelVariables
        if typeof(scalarVar.typeSpecificProperties) == RealProperties
            i_real += 1
            simulationData.modelVariables.reals[i_real] =
                RealVariable(scalarVar.typeSpecificProperties.start,
                             scalarVar.valueReference,
                             scalarVar.name)
        elseif typeof(scalarVar.typeSpecificProperties) == IntegerProperties
            i_int += 1
            simulationData.modelVariables.ints[i_int] =
                IntVariable(scalarVar.typeSpecificProperties.start,
                            scalarVar.valueReference,
                            scalarVar.name)
        elseif typeof(scalarVar.typeSpecificProperties) == BooleanProperties
            i_bool += 1
            simulationData.modelVariables.bools[i_bool] =
                BoolVariable(scalarVar.typeSpecificProperties.start,
                             scalarVar.valueReference,
                             scalarVar.name)
        elseif typeof(scalarVar.typeSpecificProperties) == StringProperties
            i_string += 1
            simulationData.modelVariables.strings[i_string] =
                StringVariable(scalarVar.typeSpecificProperties.start,
                               scalarVar.valueReference,
                               scalarVar.name)
        elseif typeof(scalarVar.typeSpecificProperties) == EnumerationProperties
            error("Enumeration variables not implemeted!")
            #i_enumertion += 1
            # TODO add enumerations here
            #simulationData.modelVariables.enumerations[i_enumertion] =
            #    EnumerationVariable()
        else
            error("Unknown scalar variable type $(typeof(scalarVar.typeSpecificProperties)).")
        end
    end

    if i_real != modelData.numberOfReals
        error("Counted number of real scalar variables $i_real didn't matched expeted $(modelData.numberOfReals + modelData.numberOfEventIndicators)")
    elseif i_int != modelData.numberOfInts
        error("Counted number of integer scalar variables $i_int didn't matched expeted $(modelData.numberOfInts)")
    elseif i_bool != modelData.numberOfBools
        error("Counted number of boolean scalar variables $i_bool didn't matched expeted $(modelData.numberOfBools)")
    elseif i_string != modelData.numberOfStrings
        error("Counted number of string scalar variables $i_string didn't matched expeted $(modelData.numberOfStrings)")
    elseif i_enumertion != 0
        error("Counted number of enumeration scalar variables $i_enumertion didn't matched expeted 0")
    end

    return simulationData
end

"""
`loadFMU(pathToFMU::String, useTemp::Bool=false, overWriteTemp::Bool=true)`

Unzips an FMU and returns handle to dynamic library containing FMI functions.

## Example calls
```
julia> fmu=loadFMU("path/to/fmu/helloWorld.fmu")
```
"""
function loadFMU(pathToFMU::String, useTemp::Bool=false, overWriteTemp::Bool=true)
    # Create uninitialized FMU
    fmu=FMU()

    # Split path
    fmu.FMUPath = pathToFMU
    name = splitext(basename(pathToFMU))[1]

    # Create temp folder
    if useTemp
        fmu.tmpFolder = joinpath(tempdir(), string("FMU_", name, "_", floor(Int, 10000*rand())))
    else
        fmu.tmpFolder = joinpath(dirname(pathToFMU), string("FMU_", name))
    end
    if isdir(fmu.tmpFolder)
        if !overWriteTemp
            error("Folder $tmpFolder already exists but overwriting of temp
            folder is prohibited.")
        end
    end

    # unzip FMU to tmp folder
    my_unzip(pathToFMU, fmu.tmpFolder)

    # parse modelDescription.xml
    fmu.modelDescription = readModelDescription(joinpath(fmu.tmpFolder, "modelDescription.xml"))
    fmu.modelName = fmu.modelDescription.modelName
    fmu.instanceName = fmu.modelDescription.modelName
    if (fmu.modelDescription.isModelExchange)
        fmu.fmuType = modelExchange
    else
        error("FMU does not support modelExchange")
    end

    # pathToDLL
    if Sys.iswindows()
        pathToDLL = joinpath(fmu.tmpFolder, "binaries", "win$(Sys.WORD_SIZE)", string(name, ".dll"))
    elseif Sys.islinux()
        pathToDLL = joinpath(fmu.tmpFolder, "binaries", "linux$(Sys.WORD_SIZE)", string(name, ".so"))
    else
        error("OS not supportet!")
    end

    if !isfile(pathToDLL)
        if Sys.iswindows()
            error("No shared library found matching $(Sys.WORD_SIZE) bit Windows.")
        elseif Sys.islinux()
            error("No shared library found matching $(Sys.WORD_SIZE) bit Linux.")
        end
    end

    # Fill model data
    fmu.modelData = modelDescriptionToModelData(fmu.modelDescription)

    # Fill Simulation Data
    fmu.simulationData = initializeSimulationData(fmu.modelDescription, fmu.modelData)

    # Set default experiment
    fmu.experimentData = deepcopy(fmu.modelDescription.defaultExperiment)

    fmu.eventInfo = EventInfo()

    # Open result and log file
    fmu.csvFile = open("$(fmu.modelName)_results.csv", "w")
    fmu.logFile = open("$(fmu.modelName).log", "w")

    # load shared library with FMU
    # TODO export DL_LOAD_PATH="/usr/lib/x86_64-linux-gnu" on unix systems
    # push!(DL_LOAD_PATH, "/usr/lib/x86_64-linux-gnu") maybe???
    fmu.libHandle = dlopen(pathToDLL)

    # Load hared library with logger function
    fmu.libLoggerHandle = dlopen(@libLogger)
    fmi2CallbacLogger_Cfunc = dlsym(fmu.libLoggerHandle, :logger)
    # fmi2CallbacLogger_funcWrapC = @cfunction(fmi2CallbackLogger, Cvoid,
    #    (Ptr{Cvoid}, Cstring, Cuint, Cstring, Tuple{Cstring}))
    fmi2AllocateMemory_funcWrapC = @cfunction(fmi2AllocateMemory, Ptr{Cvoid}, (Csize_t, Csize_t))
    fmi2FreeMemory_funcWrapC = @cfunction(fmi2FreeMemory, Cvoid, (Ptr{Cvoid},))

    fmi2Functions = CallbackFunctions(
        #fmi2CallbacLogger_funcWrapC,       # Logger in Julia
        fmi2CallbacLogger_Cfunc,            # Logger in C
        fmi2AllocateMemory_funcWrapC,
        fmi2FreeMemory_funcWrapC,
        C_NULL,
        C_NULL)

    # Fill FMU with remaining data
    fmu.fmuResourceLocation = joinpath(string("file:///", fmu.tmpFolder), "resources")
    fmu.fmuGUID = fmu.modelDescription.guid
    fmu.fmiCallbackFunctions = fmi2Functions

    fmu.modelState = modelUninstantiated

    return fmu
end

"""
```
    unloadFMU(fmu::FMU, [deleteTmpFolder=true::Bool])
```
Unload FMU and if `deleteTmpFolder=true` remove tmp files.
"""
function unloadFMU(fmu::FMU, deleteTmpFolder=true::Bool)

    # unload FMU dynamic library
    dlclose(fmu.libHandle)

    # unload C logger
    dlclose(fmu.libLoggerHandle)

    # delete tmp folder
    if deleteTmpFolder
        rm(fmu.tmpFolder, recursive=true, force=true);
    end

    # Close result and log file
    close(fmu.csvFile)
    close(fmu.logFile)
end

"""
Helper function to unzip file using installed tools.
Overwrites `destiantionDir`
"""
function my_unzip(target::String, destinationDir::String)

    if !isfile(target)
        error("Could not find file \"$target\"")
    end

    if !ispath(destinationDir)
        mkpath(destinationDir)
    end

    try
        #use unzip
        run(Cmd(`unzip -qo $target`, dir = destinationDir))
        println("Extracted FMU to $destinationDir")
    catch
        try
            #use 7-zip
            run(Cmd(`"C:\\Program Files\\7-Zip\\7z.exe" x $target -aoa -o$destinationDir`));
            println("Extracted FMU to $destinationDir")
        catch
            error("Could not unzip file \"$target\"")
        end
    end
end

"""
Gets values of all states of a FMU.
Helper function for main.
"""
function getContinuousStates!(fmu::FMU)

    if fmu.modelState != modelContinuousTimeMode
        error("Function call only allowed in modelContinuousTimeMode but model is in mode: $(fmu.modelState)")
    end

    states = Array{Float64}(undef,fmu.modelData.numberOfStates)
    fmi2GetContinuousStates!(fmu, states)

    for i=1:fmu.modelData.numberOfStates
        fmu.simulationData.modelVariables.reals[i].value = states[i]
    end
end

function getDerivatives!(fmu::FMU)

    if fmu.modelState != modelContinuousTimeMode
        error("Function call only allowed in modelContinuousTimeMode but model is in mode: $(fmu.modelState)")
    end

    derivatives = Array{Float64}(undef,fmu.modelData.numberOfStates)
    fmi2GetDerivatives!(fmu, derivatives)

    for i=1:fmu.modelData.numberOfDerivatives
        fmu.simulationData.modelVariables.reals[i+fmu.modelData.numberOfStates].value = derivatives[i]
    end
end

"""
Sets continous states in `FMU`
"""
function setContinuousStates!(fmu::FMU)

    if fmu.modelState != modelContinuousTimeMode
        error("Function call only allowed in modelContinuousTimeMode but model is in mode: $(fmu.modelState)")
    end

    states = Array{Float64}(undef,fmu.modelData.numberOfStates)

    for i=1:fmu.modelData.numberOfStates
        states[i] = fmu.simulationData.modelVariables.reals[i].value
    end

    fmi2SetContinuousStates(fmu, states)

end

"""
Get eventIndicators from FMU and update saved values.
"""
function getEventIndicators!(fmu::FMU)

    fmi2GetEventIndicators!(fmu, fmu.simulationData.eventIndicators)
end

function getVariable!(fmu::FMU, variable::RealVariable)

    fmi2GetReal!(fmu, [variable.valueReference], [variable.value])
    return variable.value
end

function getVariable!(fmu::FMU, variable::IntVariable)

    fmi2GetInteger!(fmu, [variable.valueReference], [variable.value])
    return variable.value
end

function getVariable!(fmu::FMU, variable::BoolVariable)

    fmi2GetBoolean!(fmu, [variable.valueReference], [variable.value])
    return variable.value
end

function getVariable!(fmu::FMU, variable::StringVariable)

    fmi2GetString!(fmu, [variable.valueReference], [variable.value])
    return variable.value
end

function getAllVariables!(fmu::FMU)

    for realVar in fmu.simulationData.modelVariables.reals
        realVar.value = getVariable!(fmu, realVar)
    end
    for intVar in fmu.simulationData.modelVariables.ints
        intVar.value = getVariable!(fmu, intVar)
    end
    for boolVar in fmu.simulationData.modelVariables.bools
        boolVar.value = getVariable!(fmu, boolVar)
    end
    for stringVar in fmu.simulationData.modelVariables.strings
        stringVar.value = getVariable!(fmu, stringVar)
    end
end

function setTime!(fmu::FMU, time::Float64, saveLastStepTime=true::Bool)

    if saveLastStepTime
        fmu.simulationData.lastStepTime = fmu.simulationData.time
    end
    fmu.simulationData.time = time
    fmi2SetTime(fmu, fmu.simulationData.time)
end

"""
Checks if signs of two arrays are different component-wise.
Helper function for bisection.
"""
function arrayDiffSign(array1::Array{Float64,1}, array2::Array{Float64,1})

    if length(array1) != length(array2)
        throw(DimensionMismatch("Left and right array of event indicators have different sizes."))
    end

    for i in 1:length(array1)
        if sign(array1[i]) != sign(array2[i])
            return true
        end
    end

    return false
end

"""
Find event time with bisection method for `eventIndicators` for given `fmu`.
"""
function findEvent(fmu::FMU)

    leftTime = fmu.simulationData.lastStepTime
    rightTime = fmu.simulationData.time

    leftEventIndicators = copy(fmu.simulationData.eventIndicators)
    getEventIndicators!(fmu)
    rightEventIndicators = copy(fmu.simulationData.eventIndicators)

    # Check if there are any events
    if !arrayDiffSign(leftEventIndicators, rightEventIndicators)
        return (false, 0)
    end

    steps = 0
    minimumStepSize = 1e-8      # TODO: Read mimimumStepSize from fmu experiment data
    maxSteps = ceil(log2((rightTime-leftTime)/minimumStepSize)) + 1
    centerTime = 0

    while rightTime - leftTime > minimumStepSize && steps < maxSteps
        steps += 1

        # Evaluate eventIndicators in center of intervall
        centerTime = 0.5*(rightTime - leftTime)
        setTime!(fmu, centerTime, false)
        getEventIndicators!(fmu)
        centerEventIndicators = copy(fmu.simulationData.eventIndicators)     # TODO Do I need to copy here?

        # TODO Check what happens when event is on leftTime, centerTime or rightTime
        # Check for event in first half of intervall [leftTime, centerTime]
        if arrayDiffSign(leftEventIndicators, centerEventIndicators)
            rightTime = centerTime
            rightEventIndicators = centerEventIndicators        # This does not copy memory, right?

        # Check for event in second half of intervall [centerTime, rightTime]
        else
            leftTime = centerTime
            leftEventIndicators = centerEventIndicators
        end
    end

    if steps == maxSteps
        error("Event was not found in maximum number of Steps!")
    end

    return (true, centerTime)
end

function writeNamesToCSV(fmu::FMU)

    write(fmu.csvFile, "\"time\"")
    for realVar in fmu.simulationData.modelVariables.reals
        write(fmu.csvFile, ",\"$(realVar.name)\"")
    end
    for intVar in fmu.simulationData.modelVariables.ints
        write(fmu.csvFile, ",\"$(intVar.name)\"")
    end
    for boolVar in fmu.simulationData.modelVariables.bools
        write(fmu.csvFile, ",\"$(boolVar.name)\"")
    end
    for stringVar in fmu.simulationData.modelVariables.strings
        write(fmu.csvFile, ",\"$(stringVar.name)\"")
    end
    write(fmu.csvFile, "\r\n")
end

function writeValuesToCSV(fmu::FMU)

    write(fmu.csvFile, "$(fmu.simulationData.time)")
    for realVar in fmu.simulationData.modelVariables.reals
        write(fmu.csvFile, ",$(realVar.value)")
    end
    for intVar in fmu.simulationData.modelVariables.ints
        write(fmu.csvFile, ",$(intVar.value)")
    end
    for boolVar in fmu.simulationData.modelVariables.bools
        write(fmu.csvFile, ",$(Int(boolVar.value))")
    end
    for stringVar in fmu.simulationData.modelVariables.strings
        write(fmu.csvFile, ",$(stringVar.value)")
    end
    write(fmu.csvFile,"\r\n")
end

"""
Main function to simulate a FMU
"""
function main(pathToFMU::String)
    # load FMU
    fmu = loadFMU(pathToFMU)
    writeNamesToCSV(fmu)

    try
        # Instantiate FMU
        fmi2Instantiate!(fmu)
        fmu.modelState = modelInstantiated

        # Set debug logging to true for all categories
        fmi2SetDebugLogging(fmu, true)

        # Get types platform
        typesPlatform = fmi2GetTypesPlatform(fmu)
        println("typesPlatform: $typesPlatform")

        # Get version of fmi
        fmiVersion = fmi2GetVersion(fmu)
        println("FMI version: $fmiVersion")

        # Set up experiment
        fmi2SetupExperiment(fmu, 0)

        # Set start time
        setTime!(fmu, fmu.experimentData.startTime, true)
        nextTime = fmu.experimentData.stopTime

        # Set initial variables with intial="exact" or "approx"

        # Initialize FMU
        fmi2EnterInitializationMode(fmu)

        # Exit Initialization
        fmi2ExitInitializationMode(fmu)

        # Event iteration
        fmu.eventInfo.newDiscreteStatesNeeded = true
        while fmu.eventInfo.newDiscreteStatesNeeded
            fmi2NewDiscreteStates!(fmu)
            if fmu.eventInfo.terminateSimulation
                error("FMU was terminated in Event at time $(fmu.simulationData.time)")
            end
        end
        # Initialize event indicators
        getEventIndicators!(fmu)

        # Enter Continuous time mode
        fmi2EnterContinuousTimeMode(fmu)

        # retrieve initial states
        getContinuousStates!(fmu)

        # retrive solution
        getAllVariables!(fmu)       # TODO Is not returning der(x) correctly
                                    # Needs to call fmi2GetXXX of course...
        writeValuesToCSV(fmu)

        # Iterate with explicit euler method
        k = 0
        k_max = 1000
        while (fmu.simulationData.time < fmu.experimentData.stopTime) && (k < k_max)
            k += 1
            getDerivatives!(fmu)

            # Compute next step size
            if fmu.eventInfo.nextEventTimeDefined
                h = min(fmu.experimentData.stepSize, fmu.eventInfo.nextEventTime - fmu.simulationData.time)
            else
                h = min(fmu.experimentData.stepSize, fmu.experimentData.stopTime - fmu.simulationData.time)
            end

            # Update time
            setTime!(fmu, fmu.simulationData.time + h)

            # Set states and perform euler step (x_k+1 = x_k + d/dx x_k*h)
            for i=1:fmu.modelData.numberOfStates
                fmu.simulationData.modelVariables.reals[i].value = fmu.simulationData.modelVariables.reals[i].value + h*fmu.simulationData.modelVariables.reals[i+fmu.modelData.numberOfStates].value
            end
            setContinuousStates!(fmu)

            # Get event indicators and check for events

            # Inform the model abaut an accepted step
            (enterEventMode, terminateSimulation) = fmi2CompletedIntegratorStep(fmu, true)
            if enterEventMode
                error("Should now enter Event mode...")
            end

            if terminateSimulation
                error("Solution got terminated bevore reaching end time.")
            end

            # save results
            getAllVariables!(fmu)
            writeValuesToCSV(fmu)

            # Handle events
        end

        # Terminate Simulation
        fmi2Terminate(fmu)

        # Free FMU
        # ToDo: Fix function
        #fmi2FreeInstance(fmu)
    finally
        # Unload FMU
        println("Unload FMU")
        unloadFMU(fmu)
    catch
        rethrow()
    end
end

max_stars_repo_name: cuihantao/JuliaFMI

max_issues_repo_name: cuihantao/JuliaFMI

max_forks_repo_name: cuihantao/JuliaFMI


text:
using CUDA
using Distributions: pdf
using Random
using Flux
using BSON
using AbstractTrees

RLBase.update!(p::RandomPolicy, x) = nothing

Random.rand(s::MultiContinuousSpace{<:CuArray}) = rand(CUDA.CURAND.generator(), s)

# avoid fallback silently
Flux.testmode!(p::AbstractPolicy, mode = true) =
    @error "someone forgets to implement this method!!!"

function Flux.testmode!(p::RandomStartPolicy, mode = true)
    testmode!(p.policy, mode)
end

Base.show(io::IO, p::AbstractPolicy) =
    AbstractTrees.print_tree(io, StructTree(p), get(io, :max_depth, 10))

AbstractTrees.children(t::StructTree{<:AbstractEnv}) = ()

AbstractTrees.printnode(io::IO, t::StructTree{<:AbstractEnv}) = print(
    io,
    "$(RLBase.get_name(t.x)): $(join([f(t.x) for f in RLBase.get_env_traits()], ","))",
)

function save(f::String, p::AbstractPolicy)
    policy = cpu(p)
    BSON.@save f policy
end

function load(f::String, ::Type{<:AbstractPolicy})
    BSON.@load f policy
    policy
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractEncounterModelImpl

export AbstractEncounterModel

abstract AbstractEncounterModel

end

module AbstractEncounterModelInterfaces

export
       # common
       generateEncounter,
       getInitialState,
       getNextCommand,
       getTrajectory,

       # CorrAEM
       generateEncountersToFile,
       setInitialDistributions,
       validate,
       getInitialSample

# common

function generateEncounter()
end

function getInitialState()
end

function getNextCommand()
end

function getTrajectory()
end

# CorrAEM

function generateEncountersToFile()
end

function setInitialDistributions()
end

function validate()
end

function getInitialSample()
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
@discord_object struct Activity
    name::String
    type::ActivityType.ActivityTypeEnum
    url::String
    created_at::Union{Int, DateTime}
    timestamps::ActivityTimestamps
    application_id::Snowflake
    details::String
    state::String
    emoji::ActivityEmoji
    party::ActivityParty
    assets::ActivityAssets
    secrets::ActivitySecrets
    instance::Bool
    flags::Int64
end


text:
module Types

using Base: UUID

#using Pkg.Types: Context
#uuid(ctx::Context,name::AbstractString) = ctx.env.project.deps[name]

function uuid(name::String)
    for key in keys(Base.loaded_modules)
        if key.name==name
            return key.uuid
        end
    end
end

uuid(name::Symbol) = uuid(String(name))
uuid(name::Module) = uuid(nameof(name))
#uuid(name::Union{Module,Symbol,AbstractString}) = uuid(Context(),name)

# uuid(ctx::Context,name::Symbol) = uuid(ctx,String(name))
# uuid(ctx::Context,name::Module) = uuid(ctx,nameof(name))
# uuid(name::Union{Module,Symbol,AbstractString}) = uuid(Context(),name)

const CONFIG_DIR = homedir() * "/.demenet/"

function setconfigdir(dir::String)
    global CONFIG_DIR = dir
end

demefile(uuid::UUID) = CONFIG_DIR * "/demes/$uuid"
keydir(uuid::UUID) = CONFIG_DIR * "/keys/$uuid/"
datadir(uuid::UUID) = CONFIG_DIR * "/data/$uuid/"
#datadir(uuid::UUID,app::AbstractString) = CONFIG_DIR * "/apps/$app/$uuid/"

abstract type AbstractID end

struct ID <: AbstractID
    id::Union{BigInt,Nothing}
end

function Base.string(id::ID; length=nothing, kwargs...)
    str = string(id.id; kwargs...)
    if length!=nothing
        len = Base.length(str)
        @assert len<=length
        pre = "0"^(length-len)
    else
        pre = ""
    end
    return "$pre$str"
end

ID(str::AbstractString; kwargs...) = ID(parse(BigInt,str; kwargs...))

import Base.Dict

Dict(id::ID) = Dict("id"=>string(id,base=16))

ID(dict::Dict) = ID(dict["id"],base=16)

import Base.Vector

Vector{UInt8}(id::ID; kwargs...) = Vector{UInt8}(string(id; kwargs...))
ID(bytes::Vector{UInt8}; kwargs...) = ID(String(copy(bytes)); kwargs...)

Base.:(==)(a::ID,b::ID) = a.id==b.id
Base.hash(a::ID,h::UInt) = hash(a.id,hash(:ID,h))
Base.in(a::ID,b::ID) = a==b

struct DemeID <: AbstractID
    uuid::UUID
    id::ID
end

DemeID(uuid::UUID,id::BigInt) = DemeID(uuid,ID(id))

Base.:(==)(a::DemeID,b::DemeID) = a.uuid==b.uuid && a.id==b.id
Base.hash(a::DemeID,h::UInt) = hash(a.id,hash(a.uuid,hash(:DemeID, h)))

# Do I need a base hash also here?

Base.in(a::DemeID,b::DemeID) = a==b

### The Notary and Cypher could actually be defined with package where the implementation sits in. I could introduce AbstractNotary, AbstractCypher and just add a constructor CypherSuite{UUID}(AbstractNotary) to give me the right definition.

### Could have multiple signatures as one wishes

end


text:
__precompile__(true)

module DirLang

function test()
	fields, lines = parse_dir("res/DirLang/Hello World")

	for i in lines
		run_line(fields, i)
	end
end

function parse_dir(path::String)
	types::Dict{String, Type} = Dict{String, Type}([
		"Int" 		=> Int,
		"Float" 	=> Float64,
		"String" 	=> String
	])
	t::String = readdir(joinpath(path, "__Type"))[1]
	if t == "Func"
		return parse_func(path)
	elseif haskey(types, t)
		v::String = readdir(joinpath(path, "__Val"))[1]
		t == "String" && return v
		return parse(types[t], v)
	end
end

function parse_func(path::String)
	codelines::Vector{Int},
	fields::Dict{String, Any} =
	let a::Vector{String} = readdir(path)
		sort!([parse(Int, i) for i in a
			if mapreduce(isnumeric, (x, y) -> x && y, i)]),

		Dict{String, Any}([i[2:end] => parse_dir(joinpath(path, i))
			for i in a if i[1] == '_' && i[2] != '_'])
	end

	lines::Vector{Tuple{String, Vector{String}}} = Vector{Tuple{String, Vector{String}}}(undef, 0)

	for i in codelines
		f = readdir(joinpath(path, string(i)))[1]
		argsind = let
			sort!([parse(Int, i) for i in readdir(joinpath(path, string(i), f, "__Args"))
				if mapreduce(isnumeric, (x, y) -> x && y, i)])
		end
		args = Vector{String}()
		for j in argsind
			push!(args, readdir(joinpath(path, string(i), f, "__Args", string(j)))[1])
		end
		push!(lines, (f, args))
	end

	return fields, lines
end

function run_line(fields::Dict{String, Any}, line::Tuple{String, Vector{String}})
	funcs = Dict([
		"print" => print
		"println" => println
	])
	args = [fields[i] for i in line[2]]
	funcs[line[1]](args...)
end

end


text:
# Indexing
# ========
#
# Indexing methods for biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# Getting and setting elements in a biological sequence.

function Base.checkbounds(seq::BioSequence, range::UnitRange)
    if 1 ≤ range.start && range.stop ≤ lastindex(seq)
        return true
    end
    throw(BoundsError(seq, range))
end

function Base.checkbounds(seq::BioSequence, locs::AbstractVector{Bool})
    if length(seq) == length(locs)
        return true
    end
    throw(BoundsError(seq, locs))
end

function Base.checkbounds(seq::BioSequence, locs::AbstractVector)
    for i in locs
        checkbounds(seq, i)
    end
    return true
end

function checkdimension(from::Integer, to::Integer)
    if from == to
        return true
    end
    throw(DimensionMismatch(string(
        "attempt to assign ",
        from, " elements to ",
        to,   " elements")))
end

function checkdimension(seq::BioSequence, locs::AbstractVector)
    return checkdimension(length(seq), length(locs))
end

function checkdimension(seq::BioSequence, locs::AbstractVector{Bool})
    return checkdimension(length(seq), sum(locs))
end

# creates a bit mask for given number of bits `n`
mask(n::Integer) = (UInt64(1) << n) - 1
mask(::Type{A}) where {A<:Alphabet} = mask(bitsof(A))

# assumes `i` is positive and `bitsof(A)` is a power of 2
@inline function bitindex(seq::BioSequence{A}, i::Integer) where {A}
    return BitIndex((i + first(seq.part) - 2) << trailing_zeros(bitsof(A)))
end

@inline function inbounds_getindex(seq::BioSequence{A}, i::Integer) where {A}
    j = bitindex(seq, i)
    @inbounds return decode(A, (seq.data[index(j)] >> offset(j)) & mask(A))
end

Base.getindex(seq::BioSequence, part::UnitRange) = BioSequence(seq, part)
Base.view(seq::BioSequence, part::UnitRange) = BioSequence(seq, part)

function Base.setindex!(seq::BioSequence{A},
                        other::BioSequence{A},
                        locs::AbstractVector{<:Integer}) where {A}
    checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    for (i, x) in zip(locs, other)
        unsafe_setindex!(seq, x, i)
    end
    return seq
end

function Base.setindex!(seq::BioSequence{A},
                        other::BioSequence{A},
                        locs::UnitRange{<:Integer}) where {A}
    checkbounds(seq, locs)
    checkdimension(other, locs)
    return copyto!(seq, locs.start, other, 1)
end

function Base.setindex!(seq::BioSequence{A},
                        other::BioSequence{A},
                        locs::AbstractVector{Bool}) where {A}
    checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    i = j = 0
    while true
        i = findnext(locs, i + 1)
        if i === nothing
            break
        end
        unsafe_setindex!(seq, other[j+=1], i)
    end
    return seq
end

function Base.setindex!(seq::BioSequence{A}, other::BioSequence{A}, ::Colon) where {A}
    return setindex!(seq, other, 1:lastindex(seq))
end

function Base.setindex!(seq::BioSequence, x, i::Integer)
    checkbounds(seq, i)
    orphan!(seq)
    return unsafe_setindex!(seq, x, i)
end

function Base.setindex!(seq::BioSequence{A}, x, locs::AbstractVector{<:Integer}) where {A}
    checkbounds(seq, locs)
    bin = enc64(seq, x)
    orphan!(seq)
    for i in locs
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

function Base.setindex!(seq::BioSequence{A}, x, locs::AbstractVector{Bool}) where {A}
    checkbounds(seq, locs)
    bin = enc64(seq, x)
    orphan!(seq)
    i = j = 0
    while true
        i = findnext(locs, i + 1)
        if i === nothing
            break
        end
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

Base.setindex!(seq::BioSequence{A}, x, ::Colon) where {A} = setindex!(seq, x, 1:lastindex(seq))

# this is "unsafe" because of no bounds check and no orphan! call
@inline function unsafe_setindex!(seq::BioSequence{A}, x, i::Integer) where {A}
    bin = enc64(seq, x)
    return encoded_setindex!(seq, bin, i)
end

@inline function encoded_setindex!(seq::BioSequence{A}, bin::UInt64, i::Integer) where {A}
    j, r = bitindex(seq, i)
    @inbounds seq.data[j] = (bin << r) | (seq.data[j] & ~(mask(A) << r))
    return seq
end


text:
module PassiveCable

using Unitful: uconvert, NoUnits

export cable_normalize, cable, rallpack1

# Compute solution g(x, t) to
#
#     ∂²g/∂x² - g - ∂g/∂t = 0
#
# on [0, L] × [0,∞), subject to:
#
#     g(x, 0) = 0
#     ∂g/∂x (0, t) = 1
#     ∂g/∂x (L, t) = 0
#
# Parameters:
#     x, t, L:  as described above
#   tol:  absolute error tolerance in result
#
# Return:
#     g(x, t)
#
# TODO: verify correctness when L≠1

function cable_normalized(x::Float64, t::Float64, L::Float64; tol=1e-8)
    if t<=0
        return 0.0
    else
        ginf = -cosh(L-x)/sinh(L)
        sum = exp(-t/L)
        Ltol = L*tol

        for k = Iterators.countfrom(1)
            a = k*pi/L
            b = exp(-t*(1+a^2))

            sum += 2*b*cos(a*x)/(1+a^2)
            resid_ub = b/(a^3*t)

            if resid_ub<Ltol
                break
            end
        end
        return ginf+sum/L;
     end
end

# Compute solution f(x, t) to
#
#     λ²∂²f/∂x² - f - τ∂f/∂t = 0
#
# on [0, L] x [0, ∞), subject to:
#
#     f(x, 0) = V
#     ∂f/∂x (0, t) = I·r
#     ∂f/∂x (L, t) = 0
#
# where:
#
#     λ² = 1/(r·g)   length constant
#     τ  = r·c       time constant
#
# In the physical model, the parameters correspond to the following:
#
#     L:  length of cable
#     r:  linear axial resistivity
#     g:  linear membrane conductivity
#     c:  linear membrane capacitance
#     V:  membrane reversal potential
#     I:  injected axial current on the left end (x = 0) of the cable.
#
# Note that r, g and c are specific 1-d quantities that differ from
# the cable resistivity r_L, specific membrane conductivity ḡ and
# specific membrane capacitance c_m as used elsewhere. If the
# cross-sectional area is A and cable circumference is f, then
# these quantities are related by:
#
#     r = r_L/A
#     g = ḡ·f
#     c = c_m·f
#
# Parameters:
#     x:  displacement along cable
#     t:  time
#     L, lambda, tau, r, V, I:  as described above
#   tol:  absolute error tolerance in result
#
# Return:
#     computed potential at (x,t) on cable.

function cable(x, t, L, lambda, tau, r, V, I; tol=1e-8)
    scale = I*r*lambda;
    if scale == 0
        return V
    else
        tol_n = abs(tol/scale)
        return scale*cable_normalized(uconvert(NoUnits, x/lambda), uconvert(NoUnits, t/tau), uconvert(NoUnits, L/lambda), tol=tol_n) + V
    end
end

# Rallpack 1 test
#
# One sided cable model with the following parameters:
#
#     RA = 1 Ω·m    bulk axial resistivity
#     RM = 4 Ω·m²   areal membrane resistivity
#     CM = 0.01 F/m²  areal membrane capacitance
#     d  = 1 µm     cable diameter
#     EM = -65 mV   reversal potential
#     I  = 0.1 nA   injected current
#     L  = 1 mm     cable length.
#
# (This notation aligns with that used in the Rallpacks paper.)
#
# Note that the injected current as described in the Rallpack model
# is trans-membrane, not axial. Consequently we need to swap the
# sign on I when passing to the cable function.
#
# Parameters:
#     x:  displacement along cable [m]
#     t:  time [s]
#   tol:  absolute tolerance for reported potential.
#
# Return:
#     computed potential at (x,t) on cable.

function rallpack1(x, t; tol=1e-8)
    RA = 1
    RM = 4
    CM = 1e-2
    d  = 1e-6
    EM = -65e-3
    I  = 0.1e-9
    L  = 1e-3

    r = 4*RA/(pi*d*d)
    lambda = sqrt(d/4 * RM/RA)
    tau = CM*RM

    return cable(x, t, L, lambda, tau, r, EM, -I, tol=tol)
end

end # module PassiveCable

max_stars_repo_name: kanzl/arbor

max_issues_repo_name: kanzl/arbor

max_forks_repo_name: kanzl/arbor


text:
immutable Uniform <: ContinuousUnivariateDistribution
    a::Float64
    b::Float64
    function Uniform(a::Real, b::Real)
	a < b || error("a < b required for range [a, b]")
	new(float64(a), float64(b))
    end
    Uniform() = new(0.0, 1.0)
end

@_jl_dist_2p Uniform unif

min(d::Uniform) = d.a

max(d::Uniform) = d.b

entropy(d::Uniform) = log(d.b - d.a)

insupport(d::Uniform, x::Real) = d.a <= x <= d.b

kurtosis(d::Uniform) = -6.0 / 5.0

mean(d::Uniform) = (d.a + d.b) / 2.0

median(d::Uniform) = (d.a + d.b) / 2.0

function mgf(d::Uniform, t::Real)
	a, b = d.a, d.b
	return (exp(t * b) - exp(t * a)) / (t * (b - a))
end

function cf(d::Uniform, t::Real)
	a, b = d.a, d.b
	return (exp(im * t * b) - exp(im * t * a)) / (im * t * (b - a))
end

mode(d::Uniform) = d.a
modes(d::Uniform) = error("The uniform distribution has no modes")

rand(d::Uniform) = d.a + (d.b - d.a) * rand()

skewness(d::Uniform) = 0.0

function var(d::Uniform)
	w = d.b - d.a
	return w * w / 12.0
end

# fit model

function fit_mle{T <: Real}(::Type{Uniform}, x::Vector{T})
    if isempty(x)
        throw(ArgumentError("x cannot be empty."))
    end

    xmin = xmax = x[1]
    for i = 2:length(x)
        xi = x[i]
        if xi < xmin
            xmin = xi
        elseif xi > xmax
            xmax = xi
        end
    end

    Uniform(xmin, xmax)
end


text:
@testset "add" begin
    x1 = [1,2,3]
    x2 = Var([1,2,3])
    y = x1 + x2
    @test y.data == [2,4,6]

    x = randn(3, 3)
    y = randn(3, 3)
    f = x -> x + y
    @test gradcheck(f, x)

    x = Var(randn(3, 3))
    y = randn(3)
    @test gradcheck(f, x)
end

@testset "mul" begin
    x1 = [1, 2, 3]
    x2 = Var([1, 2, 3])
    y = x1 .* x2
    @test y.data == [1, 4, 9]

    x = randn(3, 3)
    y = randn(3, 3)
    f = x -> x .* y
    @test gradcheck(f, x)

    y = randn(3, 1)
    @test gradcheck(f, x)

    y = randn(3)
    @test gradcheck(f, x)
    f = y -> x .* y
    @test gradcheck(f, x)
end

@testset "div" begin
    x1 = [1, 2, 3]
    x2 = Var([1, 2, 3])
    y = x1 / x2
    @test y.data == [1, 1, 1]

    x = randn(3, 3)
    y = randn(3, 3)
    f = x -> x ./ y
    @test gradcheck(f, x)
    y = randn(3)
    @test gradcheck(f, x)
end

@testset "pow" begin
    x = Var(2)
    @test (x^10).data == [1024]
    y = x^2
    gradient!(y)
    @test x.grad.data == [4.0]
end


text:
immutable SubSpace{DS,IT,T,DD,dim}<: Space{T,DD,dim}
    space::DS
    indexes::IT
end

SubSpace{T,DD,dim}(sp::Space{T,DD,dim},kr) =
    SubSpace{typeof(sp),typeof(kr),T,DD,dim}(sp,kr)

SubSpace(sp::SubSpace,kr) = SubSpace(sp.space,sp.indexes[kr])

domain(DS::SubSpace) = domain(DS.space)
dimension(sp::SubSpace) = length(sp.indexes)

|(sp::Space,kr::Union{AbstractCount,Range}) = SubSpace(sp,kr)

spacescompatible{DS,IT,T,DD,d}(S1::SubSpace{DS,IT,T,DD,d},S2::SubSpace{DS,IT,T,DD,d}) =
    spacescompatible(S1.space,S2.space) && S1.indexes == S2.indexes

=={DS,IT,T,DD,d}(S1::SubSpace{DS,IT,T,DD,d},S2::SubSpace{DS,IT,T,DD,d}) =
    S1.space == S2.space && S1.indexes == S2.indexes

canonicalspace(a::SubSpace) = a.space

setdomain(DS::SubSpace,d::Domain) = SubSpace(setdomain(DS.space,d),DS.indexes)

Conversion{S<:Space,IT,T,DD,d}(a::SubSpace{S,IT,T,DD,d},b::S) =
    ConcreteConversion(a,b)
bandinds{S,T,DD,d}(C::ConcreteConversion{SubSpace{S,UnitCount{Int},T,DD,d},S}) =
    1-first(domainspace(C).indexes),0

getindex{S,IT,T,DD,d}(C::ConcreteConversion{SubSpace{S,IT,T,DD,d},S},
                   k::Integer,j::Integer) =
    domainspace(C).indexes[j]==k?one(eltype(C)):zero(eltype(C))

getindex{IT,DS,T,DD,d}(E::ConcreteEvaluation{SubSpace{DS,IT,T,DD,d},Bool},k::Integer) =
    Evaluation(E.space.space,E.x,E.order)[E.space.indexes[k]]
getindex{IT,DS,T,DD,d}(E::ConcreteEvaluation{SubSpace{DS,IT,T,DD,d}},k::Integer) =
    Evaluation(E.space.space,E.x,E.order)[E.space.indexes[k]]
getindex{IT,DS,T,DD,d}(E::ConcreteEvaluation{SubSpace{DS,IT,T,DD,d}},kr::Range) =
    Evaluation(E.space.space,E.x,E.order)[E.space.indexes[kr]]

function conversion_rule(a::SubSpace,b::SubSpace)
     if a==b
        a
    else
        NoSpace()
    end
end
# return the space that has banded Conversion to the other
function conversion_rule(a::SubSpace,b::Space)
    if a.space==b
        a  # we can write droping coefficients as a banded operator
    else
        NoSpace()
    end
end

function coefficients(v::Vector,sp::SubSpace,dropsp::SubSpace)
    if sp == dropsp
        v
    else
        coefficients(v,sp,canonicalspace(sp),dropsp)
    end
end

## points

## transform
function transform(sp::SubSpace,vals::Vector)
    ret=transform(sp.space,vals)
    coefficients(ret,sp.space,sp)
end

itransform(sp::SubSpace,cfs::Vector) =
    itransform(sp.space,coefficients(cfs,sp,sp.space))

points(sp::SubSpace,n) = points(sp.space,n)

for TYP in (:SumSpace,:PiecewiseSpace,:Space) # Resolve conflict
    @eval begin
        function coefficients(v::Vector,sp::$TYP,dropsp::SubSpace)
            n=length(v)
            if sp==dropsp.space
                ret = Array(eltype(v),0)
                for k in dropsp.indexes
                    if k > n
                        return ret
                    end
                    push!(ret,v[k])
                end
            else
                coefficients(v,sp,canonicalspace(dropsp),dropsp)
            end
        end

        function coefficients(v::Vector,dropsp::SubSpace,sp::$TYP)
            if sp==dropsp.space
                ret = zeros(eltype(v),dropsp.indexes[length(v)])
                for k = eachindex(v)
                    ret[dropsp.indexes[k]] = v[k]
                end
                ret
            else
                coefficients(v,dropsp,canonicalspace(dropsp),sp)
            end
        end
    end
end

## ProductFUn

# values{S<:SubSpace,V<:SubSpace}(f::ProductFun{S,V})=values(ProductFun(f,space(f,1).space,space(f,2).space))
# values{S<:SubSpace}(f::ProductFun{S})=values(ProductFun(f,space(f,1).space,space(f,2)))
#
#
# function coefficients{n,DS,TT,DD}(f::ProductFun{SubSpace{n,1,DS,TT,DD,1}},ox::Space,oy::Space)
#     T=eltype(f)
#     m=size(f,1)
#     A=[pad!(coefficients(fx,ox),m+n) for fx in f.coefficients]
#     B=hcat(A...)::Array{T,2}
#     for k=1:size(B,1)
#         ccfs=coefficients(vec(B[k,:]),space(f,2),oy)
#         if length(ccfs)>size(B,2)
#             B=pad(B,size(B,1),length(ccfs))
#         end
#         B[k,1:length(ccfs)]=ccfs
#         #B[k,length(ccfs):1:end]=zero(T)
#     end
#
#     B
# end


text:
using Snowflakes
using Test
using JuMP
using GLPK

#INITIALIZE DATA
path = joinpath(@__DIR__,"testdata2.xlsx")
extract!(path)

@testset "Initiation" begin
    @test isequal( #all in cover_list is in keys(V)
        sort(stats().cover_list),sort(collect(keys(b().V)))
    )

    idx = rand(collect(keys(b().V))) #random point d[i,i]
    @test b().dist[idx,idx] == 999999999

    @test initStab(;slC = 500.0, suC = -500.0).slLim == abs.(b().d)
    @test initStab(;slC = 500.0, suC = -500.0).suLim == abs.(b().d)

    @test isa(root(;slC = 500.0, suC = -500.0),Snowflakes.node)
end

@testset "Optimizer Setting" begin
    set_default_optimizer!(GLPK.Optimizer)

    @test get_default_optimizer() == GLPK.Optimizer
    @test reset_default_optimizer() == nothing
end

@testset "ColGen Mechanisms" begin
    set_default_optimizer!(GLPK.Optimizer)
    test_root = root(;slC=500.0,suC=-500.0)
    test_mp = master(test_root;silent=false)

    @test has_values(test_mp) #has primal values
    @test has_duals(test_mp) #has dual values

    test_duals = getDuals(test_mp)
    test_sp = sub(test_root,test_duals;silent=false)

    @test colGen(test_root;silent=true,track=true,maxCG=Inf).status[end] == "EVALUATED"

    @test objective_value(test_sp) < 0 #first iter sub < 0
    @test updateStab!(test_root.stab,0.5).slLim == 0.5 .* b().d
    @test typeof(origin(test_root)) == Snowflakes.col
end


text:
#!/usr/bin/env julia

const USAGE = """
USAGE:
./publish.jl [--skip-tests] version

Publish NPM packages at the specified version (after running tests), update
the version in Project.toml, and upload to GitHub.
"""

function errorwithusage()
    @error USAGE
    exit(1)
end

using Pkg
using Pkg.TOML
using LibGit2: GitRepo, branch, isdirty, tag_list

const REPO = GitRepo(@__DIR__)
const PROJECT_FILENAME = normpath(joinpath(@__DIR__, "Project.toml"))
const PROJECT = TOML.parsefile(PROJECT_FILENAME)

if branch(REPO) != "master"
    @error(
        "The publish.jl script must be run on the master branch.\n"
        * "Please `git checkout master && git pull` and re-run.",
    )
    exit(1)
end

if isdirty(REPO)
    @error(
        "The Git repository is dirty. Please make sure all changes are "
        * "commmited and pushed to GitHub. To reset any changes (e.g. if "
        * "publish.jl failed and left the repo dirty), run `git reset --hard`."
    )
    exit(1)
end

const VERSION = let
    if isempty(ARGS) || startswith(last(ARGS), "-")
        errorwithusage()
    end
    VersionNumber(last(ARGS))
end

if "v$(VERSION)" in tag_list(REPO)
    @error(
        "WebIO@$(VERSION) has already been tagged. Did you forget "
        * "to increment the version before running this script?",
    )
    exit(1)
end

if !in("--skip-tests", ARGS)
    @info "Running tests..."
    # Build JS in prod mode.
    ENV["WEBIO_BUILD_PROD"] = true
    Pkg.test("WebIO")

    # Running tests should not make the repo dirty.
    @assert !isdirty(REPO)
end

@info "Publishing WebIO@$(VERSION); this will upload packages to NPM."
print("Continue with publishing [y/N]? ")
let
    confirmation = readline()
    if lowercase(confirmation) != "y"
        @info "Aborting."
        exit(1)
    end
end

const PACKAGES_DIR = normpath(joinpath(@__DIR__, "packages"))
@info "Removing previous build artifacts..."
run(`sh -c "rm -rf ./deps/bundles $(PACKAGES_DIR)/node_modules $(PACKAGES_DIR)/*/node_modules $(PACKAGES_DIR)/*/dist $(PACKAGES_DIR)/package-lock.json $(PACKAGES_DIR)/*/package-lock.json"`)

const npm = `npm -C $PACKAGES_DIR`

@info "Building production (i.e. minified) JavaScript bundles..."
run(`$npm install`)
run(`$npm run lerna -- run build-prod`)

@info "Setting NPM package versions to $VERSION..."
run(`$npm run lerna -- version --no-git-tag-version --no-push --exact --force-publish="*" $VERSION`)

@info "Setting WebIO Project.toml version to $VERSION..."
PROJECT["version"] = VERSION
open(PROJECT_FILENAME, "w") do io
    TOML.print(io, PROJECT)
end

commitmsg = "v$VERSION"
@info "Creating git commit..."
run(`git add .`)
run(`git commit -m $commitmsg`)

print("Please enter your NPM one-time-password (OTP): ")
otp = readline()
run(`$npm run lerna -- publish --otp $otp from-package`)

@info "Making sure package download works..."
download_packages_code = """
    using WebIO
    rm(WebIO.BUNDLES_PATH, force=true, recursive=true)
    WebIO.download_js_bundles()
    """
try
    run(`julia -e $download_packages_code`)
catch exc
    @error(
        "Oh no! Failed to download packages. Packages were already uploaded "
        * "to NPM and can't be overwritten, so you might need to skip this "
        * "version and bump the minor version once the issue is fixed.",
        exception=exc
    )
    rethrow()
end

@info "Pushing to GitHub..."
run(`git push`)

@info(
    "Success!\n"
    * "Comment on the latest commit "
)


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
using Sparso

set_options(SA_ENABLE, SA_VERBOSE, SA_USE_SPMP, SA_CONTEXT, SA_REORDER, SA_REPLACE_CALLS)

#include("./pcg-symgs.jl")
include("utils.jl")

# TEMPORARY: this should be removed in future once structure analysis can
# recognize that L and U are part of A
# Copied from pcg-symgs.jl to avoid affecting other tests if it were changed.
function pcg_symgs(x, A, b, tol, maxiter)
    # These are the only 3 properties we need to discover automatically.
    # We also require SA_CONST_STRUCTURED for A, but that has already been discovered.
#    set_matrix_property(:L, SA_LOWER_OF, :A)
#    set_matrix_property(:U, SA_UPPER_OF, :A)
    set_matrix_property(Dict(
        :A => SA_SYMM_STRUCTURED | SA_SYMM_VALUED,
        :U => SA_CONST_VALUED
        )
    )

    total_time = time()

    trsv_time = 0.
    spmv_time = 0.
    blas1_time = 0.

    L = tril(A)
    U  = SparseMatrixCSC{Cdouble, Cint}(spdiagm(1./diag(A))*triu(A))

    spmv_time -= time()
    r = b - A * x
    spmv_time += time()

    blas1_time -= time()
    normr0 = norm(r)
    z = copy(r)
    blas1_time += time()

    rel_err = 1

    trsv_time -= time()
    Base.SparseMatrix.fwdTriSolve!(L, z)
    Base.SparseMatrix.bwdTriSolve!(U, z)
    trsv_time += time()

    blas1_time -= time()
    p = copy(z) #NOTE: do not write "p=z"! That would make p and z aliased (the same variable)
    rz = dot(r, z)
    blas1_time += time()

    # For statistics only. We get these values so that L, U and A won't be live out of
    # the loop -- to save reverse-reordering of them.
    nnzL = nnz(L)
    nnzU = nnz(U)
    sizeL1 = size(L, 1)
    sizeL2 = size(L, 2)
    nnzA = nnz(A)
    sizeA1 = size(A, 1)
    sizeA2 = size(A, 2)

    k = 1

    while k <= maxiter
        old_rz = rz

        spmv_time -= time()
        Ap = A*p # Ap = Sparso.SpMV(A, p) # This takes most time. Compiler can reorder A to make faster
        spmv_time += time()

        blas1_time -= time()
        alpha = old_rz / dot(p, Ap)
        x += alpha * p
        r -= alpha * Ap
        rel_err = norm(r)/normr0
        blas1_time += time()

        if rel_err < tol
            break
        end

        blas1_time -= time()
        z = copy(r)
        blas1_time += time()

        trsv_time -= time()
        Base.SparseMatrix.fwdTriSolve!(L, z)
          # Could have written as z=L\z if \ is specialized for triangular
        Base.SparseMatrix.bwdTriSolve!(U, z)
          # Could have wrriten as z=U\z if \ is specialized for triangular
        trsv_time += time()

        blas1_time -= time()
        rz = dot(r, z)
        beta = rz/old_rz
        p = z + beta * p
        blas1_time += time()

        k += 1
    end

    total_time = time() - total_time
    println("total = $(total_time)s trsv_time = $(trsv_time)s ($((12.*(nnzL   + nnzU  ) + 2.*8*(sizeL1     + sizeL2    ))*k/trsv_time/1e9) gbps) spmv_time = $(spmv_time)s ($((12.*nnzA   + 8.*(sizeA1     + sizeA2    ))*(k + 1)/spmv_time/1e9) gbps) blas1_time = $blas1_time")
#   println("total = $(total_time)s trsv_time = $(trsv_time)s ($((12.*(nnz(L) + nnz(U)) + 2.*8*(size(L, 1) + size(L, 2)))*k/trsv_time/1e9) gbps) spmv_time = $(spmv_time)s ($((12.*nnz(A) + 8.*(size(A, 1) + size(A, 2)))*(k + 1)/spmv_time/1e9) gbps) blas1_time = $blas1_time")

    return x, k, rel_err
end

if length(ARGS) == 0
    m = 10
    A = generate_symmetric_nonzerodiagonal_sparse_matrix(m)
else
    A = matrix_market_read(ARGS[1], true, true)
    m = size(A, 1)
end
x       = zeros(Float64, m)
b       = ones(Float64, m)
tol     = 1e-7
maxiter = 20000

x, k, rel_err = pcg_symgs(x, A, b, tol, maxiter)
x = zeros(Float64, m)
println("Original: ")
x, k, rel_err = pcg_symgs(x, A, b, tol, maxiter)
println("\tOriginal k=", k)
println("\tOriginal rel_err=", rel_err)

A2 = copy(A) # workaround that we change A in-place
@acc x, k, rel_err = pcg_symgs(x, A2, b, tol, maxiter)
println("\nAccelerated: ")
x = zeros(Float64, m)
@acc x, k, rel_err = pcg_symgs(x, A, b, tol, maxiter)
println("\tAccelerated k=", k)
println("\tAccelerated rel_err=", rel_err)

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso